Partial Regularity for Harmonic Maps into Spheres at a Singular or Degenerate Free Boundary

We prove partial regularity of weakly stationary harmonic maps with (partially) free boundary data on manifolds where the domain metric may degenerate or become singular along the free boundary at the rate $d^\alpha$ for the distance function $d$ from the boundary.


Introduction
The regularity of harmonic mappings of Riemannian manifolds with free or constrained boundary data has been increasingly analysed in recent years in view of the connection between these maps and recent definitions of fractional harmonic mappings of Riemannian manifolds, see for example [14,15,17,21]. In order to obtain (partial) regularity for fractional harmonic maps, one is lead to the study of free boundary harmonic maps on domains where the Riemannian metric may degenerate or become singular along part of the boundary of the domain, depending on the fractional power in question. The regularity of free boundary harmonic maps, with smooth bounded metrics on the domain, was studied substantially earlier and is of independent interest. For example, free boundary harmonic maps constitute generalisations of minimal surfaces with free boundary. As pointed out in [9], they are also connected with critical points of partially linearised models of elasticity with free boundary conditions and can also arise in the theory of liquid crystals when modelling surface effects near a solid boundary interface. Motivated by their connection with fractional harmonic maps and the possible other geometric and physical applications, our goal is to establish a partial regularity theory for a class of free boundary harmonic mappings on domains where the metric is conformal to a smooth Riemannian background metric with conformal factor blowing up or vanishing at the boundary. We state some of our results for general target manifolds, but our main theorems concern sphere valued maps.
Let M be a smooth Riemannian manifold of dimension m ≥ 3 with smooth non-empty boundary, equipped with a smooth, bounded Riemannian metric g. Suppose N is a smooth, compact Riemannian manifold which we may assume is isometrically embedded in R n for some n ∈ N due to the theorem of Nash [18]. Let d := dist(·, ∂M). We are particularly interested in Riemannian metrics h on M of the form h = d α g in a neighbourhood of ∂M for a fixed α ∈ (− 2 m−2 , 2 m−2 ); we will analyse the regularity of a class of critical points v : (M, h) → R n of the Dirichlet energy on (M, h) given by Here D is the differential of v, β = α(m−2) 2 ∈ (−1, 1), and |Dv| 2 g is the energy density with respect to g. In local coordinates it is written as |Dv| 2 g = g ij ∂ i v, ∂ j v where ·, · is the Euclidean metric and (g ij ) = (g ij ) −1 is the matrix representing g −1 , and dvol g is the volume form which in local coordinates satisfies dvol g = det(g)dx. Throughout the paper, we assume the convention that repeated indices indicate summation over the appropriate range (unless the index is m, which is fixed).
Observe that when β = 0 the energy E 0 is nothing other than the Dirichlet energy on the manifold (M, g). In particular, if we were to allow m = 2 then, regardless of the power α, we have E β ≡ E 0 ; in this case all the results in this article follow from known results for free boundary harmonic maps, which we will discuss subsequently.
In [21] the second author established partial regularity up to the boundary for minimisers v : M → N of E β with respect to a free boundary condition in the case where M is a Euclidean half-space. In this article we will establish partial regularity for maps which are merely critical points of E β with respect to both outer and inner variations and still with a free boundary condition. Henceforth we refer to such maps as called weakly stationary harmonic maps with free boundary data.
If we have free boundary data on the whole natural boundary of a manifold, we may find that only constant solutions to the problem exist. As we are studying only the free boundary here, however, we simply exclude any part where the boundary data are not free. Thus we may imagine that M is part of a larger manifold M ′ and v is prescribed on M ′ \ M, but everywhere on ∂M, it is free. As a consequence, we do not assume that M is complete.
For reasons discussed subsequently, we will also focus on the case where N = S n−1 ⊂ R n is the round unit sphere. An abridged version of our partial regularity theorem states the following. Here, int(M) = M \ ∂M is the interior of M and H t is the t-dimensional Hausdorff measure. The theorem should properly be formulated for maps in a weighted Sobolev space, with weights behaving like d β near the boundary, but we leave the technical details for later.
Away from the boundary ∂M, the regularity stated in the Theorem 1.1 follows from known regularity theory for stationary harmonic maps; these maps are critical points, with respect to both outer and inner variations of the Dirichlet energy for mappings between Riemannian manifolds. Hélein established that weakly harmonic maps on two-dimensional domains are smooth [11]. However, when the dimension of the domain m satisfies m ≥ 3, harmonic maps are known to have singularities (points of discontinuity) in general and may even be discontinuous on the whole of their domain; Rivière has constructed discontinuous harmonic maps into spheres in [19]. In order to achieve even partial regularity, additional constraints on harmonic maps are required. One way to proceed is to impose geometric constraints, such as the co-domain having non-positive sectional curvature. In this case, harmonic maps are also known to be smooth, see [23] for example. In order to deal with general codomain manifolds however, one must impose a different type of constraint. Typically one considers stationary harmonic maps, or at least maps which satisfy a monotonicity inequality for the re-scaled energy. For such maps it is possible to estimate the size of the singular set; stationary harmonic maps are smooth in the interior of their domain with the possible exception of a set of vanishing (m − 2)dimensional Hausdorff measure. This was first established for sphere valued maps by Evans [5] and then extended to more general codomain manifolds N by Bethuel [2].
Near the boundary, when β = 0 the regularity stated in Theorem 1.1 follows from the existing regularity theory for stationary free boundary harmonic maps. As mentioned previously, the regularity theory for such maps was originally investigated with geometric or physical motivations in mind. The regularity of free boundary harmonic maps has been investigated by Baldes [1] and Gulliver and Jost [9] who established full regularity of free boundary harmonic maps under certain geometric constraints, such as the image of the harmonic map lies in a sufficiently small geodesic ball. Scheven later proved partial regularity of free boundary stationary harmonic maps into general target manifolds [22] using a reflection construction in the codomain.
When β = 0, Theorem 1.1 is new. There are related results of the second author [21] (studying minimisers of E β ) and of Millot, Pegon and Schikorra [14] who consider a class of stationary or minimising critical points v : Millot, Pegon and Schikorra use the maps they considered as means to establish partial regularity of a type of fractional s-harmonic maps into spheres. The maps they consider satisfy div(x β m Dv) = 0 in a half-space R m−1 × (0, ∞) together with a non-linear Neumann boundary condition. In contrast we require the maps we consider to satisfy v(M) ⊂ N , which makes the condition v(O) ⊂ N redundant. Such maps satisfy a harmonic map equation in the interior of the domain and a homogeneous Neumann-type condition on the boundary, see Section 3.
We will prove Theorem 1.1 by establishing an ε-regularity result, see Theorem 9.1, and then use a covering argument to conclude the partial regularity stated in Theorem 1.1. The key estimate required to prove Theorem 1.1 is a Caccioppoli-type inequality, see Lemma 8.2. For minimisers this kind of inequality is established using comparison maps and we used this approach to prove partial regularity for minimisers of E β in [21]. When considering stationary harmonic maps, it is no longer possible to take advantage of the minimising property and a different approach is required; we will instead take advantage of compensated compactness phenomena in the Euler-Lagrange equation.
The underlying observation is that for harmonic maps into a sphere, the expression |Dv| 2 can be written as a product with a 'div-curl' structure. This idea goes back to Hélein, who used it to show regularity of weakly harmonic maps from surfaces into spheres [12]. Evans [5] extended the results of Hélein to higher dimensional domains. The availability of such methods is the reason why we concentrate on sphere-valued maps as well. The regularity theory for the usual harmonic maps was further extended to other target manifolds by Bethuel [2], also using concentrated compactness arguments but with different technical details. There is also an alternative approach due to Rivière [20]. There is no reason a priori why these methods may not be adapted to our setting, too. However, as we wish to avoid the technical difficulties that come with greater generality and concentrate on the new challenges stemming from a Riemannian metric of the form d α g, we mostly work with the sphere here.
One of the cornerstones of the regularity theory for harmonic maps is the duality between BMO-spaces and the Hardy spaces H 1 (R m ), combined with estimates for the corresponding norms (see, e.g., the work of Coifman-Lions-Meyer-Semmes [4] or of Chanillo [3], who bypasses this duality but obtains inequalities in the same spirit nevertheless). We will require similar estimates, but adapted to weighted Sobolev spaces, which do not appear in the existing literature. The derivation of suitable inequalities will therefore occupy a significant portion of this paper.
In the aforementioned regularity theories for harmonic, and fractional harmonic, maps it is known that stationary or minimising maps are smooth away from their singular set. In the context where these regularity theories were developed, higher regularity follows from known estimates for continuous or Hölder continuous solutions to semi-linear elliptic equations, see respectively [13] and [24] for example. In [21] the second author studied minimising free boundary harmonic maps in the case where (M, g) is a half-space with the Euclidean metric. He used related ideas and some additional arguments to show smoothness in directions tangential to free boundary despite the singular factor arising from the factor d β in the energy. In this article, the anisotropy of the metric g precludes us from directly applying the theory in [21] unless g is Euclidean (in which case we may may draw conclusions on the smoothness of the maps we consider in directions tangential to the boundary as well). Once Hölder continuity is known, we expect minor modifications of the higher regularity results in [21] to yield smoothness for the stationary free boundary harmonic maps we consider here. In particular, we expect all the conclusions of Theorem 4.3 of [21] pertaining to the regularity of the maps considered there, to hold for the maps we consider here (with the exception of the reduced dimension of the singular set in the interior of the domain). However, we will not give the details since we expect all the necessary modifications to be of a purely technical nature. We furthermore observe that there are emerging regularity theories for semi-linear degenerate elliptic equations of the form we consider here, see for example [27]. We note that we could prove a version of Lemma 4.23 in [21] showing Hölder continuity implies Lipschitz continuity and then the results of [27] would apply to give C 1,γ regularity and possibly even C ∞ regularity depending on the precise assumptions required. For the aforementioned reasons, we will not discuss higher regularity any further in this article.
We now outline the organisation of the paper and the strategy for the proof, further details and references may be found at the beginning of each section. In Section 2 we introduce local coordinates such that the factor d β becomes the weight x β m near the boundary; we also introduce the function spaces we require. As discussed previously, the range of specified α implies β ∈ (−1, 1). It is then well known that the weight x β m then lies in Muckenhoupt class A 2 , see Section 7 for the definition, which allows us to apply weighted counterparts to many results, such as the Sobolev-Poincaré inequality, which hold for unweighted Sobolev spaces. In Section 3 we give the full definition of stationary harmonic maps with free boundary data and state the corresponding Euler-Lagrange equations. In Section 4 we use the criticality with respect to inner variations to establish energy monotonicity formulas for the re-scaled energy on concentric balls for maps into general target manifolds. We first do this in the coordinates in Section 2 and then show that monotonicity inequalities hold independently of the choice of coordinates. We prove control of the L 2 norm of the average difference between v and its average provided the re-scaled energy is assumed sufficiently small in Section 5, again for general target manifolds. From Section 6 onwards we specialise to the case where the target manifold is a Euclidean sphere. In this case the Euler-Lagrange equations have a similar structure to the unweighted case and we show that d β |Dv| 2 can be written in the form of a weighted 'div-curl' structure. In Section 7 we establish weighted versions of the aforementioned compensated compactness results, both globally and locally. We then use the weighted compensated compactness arguments, together with the structure of the Euler-Lagrange equations to establish a Caccioppoli-type inequality for the maps we consider in Section 8. With this in hand, it is possible to prove that the re-scaled energy at the boundary decays faster than obtained in the monotonicity formulas described in Section 4. We combine this fact with known results in the interior for stationary harmonic maps to obtain ε-regularity and partial regularity in Section 9.
Convention. We introduce the convention that unless stated otherwise we write C to denote a constant which depends only on universal factors such as m, β or other uniform constants introduced in the next section and we will not always distinguish between different such constants.

Local Coordinates and Function Spaces
Many of our arguments will make use of local coordinates. As we primarily focus on an analysis of critical points of E β near part of the boundary ∂M, it will be convenient to choose local coordinates that map a piece of ∂M to R m−1 × {0} and at the same time preserve the distance (with respect to g) to the boundary.
Consider a coordinate chart is a relatively open set. We may choose an open set U ⊂ M such that U ′ = U ∩ ∂M and such that the nearest point projection π ∂M : U → ∂M with respect to g is well-defined and smooth in U. Then the distance function d is smooth as well. We define x : U → R m by x(q) = (x ′ (π ∂M (q)), d(q)), q ∈ U i.e., x = (x ′ • π ∂M ) × d. This gives rise to local coordinates x on U. If we set U = x(U), then for x ′ ∈ R m−1 , the lines ({x ′ } × R) ∩ U correspond to geodesics in U intersecting the boundary perpendicularly, and distances with respect to g are preserved by x along these geodesics. If (g ij ) and (h ij ) denote the components of g and h, respectively, in these local coordinates, then Some of our arguments do not require anything more. However, we will sometimes need to take advantage of some additional properties satisfied by these coordinates. For example, we will prove that a variant of the well-known monotonicity formula for harmonic maps near points of the boundary holds in a uniform way.
Note that the background metric g induces a Riemannian metric g ′ on ∂M by restriction (even though h does not). Consider a point p ∈ ∂M and suppose that x ′ : U ′ → R m−1 describes normal coordinates on ∂M about p with respect to g ′ . If we use the above construction to obtain local coordinates x : U → R m with this choice of x ′ , then we automatically find that As everything is smooth, it follows that there exist ρ p > 0 and c p > 0 such that |g ij (x) − δ ij | ≤ c p |x| and |∂ k g ij | ≤ c p for i, j, k = 1, . . . , m, (2.2) as long as |x − x 0 | ≤ ρ p . Moreover, if B g r (p) denotes the ball of radius r about p with respect to the metric g, then for any r < ρ p , This gives a useful comparison between balls with respect to g on the one hand and the preimages of Euclidean balls with respect to x on the other hand. The functions p → ρ p and p → c p may be chosen continuous. Hence for any compact set K ⊂ U, there exist ρ > 0 and c > 0 such that ρ p ≥ ρ and c p ≤ c for all p ∈ K ∩ ∂M. The above inequalities are then uniform in K. Similarly, there exist C 0 , . . . , C 3 such that the inequalities for every ξ ∈ R m , and for k = 1, . . . , m, hold uniformly in K even if x is not based on normal coordinates. We are interested in the regularity of a class of critical points v : M → N of E β such that the energy is, at least locally, finite. We first observe that in the coordinates we have introduced on U the energy takes the form In order to give meaning to the notion of critical points of E β we consider Sobolev spaces on M such that this expression is, at least locally, finite. We therefore consider function spaces with respect to the weight the spaces L 2 β and W 1,2 β are Hilbert Spaces with inner products which induce the respective norms. We further define W 1,p β,0 (U ; R n ) as the closure in W 1,p β (U ; R n ) of the space of smooth compactly supported functions C ∞ 0 (U ; R n ) on U with respect to || · || W 1,p β (U ;R n ) . Since |x m | β is an A 2 weight, as defined subsequently in (7.1), the space of smooth functions in W 1,p β (U ; R n ) is dense in W 1,p β (U ; R n ) for p ≥ 2, and even for p > q ≥ 1 for some 2 > q ≥ 1, but not necessarily all p ∈ [1, ∞), see [30] Corollary 2.1.6. When β = 0, we omit the subscript in the preceding notation.
Since we are interested in manifold valued maps, we now introduce the function spaces we will use to analyse manifold valued critical points of E β . We define W 1,2 β (U ; N ) as the collection of maps in W 1,2 β (U ; R n ) with values in N for Lebesgue almost every x ∈ U . We may then define W 1,2 β,loc (M; N ) as maps such that v • x −1 ∈ W 1,2 β (K ; N ) for any local coordinates x : U → U as above and any compact set K ⊂ U , as well as v ∈ W 1,2 (K; N ) for every compact set K ⊂ int(M).
We will have occasion to consider a reflection of v ∈ W 1,2 β (U ; R n ) and of the metric g in the hyperplane We also extend g to V by the obvious reflection, giving rise tog on V , and note thatg also satisfies (2.4) for the same constants in all of V . It is worthy of note that in the coordinates specified previously in this section, the coefficients g im = g mi satisfy g im (x ′ , 0) = 0 for i < m and hence so do the coefficients g im (x ′ , 0) = g mi (x ′ , 0) for i = m. This means thatg andg −1 are continuous, and even Lipschitz continuous in compact domains where they are defined. We further note that the eigenvalues of g −1 are given by λ k (x ′ , |x m |) where λ k are eigenvalues of g −1 .

The First Variation of E β
In this section we specify the class of critical points we consider. We give the definitions in general but only calculate the Euler-Lagrange equations in a chart U as specified in Section 2 since this will suffice for the subsequent analysis.
Let π N denote the nearest point projection on N and let v ∈ W 1,2 Observe that K may intersect the boundary of M. We first consider outer variations of v given by Now let U be as in Section 2 with the given coordinates x : gives rise to an admissible outer variation of the form v t := π N (v + tψ). Following the computations of e.g [26] Chapter 2, we take the first derivative of E β at t = 0 to see that in the chosen coordinates critical points of E β with respect to outer variations satisfy This is the weak form of the equation in U with Neumann boundary condition where div g , grad g , and tr g are the divergence, gradient, and trace with respect to g, respectively.
The following observation will be useful when we choose certain test functions for (3.2).
Proof. It is well-known that |x m | β is a Muckenhoupt weight of class A 2 , see (7.1) for the definition. By Corollary 2.1.6 in [30] Chapter 2 this implies that converging to ψ can be constructed by convolution with a standard mollifying kernel. Then we also have ψ k → ψ pointwise almost everywhere (possibly after passing to a subsequence) and the sequence follows from the Dominated Convergence Theorem. The statement of the lemma now follows immediately.
We now consider inner variations of E β . Let K ⊂ M be compact and let Φ t : M → M constitute a smooth 1-parameter family of diffeomorphisms such that Φ t (p) = p for all p ∈ M \ K and all values of t. Then v t := v • Φ t is an inner variation. Note that K may include boundary points again. We say that v is weakly stationary harmonic with respect to free boundary data if it is weakly harmonic with respect to free boundary data and ∂ t | t=0 E β | K (v t ) = 0 for all inner variations of the preceding form.
Consider again local coordinates x : U → U as in Section 2.
where t is small enough to make Ψ t a diffeomorphism of U . We apply the above definition to x. An analysis of Ψ t similar to the standard theory for harmonic maps (see, e.g., [26] again) then shows that a critical point of E β with respect to inner variations satisfies where again is the divergence of ψ with respect to g. If any of the criteria defining weakly or stationary harmonic maps with respect to free boundary data are satisfied but with variations only supported away from the boundary, we refer to the maps in question as correspondingly weakly harmonic or weakly stationary harmonic, which is standard terminology.

Monotonicity Formula
Here we prove a monotonicity formula for half-balls centred on the boundary, similar to that for stationary harmonic maps, using arguments analogous to those of Große-Brauckmann [8]. Instead of using normal coordinates, we will use those discussed in Section 2. The second author proved a version of the following formula in [21], but the present situation is more general in view of the metric g. Henceforth we define The following computations take place in local coordinates x : U → U as discussed in Section 2. We now fix these local coordinates and work in U , regarding g as a Riemannian metric on U . We also fix a point x 0 ∈ U ∩ (R m−1 × {0}) and assume that g ij (x 0 ) = δ ij and that (2.4) and (2.5) are satisfied.
Lemma 4.1. There exists a constant C > 0 such that the following holds. Suppose Then for every In particular the map Proof. We may assume without loss of generality that Now let χ be a smooth cutoff function with χ ≡ 1 in [1, ∞), χ ≡ 0 in (−∞, 0] and such that 0 ≤ χ ≤ 1 and |χ ′ | ≤ 3. Define η l (x) = χ(l(ρ − |x|)). Replacing η with η l in (4.1) and letting l → ∞ using Lebesgues Dominated Convergence Theorem and Lebesgue's Differentiation Theorem together we see that for almost every ρ < R, where dS g (ρ) = det(g)dH m−1 is the surface measure with density det(g) relative to the (m − 1)-dimensional Hausdorff measure.
. Furthermore, by assumption g ij (0) = δ ij and thus we have that |g ij (x)x j − x i | ≤ Cρ 2 . It follows that for almost every ρ < R for almost every ρ < R, where C is a suitable constant. Now we multiply by (1 + Cρ) −1 ρ 1−m−β e Cρ and integrate between 0 < t < r < R to conclude the proof.
In addition to the monotonicity formula from Lemma 4.1, we can make use of the monotonicity formula derived by Große-Brauckmann [8] for balls in the interior of M. As a consequence, if we have control of the energy in some ball, then we have control of the energy in smaller balls contained therein, too. We now briefly switch our focus to geodesic balls B g r (p) in M with respect to g in order to emphasise that the corresponding inequalities are independent of the choice of local coordinates. To simplify the notation, we introduce the function Corollary 4.1. Suppose that K ⊂ M is compact. Then there exist two constants C, R > 0 such that the following holds true. Let v ∈ W 1,2 loc (M; N ) be a weakly stationary harmonic map with respect to free boundary data. Then for any p ∈ K and for 0 < t ≤ r ≤ R Proof. We observe that ω(p, r) is of order r m if p is far away from the boundary and of order r m+β for p ∈ ∂M. Thus in both cases, the scaling factors in front of the integrals are consistent with Große-Brauckmann's monotonicity formula and Lemma 4.1, respectively. More generally, we can estimate throughout K for some constant c.
If 4t > r, then the inequality from the lemma is immediate. Therefore, we assume that 4t ≤ r. Choose R so small that for every p ∈ K ∩ ∂M, the ball B 4R (p) is contained in a coordinate chart as in Section 2.
We first consider p ∈ K ∩ ∂M (unless the set is empty). We choose local coordinates x : U → U with the properties of Section 2, including (2.2), and set x 0 = x(p) ∈ U . Then Lemma 4.1 applies. Using the notation from the lemma, we conclude that x β m |Dv| 2 g dvol g for some constant C independent of p. Here we work on Euclidean balls depending on the choice of local coordinates. However, according to the discussion in Section 2, in particular by (2.3), if we choose R sufficiently small, then for any s ≤ R. Then x β m |Dv| 2 g dvol g as required.
If p ∈ K satisfies r/2 ≤ d(p) ≤ R, then we can still work in local coordinates as above, but constructed about a point p ′ ∈ ∂M. We set x 0 = x(p) again, which is now in the interior of U . We consider the Riemannian metrich(x) = r −α h(x 0 + rx). Thenh and all of its derivatives are uniformly bounded in a Euclidean ball B r 0 (0) for a fixed radius r 0 > 0 that is independent of p. The map x → v(x 0 + rx) is weakly stationary harmonic with respect to the metrich, and we can use Große-Brauckmann's monotonicity formula in B g r/4 (p) to derive the estimate For points with d(p) ≥ R, the behaviour of the metric near the boundary becomes completely irrelevant. In this case, we can use the usual arguments for harmonic maps.
We hence consider the case when d(p) < r/2. Here we choose p ′ ∈ ∂M with dist(p, p ′ ) = d(p) and we use the inequalities that we already know. Set s = d(p)/2. If 3t ≤ s, we conclude that If 3t > s, then an inequality similar to the first line holds trivially and we can still use the rest of the argument.
When we work with local coordinates, the following version of Corollary 4.1 is convenient. Here we fix the same local coordinates as at the beginning of this section again, and we regard g as a Riemannian metric on the domain U .
is a weakly stationary harmonic map with respect to free boundary data. Let y ∈ B R/2 (x 0 ) and r > 0 such that B r (y) ⊂ B R (x 0 ). If y m ≤ r, then If y m ≥ r, then Proof. This follows from Corollary 4.1 and the fact that for the given local coordinates, Euclidean balls are comparable to geodesic balls with inclusions similar to (2.3).
Using these monotonicity inequalities we can deduce an estimate in the space of bounded mean oscillation for a stationary weakly harmonic map in terms of the re-scaled energy. We recall here that for open Ω ⊂ R m , the mean oscillation of a map v : Ω → R n in Ω is We also introduce the notation |Ω| β = Ω |x m | β dx and v Ω,β = 1 |Ω| β Ω |x m | β vdx. In the following lemma, in addition to the map v on U , we also consider the even reflectionṽ on V = U ∪ U − as defined in Section 2.
is a weakly stationary harmonic map with respect to free boundary data. Then (4.4) Proof. We let B r (y) ⊂ B R/2 (x 0 ) with centre y = (y ′ , y m ) ∈ R m−1 × R. We may assume without loss of generality that y m ≥ 0. If y m ≥ r, then the Poincaré inequality, combined with Corollary 4.2, implies that For the case y m < r, we need to replace the usual Poincaré inequality with something else. Here we note instead that We apply Corollary 4.2 to complete the proof.

Control of the L 2 -norm
Here, and for most of the rest of the paper, we work in local coordinates as in Section 2 again. We establish control of the quantity for certain constants c and suitable radii r, for a solution v to the Euler-Lagrange equations for E β in terms of the energy in a larger ball, provided the latter energy is assumed small. The control we obtain will allow us to control an L 2 term which appears in the Caccioppoli type inequality in Section 8 provided the energy is sufficiently small. We establish this control for general target manifolds using a blow-up argument.
First we prove an auxiliary estimate on the energy decay for solutions to the associated linear equations with homogeneous Neumann-type boundary conditions.
. Then for any compact set K ⊂ U there exist constants C, C ′ , γ > 0 with the following properties. Let y ∈ (R m−1 ×{0}) and ρ > 0 such that for every ψ ∈ C ∞ 0 (U ; R n ). Then there exists λ ∈ R n such that We extend v to V by even reflection. We also extend G to V , using even reflection of G ij for 1 ≤ i, j ≤ m − 1 or i = j = m and odd reflection otherwise. Then v is a weak solution of the equation ρ . It is well-known that |x m | β is a Muckenhoupt weight of class A 2 , see (7.1) for the definition, which yields density of In the compact set K (and in its reflection), we have a uniform upper bound for G and a uniform, positive lower bound for the lowest eigenvalue of G. Hence the first inequality in the statement follows. For the second inequality, it now suffices to choose λ = v(y) and verify that v is Hölder continuous in K . To this end, we apply the regularity theory of [6], in particular Theorem 2.3.12 from [6], to equation (5.3) in V . Now we prove an estimate for the L 2 -norm of weak solutions of the Euler-Lagrange equations for E β . Recall the notation v Ω,β = Ω v|x m | β dx Ω |x m | β dx for Ω ⊂ R m . The following result applies to fixed local coordinates as in Section 2, as usual, and therefore we assume that U ⊂ R m−1 × [0, ∞) is relatively open and g is a Riemannian metric on U .
Proof. The proof of the lemma is based on a blow-up procedure, analogous to that of the proof of Lemma 3.5 in [16], for example. Suppose, for a contradiction, that there exist δ, C * > 0 such that the claim is false. Then for any θ 0 ∈ (0, 1 8 for every ψ ∈ C ∞ 0 (U ; R n ) and a sequence of half-balls B + R k (y k ) ⊂ K such that R k → 0 + and Consider the normalised sequence (w k ) k∈N with w k ∈ W 1,2 where C is a constant depending on the Riemannian metric and on K . Furthermore, we deduce from (5.7) that Recall that |x m | β is an A 2 weight as defined in (7.1). Hence we may apply (5.8) and the Poincaré inequality for A 2 weights to deduce that (w k ) k∈N is bounded in W 1,2 β (B + 1 (0); R n ). Using Lemma 2.5 of [21], we find a subsequence (w k l ) l∈N which converges weakly in W 1,2 β (B + 1 (0); R n ) and strongly in L 2 β (B + 1 (0); R n ) to a map w ∈ W 1,2 β (B + 1 (0); R n ). Let g k = g(R k x + y k ). Then in view of (5.6) and (5.8) we calculate Recall that K is compact. Hence, taking a convergent subsequence of (y k l ) l∈N and re-indexing, we may assume that y k l →ỹ, whereỹ m = 0. Now we observe that the continuity of g yields g k l (x) → F := g(ỹ) uniformly.
Using the weak convergence of (w k l ) l∈N to w, it follows that for every ψ ∈ C ∞ 0 (B 1 (0); R n ). Hence w satisfies (5.1) from Lemma 5.1 in B + 1 (0) for the constant function G = F −1 det(F ). We apply the lemma to see that there exist C > 0 and γ ∈ (0, 1) such that We also conclude, using strong convergence of (w k l ) l∈N to w in L 2 β (B + 1 (0); R n ) to take limits in (5.9), that Inequalities (5.10), (5.11), and (5.12) finally imply that δ ≤ Cθ 2γ 0 . If θ 0 is chosen sufficiently small, then this is a contradiction.

Maps into Spheres
All of our results in previous sections apply to (free boundary) harmonic maps into a general smooth compact manifold N . In this section and henceforth we specialise to the case where N = S n−1 = {y ∈ R n : |y| = 1} is the round n-sphere embedded isometrically in R n . In this case a weakly harmonic map v ∈ W 1,2 β (U ; S n−1 ) with respect to free boundary data satisfies for every ψ ∈ C ∞ 0 (U ; R n ). In order to establish partial regularity, as in the case for stationary harmonic maps into spheres, we take advantage of the fact that the Euler-Lagrange equations can be rewritten in a form reminiscent of a conservation law. We write (6.1) essentially in 'div-curl' form but we also have boundary conditions coming from the free boundary.
if and only if the vector fields X ab ∈ L 2 −β (U ; R m ) defined in components by for a, b = 1, . . . , n and every ψ ∈ C ∞ 0 (U ; R).
Remark. Lemma 6.1 asserts that v ∈ W 1,2 β (U ; S n−1 ) is weakly harmonic with respect to free boundary data if and only if the vector fields X ab ∈ L 2 −β (U ; R n ) satisfy the equation div(X ab ) = 0 in int(U ) with boundary condition X ab = 0 on U ∩ (R m−1 × {0}) weakly.

Hardy Space Estimate With Weights
In order to establish partial regularity for sphere-valued stationary harmonic maps we follow the approach in [16] Section 3 (based on ideas by Hélein [12] and Evans [5]). In particular, we combine the estimate given by Lemma 5.2 for the mean squared oscillation with a Caccioppoli-Type inequality to show that if the re-scaled energy is sufficiently small, then it decays faster than implied by the monotonicity formula in Lemma 4.1. In order to derive this inequality, we will take advantage of the equation derived in Lemma 6.5 satisfied by v. It implies, roughly speaking, that the product Dv · X lies, after extending by a reflection with respect to R m−1 × {0}, in a Hardy space.
We recall the definition of the Hardy space H 1 (R m ). For more details, see [28], for example. Let F be the collection of all ψ ∈ C ∞ 0 (B 1 (0); R) with sup B 1 (0) |Dψ| ≤ 1. Define ψ s (x) = s −m ψ(s −1 x) for ψ ∈ F and s > 0. For a locally integrable distribution v on R m let ψ s * v = R m ψ s (x − y)v(y)dy and define the maximal function In order to prove our Hardy space estimate, we will need to utilise the theory of A p weights, which we define and discuss presently. Suppose µ and ν are measures and let dν = wdµ for a non-negative locally integrable function w. Then w is in for every ball B ⊂ R m , see [29] Chapter I. For the A p class of weights, many counterparts to estimates in unweighted Sobolev-Spaces remain true. We recall the estimates we need here. Let dx denote the m-dimensional Lebesgue measure. We use the notation W 1,p (U ; wdx) to denote weighted Sobolev spaces comprised of classes of dx measurable functions with weak first order derivatives such that the function and all its first order derivatives are pintegrable with respect to the measure wdx. Recall that smooth functions contained in W 1,p (U ; wdx) are dense provided w is an A p weight, see Corollary 2.1.6 of [30]. We further have the following Sobolev-Poincaré inequality.
We observe that the bound m m−1 + δ may not allow the unweighted Sobolev exponent κ = m m−p , however the bound given in the Lemma is sufficient for our purpose. A more precise bound is given in Chapter 15 of [10].
We define the Hardy-Littlewood maximal operator corresponding to a measure dµ by M µ (f )(x) = sup B∋x 1 µ(B) B |f |dµ. For A p weights, the following version of the Hardy-Littlewood Maximal Theorem holds.
Lemma 7.2 ([29] Theorem 9). Let p > 1 and suppose dν = wdµ where w ∈ A p (dµ). Let M µ denote the Hardy Maximal function associated to µ. Then ||M µ (f )|| L p (dν) ≤ C p ||f || L p (dν) for every f ∈ L p (dν). Now we are in a position to prove our weighted Hardy estimates. We first establish an estimate for the product of a gradient and a quantity which is divergence free in R m . We use · to denote the inner product in R m for brevity.
Remark. We will require certain weight functions to satisfy corresponding A p conditions for the proof of Lemma 7.3 but postpone the proof of these conditions to a subsequent Lemma, see Lemma 7.5.
Proof of Lemma 7.3. The following arguments are inspired by the work of Coifman-Lions-Meyer-Semmes [4].
Let ψ s (x) = s −m ψ( x s ) for s > 0 and where ψ ∈ C ∞ 0 (B 1 (0); R) with |Dψ| ≤ 1. Then we have Now fix q ∈ (2, 2m m−1 ). Then 2 > q q−1 > 1 and 2 > q m−1 m > 1. Furthermore, we check that Hence we may choose γ = 0 such that Let v be the average of v on B s (x) with respect to |y m | γq dy. Furthermore, we write Then, using Hölder's inequality for the Lebesgue measure we see that, Note that since γ satisfies (7.3), the right hand side of (7.5) is finite. Now we observe that which is precisely the quantity in the condition that |y m | γq is an A q (dy) weight (to the power of 1 q ). Since γ satisfies (7.3), we have γ ∈ (− 1 q , 1 − 1 q ) and Lemma 7.5 implies that |y m | γq is an A q (dy) weight. Hence, We apply Lemma 7.1 with p = q m−1 m and κ = m m−1 so that κp = q, noting that our choice of γ and Lemma 7.5 imply |y m | γq is an A q m−1 m (dy) weight. This yields We combine (7.4), (7.5), (7.7) and (7.8) to see that where we use the notation M γq and M −γ q q−1 for the Maximal functions M µ 1 and M µ 2 , as defined prior to Lemma 7.2, corresponding to the measures dµ 1 = |y m | γq dy and dµ 2 = |y m | −γ q q−1 dy respectively. It follows from (7.9) and Hölder's inequality that (7.10) Now we apply Lemma 7.2 to two sets of measures. We write dµ 1 = |x m | γq dy and dν 1 = |x m | β−γq dµ 1 = |x m | β dy. It follows from Lemma 7.5 that |x m | β−γq is an Similarly we write dµ 2 = |x m | −γ q q−1 dx and dν 2 = |x m | γ q q−1 −β dµ 2 = |x m | −β dx. It follows from Lemma 7.5 that |x m | γ q q−1 −β is an A 2 q−1 q (dµ 2 ) weight and hence it follows The conclusion of the Lemma follows from (7.10), (7.11) and (7.12).
Next we prove a local version of Lemma 7.3.
and |Dη| ≤ C 0 r . Extend Dv and X by 0 to R m and let λ = R m ηdx −1 R m ηDv · Xdx. Then η(Dv · X − λ) ∈ H 1 (R m ) and, for some The constant λ is subtracted since it is a necessary condition that R m gdx = 0 for any g ∈ H 1 (R m ), see e.g. [25].
Let v ∈ R denote a constant which will later be chosen to be the average v over an appropriate set. Since X is weakly divergence free in B r (x 0 ), we see that for any x ∈ R m we have For all x ∈ R m , v ∈ R and s > 0 we calculate |v − v||X|dy. (7.14) and, using that |Dη| ≤ C 0 r −1 , we also find We further observe that We combine we combine (7.13) -(7.16) to see that for all x ∈ R m , v ∈ R and s > 0 we have To proceed we consider the cases x ∈ B 2r (x 0 ) and x ∈ R m \B 2r (x 0 ) and recall that v ∈ R is a constant we may choose. Suppose x ∈ B 2r (x 0 ) and let s > 0. For s ≤ r we have r −1 ≤ s −1 and in this case we set v = v Bs(x) in (7.17). If s > r then B s (x) ∩ B r (x 0 ) ⊂ B r (x 0 ) ⊂ B 3r (x) and we set v = v B 3r (x) and observe that |B s (x)| −1 ≤ 3 m |B 3r (x)| −1 . Then we see that for x ∈ B 2r (x 0 ) we have where t = s if s ≤ r and t = 3r if s > r. Now, choose q as in the proof of Lemma 7.3 and γ as in (7.3). Noting that X and Dv are extended by 0 outside B r (x 0 ) we conclude from (7.18) and (7.5)-(7.9) that for every x ∈ B 2r (x 0 ) and s > 0, where M γq and M −γ q q−1 are the maximal functions as described in Lemma 7.3. Now we consider x ∈ R m \B 2r (x 0 ). Then |x − x 0 | ≥ 2r. We observe that when |x − x 0 | ≥ r + s we have ψ s * (η(Dv · X − λ))(x) = 0. Hence we only need to consider the case when r < s and 2r ≤ |x − x 0 | < r + s (if s ≤ r this condition is empty). In this case we observe that s −1 < 2|x − x 0 | −1 . Hence, similarly to the proof of Proposition 1.92 of [25], we take advantage of the fact that R m η(Dv · X − λ)dy = 0, the Mean Value Theorem and the fact that |Dψ| ≤ 1 to see that Br(x 0 ) |Dv · X| + |λ|dy. (7.20) We apply Hölder's inequality to see that for every x ∈ R m \B 2r (x 0 ) and every s > 0. Now we estimate the H 1 (R m ) norm of η(Dv · X − λ). We combine (7.19) and (7.21) to see that Now we observe that |x| −(m+1) = Cr −1 dx = Cr −1 (7.23) and, using Hölder's inequality, that We use (7.10) -(7.12) from the proof of Lemma 7.3, together with the fact that Dv and X are extended by 0 in R m \B r (x 0 ) to see that We combine (7.22)-(7.25) to conclude the proof.
In the following lemma we give conditions for the measures discussed in the proof of Lemma 7.3 to satisfy their respective A p conditions. Lemma 7.5. Let m ≥ 3, q ∈ (2, 2m m−1 ) and β ∈ (−1, 1). Let dy be the Lebesgue measure on R m . Define dµ 1 = |y m | γq dy, and dµ 2 = |y m | −γ q q−1 dy. Then Proof. To establish the Lemma, we show (7.1) is satisfied in each case. In order to do this we consider balls with dist(B r (x); R m−1 × {0}) = |x m | − r ≥ r and balls with 0 ≤ |x m | < 2r separately. The subsequent calculations are similar for all four statements. We give the details for the first one and leave the rest to the reader, that is, we establish that |y m | γq is in A q (dy). Suppose first that |x m | − r ≥ r. Then Now suppose 0 ≤ |x m | < 2r. If we write x = (x ′ , x m ), then B r (x) ⊂ B 3r (x ′ , 0). Then provided γq > −1 and −γ q q−1 > −1, we have Hence we have verified (7.1).

Caccioppoli-Type Inequality and Energy Decay
With the Weighted Hardy estimate in Lemma 7.3 in hand, we may prove a Caccioppoli-Type inequality for stationary harmonic maps free boundary data. This estimate gives control of the re-scaled energy on a half-ball in terms of the mean squared oscillation and re-scaled energy on a half-ball of twice the radius. We will then use this lemma to show the re-scaled energy decays faster than implied by the monotonicity formula stated in Section 4; we will show the decay we obtain is fast enough to imply Hölder continuity in Section 9.
We need the following preliminary lemma with regard extending a BMO(B r (x 0 ))function to BMO(R m ) by multiplying by a cutoff function. The following lemma follows from well-known arguments, we give a proof based on e.g [5] Lemma 4.1 in order to elucidate the dependence of the constants on the cutoff function.
Proof. We assume with no loss of generality that v Br(x * ) = 0.
To proceed, we consider the cases x ∈ B 13 16 r (x * ) and x ∈ R m \B 13 16 r (x * ). First suppose x ∈ R m \B 13 Consequently, using the fact that |Dη| ≤ C * r −1 , we see that |v|dy. There exists a constant C > 0 such that the following holds. Suppose v ∈ W 1,2 β (U ; S n−1 ) is weakly stationary harmonic with respect to free boundary data.
Proof. The proof is analogous to the case for stationary harmonic maps, see for instance Lemma 3.8 of [16]. Let η ∈ C ∞ 0 (B 3 4 r (x 0 ); R) be a smooth cutoff function with η ≡ 1 in B r 2 (x 0 ), 0 ≤ η ≤ 1 and |Dη| ≤ C r . Henceforth we abbreviate the notationṽ Br(x 0 ),β toṽ and v B + r (x 0 ),β to v. Letg be the reflection of g. It follows from (2.4) and the fact that |ṽ| 2 = 1 almost everywhere that Now let w = η(ṽ −ṽ) ∈ W 1,2 β,0 (B r (x 0 ); R n ). We calculate for a, b = 1, . . . , n. We observe that sinceṽ, . Moreover, since the support of η is a compact subset of B r (x 0 ), we see thatṽ a w b −ṽ b w a ∈ W 1,2 β,0 (B r (x 0 ); R). We observe that the preceding discussion for w also holds true for W = ηw, including (8.8) which holds replacing η with η 2 , ∂ i η with 2η∂ i η and w with W . Lemma 6.1 implies that the vector fieldsX ab ∈ L 2 −β (V ; R n ) given in components byX i ab = |x m | β det(g(x ′ , |x m |))g ij (ṽ a ∂ jṽ b −ṽ b ∂ jṽ a ) are weakly divergence free in V . We use this fact, together with the fact thatṽ a W b −ṽ b W a ∈ W 1,2 β,0 (B r (x 0 ); R) and (8.7) and (8.8), to see that We apply Young's inequality to see that x β m |v − v| 2 dx, (8.10) where we have also used the fact that x β m |v − v| 2 dx.
monotonicity inequality, Corollary 4.2, and from [2] and [24] that there is a constant C and a γ 1 which do not depend on ρ or y such that for t ≤ ρ t 2−m Bt(y) From ε-regularity we deduce partial regularity using a covering argument. Let H t denote the t-dimensional Hausdorff measure. The following is the full version of Theorem 1.1.
Theorem 9.2. Let β ∈ (−1, 1). Suppose v ∈ W 1,2 β (M; S n−1 ) is weakly stationary harmonic with respect to free boundary data. There exists a relatively closed set Σ ⊂ M with H m−2+β (Σ ∩ intM) = 0 and H m−2 (Σ ∩ ∂M) = 0 and there exists γ ∈ (0, 1) such that v ∈ C 0,γ loc (M \ Σ; N ). Proof. It suffices to prove the statement in coordinate patches. Moreover, in the interior of M, the statement follows from the known regularity theory for harmonic maps, in particular [5]. Hence we may consider U ⊂ R m−1 × [0, ∞) as in Section 2 and regard g as a Riemannian metric on U , as we did in the last few sections. Define Σ I = {x ∈ U ∩ (R m−1 × (0, ∞)) : v is not smooth in any neighbourhood of x}.