# A weakly coupled system of $p$-Laplace type in a heat conduction problem

Morteza Fotouhi, Mohammad Safdari, Henrik Shahgholian

###### Abstract.

We study temperature distribution in a heat conducting problem, for a system
of p-Laplace equation, giving rise to a free boundary.

###### Key words and phrases:

Free boundary, Heat conduction, p-Laplacian.

###### 2010 Mathematics Subject Classification:

35R35

###### Contents

  1. 1 Introduction
    1. 1.1 Background
    2. 1.2 Structure of the paper
  2. 2 The Penalized Problem
  3. 3 Regularity of solutions to the penalized problem
  4. 4 The Original Problem
  5. 5 Regularity of the free boundary (case $p=2$)

##  1\. Introduction

###  1.1. Background

This paper delves into an extension of a classical optimization problem in
heat conduction, presenting a succinct description as follows: Given a surface
$\partial\Omega$ in $\mathbb{R}^{n}$ and positive constant functions defined
on it (representing the temperature distribution), the aim is to enclose
$\partial\Omega$ with a prescribed volume of insulating material to minimize
heat loss in a stationary scenario. Mathematically, the objective is to
discover a vector-valued function $\mathbf{u}=(u^{1},\cdots,u^{m})$ ($m\geq
1$) that corresponds to the temperature within $\Omega$. Whenever the
components of $\mathbf{u}$ are nonnegative and the volume of its support is
equal to $1$, they become $p$-harmonic. The target is to minimize the heat
flow, which can be regarded as a continuous family of convex functions
dependent on $\nabla\mathbf{u}$ along $\partial\Omega$.

Our research was inspired by a series of papers [1, 2, 3] and their
generalization presented in [12]. The initial two articles focused on studying
constant temperature distributions, specifically in the linear case where
$\Gamma(x,t)=t$. This linear setting enabled [1, 3] to reduce the quantity to
be minimized to the Dirichlet integral. However, even within the linear case,
the problem of nonconstant temperature distribution, examined in [2],
introduced various new challenges.

The main objective of our article is to explore the system version of the
nonlinear case with a nonconstant temperature distribution, wherein the
equation is governed by the $p$-Laplacian. The nonlinearity addressed in this
paper holds significant physical importance, as problems involving monotone
operators, akin to those studied in [12], arise in the optimization of domains
for electrostatic configurations.

The nonlinearity associated with $\nabla\mathbf{u}$ introduces various new
challenges. For instance, computing normal derivatives of $W^{1,p}$-functions
becomes problematic, leading to difficulties in providing a reasonable
mathematical model. In [2], this challenge was overcome by minimizing the
total mass of $\Delta u$, which can be treated as a nonnegative measure when
$u$ is subharmonic. However, in the present case, there is no integral
representation available for
$\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}\mathbf{u}(x)\big{)},d\sigma$. To
address this issue, similar to [12], we solve appropriate auxiliary
variational problems and compare them with the minimizer.

Now let us introduce the problem in mathematical framework. Let
$\Omega\subset\mathbb{R}^{n}$ ($n\geq 2$) be a bounded open set with smooth
boundary whose volume $|\Omega|>1$. Consider the $p$-Laplace differential
operator ($1<p<\infty$)

$\Delta_{p}u^{i}=\mathrm{div}\big{(}|\nabla u^{i}|^{p-2}\,\nabla
u^{i}\big{)}=\mathrm{div}\big{(}A[u^{i}]\big{)},$

where we set $A[u^{i}]=A(\nabla u^{i}):=|\nabla u^{i}|^{p-2}\,\nabla u^{i}$ to
simplify the notation.

Let $\boldsymbol{\varphi}:\partial\Omega\to\mathbb{R}^{m}$ be a $C^{1}$
function with positive components $\varphi^{i}>0$. For
$\mathbf{u}:\Omega\to\mathbb{R}^{m}$ ($m\geq 1$) satisfying

(1.1) $\begin{cases}\Delta_{p}u^{i}=0&\textrm{in }\\{|\mathbf{u}|>0\\},\\\
u^{i}=\varphi^{i}&\textrm{on }\partial\Omega,\\\
\textrm{vol}(\mathrm{spt}\,|\mathbf{u}|)=1,\end{cases}$

we want to minimize the functional

$J(\mathbf{u}):=\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}u^{1}(x),\dots,A_{\nu}u^{m}(x)\big{)}\,d\sigma(x),$

where $\nu$ is the outward normal vector on $\partial\Omega$,

$A_{\nu}u^{i}:=|\nabla u^{i}|^{p-2}\,\partial_{\nu}u^{i},$

and $\Gamma(x,\xi):\partial\Omega\times\mathbb{R}^{m}\to\mathbb{R}$ is a
continuous function that satisfies:

  1. (1)

For each fixed $x$, $\Gamma(x,\cdot)$ is a convex function.

  2. (2)

For every $i$, $\partial_{\xi_{i}}\Gamma(\cdot,\cdot)$ is positive and has a
positive lower bound on any set of the form $\\{(x,\xi):\xi_{i}\geq a\\}$. In
addition, $\partial_{\xi_{i}}\Gamma(\cdot,\cdot)$ is bounded above on any set
of the form $\\{(x,\xi):\xi_{i}\leq b\\}$. (The bounds can depend on $a,b$.)

  3. (3)

For each fixed $\xi$, $\partial_{\xi_{i}}\Gamma(\cdot,\xi)$ is a $C^{1}$
function.

Note that, as a result, for every $\xi$ we have

$\displaystyle\Gamma(x,\xi_{1},\dots,\xi_{m})$ $\displaystyle\geq\sum_{i\leq
m}\partial_{\xi_{i}}\Gamma(x,0)\xi_{i}+\Gamma(x,0)$
$\displaystyle\geq\sum_{i\leq m}\psi_{i}(x)\xi_{i}-C,$

where $\psi_{i}(x):=\partial_{\xi_{i}}\Gamma(x,0)>0$ are positive $C^{1}$
functions and $C$ is a constant. In particular we have

(1.2)
$\Gamma(x,A_{\nu}u_{1},\dots,A_{\nu}u_{m})\geq\sum_{i=1}^{m}\psi_{i}(x)A_{\nu}u^{i}-C.$

A typical example of $\Gamma$ is

$\Gamma(x,\xi)=\psi_{1}(x)\gamma_{1}(\xi_{1})+\dots+\psi_{m}(x)\gamma_{m}(\xi_{m}),$

where $\psi_{i}$’s are $C^{1}$ and positive, and $\gamma_{i}$’s are $C^{1}$
increasing convex functions with positive derivative.

###  1.2. Structure of the paper

The structure of our paper is as follows: In Section 2, we introduce the
physical problem under consideration. We then formulate a penalized version of
the variational problem for the temperature $\mathbf{u}$ and define suitable
constraint sets as part of our strategy to overcome the challenges arising
from the nonlinearity. We solve the optimization problem over weakly closed
subsets of $W^{1,p}$ (the sets $V_{\delta}$), establishing the optimal
regularity properties of the minimizers, including Lipschitz regularity. These
results are crucial for proving the existence of an optimal configuration for
the original penalized problem, as discussed in Section 3. Here we also
present fundamental geometric-measure properties of the optimal configuration,
such as linear growth away from the free boundary and uniformly positive
density. These properties allow us to establish a representation theorem
following the framework of [3].

In Section 4, we recover the original physical problem from the penalized
problem by showing that for sufficiently small $\varepsilon$, the volume of
$\\{|\mathbf{u}_{\varepsilon}|>0\\}$ automatically adjusts to be equal to $1$.

Section 5 is dedicated to the optimal regularity of the free boundary, for the
case $p=2$. We demonstrate that the normal derivative of the minimizer along
the free boundary is a Hölder continuous function, leading to the conclusion
that the free boundary is a $C^{1,\alpha}$ surface. Furthermore, using the
free boundary condition obtained during the proof of Hölder continuity, we
establish that the free boundary is an analytic surface, except for a small
singular set.

##  2\. The Penalized Problem

Let $\Omega_{\delta}:=\\{x\in\Omega:\textrm{dist}(x,\partial\Omega)<\delta\\}$
and

$\displaystyle V_{\delta}:=\\{\mathbf{u}\in W^{1,p}(\Omega;\mathbb{R}^{m}):\;$
$\displaystyle u^{i}\geq 0,\,\Delta_{p}u^{i}\geq 0,$
$\displaystyle\Delta_{p}u^{i}=0\textrm{ in
}\Omega_{\delta},\,u^{i}=\varphi^{i}\textrm{ on }\partial\Omega\\}.$

Furthermore, we set

$V:=\bigcup_{\delta>0}V_{\delta}.$

Observe that the above definition is consistent due to the assumption
$\varphi^{i}>0$ on $\partial\Omega$. Also, by $\Delta_{p}u^{i}\geq 0$ we mean
that for any test function $\zeta\in C_{c}^{\infty}(\Omega)$ with $\zeta\geq
0$ we have

$-\int_{\Omega}\nabla\zeta\cdot|\nabla u^{i}|^{p-2}\nabla u^{i}\,dx\geq 0.$

This implies that there is a Radon measure $\mu^{i}$ such that for any test
function $\zeta\in C_{c}^{\infty}(\Omega)$ we have

$\int_{\Omega}\zeta\,d\mu^{i}=-\int_{\Omega}\nabla\zeta\cdot|\nabla
u^{i}|^{p-2}\nabla u^{i}\,dx.$

To simplify the notation, we denote $\mu^{i}$ by $\Delta_{p}u^{i}$, and
$d\mu^{i}$ by $\Delta_{p}u^{i}\,dx$. (It should be noted that this notation is
not meant to imply $\mu^{i}$ is absolutely continuous with respect to the
Lebesgue measure. In fact, for the minimizer, the two measures are mutually
singular as we will see in Theorem 3.5.)

Let $f_{\varepsilon}:\mathbb{R}\to\mathbb{R}$ be

$f_{\varepsilon}(t):=\begin{cases}1+\frac{1}{\varepsilon}(t-1)&t\geq 1,\\\
1+\varepsilon(t-1)&t<1.\end{cases}$

We are interested in minimizing the penalized functional

$J_{\varepsilon}(\mathbf{u}):=\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}\mathbf{u}(x)\big{)}\,d\sigma+f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}$

over $V$. The significance of the above penalization is that it forces the
volume $\big{|}\\{|\mathbf{u}|>0\\}\big{|}$ to be $1$ for small enough
$\varepsilon$; see Theorem 4.3. (Notice that the components of $\mathbf{u}\in
V$ are $p$-harmonic near $\partial\Omega$; therefore they are smooth enough
near the boundary, and it makes sense to compute their derivatives along
$\partial\Omega$.) We first consider the minimizer of $J_{\varepsilon}$ over
$V_{\delta}$.

######  Lemma 2.1.

Let $\mathbf{u}\in V$. Then we have

(2.1)
$\int_{\Omega}\psi_{i}\Delta_{p}u^{i}\,dx+\int_{\Omega}\nabla\psi_{i}\cdot
A[u^{i}]\,dx=\int_{\partial\Omega}\psi_{i}A_{\nu}u^{i}\,d\sigma,$

where $\psi_{i}$’s are $C^{1}$ functions.

###### Proof.

Let $\phi_{k}\in C^{\infty}(\Omega)$ be such that $\phi_{k}\equiv 1$ on
$\tilde{\Omega}_{k}:=\Omega-\Omega_{1/k}$ and $\phi_{k}\equiv 0$ on
$\partial\Omega$. We know that $\mathbf{u}\in V_{\delta}$ for some $\delta$.
Suppose $k$ is large enough so that $1/k<\delta$, and thus
$\Omega_{1/k}\subset\Omega_{\delta}$. Then we have

$\displaystyle\int_{\Omega}\nabla(\phi_{k}\psi_{i})\cdot A[u^{i}]\,dx$
$\displaystyle=\int_{\Omega_{1/k}}\nabla(\phi_{k}\psi_{i})\cdot A[u^{i}]\,dx$
$\displaystyle\qquad\qquad+\int_{\tilde{\Omega}_{k}}\nabla(\underset{\equiv\,1}{\underset{\rotatebox[origin={c}]{90.0}{\big{\\{}}}{\phi_{k}}}\psi_{i})\cdot
A[u^{i}]\,dx$ $\displaystyle=\int_{\Omega_{1/k}}\nabla(\phi_{k}\psi_{i})\cdot
A[u^{i}]\,dx+\int_{\tilde{\Omega}_{k}}\nabla\psi_{i}\cdot A[u^{i}]\,dx.$

Now, noting that
$\partial\Omega_{1/k}=\partial\tilde{\Omega}_{k}\cup\partial\Omega$, and by
using the integration by parts formula proved in [5], we get

$\displaystyle\int_{\Omega_{1/k}}\nabla(\phi_{k}\psi_{i})\cdot A[u^{i}]\,dx$
$\displaystyle=\int_{\Omega_{1/k}}\nabla(\phi_{k}\psi_{i})\cdot
A[u^{i}]+\phi_{k}\psi_{i}\negthickspace\negthickspace\underset{=\,0\textrm{ on
}\Omega_{1/k}}{\underbrace{\Delta_{p}u^{i}}}\,dx$
$\displaystyle=-\int_{\partial\tilde{\Omega}_{k}}\underset{\equiv\,1}{\underset{\rotatebox[origin={c}]{90.0}{\big{\\{}}}{\phi_{k}}}\psi_{i}A_{\nu}u^{i}\,d\sigma+\int_{\partial\Omega}\underset{\equiv\,0}{\underset{\rotatebox[origin={c}]{90.0}{\big{\\{}}}{\phi_{k}}}\psi_{i}A_{\nu}u^{i}\,d\sigma$
$\displaystyle=-\int_{\partial\tilde{\Omega}_{k}}\psi_{i}A_{\nu}u^{i}\,d\sigma\underset{k\to\infty}{\longrightarrow}-\int_{\partial\Omega}\psi_{i}A_{\nu}u^{i}\,d\sigma.$

In addition, we have $\int_{\tilde{\Omega}_{k}}\nabla\psi_{i}\cdot
A[u^{i}]\,dx\underset{k\to\infty}{\longrightarrow}\int_{\Omega}\nabla\psi_{i}\cdot
A[u^{i}]\,dx$, and

$\displaystyle\int_{\Omega}\nabla(\phi_{k}\psi_{i})\cdot A[u^{i}]\,dx$
$\displaystyle=-\int_{\Omega}\phi_{k}\psi_{i}\,\Delta_{p}u^{i}\,dx$
$\displaystyle\qquad\underset{k\to\infty}{\longrightarrow}-\int_{\Omega}\psi_{i}\Delta_{p}u^{i}\,dx,$

which together give the desired result. ∎

We can similarly show that

$\int_{\Omega}\Delta_{p}u^{i}\,dx=\int_{\partial\Omega}A_{\nu}u^{i}\,d\sigma.$

In addition, note that $\int_{\Omega}u^{i}\Delta_{p}u^{i}\,dx$ is meaningful
(since $u^{i}-\varphi^{i}\in W_{0}^{1,p}$ while $\Delta_{p}u^{i}\in
W^{-1,p/(p-1)}$ and $\varphi^{i}$ is continuous) and we can similarly show
that

(2.2) $\int_{\Omega}u^{i}\Delta_{p}u^{i}+|\nabla
u^{i}|^{p}\,dx=\int_{\partial\Omega}\varphi^{i}A_{\nu}u^{i}\,d\sigma.$

Note that $u^{i}=\varphi^{i}$ on $\partial\Omega$.

######  Lemma 2.2.

For $\mathbf{u}\in V$ we have

$\sum_{i=1}^{m}\int_{\Omega}|\nabla u^{i}|^{p}\,dx\leq
C+C\int_{\partial\Omega}\sum_{i\leq m}\psi_{i}(x)A_{\nu}u^{i}\,d\sigma.$

######  Remark.

As we will see, the above inequality actually holds for each summand.
Furthermore, with a slight modification of the last part of the proof we
obtain that

$\int_{\Omega}\Delta_{p}u^{i}\,dx\leq
C+C\int_{\partial\Omega}\psi_{i}(x)A_{\nu}u^{i}\,d\sigma.$

###### Proof.

Let $\mathbf{h}_{0}$ be the vector-valued function in $\Omega$ satisfying
$\Delta_{p}h_{0}^{i}=0$, and taking the boundary values $\boldsymbol{\varphi}$
on $\partial\Omega$. Note that $h_{0}^{i}$ is $C^{1}$ and we can plug it in
(2.1). By subtracting the resulting relation from (2.2) we get

$\int_{\Omega}(u^{i}-h_{0}^{i})\Delta_{p}u^{i}\,dx+\int_{\Omega}\nabla(u^{i}-h_{0}^{i})\cdot
A[u^{i}]\,dx=0.$

Hence

$\displaystyle\int_{\Omega}|\nabla u^{i}|^{p}\,dx$
$\displaystyle=\int_{\Omega}\nabla u^{i}\cdot
A[u^{i}]\,dx=\int_{\Omega}(h_{0}^{i}-u^{i})\Delta_{p}u^{i}\,dx+\int_{\Omega}\nabla
h_{0}^{i}\cdot A[u^{i}]\,dx$
$\displaystyle\leq\int_{\Omega}h_{0}^{i}\Delta_{p}u^{i}\,dx+C\int_{\Omega}|\nabla
h_{0}^{i}|^{p}\,dx+\frac{1}{2}\int_{\Omega}|A[u^{i}]|^{\frac{p}{p-1}}\,dx$
$\displaystyle\leq
C\int_{\Omega}\Delta_{p}u^{i}\,dx+C+\frac{1}{2}\int_{\Omega}|\nabla
u^{i}|^{p}\,dx,$

where we have used the facts that $u^{i},\Delta_{p}u^{i}\geq 0$ and
$|A[u^{i}]|^{\frac{p}{p-1}}=|\nabla u^{i}|^{p}$. Thus we have

$\int_{\Omega}|\nabla u^{i}|^{p}\,dx\leq C\int_{\Omega}\Delta_{p}u^{i}\,dx.$

But since $\psi_{i},\Delta_{p}u^{i}\geq 0$ we get

$\int_{\Omega}|\nabla u^{i}|^{p}\,dx\leq C\int_{\Omega}\Delta_{p}u^{i}\,dx\leq
CC_{i}\int_{\Omega}\psi_{i}\,\Delta_{p}u^{i}\,dx,$

where $C_{i}=\max_{\overline{\Omega}}\frac{1}{\psi_{i}}>0$. Hence by (2.1) we
get

$\displaystyle\int_{\Omega}|\nabla u^{i}|^{p}\,dx$ $\displaystyle\leq
C\int_{\Omega}\psi_{i}\,\Delta_{p}u^{i}\,dx$
$\displaystyle=-C\int_{\Omega}\nabla\psi_{i}\cdot
A[u^{i}]\,dx+C\int_{\partial\Omega}\psi_{i}A_{\nu}u^{i}\,d\sigma$
$\displaystyle\leq\tilde{C}\int_{\Omega}|\nabla\psi_{i}|^{p}\,dx+\frac{1}{2}\int_{\Omega}|A[u^{i}]|^{\frac{p}{p-1}}\,dx+C\int_{\partial\Omega}\psi_{i}A_{\nu}u^{i}\,d\sigma$
$\displaystyle\leq C+\frac{1}{2}\int_{\Omega}|\nabla
u^{i}|^{p}\,dx+C\int_{\partial\Omega}\psi_{i}A_{\nu}u^{i}\,d\sigma,$

which gives the desired. ∎

######  Theorem 2.3.

There exists a minimizer $\mathbf{u}_{\varepsilon}^{\delta}\in V_{\delta}$ for
$J_{\varepsilon}$.

###### Proof.

Let $\\{\mathbf{u}_{k}\\}\subset V_{\delta}$ be a minimizing sequence. Then by
the above lemma and (1.2) we have

$\displaystyle\sum_{i=1}^{m}\int_{\Omega}|\nabla u_{k}^{i}|^{p}\,dx$
$\displaystyle\leq C+C\int_{\partial\Omega}\sum_{i\leq
m}\psi_{i}A_{\nu}u_{k}^{i}\,d\sigma$ $\displaystyle\leq
C+C\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}u_{k}^{1},\dots,A_{\nu}u_{k}^{m}\big{)}+C\,d\sigma$
$\displaystyle\leq C+CJ_{\varepsilon}(\mathbf{u}_{k}).$

Hence $\|\nabla\mathbf{u}_{k}\|_{L^{p}}$ is bounded. In addition, for the dual
exponent $q=\frac{p}{p-1}$ we can see that
$\|A[\mathbf{u}_{k}]\|_{L^{q}}=\|\nabla\mathbf{u}_{k}\|_{L^{p}}^{p-1}$ is also
bounded. Hence, up to a subsequence, we can assume that $\nabla
u_{k}^{i}\rightharpoonup\nabla u^{i}$ in $L^{p}$, $A[u_{k}^{i}]\rightharpoonup
A[u^{i}]$ in $L^{q}$, and $\mathbf{u}_{k}\to\mathbf{u}$ a.e in $\Omega$. Thus
$u^{i}\geq 0$. Also, $u^{i}=\varphi^{i}$ on $\partial\Omega$, since
$u_{k}^{i}-\varphi^{i}\in W_{0}^{1,p}(\Omega)$, which is a closed and convex
set, hence weakly closed. Finally, to see that $\Delta_{p}u^{i}$ has the
desired properties, notice that for an appropriate test function $\phi$ we
have

$\int_{\Omega}\nabla\phi\cdot
A[u^{i}]\,dx=\lim_{k\to\infty}\int_{\Omega}\nabla\phi\cdot A[u_{k}^{i}]\,dx$

due to the weak convergence of $A[\mathbf{u}_{k}]$. Therefore $\mathbf{u}\in
V_{\delta}$. Now we can repeat the proof of Lemma 3.3 in [12] to deduce the
weak lower semicontinuity of $J_{\varepsilon}$ with respect to this sequence,
and conclude the proof (the convexity of $\Gamma$ is needed here). ∎

######  Lemma 2.4 (Hopf’s lemma for $p$-harmonic functions).

Suppose $h$ is a $p$-harmonic function on $B_{1}(0)$ with nonnegative boundary
values on $\partial B_{1}$. Then we have

$h(x)\geq c(n,p)\,\mathrm{dist}(x,\partial B_{1})\sup_{B_{1/2}}h.$

###### Proof.

Consider the function $g(x)=e^{-\lambda|x|^{2}}-e^{-\lambda}$ for some
$\lambda>0$. Note that $g=0$ on $\partial B_{1}$, and $0<g<1$ on $B_{1}$. We
also have

$\partial_{i}g=-2\lambda
x_{i}e^{-\lambda|x|^{2}},\qquad\partial_{ij}g=(4\lambda^{2}x_{i}x_{j}-2\lambda\delta_{ij})e^{-\lambda|x|^{2}}.$

Now we have $\Delta g=(4\lambda^{2}|x|^{2}-2n\lambda)e^{-\lambda|x|^{2}}$, and

$\displaystyle\Delta_{\infty}g:=\sum_{i,j}\partial_{i}g\partial_{j}g\partial_{ij}g$
$\displaystyle=\sum_{i,j}4\lambda^{2}(4\lambda^{2}x_{i}^{2}x_{j}^{2}-2\lambda\delta_{ij}x_{i}x_{j})e^{-3\lambda|x|^{2}}$
$\displaystyle=4\lambda^{2}(4\lambda^{2}|x|^{4}-2\lambda|x|^{2})e^{-3\lambda|x|^{2}}.$

Therefore

$\displaystyle\Delta_{p}g$ $\displaystyle=\mathrm{div}(|\nabla g|^{p-2}\nabla
g)=|\nabla g|^{p-4}\big{(}|\nabla g|^{2}\Delta g+(p-2)\Delta_{\infty}g\big{)}$
$\displaystyle=(2\lambda|x|)^{p-4}\big{(}4\lambda^{2}|x|^{2}(4\lambda^{2}|x|^{2}-2n\lambda)$
$\displaystyle\qquad\qquad\qquad\qquad+(p-2)4\lambda^{2}(4\lambda^{2}|x|^{4}-2\lambda|x|^{2})\big{)}e^{-(p-1)\lambda|x|^{2}}$
$\displaystyle=(2\lambda)^{p-1}|x|^{p-2}\big{(}2\lambda|x|^{2}-n+(p-2)(2\lambda|x|^{2}-1)\big{)}e^{-(p-1)\lambda|x|^{2}}$
$\displaystyle=(2\lambda)^{p-1}|x|^{p-2}\big{(}2(p-1)\lambda|x|^{2}-n-p+2\big{)}e^{-(p-1)\lambda|x|^{2}}.$

Thus for $\frac{1}{2}\leq|x|\leq 1$ and large enough $\lambda$ we have

$\Delta_{p}g\geq
2\lambda^{p-1}\big{(}(p-1)\lambda/2-n-p+2\big{)}e^{-(p-1)\lambda}>0.$

Now we have $h\geq\inf_{\overline{B}_{1/2}}h>(\inf_{\overline{B}_{1/2}}h)g$ on
$\overline{B}_{1/2}$ (note that $h$ is positive on $B_{1}$ by maximum
principle), and on $B_{1}-\overline{B}_{1/2}$ we have
$\Delta_{p}h=0<\Delta_{p}g$. Also on $\partial B_{1}$ we have $h\geq 0=g$.
Hence by the maximum principle we have $h(x)\geq
g(x)(\inf_{\overline{B}_{1/2}}h)$ for $x\in B_{1}$. But by the Harnack’s
inequality we have

$\inf_{\overline{B}_{1/2}}h\geq C\sup_{\overline{B}_{1/2}}h$

for some constant $C$ which does not depend on $h$. Hence we obtain

$h(x)\geq Cg(x)\sup_{\overline{B}_{1/2}}h.$

On the other hand note that

$\displaystyle g(x)$
$\displaystyle=g(x)-g(x/|x|)=\int_{1/|x|}^{1}\frac{d}{dt}g(tx)\,dt=\int_{1/|x|}^{1}x\cdot\nabla
g(tx)\,dt$ $\displaystyle=\int_{1/|x|}^{1}-2\lambda t|x|^{2}e^{-\lambda
t^{2}|x|^{2}}\,dt=2\lambda|x|^{2}\int_{1}^{1/|x|}te^{-\lambda
t^{2}|x|^{2}}\,dt$ $\displaystyle\geq
2\lambda|x|^{2}\int_{1}^{1/|x|}te^{-\lambda}\,dt=\lambda
e^{-\lambda}|x|^{2}\Big{(}\frac{1}{|x|^{2}}-1\Big{)}=\lambda
e^{-\lambda}(1-|x|^{2})$ $\displaystyle\geq\lambda e^{-\lambda}(1-|x|)=\lambda
e^{-\lambda}\,\mathrm{dist}(x,\partial B_{1}),$

which gives the desired. ∎

If $h$ is a $p$-harmonic function on $B_{r}(x_{0})$, then
$\tilde{h}(x):=h(x_{0}+rx)$ is a $p$-harmonic function on $B_{1}(0)$. Hence we
have

$\displaystyle h(x_{0}+rx)$ $\displaystyle=\tilde{h}(x)\geq
c(n,p)\,\mathrm{dist}(x,\partial B_{1})\sup_{B_{1/2}}\tilde{h}$
$\displaystyle=c(n,p)\,(1-|x|)\sup_{B_{1/2}}\tilde{h}=c(n,p)\,\frac{(r-r|x|)}{r}\sup_{B_{r/2}(x_{0})}h$
$\displaystyle=c(n,p)\,\mathrm{dist}\big{(}x_{0}+rx,\partial
B_{r}(x_{0})\big{)}\frac{1}{r}\sup_{B_{r/2}(x_{0})}h.$

######  Lemma 2.5.

Let $w\in W^{1,p}(\Omega)$ be a nonnegative function. Then there exists $c>0$,
depending only on $p$ and the dimension, such that for any ball
$\overline{B}_{r}(x_{0})\subset\Omega$ we have

$\Big{(}\frac{1}{r}\sup_{B_{r/2}(x_{0})}h\Big{)}^{p}\cdot\big{|}B_{r}(x_{0})\cap\\{w=0\\}\big{|}\leq
c\int_{B_{r}(x_{0})}|\nabla(w-h)|^{p}\,dy,$

where $h$ satisfies $\Delta_{p}h=0$ in $B_{r}(x_{0})$ taking boundary values
equal to $w$ on $\partial B_{r}(x_{0})$.

###### Proof.

Let $\tau\in(0,1)$ be fixed. For $\xi$ with $|\xi|=1$ we set

$t_{\xi}:=\inf\\{t\in[\tau r,r]:w(x_{0}+t\xi)=0\\}$

provided that this set is nonempty. Otherwise we set $t_{\xi}:=r$. Now note
that $w-h$ and $w$ are absolutely continuous in almost every direction $\xi$;
in particular we have $w(x_{0}+t_{\xi}\xi)=0$ (note that this will not be
necessary true if we allow $\tau$ to be zero). Also $w-h$ is
$\mathcal{H}^{n-1}$-a.e. zero on $\partial B_{r}(x_{0})$ as its trace is zero
there, so $(w-h)(x_{0}+r\xi)=0$. Thus for almost every $\xi$ for which
$t_{\xi}<r$ we have

$\displaystyle h(x_{0}+t_{\xi}\xi)$
$\displaystyle=(w-h)(x_{0}+r\xi)-(w-h)(x_{0}+t_{\xi}\xi)$
$\displaystyle=\int_{t_{\xi}}^{r}\frac{d}{dt}\big{(}(w-h)(x_{0}+t\xi)\big{)}\,dt$
$\displaystyle=\int_{t_{\xi}}^{r}\nabla_{\xi}(w-h)(x_{0}+t\xi)\,dt$
$\displaystyle\leq(r-t_{\xi})^{\frac{p-1}{p}}\Big{(}\int_{t_{\xi}}^{r}\big{|}\nabla(w-h)(x_{0}+t\xi)\big{|}^{p}\,dt\Big{)}^{\frac{1}{p}}.$

On the other hand, using Hopf’s lemma we get

$\displaystyle h(x_{0}+t_{\xi}\xi)$ $\displaystyle\geq
c(n,p)\,\mathrm{dist}\big{(}x_{0}+t_{\xi}\xi,\,\partial
B_{r}(x_{0})\big{)}\frac{1}{r}\sup_{B_{r/2}(x_{0})}h$
$\displaystyle=c(n,p)(r-t_{\xi})\frac{1}{r}\sup_{B_{r/2}(x_{0})}h.$

Hence we obtain

$(r-t_{\xi})\Big{(}\frac{1}{r}\sup_{B_{r/2}(x_{0})}h\Big{)}^{p}\leq
C(n,p)\int_{t_{\xi}}^{r}\big{|}\nabla(w-h)(x_{0}+t\xi)\big{|}^{p}\,dt.$

Note that this inequality is trivially satisfied if $t_{\xi}=r$. Now by
integrating over $\xi$ we get

$\displaystyle C(n,p)\int_{B_{r}(x_{0})}\big{|}\nabla(w-h)(x)\big{|}^{p}\,dx$
$\displaystyle\geq C(n,p)\int_{\partial
B_{1}(0)}\int_{t_{\xi}}^{r}\big{|}\nabla(w-h)(x_{0}+t\xi)\big{|}^{p}\,dtd\xi$
$\displaystyle\geq\Big{(}\frac{1}{r}\sup_{B_{r/2}(x_{0})}h\Big{)}^{p}\int_{\partial
B_{1}(0)}(r-t_{\xi})\,d\xi\allowdisplaybreaks$
$\displaystyle=\Big{(}\frac{1}{r}\sup_{B_{r/2}(x_{0})}h\Big{)}^{p}\int_{\partial
B_{1}(0)}\int_{t_{\xi}}^{r}1\,dtd\xi$
$\displaystyle\geq\Big{(}\frac{1}{r}\sup_{B_{r/2}(x_{0})}h\Big{)}^{p}\int_{B_{r}(x_{0})-B_{\tau
r}(x_{0})}\chi_{\\{w=0\\}}\,dx,$

where the last inequality follows from the definition of $t_{\xi}$. Finally,
we get the desired by letting $\tau\to 0$. ∎

######  Lemma 2.6.

Let $\mathbf{u}=\mathbf{u}_{\varepsilon}^{\delta}$ be a minimizer of
$J_{\varepsilon}$ over $V_{\delta}$, and $B\subset\Omega$ be a ball. Then
there exists a unique $v^{i}\in W^{1,p}(\Omega)$ that minimizes the functional

$\int_{\Omega}|\nabla v^{i}|^{p}\,dx$

among all functions with $v^{i}=\varphi^{i}$ on $\partial\Omega$ and
$v^{i}\leq 0$ on $\\{u^{i}=0\\}-B$. The functions $v^{i}$ also satisfy

  1. (1)

$v^{i}=0$ on $\\{u^{i}=0\\}-B$,

  2. (2)

$\mathbf{v}=(v^{1},\dots,v^{m})\in V_{\delta}$,

  3. (3)

$0\leq u^{i}\leq v^{i}\leq C_{0}=\max_{\partial\Omega}|\boldsymbol{\varphi}|$,

  4. (4)

$\int_{\Omega}v^{i}\,\Delta_{p}v^{i}\,dx=0$.

######  Remark.

Instead of a ball $B$, we can also use other open subsets of $\Omega$ in the
above lemma.

###### Proof.

The proof is similar to that of Lemma 3.7 in [12] (see also the proof of
Theorem 3.6 in [2]). ∎

######  Theorem 2.7.

Let $\mathbf{u}=\mathbf{u}_{\varepsilon}^{\delta}$ be a minimizer of
$J_{\varepsilon}$ over $V_{\delta}$. There exists a constant
$M=M_{\varepsilon}$, independent of $\delta$, such that if for some $j$ we
have

$\frac{1}{r}\sup_{B_{r/2}(x)}u^{j}\geq M,$

then $B_{r}(x)\subset\\{|\mathbf{u}|>0\\}$, and $\Delta_{p}u^{i}=0$ in
$B_{r}(x)$ for every $i$.

###### Proof.

Let $\mathbf{v}\in V_{\delta}$ be the function given by Lemma 2.6 for
$B_{r}(x)$. Then we have

$J_{\varepsilon}(\mathbf{u})\leq J_{\varepsilon}(\mathbf{v}).$

Let $\mathbf{h}_{0}$ be the vector-valued function in $\Omega$ satisfying
$\Delta_{p}h_{0}^{i}=0$, and taking the boundary values $\boldsymbol{\varphi}$
on $\partial\Omega$. Since $0\leq u^{i}\leq v^{i}\leq h_{0}^{i}$, for each
$z\in\partial\Omega$ we have

$\partial_{\nu}h_{0}^{i}(z)\leq\partial_{\nu}v^{i}(z)\leq\partial_{\nu}u^{i}(z)\leq
0.$

Then by using the fact that $\mathbf{u},\mathbf{v},\mathbf{h}_{0}$ take the
same boundary values and therefore have equal tangential derivatives on
$\partial\Omega$, we deduce that

$a\leq A_{\nu}h_{0}^{i}(z)\leq A_{\nu}v^{i}(z)\leq A_{\nu}u^{i}(z),$

where $a$ is a lower bound for $A_{\nu}h_{0}^{i}$ (note that $a$ does not
depend on $\delta$).

Hence by the property (2) of $\Gamma$ we have

$\displaystyle\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}\mathbf{u}(x)\big{)}-\Gamma\big{(}x,A_{\nu}\mathbf{v}(x)\big{)}\,d\sigma$
$\displaystyle\qquad=\sum_{i=1}^{m}\int_{\partial\Omega}\Gamma(x,A_{\nu}u^{1},\dots,A_{\nu}u^{i-1},A_{\nu}u^{i},A_{\nu}v^{i+1},\dots,A_{\nu}v^{m})$
(2.3)
$\displaystyle\qquad\qquad\qquad-\Gamma(x,A_{\nu}u^{1},\dots,A_{\nu}u^{i-1},A_{\nu}v^{i},A_{\nu}v^{i+1},\dots,A_{\nu}v^{m})\,d\sigma$
$\displaystyle\qquad\geq
C_{a}\sum_{i=1}^{m}\int_{\partial\Omega}A_{\nu}u^{i}-A_{\nu}v^{i}\,d\sigma,$

where $C_{a}>0$ is the lower bound of $\partial_{\xi_{i}}\Gamma$’s on the set
$\\{(x,\xi):\xi_{i}\geq a\\}$. On the other hand, using the identity (2.2) we
get

$\displaystyle C_{0}\int_{\partial\Omega}A_{\nu}u^{i}-A_{\nu}v^{i}\,d\sigma$
$\displaystyle\geq\int_{\partial\Omega}\varphi^{i}\big{(}A_{\nu}u^{i}-A_{\nu}v^{i}\big{)}\,d\sigma$
$\displaystyle=\int_{\Omega}u^{i}\,\Delta_{p}u^{i}+|\nabla u^{i}|^{p}\,dy$
(2.4) $\displaystyle\qquad\qquad-\int_{\Omega}v^{i}\,\Delta_{p}v^{i}+|\nabla
v^{i}|^{p}\,dy$ $\displaystyle\geq\int_{\Omega}|\nabla
u^{i}|^{p}\,dy-\int_{\Omega}|\nabla v^{i}|^{p}\,dy,$

where $C_{0}=\max_{\partial\Omega}|\boldsymbol{\varphi}|$, and in the last
line we used the facts that $\int_{\Omega}v^{i}\,\Delta_{p}v^{i}\,dy=0$ and
$u^{i},\Delta_{p}u^{i}\geq 0$. Now consider the function $h^{i}$ in $B_{r}(x)$
satisfying $\Delta_{p}h^{i}=0$, and taking boundary values equal to $u^{i}$.
We extend $h^{i}$ to be equal to $u^{i}$ outside of $B_{r}(x)$. Then we have
$\mathbf{h}=(h_{1},\dots,h_{m})\in V_{\delta}$. In addition,
$h^{i}=u^{i}=\varphi^{i}$ on $\partial\Omega$ and $h^{i}=u^{i}=0$ on
$\\{u^{i}=0\\}-B_{r}(x)$. Hence due to the minimality property of $v^{i}$
given by Lemma 2.6 we have

$\int_{\Omega}|\nabla v^{i}|^{p}\,dy\leq\int_{\Omega}|\nabla h^{i}|^{p}\,dy.$

Combining this with the above inequality we get

$\displaystyle C_{0}\int_{\partial\Omega}A_{\nu}u^{i}-A_{\nu}v^{i}\,d\sigma$
$\displaystyle\geq\int_{\Omega}|\nabla u^{i}|^{p}-|\nabla v^{i}|^{p}\,dy$
$\displaystyle\geq\int_{\Omega}|\nabla u^{i}|^{p}-|\nabla h^{i}|^{p}\,dy$
$\displaystyle\geq C\int_{B_{r}(x)}|\nabla(u^{i}-h^{i})|^{p}\,dy,$

where the last inequality can be proved similarly to the proof of Lemma 3.1 in
[7]. (Note that in the last line we have also used the fact that $u^{i}=h^{i}$
outside $B_{r}(x)$.)

Summing the above inequality for each $i$, and using the facts that
$J_{\varepsilon}(\mathbf{u})\leq J_{\varepsilon}(\mathbf{v})$, and
$f_{\varepsilon}$ has Lipschitz constant equal to $\frac{1}{\varepsilon}$ we
get

$\displaystyle\frac{C_{a}}{C_{0}}\sum_{i\leq
m}\int_{B_{r}(x)}|\nabla(u^{i}-h^{i})|^{p}\,dy$ $\displaystyle\leq
C_{a}\int_{\partial\Omega}\sum_{i\leq m}(A_{\nu}u^{i}-A_{\nu}v^{i})\,d\sigma$
$\displaystyle\leq\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}\mathbf{u}(x)\big{)}-\Gamma\big{(}x,A_{\nu}\mathbf{v}(x)\big{)}\,d\sigma$
$\displaystyle\leq
f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{v}|>0\\}\big{|}\big{)}-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}$
$\displaystyle\leq\frac{1}{\varepsilon}\big{|}B_{r}(x)\cap\\{|\mathbf{u}|=0\\}\big{|},$

since $0\leq u^{i}\leq v^{i}$, and outside of $B_{r}(x)$, $|\mathbf{u}|=0$
implies $|\mathbf{v}|=0$. Therefore by Lemma 2.5 applied to $u^{j}$ we obtain

$\displaystyle\big{|}B_{r}(x)\cap\\{|\mathbf{u}|=0\\}\big{|}$
$\displaystyle\geq\frac{\varepsilon C_{a}}{C_{0}}\sum_{i\leq
m}\int_{B_{r}(x)}|\nabla(u^{i}-h^{i})|^{p}\,dy$
$\displaystyle\geq\frac{\varepsilon
C_{a}}{C_{0}}\int_{B_{r}(x)}|\nabla(u^{j}-h^{j})|^{p}\,dy$
$\displaystyle\geq\frac{\varepsilon
C_{a}}{cC_{0}}\Big{(}\frac{1}{r}\sup_{B_{r/2}(x)}h^{j}\Big{)}^{p}\cdot\big{|}B_{r}(x)\cap\\{u^{j}=0\\}\big{|}$
$\displaystyle\geq\frac{\varepsilon
C_{a}}{cC_{0}}\Big{(}\frac{1}{r}\sup_{B_{r/2}(x)}u^{j}\Big{)}^{p}\cdot\big{|}B_{r}(x)\cap\\{u^{j}=0\\}\big{|}$
$\displaystyle\geq\frac{\varepsilon
C_{a}M^{p}}{cC_{0}}\big{|}B_{r}(x)\cap\\{|\mathbf{u}|=0\\}\big{|},$

since $|\mathbf{u}|=0$ implies $u^{j}=0$, and $h^{j}\geq u^{j}$ as $u^{j}$ is
$p$-subharmonic. Hence if $M>(\frac{cC_{0}}{\varepsilon C_{a}})^{\frac{1}{p}}$
then $\big{|}B_{r}(x)\cap\\{|\mathbf{u}|=0\\}\big{|}$ must be zero, as
desired. Note that in this case the above inequality also implies that
$u^{i}=h^{i}$ in $B_{r}(x)$ for each $i$; so $u^{i}$ satisfies the equation in
$B_{r}(x)$. ∎

######  Corollary 2.8.

All minimizers $\mathbf{u}_{\varepsilon}^{\delta}$ are Lipschitz, and for
every $\Omega^{\prime}\subset\subset\Omega$ there exists a constant
$K_{\varepsilon}=K_{\varepsilon}(\Omega^{\prime})$, independent of $\delta$,
such that

$\|\mathbf{u}_{\varepsilon}^{\delta}|_{\Omega^{\prime}}\|_{\mathrm{Lip}}\leq
K_{\varepsilon}.$

In addition, $\Delta_{p}(u_{\varepsilon}^{\delta})^{i}=0$ in the open set
$\\{|\mathbf{u}_{\varepsilon}^{\delta}|>0\\}$.

###### Proof.

For simplicity we set $\mathbf{u}=\mathbf{u}_{\varepsilon}^{\delta}$. First
let us show that $\\{|\mathbf{u}|>0\\}$ is an open set. Suppose
$x\in\\{|\mathbf{u}|>0\\}$. Then $u^{j}(x)>0$ for some $j$. Then for small
enough $r$ we must have

$\frac{1}{r}\sup_{B_{r/2}(x)}u^{j}\geq\frac{1}{r}u^{j}(x)\geq M.$

Hence the previous theorem implies that $B_{r}(x)\subset\\{|\mathbf{u}|>0\\}$
and we have $\Delta_{p}u^{i}=0$ in $B_{r}(x)$.

Next note that $\nabla\mathbf{u}=0$ a.e. in $\\{|\mathbf{u}|=0\\}$. So suppose
$x\in\\{|\mathbf{u}|>0\\}\cap\Omega^{\prime}$. Let
$\Omega^{\prime}\subset\subset\tilde{\Omega}\subset\subset\Omega$, and
$B=B_{d}(x)$, where
$d=\mathrm{dist}\big{(}x,\partial(\\{|\mathbf{u}|>0\\}\cap\tilde{\Omega})\big{)}$.
If $\partial B$ touches $\partial\\{|\mathbf{u}|=0\\}$ then
$B_{d+d^{\prime}}(x)$ intersects $\\{|\mathbf{u}|=0\\}$, and by previous
theorem we have

$\frac{1}{d+d^{\prime}}\sup_{B_{(d+d^{\prime})/2}(x)}u^{i}\leq M$

for every $i$. Hence in the limit $d^{\prime}\to 0$ we get

$\frac{1}{d}\sup_{B_{d/2}(x)}u^{i}\leq M.$

Now since $u^{i}$’s are $p$-harmonic in $B$, as shown in the proof of Lemma
3.1 of [6], we have

$|\nabla u^{i}(x)|\leq C\frac{1}{d}\sup_{B_{d/2}(x)}u^{i}\leq CM,$

where the constant $C$ depends only on $p$ and the dimension $n$. On the other
hand, if $\partial B$ touches $\partial\tilde{\Omega}$ then, by the interior
derivative estimate of [9], we obtain (the dependence on $d$ follows from the
proof of this estimate; see equation (3.4) in [9])

$|\nabla u^{i}(x)|\leq\frac{C(n,p)}{d^{n}}\|\mathbf{u}\|_{W^{1,p}}\leq C,$

since $d\geq\mathrm{dist}(\Omega^{\prime},\partial\tilde{\Omega})$, and
$\|\mathbf{u}\|_{W^{1,p}}$ is bounded independently of $\delta$ as will be
shown now. Let $\Omega^{\prime}\subset\subset\Omega$ be a smooth open set with
$|\Omega-\Omega^{\prime}|=1$. Let $\mathbf{u}_{0}$ be a vector-valued function
on $\Omega-\Omega^{\prime}$ that satisfies the equation
$\Delta_{p}u_{0}^{i}=0$, and takes the boundary values $\boldsymbol{\varphi}$
on $\partial\Omega$ and $0$ on $\partial\Omega^{\prime}$. Then for every small
enough $\delta$ we have $\mathbf{u}_{0}\in V_{\delta}$. Hence (remember that
$\mathbf{u}=\mathbf{u}_{\varepsilon}^{\delta}$)

$\displaystyle C=J_{\varepsilon}(\mathbf{u}_{0})\geq
J_{\varepsilon}(\mathbf{u}_{\varepsilon}^{\delta})$
$\displaystyle\geq\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}\mathbf{u}_{\varepsilon}^{\delta}(x)\big{)}\,d\sigma$
$\displaystyle\geq\int_{\partial\Omega}\sum_{i=1}^{m}\psi_{i}(x)A_{\nu}(u_{\varepsilon}^{\delta})^{i}-C\,d\sigma,$

where we used (1.2) in the last line. Thus by Lemma 2.2 the
$\|\nabla\mathbf{u}_{\varepsilon}^{\delta}\|_{L^{p}(\Omega;\mathbb{R}^{m})}$
is bounded as $\delta\to 0$, and the boundedness of
$\|\mathbf{u}_{\varepsilon}^{\delta}\|_{W^{1,p}}$ follows from Poincare
inequality and the fact that all of $\mathbf{u}_{\varepsilon}^{\delta}$’s have
the same boundary values.

Finally, to see that $\mathbf{u}$ is Lipschitz continuous on all of $\Omega$,
note that $\mathbf{u}$ has $p$-harmonic components near the smooth boundary
$\partial\Omega$, attaining smooth boundary conditions $\boldsymbol{\varphi}$;
hence the gradient of $\mathbf{u}$ is bounded near the boundary too. ∎

######  Lemma 2.9.

There exists $\delta_{0}=\delta_{0}(\varepsilon)>0$, such that for every
$\delta$ we have $|\mathbf{u}_{\varepsilon}^{\delta}|>0$ on
$\Omega_{\delta_{0}}$.

######  Remark.

Note that as a consequence, $\Delta_{p}(u_{\varepsilon}^{\delta})^{i}=0$ on
$\Omega_{\delta_{0}}$ for every $\delta$ (by Theorem 2.7). In other words,
$\mathbf{u}_{\varepsilon}^{\delta}\in V_{\delta_{0}}$ for every $\delta$.

###### Proof.

Suppose to the contrary that there is a sequence
$\mathbf{u}_{k}=\mathbf{u}_{\varepsilon}^{\delta_{k}}$ for which we have

$2d_{k}:=\mathrm{dist}(\\{|\mathbf{u}_{k}|=0\\},\partial\Omega)\to 0.$

Then the midpoint of the closest points on $\\{|\mathbf{u}_{k}|=0\\}$ and
$\partial\Omega$, which we call $x_{k}$, has distance $d_{k}$ from both of
these sets. So the boundary of the ball $B_{d_{k}}(x_{k})$ touches both of
these sets. In addition, by Theorem 2.7, for every $t>0$ we must have

$\frac{1}{d_{k}}\sup_{B_{d_{k}/2}(x_{k}+t\nu_{k})}u_{k}^{i}\leq
M_{\varepsilon}$

for every $i$ (here $\nu_{k}$ is the direction of the line segment from
$x_{k}$ to its closest point on $\\{|\mathbf{u}_{k}|=0\\}$). So in the limit
$t\to 0$ we get

(2.5) $\sup_{B_{d_{k}/2}(x_{k})}u_{k}^{i}\leq M_{\varepsilon}d_{k}.$

We also have

$\sup_{B_{d_{k}}(x_{k})}|\mathbf{u}_{k}|\geq c_{0},$

where $c_{0}=\min_{i}\min_{\partial\Omega}\varphi^{i}>0$. Because at the point
$y_{k}\in\partial B_{d_{k}}(x_{k})\cap\partial\Omega$ we have
$u_{k}^{i}(y_{k})=\varphi^{i}(y_{k})\geq c_{0}$ (note that $u_{k}^{i}$ is
continuous up to the boundary).

Next consider the functions

$\hat{\mathbf{u}}_{k}(x):=\frac{\mathbf{u}(x_{k}+d_{k}x)}{\sup_{B_{d_{k}}(x_{k})}|\mathbf{u}_{k}|}$

on $B_{1}$. Then $\hat{u}_{k}^{i}$ is positive and $p$-harmonic on $B_{1}$,
and we have $\sup_{B_{1}}|\hat{\mathbf{u}}_{k}|=1$. In addition, by (2.5) we
have

$\sup_{B_{1/2}}\hat{u}_{k}^{i}=\frac{\sup_{B_{d_{k}/2}(x_{k})}u_{k}^{i}}{\sup_{B_{d_{k}}(x_{k})}|\mathbf{u}_{k}|}\leq\frac{M_{\varepsilon}d_{k}}{c_{0}}\underset{k\to\infty}{\longrightarrow}0.$

Furthermore, note that $\hat{u}_{k}^{i}$ is a uniformly bounded sequence of
$p$-harmonic functions on $B_{1}$, so there is $\alpha>0$ such that for all
$r<1$ the Hölder norms $\|\hat{u}_{k}^{i}\|_{C^{0,\alpha}(\overline{B}_{r})}$
are uniformly bounded (see page 251 of [8]). Hence, by a diagonal argument, we
can construct a subsequence of $\hat{u}_{k}^{i}$, which we still denote by
$\hat{u}_{k}^{i}$, that locally uniformly converges to a nonnegative
$p$-harmonic function $\hat{u}_{\infty}^{i}$ on $B_{1}$. In addition,
$\hat{u}_{\infty}^{i}$ must vanish on $B_{1/2}$ by the above estimate. Thus by
the strong maximum principle we must have $\hat{\mathbf{u}}_{\infty}\equiv 0$
on $B_{1}$.

Now for $y_{k}\in\partial B_{d_{k}}(x_{k})\cap\partial\Omega$ and $r<d_{k}$ we
have

$\underset{B_{r}(y_{k})\cap\Omega}{\mathrm{osc}}\,u_{k}^{i}\leq
C(n,p)\Big{(}r^{\alpha}+\underset{B_{r}(y_{k})\cap\partial\Omega}{\mathrm{osc}}\,\varphi^{i}\Big{)}\leq
Cr^{\alpha}$

for some $\alpha\in(0,1)$. This estimate holds by Theorem 4.19 of [11] when
$1<p\leq n$. And when $p>n$ this estimate holds due to the uniform Hölder
continuity of $u_{k}^{i}$ on $\overline{\Omega}$, since
$\|\mathbf{u}_{k}\|_{W^{1,p}(\Omega)}$ is uniformly bounded as we have seen in
the proof of Corollary 2.8. Hence for $r=d_{k}/2$ we have

$\min_{B_{d_{k}/2}(y_{k})\cap
B_{d_{k}}(x_{k})}u_{k}^{i}\geq\min_{B_{d_{k}/2}(y_{k})\cap\Omega}u_{k}^{i}\geq\frac{1}{2}c_{0},$

where $c_{0}=\min_{i}\min_{\partial\Omega}\varphi^{i}$. Therefore for
$\hat{y}_{k}=\frac{1}{d_{k}}(y_{k}-x_{k})\in\partial B_{1}$ we have

$\min_{B_{1/2}(\hat{y}_{k})\cap
B_{1}}\hat{u}_{k}^{i}=\frac{1}{\sup_{B_{d_{k}}(x_{k})}|\mathbf{u}_{k}|}\min_{B_{d_{k}/2}(y_{k})\cap
B_{d_{k}}(x_{k})}u_{k}^{i}\geq c>0,$

since $\sup_{B_{d_{k}}(x_{k})}|\mathbf{u}_{k}|\leq mC_{0}$ where
$C_{0}=\max_{\partial\Omega}|\boldsymbol{\varphi}|$. Thus $\hat{u}_{k}^{i}$
has a uniform positive lower bound on a subset of $B_{1}$ with positive volume
(where the volume is independent of $k$). So no subsequence of
$\hat{\mathbf{u}}_{k}$ can converge locally uniformly to
$\hat{\mathbf{u}}_{\infty}\equiv 0$, because otherwise they will uniformly
converge to $0$ outside a set of small volume, contradicting the uniform
boundedness from below. ∎

Now we can find a minimizer for $J_{\varepsilon}$ over $V$.

######  Theorem 2.10.

There exists a minimizer $\mathbf{u}_{\varepsilon}\in V$ for
$J_{\varepsilon}$. Moreover, $\mathbf{u}_{\varepsilon}$ is a Lipschitz
function, and $\Delta_{p}u_{\varepsilon}^{i}=0$ in the open set
$\\{|\mathbf{u}_{\varepsilon}|>0\\}$.

######  Remark.

As we will see in the following proof, $\mathbf{u}_{\varepsilon}^{\delta}\in
V_{\delta_{0}}$ for $\delta_{0}=\delta_{0}(\varepsilon)$ given by the above
lemma. So in fact $\mathbf{u}_{\varepsilon}$ is a minimizer of
$J_{\varepsilon}$ over some $V_{\delta}$, and therefore it has all the
properties of $\mathbf{u}_{\varepsilon}^{\delta}$’s that we have proved so
far. In particular, we have $|\mathbf{u}_{\varepsilon}|>0$ on
$\Omega_{\delta_{0}}$.

###### Proof.

As we have shown in the proof of Corollary 2.8,
$\|\nabla\mathbf{u}_{\varepsilon}^{\delta}\|_{L^{p}(\Omega;\mathbb{R}^{m})}$
is bounded as $\delta\to 0$. Hence there is a subsequence such that
$\mathbf{u}_{\varepsilon}^{\delta}\rightharpoonup\mathbf{u}_{\varepsilon}$
weakly in $W^{1,p}$ (and also a.e.) with
$A(\nabla(u_{\varepsilon}^{\delta})^{i})\rightharpoonup A(\nabla
u_{\varepsilon}^{i})$ in $L^{q}$ as $\delta\to 0$. So, in particular,
$u_{\varepsilon}^{i}\geq 0$, $u_{\varepsilon}^{i}$ is $p$-subharmonic, and
attains the boundary condition $\varphi^{i}$. Furthermore, by Corollary 2.8,
$\mathbf{u}_{\varepsilon}^{\delta}\to\mathbf{u}_{\varepsilon}$ uniformly on
compact subsets of $\Omega$. Hence for each ball
$\overline{B}\subset\\{|\mathbf{u}_{\varepsilon}|>0\\}$ and all small enough
$\delta$ we have
$\overline{B}\subset\\{|\mathbf{u}_{\varepsilon}^{\delta}|>0\\}$. Therefore by
using test functions with support in $B$ together with
$A(\nabla(u_{\varepsilon}^{\delta})^{i})\rightharpoonup A(\nabla
u_{\varepsilon}^{i})$ we can conclude that $u_{\varepsilon}^{i}$ is
$p$-harmonic in $B$.

The same reasoning applied to test functions with support in
$\Omega_{\delta_{0}}$, for $\delta_{0}$ given by the previous lemma, implies
that $u_{\varepsilon}^{i}$ is $p$-harmonic in $\Omega_{\delta_{0}}$, and thus
$\mathbf{u}_{\varepsilon}\in V_{\delta_{0}}\subset V$. In particular,
$u_{\varepsilon}^{i}$ is $p$-harmonic near the smooth boundary
$\partial\Omega$, attaining smooth boundary conditions $\varphi^{i}$, so it is
Lipschitz near $\partial\Omega$. Moreover, $\mathbf{u}_{\varepsilon}$ is
Lipschitz inside $\Omega$ away from its boundary, because it is the uniform
limit of a sequence of Lipschitz functions with uniformly bounded Lipschitz
constants. Hence $\mathbf{u}_{\varepsilon}$ is Lipschitz on all of $\Omega$.

Finally note that $\mathbf{u}_{\varepsilon}$ minimizes $J_{\varepsilon}$ over
$V$, since for every $\mathbf{w}\in V$ we have $\mathbf{w}\in V_{\delta}$ for
some $\delta$. Thus $J_{\varepsilon}(\mathbf{u}_{\varepsilon}^{\delta})\leq
J_{\varepsilon}(\mathbf{w})$. However,
$\mathbf{u}_{\varepsilon}^{\delta}\to\mathbf{u}_{\varepsilon}$, so we get
$J_{\varepsilon}(\mathbf{u}_{\varepsilon})\leq J_{\varepsilon}(\mathbf{w})$
due to the semicontinuity of $J_{\varepsilon}$. ∎

##  3\. Regularity of solutions to the penalized problem

To simplify the notation, throughout this section we will suppress the index
$\varepsilon$ in $\mathbf{u}_{\varepsilon}$.

######  Theorem 3.1.

For $\tau\in(0,1/4)$ there exists $m_{\varepsilon}(\tau)$ such that if for
each $i$ we have

$\frac{1}{r}\sup_{B_{r/2}(x)}u^{i}\leq m_{\varepsilon}(\tau),$

then $B_{\tau r}(x)\subset\\{|\mathbf{u}|=0\\}$.

###### Proof.

Similarly to Lemma 2.6, we can show that there is $v^{i}\in W^{1,p}(\Omega)$
that minimizes the functional $\int_{\Omega}|\nabla v^{i}|^{p}\,dx$ among all
functions with $v^{i}=\varphi^{i}$ on $\partial\Omega$ and $v^{i}\leq 0$ on
$\\{u^{i}=0\\}\cup\overline{B}_{\tau r}(x)$. The function $v^{i}$ also
satisfies $\Delta_{p}v^{i}\geq 0$,
$\int_{\Omega}v^{i}\,\Delta_{p}v^{i}\,dx=0$, and $u^{i}\geq v^{i}\geq 0$ (to
see this, note that $\Delta_{p}v^{i}\geq\Delta_{p}u^{i}$ on
$\Omega-\big{(}\\{u^{i}=0\\}\cup\overline{B}_{\tau
r}(x)\big{)}\subset\\{|\mathbf{u}|>0\\}$, and $v^{i}-u^{i}\leq 0$ on
$\\{u^{i}=0\\}\cup\overline{B}_{\tau r}(x)$ or $\partial\Omega$). In addition,
we have $\mathbf{v}=(v_{1},\dots,v_{m})\in V_{\delta_{1}}\subset V$ (where
$\delta_{1}$ is small enough so that $\overline{B}_{\tau
r}(x)\subset\Omega-\Omega_{\delta_{1}}$). Thus

$J_{\varepsilon}(\mathbf{u})\leq J_{\varepsilon}(\mathbf{v}).$

Let us assume that $\delta_{1}$ is small enough so that
$\overline{B}_{r}(x)\subset\Omega-\Omega_{\delta_{1}}$ and $\mathbf{u}\in
V_{\delta_{1}}$. Let $\mathbf{w}$ be a vector-valued $p$-harmonic function in
$\Omega_{\delta_{1}}$ with boundary values equal to $\boldsymbol{\varphi}$ on
$\partial\Omega$ and equal to 0 on
$\partial\Omega_{\delta_{1}}-\partial\Omega$. Then we have $u^{i}\geq
v^{i}\geq w^{i}\geq 0$ (since $\mathbf{u},\mathbf{v}$ are also $p$-harmonic on
$\Omega_{\delta_{1}}$, and nonnegative everywhere). Thus for each
$z\in\partial\Omega$ we have

$0\geq\partial_{\nu}w^{i}(z)\geq\partial_{\nu}v^{i}(z)\geq\partial_{\nu}u^{i}(z).$

Next using the fact that $\mathbf{u},\mathbf{v},\mathbf{w}$ take the same
boundary values on $\partial\Omega$, and therefore have equal tangential
derivatives on $\partial\Omega$, we deduce that

$0\geq A_{\nu}w^{i}(z)\geq A_{\nu}v^{i}(z)\geq A_{\nu}u^{i}(z).$

Now similar to (2) we can show that

$\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}\mathbf{v}(x)\big{)}-\Gamma\big{(}x,A_{\nu}\mathbf{u}(x)\big{)}\,d\sigma\leq
C_{1}\sum_{i=1}^{m}\int_{\partial\Omega}A_{\nu}v^{i}-A_{\nu}u^{i}\,d\sigma,$

where $C_{1}>0$ is the upper bound of $\partial_{\xi_{i}}\Gamma$’s on the set
$\\{(x,\xi):\xi_{i}\leq 0\\}$. On the other hand, using the identity (2.2) we
obtain (using the notation $c_{0}=\min_{i}\min_{\partial\Omega}\varphi^{i}$)

$\displaystyle c_{0}\int_{\partial\Omega}A_{\nu}v^{i}-A_{\nu}u^{i}\,d\sigma$
$\displaystyle\leq\int_{\partial\Omega}\varphi^{i}\big{(}A_{\nu}v^{i}-A_{\nu}u^{i}\big{)}\,d\sigma$
$\displaystyle=\int_{\Omega}v^{i}\,\Delta_{p}v^{i}+|\nabla v^{i}|^{p}\,dy$
$\displaystyle\qquad\qquad-\int_{\Omega}u^{i}\,\Delta_{p}u^{i}+|\nabla
u^{i}|^{p}\,dy$ (3.1) $\displaystyle=\int_{\Omega}|\nabla
v^{i}|^{p}\,dy-\int_{\Omega}|\nabla u^{i}|^{p}\,dy,$

where in the last line we used the facts that
$\int_{\Omega}v^{i}\,\Delta_{p}v^{i}\,dy=0$, and $\Delta_{p}u^{i}=0$ on
$\\{u^{i}\neq 0\\}\subset\\{|\mathbf{u}|>0\\}$.

Summing the above inequality for each $i$, and using the facts that
$J_{\varepsilon}(\mathbf{u})\leq J_{\varepsilon}(\mathbf{v})$, and the
derivative of $f_{\varepsilon}$ is bounded below by $\varepsilon$, we get

$\displaystyle\frac{C_{1}}{c_{0}}\sum_{i\leq m}\int_{\Omega}|\nabla
v^{i}|^{p}-|\nabla u^{i}|^{p}\,dy$ $\displaystyle\geq
C_{1}\int_{\partial\Omega}\sum_{i\leq m}(A_{\nu}v^{i}-A_{\nu}u^{i})\,d\sigma$
$\displaystyle\geq\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}\mathbf{v}(x)\big{)}-\Gamma\big{(}x,A_{\nu}\mathbf{u}(x)\big{)}\,d\sigma$
$\displaystyle\geq
f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{v}|>0\\}\big{|}\big{)}$
$\displaystyle\geq\varepsilon\big{|}\\{|\mathbf{u}|>0\\}\cap\\{|\mathbf{v}|=0\\}\big{|}$
(3.2) $\displaystyle\geq\varepsilon\big{|}\\{|\mathbf{u}|>0\\}\cap B_{\tau
r}(x)\big{|},$

since $u^{i}\geq v^{i}\geq 0$ and $v^{i}=0$ in $B_{\tau r}(x)$.

Next we define $g:(0,\infty)\to\mathbb{R}$ by

$g(t):=\begin{cases}t^{\frac{p-n}{p-1}}-(\tau r)^{\frac{p-n}{p-1}}&p>n,\\\
\log t-\log(\tau r)&p=n,\\\ (\tau
r)^{\frac{p-n}{p-1}}-t^{\frac{p-n}{p-1}}&p<n.\end{cases}$

Note that $g$ is an increasing function that vanishes at $t=\tau r$, and is
negative for $t<\tau r$. In addition, $g(|x|)$ is a $p$-harmonic function in
$\mathbb{R}^{n}-\\{0\\}$, which is negative on $B_{\tau r}(x)$ and vanishes on
$\partial B_{\tau r}(x)$. Now let us define
$h^{i}:B_{\sqrt{\tau}r}(x)\to\mathbb{R}$ by

$h^{i}(y):=\min\\{u^{i}(y),\,\frac{s_{i}}{g(\sqrt{\tau}r)}\big{(}g(|y-x|)\big{)}^{+}\\},$

where $s_{i}:=\max_{\overline{B}_{\sqrt{\tau}r}(x)}u^{i}$. We extend $h^{i}$
by $u^{i}$ outside of $B_{\sqrt{\tau}r}(x)$. Note that we have $h^{i}=0$ on
$\\{u^{i}=0\\}\cap\overline{B}_{\tau r}(x)$ and $h^{i}=u^{i}=\varphi^{i}$ on
$\partial\Omega$. Hence $h^{i}$ competes with $v^{i}$, and we have
$\int_{\Omega}|\nabla v^{i}|^{p}\,dx\leq\int_{\Omega}|\nabla h^{i}|^{p}\,dx$.
Therefore we can exchange $v^{i}$ by $h^{i}$ in the inequality (3.2) to get

$\frac{\varepsilon c_{0}}{C_{1}}\big{|}\\{|\mathbf{u}|>0\\}\cap B_{\tau
r}(x)\big{|}\leq\sum_{i\leq m}\int_{B_{\sqrt{\tau}r}(x)}|\nabla
h^{i}|^{p}-|\nabla u^{i}|^{p}\,dy.$

Now since $h^{i}=0$ on $B_{\tau r}(x)$, we can rewrite the above inequality as

$\displaystyle\frac{\varepsilon c_{0}}{C_{1}}\big{|}\\{|\mathbf{u}|>0\\}\cap
B_{\tau r}(x)\big{|}\,+$ $\displaystyle\sum_{i\leq m}\int_{B_{\tau
r}(x)}|\nabla u^{i}|^{p}\,dy$ (3.3) $\displaystyle\leq\sum_{i\leq
m}\int_{B_{\sqrt{\tau}r}(x)-B_{\tau r}(x)}|\nabla h^{i}|^{p}-|\nabla
u^{i}|^{p}\,dy.$

But

$|\nabla h^{i}|^{p}-|\nabla u^{i}|^{p}\leq-p|\nabla h^{i}|^{p-2}\nabla
h^{i}\cdot\nabla(u^{i}-h^{i}),$

since for two vectors $a,b$ we have
$|a|^{p}-|b|^{p}\leq-p|a|^{p-2}a\cdot(b-a)$. So we can estimate the right hand
side of (3.3) as follows (using integration by parts, and the facts that
$\Delta_{p}h^{i}=0$ on $\\{u^{i}>h^{i}\\}$, $h^{i}=0$ on $\partial B_{\tau
r}(x)$, and $h^{i}=u^{i}$ on $\partial B_{\sqrt{\tau}r}(x)$):

$\displaystyle\int_{B_{\sqrt{\tau}r}(x)-B_{\tau r}(x)}$ $\displaystyle|\nabla
h^{i}|^{p}-|\nabla u^{i}|^{p}\,dy$
$\displaystyle\leq-p\int_{B_{\sqrt{\tau}r}(x)-B_{\tau r}(x)}|\nabla
h^{i}|^{p-2}\nabla(u^{i}-h^{i})\cdot\nabla h^{i}\,dy$
$\displaystyle=p\int_{\partial B_{\tau r}(x)}(u^{i}-h^{i})|\nabla
h^{i}|^{p-2}\nabla h^{i}\cdot\nu\,d\sigma$
$\displaystyle\qquad-p\int_{\partial B_{\sqrt{\tau}r}(x)}(u^{i}-h^{i})|\nabla
h^{i}|^{p-2}\nabla h^{i}\cdot\nu\,d\sigma$ $\displaystyle=p\int_{\partial
B_{\tau r}(x)}u^{i}|\nabla h^{i}|^{p-2}\nabla
h^{i}\cdot\nu\,d\sigma=C(n,p,\tau)\frac{s_{i}^{p-1}}{r^{p-1}}\int_{\partial
B_{\tau r}(x)}u^{i}\,d\sigma,$

where the last equality is calculated using the fact
$h^{i}(y)=\frac{s_{i}}{g(\sqrt{\tau}r)}\big{(}g(|y-x|)\big{)}^{+}=0$ on
$\overline{B}_{\tau r}(x)$; hence on $\partial B_{\tau r}(x)$ we have

$\nabla
h^{i}=C(\tau)s_{i}r^{\frac{n-p}{p-1}}\begin{cases}\frac{2|p-n|}{p-1}|y-x|^{\frac{2-n-p}{p-1}}(y-x)&p\neq
n,\\\ |y-x|^{-2}(y-x)&p=n,\end{cases}$

and thus

$\displaystyle|\nabla h^{i}|^{p-2}\nabla h^{i}\cdot\nu$
$\displaystyle=C(n,p,\tau)s_{i}^{p-1}r^{n-p}|y-x|^{2-n-p}|y-x|^{p-2}(y-x)\cdot\frac{(y-x)}{\tau
r}$ $\displaystyle=C(n,p,\tau)s_{i}^{p-1}r^{n-p}|y-x|^{-n+2}\frac{1}{\tau
r}=C(n,p,\tau)\frac{s_{i}^{p-1}}{r^{p-1}}.$

Hence (3.3) becomes

$\displaystyle\frac{\varepsilon c_{0}}{C_{1}}\big{|}\\{|\mathbf{u}|>0\\}\cap
B_{\tau r}(x)\big{|}\,+$ $\displaystyle\sum_{i\leq m}\int_{B_{\tau
r}(x)}|\nabla u^{i}|^{p}\,dy$ (3.4) $\displaystyle\leq C(n,p,\tau)\sum_{i\leq
m}\frac{s_{i}^{p-1}}{r^{p-1}}\int_{\partial B_{\tau r}(x)}u^{i}\,d\sigma.$

On the other hand we have

$\displaystyle\int_{\partial B_{\tau r}(x)}u^{i}\,d\sigma$ $\displaystyle\leq
c(n,\tau)\Big{(}\int_{B_{\tau r}(x)}u^{i}\,dy+\int_{B_{\tau r}(x)}|\nabla
u^{i}|\,dy\Big{)}$ (3.5) $\displaystyle\leq
c(n,\tau)\Big{(}(s_{i}+1)\cdot\big{|}\\{|\mathbf{u}|>0\\}\cap B_{\tau
r}(x)\big{|}+\int_{B_{\tau r}(x)}|\nabla u^{i}|^{p}\,dy\Big{)},$

where in the last line we estimated $u^{i},|\nabla u^{i}|$ from above by
$s_{i},1+|\nabla u^{i}|^{p}$ on the set
$\\{u^{i}>0\\}\subset\\{|\mathbf{u}|>0\\}$. Next note that

(3.6)
$s_{i}=\max_{\overline{B}_{\sqrt{\tau}r}(x)}u^{i}\leq\sup_{B_{r/2}(x)}u^{i}\leq
rm_{\varepsilon}(\tau),$

since $\sqrt{\tau}<1/2$. Combining the inequalities (3.4), (3.5), and (3.6) we
get

$\displaystyle\frac{\varepsilon c_{0}}{C_{1}}\big{|}\\{|\mathbf{u}|>0\\}\cap
B_{\tau r}(x)\big{|}+\sum_{i\leq m}\int_{B_{\tau r}(x)}|\nabla u^{i}|^{p}\,dy$
$\displaystyle\leq cC\sum_{i\leq
m}\frac{s_{i}^{p-1}}{r^{p-1}}\Big{(}(s_{i}+1)\cdot\big{|}\\{|\mathbf{u}|>0\\}\cap
B_{\tau r}(x)\big{|}+\int_{B_{\tau r}(x)}|\nabla u^{i}|^{p}\,dy\Big{)}$
$\displaystyle\leq
cC\,m_{\varepsilon}^{p-1}(\tau)\Big{(}\big{|}\\{|\mathbf{u}|>0\\}\cap B_{\tau
r}(x)\big{|}\sum_{i\leq m}(s_{i}+1)+\sum_{i\leq m}\int_{B_{\tau r}(x)}|\nabla
u^{i}|^{p}\,dy\Big{)}.$

Now if $m_{\varepsilon}(\tau)$ is small enough, we must necessarily have
$|\mathbf{u}|=0$ on $B_{\tau r}(x)$, as desired. ∎

Now let us set

$\displaystyle U$ $\displaystyle:=\\{x\in\Omega:|\mathbf{u}(x)|>0\\},$
$\displaystyle E$ $\displaystyle:=\\{x\in\Omega:|\mathbf{u}(x)|=0\\}.$

######  Lemma 3.2.

For every $i$ we have

$U=\\{x\in\Omega:u^{i}(x)>0\\},\qquad E=\\{x\in\Omega:u^{i}(x)=0\\}.$

###### Proof.

By Theorem 2.10, each $u^{i}$ is $p$-harmonic in the open set $U$. So in each
component of $U$ either $u^{i}>0$ or $u^{i}\equiv 0$ (by the strong maximum
principle). Now consider a component of $U$, say $U_{1}$. If $\partial U_{1}$
does not intersect $\partial\Omega$ then it must be a subset of $E$. Therefore
every $u^{i}$ vanishes on $\partial U_{1}$, and hence every $u^{i}$ vanishes
on $U_{1}$ by the maximum principle. So we would have $U_{1}\subset E$, which
is a contradiction. Thus $\partial U_{1}$ must intersect $\partial\Omega$.
Hence each $u^{i}>0$ on $U_{1}$, since they are positive on $\partial\Omega$.
Therefore each $u^{i}$ is positive on every component of $U$, as desired. ∎

######  Corollary 3.3.

There are $c,C>0$ such that for $x\in U$ near $\partial E$ we have

$c\cdot\mathrm{dist}(x,\partial E)\leq|\mathbf{u}(x)|\leq
C\cdot\mathrm{dist}(x,\partial E).$

###### Proof.

The right hand side inequality holds according to the Lipschitz regularity of
the solutions, Theorem 2.10. To see the left hand side inequality, we argue
indirectly. Assume to the contrary that there exists a sequence $x_{k}\in U$
such that

(3.7) $|\mathbf{u}(x_{k})|\leq\frac{1}{k}\mathrm{dist}(x_{k},\partial E).$

Let $r_{k}=\mathrm{dist}(x_{k},\partial E)$ and define

$\mathbf{u}_{k}(x)=\frac{\mathbf{u}(x_{k}+r_{k}x)}{r_{k}}.$

The sequence $\mathbf{u}_{k}$ is uniformly bounded and uniformly Lipschitz in
$B_{1}$ due to Lipschitz regularity of $\mathbf{u}$ and assumption (3.7).

Recall that $\Delta_{p}u_{k}^{i}=0$ in $U$, then we may choose a converging
subsequence $\mathbf{u}_{k}\to\mathbf{u}_{0}$ such that $u_{0}^{i}$ is also
$p$-harmonic. Furthermore, by Theorem 3.1 we get that

$\sup_{B_{1/2}(0)}|\mathbf{u}_{0}|=\lim_{k\to\infty}\sup_{B_{1/2}(0)}|\mathbf{u}_{k}|\geq
m_{\varepsilon}>0,$

since $|\mathbf{u}_{k}(0)|>0$. Also, (3.7) yields that $\mathbf{u}_{0}(0)=0$,
which contradicts the maximum (minimum) principle; remember that each
component of $\mathbf{u}_{0}$ is nonnegative. ∎

######  Corollary 3.4.

There exists $c=c_{\varepsilon}\in(0,1)$ such that for any $x\in\partial U$
and small enough $r$ we have

(3.8) $c\leq\frac{\big{|}E\cap B_{r}(x)\big{|}}{\big{|}B_{r}(x)\big{|}}\leq
1-c.$

###### Proof.

The proof is similar to the proof of Theorem 4.2 in [7]. By Theorem 3.1, there
exists $z\in B_{r/2}(x)$ such that $|\mathbf{u}(z)|\geq m_{\varepsilon}r>0$.
Now for any $y\in B_{\tau r}(z)$ we have

$|\mathbf{u}(y)-\mathbf{u}(z)|\leq\mathrm{Lip}(\mathbf{u})|y-z|<\mathrm{Lip}(\mathbf{u})\tau
r<\frac{m_{\varepsilon}r}{2},$

provided that $\tau$ is small enough. Hence we must have
$|\mathbf{u}(y)|>\frac{m_{\varepsilon}r}{2}>0$. This gives the upper estimate
in (3.8).

To prove the estimate from below, suppose to the contrary that there exists a
sequence of points $x_{k}\in\partial U$ and radii $r_{k}\to 0$ such that

$\big{|}\\{|\mathbf{u}|=0\\}\cap
B_{r_{k}}(x_{k})\big{|}<\frac{1}{k}|B_{r_{k}}(x_{k})|=\frac{1}{k}r_{k}^{n}\,|B_{1}|.$

Now let us define

$\mathbf{u}_{k}(x)=\frac{\mathbf{u}(x_{k}+r_{k}x)}{r_{k}}.$

Note that $\mathbf{u}_{k}(0)=\mathbf{u}(x_{k})=0$, and thus $\mathbf{u}_{k}$
is uniformly bounded and uniformly Lipschitz in $B_{1}=B_{1}(0)$ due to
Lipschitz regularity of $\mathbf{u}$. Also

$\big{|}\\{|\mathbf{u}_{k}|=0\\}\cap
B_{1}\big{|}=\frac{1}{r_{k}^{n}}\big{|}\\{|\mathbf{u}|=0\\}\cap
B_{r_{k}}(x_{k})\big{|}\underset{k\to\infty}{\longrightarrow}0.$

Let $v_{k}^{i}$ be a $p$-harmonic function in $B_{1/2}$ with boundary data
$v_{k}^{i}=u_{k}^{i}$ on $\partial B_{1/2}$. Then
$h_{k}^{i}(x)=r_{k}v_{k}^{i}(\frac{x-x_{k}}{r_{k}})$ is a $p$-harmonic
function in $B_{r_{k}/2}(x_{k})$ with boundary data $h_{k}^{i}=u^{i}$ on
$\partial B_{r_{k}/2}(x_{k})$. Now, similarly to the proof of Theorem 2.7, we
can show that

(3.9) $\displaystyle\int_{B_{1/2}}|\nabla(u_{k}^{i}-v_{k}^{i})|^{p}\,dx$
$\displaystyle=\frac{1}{r_{k}^{n}}\int_{B_{r_{k}/2}(x_{k})}|\nabla(u^{i}-h_{k}^{i})|^{p}\,dx$
$\displaystyle\leq\frac{C}{\varepsilon}\frac{1}{r_{k}^{n}}\big{|}\\{|\mathbf{u}|=0\\}\cap
B_{r_{k}}(x_{k})\big{|}\underset{k\to\infty}{\longrightarrow}0.$

(Note that the constant $C$ does not depend on the radius $r_{k}$ or the point
$x_{k}$.)

Since $u_{k}^{i}$ and therefore $v_{k}^{i}$ are uniformly Lipschitz in
$B_{1/4}$, we may assume that $u_{k}^{i}\to u_{0}^{i}$ and $v_{k}^{i}\to
v_{0}^{i}$ uniformly in $B_{1/4}$. Observe that $\Delta_{p}v_{0}^{i}=0$, and
(3.9) implies that $u_{0}^{i}=v_{0}^{i}+C$ for some constant $C$. Thus
$\Delta_{p}u_{0}^{i}=0$ in $B_{1/4}$ and from the strong minimum principle it
follows $u_{0}^{i}\equiv 0$ in $B_{1/4}$, since $u_{0}^{i}\geq 0$ and
$u_{0}^{i}(0)=\lim u_{k}^{i}(0)=0$. On the other hand the nondegeneracy
property, Theorem 3.1, implies that (since $x_{k}$ is not in the interior of
$\\{|\mathbf{u}|=0\\}$)

$\|\mathbf{u}_{k}\|_{L^{\infty}(B_{1/4})}=\frac{1}{r_{k}}\|\mathbf{u}\|_{L^{\infty}(B_{r_{k}/4}(x_{k}))}\geq
m_{\varepsilon}/2>0.$

Therefore we get $\|\mathbf{u}_{0}\|_{L^{\infty}(B_{1/4})}\geq
m_{\varepsilon}/2$, which is a contradiction. ∎

Hence we can apply the results in section 4 of [3] and in section 3 of [4] to
conclude (see also sections 5 and 6 of [7])

######  Theorem 3.5.

Let $\mathbf{u}=\mathbf{u}_{\varepsilon}$ be a minimizer of $J_{\varepsilon}$
over $V$. Then we have

  1. (1)

The $(n-1)$-dimensional Hausdorff measure of $\partial E$ is locally finite,
i.e. $\mathcal{H}^{n-1}(\Omega^{\prime}\cap\partial E)<\infty$ for every
$\Omega^{\prime}\subset\subset\Omega$. Moreover, there exist positive
constants $c_{\varepsilon},C_{\varepsilon}$, depending on
$n,p,\Omega,\Omega^{\prime},\varepsilon$, such that for each ball
$B_{r}(x)\subset\Omega^{\prime}$ with $x\in\partial E$ we have

$c_{\varepsilon}r^{n-1}\leq\mathcal{H}^{n-1}(B_{r}(x)\cap\partial E)\leq
C_{\varepsilon}r^{n-1}.$

  2. (2)

There exist Borel functions $q^{i}=q_{\varepsilon}^{i}$ such that

$\Delta_{p}u^{i}=q^{i}\,\mathcal{H}^{n-1}\mathbin{\vrule
height=6.88889pt,depth=0.0pt,width=0.55974pt\vrule
height=0.55974pt,depth=0.0pt,width=5.59721pt}\partial E,$

that is, for any $\zeta\in C_{0}^{\infty}(\Omega)$ we have

$-\int_{\Omega}A[u^{i}]\cdot\nabla\zeta\,dy=\int_{\partial E}\zeta
q^{i}\,d\mathcal{H}^{n-1}.$

  3. (3)

For $\mathcal{H}^{n-1}$-a.e. points $x\in\partial E$ we have

$c_{\varepsilon}\leq\sum_{i=1}^{m}q^{i}(x)\leq C_{\varepsilon}.$

  4. (4)

For $\mathcal{H}^{n-1}$-a.e. points $x\in\partial E$ an outward unit normal
$\nu=\nu_{E}(x)$ is defined, and

$u^{i}(x+y)=(q^{i}(x))^{\frac{1}{p-1}}(y\cdot\nu)^{+}+o(|y|),$

which allows us to define $A_{\nu}u^{i}(x)=q^{i}(x)$ at those points.

  5. (5)

The reduced boundary $\partial_{\mathrm{red}}E$ satisfies
$\mathcal{H}^{n-1}(\partial E-\partial_{\mathrm{red}}E)=0$.

##  4\. The Original Problem

In this section we will show that for $\varepsilon>0$ small enough, a
minimizer of $J_{\varepsilon}$ over $V$ satisfies
$\big{|}\\{|\mathbf{u}_{\varepsilon}|>0\\}\big{|}=1$, and hence it can be
regarded as a solution to our original problem (1.1). Remember that

$U=U_{\varepsilon}=\\{|\mathbf{u}_{\varepsilon}|>0\\},\qquad
E=E_{\varepsilon}=\\{|\mathbf{u}_{\varepsilon}|=0\\}.$

Note that by Lemma 2.9, the free boundary $\partial E$ has a positive distance
from the fixed boundary $\partial\Omega$. We say $x\in\partial E$ is a regular
point of the free boundary if it satisfies (3) and (4) in Theorem 3.5. The set
of such regular points of the free boundary will be denoted by
$\mathcal{R}=\mathcal{R}_{\varepsilon}$; Theorem 3.5 shows that
$\mathcal{H}^{n-1}(\partial E-\mathcal{R})=0$.

######  Lemma 4.1.

There is a constant $C>0$, independent of $\varepsilon$, such that

$\inf_{\mathcal{R}_{\varepsilon}}\Big{(}\sum_{i\leq
m}q_{\varepsilon}^{i}\Big{)}\leq C.$

######  Remark.

Note that $\sum_{i\leq m}q_{\varepsilon}^{i}\geq c_{\varepsilon}>0$ by Theorem
3.5.

###### Proof.

Let $\Omega^{\prime}\subset\subset\Omega$ be a smooth open set with
$|\Omega-\Omega^{\prime}|=1$. Let $\mathbf{u}_{0}$ be a vector-valued function
on $\Omega-\Omega^{\prime}$ that satisfies the equation
$\Delta_{p}u_{0}^{i}=0$, and takes the boundary values $\boldsymbol{\varphi}$
on $\partial\Omega$ and $0$ on $\partial\Omega^{\prime}$. Then for some small
enough $\delta_{0}$ we have $\mathbf{u}_{0}\in V_{\delta_{0}}\subset V$; hence

$\displaystyle
C=\int_{\partial\Omega}\Gamma(x,A_{\nu}\mathbf{u}_{0})\,d\sigma+1$
$\displaystyle=J_{\varepsilon}(\mathbf{u}_{0})\geq
J_{\varepsilon}(\mathbf{u}_{\varepsilon})$
$\displaystyle=\int_{\partial\Omega}\Gamma(x,A_{\nu}\mathbf{u}_{\varepsilon})\,d\sigma+f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}_{\varepsilon}|>0\\}\big{|}\big{)}$
$\displaystyle\geq\int_{\partial\Omega}\sum_{i=1}^{m}\psi_{i}A_{\nu}u_{\varepsilon}^{i}-C\,d\sigma+f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}_{\varepsilon}|>0\\}\big{|}\big{)}$
$\displaystyle\geq C\sum_{i=1}^{m}\int_{\Omega}|\nabla
u_{\varepsilon}^{i}|^{p}\,dx-C+f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}_{\varepsilon}|>0\\}\big{|}\big{)}$
$\displaystyle\geq-C+\frac{1}{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}_{\varepsilon}|>0\\}\big{|}-1\big{)},$

where we have used (1.2) and Lemma 2.2. Thus we get the bound

$|U|=\big{|}\\{|\mathbf{u}_{\varepsilon}|>0\\}\big{|}\leq 1+C\varepsilon.$

Note that $J_{\varepsilon}(\mathbf{u}_{0})$, and thus $C$, does not depend on
$\varepsilon$ due to the definition of $f_{\varepsilon}$. As a result, we have
a lower bound for the volume of $E$. Hence, by the isoperimetric inequality,
we have a lower bound for $\mathcal{H}^{n-1}(\partial E)$, independent of
$\varepsilon$. Now note that (keep in mind that $\nu_{E}$ points to the
interior of $U$)

$\int_{\partial\Omega}A_{\nu}u_{\varepsilon}^{i}\,d\sigma-\int_{\partial
E}A_{\nu}u_{\varepsilon}^{i}\,d\mathcal{H}^{n-1}=\int_{U}\Delta_{p}u_{\varepsilon}^{i}\,dx=0.$

Therefore we get

$\displaystyle\int_{\partial
E}A_{\nu}u_{\varepsilon}^{i}\,d\mathcal{H}^{n-1}=\int_{\partial\Omega}A_{\nu}u_{\varepsilon}^{i}\,d\sigma$
$\displaystyle=\int_{\Omega}\Delta_{p}u_{\varepsilon}^{i}\,dx$
$\displaystyle\leq
C+C\int_{\partial\Omega}\psi_{i}A_{\nu}u_{\varepsilon}^{i}\,d\sigma,$

where the last inequality follows from the remark below Lemma 2.2. Thus we
have

$\displaystyle\inf_{\mathcal{R}_{\varepsilon}}\Big{(}\sum_{i\leq
m}A_{\nu}u_{\varepsilon}^{i}\Big{)}\mathcal{H}^{n-1}(\partial E)$
$\displaystyle\leq\int_{\partial E}\sum_{i\leq
m}A_{\nu}u_{\varepsilon}^{i}\,d\mathcal{H}^{n-1}$ $\displaystyle\leq
C+C\int_{\partial\Omega}\sum_{i\leq
m}\psi_{i}A_{\nu}u_{\varepsilon}^{i}\,d\sigma$ (by (1.2)) $\displaystyle\leq
C+C\int_{\partial\Omega}\Gamma\big{(}x,A_{\nu}\mathbf{u}_{\varepsilon}\big{)}\,d\sigma$
$\displaystyle\leq C+CJ_{\varepsilon}(\mathbf{u}_{\varepsilon})\leq
C+CJ_{\varepsilon}(\mathbf{u}_{0})\leq C,$

which gives the desired (noting that $q^{i}=A_{\nu}u_{\varepsilon}^{i}$ by
Theorem 3.5). ∎

######  Lemma 4.2.

For small enough $\varepsilon$ we have

$\big{|}\\{|\mathbf{u}_{\varepsilon}|>0\\}\big{|}\geq 1.$

###### Proof.

Consider a point $z_{0}\in\Omega^{c}$ which has distance $\delta_{0}$ from
$\partial\Omega$. Then the ball $B_{\delta_{0}}(z_{0})$ is an exterior tangent
ball to $\partial\Omega$. Let $t=t(\varepsilon)$ be the first time at which
$\partial B_{\delta_{0}+t}(z_{0})$ intersects
$\partial\\{|\mathbf{u}_{\varepsilon}|=0\\}$, at a point
$x_{0}=x_{0}(\varepsilon)$. Now let $v$ be a $p$-harmonic function in
$B_{\delta_{0}+t}(z_{0})-\overline{B}_{\delta_{0}}(z_{0})$ with boundary
values $0$ on $\partial B_{\delta_{0}+t}(z_{0})$ and $c_{0}$ on $\partial
B_{\delta_{0}}(z_{0})$, where
$c_{0}=\min_{i}\min_{\partial\Omega}\varphi^{i}>0$. Then on
$\partial\big{(}\Omega\cap B_{\delta_{0}+t}(z_{0})\big{)}$ we have $v\leq
u^{i}$; so by the maximum principle we have $v\leq u^{i}$ in $\Omega\cap
B_{\delta_{0}+t}(z_{0})$. However, by an easy modification of the proof of
Hopf’s lemma (Lemma 2.4), we can see that

$v(x)\geq cc_{0}\,\mathrm{dist}\big{(}x,\partial
B_{\delta_{0}+t}(z_{0})\big{)},$

where the constant $c$ only depends on $n,p,\delta_{0}$. Therefore, for points
$x$ in the line segment between $x_{0},z_{0}$ we have

$u^{i}(x)\geq v(x)\geq cc_{0}\,\mathrm{dist}\big{(}x,\partial
B_{\delta_{0}+t}(z_{0})\big{)}=cc_{0}|x-x_{0}|.$

Now consider the ball $B_{r}(x_{0})$ for small enough $r$. Then we have

$\frac{1}{r}\sup_{B_{r/2}(x_{0})}u^{i}\geq\frac{1}{r}cc_{0}\frac{r}{2}=\frac{cc_{0}}{2},$

independently of $\varepsilon$.

Let $\mathbf{h}$ be the vector-valued function which satisfies
$\Delta_{p}h^{i}=0$ in $B_{r}(x_{0})$, and is equal to $\mathbf{u}$ in
$\Omega-B_{r}(x_{0})$. By Lemma 2.5 and the fact that $h^{i}\geq u^{i}$ we
have

$\displaystyle\int_{B_{r}(x_{0})}|\nabla(u^{i}-h^{i})|^{p}\,dy$
$\displaystyle\geq
C\Big{(}\frac{1}{r}\sup_{B_{r/2}(x_{0})}h^{i}\Big{)}^{p}\cdot\big{|}B_{r}(x_{0})\cap\\{u^{i}=0\\}\big{|}$
$\displaystyle\geq
C\Big{(}\frac{1}{r}\sup_{B_{r/2}(x_{0})}u^{i}\Big{)}^{p}\cdot\big{|}B_{r}(x_{0})\cap\\{u^{i}=0\\}\big{|}$
$\displaystyle\geq C\big{|}B_{r}(x_{0})\cap\\{u^{i}=0\\}\big{|}\geq
C\big{|}B_{r}(x_{0})\cap\\{|\mathbf{u}|=0\\}\big{|}.$

Next let $\mathbf{v}$ be the function given by Lemma 2.6 for $B_{r}(x_{0})$.
We know that $J_{\varepsilon}(\mathbf{u})\leq J_{\varepsilon}(\mathbf{v})$.
Then similarly to the proof of Theorem 2.7 we can see that

$\displaystyle C\sum_{i\leq
m}\int_{B_{r}(x_{0})}|\nabla(u^{i}-h^{i})|^{p}\,dy$
$\displaystyle\leq\int_{\partial\Omega}\Gamma(x,A_{\nu}\mathbf{u})-\Gamma(x,A_{\nu}\mathbf{v})\,d\sigma$
$\displaystyle\leq
f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{v}|>0\\}\big{|}\big{)}-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}.$

(A closer inspection of the proof of Theorem 2.7 reveals that the constant $C$
in the above estimate only depends on
$n,p,\Omega,\boldsymbol{\varphi},\Gamma$.)

Now suppose to the contrary that $\big{|}\\{|\mathbf{u}|>0\\}\big{|}<1$. Then,
since $0\leq u^{i}\leq v^{i}$, and outside of $B_{r}(x_{0})$, $|\mathbf{u}|=0$
implies $|\mathbf{v}|=0$, we have

$\big{|}\\{|\mathbf{v}|>0\\}\big{|}\leq\big{|}\\{|\mathbf{u}|>0\\}\big{|}+\big{|}B_{r}(x_{0})\cap\\{|\mathbf{u}|=0\\}\big{|}<1$

for small enough $r$. Hence (using the monotonicity of $f_{\varepsilon}$) we
have

$\displaystyle
f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{v}|>0\\}\big{|}\big{)}-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}$
$\displaystyle\qquad\qquad\leq
f_{\varepsilon}\Big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}+\big{|}B_{r}(x_{0})\cap\\{|\mathbf{u}|=0\\}\big{|}\Big{)}-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}$
$\displaystyle\qquad\qquad=\varepsilon\big{|}B_{r}(x_{0})\cap\\{|\mathbf{u}|=0\\}\big{|}.$

Combining this estimate with the estimates of the above paragraph, and using
(3.8), we obtain

$0<C\big{|}B_{r}(x_{0})\cap\\{|\mathbf{u}|=0\\}\big{|}\leq\varepsilon\big{|}B_{r}(x_{0})\cap\\{|\mathbf{u}|=0\\}\big{|},$

which gives a positive lower bound for $\varepsilon$, and results in a
contradiction. ∎

######  Theorem 4.3.

When $\varepsilon$ is small enough we have

$\big{|}\\{|\mathbf{u}_{\varepsilon}|>0\\}\big{|}=1.$

###### Proof.

By the above lemma we only need to show that
$\big{|}\\{|\mathbf{u}_{\varepsilon}|>0\\}\big{|}\leq 1$. To this end, we will
compare $\mathbf{u}_{\varepsilon}$ with a suitable perturbation of itself. Let
$x_{0}\in\mathcal{R}$, and let $\rho:\mathbb{R}\to\mathbb{R}$ be a nonnegative
smooth function supported in $(0,1)$. For small enough $r,\lambda>0$ we
consider the vector field

$T_{r}(x):=\begin{cases}x+r\lambda\rho(|x-x_{0}|/r)\nu(x_{0})&\textrm{if }x\in
B_{r}(x_{0}),\\\ x&\text{elsewhere}.\end{cases}$

Here, $\nu(x_{0})$ is the outward normal vector provided in (4) of Theorem
3.5. We can easily see that for $x$ in $B_{r}(x_{0})$ we have

(4.1)
$DT_{r}(x)\,\cdot=I\cdot+\,\lambda\rho^{\prime}(|x-x_{0}|/r)\frac{\langle
x-x_{0},\cdot\,\rangle}{|x-x_{0}|}\nu(x_{0}),$

where $I$ is the identity matrix. Hence, if $\lambda$ is small enough, $T_{r}$
is a diffeomorphism that maps $B_{r}(x_{0})$ onto itself.

Now consider

$\mathbf{v}_{r}(x):=\mathbf{u}(T_{r}^{-1}(x))$

for $r>0$ small enough. Similarly to the proof of Theorem 3.1, we consider the
vector-valued function $\mathbf{w}$ whose components minimize the Dirichlet
$p$-energy subject to the condition

$w^{i}\leq 0\;\textrm{ on
}\;\\{\mathbf{u}=0\\}\cup\big{(}\overline{B}_{r}(x_{0})\cap\\{\mathbf{v}_{r}=0\\}\big{)}.$

With a calculation similar to (3) and (2) we get

$\displaystyle 0\leq J_{\varepsilon}(\mathbf{w})-J_{\varepsilon}(\mathbf{u})$
$\displaystyle\leq C\sum_{i=1}^{m}\int_{\Omega}|\nabla w^{i}|^{p}-|\nabla
u^{i}|^{p}\,dx$
$\displaystyle\qquad\qquad+f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{w}|>0\\}\big{|}\big{)}-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}$
(4.2) $\displaystyle\leq C\sum_{i=1}^{m}\int_{B_{r}(x_{0})}|\nabla
v_{r}^{i}|^{p}-|\nabla u^{i}|^{p}\,dx$
$\displaystyle\qquad\qquad+f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{w}|>0\\}\big{|}\big{)}-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)},$

where in the last inequality we have compared the Dirichlet $p$-energy of
$\mathbf{w}$ with that of
$\mathbf{v}_{r}\chi_{B_{r}(x_{0})}+\mathbf{u}\chi_{\Omega-B_{r}(x_{0})}$.

Now notice that

$\displaystyle\int_{B_{r}(x_{0})}|\nabla v_{r}^{i}|^{p}\,dx$
$\displaystyle=\int_{B_{r}(x_{0})}\big{|}DT_{r}(T_{r}^{-1}(x))^{-1}\nabla
u^{i}(T_{r}^{-1}(x))\big{|}^{p}\,dx$
$\displaystyle=\int_{B_{r}(x_{0})}\big{|}DT_{r}(y)^{-1}\nabla
u^{i}(y)\big{|}^{p}\,|\det DT_{r}(y)|\,dy$ ($z=\frac{y-x_{0}}{r}$)
$\displaystyle=r^{n}\int_{B_{1}}\big{|}DT_{r}(y)^{-1}\nabla
u^{i}(y)\big{|}^{p}\,|\det DT_{r}(y)|\,dz.$

From (4.1), for small enough $\lambda$ we can write

$\displaystyle DT_{r}(y)^{-1}$
$\displaystyle=I+\Big{(}\sum_{k=1}^{\infty}(-1)^{k}\lambda^{k}\rho^{\prime}(|z|)^{k}\frac{\langle
z,\nu\rangle^{k-1}}{|z|^{k-1}}\Big{)}\frac{\langle
z,\cdot\,\rangle}{|z|}\nu(x_{0})$ (4.3)
$\displaystyle=I-\lambda\rho^{\prime}(|z|)\frac{\langle
z,\cdot\,\rangle}{|z|}\nu(x_{0})+\lambda^{2}g(\lambda,z)\frac{\langle
z,\cdot\,\rangle}{|z|}\nu(x_{0})$

for some $g$. Hence we have

$DT_{r}(y)^{-1}\nabla u^{i}(y)=\nabla
u^{i}(y)-\lambda\rho^{\prime}(|z|)\frac{\langle z,\nabla
u^{i}(y)\rangle}{|z|}\nu(x_{0})+O(\lambda^{2}).$

Thus

$\big{|}DT_{r}(y)^{-1}\nabla u^{i}(y)\big{|}^{2}=|\nabla
u^{i}(y)|^{2}-2\lambda\rho^{\prime}(|z|)\frac{\langle z,\nabla
u^{i}(y)\rangle}{|z|}\langle\nu(x_{0}),\nabla u^{i}(y)\rangle+O(\lambda^{2}),$

and therefore

$\displaystyle\big{|}DT_{r}(y)^{-1}\nabla u^{i}(y)\big{|}^{p}$
$\displaystyle=|\nabla
u^{i}(y)|^{p}\Big{(}1-p\lambda\rho^{\prime}(|z|)\frac{\langle z,\nabla
u^{i}(y)\rangle}{|z||\nabla u^{i}(y)|^{2}}\langle\nu(x_{0}),\nabla
u^{i}(y)\rangle\Big{)}$ $\displaystyle\qquad+O(\lambda^{2}).$

Also, we have (noting that $DT_{r}$ is the identity matrix plus a rank $1$
matrix)

$|\det DT_{r}(y)|=1+\lambda\rho^{\prime}(|z|)\frac{\langle
z,\nu(x_{0})\rangle}{|z|}.$

All these together, we obtain (remember that $y=x_{0}+rz$)

$\displaystyle r^{-n}\int_{B_{r}(x_{0})}|\nabla v_{r}^{i}|^{p}-|\nabla
u^{i}|^{p}\,dx$ $\displaystyle\qquad=\lambda\int_{B_{1}}|\nabla
u^{i}(y)|^{p}\rho^{\prime}(|z|)\left(\frac{\langle
z,\nu(x_{0})\rangle}{|z|}-p\frac{\langle z,\nabla u^{i}(y)\rangle\langle\nabla
u^{i}(y),\nu(x_{0})\rangle}{|z||\nabla u^{i}(y)|^{2}}\right)dz$
$\displaystyle\qquad\quad+O(\lambda^{2}).$

Now consider the blowup sequence $\mathbf{u}_{r}(z):=\mathbf{u}(x_{0}+rz)/r$.
We know that as $r\to 0$ (see [4])

$\begin{split}&\\{u_{r}^{i}>0\\}\cap B_{1}\to\\{z:z\cdot\nu(x_{0})>0\\}\cap
B_{1},\\\ &\nabla u^{i}(y)=\nabla
u_{r}^{i}(z)\to(q^{i}(x_{0}))^{\frac{1}{p-1}}\nu(x_{0})\chi_{\\{z\cdot\nu(x_{0})>0\\}},\qquad\text{a.e.
in }B_{1}.\end{split}$

Therefore we get

$\displaystyle r^{-n}\int_{B_{r}(x_{0})}|\nabla v_{r}^{i}|^{p}-|\nabla
u^{i}|^{p}\,dx$ $\displaystyle\qquad\underset{r\to
0}{\longrightarrow}-(p-1)\lambda|q^{i}(x_{0})|^{\frac{p}{p-1}}\int_{B_{1}\cap\\{z\cdot\nu(x_{0})>0\\}}\rho^{\prime}(|z|)\frac{\langle
z,\nu(x_{0})\rangle}{|z|}\,dz+O(\lambda^{2}).$

Note that the formula (4.3) for $(DT_{r})^{-1}$ does not depend on $r$, and
the function $\cdot\mapsto|\cdot|^{p}$ is continuous; so the $O(\lambda^{2})$
term converges to an $O(\lambda^{2})$ term as $r\to 0$. Next note that

$\mathrm{div}(\rho(|z|)\nu)=\frac{\rho^{\prime}(|z|)}{|z|}\langle
z,\nu\rangle.$

Thus (noting that $\rho(|z|)$ is zero near $\partial B_{1}$)

$\displaystyle\int_{B_{1}\cap\\{z\cdot\nu(x_{0})>0\\}}\rho^{\prime}(|z|)\frac{\langle
z,\nu(x_{0})\rangle}{|z|}\,dz$
$\displaystyle=-\int_{B_{1}\cap\\{z\cdot\nu(x_{0})=0\\}}\rho(|z|)\,dz$
$\displaystyle=-\omega_{n-1}\int_{0}^{1}\rho(t)t^{n-1}\,dt=-C_{\rho}\omega_{n-1},$

where $\omega_{n-1}$ is the volume of the $(n-1)$-dimensional ball of radius
$1$, and $C_{\rho}$ depends only on $\rho$. Hence we can write

$\displaystyle\int_{B_{r}(x_{0})}|\nabla v_{r}^{i}|^{p}-|\nabla
u^{i}|^{p}\,dx$ $\displaystyle\qquad=\big{[}(p-1)\lambda
C_{\rho}\omega_{n-1}|q^{i}(x_{0})|^{\frac{p}{p-1}}+O(\lambda^{2})\big{]}r^{n}+o(r^{n}).$

On the other hand,

$\begin{split}&\lim_{r\to
0}r^{-n}\big{|}B_{r}(x_{0})\cap\\{|\mathbf{v}_{r}|>0\\}\big{|}=\lim_{r\to
0}r^{-n}\int_{\\{|\mathbf{v}_{r}|>0\\}\cap B_{r}(x_{0})}dx\\\
&\qquad=\lim_{r\to 0}r^{-n}\int_{\\{|\mathbf{u}|>0\\}\cap B_{r}(x_{0})}|\det
DT_{r}(y)|\,dy\\\
&\qquad=\int_{B_{1}\cap\\{z\cdot\nu(x_{0})>0\\}}1+\lambda\rho^{\prime}(|z|)\frac{\langle
z,\nu(x_{0})\rangle}{|z|}\,dz\\\
&\qquad=\frac{1}{2}\omega_{n}-\lambda\omega_{n-1}\int_{0}^{1}\rho(t)t^{n-1}\,dt=\frac{1}{2}\omega_{n}-\lambda
C_{\rho}\omega_{n-1}.\end{split}$

Thus for
$A_{0}:=\big{(}\\{|\mathbf{u}|>0\\}-B_{r}(x_{0})\big{)}\cup\big{(}\\{|\mathbf{v}_{r}|>0\\}\cap
B_{r}(x_{0})\big{)}$ we have

$\displaystyle|A_{0}|-\big{|}\\{|\mathbf{u}|>0\\}\big{|}$
$\displaystyle=\big{|}B_{r}(x_{0})\cap\\{|\mathbf{v}_{r}|>0\\}\big{|}-\big{|}B_{r}(x_{0})\cap\\{|\mathbf{u}|>0\\}\big{|}$
$\displaystyle=-\lambda C_{\rho}\omega_{n-1}r^{n}+o(r^{n}).$

In addition, it is easy to see that $\\{|\mathbf{w}|>0\\}\subset A_{0}$.

Now suppose to the contrary that $\big{|}\\{|\mathbf{u}|>0\\}\big{|}>1$. Then
we can choose $r$ small enough so that

$|A_{0}|=\big{|}\\{|\mathbf{u}|>0\\}\big{|}-\lambda
C_{\rho}\omega_{n-1}r^{n}+o(r^{n})>1.$

Therefore, using the monotonicity of $f_{\varepsilon}$ we get

$\begin{split}f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{w}|>0\\}\big{|}\big{)}&-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}\leq
f_{\varepsilon}(|A_{0}|)-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}\\\
&=\frac{1}{\varepsilon}\big{(}|A_{0}|-\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}=-\frac{1}{\varepsilon}\lambda
C_{\rho}\omega_{n-1}r^{n}+o(r^{n}).\end{split}$

Finally, by putting all these estimates in (4.2), we obtain

$\displaystyle 0$ $\displaystyle\leq C\sum_{i=1}^{m}\int_{B_{r}(x_{0})}|\nabla
v_{r}^{i}|^{p}-|\nabla
u^{i}|^{p}\,dx+f_{\varepsilon}(|A_{0}|)-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}$
$\displaystyle=\big{[}(p-1)\lambda
C_{\rho}\omega_{n-1}\sum_{i=1}^{m}|q^{i}(x_{0})|^{\frac{p}{p-1}}+O(\lambda^{2})\big{]}r^{n}-\frac{1}{\varepsilon}\lambda
C_{\rho}\omega_{n-1}r^{n}+o(r^{n}).$

Dividing by $r^{n}$ and letting $r\to 0$, and then dividing by $\lambda$ and
letting $\lambda\to 0$, we get

$\frac{1}{\varepsilon}\leq(p-1)\sum_{i=1}^{m}|q^{i}(x_{0})|^{\frac{p}{p-1}}.$

Now if we choose $x_{0}$ such that

$\sum_{i\leq
m}q^{i}(x_{0})\leq\inf_{\mathcal{R}_{\varepsilon}}\Big{(}\sum_{i\leq
m}q^{i}\Big{)}+1,$

then by Lemma 4.1 (and the equivalence of all norms on the finite-dimensional
space $\mathbb{R}^{m}$) we have $\sum_{i\leq
m}|q^{i}(x_{0})|^{\frac{p}{p-1}}\leq C$, independently of $\varepsilon$.
However, this implies that $\varepsilon$ has a positive lower bound, which is
a contradiction. ∎

##  5\. Regularity of the free boundary (case $p=2$)

We are going to show that $\mathcal{R}$ is an analytic hypersurface when
$p=2$. To see this, we first derive the free boundary condition, also known as
the optimality condition, in the following lemma. We perturb the optimal set
$\Omega$ and compute the first variation of the energy functional
$J_{\varepsilon}$. To perform this computation, it is crucial to ensure that
the $p$-harmonic solution within the perturbed domain is differentiable with
respect to the perturbation parameter. When $p=2$, this can be established
through the implicit function theorem . However, it is noteworthy that for
$p\neq 2$ the proof does not hold, primarily due to the ill-posedness of the
derivative of the map $u\mapsto\Delta_{p}u$.

######  Lemma 5.1.

Let $\mathbf{u}$ be a solution of the minimization problem (1.1) for $p=2$.
Let $h^{i}$ be the solution of

$\begin{cases}\Delta h^{i}=0&\text{in }\Omega-E,\\\ h^{i}=0&\text{on }E,\\\
h^{i}=\partial_{\xi_{i}}\Gamma(x,\partial_{\nu}\mathbf{u})&\text{on
}\partial\Omega.\end{cases}$

Then, on the regular part of the free boundary we have

(5.1) $\sum_{i=1}^{m}\partial_{\nu}h^{i}\partial_{\nu}u^{i}=C$

for some positive constant $C$.

###### Proof.

Let $x_{1}$ and $x_{2}$ be two regular points in $\mathcal{R}$ with
corresponding unit normal vectors $\nu(x_{1})$ and $\nu(x_{2})$. Also, let
$\rho:\mathbb{R}\to\mathbb{R}$ be a nonnegative smooth function supported in
$(0,1)$. Similarly to the proof of Theorem 4.3 we define the vector field

$T_{r,\lambda}(x):=\begin{cases}x-r\lambda\rho(|x-x_{1}|/r)\nu(x_{1})&\textrm{if
}x\in B_{r}(x_{1}),\\\ x+r\lambda\rho(|x-x_{2}|/r)\nu(x_{2})&\textrm{if }x\in
B_{r}(x_{2}),\\\ x&\text{elsewhere},\end{cases}$

for small enough $r,\lambda>0$ (which makes $T_{r,\lambda}$ a diffeomorphism
from $B_{r}(x_{a})$ onto itself for $a=1,2$).

Now for some fixed $r>0$ let $E_{\lambda}=T_{r,\lambda}^{-1}(E)$, and assume
that $\mathbf{w}_{\lambda}$ solves

$\begin{cases}\Delta w_{\lambda}^{i}=0&\text{in }\Omega-E_{\lambda},\\\
w_{\lambda}^{i}=\varphi^{i}&\text{on }\partial\Omega,\\\
w_{\lambda}^{i}=0&\text{on }\partial E_{\lambda}.\end{cases}$

Define $\mathbf{v}_{\lambda}(y):=\mathbf{w}_{\lambda}(T_{r,\lambda}^{-1}(y))$.
We are going to show that $\lambda\mapsto\mathbf{v}_{\lambda}$ is a $C^{1}$
map from a neighborhood of $\lambda=0$ into $W^{1,2}(\Omega-E)$. We know that
each $v_{\lambda}^{i}$ satisfies an elliptic PDE of the form

$F[v,\lambda]=F(D_{y}^{2}v,\nabla_{y}v,y,\lambda)=0,\qquad\text{ in
}\;U=\Omega-E.$

We also know that $F=\Delta$ when $y\notin B_{r}(x_{1})\cup B_{r}(x_{2})$ or
when $\lambda=0$. In addition, we can consider $F$ as a $C^{1}$ map

$\displaystyle F:W^{1,2}(U)\times\mathbb{R}$ $\displaystyle\to W^{-1,2}(U),$
$\displaystyle(v,\lambda)$ $\displaystyle\mapsto F[v,\lambda]$

where $U=\Omega-E$.

Now we employ the implicit function theorem to show that
$\lambda\mapsto\mathbf{v}_{\lambda}$ is $C^{1}$. This can be readily deduced
from the fact that

$\partial_{v}F|_{\lambda=0}:W_{0}^{1,2}(U)\to W^{-1,2}(U)$

is invertible, since we have

$\partial_{v}F|_{\lambda=0}\;\cdot=\frac{d}{ds}\Big{|}_{s=0}F[v+s\,\cdot\,,0]=\frac{d}{ds}\Big{|}_{s=0}\Delta(v+s\,\cdot\,)=\Delta\cdot.$

Therefore, $\mathbf{v}_{\lambda}=\mathbf{u}+\lambda\mathbf{u}_{0}+o(\lambda)$
in $W^{1,2}(U)$, where $\mathbf{u}_{0}\in W_{0}^{1,2}(U)$ solves

$0=\frac{d}{d\lambda}\Big{|}_{\lambda=0}F[v_{\lambda}^{i},\lambda]=\partial_{v}Fu_{0}^{i}+\partial_{\lambda}F.$

In other words

$\Delta u_{0}^{i}=-\partial_{\lambda}F|_{v=u^{i},\,\lambda=0}.$

Note that we also have
$\nabla\mathbf{v}_{\lambda}=\nabla\mathbf{u}+\lambda\nabla\mathbf{u}_{0}+o(\lambda)$,
since $\lambda\mapsto\mathbf{v}_{\lambda}$ is a $C^{1}$ map into $W^{1,2}(U)$;
so $\lambda\mapsto\nabla\mathbf{v}_{\lambda}$ is a $C^{1}$ map into
$L^{2}(U)$.

Now let $h^{i}$ be the solution of $\Delta h^{i}=0$ in $U=\Omega-E$ with
boundary data $h^{i}=\partial_{\xi_{i}}\Gamma(x,\partial_{\nu}\mathbf{u})$ on
$\partial\Omega$ and $h^{i}=0$ on $\partial E$. Then for small $\lambda>0$ we
have (note that for $p=2$ we have $A_{\nu}=\partial_{\nu}$)

$\displaystyle\int_{\partial\Omega}\Gamma(x,\partial_{\nu}\mathbf{v}_{\lambda})-\Gamma(x,\partial_{\nu}\mathbf{u})\,d\sigma$
$\displaystyle\hskip
28.45274pt=\int_{\partial\Omega}\sum_{i}\partial_{i}\Gamma(x,\partial_{\nu}\mathbf{u})(\partial_{\nu}v_{\lambda}^{i}-\partial_{\nu}u^{i})\,d\sigma+o(\lambda)$
$\displaystyle\hskip
28.45274pt=\lambda\int_{\partial\Omega}\sum_{i}\partial_{i}\Gamma(x,\partial_{\nu}\mathbf{u})\partial_{\nu}u_{0}^{i}\,d\sigma+o(\lambda)$
$\displaystyle\hskip
28.45274pt=\lambda\int_{\partial\Omega}\sum_{i}h^{i}\partial_{\nu}u_{0}^{i}\,d\sigma+o(\lambda)$
$\displaystyle\hskip 28.45274pt=\lambda\sum_{i}\int_{U}\nabla h^{i}\cdot\nabla
u_{0}^{i}+h^{i}\Delta u_{0}^{i}\,dx+o(\lambda)$ $\displaystyle\hskip
28.45274pt=-\lambda\sum_{i}\int_{{\textstyle(}B_{r}(x_{1})\cup
B_{r}(x_{2}){\textstyle)}-E}h^{i}\partial_{\lambda}F|_{v=u^{i},\,\lambda=0}\,dx+o(\lambda).$

Note that in the last line we have used the facts that $\Delta h^{i}=0$ in $U$
and $u_{0}^{i}=0$ on $\partial U=\partial\Omega\cup\partial E$. Also,
$\partial_{\lambda}F|_{v=u^{i},\,\lambda=0}=0$ outside $B_{r}(x_{1})\cup
B_{r}(x_{2})$, because in that region $F=\Delta$ for all $\lambda$.

Now let us extend $\mathbf{w}_{\lambda}$ to all of $\Omega$ by setting it
equal to $0$ on $E_{\lambda}$. Note that $w_{\lambda}^{i}$ is positive on
$\Omega-E_{\lambda}$ by the maximum principle. Hence

$\\{|\mathbf{w}_{\lambda}|>0\\}=\Omega-E_{\lambda}.$

Furthermore, similarly to the proof of Theorem 4.3, we obtain

$\displaystyle
f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{w}_{\lambda}|>0\\}\big{|}\big{)}$
$\displaystyle-
f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}\leq\frac{1}{\varepsilon}(|E|-|E_{\lambda}|)$
$\displaystyle=\frac{\lambda}{\varepsilon}\Big{(}\int_{B_{r}(x_{2})\cap\\{|\mathbf{u}|>0\\}}\rho^{\prime}(|x-x_{2}|)\frac{\langle
x-x_{2},\nu(x_{2})\rangle}{|x-x_{2}|}\,dx$
$\displaystyle\qquad\;-\int_{B_{r}(x_{1})\cap\\{|\mathbf{u}|>0\\}}\rho^{\prime}(|x-x_{1}|)\frac{\langle
x-x_{1},\nu(x_{1})\rangle}{|x-x_{1}|}\,dx\Big{)}=\frac{\lambda}{\varepsilon}o(r^{n}).$

Therefore if we compare the energy of $\mathbf{u}$ with $\mathbf{w}_{\lambda}$
(it is easy to see that $\mathbf{w}_{\lambda}\in V$) we get (in the second
equality below we use the fact that
$\mathbf{v}_{\lambda}=\mathbf{w}_{\lambda}$ near $\partial\Omega$)

$\displaystyle 0\leq
J_{\varepsilon}(\mathbf{w}_{\lambda})-J_{\varepsilon}(\mathbf{u})$
$\displaystyle=\int_{\partial\Omega}\Gamma(x,\partial_{\nu}\mathbf{w}_{\lambda})-\Gamma(x,\partial_{\nu}\mathbf{u})\,d\sigma$
$\displaystyle\qquad\qquad+f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{w}_{\lambda}|>0\\}\big{|}\big{)}-f_{\varepsilon}\big{(}\big{|}\\{|\mathbf{u}|>0\\}\big{|}\big{)}$
$\displaystyle=\int_{\partial\Omega}\Gamma(x,\partial_{\nu}\mathbf{v}_{\lambda})-\Gamma(x,\partial_{\nu}\mathbf{u})\,d\sigma+\frac{\lambda}{\varepsilon}o(r^{n})$
$\displaystyle=-\lambda\sum_{i}\int_{{\textstyle(}B_{r}(x_{1})\cup
B_{r}(x_{2}){\textstyle)}-E}h^{i}\partial_{\lambda}F|_{v=u^{i},\,\lambda=0}\,dx+o(\lambda)+\frac{\lambda}{\varepsilon}o(r^{n}).$

Hence if we divide by $\lambda$ and let $\lambda\to 0$ we obtain

(5.2) $0\leq-\sum_{i}\int_{{\textstyle(}B_{r}(x_{1})\cup
B_{r}(x_{2}){\textstyle)}-E}h^{i}\partial_{\lambda}F|_{v=u^{i},\,\lambda=0}\,dx+o(r^{n}).$

So we need to compute $\partial_{\lambda}F|_{v=u^{i},\,\lambda=0}$.

Next let us compute $F$ explicitly. Set $x=T_{r,\lambda}^{-1}(y)$ so that
$y=T_{r,\lambda}(x)$. To simplify the notation we suppress the $\lambda$ or
$r$ in the indices. We have $v^{i}(T(x))=v^{i}(y)=w^{i}(x)$. Hence

$\displaystyle\partial_{x_{k}}w^{i}$
$\displaystyle=\sum_{j}\partial_{y_{j}}v^{i}\partial_{x_{k}}T^{j},$
$\displaystyle\partial_{x_{k}x_{k}}^{2}w^{i}$
$\displaystyle=\sum_{j}\partial_{x_{k}}\big{(}\partial_{y_{j}}v^{i}\partial_{x_{k}}T^{j}\big{)}$
$\displaystyle=\sum_{j,\ell}\partial_{y_{j}y_{\ell}}^{2}v^{i}\partial_{x_{k}}T^{j}\partial_{x_{k}}T^{\ell}+\sum_{j}\partial_{y_{j}}v^{i}\partial_{x_{k}x_{k}}^{2}T^{j}.$

Therefore

$0=\Delta
w^{i}=\sum_{j,\ell,k}\partial_{y_{j}y_{\ell}}^{2}v^{i}\partial_{x_{k}}T^{j}\partial_{x_{k}}T^{\ell}+\sum_{j,k}\partial_{y_{j}}v^{i}\partial_{x_{k}x_{k}}^{2}T^{j}.$

It is easy to see that inside $B_{r}(x_{a})$ ($a=1,2$) we have

$\displaystyle\partial_{x_{k}}T^{j}$
$\displaystyle=\delta_{jk}+(-1)^{a}\lambda\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\nu^{j}(x_{a}),\qquad\big{(}z=\frac{x-x_{a}}{r}\big{)},$
$\displaystyle\partial_{x_{k}x_{k}}^{2}T^{j}$
$\displaystyle=(-1)^{a}\lambda\partial_{x_{k}}\big{(}\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\big{)}\nu^{j}(x_{a}).$

Thus

$\displaystyle F[v,\lambda]$
$\displaystyle=\sum_{j,\ell,k}\partial_{y_{j}y_{\ell}}^{2}v\partial_{x_{k}}T^{j}\partial_{x_{k}}T^{\ell}+\sum_{j,k}\partial_{y_{j}}v\partial_{x_{k}x_{k}}^{2}T^{j}$
$\displaystyle=\sum_{j,\ell,k}\big{[}\delta_{jk}+(-1)^{a}\lambda\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\nu^{j}(x_{a})\big{]}\big{[}\delta_{\ell
k}+(-1)^{a}\lambda\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\nu^{\ell}(x_{a})\big{]}\partial_{y_{j}y_{\ell}}^{2}v$
$\displaystyle\qquad\qquad+\sum_{j,k}\big{[}(-1)^{a}\lambda\partial_{x_{k}}\big{(}\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\big{)}\nu^{j}(x_{a})\big{]}\partial_{y_{j}}v$

in $B_{r}(x_{a})$ for $a=1,2$, and $F[v,\lambda]=\Delta v$ elsewhere. Now note
that

$\sum_{k}\partial_{x_{k}}\big{(}\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\big{)}=\sum_{k}\big{(}\rho^{\prime\prime}(|z|)\frac{z_{k}^{2}}{r|z|^{2}}+\rho^{\prime}(|z|)\frac{1}{r|z|}-\rho^{\prime}(|z|)\frac{z_{k}^{2}}{r|z|^{3}}\big{)}=\frac{1}{r}\rho^{\prime\prime}(|z|).$

Hence we get

$\partial_{\lambda}F|_{v=u^{i},\,\lambda=0}=(-1)^{a}\Big{(}2\rho^{\prime}(|z|)\sum_{j,k}\frac{z_{k}}{|z|}\nu^{j}(x_{a})\partial_{jk}^{2}u^{i}+\frac{1}{r}\rho^{\prime\prime}(|z|)\sum_{j}\nu^{j}(x_{a})\partial_{j}u^{i}\Big{)}$

in $B_{r}(x_{a})$ for $a=1,2$. Note that although a priori $z,u^{i}$ in the
above equation are functions of $y$, at $\lambda=0$ we have $y=x$, and thus we
can regard them as functions of $x$ too.

Let
$\mathbf{u}_{r}(z)=\frac{1}{r}\mathbf{u}(x_{a}+rz)=\frac{1}{r}\mathbf{u}(x)$
and $h_{r}^{i}(z)=\frac{1}{r}h^{i}(x_{a}+rz)=\frac{1}{r}h^{i}(x)$. Putting all
these in (5.2) we get (note that in the following integration by parts the
boundary term is zero, since $\rho$ is $0$ for $z$ near $\partial B_{1}$ and
$h^{i}$ is $0$ on $\partial E$)

$\displaystyle 0$
$\displaystyle\leq-\sum_{i}\int_{{\textstyle(}B_{r}(x_{1})\cup
B_{r}(x_{2}){\textstyle)}-E}h^{i}\partial_{\lambda}F|_{v=u^{i},\,\lambda=0}\,dx+o(r^{n})$
$\displaystyle=\sum_{a,i}(-1)^{a+1}\int_{B_{r}(x_{a})-E}h^{i}\Big{(}2\rho^{\prime}(|z|)\sum_{j,k}\frac{z_{k}}{|z|}\nu^{j}(x_{a})\partial_{jk}^{2}u^{i}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\frac{1}{r}\rho^{\prime\prime}(|z|)\sum_{j}\nu^{j}(x_{a})\partial_{j}u^{i}\Big{)}dx+o(r^{n})\allowdisplaybreaks$
$\displaystyle=\sum_{a,i}(-1)^{a+1}\int_{B_{r}(x_{a})-E}\Bigl{(}-2\sum_{k}\partial_{k}\Big{[}h^{i}\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\Big{]}\sum_{j}\nu^{j}(x_{a})\partial_{j}u^{i}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\frac{1}{r}h^{i}\rho^{\prime\prime}(|z|)\sum_{j}\nu^{j}(x_{a})\partial_{j}u^{i}\Big{)}dx+o(r^{n})\allowdisplaybreaks$
$\displaystyle=\sum_{a,i}(-1)^{a+1}\int_{B_{r}(x_{a})-E}\Bigl{(}-2\sum_{k}\Big{[}\partial_{k}h^{i}\rho^{\prime}(|z|)\frac{z_{k}}{|z|}+h^{i}\partial_{k}\big{(}\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\big{)}\Big{]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\frac{1}{r}h^{i}\rho^{\prime\prime}(|z|)\Big{)}\sum_{j}\nu^{j}(x_{a})\partial_{j}u^{i}\,dx+o(r^{n})\allowdisplaybreaks$
$\displaystyle=\sum_{a,i}(-1)^{a+1}\int_{B_{r}(x_{a})-E}\Bigl{(}-2\sum_{k}\Big{[}\partial_{k}h^{i}\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\Big{]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\frac{1}{r}h^{i}\rho^{\prime\prime}(|z|)\Big{)}\sum_{j}\nu^{j}(x_{a})\partial_{j}u^{i}\,dx+o(r^{n})\allowdisplaybreaks$
$\displaystyle=\sum_{a,i}(-1)^{a+1}r^{n}\int_{B_{1}\cap\\{|\mathbf{u}_{r}|>0\\}}\Bigl{(}-2\sum_{k}\Big{[}\partial_{k}h_{r}^{i}\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\Big{]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\frac{1}{r}rh_{r}^{i}\rho^{\prime\prime}(|z|)\Big{)}\sum_{j}\nu^{j}(x_{a})\partial_{j}u_{r}^{i}\,dz+o(r^{n}).$

Now note that $\partial_{j}u_{r}^{i}(z)\to
q^{i}(x_{a})\nu^{j}(x_{a})=\partial_{j}u^{i}(x_{a})$ when $z\cdot\nu(x_{a})>0$
by the results of [4]. Next note that $h^{i}$ is Lipschitz continuous, since
$u^{i}$ is Lipschitz and we have $0\leq h^{i}\leq cu^{i}$ for some constant
$c$. To see this note that the function
$\partial_{\xi_{i}}\Gamma(x,\partial_{\nu}\mathbf{u})$ is positive and
continuous on the compact set $\partial\Omega$, so it is bounded there, and
thus for some $c>0$ we have
$h^{i}=\partial_{\xi_{i}}\Gamma(x,\partial_{\nu}\mathbf{u})\leq
c\varphi^{i}=cu^{i}$ on $\partial\Omega$. Hence the claim follows by the
maximum principle. Therefore, by Lemma B.1 in [7], we also have
$\partial_{k}h_{r}^{i}(z)\to
p^{i}(x_{a})\nu^{k}(x_{a})=\partial_{k}h^{i}(x_{a})$ for some function
$p^{i}$, and $h_{r}^{i}(z)\to\nabla h^{i}(x_{a})\cdot z$ as $h^{i}(x_{a})=0$.
Thus if we divide the above expression by $r^{n}$ and let $r\to 0$ we obtain

$\displaystyle 0$
$\displaystyle\leq\sum_{a,i}(-1)^{a+1}\int_{B_{1}\cap\\{z\cdot\nu(x_{a})>0\\}}\Bigl{(}-2\sum_{k}\Big{[}\partial_{k}h^{i}(x_{a})\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\Big{]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\big{(}\nabla
h^{i}(x_{a})\cdot
z\big{)}\rho^{\prime\prime}(|z|)\Big{)}\sum_{j}\nu^{j}(x_{a})\partial_{j}u^{i}(x_{a})\,dz\allowdisplaybreaks$
$\displaystyle=\sum_{a,i}(-1)^{a+1}\int_{B_{1}\cap\\{z\cdot\nu(x_{a})>0\\}}\Bigl{(}-2\sum_{k}\Big{[}p^{i}(x_{a})\nu^{k}(x_{a})\rho^{\prime}(|z|)\frac{z_{k}}{|z|}\Big{]}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-p^{i}(x_{a})\big{(}\nu(x_{a})\cdot
z\big{)}\rho^{\prime\prime}(|z|)\Big{)}\partial_{\nu}u^{i}(x_{a})\,dz\allowdisplaybreaks$
$\displaystyle=\sum_{a,i}(-1)^{a}\int_{B_{1}\cap\\{z\cdot\nu(x_{a})>0\\}}\Bigl{(}\frac{2}{|z|}\rho^{\prime}(|z|)+\rho^{\prime\prime}(|z|)\Big{)}\big{(}\nu(x_{a})\cdot
z\big{)}p^{i}(x_{a})\partial_{\nu}u^{i}(x_{a})\,dz$
$\displaystyle=\sum_{a,i}(-1)^{a}\partial_{\nu}h^{i}(x_{a})\partial_{\nu}u^{i}(x_{a})\int_{B_{1}\cap\\{z\cdot\nu(x_{a})>0\\}}\Bigl{(}\frac{2}{|z|}\rho^{\prime}(|z|)+\rho^{\prime\prime}(|z|)\Big{)}\big{(}\nu(x_{a})\cdot
z\big{)}dz$
$\displaystyle=C_{\rho}\Big{(}\sum_{i}\partial_{\nu}h^{i}(x_{2})\partial_{\nu}u^{i}(x_{2})-\sum_{i}\partial_{\nu}h^{i}(x_{1})\partial_{\nu}u^{i}(x_{1})\Big{)},$

where
$C_{\rho}=\int_{B_{1}\cap\\{z\cdot\nu(x_{a})>0\\}}\bigl{(}\frac{2}{|z|}\rho^{\prime}(|z|)+\rho^{\prime\prime}(|z|)\big{)}\big{(}\nu(x_{a})\cdot
z\big{)}dz$ does not depend on $x_{a}$; we have also used the fact that
$p^{i}(x_{a})=\partial_{\nu}h^{i}(x_{a})$. By switching the role of
$x_{1},x_{2}$ we conclude that

$\sum_{i}\partial_{\nu}h^{i}(x_{2})\partial_{\nu}u^{i}(x_{2})-\sum_{i}\partial_{\nu}h^{i}(x_{1})\partial_{\nu}u^{i}(x_{1})$

must be zero, as desired. ∎

The main idea to show the regularity of the free boundary lies in utilizing
the boundary Harnack principle, which allows us to reduce the system into a
scalar problem. The key tool in employing this approach is non-tangential
accessibility of the domain; for the definition of non-tangentially accessible
(NTA) domains we refer to [2].

######  Lemma 5.2.

Let $\mathbf{u}$ be a solution of the minimization problem (1.1) for $p=2$.
Then $U=\\{x:|\mathbf{u}(x)|>0\\}$ is a non-tangentially accessible domain.

###### Proof.

This result follows from the same analysis as of Theorem 4.8 in [2] for the
function ${\scriptstyle\mathcal{U}}=u^{1}+\cdots+u^{m}$. Note that
${\scriptstyle\mathcal{U}}$ is harmonic in
$\\{|\mathbf{u}|>0\\}=\\{{\scriptstyle\mathcal{U}}>0\\}$ (these two sets are
equal due to Lemma 3.2), and the function ${\scriptstyle\mathcal{U}}$ is also
Lipschitz continuous and satisfies the nondegeneracy property by Corollary
3.3. ∎

######  Theorem 5.3.

Let $x_{0}\in\mathcal{R}$ be a regular point of the free boundary. Then there
is $r>0$ such that $B_{r}(x_{0})\cap\partial\\{|\mathbf{u}|>0\\}$ is a
$C^{1,\alpha}$ hypersurface for some $\alpha>0$.

###### Proof.

We may assume that $u^{1}>0$ in $B_{r_{0}}(x_{0})\cap\\{|\mathbf{u}|>0\\}$ for
some $r_{0}>0$. First we show that for some $0<r\leq r_{0}$ there is a Hölder
continuous function $g$ defined on
$B_{r}(x_{0})\cap\partial\\{|\mathbf{u}|>0\\}$, such that in the viscosity
sense we have

$\partial_{\nu}h^{1}\partial_{\nu}u^{1}=g\quad\text{ on }\partial E,$

where $h^{1}$ is defined in Lemma 5.1. Since
$B_{r_{0}}(x_{0})\cap\\{|\mathbf{u}|>0\\}$ is an NTA domain, the boundary
Harnack inequality implies that $G^{i}:=u^{i}/u^{1}$ and $H^{i}=h^{i}/h^{1}$
are Hölder continuous functions in
$B_{r}(x_{0})\cap\overline{\\{|\mathbf{u}|>0\\}}$ for some $0<r\leq r_{0}$.
Now if we consider a one-sided tangent ball at some point $y\in
B_{r}(x_{0})\cap\partial\\{|\mathbf{u}|>0\\}$, we have asymptotic developments
(see Lemma B.1 in [7], noting that $h^{i}$ is Lipschitz as we have shown in
the proof of Lemma 5.1)

$\begin{split}u^{i}(y+x)&=q^{i}(y)(x\cdot\nu(y))^{+}+o(|x|),\\\
h^{i}(y+x)&=p^{i}(y)(x\cdot\nu(y))^{+}+o(|x|).\end{split}$

Therefore $G^{i}(y)=q^{i}(y)/q^{1}(y)$ and $H^{i}(y)=p^{i}(y)/p^{1}(y)$. Thus
from (5.1) we can infer that

$p^{1}(y)q^{1}(y)\Big{(}1+\sum_{i>1}G^{i}(y)H^{i}(y)\Big{)}=\sum_{i}p^{i}(y)q^{i}(y)=\sum_{i}\partial_{\nu}h^{i}\partial_{\nu}u^{i}$

is constant for every $y\in B_{r}(x_{0})\cap\partial\\{|\mathbf{u}|>0\\}$.
Note that $G^{i},H^{i}>0$ at $y$ as $p^{i},q^{i}>0$. Hence by applying Theorem
3.1 in [10] we get the desired result. ∎

######  Corollary 5.4.

Let $\mathbf{u}$ be a solution of the minimization problem (1.1) for $p=2$.
Then the regular part of the free boundary, $\mathcal{R}$, is analytic.

###### Proof.

Suppose $0\in\mathcal{R}$ and $u^{1}>0$ in $B_{r}\cap\\{|\mathbf{u}|>0\\}$.
Then we apply the hodograph-Legendre transformation $x\mapsto
y=(x_{1},\dots,x_{n-1},u^{1})$. Next we define the partial Legendre functions

$\displaystyle v^{1}(y):=x_{n},\quad v^{i}(y):=u^{i}(x),$
$\displaystyle\quad\text{for }i=2,\cdots,m,$ $\displaystyle
w^{i}(y):=h^{i}(x),$ $\displaystyle\quad\text{for }i=1,\cdots,m.$

As $\mathcal{R}$ is $C^{1,\alpha}$, it follows that $u^{i}$ and $h^{i}$ are in
$C^{1,\alpha}(\overline{B_{r}\cap\\{|\mathbf{u}|>0\\}})$. So, $v^{i}$ and
$w^{i}$ are $C^{1,\alpha}$ in a neighborhood of the origin in $\\{y_{n}\geq
0\\}$. Now we have verified all the hypothesis of Theorem 7.1 in [2], and
through a similar argument we can obtain the analyticity of $\mathcal{R}$. ∎

## Acknowledgements

This project was initiated while the authors stayed at Institute Mittag
Leffler (Sweden), during the program Geometric aspects of nonlinear PDE. H.
Shahgholian was supported by Swedish research Council.

## Declarations

Data availability statement: All data needed are contained in the manuscript.

Funding and/or Conflicts of interests/Competing interests: The authors declare
that there are no financial, competing or conflict of interests.

## References

  * [1] N. E. Aguilera, H. W. Alt, and L. A. Caffarelli.  An optimization problem with volume constraint.  SIAM journal on control and optimization, 24(2):191–198, 1986.
  * [2] N. E. Aguilera, L. A. Caffarelli, and J. Spruck.  An optimization problem in heat conduction.  Annali della Scuola Normale Superiore di Pisa-Classe di Scienze, 14(3):355–387, 1987.
  * [3] H. W. Alt and L. A. Caffarelli.  Existence and regularity for a minimum problem with free boundary.  Journal für die reine und angewandte Mathematik, 325:105–144, 1981.
  * [4] H. W. Alt, L. A. Caffarelli, and A. Friedman.  A free boundary problem for quasi-linear elliptic equations.  Annali della Scuola Normale Superiore di Pisa-Classe di Scienze, 11(1):1–44, 1984.
  * [5] G. Anzellotti.  Pairings between measures and bounded functions and compensated compactness.  Annali di Matematica pura ed applicata, 135:293–318, 1983.
  * [6] D. Danielli and A. Petrosyan.  A minimum problem with free boundary for a degenerate quasilinear operator.  Calculus of Variations and Partial Differential Equations, 23(1):97–124, 2005.
  * [7] M. Fotouhi and H. Shahgholian.  A minimization problem with free boundary for $p$-laplacian weakly coupled system.  arXiv preprint arXiv:2301.02236, 2023.
  * [8] O. A. Ladyzhenskaya and N. N. Ural’tseva.  Linear and quasilinear elliptic equations.  Translated from the Russian by Scripta Technica, Inc. Translation editor: Leon Ehrenpreis. Academic Press, New York-London, 1968.
  * [9] J. L. Lewis.  Regularity of the derivatives of solutions to certain degenerate elliptic equations.  Indiana University Mathematics Journal, 32(6):849–858, 1983.
  * [10] F. P. Maiale, G. Tortone, and B. Velichkov.  Epsilon-regularity for the solutions of a free boundary system.  arXiv preprint arXiv:2108.03606, 2021.
  * [11] J. Malỳ and W. P. Ziemer.  Fine regularity of solutions of elliptic partial differential equations.  Number 51. American Mathematical Soc., 1997.
  * [12] E. V. Teixeira.  The nonlinear optimization problem in heat conduction.  Calculus of Variations and Partial Differential Equations, 24(1):21–46, 2005.