y}(ih,jh)\approx\begin{dcases}\frac{3U_{i,j}-4U_{i,j-1}+U_{i,j-2}}{2h}&\text{if
$n_{y}(ih,jh)\geq 0$}\\\ \frac{-U_{i,j+2}+4U_{i,j+1}-3U_{i,j}}{2h}&\text{if
$n_{y}(ih,jh)<0$}.\end{dcases}$ (104)

After introducing, as for the disk case, a vector
$X_{h}\in\mathbb{R}^{(N+1)^{2}}$ such that $U_{ij}=X_{j+(N+1)i}=X_{n}$, we
obtain the following linear system

$A_{h}X_{h}=B_{h},$ (105)

where

  * •

if $x_{n}=0$ or $1$, or $y_{n}=0$ or $1$ then
$A_{n,n}=1\text{ and }A_{n,m}=0\text{ for $m\neq n$},\text{ and }B_{n}=0\text{
(Dirichlet boundary conditions)}$

  * •

else the nonzero entries of the matrix $A_{h}$ are summarized in Table 3,

  * •

if $x_{n}\notin\\{0,1\\}$ and $x_{n}\notin\\{0,1\\}$ then
$B_{n}=(1-\chi_{n})f(x_{n},y_{n})+\frac{\chi_{n}}{\varepsilon}g(x_{n},y_{n})$,

using the same notations as for the disk case:

  * •

$x_{n}=ih$ and $y_{n}=jh$ where $n=(N+1)i+j$

  * •

$F_{n}=F(x_{n},y_{n})$ for a function
$F\colon{\mathbb{R}}^{2}\to{\mathbb{R}}$.

$x_{n}$ | $\geq 0.5+R$ | $\leq 0.5-R$ | $>0.5-R$ | else
---|---|---|---|---
$y_{n}$ | $>0.5-R$ | $<0.5+R$ | $\leq 0.5-R$ | else
$\implies(n_{x})_{n}$ | $\geq 0$ | $\leq 0$ | $\geq 0$ | $\leq 0$
$\implies(n_{y})_{n}$ | $\geq 0$ | $\leq 0$ | $\leq 0$ | $\geq 0$
$A_{n,n-2N-2}$ | $\dfrac{\chi_{n}(n_{x})_{n}}{2\varepsilon h}$ | $0$ | $\dfrac{\chi_{n}(n_{x})_{n}}{2\varepsilon h}$ | $0$
$A_{n,n-N-1}$ | $-\dfrac{1}{h^{2}}-\dfrac{2\chi_{n}(n_{x})_{n}}{\varepsilon h}$ | $-\dfrac{1}{h^{2}}$ | $-\dfrac{1}{h^{2}}-\dfrac{2\chi_{n}(n_{x})_{n}}{\varepsilon h}$ | $-\dfrac{1}{h^{2}}$
$A_{n,n-2}$ | $\dfrac{\chi_{n}(n_{y})_{n}}{2\varepsilon h}$ | $0$ | $0$ | $\dfrac{\chi_{n}(n_{y})_{n}}{2\varepsilon h}$
$A_{n,n-1}$ | $-\dfrac{1}{h^{2}}-\dfrac{2\chi_{n}(n_{y})_{n}}{\varepsilon h}$ | $-\dfrac{1}{h^{2}}$ | $-\dfrac{1}{h^{2}}$ | $-\frac{1}{h^{2}}-\dfrac{2\chi_{n}(n_{y})_{n}}{\varepsilon h}$
$A_{n,n}$ | $\dfrac{4}{h^{2}}+1+\dfrac{3\chi_{n}}{2\varepsilon h}(|(n_{x})_{n}|+|(n_{y})_{n}|)+\dfrac{\alpha\chi_{n}}{\varepsilon}$
$A_{n,n+1}$ | $-\dfrac{1}{h^{2}}$ | $-\dfrac{1}{h^{2}}+\dfrac{2\chi_{n}(n_{y})_{n}}{\varepsilon h}$ | $-\dfrac{1}{h^{2}}+\dfrac{2\chi_{n}(n_{y})_{n}}{\varepsilon h}$ | $-\frac{1}{h^{2}}$
$A_{n,n+2}$ | $0$ | $-\dfrac{\chi_{n}(n_{y})_{n}}{2\varepsilon h}$ | $-\dfrac{\chi_{n}(n_{y})_{n}}{2\varepsilon h}$ | $0$
$A_{n,n+N+1}$ | $-\dfrac{1}{h^{2}}$ | $-\dfrac{1}{h^{2}}+\dfrac{2\chi_{n}(n_{x})_{n}}{\varepsilon h}$ | $-\dfrac{1}{h^{2}}$ | $-\dfrac{1}{h^{2}}+\frac{2\chi_{n}(n_{x})_{n}}{\varepsilon h}$
$A_{n,n+2N+2}$ | $0$ | $-\dfrac{\chi_{n}(n_{x})_{n}}{2\varepsilon h}$ | $0$ | $-\dfrac{\chi_{n}(n_{x})_{n}}{2\varepsilon h}$
Table 3: Nonzero entries of the matrix $A_{h}$ of the finite difference scheme
(105).

##### Special treatment of discontinuities:

We recall that

$\chi(x,y)=\chi_{\mathcal{\omega}}(x,y)$ (106)

is the characteristic function of $\omega=\Omega\in\overline{\mathcal{U}}$ and
hence is a discontinuous function, and that $n$ and $g$ are extensions to
$\omega$ of the data on the boundary of the initial domain, and have also
discontinuities. Since we are using finite differences, we have to handle
these discontinuities by making suitable choices for the values at
discontinuity points. Following the numerical tests carried out in the one-
dimensional case [7] where second order convergence was recovered when we
impose $\chi$ to be $1$ on the boundary of the initial domain, we do the same
in the two-dimensional case by forcing $\chi$ to be $1$ on the boundary of
$\mathcal{U}$. Also, special treatment has to be done for the corners of the
initial square since the extensions $n$ and $g$ have discontinuities at these
points.

Consider for example the upper left corner $(x_{c},y_{c})=(0.5-R,0.5+R)$. At
this point, the continuous equation reads

$\displaystyle-\Delta
u_{\varepsilon}(x_{c},y_{c})+u_{\varepsilon}(x_{c},y_{c})+\frac{\chi(x_{c},y_{c})}{\varepsilon}\Bigl{(}\frac{\partial
u_{\varepsilon}}{\partial x}(x_{c},y_{c})n_{x}(x_{c},y_{c})+\frac{\partial
u_{\varepsilon}}{\partial y}(x_{c},y_{c})n_{y}(x_{c},y_{c})+\alpha
u_{\varepsilon}(x_{c},y_{c})\Bigr{)}$
$\displaystyle\quad=(1-\chi(x_{c},y_{c}))f(x_{c},y_{c})+\frac{\chi(x_{c},y_{c})}{\varepsilon}g(x_{c},y_{c})$

thus, if we impose $\chi(x_{c},y_{c})=1$,

$\displaystyle-\Delta
u_{\varepsilon}(x_{c},y_{c})+u_{\varepsilon}(x_{c},y_{c})+\frac{1}{\varepsilon}\Bigl{(}\frac{\partial
u_{\varepsilon}}{\partial x}(x_{c},y_{c})n_{x}(x_{c},y_{c})+\frac{\partial
u_{\varepsilon}}{\partial y}(x_{c},y_{c})n_{y}(x_{c},y_{c})+\alpha
u_{\varepsilon}(x_{c},y_{c})\Bigr{)}$
$\displaystyle\quad=\frac{1}{\varepsilon}g(x_{c},y_{c}).$

Since the aim of penalization is to recover the boundary conditions of the
initial problem for $\varepsilon$ small enough, we propose to choose $n_{x}$,
$n_{y}$ and $g$ at the corner so that the equality

$\frac{\partial u}{\partial x}(x_{c},y_{c})n_{x}(x_{c},y_{c})+\frac{\partial
u}{\partial y}(x_{c},y_{c})n_{y}(x_{c},y_{c})+\alpha
u(x_{c},y_{c})=g(x_{c},y_{c})$ (107)

is satisfied for the exact solution $u$. For sufficiently smooth $u\in
C^{2}([0.5-R,0.5+R]^{2})$ we have the following equations satisfied pointwise
on the sides of the square

$\begin{dcases}\frac{\partial u}{\partial
x}(x,y_{c})n_{x}(x,y_{c})+\frac{\partial u}{\partial
y}(x,y_{c})n_{y}(x,y_{c})+\alpha
u(x,y_{c})=g(x,y_{c})&x\in{]}x_{c},x_{c}+2R[\\\ \frac{\partial u}{\partial
x}(x_{c},y)n_{x}(x_{c},y)+\frac{\partial u}{\partial
y}(x_{c},y)n_{y}(x_{c},y)+\alpha
u(x_{c},y)=g(x_{c},y)&y\in{]}y_{c}-2R,y_{c}[.\end{dcases}$ (108)

Thus, it is sufficient for example to take

$\begin{dcases}n_{x}(x_{c},y_{c})=\lim\limits_{x\to
x_{c}^{+}}n_{x}(x,y_{c})\\\ n_{y}(x_{c},y_{c})=\lim\limits_{x\to
x_{c}^{+}}n_{y}(x,y_{c})\\\ g(x_{c},y_{c})=\lim\limits_{x\to
x_{c}^{+}}g(x,y_{c})\end{dcases}$ (109)

in order to satisfy (107). Other choices are also possible by taking instead
$\lim\limits_{y\to y_{c}^{-}}$, or even the sum of the two equations and then
the mean:

$\begin{dcases}n_{x}(x_{c},y_{c})=\frac{1}{2}\Bigl{(}\lim\limits_{x\to
x_{c}^{+}}n_{x}(x,y_{c})+\lim\limits_{y\to y_{c}^{-}}n_{x}(x_{c},y)\Bigr{)}\\\
n_{y}(x_{c},y_{c})=\frac{1}{2}\Bigl{(}\lim\limits_{x\to
x_{c}^{+}}n_{y}(x,y_{c})+\lim\limits_{y\to y_{c}^{-}}n_{y}(x_{c},y)\Bigr{)}\\\
g(x_{c},y_{c})=\frac{1}{2}\Bigl{(}\lim\limits_{x\to
x_{c}^{+}}g(x,y_{c})+\lim\limits_{y\to
y_{c}^{-}}g(x_{c},y)\Bigr{)}.\end{dcases}$ (110)

####  5.2.3 Numerical tests

We consider the following test example, from a manufactured solution:

  * •

$u(x,y)=\sin(5(x+y-1))$ the exact solution of the initial problem

  * •

$f$ and $\tilde{g}$ accordingly: $f=-\Delta u+u$ in $\mathcal{U}$,
$\tilde{g}=\nabla u\cdot\tilde{n}+\alpha u$ on $\partial\mathcal{U}$.

In the following tests, we will consider $\alpha=2$ and $R=0.3$.

#####  5.2.3.1 Error inside the fluid domain

We compare the numerical solution of the penalized problem using finite
differences to the exact solution of the initial problem inside the fluid
domain. The results are reported on Figures 8 and 9, and Tables 4 and 5.

We consider here the following errors:

  * •

$L^{2}$-norm error:

$\mathsf{Er_{L^{2}}}=\sum_{i=0}^{N}\sum_{j=0}^{N}h^{2}(u(ih,jh)-U_{ij})^{2}\times\mathds{1}_{C}(ih,jh),$

  * •

$L^{\infty}$-norm error:

$\mathsf{Er_{L^{\infty}}}=\max\limits_{0\leq i,j\leq
N}|u(ih,jh)-U_{ij}|\times\mathds{1}_{C}(ih,jh),$

  * •

$H^{1}$-norm error:

$\mathsf{Er_{H^{1}}}=\sum_{i=0}^{N}\sum_{j=0}^{N}h^{2}\Bigl{[}(u(ih,jh)-U_{ij})^{2}+\Bigl{(}\frac{\partial
u}{\partial x}(ih,jh)-U^{x}_{ij}\Bigr{)}^{2}+\Bigl{(}\frac{\partial
u}{\partial
y}(ih,jh)-U^{y}_{ij}\Bigr{)}^{2}\Bigr{]}\times\mathds{1}_{C}(ih,jh),$

where $U^{x}_{ij}$ and $U^{y}_{ij}$ are second-order finite difference
approximations of $u_{\varepsilon}(ih,jh)\approx U_{ij}$ (decentered at the
boundary of $\mathcal{U}$). In the tables, the errors correspond to
$C=[0.5-R,0.5+R]^{2}$, errors (no boundary) correspond to
$C=]0.5-R,0.5+R[^{2}$, while the interior errors correspond to
$C=[0.5-R+S,0.5+R-S]^{2}$ where $0<S<R$. Fox example, for $R=0.3$, we took
$S=0.1$.

(a) Numerical solution

(b) Exact solution

Figure 8: Comparison between the numerical solution of the penalized problem
and thee exact solution of the initial problem for $\varepsilon=10^{-10}$ and
$N=150$.

(a)

(b)

Figure 9: Square embedded in a square. Error (inside the fluid domain) between the numerical solution of the penalized problem and the exact solution of the initial problem with respect to: (a) $\varepsilon$ for $N=150$, (b) $N$ for $\varepsilon=10^{-10}$. $N$ | $L^{\infty}$ | $L^{2}$ | $H^{1}$ |  $L^{\infty}$ (no |  $L^{2}$ (no |  $H^{1}$ (no | $L^{\infty}$ | $L^{2}$ | $H^{1}$
---|---|---|---|---|---|---|---|---|---
|  |  |  | boundary) | boundary) | boundary) | (interior) | (interior) | (interior)
50 | – | – | – |  |  |  | – | – | –
100 | -1.969 | -1.995 | -2.010 | -1.908 | -1.899 | -1.939 | -1.872 | -1.873 | -1.919
150 | -1.982 | -1.997 | -2.005 | -1.949 | -1.943 | -1.965 | -1.843 | -1.849 | -1.833
200 | -1.987 | -1.998 | -2.003 | -1.965 | -1.960 | -1.975 | -2.072 | -2.063 | -2.132
Table 4: Square embedded in a square. Convergence order with respect to the mesh size $h=\frac{1}{N}$ for $\varepsilon=10^{-10}$. $\varepsilon$ | $L^{\infty}$ | $L^{2}$ | $H^{1}$ |  $L^{\infty}$ (no |  $L^{2}$ (no |  $H^{1}$ (no | $L^{\infty}$ | $L^{2}$ | $H^{1}$
---|---|---|---|---|---|---|---|---|---
|  |  |  | boundary) | boundary) | boundary) | (interior) | (interior) | (interior)
$10^{-1}$ | – | – | – | – | – | – | – | – | –
$10^{-2}$ | 0.862 | 0.793 | 0.915 | 0.857 | 0.789 | 0.912 | 0.749 | 0.766 | 0.825
$10^{-3}$ | 0.922 | 0.900 | 0.934 | 0.921 | 0.900 | 0.927 | 0.863 | 0.905 | 0.895
$10^{-4}$ | 1.071 | 1.208 | 0.869 | 1.082 | 1.216 | 0.868 | 1.325 | 1.304 | 0.938
Table 5: Square embedded in a square. Convergence order with respect to the
penalization parameter $\varepsilon$ for $N=150$.

The numerical results suggest as expected the convergence of $u_{\varepsilon}$
towards $u$ inside the fluid domain, with an order of convergence of $1$ with
respect to the penalization parameter as in the one-dimensional case [7], and
as obtained in the theoretical study assuming more regularity (Section 4). We
recover a second-order convergence with respect to the mesh size, if special
care is taken in the choice of the values at discontinuity points as explained
earlier (in addition to $N$ chosen so that the Cartesian grid contains the
points on the boundary of the initial square).

#####  5.2.3.2 Error inside the obstacle domain and investigation of the
boundary layer

In this section, we propose to study also the error inside the obstacle domain
as well as the presence of a boundary layer, that is expected theoretically
and observed numerically at the boundary $\partial\Omega$.

Indeed, a cross-section of the numerical solution of the penalized problem for
$y=0.5$ shows that $u_{\varepsilon}(\cdot,0.5)$ converges to
$u_{\lim}(\cdot,0.5)$ far from the boundary, with a boundary layer near the
boundary as the solution $u_{\varepsilon}$ has to satisfy homogeneous
Dirichlet boundary conditions on the boundary, see Figure 10.

Theoretically (Theorem 10), the expression of $u_{\lim}$ is the following in
$[0,1]\times]0.2,0.8[$

$u_{\lim}(x,y)=\begin{dcases}V^{0}(x,y)=u(x,y)&\text{if $x\in[0.2,0.8]$}\\\
\overline{W}^{0}(x,y)&\text{if $x\in[0,0.2]\cup[0.8,1]$},\end{dcases}$ (111)

and we obtain here explicitly (by solving (48)), again for $y\in{]}0.2,0.8[$,

$\overline{W}^{0}(x,y)=\begin{dcases}\Bigl{(}u(0.2,y)-\frac{\tilde{g}(0.2,y)}{\alpha}\Bigr{)}e^{-\alpha(0.2-x)}+\frac{\tilde{g}(0.2,y)}{\alpha}&\text{if
$x\in[0,0.2]$}\\\
\Bigl{(}u(0.8,y)-\frac{\tilde{g}(0.8,y)}{\alpha}\Bigr{)}e^{-\alpha(x-0.8)}+\frac{\tilde{g}(0.8,y)}{\alpha}&\text{if
$x\in[0.8,1]$}.\end{dcases}$ (112)

(a) Plots of $u_{\varepsilon}$ and $u_{\lim}$

(b) Plots of the ratio $\frac{u_{\varepsilon}}{u_{\lim}}$

Figure 10: Comparison between the numerical solution of the penalized problem
$u_{\varepsilon}$ for $\varepsilon\in\\{10^{-1},10^{-2},10^{-3},10^{-4}\\}$
and (red curve) the limit solution $u_{\lim}$ ($V^{0}=u$ in the fluid domain,
$\overline{W}^{0}$ in the obstacle domain) for $y=0.5$. Here, $N=150$.

In addition, Table 6 shows the convergence order of the error between the
numerical solution of the penalized problem and the limit solution inside the
obstacle domain far from the boundary: more precisely, we consider the error
inside the strip $[0.1,0.2]\times[0.5,0.5+h]$. An order of convergence of 1 is
recovered in the $H^{1}$-norm (as expected theoretically (cf. Theorem 10)) if
the mesh size is chosen small enough (we took $N=1000$, hence $h=10^{-3}$) in
order to witness this order for small values of $\varepsilon$.

$\varepsilon$ | $L^{\infty}$ | $L^{2}$ | $H^{1}$
---|---|---|---
$10^{-1}$ | – | – | –
$10^{-2}$ | 0.991 | 0.949 | 1.142
$10^{-3}$ | 0.929 | 0.915 | 1.012
$10^{-4}$ | 0.974 | 0.968 | 1.011
Table 6: Convergence order of $\|u_{\varepsilon}-\overline{W}^{0}\|$ in the
strip $[0.1,0.2]\times[0.5,0.5+h]$ with respect to the penalization parameter
$\varepsilon$ for $N=1000$.

In order to check the thickness of the boundary layer, which theoretically
should be of order $\varepsilon$, we recall that we expect, using (56) and
(59),

$u(x,0.5)=\overline{W}^{0}(x,0.5)+\theta(x)\widetilde{W}^{0}\Bigl{(}x,0.5,\frac{\varphi(x,0.5)}{\varepsilon}\Bigr{)}+O(\varepsilon),$
(113)

and we obtain here explicitly (by solving (49))

$\widetilde{W}^{0}\Bigl{(}x,0.5,\frac{\varphi(x,0.5)}{\varepsilon}\Bigr{)}=\begin{dcases}-\overline{W}^{0}(x,0.5)e^{-\frac{x}{\varepsilon}}&\text{if
$x\in[0,0.2]$}\\\
-\overline{W}^{0}(x,0.5)e^{-\frac{1-x}{\varepsilon}}&\text{if
$x\in[0.8,1]$}.\end{dcases}$ (114)

where $\theta$ is $C^{\infty}$, localized near the boundary $x=0$ or $x=1$ and
$\theta\equiv 1$ near this same boundary. Thus, near the boundary $x=0$ for
instance, we should have, for $\varepsilon$ small enough,

$u_{\varepsilon}(x,0.5)\approx u_{\lim}(x,0.5)(1-e^{-\frac{x}{\varepsilon}}).$
(115)

The thickness of the boundary layer near the boundary $x=0$ can be defined
(conventionally) by $E_{\varepsilon}$ such that

$\frac{u_{\varepsilon}(E_{\varepsilon},0.5)}{u_{\lim}(E_{\varepsilon},0.5)}=0.95\quad\implies\quad
1-e^{-\frac{E_{\varepsilon}}{\varepsilon}}\approx 0.95\quad\implies\quad-
E_{\varepsilon}\approx\log(0.05)\,\varepsilon.$ (116)

As estimating $E_{\varepsilon}$ numerically using its initial definition is
not easy given that the approximated values of $u_{\varepsilon}(x,0.5)$ are
discrete, and since we are more interested in the order of convergence of the
thickness of the boundary layer rather than the value of the thickness itself,
we suggest to estimate the thickness $BL_{\varepsilon}=\varepsilon$ instead
(this corresponds to the definition $E_{\varepsilon}$ such that
$\frac{u_{\varepsilon}(E_{\varepsilon},0.5)}{u_{\lim}(E_{\varepsilon},0.5)}\approx
1-e^{-1}\approx 0.63$). This can be obtained, if we denote
$RU_{\varepsilon}(x)\approx\frac{u_{\varepsilon}(x,0.5)}{u_{\lim}(x,0.5)}$, by

$BL^{1}_{\varepsilon}\approx\frac{1}{\frac{\partial RU_{\varepsilon}}{\partial
x}(0)}\approx\frac{h}{RU_{\varepsilon}(h)}$ (117)

or

$BL^{2}_{\varepsilon}\approx-\frac{h}{\log(1-RU_{\varepsilon}(h))}$ (118)

where $h=\frac{1}{N}$ is the meshsize. The results are reported on Table 7. We
recover the first-order convergence of the length of the boundary layer with
respect to $\epsilon$, at least for values of $\varepsilon\geq 10^{-3}$. For
smaller values of $\varepsilon$ ($10^{-4}$ for instance), we are limited by
the discretization step $h=10^{-3}$ which does not allow to capture properly
the length of the boundary layer, and one needs to refine the mesh further for
this purpose. In addition, it seems that the estimate $BL^{2}_{\varepsilon}$
is more accurate, which makes sense since the estimate $BL_{\varepsilon}^{1}$
contains an additional error with respect to $h$ related to the approximation
of the derivative.

$\varepsilon$ | Value of $BL^{1}_{\varepsilon}$ | Order of convergence | Value of $BL^{2}_{\varepsilon}$ | Order of convergence
---|---|---|---|---
$10^{-1}$ | $0.116$ | – | $0.116$ | –
$10^{-2}$ | $1.10\times 10^{-2}$ | 1.024 | $1.05\times 10^{-2}$ | 1.043
$10^{-3}$ | $1.79\times 10^{-3}$ | 0.789 | $1.22\times 10^{-3}$ | 0.933
$10^{-4}$ | $1.06\times 10^{-3}$ | 0.224 | $3.64\times 10^{-4}$ | 0.525
Table 7: Estimate and convergence order of the thickness of the boundary layer
of $u_{\varepsilon}(x,0.5)$ near the boundary $x=0$ with respect to the
penalization parameter $\varepsilon$ for $N=1000$.

##  6 Generalization

Using the same approach as above, the main result of this article stating the
convergence of the penalization method holds in higher dimensions $d\geq 2$,
assuming enough regularity on the domains and the data. More precisely,
Theorems 10 and 11 hold assuming the following regularity hypotheses:

Regularity hypotheses in dimension $d$:

  * •

$\Omega$ and $\mathcal{U}$ are supposed to be of class
$C^{k+3\lfloor\frac{d}{2}\rfloor+5,1}$ in dimension $d$ for some $k\geq 2$

  * •

$f\in H^{k+3\lfloor\frac{d}{2}\rfloor+4}(\mathcal{U})$

  * •

$\tilde{g}\in
H^{k+3\lfloor\frac{d}{2}\rfloor+4+\frac{1}{2}}(\partial\mathcal{U})$

  * •

$n=\nabla\psi$ with

    * –

$\psi\in C^{k+2\lfloor\frac{d}{2}\rfloor+5}(\overline{\omega})$

    * –

$\psi=0$ on $\partial\mathcal{U}$

    * –

$\psi>0$ in $\omega$

    * –

$|\nabla\psi|\geq\lambda>0$ in $\omega$

    * –

$n\cdot\nu_{\Omega}\geq\lambda_{0}>0$ on $\partial\Omega$, where
$\nu_{\Omega}$ is the outward normal vector to $\Omega$

##  7 Conclusion

In this paper, we have investigated the extension to higher dimensions of a
previously suggested penalization method for Neumann or Robin boundary
conditions [7]. Our results include the existence and uniqueness of the
solution of the penalized problem and the convergence with respect to the
penalization parameter towards the solution. The existence and uniqueness have
been established using the method of Droniou for non-coercive linear elliptic
problems [14]. To prove the convergence, we developed a boundary layer
approach, assuming enough regularity, for Neumann or Robin boundary
conditions, adapted from the Dirichlet approach [9]. The primary challenge in
applying the boundary layer approach to Neumann or Robin boundary conditions
is that, in contrast to Dirichlet boundary conditions, the penalized problem
is not coercive, making it more challenging to derive essential estimates on
the remainders. To this aim, in order to derive such estimates, a thorough
examination of the dual problem turned out to be pertinent: the convergence
result was reduced to finding suitable supersolutions to the dual problem.
While the construction of these supersolutions might be explicit in special
cases: one-dimensional case, spherical symmetric case, it is more tricky in
general settings. We proposed to conduct an additional formal boundary layer
approach on the dual problem in order to construct theses supersolutions as
approximated solutions of the dual problem. Finally, we have derived numerical
experiments in two dimensions, using adequate finite difference schemes
suitable for advection-dominated problems. Even if finite differences are not
the best choice of discretization given the discontinuity of the coefficients
introduced by the penalization, these numerical simulations allowed to
illustrate the theoretical results and also suggest that the convergence
result may be valid in less regular settings. It would thus be interesting to
weaken the regularity assumptions (ours are not necessarily sharp) in the
context of the boundary layer approach or using weak convergence approaches.

Perspectives of this work include different directions:

  * •

extension of the proposed method to other model problems of interest in fluid
mechanics: Stokes, Navier-Stokes equations, …as it has been performed for
penalized Dirichlet boundary conditions [9],

  * •

practical extension of the various data on the boundary to the complementary
domain,

  * •

improvement of numerical approximations using stabilized finite elements
suitable for singular perturbation problems [19, 24], more suitable for our
penalized problem,

  * •

application to problems raised on moving domains such as those associated with
the simulation of population dynamics under climate change [8, 28].

## Acknowledgments

B. Bensiali thanks Centrale Marseille/I2M for the hospitality and assistance
during her research stays. She also thanks Prof. Stefan Sauter for useful
discussions on this work at the 11th Zurich Summer School 2021 and for the
invitation to the “Institut für Mathematik” at the University of Zurich, which
indirectly allowed her to treat the spherical symmetric case using Bessel
functions (Appendix C.2).

This research was supported through computational resources of HPC-MARWAN
(hpc.marwan.ma) provided by the National Center for Scientific and Technical
Research (CNRST), Rabat, Morocco.

## References

  * [1] M Abramovitz and IA Stegun.  Handbook of mathematical functions, vol. 55.  National Bureau of Standards Applied Mathematics Series, Washington, DC, 1972.
  * [2] Robert A Adams and John JF Fournier.  Sobolev spaces.  Elsevier, 2003.
  * [3] G. Allaire and A. ERN.  Optimisation et contrôle.  Cours de l’Ecole Polytechnique. 2023.
  * [4] P. Angot, C.H. Bruneau, and P. Fabrie.  A penalization method to take into account obstacles in incompressible viscous flows.  Numer. Math., 81:497–520, 1999.
  * [5] Philippe Angot.  Analysis of singular perturbations on the brinkman problem for fictitious domain models of viscous flows.  Mathematical methods in the applied sciences, 22(16):1395–1412, 1999.
  * [6] Claude Bardos.  Problèmes aux limites pour les équations aux dérivées partielles du premier ordre à coefficients réels; théorèmes d’approximation; application à l’équation de transport.  In Annales scientifiques de l’École normale supérieure, volume 3, pages 185–233, 1970.
  * [7] B. Bensiali, G. Chiavassa, and J. Liandrat.  Penalization of Robin boundary conditions.  Appl. Numer. Math., 96:134–152, 2015.
  * [8] Henri Berestycki, Odo Diekmann, Cornelis J Nagelkerke, and Paul A Zegeling.  Can a species keep pace with a shifting climate?  Bulletin of mathematical biology, 71(2):399–429, 2009.
  * [9] F. Boyer and P. Fabrie.  Mathematical tools for the study of incompressible Navier-Stokes equations and related models, volume 183 of Applied Mathematical Sciences.  Springer, 2010.
  * [10] H. Brezis.  Functional analysis, Sobolev spaces and partial differential equations.  Universitext. Springer, 2010.
  * [11] P.G. Ciarlet.  Introduction à l’analyse numérique matricielle et à l’optimisation.  Collection Mathématiques appliquées pour la maîtrise. Dunod, 1998.
  * [12] Anne-Laure Dalibard and Laure Saint-Raymond.  Mathematical study of degenerate boundary layers: A large scale ocean circulation problem, volume 253.  American Mathematical Society, 2018.
  * [13] K. Deimling.  Nonlinear functional analysis (springer, berlin, 1985).
  * [14] J. Droniou.  Non-coercive linear elliptic problems.  Potential Analysis, 2002.
  * [15] J. Droniou and C. Imbert.  Solutions de viscosité et solutions variationnelles pour EDP non-linéaires.  Cours de DEA, Département de Mathématiques, Montpellier II, 2002\.
  * [16] Alexandre Ern and Jean-Luc Guermond.  Finite Elements III: First-Order and Time-Dependent PDEs, volume 74.  Springer Nature, 2021.
  * [17] Lawrence C Evans.  Partial differential equations, volume 19.  American Mathematical Society, 2022.
  * [18] L.C. Evans.  Partial differential equations.  American Mathematical Society, RI. second edition 2010.
  * [19] Klaus Gerdes, Jens Markus Melenk, Christoph Schwab, and Dominik Schötzau.  The hp-version of the streamline diffusion finite element method in two space dimensions.  Mathematical Models and Methods in Applied Sciences, 11(02):301–337, 2001.
  * [20] Gung-Min Gie, Makram Hamouda, Chang-Yeol Jung, and Roger M Temam.  Singular perturbations and boundary layers.  Springer, 2018.
  * [21] Philip Hartman.  Ordinary differential equations.  SIAM, 2002.
  * [22] H. Le Dret.  Fixed point theorems and applications.  In Nonlinear Elliptic Partial Differential Equations, pages 47–68. Springer, 2018.
  * [23] Jacques-Louis Lions.  Perturbations singulières dans les problèmes aux limites et en contrôle optimal, volume 323.  Springer, 2006.
  * [24] Jens M Melenk.  hp-Finite element methods for singular perturbations.  Number 1796. Springer Science & Business Media, 2002.
  * [25] R. Mittal and G. Iaccarino.  Immersed boundary methods.  Annu. Rev. Fluid Mech., 37:239–261, 2005.
  * [26] C. Peskin.  The immersed boundary method.  Acta Numerica, 11:479–517, 2002.
  * [27] I. Ramiere, P. Angot, and M. Belliard.  A fictitious domain approach with spread interface for elliptic problems with general boundary conditions.  Comput. Methods Appl. Mech. Eng., 196:766–781, 2007.
  * [28] Lionel Roques, Alain Roques, Henri Berestycki, and André Kretzschmar.  A population facing climate change: joint influences of Allee effects and environmental boundary geometry.  Population Ecology, 50(2):215–225, 2008.
  * [29] François Rouvière.  Petit guide de calcul différentiel: à l’usage de la licence et de l’agrégation.  Cassini, 2009.
  * [30] Teluo Sakurai, Katsunori Yoshimatsu, Naoya Okamoto, and Kai Schneider.  Volume penalization for inhomogeneous neumann boundary conditions modeling scalar flux in complicated geometry.  Journal of Computational Physics, 390:452–469, 2019.
  * [31] Stefan Sauter and Céline Torres.  The heterogeneous Helmholtz problem with spherical symmetry: Green’s operator and stability estimates.  Asymptotic Analysis, 125(3-4):289–325, 2021.
  * [32] K. Schneider.  Immersed boundary methods for numerical simulation of confined fluid and plasma turbulence in complex geometries: a review.  J. Plasma Phys., 6, 2015.
  * [33] Ramakrishnan Thirumalaisamy, Neelesh A Patankar, and Amneet Pal Singh Bhalla.  Handling Neumann and Robin boundary conditions in a fictitious domain volume penalization framework.  Journal of Computational Physics, 448:110726, 2022.

##  Appendix A Regularity for advection-reaction equation using the method of
characteristics

We give here the proof of Theorem 8.

###### Proof.

We use the method of characteristics to solve (12), starting by considering
the constant case $\beta(x)=\alpha$, and follow the following steps.

##### Step 1. Method of characteristics

For simplicity, we assume $\psi$ defined in all
${\mathbb{R}}^{2}\setminus\mathcal{U}$ and embed the domain $\omega$ inside
the domain delimited by $\psi=0$ ($\partial\mathcal{U}$) and
$\psi=c_{2}>\max_{\overline{\omega}}\psi$. The previous set is not empty (see
below). In addition, the assumption $\psi\in
C^{k+1}({\mathbb{R}}^{2}\setminus\mathcal{U})$ and $|\nabla\psi|\geq\lambda>0$
imply $\frac{\partial\psi}{\partial x}(x,y)\neq 0$ or
$\frac{\partial\psi}{\partial y}(x,y)\neq 0$, thus using the theorem of
implicit functions, we deduce that $\psi=c$ is locally a curve of class
$C^{k+1}$. For simplisity, we will also assume that the curves $\psi=c$ are
closed.. For $(x,y)\in\omega$, we introduce

$\begin{dcases}\dot{X}(s)=n(X(s))=\nabla\psi(X(s))\\\ X(t)=(x,y).\end{dcases}$
(119)

The conditions on $\nabla\psi$ (continuous and globally Lipschitzian) enable
to use the global Cauchy-Lipschitz theorem [29] that ensures the existence and
uniqueness in $C^{1}$ of $s\mapsto X(s)$ for all $s$, then denoting
$\Psi(s)=\psi(X(s))$, one has

$\displaystyle\frac{d\Psi}{ds}(s)$ $\displaystyle=\frac{\partial\psi}{\partial
x}(X(s))\,\dot{X_{1}}(s)+\frac{\partial\psi}{\partial
y}(X(s))\,\dot{X_{2}}(s)$ $\displaystyle=|\nabla\psi(X(s))|^{2}.$

Using the hypothesis $|\nabla\psi|\geq\lambda>0$ and integrating we get for
$t\leq T$, $\Psi(T)\geq\Psi(t)+(T-t)\lambda^{2}$ and, for $t>T$,
$\Psi(T)\leq\Psi(t)-(t-T)\lambda^{2}.$

Thus for all $(x,y)\in\omega$ (embeded between $\psi=0$ and $\psi=c_{2}$),
there exists a unique $C^{1}$ curve $s\mapsto X(s)=(X_{1}(s),X_{2}(s))$
solution to (119) that passes through $(x,y)$ and cuts the curve $\psi=0$
($\partial\mathcal{U}$) and the curve $\psi=c_{2}$. We deduce that for all
$(x,y)\in\omega$ there exists a unique $C^{1}$ curve $s\mapsto
X(s)=(X_{1}(s),X_{2}(s))$ solution to (119), which cuts the curve $\psi=0$
($\partial\mathcal{U}$) and $\partial\Omega$.

We now introduce

$\begin{dcases}\frac{\partial X}{\partial
t}(\xi,t)=n(X(\xi,t))=\nabla\psi(X(\xi,t))\\\
X(\xi,0)=\gamma(\xi)\in\partial\mathcal{U}.\end{dcases}$ (120)

where $t\geq 0$ and $\xi\in[0,\Xi]$ is a parameterization of
$\partial\mathcal{U}$ ($\psi=0$), i.e. $\partial\mathcal{U}=\\{\gamma(\xi),\
\xi\in[0,\Xi]\\}$. Given the assumption on $\partial\mathcal{U}$, $\gamma\in
C^{k}([0,\Xi])$. We know, from the previous reasoning, that for all
$\xi\in[0,\Xi]$ the curve $t\mapsto X(\xi,t)$ will cut the curve
$\psi=c_{2}>0$ at a time $T_{2,\xi}$ such that

$\displaystyle c_{2}=\Psi(T_{2,\xi})$
$\displaystyle\geq\Psi(0)+T_{2,\xi}\,\lambda^{2}$ $\displaystyle T_{2,\xi}$
$\displaystyle\leq\frac{c_{2}}{\lambda^{2}}.$

Thus there exists a time $T_{f}<\infty$ (since
$T_{f}\leq\frac{c_{2}}{\lambda^{2}}$) such that all the characteristics
starting from $\partial\mathcal{U}$ (at $t=0$) arrive at the curve
$(\psi=c_{2})$. We deduce that for all $\xi\in[0,\Xi]$, there exists a time
$T_{\xi}<\infty$ such that the curve $t\mapsto X(\xi,t)$ arrives at
$\partial\Omega$. Since we assume $n\cdot\nu_{\Omega}>0$, the curve $t\mapsto
X(\xi,t)$ cannot come back inside of $\omega$ after reaching it for the first
time. We thus know that by solving (120) for $\xi\in[0,\Xi]$ and
$t\in[0,T_{\xi}]$, the union of all the solutions will fill exactly
$\overline{\omega}$.

Then, by denoting

$\displaystyle\mathcal{G}(\xi,t)$ $\displaystyle=g(X(\xi,t))$
$\displaystyle\mathcal{W}(\xi,t)$ $\displaystyle=W(X(\xi,t))$
$\displaystyle\mathcal{V}(\xi)$ $\displaystyle=V(X(\xi,0))=V(\gamma(\xi)),$

we obtain from the equation on $W$ (12) that $\mathcal{W}(\xi,t)$ satisfies

$\displaystyle\dot{\mathcal{W}}(\xi,t)$ $\displaystyle=\frac{\partial
W}{\partial x}(X(\xi,t))\,\dot{X_{1}}(\xi,t)+\frac{\partial W}{\partial
y}(X(\xi,t))\,\dot{X_{2}}(\xi,t)$ $\displaystyle=\frac{\partial W}{\partial
x}(X(\xi,t))\,n_{1}(X(\xi,t))+\frac{\partial W}{\partial
y}(X(\xi,t))\,n_{2}(X(\xi,t))$ $\displaystyle=\nabla W(X(\xi,t))\cdot
n(X(\xi,t))$ $\displaystyle=g(X(\xi,t))-\alpha W(X(\xi,t))$
$\displaystyle=\mathcal{G}(\xi,t)-\alpha\mathcal{W}(\xi,t)$
$\displaystyle\dot{\mathcal{W}}(\xi,t)+\alpha\mathcal{W}(\xi,t)$
$\displaystyle=\mathcal{G}(\xi,t),$ (121)

whose solution is given by

$\mathcal{W}(\xi,t)=k(\xi)e^{-\alpha t}+\mathcal{W}_{\mathrm{part}}(\xi,t),$
(122)

where $\mathcal{W}_{\mathrm{part}}(\xi,t)$ is a particular solution that we
can obtain using the method of variation of constants: we look for a
particular solution of the form

$\mathcal{W}_{\mathrm{part}}(\xi,t)=k(\xi,t)e^{-\alpha t}$ (123)

then

$\displaystyle\dot{\mathcal{W}}_{\mathrm{part}}(\xi,t)$
$\displaystyle=\dot{k}(\xi,t)e^{-\alpha t}-k(\xi,t)\alpha e^{-\alpha t}$ (124)
$\displaystyle=\dot{k}(\xi,t)e^{-\alpha
t}-\alpha\mathcal{W}_{\mathrm{part}}(\xi,t).$ (125)

We want $\mathcal{W}_{\mathrm{part}}(\xi,t)$ to be a solution to (A), i.e.

$\displaystyle\dot{k}(\xi,t)=\mathcal{G}(\xi,t)e^{\alpha t}.$ (126)

We take

$k(\xi,t)=\int_{0}^{t}\mathcal{G}(\xi,s)e^{\alpha s}\,ds.$ (127)

Finally,

$\mathcal{W}(\xi,t)=(k(\xi)+k(\xi,t))e^{-\alpha t},$ (128)

and using the boundary condition in (12),
$\mathcal{W}(\xi,0)=W(X(\xi,0))=W(\gamma(\xi))=V(\gamma(\xi))=\mathcal{V}(\xi)$
and we obtain

$\mathcal{W}(\xi,t)=\Bigl{(}\mathcal{V}(\xi)+\int_{0}^{t}\mathcal{G}(\xi,s)e^{\alpha
s}\,ds\Bigr{)}e^{-\alpha t},$ (129)

thus

$W(X(\xi,t))=\Bigl{(}V(\gamma(\xi))+\int_{0}^{t}g(X(\xi,s))e^{\alpha
s}\,ds\Bigr{)}e^{-\alpha t}.$ (130)

Because the characteristics fill all $\overline{\omega}$, the previous formula
allows formally to define $W(x,y)$ at each point $(x,y)\in\overline{\omega}$.
Since we are interested in the regularity of the solution $W$ as a function of
the space variables, we need to express $W(x,y)$.

##### Step 2. Regularity of $X(\xi,t)$ with respect to $\xi$ and $t$

We denote by $G$ the function

$G\colon(\xi,t)\mapsto X(\xi,t)$ (131)

where for every $\xi\in[0,\Xi]$, $t\mapsto X(\xi,t)$ is the unique solution to
the Cauchy problem (120).

The regularity of $G$ with respect to $\xi$ is related to the differentiable
dependence on initial conditions for a Cauchy problem. Indeed, let’s introduce

$\begin{dcases}\frac{\partial Z}{\partial
t}(z_{0},t)=n(Z(z_{0},t))=\nabla\psi(Z(z_{0},t))\\\
Z(z_{0},0)=z_{0}\in{\mathbb{R}}^{2}\setminus\mathcal{U},\end{dcases}$ (132)

then since $\psi\in C^{k+1}(\overline{\omega})$, hence $n\in
C^{k}(\overline{\omega})$, we deduce that the solution $(z_{0},t)\mapsto
Z(z_{0},t)$ is $C^{k}$ on its domain of existence ([21], Corollary 4.1).

Using the fact that $G(\xi,t)=X(\xi,t)=Z(\gamma(\xi),t)$, we deduce, since
$\gamma\in C^{k}([0,\Xi])$, that $G\in C^{k}([0,\Xi]\times{\mathbb{R}}^{+})$.

##### Step 3. From local to global inverse function theorem

Here we show that $G$ is invertible, and that $H=G^{-1}$ is
$C^{k}(\overline{\omega})$.

$H\colon(x,y)\mapsto(\xi,t)\text{ such that }G(\xi,t)=(x,y).$ (133)

We make use of the inverse function theorem.

######  Theorem 12 (Inverse function theorem).

Let $f\colon U\to V$ be a map between open subsets of ${\mathbb{R}}^{n}$.
Assume $f$ is $C^{k}$. If $f$ is injective on a closed subset $A\subset U$ and
if the jacobian matrix of $f$ is invertible at each point of $A$, then $f$ is
injective in a neighborhood $A^{\prime}$ of $A$ and $f^{-1}\colon
f(A^{\prime})\to A^{\prime}$ is $C^{k}$.

We fix $\xi_{0}\in[0,\Xi]$. By definition $G(\xi,t)=X(\xi,t)$. Therefore, the
Jacobian of $G$ is given by

$\displaystyle JG(\xi,t)$
$\displaystyle=\begin{vmatrix}\frac{\partial{X_{1}}}{\partial\xi}(\xi,t)&\frac{\partial{X_{1}}}{\partial
t}(\xi,t)\\\
\frac{\partial{X_{2}}}{\partial\xi}(\xi,t)&\frac{\partial{X_{2}}}{\partial
t}(\xi,t)\end{vmatrix}$
$\displaystyle=\frac{\partial{X_{1}}}{\partial\xi}\frac{\partial{X_{2}}}{\partial
t}-\frac{\partial{X_{1}}}{\partial t}\frac{\partial{X_{2}}}{\partial\xi}.$

At the point $(\xi_{0},0)$,

$\displaystyle JG(\xi_{0},0)$
$\displaystyle=\frac{\partial{X_{1}}}{\partial\xi}(\xi_{0},0)\frac{\partial{X_{2}}}{\partial
t}(\xi_{0},0)-\frac{\partial{X_{1}}}{\partial
t}(\xi_{0},0)\frac{\partial{X_{2}}}{\partial\xi}(\xi_{0},0)$
$\displaystyle=\gamma_{1}^{\prime}(\xi_{0})n_{2}(\gamma(\xi_{0}))-n_{1}(\gamma(\xi_{0}))\gamma_{2}^{\prime}(\xi_{0})$
$\displaystyle=n(\gamma(\xi_{0}))\cdot(\gamma_{1}^{\prime}(\xi_{0}),-\gamma_{2}^{\prime}(\xi_{0})),$

since, by definition, $X(\xi,t)$ satisfies the characteristic ODE (120).
Forall $\xi_{0}\in[0,\Xi]$, $JG(\xi_{0},0)\neq 0$ since $n$ is not tangent to
$\partial\mathcal{U}$ (we say that $\partial\mathcal{U}$ is
noncharacteristic). We already know that $G$ is injective on
$[0,\Xi]\times{\mathbb{R}}^{+}$ (by the global Cauchy-Lipschitz theorem). By
the inverse function theorem, $G^{-1}\colon G(A^{\prime})\to A^{\prime}$ is
$C^{k}$ where $A^{\prime}$ is a neighborhood of $[0,\Xi]\times\\{0\\}$. We
thus obtain that $G^{-1}$ is $C^{k}$ on a neighborhood $B$ of
$\partial\mathcal{U}$.

Our goal is to establish a global differentiability, that is $G^{-1}\in
C^{k}(\overline{\omega})$. To this end, we will use a step by sted advancing
procedure. Let us explain the first step. We denote by $G_{I}^{-1}\colon
B_{I}\mapsto[0,\Xi]\times{\mathbb{R}}^{+}$ the restriction of $G^{-1}$ to
$B_{I}$ a neighborhood of $\partial\mathcal{U}$. We denote by
$c=\min_{\partial B_{I}\setminus\partial\mathcal{U}}\psi$, then we take
$\Gamma_{II}$ to be defined by $\psi=c_{II}$ where $0>c_{II}<c$. We recall
that we assumed that the curves $\psi=c$ are closed. $\Gamma_{II}\subset
B_{I}$ and we consider $\Gamma_{II}=\\{\gamma_{II}(\xi_{II}),\
\xi_{II}\in[0,\Xi_{II}]\\}$ with $\gamma_{II}\in C^{k}$ (this follows from the
regularity of $\psi$, the assumption on $\nabla\psi$ and using the theorem of
implicit functions, as stated previously). We introduce
$G_{II}(\xi_{II},t)=Y(\xi_{II},t)$, where $Y(\xi_{II},t)$ is the unique
solution to

$\begin{dcases}\frac{\partial Y}{\partial
t}(\xi_{II},t)=n(Y(\xi_{II},t))=\nabla\psi(Y(\xi_{II},t))\\\
Y(\xi_{II},0)=\gamma_{II}(\xi_{II})\in\Gamma_{II}.\end{dcases}$ (134)

Since $\Gamma_{II}$ is noncharacteristic ($n=\nabla\psi$ is normal to
$\Gamma_{II}$), we obtain as previously using the inverse function theorem
that $G_{II}^{-1}$ is $C^{k}$ on a neighborhood $B_{II}$ of $\Gamma_{II}$. For
$(x,y)\in B_{II}\setminus B_{I}$, since

$\displaystyle G_{II}^{-1}(x,y)=(\xi_{II},t_{II})\text{ such that
$G_{II}(\xi_{II},t_{II})=(x,y)$}$

and

$\displaystyle G_{I}^{-1}(\gamma_{II}(\xi_{II}))=(\xi,t)\text{ such that
$G_{I}(\xi,t)=\gamma_{II}(\xi_{II})$},$

we obtain

$\displaystyle
G^{-1}(x,y)=\Bigl{(}[G^{-1}_{I}\bigl{(}\gamma_{II}([G_{II}^{-1}(x,y)]_{1})\bigr{)}]_{1},[G^{-1}_{II}(x,y)]_{2}+[G^{-1}_{I}\bigl{(}\gamma_{II}([G_{II}^{-1}(x,y)]_{1})\bigr{)}]_{2}\Bigr{)}$
(135)

which allows to show by composition, that $G^{-1}$ is $C^{k}$ on $B_{I}\cup
B_{II}$.

It remains to show that this process allows to visit all $\overline{\omega}$.
Assume, by contradiction, that there exist unreachable curves
$\Gamma_{\hat{\beta}}=\\{x,\ \psi(x)=\hat{\beta}\\}$ and let $\beta_{*}$ be
the infimum of the $\hat{\beta}$ and note $\gamma_{*}$ the parameterization of
$\Gamma_{\beta_{*}}$. If $\Gamma_{*}$ can be reached, the procedure above
provides and extension and then a contradiction. If $\Gamma_{*}$ can not be
reached, we introduce $\tilde{G}_{*}(\xi_{*},t)=Z(\xi_{II},t)$, where
$Z(\xi_{II},t)$ is the unique solution to the backward Cauchy problem

$\begin{dcases}\frac{\partial Z}{\partial
t}(\xi_{*},t)=-n(Z(\xi_{*},t))=-\nabla\psi(Z(\xi_{*},t))\\\
Z(\xi_{*},0)=\gamma_{*}(\xi_{*})\in\Gamma_{*},\end{dcases}$ (136)

then $\tilde{G}_{*}$ is $C^{k}$ on its domain of existence.

Let us fix $t_{2}>0$ small enough such that, for all $\xi_{*}\in[0,\Xi_{*}]$,
$\tilde{G}_{*}(\xi_{*},t_{2})\in B^{*}$, where $B_{*}$ is a neighborhood of
$\Gamma_{*}$, then
$\tilde{G}_{*}(\xi_{*},t_{2})=\gamma_{\beta^{\prime}}(\xi_{\beta^{\prime}})$
($\tilde{G}_{*}(\xi_{*},t_{*})$ belongs to a curve $(\psi=\beta^{\prime})$
included in $B_{*}$ with $\beta^{\prime}<\beta_{*}$). Since we know that
$G^{-1}$ is $C^{k}$ on $\\{x,\ \psi(x)=\beta^{\prime}\\}$, and

$\displaystyle G^{-1}(\gamma_{*}(\xi_{*}))$ $\displaystyle=(\xi,t)\text{ such
that $G(\xi,t)=\gamma_{*}(\xi_{*})$}$
$\displaystyle\tilde{G}_{*}(\xi_{*},t_{2})$
$\displaystyle=\gamma_{\beta^{\prime}}(\xi_{\beta^{\prime}})$ $\displaystyle
G^{-1}(\gamma_{\beta^{\prime}}(\xi_{\beta^{\prime}}))$
$\displaystyle=(\xi,t_{1})\text{ such that
$G(\xi,t_{1})=\gamma_{\beta^{\prime}}(\xi_{\beta^{\prime}})$},$

we can extend $G^{-1}$ up to $\Gamma_{*}$ forming

$\displaystyle
G^{-1}(\gamma_{*}(\xi_{*}))=\Bigl{(}[G^{-1}\bigl{(}\tilde{G}_{*}(\xi_{*},t_{2})\bigr{)}]_{1},[G^{-1}\bigl{(}\tilde{G}_{*}(\xi_{*},t_{2})\bigr{)}]_{2}+t_{2}\Bigr{)}$
(137)

which is $C^{k}$ ($\gamma_{*}^{-1}$ is $C^{k}$) and provides the
contradiction.

##### Step 4. Regularity of $W$

We are now able to discuss the regularity of $W$. By inversion, (130) writes

$W(x,y)=\Bigl{(}V(\gamma(H_{1}(x,y)))+\int_{0}^{H_{2}(x,y)}g(G(H_{1}(x,y),s))e^{\alpha
s}\,ds\Bigr{)}e^{-\alpha H_{2}(x,y)},$ (138)

where $G$ and $H$ are given by (131) and (133) respectively.

Since $V$ is $C^{k}(\partial\mathcal{U})$, $\gamma\in C^{k}([0,\Xi])$ and
$H_{1}\in C^{k}(\overline{\omega})$, we deduce $(x,y)\mapsto
V(\gamma(H_{1}(x,y)))$ is $C^{k}(\overline{\omega})$. Similarly $(x,y)\mapsto
e^{-\alpha H_{2}(x,y)}$ is $C^{k}(\overline{\omega})$ by composition.

It then remains to study the regularity of the function

$(x,y)\mapsto L(x,y)=\int_{0}^{H_{2}(x,y)}g(G(H_{1}(x,y),s))e^{\alpha s}\,ds.$
(139)

Let us introduce

$(x,y,t)\mapsto K(x,y,t)=\int_{0}^{t}g(G(H_{1}(x,y),s))e^{\alpha
s}\,ds=\int_{0}^{t}h(x,y,s)\,ds.$ (140)

then

$L(x,y)=K(x,y,H_{2}(x,y)).$ (141)

Since $g\in C^{k}(\overline{\omega})$, $G\in
C^{k}([0,\Xi]\times{\mathbb{R}}^{+})$ and $H\in C^{k}(\overline{\omega})$ then
$(x,y,s)\mapsto h(x,y,s)$ is $C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$.

We prove by induction on $k\geq 1$ that if $h$ is
$C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$ then $K$ is
$C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$.

######  Lemma 5.

Let $\omega$ be a bounded domain of ${\mathbb{R}}^{2}$. For all $k\geq 1$, if
$h\in C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$ then $(x,y,t)\mapsto
K(x,y,t)=\int_{0}^{t}h(x,y,s)\,ds$ is also in
$C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$ and

$\displaystyle\frac{\partial^{k}K}{\partial x^{k}}(x,y,t)$
$\displaystyle=\int_{0}^{t}\frac{\partial^{k}h}{\partial x^{k}}(x,y,s)\,ds,$
$\displaystyle\frac{\partial^{k}K}{\partial y^{k}}(x,y,t)$
$\displaystyle=\int_{0}^{t}\frac{\partial^{k}h}{\partial y^{k}}(x,y,s)\,ds,$
$\displaystyle\frac{\partial^{k}K}{\partial t^{k}}(x,y,t)$
$\displaystyle=\frac{\partial^{k-1}h}{\partial t^{k-1}}(x,y,t).$

###### Proof.

We use a proof by induction on $k\geq 1$.

  * •

We first prove the result for $k=1$: Since $h\in
C^{1}(\overline{\omega}\times{\mathbb{R}}^{+})\subset
C^{0}(\overline{\omega}\times{\mathbb{R}}^{+})$, we obtain

$\frac{\partial K}{\partial t}(x,y,t)=h(x,y,t),$

and $\frac{\partial K}{\partial t}\in
C(\overline{\omega}\times{\mathbb{R}}^{+})$. Now we fix $t\in{\mathbb{R}}^{+}$
and consider the function $\widetilde{K}_{t}$ defined on $\overline{\omega}$
by

$\widetilde{K}_{t}(x,y)=K(x,y,t).$

We use the theorem of differentiability of a function defined by in integral
to show that $\widetilde{K}_{t}\in C^{1}(\overline{\omega})$. Indeed,

    * –

For all $s\in{\mathbb{R}}^{+}$, $h$ is of class $C^{1}$ in the variables
$(x,y)\in\overline{\omega}$.

    * –

Since $\frac{\partial h}{\partial x}$ and $\frac{\partial h}{\partial y}$ are
continuous on $\overline{\omega}\times{\mathbb{R}}^{+}$, there exists a
constant $M$ such that

$\displaystyle\forall(x,y)\in\overline{\omega},\ \forall s\in[0,t],\quad$
$\displaystyle\biggl{|}\frac{\partial h}{\partial x}(x,y,s)\biggr{|}\leq M,$
$\displaystyle\biggl{|}\frac{\partial h}{\partial y}(x,y,s)\biggr{|}\leq M,$

and the function $s\mapsto M$ is integrable on $[0,t]$.

In conclusion, for all $t\in{\mathbb{R}}^{+}$, $\widetilde{K}_{t}\in
C^{1}(\overline{\omega})$ and

$\displaystyle\frac{\partial K}{\partial
x}(x,y,t)=\frac{\partial\widetilde{K}_{t}}{\partial
x}(x,y)=\int_{0}^{t}\frac{\partial h}{\partial x}(x,y,s)\,ds,$
$\displaystyle\frac{\partial K}{\partial
y}(x,y,t)=\frac{\partial\widetilde{K}_{t}}{\partial
y}(x,y)=\int_{0}^{t}\frac{\partial h}{\partial y}(x,y,s)\,ds.$

We now use the theorem of continuity of a function defined by an integral to
show that $\frac{\partial K}{\partial x}$ and $\frac{\partial K}{\partial y}$
are in $C(\overline{\omega}\times{\mathbb{R}}^{+})$. Indeed, let’s detail the
proof for $\frac{\partial K}{\partial x}$. We first write

$\frac{\partial K}{\partial x}(x,y,t)=\int_{0}^{T_{f}}\frac{\partial
h}{\partial
x}(x,y,s)\mathds{1}_{[0,t]}(s)\,ds=\int_{0}^{T_{f}}i(x,y,t,s)\,ds.$ (142)

Since $\frac{\partial h}{\partial x}\in
C(\overline{\omega}\times{\mathbb{R}}^{+})$, we obtain

    * –

$i$ is continuous in the variables $(x,y,t)$ for almost every $s\in[0,T_{f}]$,
i.e.

$i(x,y,t,s)\underset{(x,y,t)\to(x_{0},y_{0},t_{0})}{\longrightarrow}i(x_{0},y_{0},t_{0},s)\quad\text{for
almost every $s\in[0,T_{f}]$}.$

    * –

There exists a constant $M$ such that

$\displaystyle\forall(x,y,t)\in\overline{\omega},\quad$
$\displaystyle|i(x,y,t,s)|\leq M,\quad\text{for almost every
$s\in[0,T_{f}]$,}$

and the function $s\mapsto M$ is integrable on $[0,T_{f}]$.

In conclusion, $K\in C^{1}(\overline{\omega}\times{\mathbb{R}}^{+})$ since the
partial derivatives exist and are continuous on
$\overline{\omega}\times{\mathbb{R}}^{+}$.

  * •

We assume the proposition true for a given $k\geq 1$ and prove it holds for
$k+1$: If $h\in C^{k+1}(\overline{\omega}\times{\mathbb{R}}^{+})$ then $h\in
C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$ and using the induction
hypothesis, it follows $K\in C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$
and

$\displaystyle\frac{\partial K}{\partial x}(x,y,t)$
$\displaystyle=\int_{0}^{t}\frac{\partial h}{\partial x}(x,y,s)\,ds,$
$\displaystyle\frac{\partial K}{\partial y}(x,y,t)$
$\displaystyle=\int_{0}^{t}\frac{\partial h}{\partial y}(x,y,s)\,ds,$
$\displaystyle\frac{\partial K}{\partial t}(x,y,t)$ $\displaystyle=h(x,y,t).$

Since $\frac{\partial h}{\partial x}\in
C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$ and $\frac{\partial
h}{\partial y}\in C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$, then using
again the induction hypothesis, we deduce $\frac{\partial K}{\partial x}\in
C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$ and $\frac{\partial
K}{\partial y}\in C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$. In addition
$\frac{\partial K}{\partial t}=h\in
C^{k}(\overline{\omega}\times{\mathbb{R}}^{+})$. In conclusion, $K\in
C^{k+1}(\overline{\omega}\times{\mathbb{R}}^{+})$.

  * •

The proposition thus holds for all $k\geq 1$, by induction.

∎

By composition, we deduce from its expression (141) that $L\in
C^{k}(\overline{\omega})$. In conclusion $W\in C^{k}(\overline{\omega})$.

##### Step 5. Genral case ($\beta\in C^{k}(\overline{\omega})$)

Using the same approach and notations as in the proof of the constant case
$\beta(x)=\alpha$ (see Appendix A), we find by the method of characteristics
that the solution of (13) is given along the characteristics by

$W(X(\xi,t))=\Bigl{(}V(\gamma(\xi))+\int_{0}^{t}g(X(\xi,s))e^{\int_{0}^{s}\beta(X(\xi,u))\,du}\,ds\Bigr{)}e^{-\int_{0}^{t}\beta(X(\xi,s))\,ds},$
(143)

and by inversion

$W(x,y)=\Bigl{(}V\bigl{(}\gamma(H_{2}(x,y))\bigr{)}+\int_{0}^{H_{2}(x,y)}g(G(H_{1}(x,y),s))e^{\int_{0}^{s}\beta(G(H_{1}(x,y),u))\,du}\,ds\Bigr{)}e^{-\int_{0}^{H_{2}(x,y)}\beta(G(H_{1}(x,y),s))\,ds}.$
(144)

The proof of the regularity follows the same ideas as in the previous case
($\beta(x)=\alpha$).

∎

##  Appendix B Hypotheses checking in Droniou’s article [14]

###### Proof.

We first recall the notation and hypotheses used in [14]. Note that these
notations are local to this proof.

##### Notation and hypotheses in [14]

Let $\Omega$ be a bounded domain in $\mathbb{R}^{N}(N\geqslant 2)$ with a
Lipschitz continuous boundary. We denote by $\mathbf{n}$ the unit normal to
$\partial\Omega$ outward to $\Omega$ and by $\sigma$ the measure on
$\partial\Omega$.

If $\Gamma$ is a measurable subset of
$\partial\Omega,W_{\Gamma}^{1,q}(\Omega)$ is the space of all functions in
$W^{1,q}(\Omega)$ (the usual Sobolev space) the trace of which is null on
$\Gamma;$ it is endowed with the same norm as $W^{1,q}(\Omega)$, that is to
say $\|v\|_{W^{1,q}(\Omega)}=\|v\|_{L^{q}(\Omega)}+\||\nabla
v|\|_{L^{q}(\Omega)}$. When $q=2$, we denote as usual $W^{1,2}=H^{1}$.

We take, $N_{*}=N$ when $N\geqslant 3$, and, $N_{*}\in{]}2,\infty[$ when
$N=2$.

We make the following hypotheses on the data.

  * •

$\Gamma_{d}$ and $\Gamma_{f}$ are measurable subsets of $\partial\Omega$ such
that $\sigma\left(\Gamma_{d}\cap\Gamma_{f}\right)=0$ and
$\partial\Omega=\Gamma_{d}\cup\Gamma_{f}$, [Hyp (8)]

  * •

$A\colon\Omega\rightarrow M_{N}(\mathbb{R})$ is a measurable matrix-valued
function which satisfies:
$\exists\alpha_{A}>0$ s.t. $A(x)\xi\cdot\xi\geqslant\alpha_{A}|\xi|^{2}\quad$
for a.e. $x\in\Omega$, for all $\xi\in\mathbb{R}^{N}$,
$\exists\Lambda_{A}>0$ s.t. $\|A(x)\|\leqslant\Lambda_{A}\quad$ for a.e.
$x\in\Omega$
(where, for $\left.M\in
M_{N}(\mathbb{R}),\|M\|:=\sup\left\\{|M\xi|,\xi\in\mathbb{R}^{N},|\xi|=1\right\\}\right)$,
[Hyp (9)]

  * •

$b\in L^{N_{*}/2}(\Omega),\quad b\geqslant 0\text{ a.e. on }\Omega$, [Hyp
(10)]

  * •

$\lambda\in L^{N_{*}-1}(\partial\Omega),\quad\lambda\geqslant 0\
\sigma\text{-a.e. on }\partial\Omega$, [Hyp (11)]

  * •

$\mathbf{v}\in\left(L^{N_{*}}(\Omega)\right)^{N}$, [Hyp (12)]

  * •

$L\in\left(H_{\Gamma_{d}}^{1}(\Omega)\right)^{\prime}$ [Hyp (13)]
(recall that $N_{*}=N$ when $N\geqslant 3$ and that $N_{*}\in{]}2,\infty[$
when $N=2$ ).

  * •

The non-convection parts of equation (26) (Eq. (4) in Droniou’s article) are
supposed to be coercive, that is to say:
$\exists b_{0}>0,\exists E\subset\Omega$ such that $b\geqslant b_{0}$ on $E$,
$\exists\lambda_{0}>0,\exists S\subset\Gamma_{f}$ such that
$\lambda\geqslant\lambda_{0}$ on $S$ and either
$\sigma\left(\Gamma_{d}\right)>0$ or $|E|>0$ or $\sigma(S)>0$. [Hyp (14)]

##### Hypotheses checking

The penalized problem (14) is of the form (26) and the dual problem (25) of
the form (27) by taking $N=d$ and

  * •

$\Gamma_{d}=\partial\Omega$, $\Gamma_{f}=\emptyset$. Thus
$\Gamma_{d}\cap\Gamma_{f}=\emptyset$ and
$\Gamma_{d}\cup\Gamma_{f}=\partial\Omega$. [Hyp (8) OK]

  * •

$A=Id$. [Hyp (9) OK]

  * •

$b=1+\frac{\chi}{\varepsilon}\alpha$. Thus $b\in L^{\infty}(\Omega)\subset
L^{N_{*}/2}(\Omega)$, and $b\geq 0$. [Hyp (10) OK]

  * •

$\lambda=0$. Thus $\lambda\in L^{N_{*}-1}(\partial\Omega)$, and $\lambda\geq
0$ on $\partial\Omega$. [Hyp (11) OK]

  * •

For the penalized problem, $L=(1-\chi)f+\frac{\chi}{\varepsilon}g$. Thus $L\in
L^{2}(\Omega)\hookrightarrow H^{-1}(\Omega)=(H^{1}_{0}(\Omega))^{\prime}$.
For the dual problem, $L$ is assumed in $H^{-1}(\Omega)$. [Hyp (13) OK]

  * •

${\bf v}=\frac{\chi}{\varepsilon}n$

where $N_{*}=N$ (dimension) when $N\geq 3$ and $N_{*}\in{]}2,\infty[$ when
$N=2$.

Also, the non-convection parts of the problem are coercive ($E=\Omega$,
$b_{0}=1$ and $S=\emptyset$). [Hyp (14) OK]

Finally, the remaining hypothesis to be checked is ${\bf
v}\in(L^{N_{*}}(\Omega))^{N}$. We have ${\bf v}=0$ in $\mathcal{U}$ and ${\bf
v}=\frac{n}{\varepsilon}$ in $\omega$, where $n$ is an extension of the normal
unit vector $\tilde{n}$ defined on $\partial\mathcal{U}$. $n$ is assumed to be
in $H^{1}(\omega)$.

In practice, we are interested in dimension $N=2$ or $3$. Because $\omega$ is
a bounded Lipschiz domain, Sobolev embeddings (see Subsection 2.2.2) give the
following:

  * •

If $N=2$, then $n\in H^{1}(\omega)\hookrightarrow L^{3}(\omega)$ [Hyp (12) OK]

  * •

If $N=3$, then $n\in H^{1}(\omega)\hookrightarrow L^{p*}(\omega)$ where
$\frac{1}{p^{*}}=\frac{1}{2}-\frac{1}{3}=\frac{1}{6}$. Thus $n\in
L^{6}(\omega)\subset L^{3}(\omega)$ [Hyp (12) OK]

In conclusion, there exists a unique solution $u_{\varepsilon}$ to the
penalized problem (14) and a unique solution $w_{\varepsilon}$ to the dual
problem (25) in dimension $2$ or $3$ if $n\in H^{1}(\omega)$. If we can choose
the lifting $n\in L^{d}(\omega)$, then the hypothesis (12) holds also in any
dimension $N=d\geq 4$.

The proof in [14] starts by proving an existence result for (27), using an
approximated solution and fixed point arguments (Leray-Schauder topological
degree [13, 22]). This existence result is then used to prove an a priori
estimate on the solution of (26) that leads to an existence result for (26).
Finally, the uniqueness results for both problems are obtained using the
linearity of these equations and a duality argument. A simplified proof of the
existence of the solution of the dual problem when $n\in L^{\infty}(\omega)$
can be found in [15].

∎

##  Appendix C Explicit construction of supersolutions in special cases

###  C.1 One-dimensional case

The one-dimensional case was studied in [7]. Here, we present an alternative
proof for the convergence of the penalization method based on the boundary
layer approach presented in this paper for the general case. The advantage of
the boundary layer approach is that it is generalizable in higher dimension,
unlike the approach used in [7] which was based on the explicit computation of
the solution.

We thus consider the following one-dimensional problem with Neumann or Robin
boundary conditions at $x=1$

$\begin{dcases}-u^{\prime\prime}+u=f\quad\text{in $\mathcal{U}={]}1,2[$}\\\
-u^{\prime}(1)+\alpha u(1)=g(1),\quad u(2)=0,\end{dcases}$ (145)

where $f\in L^{2}(]1,2[)$, $g(1)\in{\mathbb{R}}$ and $\alpha\geq 0$ are given.
The corresponding penalized problem reads

$\begin{dcases}-u_{\varepsilon}^{\prime\prime}+u_{\varepsilon}+\frac{\chi}{\varepsilon}(-u_{\varepsilon}^{\prime}+\alpha
u_{\varepsilon}-g(1))=(1-\chi)f\quad\text{in $\Omega={]}0,2[$}\\\
u_{\varepsilon}(0)=0,\quad u_{\varepsilon}(2)=0,\end{dcases}$ (146)

where $\varepsilon>0$ is a small parameter (the penalization parameter), and
$\chi$ is the characteristic function of the obstacle $\omega={]}0,1[$. Using
the notations of this paper, here $\varphi(x)=x$ and $n(x)=-1=\nabla\psi(x)$
in $\omega$ where $\psi(x)=1-x$. Following the approach used in Section 4, the
question of convergence of the asymptotic expansion reduces to finding
suitable weight functions $p_{\varepsilon}$ and $q_{\varepsilon}$ satisfying
(70) which reduces in the considered one-dimensional case to

$\begin{dcases}-q_{\varepsilon}^{\prime\prime}+\frac{q_{\varepsilon}^{\prime}}{\varepsilon}+q_{\varepsilon}=b_{\varepsilon}\geq
0&\text{in $\omega=(0,1)$}\\\
-p_{\varepsilon}^{\prime\prime}+p_{\varepsilon}=a_{\varepsilon}\geq 0&\text{in
$\mathcal{U}=(1,2)$}\\\ p_{\varepsilon}(1)=q_{\varepsilon}(1)\\\
p_{\varepsilon}^{\prime}(1)=q_{\varepsilon}^{\prime}(1)-\frac{q_{\varepsilon}(1)}{\varepsilon}\\\
p_{\varepsilon}\geq\beta>0&\text{in $\mathcal{U}=(1,2)$}\\\
q_{\varepsilon}\geq\beta>0&\text{in $\omega=(0,1)$}\\\
\|q_{\varepsilon}\|_{L^{\infty}(0,\delta)}\leq C\text{ and
}\|q_{\varepsilon}\|_{L^{\infty}(\delta,1)}\leq\frac{C}{\varepsilon},\end{dcases}$
(147)

where $0<\delta<1$.

In this one-dimensional case, we can exhibit suitable supersolutions
satisfying (147): we take

$q_{\varepsilon}(x)=1+\frac{1}{\varepsilon}e^{\frac{x-1}{\varepsilon}}-\frac{1}{\varepsilon}e^{-\frac{1}{\varepsilon}}\quad\text{in
$\omega=(0,1)$},$ (148)

so that

$q_{\varepsilon}(x)\geq 1-\frac{1}{\varepsilon}e^{-\frac{1}{\varepsilon}}\to
1\ \text{when $\varepsilon\to 0$}$ (149)

and

$\displaystyle-
q_{\varepsilon}^{\prime\prime}+\frac{q_{\varepsilon}^{\prime}}{\varepsilon}+q_{\varepsilon}=q_{\varepsilon}\geq
0.$ (150)

To satisfy the transmission conditions, we choose

$p_{\varepsilon}(x)=q_{\varepsilon}(1)+\Bigl{(}q_{\varepsilon}^{\prime}(1)-\frac{q_{\varepsilon}(1)}{\varepsilon}\Bigr{)}(x-1)\quad\text{in
$\mathcal{U}=(1,2)$},$ (151)

so that

$p_{\varepsilon}(2)=q_{\varepsilon}(1)+q_{\varepsilon}^{\prime}(1)-\frac{q_{\varepsilon}(1)}{\varepsilon}=\Bigl{(}1-\frac{1}{\varepsilon}\Bigr{)}\Bigl{(}1+\frac{1}{\varepsilon}-\frac{1}{\varepsilon}e^{-\frac{1}{\varepsilon}}\Bigr{)}+\frac{1}{\varepsilon^{2}}=1+\frac{1}{\varepsilon^{2}}e^{-\frac{1}{\varepsilon}}\to
1\ \text{when $\varepsilon\to 0$}$ (152)

and

$-p_{\varepsilon}^{\prime\prime}+p_{\varepsilon}=p_{\varepsilon}\geq\min(p_{\varepsilon}(1),p_{\varepsilon}(2))=\min(q_{\varepsilon}(1),p_{\varepsilon}(2))\geq
0\quad\text{for $\varepsilon$ small enough}.$ (153)

Thus the strict positivity of $p_{\varepsilon}$ and $q_{\varepsilon}$ for
sufficiently small $\varepsilon>0$. Finally, for all $x\in(0,\delta)$ with
$\delta<1$,

$\displaystyle|q_{\varepsilon}(x)|\leq
1+\frac{1}{\varepsilon}e^{\frac{\delta-1}{\varepsilon}}-\frac{1}{\varepsilon}e^{-\frac{1}{\varepsilon}}\leq
C$

for sufficiently small $\varepsilon>0$, and, for $x\in(\delta,1)$,

$\displaystyle|q_{\varepsilon}(x)|\leq
1+\frac{1}{\varepsilon}-\frac{1}{\varepsilon}e^{-\frac{1}{\varepsilon}}\leq\frac{C}{\varepsilon}$

for sufficiently small $\varepsilon$, then the desired estimates on
$\|q_{\varepsilon}\|_{L^{\infty}}$. This completes the existence of suitable
weight functions $p_{\varepsilon}$, $q_{\varepsilon}$ in the one-dimensional
case.

######  Remark 4.

We thus recover the results obtained in [7] using a boundary layer approach.
In fact, we obtain a finer description of what happens in the obstacle domain,
which was not carried out in [7].

Figure 11: Plot of suitable supersolutions $q_{\varepsilon}$ (149) in $(0,1)$
and $p_{\varepsilon}$ (151) in $(1,2)$ satisfying (147) in the one-dimensional
case, for different values of $\varepsilon$.

###  C.2 Spherical symmetric case (in 2D)

We consider here $\mathcal{U}=B(0,1)$ and $\Omega=B(0,2)$ in 2D. In this case,
we can choose $n=\vec{e}_{r}=\nabla\psi$ where
$\psi(x,y)=\sqrt{x^{2}+y^{2}}-1=r-1$ and look for weight functions
$p_{\varepsilon}$, $q_{\varepsilon}$ satisfying (70) that are radial

$\displaystyle p_{\varepsilon}(x,y)=p_{\varepsilon}(r)\quad\text{and}\quad
q_{\varepsilon}(x,y)=q_{\varepsilon}(r).$

Using polar coordinates,

$\begin{dcases}\Delta f=\frac{1}{r}\frac{\partial}{\partial
r}\Bigl{(}r\frac{\partial f}{\partial
r}\Bigr{)}+\frac{1}{r^{2}}\frac{\partial^{2}f}{\partial\theta^{2}}=\frac{\partial^{2}f}{\partial
r^{2}}+\frac{1}{r}\frac{\partial f}{\partial
r}+\frac{1}{r^{2}}\frac{\partial^{2}f}{\partial\theta^{2}}\\\
\mathrm{div}(A)=\frac{1}{r}\frac{\partial}{\partial
r}\Bigl{(}rA_{r}\Bigr{)}+\frac{1}{r}\frac{\partial
A_{\theta}}{\partial\theta}=\frac{\partial A_{r}}{\partial
r}+\frac{A_{r}}{r}+\frac{1}{r}\frac{\partial
A_{\theta}}{\partial\theta}.\end{dcases}$ (154)

Thus, (70) rewrites in the radial symmetric case (see for instance [31])

$\begin{dcases}-p_{\varepsilon}^{\prime\prime}-\frac{1}{r}p_{\varepsilon}^{\prime}+p_{\varepsilon}=a_{\varepsilon}(r)\geq
0&\text{in $]0,1[$}\\\
-q_{\varepsilon}^{\prime\prime}-\Bigl{(}\frac{1}{r}+\frac{1}{\varepsilon}\Bigr{)}q_{\varepsilon}^{\prime}+\Bigl{(}1-\frac{1}{\varepsilon
r}\Bigr{)}q_{\varepsilon}=b_{\varepsilon}(r)\geq 0&\text{in $]1,2[$}\\\
p_{\varepsilon}^{\prime}(0)=0&\\\ p_{\varepsilon}(1)=q_{\varepsilon}(1)&\\\
p_{\varepsilon}^{\prime}(1)=q_{\varepsilon}^{\prime}(1)+\frac{q_{\varepsilon}(1)}{\varepsilon}&\\\
p_{\varepsilon}\geq\beta>0&\text{in $]0,1[$}\\\
q_{\varepsilon}\geq\beta>0&\text{in $]0,2[$}\\\
\|q_{\varepsilon}\|_{L^{\infty}(1+\delta,2)}\leq C\text{ and
}\|q_{\varepsilon}\|_{L^{\infty}(1,1+\delta)}\leq\frac{C}{\varepsilon}.\end{dcases}$
(155)

We construct by hand suitable weight functions satisfying the previous
conditions (155): we take

$p_{\varepsilon}(r)=\frac{1}{\varepsilon}J_{0}(ir)$ (156)

where $J_{0}$ is the Bessel function of the first kind, then $p_{\varepsilon}$
satisfies

$\begin{dcases}-p_{\varepsilon}^{\prime\prime}-\frac{1}{r}p_{\varepsilon}^{\prime}+p_{\varepsilon}=0&\text{in
$]0,1[$}\\\ p_{\varepsilon}^{\prime}(0)=0\\\
p_{\varepsilon}(r)\geq\frac{1}{\varepsilon}J_{0}(0)=\frac{1}{\varepsilon}>0&\text{in
$]0,1[$}\end{dcases}$ (157)

using the properties of the Bessel functions [1]. And we take
$q_{\varepsilon}$ of the form

$q_{\varepsilon}(r)=\frac{d_{\varepsilon}}{r}+\frac{e_{\varepsilon}}{\varepsilon}e^{-\frac{r-1}{\varepsilon}}.$
(158)

The transmission conditions at $r=1$ give, since
$\frac{d}{dr}\bigl{(}J_{0}(ir)\bigr{)}=-iJ_{1}(ir)$,

$\begin{dcases}d_{\varepsilon}=\frac{iJ_{1}(i)}{\varepsilon-1}\to-
iJ_{1}(i)>0\\\ e_{\varepsilon}=J_{0}(i)+\frac{\varepsilon
iJ_{1}(i)}{\varepsilon-1}\to J_{0}(i)>0,\end{dcases}$ (159)

when $\varepsilon\to 0^{+}$. Therefore, for sufficiently small $\varepsilon$,

$q_{\varepsilon}(r)\geq\frac{d_{\varepsilon}}{r}\geq\frac{C}{2}>0\qquad\text{in
$]1,2[$}.$ (160)

One then has

$\displaystyle-
q_{\varepsilon}^{\prime\prime}-\Bigl{(}\frac{1}{r}+\frac{1}{\varepsilon}\Bigr{)}q_{\varepsilon}^{\prime}+\Bigl{(}1-\frac{1}{\varepsilon
r}\Bigr{)}q_{\varepsilon}$
$\displaystyle=q_{\varepsilon}+d_{\varepsilon}\Bigl{(}-\Bigl{(}\frac{1}{r}\Bigr{)}^{\prime\prime}-\frac{1}{r}\Bigl{(}\frac{1}{r}\Bigr{)}^{\prime}\Bigr{)}$
(161)
$\displaystyle=q_{\varepsilon}+d_{\varepsilon}\Bigl{(}\frac{-1}{r^{3}}\Bigr{)}$
(162)
$\displaystyle=\frac{e_{\varepsilon}}{\varepsilon}e^{-\frac{r-1}{\varepsilon}}+d_{\varepsilon}\frac{r^{2}-1}{r^{3}}\geq
0$ in $]1,2[$ (163)

for sufficiently small $\varepsilon$. The previous calculations follow from
the fact that $r\mapsto\frac{1}{r}$ solves

$-u^{\prime}-\frac{1}{r}u=0$

and $r\mapsto e^{-\frac{r-1}{\varepsilon}}$ solves

$-u^{\prime\prime}-\Bigl{(}\frac{1}{r}+\frac{1}{\varepsilon}\Bigr{)}u^{\prime}-\frac{1}{\varepsilon
r}u=0.$

Finally, the estimates on $q_{\varepsilon}$ in (155) are also satisfied since,
for sufficiently small $\varepsilon$,

$|q_{\varepsilon}(r)|\leq
C+C^{\prime}\frac{1}{\varepsilon}e^{-\frac{\delta}{\varepsilon}}\to
0,\quad\text{for $r\in{]1+\delta},2[$}$ (164)

and

$|q_{\varepsilon}(r)|\leq
C+C^{\prime}\frac{1}{\varepsilon}\leq\frac{C}{\varepsilon},\quad\text{for
$r\in{]1},1+\delta[$}.$ (165)

######  Remark 5.

The approach used here could also be used in the one-dimensional case. The
approach could also be adapted to the spherical symmetric case in $n$
dimensions.

######  Remark 6.

The choices made for $q_{\varepsilon}$ in the special cases (1D or spherical
symmetric case) will be understood in the following appendix where we use a
formal boundary layer approach for the dual problem.

(a) Plots of $p_{\varepsilon}(r)$ and $q_{\varepsilon}(r)$

(b) Plots of $a_{\varepsilon}(r)$ and $b_{\varepsilon}(r)$

Figure 12: Plot of suitable supersolutions $p_{\varepsilon}(r)$ (156) in
$(0,1)$ and $q_{\varepsilon}(r)$ (158) in $(1,2)$ satisfying (155) in the
spherical symmetric case (2D), for different values of $\varepsilon$.

##  Appendix D Formal boundary layer approach for the dual problem

We are going to use again a boundary layer approach for the dual problem

$\displaystyle-\Delta
p_{\varepsilon}+p_{\varepsilon}=a_{\varepsilon}=\frac{a^{-1}(x)}{\varepsilon}+a^{0}(x)+\ldots\geq
0$ in $\mathcal{U}$ (166a) $\displaystyle-\Delta
q_{\varepsilon}-\frac{1}{\varepsilon}\mathrm{div}(q_{\varepsilon}n)+q_{\varepsilon}=b_{\varepsilon}=\frac{b^{-1}(x)}{\varepsilon}+b^{0}(x)+\ldots\geq
0$ in $\omega$ (166b) $\displaystyle p_{\varepsilon}=q_{\varepsilon}$ on
$\partial\mathcal{U}$ (166c) $\displaystyle\frac{\partial
p_{\varepsilon}}{\partial\nu}=\frac{\partial
q_{\varepsilon}}{\partial\nu}+\frac{q_{\varepsilon}}{\varepsilon}$ on
$\partial\mathcal{U}$ (166d) $\displaystyle
q_{\varepsilon}=q_{0}(x)+\varepsilon q_{1}(x)+\ldots$ on $\partial\Omega$
(166e)

The reader might think that we are going around in circles, but the goal here
is not to estabish a convergence for the dual problem, but simply to use the
boundary layer approach in order to construct an approximate solution which we
wish to be a supersolution and satisfy the expected properties in (70).

Based on numerical experiments, we consider the following ansatz for the
solution (boundary layer near $\partial\mathcal{U}$):

$\begin{dcases}p_{\varepsilon}(x)=\frac{1}{\varepsilon}P^{-1}(x)+P^{0}(x)+\varepsilon
P^{1}(x)+\ldots,\\\
q_{\varepsilon}(x)=\frac{1}{\varepsilon}Q^{-1}\Bigl{(}x,\frac{\psi(x)}{\varepsilon}\Bigr{)}+Q^{0}\Bigl{(}x,\frac{\psi(x)}{\varepsilon}\Bigr{)}+\varepsilon
Q^{1}\Bigl{(}x,\frac{\psi(x)}{\varepsilon}\Bigr{)}+\ldots\end{dcases}$ (167)

where $\psi$ is somehow the distance to $\partial\mathcal{U}$, and the profile
terms inside $\omega$ have the form

$\begin{dcases}Q^{i}(x,z)=\overline{Q}^{i}(x)+\widetilde{Q}^{i}(x,z)\\\
\text{where}\ \forall i,k\geq 0,\
\partial_{z}^{k}\widetilde{Q}^{i}\xrightarrow[z\to+\infty]{}0,\quad\forall
x\in\omega,\end{dcases}$ (168)

Using the same notations as before, for a function $(x,z)\mapsto Q(x,z)$, the
function $x\mapsto q(x)=Q\bigl{(}x,\frac{\psi(x)}{\varepsilon}\bigr{)}$
satisfies

$\begin{dcases}\nabla q=\nabla Q+\frac{1}{\varepsilon}Q_{z}\nabla\psi\\\
\Delta q=\Delta Q+\frac{2}{\varepsilon}\nabla
Q_{z}\cdot\nabla\psi+\frac{1}{\varepsilon}\Delta\psi
Q_{z}+\frac{1}{\varepsilon^{2}}|\nabla\psi|^{2}Q_{zz},\end{dcases}$ (169)

where the derivatives of $q$ are evaluated at $x$ while the derivatives of $Q$
are evaluated at $\bigl{(}x,\frac{\psi(x)}{\varepsilon}\bigr{)}$.

For simplicity, we assume $\psi$ to be sufficiently regular and $n=\nabla\psi$
in $\omega$, we also assume $\nabla\psi\neq 0$ in $\omega$. In this cas, using

$\displaystyle\mathrm{div}(q_{\varepsilon}n)$
$\displaystyle=q_{\varepsilon}\mathrm{div}(n)+n\cdot\nabla q_{\varepsilon}$
$\displaystyle=q_{\varepsilon}\Delta\psi+\nabla
q_{\varepsilon}\cdot\nabla\psi,$

Eq. (166b) reduces to

$-\Delta q_{\varepsilon}-\frac{1}{\varepsilon}\nabla
q_{\varepsilon}\cdot\nabla\psi-\frac{1}{\varepsilon}\Delta\psi
q_{\varepsilon}+q_{\varepsilon}=b_{\varepsilon}=\frac{b^{-1}(x)}{\varepsilon}+b^{0}(x)+\ldots\geq
0\quad\text{in $\omega$}.$ (170)

###  D.1 Determination of profiles

####  D.1.1 Asymptotic expansion of Eq. (170) inside the obstacle $\omega$

  * •

Order $\varepsilon^{-3}$: identifying the terms corresponding to the power
$-3$ of $\varepsilon$ leads to

$-|\nabla\psi|^{2}Q^{-1}_{zz}-|\nabla\psi|^{2}Q^{-1}_{z}=0,$

i.e.

$-Q^{-1}_{zz}-Q^{-1}_{z}=0.$

Using the decomposition (168) one obtains

$-\widetilde{Q}_{zz}^{-1}-\widetilde{Q}_{z}^{-1}={0}\quad\text{in
$\omega\times{\mathbb{R}}^{+}$}.$ (171)

  * •

Order $\varepsilon^{-2}$:

$-|\nabla\psi|^{2}Q^{0}_{zz}-\Delta\psi Q^{-1}_{z}-2\nabla
Q^{-1}_{z}\cdot\nabla\psi-Q^{0}_{z}|\nabla\psi|^{2}-\nabla
Q^{-1}\cdot\nabla\psi-\Delta\psi Q^{-1}=0,$

which, by using hypothesis (168) and taking the limit when $z\to+\infty$,
leads to

$-\nabla\overline{Q}^{-1}\cdot\nabla\psi-\Delta\psi\overline{Q}^{-1}=0\quad\text{in
$\omega$},$

i.e.

$\mathrm{div}(\overline{Q}^{-1}n)=0\quad\text{in $\omega$},$ (172)

and by difference, we obtain

$-|\nabla\psi|^{2}\widetilde{Q}_{zz}^{0}-\Delta\psi\widetilde{Q}_{z}^{-1}-2\nabla\widetilde{Q}^{-1}_{z}\cdot\nabla\psi-\widetilde{Q}^{0}_{z}|\nabla\psi|^{2}-\nabla\widetilde{Q}_{z}^{-1}\cdot\nabla\psi-\Delta\psi\widetilde{Q}^{-1}=0\quad\text{in
$\omega\times{\mathbb{R}}^{+}$}.$ (173)

  * •

Order $\varepsilon^{-1}$:

$-|\nabla\psi|^{2}Q^{1}_{zz}-\Delta\psi Q^{0}_{z}-2\nabla
Q^{0}_{z}\cdot\nabla\psi-\Delta Q^{-1}-Q^{1}_{z}|\nabla\psi|^{2}-\nabla
Q^{0}\cdot\nabla\psi-\Delta\psi Q^{0}+Q^{-1}=b^{-1},$

and again using the same reasoning as above, by using hypothesis (168) and
$z\to+\infty$, we obtain

$-\Delta\overline{Q}^{-1}-\nabla\overline{Q}^{0}\cdot\nabla\psi-\Delta\psi\overline{Q}^{0}+\overline{Q}^{-1}=b^{-1}\quad\text{in
$\omega$},$

i.e.

$\nabla\overline{Q}^{0}\cdot\nabla\psi+\Delta\psi\overline{Q}^{0}=\mathrm{div}(\overline{Q}^{0}n)=-\Delta\overline{Q}^{-1}+\overline{Q}^{-1}-b^{-1}\quad\text{in
$\omega$},$ (174)

and by difference,

$-|\nabla\psi|^{2}\widetilde{Q}_{zz}^{1}-\Delta\psi\widetilde{Q}_{z}^{0}-2\nabla\widetilde{Q}^{0}_{z}\cdot\nabla\psi-\Delta\widetilde{Q}^{-1}-\widetilde{Q}^{1}_{z}|\nabla\psi|^{2}-\nabla\widetilde{Q}^{0}\cdot\nabla\psi-\Delta\psi\widetilde{Q}^{0}+\widetilde{Q}^{-1}=0\quad\text{in
$\omega\times{\mathbb{R}}^{+}$}.$ (175)

  * •

Order $\varepsilon^{0}$: from

$-|\nabla\psi|^{2}Q^{2}_{zz}-\Delta\psi Q^{1}_{z}-2\nabla
Q^{1}_{z}\cdot\nabla\psi-\Delta Q^{0}-Q^{2}_{z}|\nabla\psi|^{2}-\nabla
Q^{1}\cdot\nabla\psi-\Delta\psi Q^{1}+Q^{0}=b^{0},$

we deduce once again as before,

$-\Delta\overline{Q}^{0}-\nabla\overline{Q}^{1}\cdot\nabla\psi-\Delta\psi\overline{Q}^{1}+\overline{Q}^{0}=b^{0}\quad\text{in
$\omega$},$

i.e.

$\nabla\overline{Q}^{1}\cdot\nabla\psi+\Delta\psi\overline{Q}^{1}=\mathrm{div}(\overline{Q}^{1}n)=-\Delta\overline{Q}^{0}+\overline{Q}^{0}-b^{0}\quad\text{in
$\omega$},$ (176)

and by difference,

$-|\nabla\psi|^{2}\widetilde{Q}_{zz}^{2}-\Delta\psi\widetilde{Q}_{z}^{1}-2\nabla\widetilde{Q}^{1}_{z}\cdot\nabla\psi-\Delta\widetilde{Q}^{0}-\widetilde{Q}^{2}_{z}|\nabla\psi|^{2}-\nabla\widetilde{Q}^{1}\cdot\nabla\psi-\Delta\psi\widetilde{Q}^{1}+\widetilde{Q}^{0}=0\quad\text{in
$\omega\times{\mathbb{R}}^{+}$}.$ (177)

  * •

Order $\varepsilon^{1}$: from

$-|\nabla\psi|^{2}Q^{3}_{zz}-\Delta\psi Q^{2}_{z}-2\nabla
Q^{2}_{z}\cdot\nabla\psi-\Delta Q^{1}-Q^{3}_{z}|\nabla\psi|^{2}-\nabla
Q^{2}\cdot\nabla\psi-\Delta\psi Q^{2}+Q^{1}=b^{1},$

we deduce once again as before,

$-\Delta\overline{Q}^{1}-\nabla\overline{Q}^{2}\cdot\nabla\psi-\Delta\psi\overline{Q}^{2}+\overline{Q}^{2}=b^{1}\quad\text{in
$\omega$},$

i.e.

$\nabla\overline{Q}^{2}\cdot\nabla\psi+\Delta\psi\overline{Q}^{2}=\mathrm{div}(\overline{Q}^{2}n)=-\Delta\overline{Q}^{1}+\overline{Q}^{1}-b^{1}\quad\text{in
$\omega$},$ (178)

and by difference,

$-|\nabla\psi|^{2}\widetilde{Q}_{zz}^{3}-\Delta\psi\widetilde{Q}_{z}^{2}-2\nabla\widetilde{Q}^{2}_{z}\cdot\nabla\psi-\Delta\widetilde{Q}^{1}-\widetilde{Q}^{3}_{z}|\nabla\psi|^{2}-\nabla\widetilde{Q}^{2}\cdot\nabla\psi-\Delta\psi\widetilde{Q}^{2}+\widetilde{Q}^{1}=0\quad\text{in
$\omega\times{\mathbb{R}}^{+}$}.$ (179)

####  D.1.2 Asymptotic expansion of Eq. (166b) inside the fluid domain
$\mathcal{U}$

  * •

Order $\varepsilon^{-1}$:

$-\Delta P^{-1}+P^{-1}=a^{-1}(x)\quad\text{in $\mathcal{U}$}.$ (180)

  * •

Order $\varepsilon^{0}$:

$-\Delta P^{0}+P^{0}=a^{0}(x)\quad\text{in $\mathcal{U}$}.$ (181)

  * •

Order $\varepsilon^{1}$:

$-\Delta P^{1}+P^{1}=a^{1}(x)\quad\text{in $\mathcal{U}$}.$ (182)

####  D.1.3 Asymptotic expansion of Eq. (166c) (transmission condition on
$\partial\mathcal{U}$)

We obtain simply

$\displaystyle Q^{-1}(x,0)$ $\displaystyle=P^{-1}(x)$ $\displaystyle
Q^{0}(x,0)$ $\displaystyle=P^{0}(x)\qquad\text{on $\partial\mathcal{U}$}$
$\displaystyle Q^{1}(x,0)$ $\displaystyle=P^{1}(x)$

i.e.

$\displaystyle\overline{Q}^{-1}(x)+\widetilde{Q}^{-1}(x,0)$
$\displaystyle=P^{-1}(x)$ (183)
$\displaystyle\overline{Q}^{0}(x,0)+\widetilde{Q}^{0}(x,0)$
$\displaystyle=P^{0}(x)\qquad\text{on $\partial\mathcal{U}$}.$ (184)
$\displaystyle\overline{Q}^{1}(x,0)+\widetilde{Q}^{1}(x,0)$
$\displaystyle=P^{1}(x)$ (185)

####  D.1.4 Asymptotic expansion of Eq. (166d) (transmission condition on
$\partial\mathcal{U}$)

We have

$\nabla p_{\varepsilon}\cdot\nu=\nabla
q_{\varepsilon}\cdot\nu+\frac{q_{\varepsilon}}{\varepsilon}=\nabla
q_{\varepsilon}\cdot\nu+\frac{p_{\varepsilon}}{\varepsilon}\quad\text{on
$\partial\mathcal{U}$}.$

  * •

Order $\varepsilon^{-2}$

$0=Q^{-1}_{z}\nabla\psi\cdot\nu+P^{-1}=Q^{-1}_{z}+P^{-1}$

  * •

Order $\varepsilon^{-1}$

$\nabla P^{-1}\cdot\nu=Q^{0}_{z}\nabla\psi\cdot\nu+\nabla
Q^{-1}\cdot\nu+P^{0}=Q^{0}_{z}+\nabla Q^{-1}\cdot\nu+P^{0}$

  * •

Order $\varepsilon^{0}$

$\nabla P^{0}\cdot\nu=Q^{1}_{z}\nabla\psi\cdot\nu+\nabla
Q^{0}\cdot\nu+P^{1}=Q^{1}_{z}+\nabla Q^{0}\cdot\nu+P^{1}.$

Thus

$\displaystyle\widetilde{Q}^{-1}_{z}+\overline{Q}^{-1}+\widetilde{Q}^{-1}=0$
(186) $\displaystyle\nabla
P^{-1}\cdot\nu=\widetilde{Q}^{0}_{z}+\nabla\overline{Q}^{-1}\cdot\nu+\nabla\widetilde{Q}^{-1}\cdot\nu+\overline{Q}^{0}+\widetilde{Q}^{0}$
on $\partial\mathcal{U}\times\\{z=0\\}$. (187) $\displaystyle\nabla
P^{0}\cdot\nu=\widetilde{Q}^{1}_{z}+\nabla\overline{Q}^{0}\cdot\nu+\nabla\widetilde{Q}^{0}\cdot\nu+\overline{Q}^{1}+\widetilde{Q}^{1}$
(188)

####  D.1.5 Asymptotic expansion of Eq. (166e) (Dirichlet boundary conditions
on $\partial\Omega$)

By neglecting
$\widetilde{Q}^{i}\bigl{(}x,\frac{\psi(x)}{\varepsilon}\bigr{)}\to 0$ when
$\varepsilon\to 0$ (we will find later an exponential decrease), we obtain

$\begin{dcases}\overline{Q}^{-1}(x)=0&\\\ \overline{Q}^{0}(x)=q_{0}&\text{on
$\partial\Omega$}.\\\ \overline{Q}^{1}(x)=q_{1}&\end{dcases}$ (189)

###  D.2 Formal resolution of the profile equations

####  D.2.1 Determination of $Q^{-1}$ and $P^{-1}$

From (172) and (189), we obtain

$\overline{Q}^{-1}=0\quad\text{in $\omega$}.$ (190)

From (171) and the fact that $\widetilde{Q}^{-1}$ tends to $0$ as
$z\to+\infty$ , we obtain

$\widetilde{Q}^{-1}(x,z)=A(x)e^{-z}\quad\text{in
$\omega\times{\mathbb{R}}^{+}$}.$

The condition (186) rewrites

$\widetilde{Q}^{-1}_{z}(x,0)+\widetilde{Q}^{-1}(x,0)=0\quad\text{on
$\partial\mathcal{U}$}$

and is satisfied since $\widetilde{Q}^{-1}$ satisfy everywhere

$\widetilde{Q}^{-1}_{z}(x,z)=-\widetilde{Q}^{-1}(x,z)\quad\text{in
$\omega\times{\mathbb{R}}^{+}$}.$ (191)

From (180) and (183), we obtain that $P^{-1}$ is solution to

$\begin{dcases}-\Delta P^{-1}+P^{-1}=a^{-1}&\text{in $\mathcal{U}$}\\\
P^{-1}(x)=\widetilde{Q}^{-1}(x,0)=A(x)&\text{on
$\partial\mathcal{U}$},\end{dcases}$ (192)

and a possible choice for $\widetilde{Q}^{-1}$ is given by

$\widetilde{Q}^{-1}(x,z)=\underline{P}^{-1}(x)e^{-z}\quad\text{in
$\omega\times{\mathbb{R}}^{+}$},$ (193)

where $\underline{P}^{-1}$ is a simultaneous lifting in $\omega$ of the trace
of $P^{-1}$ and of that of $\frac{\partial P^{-1}}{\partial\nu}$ on
$\partial\mathcal{U}$ and having the same regularity as $P^{-1}$.

Notice that $A(x)$ on $\partial\mathcal{U}$ is not yet determined at this
stage. It will be using the equations on $P^{0},Q^{0}$.

####  D.2.2 Determination of $Q^{0}$ and $P^{0}$

Similarly, from (174) and (189), we obtain that $\overline{Q}^{0}$ is solution
to

$\begin{dcases}\mathrm{div}(\overline{Q}^{0}n)=-\Delta\overline{Q}^{-1}+\overline{Q}^{-1}-b^{-1}=-b^{-1}&\text{in
$\omega$}\\\ \overline{Q}^{0}(x)=q_{0}(x)&\text{on
$\partial\Omega$}.\end{dcases}$ (194)

From (173), we have

$\displaystyle-|\nabla\psi|^{2}\widetilde{Q}_{zz}^{0}-|\nabla\psi|^{2}\widetilde{Q}^{0}_{z}$
$\displaystyle=\Delta\psi\widetilde{Q}_{z}^{-1}+2\nabla\widetilde{Q}^{-1}_{z}\cdot\nabla\psi+\nabla\widetilde{Q}_{z}^{-1}\cdot\nabla\psi+\Delta\psi\widetilde{Q}^{-1}$
$\displaystyle=-\nabla\widetilde{Q}_{z}^{-1}\cdot\nabla\psi$
$\displaystyle=-\nabla\underline{P}^{-1}\cdot\nabla\psi\,e^{-z}$
$\displaystyle-\widetilde{Q}_{zz}^{0}-\widetilde{Q}^{0}_{z}$
$\displaystyle=-\frac{\nabla\underline{P}^{-1}\cdot\nabla\psi}{|\nabla\psi|^{2}}\,e^{-z}=F(x)\,e^{-z},$

using (191) and the obtained expression of $\widetilde{Q}^{-1}$ (193). The
solution which tends to $0$ as $z\to+\infty$ is given by

$\widetilde{Q}^{0}(x,z)=B(x)e^{-z}+F(x)ze^{-z}\quad\text{in
$\omega\times{\mathbb{R}}^{+}$},$ (195)

where $B(x)$ is to be determined.

Using (186), we have

$\displaystyle\nabla
P^{-1}\cdot\nu=\widetilde{Q}^{0}_{z}(x,0)+\nabla\overline{Q}^{-1}\cdot\nu+\nabla\widetilde{Q}^{-1}\cdot\nu+\overline{Q}^{0}(x)+\widetilde{Q}^{0}(x,0)$
on $\partial\mathcal{U}$.

i.e.

$\displaystyle\widetilde{Q}^{0}_{z}(x,0)+\widetilde{Q}^{0}(x,0)$
$\displaystyle=\nabla
P^{-1}\cdot\nu-\nabla\overline{Q}^{-1}\cdot\nu-\overline{Q}^{0}(x)$ on
$\partial\mathcal{U}$ $\displaystyle=(\nabla
P^{-1}\cdot\nu)_{|\mathcal{U}}-(\nabla\underline{P}^{-1}\cdot\nu)_{|\omega}-\overline{Q}^{0}(x)$
on $\partial\mathcal{U}$ $\displaystyle=-\overline{Q}^{0}(x)$
$\displaystyle\text{on $\partial\mathcal{U}$}.$

Using the expression found for $\widetilde{Q}^{0}$, we have

$\displaystyle\widetilde{Q}^{0}_{z}(x,z)+\widetilde{Q}^{0}(x,z)$
$\displaystyle=-B(x)e^{-z}+F(x)e^{-z}-F(x)ze^{-z}+B(x)e^{-z}+F(x)ze^{-z}$
$\displaystyle=F(x)\,e^{-z}.$

Thus,

$\displaystyle F(x)$ $\displaystyle=-\overline{Q}^{0}(x)$ on
$\partial\mathcal{U}$
$\displaystyle-\frac{\nabla\underline{P}^{-1}\cdot\nabla\psi}{|\nabla\psi|^{2}}$
$\displaystyle=-\overline{Q}^{0}(x)$ on $\partial\mathcal{U}$
$\displaystyle\frac{\partial P^{-1}}{\partial\nu}$
$\displaystyle=\overline{Q}^{0}$ $\displaystyle\text{on
$\partial\mathcal{U}$}.$

since $\nabla\psi=\nu$ on $\partial\mathcal{U}$. This fixes $A(x)$ on
$\partial\mathcal{U}$, since we obtain that $P^{-1}$ is solution to

$\begin{dcases}-\Delta P^{-1}+P^{-1}=a^{-1}&\text{in $\mathcal{U}$}\\\
\frac{\partial P^{-1}}{\partial\nu}=\overline{Q}^{0}&\text{on
$\partial\mathcal{U}$},\end{dcases}$ (196)

and $A(x)$ is then the trace of $P^{-1}$ on $\partial\mathcal{U}$.

From (181) and (183), we obtain that $P^{0}$ is solution to

$\begin{dcases}-\Delta P^{0}+P^{0}=a^{0}&\text{in $\mathcal{U}$}\\\
P^{0}(x)=\widetilde{Q}^{0}(x,0)+\overline{Q}^{0}(x)=B(x)+\overline{Q}^{0}(x)&\text{on
$\partial\mathcal{U}$}.\end{dcases}$ (197)

Thus we can take $B(x)=\underline{P}^{0}-\overline{Q}^{0}$ in $\omega$, where
$\underline{P}^{0}$ is a simultaneous lifting in $\omega$ of the trace of
$P^{0}$ and of that of $\frac{\partial P^{0}}{\partial\nu}$ on
$\partial\mathcal{U}$ and having the same regularity as $P^{0}$. The equation
satisfied by $P^{0}$ will be determined using the equations on $P^{1},Q^{1}$.

####  D.2.3 Determination of $Q^{1}$ and $P^{1}$

From (176) and (189), we obtain that $\overline{Q}^{1}$ is solution to

$\begin{dcases}\mathrm{div}(\overline{Q}^{1}n)=-\Delta\overline{Q}^{0}+\overline{Q}^{0}-b^{0}&\text{in
$\omega$}\\\ \overline{Q}^{1}(x)=q_{1}(x)&\text{on
$\partial\Omega$}.\end{dcases}$ (198)

From (175) and the obtained expression of $\widetilde{Q}^{0}$ (195), we have

$\displaystyle-|\nabla\psi|^{2}\widetilde{Q}_{zz}^{1}-|\nabla\psi|^{2}\widetilde{Q}^{1}_{z}$
$\displaystyle=\Delta\psi\widetilde{Q}_{z}^{0}+2\nabla\widetilde{Q}^{0}_{z}\cdot\nabla\psi+\Delta\widetilde{Q}^{-1}+\nabla\widetilde{Q}^{0}\cdot\nabla\psi+\Delta\psi\widetilde{Q}^{0}-\widetilde{Q}^{-1}$
$\displaystyle=\Delta\psi
F(x)\,e^{-z}+(\Delta\underline{P}^{-1}-\underline{P}^{-1})e^{-z}+2\nabla(F(x)e^{-z}-\widetilde{Q}^{0})\cdot\nabla\psi+\nabla\widetilde{Q}^{0}\cdot\nabla\psi$
$\displaystyle=\Delta\psi
F(x)\,e^{-z}+(\Delta\underline{P}^{-1}-\underline{P}^{-1})e^{-z}-\nabla\widetilde{Q}^{0}\cdot\nabla\psi+2\nabla
F(x)\cdot\nabla\psi e^{-z}$ $\displaystyle=\Delta\psi
F\,e^{-z}+(\Delta\underline{P}^{-1}-\underline{P}^{-1})e^{-z}-\nabla
B\cdot\nabla\psi\,e^{-z}-\nabla F\cdot\nabla\psi\,ze^{-z}+2\nabla
F\cdot\nabla\psi\,e^{-z}$
$\displaystyle-\widetilde{Q}_{zz}^{1}-\widetilde{Q}^{1}_{z}$
$\displaystyle=\frac{\Delta\psi
F+\Delta\underline{P}^{-1}-\underline{P}^{-1}-\nabla B\cdot\nabla\psi+2\nabla
F\cdot\nabla\psi}{|\nabla\psi|^{2}}\,e^{-z}-\frac{\nabla
F\cdot\nabla\psi}{|\nabla\psi|^{2}}\,ze^{-z}$
$\displaystyle=G(x)\,e^{-z}+H(x)\,ze^{-z}.$

The solution which tends to $0$ as $z\to+\infty$ is given by

$\widetilde{Q}^{1}(x,z)=C(x)e^{-z}+(G(x)+H(x))ze^{-z}+\frac{H(x)}{2}z^{2}e^{-z}\quad\text{in
$\omega\times{\mathbb{R}}^{+}$},$

where $C(x)$ is to be determined.

Using (186), we have

$\displaystyle\nabla
P^{0}\cdot\nu=\widetilde{Q}^{1}_{z}+\nabla\overline{Q}^{0}\cdot\nu+\nabla\widetilde{Q}^{0}\cdot\nu+\overline{Q}^{1}+\widetilde{Q}^{1}$
on $\partial\mathcal{U}$.

i.e.

$\displaystyle\widetilde{Q}^{1}_{z}(x,0)+\widetilde{Q}^{1}(x,0)$
$\displaystyle=\nabla
P^{0}\cdot\nu-\nabla\overline{Q}^{0}\cdot\nu-\nabla\widetilde{Q}^{0}(x,0)\cdot\nu-\overline{Q}^{1}$
on $\partial\mathcal{U}$ $\displaystyle=\nabla
P^{0}\cdot\nu-\nabla\overline{Q}^{0}\cdot\nu-\nabla
B\cdot\nu-\overline{Q}^{1}$ $\displaystyle\text{on $\partial\mathcal{U}$}.$

Using the expression found for $\widetilde{Q}^{1}$, we have

$\displaystyle[\widetilde{Q}^{1}_{z}(x,z)+\widetilde{Q}^{1}(x,z)]_{z=0}$
$\displaystyle=-C(x)+(G(x)+H(x))+C(x)$ $\displaystyle=G(x)+H(x).$

Thus

$\displaystyle G(x)+H(x)$ $\displaystyle=\frac{\partial
P^{0}}{\partial\nu}-\frac{\partial\overline{Q}^{0}}{\partial\nu}-\nabla
B\cdot\nu-\overline{Q}^{1}$ on $\partial\mathcal{U}$ $\displaystyle\Delta\psi
F+\Delta\underline{P}^{-1}-\underline{P}^{-1}-\nabla B\cdot\nu+\nabla
F\cdot\nu$ $\displaystyle=\frac{\partial
P^{0}}{\partial\nu}-\frac{\partial\overline{Q}^{0}}{\partial\nu}-\nabla
B\cdot\nu-\overline{Q}^{1}$ on $\partial\mathcal{U}$

and it follows

$\displaystyle\frac{\partial P^{0}}{\partial\nu}$
$\displaystyle=\Delta\psi(-\overline{Q}^{0})+\Delta\underline{P}^{-1}-\underline{P}^{-1}+\nabla
F\cdot\nu+\frac{\partial\overline{Q}^{0}}{\partial\nu}+\overline{Q}^{1}$ on
$\partial\mathcal{U}$
$\displaystyle=2\frac{\partial\overline{Q}^{0}}{\partial\nu}+b^{-1}+\Delta\underline{P}^{-1}-\underline{P}^{-1}+\nabla
F\cdot\nu+\overline{Q}^{1}$ on $\partial\mathcal{U}$,

since $\overline{Q}^{0}$ satisfies
$\overline{Q}^{0}\Delta\psi+\nabla\overline{Q}^{0}\cdot\nabla\psi=-b^{-1}$,
where we recall that

$F(x)=-\frac{\nabla\underline{P}^{-1}\cdot\nabla\psi}{|\nabla\psi|^{2}}\quad\text{in
$\omega$}.$

Therefore, we obtain that $P^{0}$ is solution to

$\begin{dcases}-\Delta P^{0}+P^{0}=a^{0}&\text{in $\mathcal{U}$}\\\
\frac{\partial
P^{0}}{\partial\nu}=2\frac{\partial\overline{Q}^{0}}{\partial\nu}+b^{-1}+\Delta\underline{P}^{-1}-\underline{P}^{-1}+\nabla
F\cdot\nu+\overline{Q}^{1}&\text{on $\partial\mathcal{U}$},\end{dcases}$ (199)

and $B(x)$ is then well defined by $\underline{P}^{0}-\overline{Q}^{0}$.