# On the diffusion approximation of the stationary radiative transfer equation
with absorption and emission

Elena Demattè, Juan J.L. Velázquez Institute for Applied Mathematics,
University of Bonn, 53115 Bonn, Germany.E-mail<EMAIL_ADDRESS>bonn.deInstitute for Applied Mathematics, University of Bonn, 53115 Bonn,
Germany.E-mail<EMAIL_ADDRESS>
###### Abstract

We study the situation in which the distribution of temperature a body is due
to its interaction with radiation. We consider the boundary value problem for
the stationary radiative transfer equation under the assumption of the local
thermodynamic equilibrium. We study the diffusion equilibrium approximation in
the absence of scattering. We consider absorption coefficient independent of
the frequency $\nu$ (the so-called Grey approximation) and the limit when the
photons’ mean free path tends to zero, i.e. the absorption coefficient tends
to infinity. We show that the densitive of radiative energy $u$, which is
proportional to the fourth power of the temperature due to the Stefan-
Boltzmann law, solves in the limit an elliptic equation where the boundary
value can be determined uniquely in terms of the original boundary condition.
We derive formally with the method of matched asymptotic expansions the
boundary condition for the limit problem and we prove rigorously the
convergence to the solution of the limit problem with a careful analysis of
some non-local integral operators. The method developed here allows to prove
all the results using only maximum principle arguments.

Acknowledgments: The authors gratefully acknowledge the financial support of
the collaborative research centre The mathematics of emerging effects (CRC
1060, Project-ID 211504053) and Bonn International Graduate School of
Mathematics (BIGS) at the Hausdorff Center for Mathematics founded through the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC-2047/1 – 390685813.

Keywords: Radiative transfer equation, diffusion approximation, stationary
solution, maximum principle, boundary layers.

Statements and Declarations: The authors have no relevant financial or non-
financial interests to disclose.

Data availability: Data sharing not applicable to this article as no datasets
were generated or analysed during the current study.

###### Contents

  1. 1 Introduction
    1. 1.1 Motivation and previous results
    2. 1.2 Structure of the paper and notation
  2. 2 Derivation of the limit problem
    1. 2.1 Formal derivation of the limit problem in the diffusive equilibrium approximation
    2. 2.2 Formal derivation of boundary condition for the limit problem in the diffusive equilibrium approximation
    3. 2.3 Some properties of the kernel
  3. 3 The boundary condition for the limit problem
    1. 3.1 The homogeneous equation
    2. 3.2 Well-posedness theory for the inhomogeneous equation
    3. 3.3 Asymptotic behavior of the bounded solution of the inhomogeneous equation
  4. 4 Rigorous proof of the diffusion equilibrium approximation for constant absorption coefficient
    1. 4.1 Derivation of the equation for $u^{\varepsilon}$
    2. 4.2 Uniform boundedness of $u^{\varepsilon}$
    3. 4.3 Estimates of $u^{\varepsilon}-\overline{u}$ near the boundary $\partial\Omega$
    4. 4.4 Convergence of $u^{\varepsilon}$ to the solution of the new boundary value problem
  5. 5 Diffusion approximation for space dependent absorption coefficient
    1. 5.1 The limit problem and the boundary layer equation
    2. 5.2 Rigorous proof of the convergence: equation for $u^{\varepsilon}$ and properties of the kernel
    3. 5.3 Rigorous proof of the convergence: uniform boundedness of $u^{\varepsilon}$
    4. 5.4 Rigorous proof of the convergence: estimates of $u^{\varepsilon}-\overline{u}$ near the boundary $\partial\Omega$
    5. 5.5 Rigorous proof of the convergence of $u^{\varepsilon}$ to the solution of the new boundary value problem

##  1 Introduction

The radiative transfer equation is the kinetic equation which describes the
distribution of energy and direction of motions of a set of photons, which can
be absorbed and scattered by a medium. This equation can be used to describe
the transfer of heat in a material due to radiative processes. The radiative
transfer equation can be written in its more general form as

$\frac{1}{c}\partial_{t}I_{\nu}(x,n,t)+n\cdot\nabla_{x}I_{\nu}(x,n,t)=\alpha_{\nu}^{e}-\alpha_{\nu}^{a}I_{\nu}(x,n,t)-\alpha_{\nu}^{s}I_{\nu}(x,n,t)+\alpha_{\nu}^{s}\int_{\mathbb{S}^{2}}K(n,n^{\prime})I_{\nu}(x,n^{\prime},t)dn^{\prime}.$
(1.1)

We denote by $I_{\nu}(x,n,t)$ the intensity of radiation (i.e. radiating
energy) of frequency $\nu$ at position $x\in\Omega$ and in direction
$n\in\mathbb{S}^{2}$ and at time $t\geq 0$. The coefficients
$\alpha_{\nu}^{a}$, $\alpha_{\nu}^{e}$ and $\alpha_{\nu}^{s}$ are respectively
the absorption, the emission and the scattering coefficient. In the scattering
term the kernel is normalized such that
$\int_{\mathbb{S}^{2}}K(n,n^{\prime})dn^{\prime}=1$. The speed of light is
indicated by $c$.
In this paper we focus on the stationary problem and on processes, where the
scattering is negligible. Therefore the equation we will study reduces to

$n\cdot\nabla_{x}I_{\nu}\left(x,n\right)=\alpha_{\nu}^{e}-\alpha^{a}_{\nu}I_{\nu}\left(x,n\right).$
(1.2)

In this article we consider the situation of local thermal equilibrium (LTE),
which means that at every point $x\in\Omega$ there is a well-defined
temperature $T(x)\geq 0$. This yields according to the Kirchhoff’s law (cf.
[37]) the following relation for the absorption and emission coefficient

$\alpha_{\nu}^{e}(x)=\alpha_{\nu}^{a}(x)B_{\nu}(T(x)),$

where $B_{\nu}(T(x))=\frac{2h\nu^{3}}{c^{2}}\frac{1}{e^{\frac{h\nu}{kT}}-1}$
is the Plank emission of a black body at temperature $T(x)$ and $k$ the
Boltzmann constant. Moreover, it is well-known that

$\int_{0}^{\infty}B_{\nu}\left(T(x)\right)\;d\nu=\sigma T^{4}(x),$ (1.3)

where $\sigma=\frac{2\pi^{4}k^{4}}{15h^{3}c^{2}}$ is the Stefan-Boltzmann
constant. We will denote for simplicity from now on as the absorption
coefficient $\alpha_{\nu}^{a}$ as $\alpha_{\nu}$.

The solution $I_{\nu}(x,n)$ of (1.2) can be used to compute the flux of energy
at each point $x\in\Omega$ of the material, which is given by

$\mathcal{F}(x):=\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;n\>I_{\nu}\left(x,n\right).$

In the stationary case, if the temperature is independent of time at every
point, the total incoming and outgoing flux of energy should balance. In
mathematical terms this can be formulated by the condition for the flux of
energy to be divergence free, i.e.

$\nabla_{x}\cdot\mathcal{F}(x)=0.$

This situation is denoted in the physical literature by pointwise radiative
equilibrium.

We study the situation when the radiation is coming from a very far source of
infinite distance. This can be formalized in mathematical terms by means of
the boundary condition

$I_{\nu}\left(x,n\right)=g_{\nu}\left(n\right)\geq 0$ (1.4)

if $x\in\partial\Omega$ and $n\cdot N_{x}<0$ for $N_{x}$ the outer normal
vector of the boundary at point $x\in\partial\Omega$. Throughout this paper we
will consider $\Omega\subset\mathbb{R}^{3}$ to be a bounded convex domain with
$C^{3}$-boundary.

We are concerned in this paper with the study of the diffusion approximation
that arises in optically thick media. This means, that we consider the case
when the optical depth is very large compared to the characteristic length of
the system. Hence, we rescale and for $\varepsilon\ll 1$ we consider the
following boundary value problem

$\begin{cases}n\cdot\nabla_{x}I_{\nu}\left(x,n\right)=\frac{\alpha_{\nu}(x)}{\varepsilon}\left(B_{\nu}\left(T\left(x\right)\right)-I_{\nu}\left(x,n\right)\right)&x\in\Omega,\\\
\nabla_{x}\cdot\mathcal{F}=0&x\in\Omega,\\\
I_{\nu}\left(x,n\right)=g_{\nu}\left(n\right)&x\in\partial\Omega\text{ and
}n\cdot N_{x}<0.\end{cases}$ (1.5)

For the solution to this equation we will prove that the intensity of
radiation $I_{\nu}(x,n)$ is approximately the Plank distribution
$B_{\nu}(T(x))$ with the local temperature at each point $x\in\Omega$, i.e. we
will show

$I_{\nu}^{\varepsilon}(x,n)\to B_{\nu}(T(x))\;\;\;\;\;\;\;\;\text{ as
}\varepsilon\to 0.$ (1.6)

Notice however, that this approximation cannot be expected for points $x$ that
are close to the boundary $\partial\Omega$. The situation in which (1.6) holds
is denoted in the physical literature as the diffusion equilibrium
approximation (see e.g. [24] and [37]). More precisely, we will consider the
limit problem when $\varepsilon\to 0$ and we will rigorously prove that it is
given by a Dirichlet problem for the heat equation of the temperature with
boundary value uniquely determined by the incoming source $g_{\nu}(n)$ and the
outer normal $N_{x}$ for $x\in\partial\Omega$. The main result we will prove
in this paper is for the so called Grey approximation, i.e. the case when the
absorption coefficient is independent of the frequency $\nu$. The main reason
for that is that some of the estimates are already in this case very
technical. Hopefully, the type of methods we are developing in this paper can
be extended to the non-Grey case.

######  Theorem 1.1.

Let $\alpha_{\nu}(x)=\alpha(x)$ independent of $\nu$, $\alpha\in
C^{3}\left(\Omega\right)$, $g_{\nu}\geq 0$ with
$\int_{0}^{\infty}g_{\nu}(n)\;d\nu\in L^{\infty}\left(\mathbb{S}\right)$ in
(1.5), $\Omega$ bounded convex with $C^{3}$-boundary and strictly positive
curvature. Let $T_{\varepsilon}$ be the temperature associated to the
intensity $I_{\nu}$ solution to the initial value problem (1.5). Then there
exists a functional
$T_{\Omega}:L^{\infty}\left(\mathbb{S}^{2},L^{1}\left(\mathbb{R}_{+}\right)\right)\to
C\left(\partial\Omega\right)$ which maps $g_{\nu}$ to a continuous function
$T_{\Omega}[g_{\nu}](p)$ on the boundary $p\in\partial\Omega$ such that

$T_{\varepsilon}(x)\to T(x)$

uniformly in every compact subset of $\Omega$, where $T$ is the solution to
the Dirichlet problem

$\begin{cases}-\operatorname*{div}\left(\frac{\sigma 4T^{3}}{\alpha}\nabla
T\right)=0&x\in\Omega,\\\
T(p)=T_{\Omega}[g_{\nu}](p)&p\in\partial\Omega.\end{cases}$

###  1.1 Motivation and previous results

The computation of the distribution of temperature of matter interacting with
radiation is an important issue in many physical application and in addition
it rises interesting mathematical questions. The kinetic equation describing
the interaction of matter with radiation is the radiative transfer equation. A
detailed explanation of its derivation and its main properties can be found in
[5, 24, 26, 31, 37]. In particular, in [24, 37] there is an extensive
discussion about the diffusion equilibrium approximation and the situations
where this can be expected or not.

Since the earlier result by Compton [6] in 1922 the interaction of a gas with
radiation has been extensively studied. Milne for example studied a simplified
model, where the radiation is monochromatic and the gas density depends only
on one space variable (cf. [25]).

A question which has been much studied in mathematical literature is the
situation in which $\alpha_{\nu}^{e}=\alpha_{\nu}^{a}=0$ in (1.1), i.e. the
interaction between matter and radiation is due to scattering only. In this
case the problem reduces to

$\frac{1}{c}\partial_{t}I_{\nu}(x,n,t)+n\cdot\nabla_{x}I_{\nu}(x,n,t)=-\alpha_{\nu}^{s}(x)I_{\nu}(x,n,t)+\alpha_{\nu}^{s}(x)\int_{\mathbb{S}^{2}}K(n,n^{\prime})I_{\nu}(x,n^{\prime},t)dn^{\prime}.$
(1.7)

The same equation arises also in the study of neutron transport, a problem
which has been extensively studied in mathematics.

It turns out that in the Grey approximation the problem (1.5) can be reduced
exactly to the study of a particular neutron transport equation, namely the
case when the kernel $K$ is constant $1$. Indeed, denoting by
$J(x,n)=\int_{0}^{\infty}I_{\nu}(x,n)\;d\nu$ and combining the first two
equations of (1.5) we obtain
$\int_{0}^{\infty}B_{\nu}(T(x))\;d\nu=\fint_{\mathbb{S}^{2}}J(x,n)\;dn=\frac{1}{4\pi}\int_{\mathbb{S}^{2}}J(x,n)\;dn$.
Hence, equation (1.5) is equivalent to the study of

$\begin{cases}n\cdot\nabla_{x}J(x,n)=\frac{\alpha(x)}{\varepsilon}\left(\fint_{\mathbb{S}^{2}}J(x,n)\;dn-J(x,n)\right)&\text{
if }x\in\Omega,\\\ J(x,n)=\int_{0}^{\infty}g_{\nu}(n)\;d\nu&\text{ if
}x\in\partial\Omega,\;n\cdot N_{x}<0.\end{cases}$ (1.8)

However, the equivalence between (1.5) and (1.8) does not hold in the non-Grey
case. The properties of equation (1.8) as well as the diffusion approximation
limit have been studied for a long time, starting with the seminal paper [4]
of 1979, where the stationary version of (1.7) was studied. In that work the
authors proved the diffusion approximation for the neutron transport equation
using a stochastic method. The result they obtained for $J$ would imply in
particular our main Theorem 1.1.

More recently, in a series of papers [16, 33, 34, 35, 36] Yan Guo and Lei Wu
have studied the diffusion approximation of both the stationary and the time
dependent neutron transport equation (1.7) when $K\equiv 1$ and
$\alpha^{s}_{\nu}(x)\equiv 1$, independent of $x$, for different classes of
boundary conditions in $2$ and $3$ dimensions, in bounded convex domains or
annuli (in 2D). In particular the result in paper [34] imply again the main
Theorem 1.1 when $\alpha\equiv 1$. Their proof relies on PDE methods and not
on a stochastic approach. Moreover, they also computed the geometric
approximation in the structure of the boundary layer.

The main goal of this paper is to develop a method which allows to obtain
diffusive limit approximations like the one in Theorem 1.1 for the radiative
transfer equation (1.1) using PDE methods that rely only in maximum principle
tools. This tools are different from those used by Guo and Wu. Specifically,
the method in [16, 33, 34, 35, 36] relies on the $L^{2}$-$L^{p}$-$L^{\infty}$
estimates that were introduced for the analysis of kinetic equations by Yan
Guo in [15]. In particular, the method is based on the estimates of the
velocity distribution $J$. Our approach is based on the direct derivation of
estimates for the temperature $T(x)$ associated to a given distribution of
radiation $I_{\nu}(x,n)$. More precisely, equation (1.5) can be reformulated
as a non-local integral equation for the temperature (c.f. [21]). In the case
of the Grey approximation we have the following equation for $u(x)=4\pi\sigma
T^{4}(x)$

$u(x)-\int_{\Omega}K_{\varepsilon}(x,\eta)u(\eta)\;d\eta=S(x),$ (1.9)

where the precise form of the kernel $K_{\varepsilon}$ and of the source
$S(x)$ are discussed in Sections 4 and 5.

Equation (1.9) can be thought as a non-local elliptic equation which in
particular satisfies good properties, such as the maximum principle.
Specifically, our proof relies only in finding suitable supersolutions and
applying the maximum principle. The way in which we constructed these
supersolutions is mimicking particular solutions of elliptic equations with
constant coefficients. These supersolutions give also an insight of the
behavior of the solution near the boundary $\partial\Omega$. Our hope is that
the method developed in this paper could be extended to the non-Grey case, at
least for some suitable choice of $\alpha_{\nu}(x)$. One reason why this
should be possible is that [21] shows how to solve the non-local equation
(1.9) for some class of non-Grey problems.

Another type of diffusion approximation for (1.1) is the one in [13, 14] in
which it has been considered the situation when $\alpha_{\nu}^{s}\to\infty$
while $\alpha_{\nu}^{e}$ and $\alpha_{\nu}^{a}$ remain bounded combined with
the equation for balancing the energy either in the one dimensional case or in
the whole space.

The well-posedness and the diffusion approximation of the time dependent
problem (1.7) in the frame work of $L^{1}$-functions using the theory of
$m$-accretive operators has ben studied in a series of papers [2, 3].
Seemingly, although the techniques in these papers allow to develop a theory
for the time dependent problem, they do not provide information about the
stationary solution.

Some versions of the stationary problem involving the radiative transfer
equation can be found in [22, 23, 27, 32]. The problems studied in these
papers include also heat conduction and different type of boundary condition
of our model (for a more detailed discussion see [21]).

It is important to emphasize that equation (1.5) is very different in the non-
Grey case from the scattering problem (1.7), in the sense that the system
(1.5) provides an equation for the temperature. Specifically, the equation
$\nabla_{x}\cdot\mathcal{F}=0$ is automatic satisfied in the stationary
version of (1.7). Physically, this is due to the fact that the radiation
arriving at every point is just deflected. Equation (1.5) plays the same role
as the Laplace equation in order to describe the stationary distribution of
temperature in systems where the energy is transported by means of heat
conduction. In the case of (1.5) the energy is transported by means of
radiation which results in non-locality for determining the temperature
distribution. The fact that the determination of the temperature in a body
where the energy is transported by radiation is non-local was first formulated
in [18]. Since the approximation (1.6) fails at the boundary, some boundary
layers appears for which the intensity of radiation $I_{\nu}^{\varepsilon}$
differs from the Plank distribution $B_{\nu}(T)$. Hence, a careful analysis
must be made for these boundary layers where the radiation is out of
equilibrium. This will be essential in order to determine the functional
$T_{\Omega}$ in Theorem 1.1, which defines the temperature at every point of
the boundary.

Finally, we mention that one can consider more complicated interactions
between radiation and matter. For instance when the matter that interacts with
radiation is a moving fluid. (cf. [11, 12, 24, 37]). The case when the
interacting medium is a Boltzmann gas whose molecules can be in different
energy layers has been considered in [8, 20, 26, 30].

###  1.2 Structure of the paper and notation

We aim to prove Theorem 1.1. In Section 2 we will formally derive the limit
problem using ideas from the theory of matched asymptotic expansions and from
boundary layers theory. Sections 3 and 4 deal with the case of constant
absorption coefficient $\alpha\equiv 1$, while Section 5 shows the diffusion
approximation in the case of space dependent coefficient $\alpha\in
C^{3}(\Omega)$. In Section 3 we will study the properties of the solution for
the non-local integral equation describing the distribution of energy at the
boundary layer. In particular, this will allow us to construct in the case of
$\alpha\equiv 1$ the functional $T_{\Omega}$ of Theorem 1.1 which assigns the
boundary value of the limit solution at every point of the boundary. Important
tools we will use for the well-posedness are the maximum principle and a
combination of Fourier methods with the classical tools of sub- and super-
solution which resembles the Perron method for the Laplace equation. For the
asymptotic theory we use the theory of Fourier transform for distributions.
Section 4 deals with the rigorous proof of the convergence to the diffusion
equilibrium approximation for the limit problem in the constant absorption
coefficient case. Here the main tool is the maximum principle for the non-
local integral operator we can construct for the boundary value problem (1.5).
Finally, in Section 5 we consider the more general Grey approximation in which
$\alpha\in C^{3}(\overline{\Omega})$ is not constant. We will derive formally
the limit problem and then we will extend the theory developed in Sections 3
and 4 for this case. We will hence prove again by means of the maximum
principle and suitable supersolutions the convergence to the diffusion
equilibrium approximation.

We introduce here some notation we will use throughout this paper. First of
all, $\Omega\subset\mathbb{R}^{3}$ is an open bounded convex domain with
$C^{3}$-boundary and strictly positive curvature. In order to avoid
meaningless constants we assume without loss of generality that
$0\in\overline{\Omega}$. $N_{x}$ indicates always the outer normal vector for
a point $x\in\partial\Omega$.

We assume $\Omega$ to be convex in order to simplify some geometrical
argument. First of all this assumption implies that for every point
$p\in\partial\Omega$ the tangent plane to the boundary at $p$ divided the
space $\mathbb{R}^{3}$ in two disjoint half-spaces, one of them containing the
whole domain $\Omega$. This will be used several times in the definition for
every point $p\in\partial\Omega$ of the isometric transformation mapping $p$
to $0$ and $\Omega$ in the positive half-space
$\mathbb{R}_{+}\times\mathbb{R}^{2}$. The assumption of convexity can be
relaxed and the geometrical estimates should still hold, but we would need a
more careful analysis of the geometry of the problem.
Moreover, for $g_{\nu}(n)\geq 0$ with $\int_{0}^{\infty}g_{\nu}(n)\;d\nu\in
L^{\infty}\left(\mathbb{S}^{2}\right)$ we define the norms

$\Arrowvert
g\Arrowvert_{1}:=\int_{0}^{\infty}\int_{\mathbb{S}^{2}}g_{\nu}(n)\;d\nu\;dn$
(1.10)

and

$\Arrowvert
g\Arrowvert_{\infty}:=\sup\limits_{n\in\mathbb{S}^{2}}\left(\int_{0}^{\infty}g_{\nu}(n)\;d\nu\right).$
(1.11)

######  Remark.

The reason why we are assuming the seemingly restrictive boundary condition
(1.4) is because we are supposing that the source of radiation is placed at
infinity. We can obtain analogous results to the one of the paper if we
consider the more general boundary condition $g_{\nu}(n,x)$ depending also on
$x\in\partial\Omega$. In addition to the assumption above we need to require
$g_{\nu}(n,x)$ to be a $C^{1}$-function with respect to $x\in\partial\Omega$.

Figure 1: Representation of the boundary value problem.

For any point $p\in\partial\Omega$ we choose a fixed isometry mapping $p$ to
$0$ and the vector $N_{p}$ to $-e_{1}$. We will denote this rigid motion by
$\mathcal{R}_{p}:\mathbb{R}^{3}\to\mathbb{R}^{3}$ with the following
properties

$\mathcal{R}_{p}(p)=0\;\;\;\;\;\;\;\text{ and
}\;\;\;\;\;\;\;\mathcal{R}_{p}(N_{p})=-e_{1}.$ (1.12)

Finally, we define by

$\begin{split}\pi_{\partial\Omega}:\left\\{x\in\mathbb{R}^{3}:\text{dist}(x,\partial\Omega)<\delta\right\\}&\to\partial\Omega\\\
x&\mapsto\pi_{\partial\Omega}(x)\end{split}$ (1.13)

the projection to the unique closest point in the boundary $\partial\Omega$.
This function is continuous and well-defined in small neighborhood of
$\partial\Omega$, i.e. for $\delta>0$ small enough.

##  2 Derivation of the limit problem

###  2.1 Formal derivation of the limit problem in the diffusive equilibrium
approximation

We first remind how to obtain formally the equation for the interior in the
limit problem. First of all we expand the intensity of radiation

$I_{\nu}\left(x,n\right)=f_{\nu}^{0}\left(x,n\right)+\varepsilon
f_{\nu}^{1}\left(x,n\right)+\varepsilon^{2}...$ (2.1)

Substituting it in the first equation of (1.5) and identifying the terms
containing $\varepsilon^{-1}$ and $\varepsilon^{0}$ we see

$f_{\nu}^{0}\left(x,n\right)=B_{\nu}\left(T\left(x\right)\right)$

and

$f_{\nu}^{1}\left(x,n\right)=-\frac{1}{\alpha_{\nu}(x)}n\cdot\nabla_{x}B_{\nu}\left(T\left(x\right)\right)$

Using the second equation in (1.5) and the expansion in (2.1) we deduce

$\begin{split}0&=\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;n\cdot\nabla_{x}I_{\nu}\left(x,n\right)\\\
&=\operatorname*{div}\left[\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;nB_{\nu}\left(T\left(x\right)\right)\right]-\varepsilon\operatorname*{div}\left[\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\>\left(n\otimes
n\right)\frac{1}{\alpha_{\nu}(x)}\nabla_{x}B_{\nu}\left(T\left(x\right)\right)\right]\\\
&=-\varepsilon\frac{4}{3}\pi\operatorname*{div}\left(\left(\int_{0}^{\infty}d\nu\frac{1}{\alpha_{\nu}(x)}\nabla_{x}B_{\nu}\left(T\left(x\right)\right)\right)\right),\end{split}$

where we used

$\int_{\mathbb{S}^{2}}dn(n\otimes n)=\frac{4}{3}\pi\text{Id}\;\;\;\;\;\text{
and }\;\;\;\;\;\int_{\ss^{2}}dn\;n=0.$

Therefore,

$\operatorname*{div}\left(\kappa\left(T\right)\nabla_{x}T\right)=0,$ (2.2)

where
$\kappa\left(T\right):=\int_{0}^{\infty}d\nu\frac{\partial_{T}B_{\nu}\left(T\left(x\right)\right)}{\alpha_{\nu}(x)}$.
In the particular case of the Grey approximation when $\alpha_{\nu}(x)=1$ we
have $\kappa(T)=4\sigma T^{3}(x)$. Then defining $u(x):=4\pi\sigma T^{4}(x)$
we obtain the following equation

$\Delta u=0.$

This is the limit problem we will study.

###  2.2 Formal derivation of boundary condition for the limit problem in the
diffusive equilibrium approximation

In order to obtain the intensity of radiation closed to the boundary of
$\Omega$ we derive a boundary layer equation, whose solution will be used to
determine the value of the temperature at the boundary by means of a matched
argument. Suppose that $x_{0}\in\partial\Omega$, without loss of generality we
can assume $x_{0}=0$ and $N_{x_{0}}=N=-e_{1}$ using the rigid motion
$\mathcal{R}_{x_{0}}$ defined in (1.12) and putting
$\overline{g}_{\nu}(n):=g_{\nu}\left(\mathcal{R}^{-1}_{x_{0}}(n)\right)$. We
rescale $x=\varepsilon y$, where $y\in\frac{1}{\varepsilon}\Omega$. Thus, at
the leading order as $\varepsilon\to 0$ we obtain
$\alpha_{\nu}(x)=\alpha_{\nu}(\varepsilon
y)=\alpha_{\nu}(0)+\mathcal{O}(\varepsilon)$. Taking $\varepsilon\to 0$ we
obtain that the intensity of radiation satisfies

$\begin{cases}n\cdot\nabla_{y}I_{\nu}\left(y,n\right)=\alpha_{\nu}(0)\left(B_{\nu}\left(T\left(y\right)\right)-I_{\nu}\left(y,n\right)\right)&y\in\mathbb{R}_{+}\times\mathbb{R}^{2}\\\
\nabla_{y}\cdot\mathcal{F}=0&y\in\mathbb{R}_{+}\times\mathbb{R}^{2}\\\
I_{\nu}\left(y,n\right)=\overline{g}_{\nu}\left(n\right)&y\in\\{0\\}\times\mathbb{R}^{2}\text{
and }n\cdot N<0\end{cases}$ (2.3)

The first equation can be solved for $I_{\nu}$ using the method of
characteristics. Given $y\in\mathbb{R}_{+}\times\mathbb{R}^{2}$ and
$n\in\mathbb{S}^{2}$ with $n\cdot N<0$ we call $Y(y,n)$ the unique point
belonging to
$\partial\left(\mathbb{R}_{+}\times\mathbb{R}^{2}\right)=\\{0\\}\times\mathbb{R}^{2}$
such that

$y=Y(y,n)+s(y,n)n,$

where $s(y,n)=\left|y-Y(y,n)\right|$. Notice that $s(y,n)$ is the distance to
the first intersection point of the boundary $\\{0\\}\times\mathbb{R}^{2}$
with the half line $\\{y-tn:t>0\\}$. For $n\cdot N\geq 0$ we define
$s(y,n)=\infty$. Solving the equation by characteristics we obtain

$\begin{split}I_{\nu}\left(y,n\right)=&\overline{g}_{\nu}(n)e^{-\alpha_{\nu}(0)s(y,n)}\raisebox{2.0pt}{$\chi$}_{n\cdot
N<0}+\int_{0}^{s(y,n)}\;e^{-\alpha_{\nu}(0)t}\alpha_{\nu}(0)B_{\nu}\left(T\left({y-tn}\right)\right)dt.\end{split}$

Using the second equation in the rescaled problem (2.3) we calculate

$\begin{split}0=&\operatorname*{div}\left[\int_{0}^{\infty}d\nu\int_{n\cdot
N<0}dn\;n\overline{g}_{\nu}(n)e^{-\alpha_{\nu}(0)s(y,n)}+\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\int_{0}^{s(y,n)}dt\;ne^{-\alpha_{\nu}(0)t}\alpha_{\nu}B_{\nu}\left(T\left({y-tn}\right)\right)\right]\\\
=&-\int_{0}^{\infty}d\nu\int_{n\cdot
N<0}dn\;\overline{g}_{\nu}(n)\alpha_{\nu}n\cdot\nabla_{y}s(y,n)e^{-\alpha_{\nu}(0)s(y,n)}\\\
&+\operatorname*{div}\left(\int_{0}^{\infty}d\nu\int_{\mathbb{R}_{+}\times\mathbb{R}^{2}}d\eta\;\frac{y-\eta}{\left|y-\eta\right|^{3}}e^{-\alpha_{\nu}(0)\left|y-\eta\right|}\alpha_{\nu}(0)B_{\nu}\left(T\left({\eta}\right)\right)\right)\\\
=&-\int_{0}^{\infty}d\nu\int_{n\cdot
N<0}dn\;\overline{g}_{\nu}(n)\alpha_{\nu}(0)e^{-\alpha_{\nu}(0)s(y,n)}+4\pi\int_{0}^{\infty}d\nu(0)\;\alpha_{\nu}B_{\nu}\left(T\left(y\right)\right)\\\
&-\int_{0}^{\infty}d\nu\int_{\mathbb{R}_{+}\times\mathbb{R}^{2}}d\eta\;\frac{\alpha_{\nu}^{2}(0)}{\left|y-\eta\right|^{2}}e^{-\alpha_{\nu}(0)\left|y-\eta\right|}B_{\nu}\left(T\left({\eta}\right)\right).\end{split}$
(2.4)

The second equality holds via the spherical change of variable

$\begin{split}\mathbb{S}^{2}\times\mathbb{R}_{+}&\to\mathbb{R}_{+}\times\mathbb{R}^{2}\\\
(n,t)&\mapsto\eta=y-tn\end{split}$

so that $n=\frac{y-\eta}{\left|y-\eta\right|}$. For the third inequality we
use on the one hand that
$\operatorname*{div}_{y}\left(\frac{y-\eta}{\left|y-\eta\right|^{3}}\right)=4\pi\delta(y-\eta)$
and on the other hand that $n\cdot\nabla_{y}s(y,n)=1$. The latter can be seen
by the fact that

$Y(y,n)+s(y+tn,n)n=y+tn=Y(y+tn,n)+s(y+tn,n)n.$

This implies that $Y(y+tn,n)$ is $t$-constant and therefore
$1=\partial_{t}s(y+tn,n)=\left(\nabla_{y}s(y+tn,n)\right)\cdot n$. We assume
now that the temperature depends only on the first variable. This can be
expected because we are considering limits for $\varepsilon\ll 1$ and hence
the temperature can be considered to depend only on the distance to the point
$x_{0}$, which is approximated by the first variable in this setting. After
the change of variables $\xi=(y_{2}+\eta_{2},y_{3}-\eta_{3})$ and calling
$y-\eta:=y_{1}-\eta_{1}$ the last integral in (2.4) can be written as

$\int_{0}^{\infty}d\nu\int_{\mathbb{R}_{+}}d\eta\int_{\mathbb{R}^{2}}d\xi\;\alpha_{\nu}^{2}(0)\frac{e^{-\alpha_{\nu}(0)\sqrt{(y-\eta)^{2}+|\xi|^{2}}}}{(y-\eta)^{2}+|\xi|^{2}}B_{\nu}\left(T\left(\eta\right)\right).$

Using polar coordinates we obtain

$\begin{split}\int_{\mathbb{R}^{2}}d\xi\;\frac{e^{-\alpha_{\nu}(0)\sqrt{(y-\eta)^{2}+|\xi|^{2}}}}{(y-\eta)^{2}+|\xi|^{2}}=&\pi\int_{|y-\eta|^{2}}^{\infty}dx\;\frac{e^{-\alpha_{\nu}(0)\sqrt{x}}}{x}=2\pi\int_{\alpha_{\nu}(0)|y-\eta|}^{\infty}dt\;\frac{e^{-t}}{t}=4\pi
K(\alpha_{\nu}(0)|y-\eta|),\end{split}$ (2.5)

where we will denote $K(x)=\frac{1}{2}\int_{|x|}^{\infty}dt\;\frac{e^{-t}}{t}$
as the normalized exponential integral.

Notice that $s(y,n)=\frac{y_{1}}{\left|n\cdot N\right|}$ if $n\cdot N<0$. We
can summarize the equation the temperature satisfies in the non-Grey
approximation as follows

$\begin{split}&\int_{0}^{\infty}d\nu\;\alpha_{\nu}(0)B_{\nu}\left(T(y_{1})\right)-\int_{0}^{\infty}d\nu\int_{0}^{\infty}d\eta\;\alpha_{\nu}^{2}(0)K\left(\alpha_{\nu}(0)\left|y_{1}-\eta_{1}\right|\right)B_{\nu}\left(T\left(\eta_{1}\right)\right)\\\
=&\int_{0}^{\infty}d\nu\int_{n\cdot
N<0}dn\;\overline{g}_{\nu}(n)\alpha_{\nu}(0)e^{-\alpha_{\nu}(0)\frac{y_{1}}{\left|n\cdot
N\right|}}.\end{split}$ (2.6)

In the particular case of the Grey approximation when $\alpha\equiv 1$ using
that $u(y)=4\pi\sigma T^{4}(y)$ we can simplify equation (2.6) by property
(1.3)

$u(y_{1})-\int_{0}^{\infty}d\eta\;K(y_{1}-\eta)u(\eta)=\int_{0}^{\infty}d\nu\int_{n\cdot
N<0}dn\;\overline{g}_{\nu}(n)e^{-\frac{y_{1}}{\left|n\cdot N\right|}}.$ (2.7)

In some occasions, when the dependence of the boundary layer function $u$ on
the point $p\in\partial\Omega$ is needed, we will use the notation
$\overline{u}(y_{1},p)$, where this function solves according to the rigid
motion $\mathcal{R}_{p}$ in (1.12)

$\overline{u}(y_{1},p)-\int_{0}^{\infty}d\eta\;K(y_{1}-\eta)u(\eta,p)=\int_{0}^{\infty}d\nu\int_{n\cdot
N_{p}<0}dn\;g_{\nu}(n)e^{-\frac{y_{1}}{\left|n\cdot N_{p}\right|}}.$ (2.8)

For the rest of Section 2 and Section 3 we will focus on the study of
$\overline{u}(y_{1},p)$ for an arbitrary given $p\in\partial\Omega$, hence we
will call $u(y_{1})=\overline{u}(y_{1},p)$ and $N=N_{p}$. In order to simplify
the reading from now on we set $G(x)=\int_{0}^{\infty}d\nu\int_{n\cdot
N<0}dn\;\overline{g}_{\nu}(n)e^{-\frac{x}{\left|n\cdot
N\right|}}\raisebox{2.0pt}{$\chi$}_{\\{x>0\\}}$ and if we want to stress out
the dependence on $p\in\partial\Omega$ we write
$G_{p}(x)=\int_{0}^{\infty}d\nu\int_{n\cdot
N_{p}<0}\;g_{\nu}(n)e^{-\frac{x}{\left|n\cdot
N_{p}\right|}}\raisebox{2.0pt}{$\chi$}_{\\{x>0\\}}$.

From now on until Section 5 we consider the case of constant absorption
coefficient $\alpha\equiv 1$.

###  2.3 Some properties of the kernel

We consider the kernel $K$ introduced in Section 2.2. We remark that
$K(x)=\frac{1}{2}E_{1}(|x|)$, where $E_{1}$ is the standard exponential
integral function. See [1]. We collect some properties of the normalized
exponential integral.

######  Proposition 2.1.

The function $K$ satisfies $\int_{-\infty}^{\infty}dx\;K(x)=1$, $K\in
L^{{1}}\left(\mathbb{R}\right)\cap L^{{2}}\left(\mathbb{R}\right)$ and the
following estimate holds

$\frac{1}{4}e^{-|x|}\ln(1+\frac{2}{|x|})\leq
K(x)\leq\frac{1}{2}e^{-|x|}\ln(1+\frac{1}{|x|}).$

Moreover, the Fourier transform of $K$ is
$\hat{K}(\xi)=\frac{1}{\sqrt{2\pi}}\frac{\arctan(\xi)}{\xi}$.

###### Proof.

Since $K$ is even and non negative we can calculate, applying Tonelli’s
Theorem

$\int_{-\infty}^{\infty}K(s)\;ds=2\int_{0}^{\infty}K(s)\;ds=\int_{0}^{\infty}\int_{s}^{\infty}\frac{e^{-t}}{t}\;dt\;ds=\int_{0}^{\infty}\frac{e^{-t}}{t}\int_{0}^{t}\;ds\;dt=\int_{0}^{\infty}e^{-t}\;dt=1.$

This proves also that $K\in L^{{1}}\left(\mathbb{R}\right)$.

For the square integrability we refer to equation 5.1.33 in [1] and see
$\int_{\mathbb{R}}\left|K(x)\right|^{2}\;dx=\ln(2)$. Estimate 5.1.20 in [1]
also implies $\frac{1}{4}e^{-|x|}\ln(1+\frac{2}{|x|})\leq
K(x)\leq\frac{1}{2}e^{-|x|}\ln(1+\frac{1}{|x|})$.

We now move to the computation of the Fourier transform of the kernel $K$. The
kernel is an even function, hence we compute

$\begin{split}\hat{K}(\xi)=&\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{-i\xi
x}K(x)\;dx=\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}\cos\left(\xi
x\right)\int_{x}^{\infty}\frac{e^{-t}}{t}\;dt\;dx\\\
=&\frac{1}{\sqrt{2\pi}}\frac{1}{\xi}\int_{0}^{\infty}\frac{e^{-t}}{t}\sin(\xi
t)\;dt=\frac{1}{\sqrt{2\pi}}\frac{\arctan(\xi)}{\xi}.\end{split}$

The last identity can be justified noticing that
$F(\xi)=\int_{0}^{\infty}\frac{e^{-t}}{t}\sin(\xi t)\;dt$ has derivative
$F^{\prime}(\xi)=\frac{1}{\xi^{2}+1}$. ∎

The following calculation will also be very useful in the next section.

######  Proposition 2.2.

Let $x>0$. Then we can compute

$\int_{-x}^{\infty}K(s)\;ds=1-\frac{e^{-x}}{2}+xK(x);$ (2.9)
$\int_{x}^{\infty}K(s)\;ds=\frac{e^{-x}}{2}-xK(x);$ (2.10)
$\int_{-x}^{\infty}sK(s)\;ds=\int_{x}^{\infty}sK(s)\;ds=\frac{xe^{-x}}{4}+\frac{e^{-x}}{4}-\frac{x^{2}}{2}K(x);$
(2.11)

###### Proof.

The proof relies on basic integral computations. We have to compute several
integrals changing the order of integration applying Tonelli’s Theorem and
integrating by parts. We assume $x>0$. We prove only (2.9), since all other
formulas can be obtained in a similar way.

$\begin{split}\int_{-x}^{\infty}K(s)\;ds=&\frac{1}{2}\int_{-x}^{0}\int_{|s|}^{\infty}\frac{e^{-t}}{t}\;dt\;ds+\frac{1}{2}\int_{0}^{\infty}\int_{s}^{\infty}\frac{e^{-t}}{t}\;dt\;ds\\\
=&\frac{1}{2}\int_{0}^{x}\frac{e^{-t}}{t}\int_{0}^{t}\;ds\;dt+\frac{1}{2}\int_{x}^{\infty}\frac{e^{-t}}{t}\int_{0}^{x}\;ds\;dt+\frac{1}{2}\\\
=&1-\frac{e^{-x}}{2}+xK(x).\end{split}$

∎

##  3 The boundary condition for the limit problem

We now start with the boundary layer analysis. This boundary layer problem,
know in the literature as Milne problem, was studied with different
approaches, e.g. [2, 3, 7, 10, 19]. We will present another proof of the
boundary layer analysis for the equation (2.7) of the temperature. The proof
uses a combination of comparison arguments and Fourier analysis. In addition,
instead of considering the intensity of radiation the analysis is made
directly for the temperature.

Our aim is now to solve equation (2.7). Indeed, according to the method of
matched asymptotic expansions we expect the boundary condition for the limit
problem to be the limit of $u$ as $y\to\infty$ for every point
$x\in\partial\Omega$. In order to simplify the notation we call
$\mathcal{L}\left(u\right)(x):=u(x)-\int_{0}^{\infty}dy\;K\left(x-y\right)u(y)$
and
$\overline{\mathcal{L}}\left(u\right)(x):=u(x)-\int_{-\infty}^{\infty}dy\;K\left(x-y\right)u(y)$.

###  3.1 The homogeneous equation

We start with the study of the homogeneous equation, i.e. (2.7) with
$G(x)\equiv 0$. We will show using maximum principle that any bounded solution
is the trivial solution $u\equiv 0$. We will use the following version of the
maximum principle for the non-local operator $\mathcal{L}$.

######  Lemma 3.1.

Let $\overline{u}\in C\left([0,\infty)\right)$ with
$\lim\limits_{x\to\infty}\overline{u}(x)\in[0,\infty]$ be a supersolution of
(2.7), i.e.

$\begin{cases}\overline{u}(x)-\int_{0}^{\infty}dy\;K\left(x-y\right)\overline{u}(y)\geq
0&x>0\\\ \overline{u}(x)=0&x<0\end{cases}$

Then $u\geq 0$ for all $x\geq 0$.

###### Proof.

Let us assume the contrary, i.e. that there exists some $x\in[0,\infty]$ such
that $\overline{u}(x)<0$. By assumption $x\in[0,\infty)$. Since $\overline{u}$
is continuous in $[0,\infty)$ and it has non-negative limit at infinity which
is bounded or infinity, $u$ attains its global minimum in $[0,\infty)$, i.e.
there exists some $x_{0}\in[0,\infty)$ such that
$\overline{u}(x_{0})=\inf_{x\in[0,\infty)}\overline{u}(x)<0$. Since
$\overline{u}$ is a super solution we can calculate

$\begin{split}0\leq&\mathcal{L}(\overline{u})(x_{0})=\overline{u}(x_{0})-\int_{0}^{\infty}dy\;K\left(x_{0}-y\right)\overline{u}(y)\\\
=&\int_{-\infty}^{\infty}dy\;K\left(x_{0}-y\right)\overline{u}(x_{0})-\int_{0}^{\infty}dy\;K\left(x_{0}-y\right)\overline{u}(y)\\\
=&\int_{-\infty}^{0}dy\;K\left(x_{0}-y\right)\overline{u}(x_{0})+\int_{0}^{\infty}dy\;K\left(x_{0}-y\right)\left(\overline{u}(x_{0})-\overline{u}(y)\right)\\\
<&0,\end{split}$

where we used the positivity of $K\left(x_{0}-y\right)$, the fact that the
integral of the kernel $K$ is $1$ and the fact that $\overline{u}(x_{0})$ is
the minimum of $\overline{u}$ and it is strictly negative. This leads to a
contradiction and thus we conclude the proof. ∎

With the maximum principle we can now show the following theorem on the
triviality of the solution to the homogeneous equation.

######  Theorem 3.1.

Assume $u$ is a bounded solution to

$\overline{\mathcal{L}}(u)(x)=0$ (3.1)

with $u(x)\equiv 0$ for $x<0$. Then $u=0$ for almost every $x\in\mathbb{R}$.

###### Proof.

We will construct a supersolution $\overline{u}$ which converges to infinity
and we will apply Lemma 3.1 to the supersolutions $\overline{u}-u$ and
$u+\overline{u}$. First of all we see that for $x>0$ the bounded solution $u$
is continuous, indeed $u(x)=K*u(x)$. Since $K\in
L^{{1}}\left(\mathbb{R}\right)$ and $u\in L^{{\infty}}\left(\mathbb{R}\right)$
then the convolution is a continuous bounded function. Moreover we can extend
continuously $u$ in $0$. Indeed, we define

$u(0)=\lim\limits_{x\to
0}\left[G(x)+\int_{0}^{\infty}dy\;K\left(x-y\right)u(y)\right].$

This limit exists because $G$ is continuous in $[0,\infty)$ and for the
integral term we can apply the generalized dominated convergence theorem using
that the sequence $K\left(x-y\right)\to K\left(y\right)$ as $x\to 0$ pointwise
and in $L^{{1}}\left(\mathbb{R}\right)$.

We consider now the function

$\overline{u}(x)=\begin{cases}1+x&x\geq 0\\\ 0&x<0\end{cases}$

$\overline{u}$ is a supersolution. It is indeed possible to calculate
$\mathcal{L}(\overline{u})(x)$. Let $x\geq 0$. Then
$\mathcal{L}(\overline{u})=\mathcal{L}(Id)+\mathcal{L}(1)$. By a simple
calculation we get on the one hand

$\begin{split}\mathcal{L}(Id)(x)=&x-\int_{0}^{\infty}dy\;K\left(x-y\right)y=x-\int_{-x}^{\infty}dy\;K\left(y\right)(x+y)=\frac{x}{4}e^{-x}-\frac{e^{-x}}{4}-\frac{x^{2}}{2}K(x)\end{split}$

and on the other hand

$\begin{split}\mathcal{L}(1)(x)=&1-\int_{0}^{\infty}dy\;K\left(x-y\right)=1-\int_{-x}^{\infty}dy\;K\left(y\right)=\frac{e^{-x}}{2}-xK(x).\end{split}$

Therefore we want to show that the function
$f(x):=\mathcal{L}(\overline{u})(x)=\frac{e^{-x}}{4}(1+x)-\frac{x}{2}K(x)(2+x)$
is non-negative for all $x\geq 0$. It is not difficult to see that
$f(0)=\frac{1}{4}>0$ and that $\lim\limits_{x\to\infty}f(x)=0$. Moreover, we
can consider the derivative

$\begin{split}f^{\prime}(x)=&\frac{1}{2}\left(e^{-x}-K(x)(2x+2)\right)\leq\frac{1}{2}\left(e^{-x}-\frac{e^{-x}}{2}\ln\left(1+\frac{2}{x}\right)(x+1)\right)\leq
0.\end{split}$

The first inequality is given by the estimate of Proposition 2.1 and the
second one is due to the well-know estimate $\ln(1+x)\geq\frac{2x}{2+x}$. The
non-positivity of the derivative implies that $f$ is monotonously decreasing,
and therefore $\mathcal{L}(\overline{u})(x)=f(x)\geq 0$ for all $x\geq 0$.

Let now $\varepsilon>0$ arbitrary. We know that $u$ is bounded and
$\overline{u}$ converges to infinity, moreover both $u$ and $\overline{u}$ are
continuous in $[0,\infty)$. Also $u$ is a homogeneous solution of (2.7) and
the operator $\mathcal{L}$ is linear. Therefore we can apply Lemma 3.1 to the
supersolutions $\varepsilon\overline{u}-u$ and $u+\varepsilon\overline{u}$ and
get that the
$\inf_{x\in[0,\infty)}\left[\varepsilon\overline{u}(x)-u(x)\right]\geq 0$ and
$\inf_{x\in[0,\infty)}\left[\varepsilon\overline{u}(x)+u(x)\right]\geq 0$.
This implies that for any $x\in\mathbb{R}$ the following holds

$-\varepsilon\overline{u}(x)\leq u(x)\leq\varepsilon\overline{u}(x)$

Since $\varepsilon$ was arbitrary we conclude $u(x)=0$ for all
$x\in\mathbb{R}$. ∎

###  3.2 Well-posedness theory for the inhomogeneous equation

We can now move to the well-posedness theory for the inhomogeneous equation,
for which the next theorem is the main result.

######  Theorem 3.2.

Let $H:\mathbb{R}_{+}\to\mathbb{R}_{+}$ be a continuous function bounded by an
exponential function, i.e. $|H(x)|\leq
Ce^{-Ax}\raisebox{2.0pt}{$\chi$}_{\\{x>0\\}}$ for $C,A>0$. Then there exists a
unique bounded solution to the equation

$\begin{cases}u(x)-\int_{0}^{\infty}dy\;K\left(x-y\right)u(y)=H(x)&x>0,\\\
u(x)=0&x<0.\end{cases}$ (3.2)

Moreover, $u$ is continuous on $(0,\infty)$.

###### Proof.

The assumption on the exponential decay of $H$ yields $H\in
L^{1}\left(\mathbb{R}\right)\cap L^{2}\left(\mathbb{R}\right)\cap
L^{\infty}\left(\mathbb{R}\right)$. In order to find a bounded solution for
(3.2) we will follow several steps. We will look for functions $\tilde{u}$ and
$v$ solutions of the following equations

$\tilde{u}(x)-\int_{-\infty}^{\infty}K(x-y)\tilde{u}(y)\;dy=\overline{H}(x):=H(x)-H(-x)\;\;\;\;\;\;\;\;x\in\mathbb{R}$

and

$\begin{cases}v(x)-\int_{-\infty}^{\infty}K(x-y)v(v)\;dy=0,&x>0\\\
v(x)=-\tilde{u}(x)&x<0.\end{cases}$ (3.3)

Then $u=\tilde{u}+v$ will be the desired solution.

Step 1: Construction of $\tilde{u}$.
We can construct the solution $\tilde{u}$ via Fourier method. First of all we
notice that any affine function is a solution to the homogeneous equation in
the whole space $\mathbb{R}$. This is because
$\int_{-\infty}^{\infty}K(x)\;dx=1$ and $\int_{-\infty}^{\infty}xK(x)\;dx=0$.
Since by assumption $H\in L^{{2}}\left(\mathbb{R}\right)$ also
$\overline{H}\in L^{{2}}\left(\mathbb{R}\right)$. We define for an integrable
function $f$ the kth-moment as
$m_{k}\left(f\right)=\int_{-\infty}^{\infty}x^{k}f(x)\;dx$ assuming it exists.
Then clearly by construction $m_{0}\left(\overline{H}\right)=0$ while
$m_{1}\left(\overline{H}\right)=\frac{2C}{A}>0$. Moreover, since
$\overline{H}$ has exponential decay, all moments
$m_{k}\left(\overline{H}\right)<\infty$ are bounded.

We define also the function
$F(x)=\overline{\mathcal{L}}\left(\frac{3}{2}\text{sgn}\right)(x)$. It can be
compute that

$F(x)=\frac{3}{2}\left(\text{sgn}(x)-\int_{-\infty}^{\infty}K(x-y)\text{sgn}(y)\;dy\right)=\frac{3}{2}\left(\text{sgn}(x)-\int_{-x}^{x}K(y)\;dy\right).$

It is not difficult to see that $F(0)=0$, $\lim\limits_{|x|\to\infty}F(x)=0$
and that $F$ is a stepwise continuous function with the discontinuity in $0$.
Therefore $F(x)$ is bounded. We proceed with the construction of $\tilde{u}$.
We can write it as $\tilde{u}=u^{(1)}+u^{(2)}+a+bx$, where
$u^{(1)}(x)=m_{1}\left(\overline{H}\right)\frac{3}{2}\text{sgn}(x)$ solves the
equation

$\mathcal{L}\left(u^{(1)}\right)(x)=m_{1}\left(\overline{H}\right)F(x)\;\;\;\;\;\;\;\;\;x\in\mathbb{R}$

and $u^{(2)}$ solves

$\mathcal{L}\left(u^{(2)}\right)(x)=\overline{H}(x)-m_{1}\left(\overline{H}\right)F(x)\;\;\;\;\;\;\;\;\;x\in\mathbb{R}.$
(3.4)

applying now the Fourier transform to the equation (3.4), recalling the
convolution rule and the Fourier transforms of the kernel $K$ and of the sgn
function we get first in distributional sense

$\hat{u}^{(2)}(s)\left(\frac{s-\arctan(s)}{s}\right)=\mathcal{F}(\overline{H})(s)+\frac{3m_{1}\left(\overline{H}\right)}{\sqrt{2\pi}}\frac{i}{s}\frac{s-\arctan(s)}{s}.$
(3.5)

The Fourier transform of $\overline{H}$ is in $C^{\infty}$, since
$\overline{H}$ has exponential decay and therefore it has all kth-moment
finite. Therefore there exists a function $\tilde{H}$ such that
$\tilde{H}(0)=\tilde{H}^{\prime}(0)=0$ such that
$\mathcal{F}(\overline{H})(s)=-\frac{i}{\sqrt{2\pi}}m_{1}\left(\overline{H}\right)s+\tilde{H}(s)$,
since $m_{0}\left(\overline{H}\right)=0$ and by definition
$\mathcal{F}(\overline{H})^{\prime}(s)\big{|}_{s=0}=-\frac{i}{\sqrt{2\pi}}m_{1}\left(\overline{H}\right)$.
We can therefore find first formally $u^{(2)}$ analyzing its Fourier transform

$\begin{split}\hat{u}^{(2)}(s)=&\frac{s}{s-\arctan(s)}\mathcal{F}(\overline{H})(s)+\frac{3m_{1}\left(\overline{H}\right)}{\sqrt{2\pi}}\frac{i}{s}\\\
=&-\frac{is^{2}}{s-\arctan(s)}\frac{m_{1}\left(\overline{H}\right)}{\sqrt{2\pi}}+\frac{3m_{1}\left(\overline{H}\right)}{\sqrt{2\pi}}\frac{i}{s}+\tilde{H}(s)\frac{s}{s-\arctan(s)}\\\
=&\mathbb{H}(s).\end{split}$ (3.6)

It is important to notice that $\lim\limits_{s\to
0}\frac{s^{2}}{s-\arctan(s)}-\frac{3}{s}=0$, since
$\frac{s}{s-\arctan(s)}=\frac{3}{s^{2}}+\frac{9}{5}+O(s^{2})$ near zero. Using
L’Hôpital rule we see also that $\lim\limits_{s\to
0}\tilde{H}(s)\frac{s}{s-\arctan(s)}$ is finite. On the other hand
$\frac{s}{s-\arctan(s)}$ is bounded for $|s|>1$. Since
$\mathcal{F}(\overline{H})(s)$ and $\frac{1}{s}$ are both square integrable
functions and since $\mathbb{H}$ is bounded near $0$ we conclude that
$\mathbb{H}\in L^{{2}}\left(\mathbb{R}\right)$. Therefore also the in (3.6)
defined $\hat{u}^{(2)}$ is square integrable. We can hence invert it

$u^{(2)}(x):=\mathcal{F}^{-1}\left(\mathbb{H}\right)(x)\in
L^{{2}}\left(\mathbb{R}\right).$

Since this function solves (3.5) not only in distributional sense but also
pointwise almost everywhere, we can conclude rigorously that indeed the
function in (3.6) is the desired $u^{(2)}$ solving (3.4). Moreover,
$u^{(2)}=K*u^{(2)}+\overline{H}-F$ and since both $K$ and $u^{(2)}$ itself are
square integrable and both $H$ and $F$ are bounded, then also $u^{(2)}$ is
bounded. We can conclude this step therefore defining

$\tilde{u}(x)=\frac{3}{2}m_{1}\left(\overline{H}\right)\text{sgn}(x)+a+bx+u^{(2)}(x).$
(3.7)

Step 2: Construction of $v$.

We recall that the equation $v$ shall solve (3.3) As we found out in the first
step,
$\tilde{u}=\frac{3}{2}m_{1}\left(\overline{H}\right)\text{sgn}(x)+a+bx+u^{(2)}(x)$.
As we already pointed out, affine solutions are always solution of the
homogeneous equation in the whole space $\mathbb{R}$. Therefore, we shall look
for a function of the form

$v(x)=\frac{3}{2}m_{1}\left(\overline{H}\right)-a-bx+v^{(2)}(x)$ (3.8)

where $v^{(2)}$ solves similarly as above

$\begin{cases}v^{(2)}(x)-\int_{-\infty}^{\infty}K(x-y)v(v)\;dy=0&x>0,\\\
v^{(2)}(x)=-u^{(2)}(x)&x<0.\end{cases}$ (3.9)

We proceed now iteratively constructing the desired solution. We call $B>0$
the constant such that $\left\Arrowvert u^{(2)}\right\Arrowvert_{\infty}\leq
B$ and we define $\overline{v}=B$ and $\underline{v}=-B$. Inductively we
define $v_{0}:=\underline{v}$ and for $k\geq 1$ we set

$v_{k}(x)=\begin{cases}-u^{(2)}(x)&x<0,\\\
\int_{-\infty}^{\infty}K\left(x-y\right)v_{k-1}(y)\;dy&x>0.\\\ \end{cases}$

We claim that $\underline{v}=v_{0}\leq v_{1}\leq v_{2}\leq...\leq v_{k}\leq
v_{k+1}\leq...$ and that $v_{k}\leq\underline{v}$ for all $k\in\mathbb{N}$.
Clearly for $k=0$ both statements hold. On the one hand since
$\int_{-\infty}^{\infty}K\left(x-y\right)v_{0}(y)=-B$ we see that

$v_{1}(x)-v_{0}(x)=\begin{cases}-u^{(2)}(x)+B\geq 0&x<0,\\\ 0&x>0,\\\
\end{cases}$

on the other hand per definition we have $\overline{v}-\underline{v}=2B\geq
0$. We see also that $\overline{v}-v_{1}\geq 0$, indeed

$\overline{v}(x)-v_{1}(x)=\begin{cases}B+u^{(2)}(x)\geq 0&x<0,\\\ 2B&x>0.\\\
\end{cases}$

We now prove inductively that $v_{k}\geq v_{k-1}$ and $\overline{v}\geq
v_{k}$. Hence, we assume that these inequalities are satisfied for $k$ and we
prove them for $k+1$. Indeed this just follows from the identities

$v_{k+1}(x)-v_{k}(x)=\begin{cases}0&x<0,\\\
\int_{-\infty}^{\infty}K\left(x-y\right)\left(v_{k}(y)-v_{k-1}(y)\right)\geq
0&x>0,\\\ \end{cases}$
$\overline{v}(x)-v_{k+1}(x)=\begin{cases}B+u^{(2)}(x)\geq 0&x<0,\\\
\int_{-\infty}^{\infty}K\left(x-y\right)\left(B-v_{k}(y)\right)\geq 0&x>0,\\\
\end{cases}$

where we used again that the integral in the whole line of the kernel $K$ is
$1$. Therefore the sequence $v_{k}(x)$ is increasing and bounded. This means
that there exists a pointwise limit. By the dominated convergence theorem and
by construction this will be also the desired solution of (3.9), i.e.

$v^{(2)}(x):=\lim\limits_{k\to\infty}v_{k}(x)$

solves the equation (3.9) and it is by construction bounded.

Step 3: Properties of $u$.

Now we are ready to write down the whole solution. As we remarked at the
beginning $u=\tilde{u}+v$, where $\tilde{u}$ solves as in Step 1 (3.2) and $v$
solves as in Step 2 (3.3). Therefore by (3.8) and by (3.7)

$u(x)=\begin{cases}6m_{1}\left(H\right)+u^{(2)}(x)+v^{(2)}(x)&x>0,\\\
0&x<0,\end{cases}$

solves the initial problem (3.2) and it is by construction bounded. Moreover,
since $K$ is integrable and $H$ is continuous in $[0,\infty)$ also $u=K*u+H$
is continuous in $[0,\infty)$.

Step 4: Uniqueness.

Let us assume that $u_{1}$ and $u_{2}$ are two bounded solution to the problem
(3.2). Then $u_{1}-u_{2}$ will be a bounded continuous solution to the
homogeneous problem (3.1). Therefore by Theorem 3.1 $u_{1}-u_{2}=0$. Hence,
there exists a unique bounded solution $u$ to the inhomogeneous problem (3.2).
This concludes the proof. ∎

######  Corollary 3.1.

Let $p\in\partial\Omega$ and $G_{p}(x)$ as defined in (2.8). Let
$g_{\nu}(n)\geq 0$ and assume $\int_{0}^{\infty}d\nu\;g_{\nu}(n)\in
L^{\infty}\left(\mathbb{S}\right)$. Then there exists a unique bounded
solution to the equation

$\begin{cases}u(x)-\int_{0}^{\infty}dy\;K\left(x-y\right)u(y)=G_{p}(x)&x>0,\\\
u(x)=0&x<0.\end{cases}$ (3.10)

Moreover, $u$ is continuous on $(0,\infty)$.

###### Proof.

By assumption $G_{p}$ is continuous for $x>0$ and $|G_{p}(x)|\leq\Arrowvert
g\Arrowvert_{1}e^{-y}\raisebox{2.0pt}{$\chi$}_{\\{x>0\\}}$. Hence we can apply
Theorem 3.2. ∎

It is also possible to show, that the bounded solution $u$ is non-negative

######  Lemma 3.2.

Let $u$ be the unique bounded solution to (3.10). Then $u(x)\geq 0$ for all
$x\in\mathbb{R}$.

###### Proof.

The proof is very similar to the proof of Theorem 3.1. We consider the
supersolution

$\overline{u}(x)=\begin{cases}1+x&x\geq 0\\\ 0&x<0\end{cases}.$

As we have seen before, $u=K*u+G$ is continuous in $[0,\infty)$. Moreover,
since $G>0$ as $x\geq 0$, $u$ is a supersolution too. Let now $\varepsilon>0$
be arbitrary. Let us consider the supersolution $\varepsilon\overline{u}+u$.
This is continuous in $[0,\infty)$ and since $u$ is bounded it converges to
infinity as $x\to\infty$. Therefore Lemma 3.1 implies that there exists no
$x_{0}\in[0,\infty)$ such that

$\inf_{x\in[0,\infty)}\left(\varepsilon\overline{u}(x)+u(x)\right)=\varepsilon\overline{u}(x_{0})+u(x_{0})<0.$

Hence $u\geq-\varepsilon\overline{u}$ and since $\varepsilon>0$ was arbitrary
we conclude $u\geq 0$. ∎

######  Remark.

Theorem 3.2 can be proved also using the Wiener-Hopf solution formula for the
problem (2.7) as given in [17]. It is true that in this way one obtains an
explicit formula, which not only assures the well-posedness of the planar
problem we are studying but also directly shows the existence of a limit for
the solution $u$ when $x\to\infty$. However, the Wiener-Hopf method produces a
complicate formula which requires a careful analysis with complex variables in
order to be understood. We have preferred to use this soft method approach
which in particular allows us to prove some relevant properties of the
solution, such as the positivity.

###  3.3 Asymptotic behavior of the bounded solution of the inhomogeneous
equation

We were able to show that the equation for the boundary value in the Grey
approximation has a unique bounded solution which is positive whenever $G>0$.
As we anticipated at the beginning of this section, we would like to study the
limit as $x\to\infty$ of the solution $u(x)$. We will show, that such limit
exists and is uniquely characterized by $g_{\nu}(n)$ and $N$. To this end we
first prove that the function $u$ is uniformly continuous.

######  Lemma 3.3.

Let $u$ be the unique bounded solution to the problem (2.7). The $u$ is
uniform continuous on $[0,\infty)$ and it satisfies for $x,y\in[0,\infty)$

$\begin{split}\left|u(x)-u(y)\right|\leq&\left|G(x)-G(y)\right|\\\
+&\left\Arrowvert
u\right\Arrowvert_{\infty}\left[\frac{\left|e^{-x}-e^{-y}\right|}{2}+2\left(1-e^{\frac{\left|x-y\right|}{2}}\right)+4\left|\frac{y-x}{2}\right|K\left(\frac{y-x}{2}\right)+\left|xK(x)-yK(y)\right|\right].\end{split}$
(3.11)

###### Proof.

This is a consequence of the uniform continuity of $G$ and $xK(x)$. Clearly,
since $u$ solves the problem (3.10), we have the estimate

$\left|u(x)-u(y)\right|\leq\left|G(x)-G(y)\right|+\int_{0}^{\infty}\left|K\left(\eta-x\right)-K\left(\eta-y\right)\right|u(\eta)\;d\eta.$
(3.12)

Since $G$ is continuous on $[0,1]$, and therefore uniformly continuous on
$[0,1]$ and since $G$ is Lipschitz continuous in $[1,\infty)$, $G$ is uniform
continuous in $[0,\infty)$. The latter affirmation is true, since

$\sup_{x\geq 1}\left|G^{\prime}(x)\right|\leq\int_{0}^{\infty}d\nu\int_{n\cdot
N<0}dn\;g_{\nu}(n)\frac{e^{-\frac{1}{|n\cdot N|}}}{|n\cdot N|}<\infty,$

where the finiteness is due to the fact that $\lim\limits_{|n\cdot N|\to
0}\frac{e^{-\frac{1}{|n\cdot N|}}}{|n\cdot N|}=0$.

For the integral term in (3.12) we assume that $x<y$. Then we can calculate
using the fact that for positive arguments the kernel $K$ is decreasing

$\begin{split}\int_{0}^{\infty}&\left|K\left(\eta-x\right)-K\left(\eta-y\right)\right|u(\eta)\;d\eta\\\
=&\int_{0}^{\frac{x+y}{2}}\left(K\left(\eta-x\right)-K\left(\eta-y\right)\right)u(\eta)\;d\eta+\int_{\frac{x+y}{2}}^{\infty}\left(K\left(\eta-y\right)-K\left(\eta-x\right)\right)u(\eta)\;d\eta\\\
\leq&\left\Arrowvert
u\right\Arrowvert_{\infty}\left[\int_{0}^{\frac{x+y}{2}}\left(K\left(\eta-x\right)-K\left(\eta-y\right)\right)\;d\eta+\int_{\frac{x+y}{2}}^{\infty}\left(K\left(\eta-y\right)-K\left(\eta-x\right)\right)\;d\eta\right]\\\
\end{split}$

We can calculate explicitly the last two integrals using the result of
Proposition 2.2, indeed by a change of variable

$\begin{split}\int_{0}^{\infty}&\left|K\left(\eta-x\right)-K\left(\eta-y\right)\right|u(\eta)\;d\eta\\\
\leq&\left\Arrowvert
u\right\Arrowvert_{\infty}\left[\int_{-x}^{\frac{y-x}{2}}K\left(\eta\right)\;d\eta-\int_{-y}^{\frac{x-y}{2}}K\left(\eta\right)\;d\eta+\int_{\frac{x-y}{2}}^{\infty}K\left(\eta\right)\;d\eta-\int_{\frac{y-x}{2}}^{\infty}K\left(\eta\right)\;d\eta\right]\\\
=&\left\Arrowvert
u\right\Arrowvert_{\infty}\left[\frac{e^{-y}-e^{-x}}{2}+2\left(1-e^{\frac{x-y}{2}}\right)+4\frac{y-x}{2}K\left(\frac{y-x}{2}\right)+xK(x)-yK(y)\right].\end{split}$

Recalling that $x<y$ we get the estimate (3.11). From the well-known estimates
$\left|e^{-x}-e^{-y}\right|\leq\left|x-y\right|$ and
$\left|1-e^{\frac{x-y}{2}}\right|\leq\frac{\left|x-y\right|}{2}$ we see that
we shall only consider the function $f(x)=xK(x)$. Since $f(0)=0$ and $f$ is
continuous, $f$ is uniformly continuous on $[0,1]$, on the other hand $f$ is
Lipschitz continuous on $[1,\infty]$. This is because

$\sup\limits_{x\geq 1}\left|f^{\prime}(x)\right|=\sup\limits_{x\geq
1}\left|K(x)-e^{-x}\right|\leq\frac{1}{e}+K(1)<\infty.$

Therefore $f$ is uniform continuous on $[0,\infty)$. By the continuity of $f$
in $0$ we also now that given an $\varepsilon>0$ there exists some $\delta$
such that $\frac{y-x}{2}K\left(\frac{y-x}{2}\right)<\varepsilon$ for all
$\left|x-y\right|<\delta$. Hence, we conclude that $u$ is uniform continuous.
∎

We want now to show that the limit $\lim\limits_{y\to\infty}u(y)$ exists. To
this end we proceed again using Fourier methods.

######  Theorem 3.3.

Let $u$ be the unique bounded solution to the problem (3.10). Then
$\lim\limits_{x\to\infty}u(x)$ exists and it is uniquely determined by $G$ and
$u$ itself. Moreover, the limit is positive if $\left\\{n\in\mathbb{S}:n\cdot
N<0\text{ and }\int_{0}^{\infty}d\nu\overline{g}_{\nu}(n)\not\equiv
0\right\\}$ is not a zero measure set.

###### Proof.

Since $u$ is the unique bounded solution, $u$ solves for all $x\in\mathbb{R}$

$u(x)-\int_{-\infty}^{\infty}K(y-x)u(y)\;dy=G(x)\;\raisebox{2.0pt}{$\chi$}_{\\{x>0\\}}-\int_{0}^{\infty}K(y-x)u(y)\;dy\;\raisebox{2.0pt}{$\chi$}_{\\{x<0\\}}\equiv
W(x).$ (3.13)

Indeed, this is equivalent to (3.10). This can be seen easily, since $u$
solves for $x<0$

$u(x)-\int_{-\infty}^{0}K(y-x)u(y)\;dy=0$

and since $u=0$ for $x<0$ is a possible solution, by uniqueness, this is the
only possible solution. It is not only true that $W\in
L^{1}\left(\mathbb{R}\right)\cap L^{2}\left(\mathbb{R}\right)$ but also that
$W$ has all moments bounded. This follows from the similar property of $G$
(cf. Step 1 in Theorem 3.2) as well as from the inequality
$0\leq\int_{0}^{\infty}K(y-x)u(y)\;dy\raisebox{2.0pt}{$\chi$}_{\\{x<0\\}}\leq\Arrowvert
u\Arrowvert_{\infty}\;\raisebox{2.0pt}{$\chi$}_{\\{x<0\\}}\left(\frac{e^{-|x|}}{2}-|x|K(x)\right)$.
Notice that $|x|K(x)\leq\frac{e^{-}|x|}{2}$. Hence, finite moments and
Riemann-Lebesgue Theorem imply that $W$ has a Fourier transform $\hat{W}\in
C_{0}\left(\mathbb{R}\right)\cap C^{\infty}\left(\mathbb{R}\right)\cap
L^{2}\left(\mathbb{R}\right)$. Moreover, looking at the left hand side of
(3.13) we recall as in [28] that in distributional sense for all
$\phi\in\mathcal{S}\left(\mathbb{R}\right)$

$\langle\hat{u}-\mathcal{F}\left(u*K\right),\phi\rangle:=\langle
u-u*K,\hat{\phi}\rangle=\langle
u,\mathcal{F}\left((1-\sqrt{2\pi}\hat{K})\phi\right)\rangle,$

where the last equality is due to an elementary calculation involving the
convolution and we define $\langle f,g\rangle=\int_{\mathbb{R}}f(x)g(x)\;dx$.
We recall also that
$1-\sqrt{2\pi}\hat{K}(\xi)=\frac{\xi-\arctan(\xi)}{\xi}:=F(\xi)$. Hence, for
all $\phi\in\mathcal{S}\left(\mathbb{R}\right)$ we have

$\langle u,\mathcal{F}(\phi F)\rangle=\langle\hat{W},\phi\rangle.$ (3.14)

Now we consider for $\varepsilon>0$ the sequence of standard mollifiers
$\phi_{\varepsilon}(\xi):=\frac{1}{\varepsilon}\phi\left(\frac{\xi}{\varepsilon}\right)\in
C_{c}^{\infty}\left(\mathbb{R}\right)\subset\mathcal{S}\left(\mathbb{R}\right)$
such that in distributional sense $\phi_{\varepsilon}\rightharpoonup\delta$.
The smoothness of $\hat{W}$ implies $\langle\hat{W},\phi\rangle\to\hat{W}(0)$
as $\varepsilon\to 0$. It is our first aim to show that $\hat{W}(0)$ is zero.
To this end we study the left hand side of (3.14). We calculate

$\begin{split}\langle
u,&\mathcal{F}(\phi_{\varepsilon}F)\rangle=\frac{1}{\sqrt{2\pi}}\int_{0}^{\infty}dx\;u(x)\int_{\mathbb{R}}d\xi\;\phi_{\varepsilon}(\xi)F(\xi)e^{-i\xi
x}\\\
=&\frac{1}{\sqrt{2\pi}}\int_{0}^{1}dx\;u(x)\int_{\mathbb{R}}d\xi\;\phi_{\varepsilon}(\xi)F(\xi)e^{-i\xi
x}-\frac{1}{\sqrt{2\pi}}\int_{1}^{\infty}dx\;\frac{u(x)}{x^{2}}\int_{\mathbb{R}}d\xi\;\left(\phi_{\varepsilon}(\xi)F(\xi)\right)^{\prime\prime}e^{-i\xi
x},\end{split}$

where for the last equality we integrated twice by parts in $\xi$. By a change
of coordinates and the dominated convergence theorem, since $F(0)=0$ and
$|F(\varepsilon\xi)\phi(\xi)|\leq|\phi(\xi)|$ we see for the first term as
$\varepsilon\to 0$

$\left|\frac{1}{\sqrt{2\pi}}\int_{0}^{1}dx\;u(x)\int_{\mathbb{R}}d\xi\;\phi_{\varepsilon}(\xi)F(\xi)e^{-i\xi
x}\right|\leq\frac{1}{\sqrt{2\pi}}\int_{0}^{1}dx\;u(x)\int_{\mathbb{R}}d\xi\;|F(\varepsilon\xi)\phi(\xi)|\to
0.$

Thus, we shall consider only the second term. We use the following well-known
estimate $\left|e^{-i\xi x}-1\right|\leq 2|\xi|^{\delta}|x|^{\delta}$ for
$0<\delta<1$ and $x\in\mathbb{R}$. Then using
$\int_{\mathbb{R}}\left(\phi_{\varepsilon}F\right)^{\prime\prime}=0$

$\begin{split}\frac{1}{\sqrt{2\pi}}\int_{1}^{\infty}dx\;\frac{u(x)}{x^{2}}\int_{\mathbb{R}}d\xi\;\left(\phi_{\varepsilon}(\xi)F(\xi)\right)^{\prime\prime}e^{-i\xi
x}=&\frac{1}{\sqrt{2\pi}}\int_{1}^{\infty}dx\;\frac{u(x)}{x^{2}}\int_{\mathbb{R}}d\xi\;\left(\phi_{\varepsilon}(\xi)F(\xi)\right)^{\prime\prime}\left(e^{-i\xi
x}-1\right),\end{split}$

and hence

$\begin{split}\left|\frac{1}{\sqrt{2\pi}}\int_{1}^{\infty}dx\;\frac{u(x)}{x^{2}}\int_{\mathbb{R}}d\xi\;\left(\phi_{\varepsilon}(\xi)F(\xi)\right)^{\prime\prime}e^{-i\xi
x}\right|\leq\frac{1}{\sqrt{2\pi}}\int_{1}^{\infty}dx\;\frac{u(x)}{x^{2-\delta}}\int_{\mathbb{R}}d\xi\;\left|\left(\phi_{\varepsilon}(\xi)F(\xi)\right)^{\prime\prime}\right|2|\xi|^{\delta}.\end{split}$

Now we notice that $\int_{1}^{\infty}dx\;\frac{u(x)}{x^{2-\delta}}<\infty$ and
also we see that $F(\xi)\simeq\frac{\xi^{2}}{3}$ as $x\to 0$, similarly as
$\xi\to 0$ also $F^{\prime}(\xi)\simeq\frac{2}{3}\xi$ and
$F^{\prime\prime}(\xi)\simeq\frac{2}{3}$. Hence, with a change of variables we
see that

$\begin{split}\int_{\mathbb{R}}d\xi\;&\left|\left(\phi_{\varepsilon}(\xi)F(\xi)\right)^{\prime\prime}\right||\xi|^{\delta}\\\
\leq&\int_{\mathbb{R}}d\xi\;\left[|\phi(\xi)||F^{\prime\prime}(\varepsilon\xi)|\varepsilon^{\delta}|\xi|^{\delta}+2|\phi^{\prime}(\xi)|\frac{|F^{\prime}(\varepsilon\xi)||\xi|^{\delta}}{\varepsilon^{1-\delta}}+|\phi^{\prime\prime}(\xi)|\frac{|F(\varepsilon\xi)||\xi|^{\delta}}{\varepsilon^{2-\delta}}\right].\end{split}$

With the consideration above about $F$ and since $\phi\in
C_{c}^{\infty}\left(\mathbb{R}\right)$ we see that there exists a constant
$C=2\Arrowvert\phi\Arrowvert_{C_{c}^{\infty}\left(\mathbb{R}\right)}\left(\max\limits_{\text{supp}\phi}|\xi|\right)^{2+\delta}<\infty$
such that

$|\phi(\xi)||F^{\prime\prime}(\varepsilon\xi)|\varepsilon^{\delta}|\xi|^{\delta}+|\phi^{\prime}(\xi)|\frac{|F^{\prime}(\varepsilon\xi)||\xi|^{\delta}}{\varepsilon^{1-\delta}}+|\phi^{\prime\prime}(\xi)|\frac{|F(\varepsilon\xi)||\xi|^{\delta}}{\varepsilon^{2-\delta}}\leq
C\varepsilon^{\delta}$

for any $\xi\in\text{supp}(\phi)$. Thus, again with the dominated convergence
theorem we conclude

$\left|\frac{1}{\sqrt{2\pi}}\int_{1}^{\infty}dx\;\frac{u(x)}{x^{2}}\int_{\mathbb{R}}d\xi\;\left(\phi_{\varepsilon}(\xi)F(\xi)\right)^{\prime\prime}e^{-i\xi
x}\right|\to 0,$

which implies the first claim, namely $\hat{W}(0)=0$.

As next step we prove that the limit $\lim\limits_{x\to\infty}u(x)$ exists.
First of all we know that in distributional sense $\hat{u}$ solves the
equation

$F\hat{u}\overset{\mathcal{S}^{\prime}}{=}\hat{W}.$ (3.15)

Given any distributional solution $\hat{u}$ to (3.15) also
$\hat{u}+\hat{u}_{h}$ is a solution, where $\hat{u}_{h}$ is the homogeneous
solution to $F\hat{u}_{h}\overset{\mathcal{S}^{\prime}}{=}0$. Let us consider
the tempered distribution given by $\hat{u}_{h}$ and let
$\varphi\in\mathcal{S}(\mathbb{R})$ be any testfunction with support away from
zero, i.e. $\text{supp}(\varphi)\subset\mathbb{R}\setminus\\{0\\}$. Since
$F(\xi)=0$ if and only $\xi=0$ and since it is bounded, the function
$\frac{\varphi}{F}\in\mathcal{S}\left(\mathbb{R}\right)$. Hence,
$\int_{\mathbb{R}}\hat{u}_{h}\varphi=0$. This implies (see [29]) that
$\hat{u}_{h}\overset{\mathcal{S}^{\prime}}{=}\sum\limits_{0\leq\alpha<m}c_{\alpha}(D^{\alpha}\delta)$,
for $c_{\alpha}$ constants and a suitable $m\in\mathbb{N}$. Since
$c_{\alpha}F(D^{\alpha}\delta)\not\equiv 0$ for any $\alpha\geq 2$ we conclude

$\hat{u}_{h}=c_{0}\delta+c_{1}\delta^{\prime}$

for suitable constants $c_{0},c_{1}$. Using the smoothness of $\hat{W}$ we can
write $\hat{W}(\xi)=\hat{W}^{\prime}(0)\xi+H(\xi)$ where
$\hat{W}^{\prime}(0)=\frac{m_{1}(W)}{\sqrt{2\pi}i}$ and $H\in
C^{\infty}\left(\mathbb{R}\right)$ with $H(0)=H^{\prime}(0)=0$. Let us
consider the behavior of $F$

$F(\xi)\simeq\begin{cases}\frac{\xi^{2}}{3}-\frac{\xi^{4}}{5}+\mathcal{O}(\xi^{6})&\xi\to.0,\\\
1-\frac{\pi}{2\xi}+\mathcal{O}\left(\frac{1}{\xi^{2}}\right)&\xi\to\infty\end{cases}$
(3.16)

Hence,

$f(\xi):=\hat{W}(\xi)-\frac{3m_{1}(W)}{\sqrt{2\pi}i}\frac{F(\xi)}{\xi}\in
L^{2}\left(\mathbb{R}\right)$ (3.17)

and it also satisfies

$f(\xi)\simeq
H^{\prime\prime}(0)\xi^{2}+\mathcal{O}(\xi^{3})\;\;\;\;\;\;\text{ as }\xi\to
0.$ (3.18)

By the boundedness of $F$ and given its behavior as in (3.16) we conclude that
the function $\hat{h}:=\frac{f}{F}\in L^{2}\left(\mathbb{R}\right)$, in
particular $\hat{h}$ is well-defined in zero. It is easy to see that $\hat{u}$
solves

$F(\xi)\hat{u}(\xi)\overset{\mathcal{S}^{\prime}}{=}\frac{3m_{1}(W)}{\sqrt{2\pi}i}\frac{F(\xi)}{\xi}+f(\xi).$
(3.19)

Therefore, since $\hat{h}\in L^{2}(\mathbb{R})$ we have that
$\hat{u}(\xi)=\frac{3m_{1}(W)}{\sqrt{2\pi}i}PV\left(\frac{1}{\xi}\right)+\hat{h}(\xi)$
is a solution to (3.19). We denote by $PV(\cdot)$ the principal value. Thus,
adding the homogeneous solution we conclude

$\hat{u}(\xi)\overset{\mathcal{S}^{\prime}}{=}c_{0}\delta+c_{1}\delta^{\prime}+\frac{3}{2i}m_{1}(W)\sqrt{\frac{2}{\pi}}PV\left(\frac{1}{\xi}\right)+\hat{h}(\xi),$

which yields

$u(x)\overset{\mathcal{S}^{\prime}}{=}\frac{c_{0}}{\sqrt{2\pi}}-\frac{c_{1}i}{\sqrt{\pi}}x+\frac{3}{2}m_{1}(W)\text{sgn}(x)+h(x),$

where $h\in L^{2}(\mathbb{R})$ is the inverse transform of $\hat{h}$. Since
$u$ is bounded and satisfies $u(x)=0$ for all $x<0$, we have in distributional
sense

$u(x)=\frac{3}{2}m_{1}(W)+\frac{3}{2}m_{1}(W)\text{sgn}(x)+h(x).$

Hence for $x>0$ also $u(x)=3m_{1}(W)+h(x)$ pointwise. Lemma 3.3 implies also
that $h$ is uniformly continuous in the positive real line. Hence,
$\lim\limits_{x\to\infty}h(x)=0$ and therefore the limit of $u$ as
$x\to\infty$ exists and is uniquely determined by $g_{\nu}(n)$ and $N$. This
is true since

$\lim\limits_{y\to\infty}u(y)=3m_{1}(W)=3\left(\int_{0}^{\infty}dx\;xG(x)-\int_{-\infty}^{0}dx\;x\int_{0}^{\infty}dy\;K(y-x)u(y)\right)\geq
0.$

Also the positivity of the limit is guaranteed when
$\left\\{n\in\mathbb{S}:n\cdot N<0\text{ and
}\int_{0}^{\infty}d\nu\overline{g}_{\nu}(n)\not\equiv 0\right\\}$ is not a
zero measure set. ∎

We will define
$\overline{u}_{\infty}(p):=\lim\limits_{y\to\infty}\overline{u}(y,p)$ for
$p\in\partial\Omega$. We can also show that $\overline{u}$ converges to
$\overline{u}_{\infty}$ with exponential rate.

######  Lemma 3.4.

Let $u$ be the unique bounded solution to the problem (3.10) and
$u_{\infty}=\lim\limits_{x\to\infty}u(x)$. Then there exists a constant $C>0$
such that

$|u-u_{\infty}|\leq Ce^{-\frac{|x|}{2}}.$

###### Proof.

We use the same notation as in Theorem 3.3. Hence, we know that

$\hat{u}(\xi)\overset{\mathcal{S}^{\prime}}{=}\frac{u_{\infty}}{2}\sqrt{2\pi}\delta+\frac{u_{\infty}}{2}\sqrt{\frac{2}{\pi}}PV\left(\frac{1}{\xi}\right)+\hat{h}(\xi),$
(3.20)

with
$F(\xi)\hat{h}(\xi)=\hat{W}(\xi)-\frac{3m_{1}(W)}{\sqrt{2\pi}i}\frac{F(\xi)}{\xi}$.
By the definition of $W$ we see

$\lim\limits_{x\nearrow 0}W(x)-\lim\limits_{x\searrow
0}W(x)=W(0^{+})-W(0^{-})=u(0).$ (3.21)

We recall that $W$ has exactly one discontinuity in $x=0$ and that
$W\raisebox{2.0pt}{$\chi$}_{\\{x<0\\}}\in
C^{\infty}\left(\mathbb{R}_{-}\right)$ and
$W\raisebox{2.0pt}{$\chi$}_{\\{x>0\\}}\in
C^{\infty}\left(\mathbb{R}_{+}\right)$. By the monotonicity of the two
functions $W\raisebox{2.0pt}{$\chi$}_{\\{x<0\\}}$ and
$W\raisebox{2.0pt}{$\chi$}_{\\{x>0\\}}$ and since $W\in
L^{\infty}\left(\mathbb{R}\right)$ we see that
$W^{\prime}\raisebox{2.0pt}{$\chi$}_{\\{x<0\\}}\in L^{1}(\mathbb{R}_{-})$ and
$W^{\prime}\raisebox{2.0pt}{$\chi$}_{\\{x>0\\}}\in L^{1}(\mathbb{R}_{+})$.
Moreoevr, we have the asymptotics
$\hat{W}(\xi)\simeq\frac{u(0)}{\sqrt{2\pi}i\xi}+\mathcal{O}\left(\frac{1}{\xi^{1+\delta}}\right)$
as $|\xi|\to\infty$ for $0<\delta<1$. Indeed, integrating by parts and using
that $\lim\limits_{|x|\to\infty}W(x)=0$ we compute

$\begin{split}\sqrt{2\pi}\hat{W}(\xi)=&\int_{-\infty}^{0}W(x)e^{-i\xi
x}\;dx+\int_{0}^{\infty}W(x)e^{-i\xi x}\;dx\\\
=&\frac{u(0)}{i\xi}+\frac{1}{i\xi}\left(\int_{-\infty}^{0}W^{\prime}(x)e^{-i\xi
x}\;dx+\int_{0}^{\infty}W^{\prime}(x)e^{-i\xi x}\;dx\right)\\\
=&\frac{u(0)}{i\xi}-\frac{1}{i\xi}\left(\int_{-\infty}^{-1}dx\int_{0}^{\infty}dy\;\frac{e^{-(y-x)}u(y)}{2(y-x)}e^{-i\xi
x}+\int_{-\infty}^{0}dx\int_{1}^{\infty}dy\;\frac{e^{-(y-x)}u(y)}{2(y-x)}e^{-i\xi
x}\right)\\\
&-\frac{1}{i\xi}\left(\int_{-1}^{0}dx\int_{0}^{1}dy\;\frac{e^{-(y-x)}u(y)}{2(y-x)}\frac{d}{dx}\frac{e^{-i\xi
x}-1}{-i\xi}\right)\\\
&+\frac{1}{i\xi}\left(\int_{1}^{\infty}G^{\prime}(x)e^{-i\xi
x}\;dx+\int_{0}^{1}G^{\prime}(x)\frac{d}{dx}\frac{e^{-i\xi
x}-1}{-i\xi}\right)\end{split}$ (3.22)

We conclude integrating by parts and applying the Riemann-Lebesgue Theorem in
the following way. First of all, the function
$\partial_{x}\frac{e^{-(y-x)}}{(y-x)}$ is integrable on
$(-\infty,1)\times\mathbb{R}_{+}\cup\mathbb{R}_{-}\times(1,\infty)$ and also
$G^{\prime\prime}(x)$ is integrable in $(1,\infty)$. Moreover, using
$\left|e^{-i\xi x}-1\right|\leq 2|\xi|^{\delta}|x|^{\delta}$ for $0<\delta<1$
we have

$\int_{-1}^{0}dx\int_{0}^{1}dy\;\frac{e^{-(y-x)}u(y)}{(y-x)}|x|^{\delta}\leq
C\int_{-1}^{0}dx\;\left(|x|^{\delta-1}\right)<\infty$

and

$\int_{0}^{1}dx\;G^{\prime\prime}(x)|x|^{\delta}\leq
C\int_{0}^{1}\frac{e^{-x}}{x^{1-\delta}}\;dx<\infty.$

For this last estimate we also used that
$\frac{d}{d\theta}e^{\frac{x}{\cos(\theta)}}=x\frac{e^{\frac{x}{cos(\theta)}}}{\cos^{2}(\theta)}\sin(\theta)$,
which implies $\left|G^{\prime\prime}(x)\right|\leq 2\pi\Arrowvert
g\Arrowvert_{\infty}\frac{e^{-x}}{x}$. Thus, by the definition of $\hat{h}$
and using (3.16) we have

$\hat{h}(\xi)\simeq\begin{cases}\mathcal{O}(1)&|\xi|\to 0,\\\
\frac{u(0)}{\sqrt{2\pi}}\frac{1}{i\xi}-\frac{u_{\infty}}{\sqrt{2\pi}}\frac{1}{i\xi}+\mathcal{O}\left(\frac{1}{\xi^{1+\delta}}\right)&|\xi|\to\infty.\end{cases}$
(3.23)

By the definition of $\hat{u}$ in (3.20) we see

$\hat{v}(\xi):=\hat{u}(\xi)-\frac{u_{\infty}}{2}\sqrt{2\pi}\delta-
PV\left(\frac{1}{i\xi}\right)\left(\frac{u_{\infty}}{2}\sqrt{\frac{2}{\pi}}\frac{1}{1+\xi^{2}}+\frac{u(0)}{\sqrt{2\pi}}\frac{\xi^{2}}{1+\xi^{2}}\right)\in
L^{2}(\mathbb{R}).$ (3.24)

We claim that

  1. (i)

$\hat{v}$ is analytic in the strip
$S=\\{z\in\mathbb{C}:|\Im(z)|<\frac{3}{4}\\}$;

  2. (ii)

$|\hat{v}(\xi)|\leq\frac{C}{|1+\xi^{1+\delta}|}$;

  3. (iii)

$v(x)=u(x)-u_{\infty}+\frac{e^{-|x|}}{2}\left(u_{\infty}-u(0)\right)$ for
$x>0$ and $v(x)=\mathcal{F}^{-1}(\hat{v})(x)$.

A contour integral implies then the lemma. Indeed for $x>0$ we can compute

$\begin{split}\sqrt{2\pi}|v(x)|=&\lim\limits_{R\to\infty}\left|\int_{-R}^{R}\hat{v}(\xi)e^{i\xi
x}\;d\xi\right|\\\
\leq&\lim\limits_{R\to\infty}\left|i\int_{0}^{\frac{1}{2}}\hat{v}(R+it)e^{iRx}e^{-tx}\;dt\right|+\lim\limits_{R\to\infty}\left|i\int_{0}^{\frac{1}{2}}\hat{v}(-R+it)e^{-iRx}e^{-tx}\;dt\right|\\\
&+\lim\limits_{R\to\infty}\left|\int_{-R}^{R}\hat{v}\left(t+\frac{1}{2}i\right)e^{itx}e^{-\frac{x}{2}}\;dt\right|\\\
\leq&e^{-\frac{x}{2}}\lim\limits_{R\to\infty}\int_{-R}^{R}\frac{C}{\frac{1}{2}+t^{1+\delta}}\;dt=\overline{C}e^{-\frac{x}{2}},\end{split}$
(3.25)

where for the first inequality we used the triangle inequality and the
analycity of $\hat{v}$ by (i), the second inequality is due to dominate
convergence and the claim (ii), finally the last integral is finite. Equation
(3.25) and claim (iii) imply $|u(x)-u_{\infty}|\leq Ce^{-\frac{x}{2}}$ for
$x>0$.

We prove now the claims. To prove claim (i) it is enough to show that
$\hat{h}$ is analytic in $S$. Then, (3.20) and (3.24) implies (i). First of
all we recall that $W$ has an exponential decay like $|W(x)|\leq Ce^{-|x|}$,
hence $|W(x)|e^{\frac{3}{4}|x|}\in L^{1}\left(\mathbb{R}\right)$ and therefore
Paley-Wiener Theorem implies that $\hat{W}$ is analytic in $S$. Since
$\arctan(z)=\frac{1}{2i}\ln(\frac{1+iz}{1-iz})$ is analytic in
$\\{z\in\mathbb{C}:|\Im(z)|<1\\}$ and since $F(z)=\frac{z-\arctan(z)}{z}$ has
exactly one zero in $z=0$, which is of degree $2$, the definition of
$\hat{h}=\frac{f}{F}$ together with (3.17) implies that $\hat{h}$ is analytic
in $S$ since (3.18) implies that $0$ is a removable singularity.

For claim (ii) we just put together equations (3.20), (3.23) and (3.24). We
notice also that the constant $C>0$ of claim (ii) depends only on $\hat{W}$.

Claim (iii) is more involved. We have to consider again two different contour
integrals in order to compute the inverse Fourier transform of $\hat{v}$. We
start with considering the function
$PV\left(f(\xi)\right)=PV\left(\frac{1}{i\xi(1+\xi^{2})}\right)$. Let first of
all $x>0$ and let $\gamma^{+}_{1}$ the path around $i$ given as in the
following picture.

Figure 2: sketch of $\gamma^{+}_{1}$.

Hence, we compute

$\begin{split}\mathcal{F}^{-1}&\left(PV\left(\frac{1}{i\xi}\frac{1}{1+\xi^{2}}\right)\right)(x)=\frac{1}{\sqrt{2\pi}}\lim\limits_{R\to\infty}\left(\int_{-R}^{-\frac{1}{R}}f(\xi)e^{i\xi
x}\;d\xi+\int_{\frac{1}{R}}^{R}f(\xi)e^{i\xi x}\;d\xi\right)\\\
=&\frac{1}{\sqrt{2\pi}}\lim\limits_{R\to\infty}\left(\int_{\gamma^{+}_{1}}f(\xi)e^{i\xi
x}d\xi\right)+\frac{1}{\sqrt{2\pi}}\lim\limits_{R\to\infty}\left(\int_{0}^{\pi}f\left(\frac{e^{i\theta}}{R}\right)\frac{ie^{i\theta}}{R}e^{-\frac{\sin(\theta)x}{R}}e^{\frac{i\cos(\theta)x}{R}}\;d\theta\right)\\\
&-\frac{1}{\sqrt{2\pi}}\lim\limits_{R\to\infty}\left(\int_{0}^{\pi}f\left(Re^{i\theta}\right)Rie^{i\theta}e^{-R\sin(\theta)x}e^{iR\cos(\theta)x}\;d\theta\right)\\\
=&\sqrt{\frac{\pi}{2}}\left(1-e^{-x}\right).\end{split}$

For the computation of these integrals we used the Cauchy’s residue theorem
and $\text{Res}_{i}f(\xi)e^{i\xi x}=\frac{ie^{-x}}{2}$, the second integral
converges to $\pi$ as $R\to\infty$ and the third converges to zero, both limit
are due to the Lebesgue dominated convergence theorem. Denoting by
$\gamma^{-}_{1}$ the mirrored path to $\gamma^{+}_{1}$ with respect to the
real axis and arguing similarly we also get that for $x<0$ the inverse Fourier
transformation is
$\mathcal{F}^{-1}\left(PV\left(\frac{1}{i\xi}\frac{1}{1+\xi^{2}}\right)\right)(x)=-\sqrt{\frac{\pi}{2}}\left(1-e^{-|x|}\right)$.
Hence,

$\mathcal{F}^{-1}\left(PV\left(\frac{1}{i\xi}\frac{1}{1+\xi^{2}}\right)\right)(x)=\text{sgn}(x)\sqrt{\frac{\pi}{2}}\left(1-e^{-|x|}\right).$
(3.26)

For the function $g(x)(\xi)=\frac{\xi}{i(1+\xi^{2})}$ we consider again first
of all $x>0$ and the path $\gamma^{+}_{2}$ around $i$ given by

Figure 3: sketch of $\gamma^{+}_{2}$.

Hence, the Cauchy’s residue theorem and the dominated convergence imply

$\begin{split}\mathcal{F}^{-1}&\left(\frac{\xi}{i(1+\xi^{2})}\right)(x)=\frac{1}{\sqrt{2\pi}}\lim\limits_{R\to\infty}\int_{-R}^{R}g(\xi)e^{i\xi
x}\;d\xi\\\
=&\frac{1}{\sqrt{2\pi}}\lim\limits_{R\to\infty}\left(\int_{\gamma^{+}_{2}}g(\xi)e^{i\xi
x}d\xi\right)-\frac{1}{\sqrt{2\pi}}\lim\limits_{R\to\infty}\left(\int_{0}^{\pi}g\left(Re^{i\theta}\right)Rie^{i\theta}e^{-R\sin(\theta)x}e^{iR\cos(\theta)x}\;d\theta\right)\\\
=&\sqrt{\frac{\pi}{2}}e^{-x},\end{split}$

where we also used that $\text{Res}_{i}g(\xi)e^{-\xi x}=\frac{e^{-x}}{2i}$.
Denoting similarly as before by $\gamma^{-}_{2}$ the mirrored path to
$\gamma^{+}_{2}$ with respect to the real axis we obtain
$\mathcal{F}^{-1}(g)(x)=-\sqrt{\frac{\pi}{2}}e^{-|x|}$for $x<0$ and thus

$\mathcal{F}^{-1}\left(\frac{\xi}{i(1+\xi^{2})}\right)(x)=\text{sgn}(x)\sqrt{\frac{\pi}{2}}e^{-|x|}.$
(3.27)

Hence, the definition of $u$ and equations (3.26), (3.27) imply claim (iii)
for $x>0$

$v(x)=u(x)-u_{\infty}+\frac{e^{-|x|}}{2}\left(u_{\infty}-u(0)\right).$

∎

There are still two important properties of $\overline{u}(y,p)$ we will need
for the rest of the paper and which are explained in the next two Lemmas.
First of all $\overline{u}(y,p)$ is uniformly bounded in both variables.

######  Lemma 3.5.

Let $\overline{u}(y,p)$ be the non-negative bounded solution to the problem
(2.8) for $g_{\nu}(n)$ satisfying the assumption as in Theorem 3.2. Then there
exists a constant $C$ such that

$\sup\limits_{y\in\mathbb{R},\;p\in\partial\Omega}\overline{u}(y,p)\leq
C<\infty.$

###### Proof.

By definition $\overline{u}$ satisfies $\mathcal{L}(\overline{u})(y)=G_{p}(y)$
for $y>0$ and $\overline{u}(y,p)=0$ for $y<0$. Moreover, recalling the norm as
in (1.10) the source can be estimated by

$0\leq G_{p}(y)\leq\Arrowvert g\Arrowvert_{1}e^{-y},$

since $|n\cdot N_{p}|\leq 1$.

Theorem 3.2 assures us the existence of a unique bounded continuous (in the
positive line) solution $v$ of $\mathcal{L}(v)(y)=\Arrowvert
g\Arrowvert_{1}e^{-y}$ for $y>0$ and $v(y)=0$ for $y<0$. Hence, we can apply
the maximum principle of Theorem 3.1 as we did in Lemma 3.2 to the function
$v-\overline{u}(\cdot,p)\in C\left([0,\infty]\right)$ and we conclude

$0\leq\overline{u}(y,p)\leq v(y)\leq\Arrowvert v\Arrowvert_{\infty}:=C<\infty$

for all $y\in\mathbb{R}$ and $p\in\partial\Omega$. ∎

Also, the rate of convergence of $\overline{u}(y,p)$ to
$\overline{u}_{\infty}(p)$ can be bounded independently of
$p\in\partial\Omega$.

######  Corollary 3.2.

There exists a constant $C>0$ independent of $p\in\partial\Omega$ such that

$|\overline{u}(y,p)-\overline{u}_{\infty}(p)|\leq Ce^{-\frac{y}{2}}$

.

###### Proof.

This is a consequence of Lemma 3.4 and Lemma 3.5. From Lemma 3.5 we know that
there exists a constant $C>0$ independent of $p\in\partial\Omega$ such that

$\left|W(x)\right|\leq
C\left(e^{-|x|}+|x|K(x)\raisebox{2.0pt}{$\chi$}_{\\{x<0\\}}\right)\in
L^{1}(\mathbb{R})\cap L^{2}(\mathbb{R})\cap L^{\infty}(\mathbb{R}),$

where $W$ is the function defined in (3.13). Since
$|x|K(x)\leq\frac{e^{-|x|}}{2}$ all moments of $W$ are finite and for any
$n\in\mathbb{N}$ there exists a constant $C_{n}>0$ independent of
$p\in\partial\Omega$ such that

$\left|m_{n}\left(W\right)\right|\leq C_{n}<\infty.$

Hence, $\hat{W}\in C_{0}(\mathbb{R})\cap C^{\infty}(\mathbb{R})\cap
L^{2}(\mathbb{R})$ and also all derivatives are uniformly bounded in
$p\in\partial\Omega$ since
$\left|\hat{W}^{(n)}(\xi)\right|\leq\frac{C_{n}}{\sqrt{2\pi}}$. Thus, the
function $\hat{h}$ in (3.23) defined using (3.18) can be bounded independently
of $p\in\partial\Omega$.

Moreover, we notice that in (3.22) as $|\xi|\to\infty$ we can bound
$\left|\hat{W}(\xi)-\frac{u(0)}{\sqrt{(}2\pi)i\xi}\right|$ by
$\frac{C}{\left|\xi^{1+\delta}\right|}$ with a constant $C>0$ independent of
$p\in\partial\Omega$. Indeed, as we have seen in Lemma 3.4 we have
$\left|G^{\prime\prime}(x)\right|\leq 2\pi\Arrowvert
g\Arrowvert_{\infty}\frac{e^{-x}}{x}$ and by Lemma 3.5 we have also
$\left|\overline{u}(y,p)\right|\leq C$.

Hence, we conclude as we did in Lemma 3.4 that there exists a constant $C>0$
independent of $p\in\partial\Omega$ such that
$\left|\hat{v}(\xi)\right|\leq\frac{C}{|1+\xi^{1+\delta}|}$, where $\hat{v}$
was defined in (3.24).

Arguing now exactly as in Lemma 3.4 using also Lemma 3.5 we conclude that
there exists a constant $C>0$ independent of $p\in\partial\Omega$ such that
$|\overline{u}(y,p)-\overline{u}_{\infty}(p)|\leq Ce^{-\frac{y}{2}}$. ∎

Next, using again the maximum principle we can also show that
$\overline{u}(y,p)$ is Lipschitz continuous with respect to
$p\in\partial\Omega$ uniformly in $y$.

######  Lemma 3.6.

Let $g_{\nu}(n)$ be as in Theorem 3.2 and let $\overline{u}$ be the unique
bounded solution to (2.8). Then $\overline{u}$ is uniformly continuous with
respect the variable $p\in\partial\Omega$ uniformly in $y$. More precisely, it
is Lipschitz continuous, i.e. there exists a constant $C>0$ such that for
every $p,q\in\partial\Omega$

$\sup_{y\geq 0}\left|\overline{u}(y,p)-\overline{u}(y,q)\right|\leq
C|p-q|:=\omega_{1}(|p-q|).$

###### Proof.

The proof is based on the maximum principle. We start taking
$0<\tilde{\delta}<1$ sufficiently small and we consider $p,q\in\partial\Omega$
with $|p-q|<\tilde{\delta}$. We denote by $S_{p}(q)$ the plane defined by the
vector $\overset{\rightharpoonup}{pq}$ and the unit vector $N_{p}$. Given that
$\partial\Omega$ is a $C^{3}$-surface we can define $\rho_{p}$ to be the
radius of curvature of the curve $C_{p}(q):=S_{p}(q)\cap\partial\Omega$ at
$p$. Since by assumption the curvature of $\partial\Omega$ is bounded from
below by a positive constant, for $\tilde{\delta}$ small enough we can
estimate

$\frac{1}{2}\rho_{p}\theta_{pq}\leq|p-q|\leq 2\rho_{p}\theta_{pq},$ (3.28)

where $\theta_{pq}$ is the angle between $N_{p}$ and $N_{q}$. This is true,
because for $\tilde{\delta}$ sufficiently small the angle $\theta_{pq}$ is not
zero and it is approximately the central angle between the rays connecting $p$
and $q$ with the center of the circle with radius $\rho_{p}$ tangent to $p$.
We denote by $R$ the minimal radius of curvature of $\partial\Omega$, hence
$\rho_{p}\geq R$. Now we consider the operator $\mathcal{L}$ acting on the
difference $\overline{u}(y,p)-\overline{u}(y,q)$. We can estimate its absolute
value by the sum of the following six terms

$\begin{split}\left|\mathcal{L}\left(\right.\right.&\left.\left.\overline{u}(y,p)-\overline{u}(y,q)\right)\right|\leq\int_{A_{1}}\int_{0}^{\infty}g_{\nu}(n)e^{-\frac{y}{|n\cdot
N_{p}|}}\;d\nu\;dn+\int_{A_{2}}\int_{0}^{\infty}g_{\nu}(n)e^{-\frac{y}{|n\cdot
N_{q}|}}\;d\nu\;dn\\\
&+\int_{A_{3}}\int_{0}^{\infty}g_{\nu}(n)\left|e^{-\frac{y}{|n\cdot
N_{p}|}}-e^{-\frac{y}{|n\cdot
N_{q}|}}\right|\;d\nu\;dn+\int_{A_{4}}\int_{0}^{\infty}g_{\nu}(n)\left|e^{-\frac{y}{|n\cdot
N_{p}|}}-e^{-\frac{y}{|n\cdot N_{q}|}}\right|\;d\nu\;dn\\\
&+\int_{A_{5}}\int_{0}^{\infty}g_{\nu}(n)\left|e^{-\frac{y}{|n\cdot
N_{p}|}}-e^{-\frac{y}{|n\cdot
N_{q}|}}\right|\;d\nu\;dn+\int_{A_{6}}\int_{0}^{\infty}g_{\nu}(n)\left|e^{-\frac{y}{|n\cdot
N_{p}|}}-e^{-\frac{y}{|n\cdot N_{q}|}}\right|\;d\nu\;dn,\end{split}$ (3.29)

where we denote by $A_{i}$ the following sets

$\begin{split}A_{1}&:=\left\\{n\in\mathbb{S}^{2}:n\cdot N_{p}<0,\;n\cdot
N_{q}\geq 0\right\\},\;\;\;\;A_{2}:=\left\\{n\in\mathbb{S}^{2}:n\cdot
N_{p}\geq 0,\;n\cdot N_{q}<0\right\\},\\\
A_{3}&:=\left\\{n\in\mathbb{S}^{2}:n\cdot N_{p}<0,\;n\cdot N_{q}<0,\;|n\cdot
N_{p}|\geq|n\cdot N_{q}|,|;|n\cdot N_{p}|>\frac{4}{R}|p-q|\right\\},\\\
A_{4}&:=\left\\{n\in\mathbb{S}^{2}:n\cdot N_{p}<0,\;n\cdot N_{q}<0,\;|n\cdot
N_{p}|\geq|n\cdot N_{q}|,|;|n\cdot N_{p}|\leq\frac{4}{R}|p-q|\right\\},\\\
A_{5}&:=\left\\{n\in\mathbb{S}^{2}:n\cdot N_{p}<0,\;n\cdot N_{q}<0,\;|n\cdot
N_{q}|\geq|n\cdot N_{p}|,|;|n\cdot N_{q}|>\frac{4}{R}|p-q|\right\\}\text{ and
}\\\ A_{6}&:=\left\\{n\in\mathbb{S}^{2}:n\cdot N_{p}<0,\;n\cdot
N_{q}<0,\;|n\cdot N_{q}|\geq|n\cdot N_{p}|,|;|n\cdot
N_{q}|\leq\frac{4}{R}|p-q|\right\\}.\\\ \end{split}$

By symmetry, we need to estimate only the first, the third and the fourth
terms. We start with the first line of equation (3.29). The set $A_{1}$ is
contained by the set given by all the $n$ such that their angle with $N_{p}$
is in the interval $(\frac{\pi}{2},\frac{\pi}{2}+\theta_{pq})$. Using the fact
that $\frac{y}{|n\cdot N_{p}|}>y$, we estimate the exponential by $e^{-y}$ and
hence we see

$\int_{A_{1}}\int_{0}^{\infty}g_{\nu}(n)e^{-\frac{y}{|n\cdot
N_{p}|}}\;d\nu\;dn\leq\Arrowvert
g\Arrowvert_{\infty}2\pi\theta_{pq}e^{-y}\leq\frac{4\pi}{R}\Arrowvert
g\Arrowvert_{\infty}e^{-y}.$ (3.30)

The second term in (3.29) is estimated similarly. For the third term of
equation (3.29) we estimate the difference of the exponential as follows,
assuming $|n\cdot N_{p}|\geq|n\cdot N_{q}|$

$\left|e^{-\frac{y}{|n\cdot N_{p}|}}-e^{-\frac{y}{|n\cdot N_{q}|}}\right|\leq
e^{-\frac{y}{|n\cdot N_{p}|}}y\left|\frac{1}{|n\cdot N_{q}|}-\frac{1}{|n\cdot
N_{p}|}\right|\leq e^{-\frac{y}{|n\cdot N_{p}|}}y\left|\frac{|n\cdot
N_{p}|-|n\cdot N_{q}|}{|n\cdot N_{q}||n\cdot N_{p}|}\right|,$

where we used for $x>0$ the inequality $1-e^{-x}\leq x$. By definition
$|n\cdot(N_{p}-N_{q})|\leq\theta_{pq}\leq\frac{2}{R}|p-q|$ which implies

$0\leq|n\cdot N_{p}|-|n\cdot
N_{q}|=|n\cdot(N_{q}-N_{p})|\leq\frac{2}{R}|p-q|.$

Since $|n\cdot N_{p}|>\frac{4}{R}|p-q|$ we see also that

$|n\cdot N_{q}|\geq|n\cdot N_{p}|-\frac{2}{R}|p-q|\geq\frac{|n\cdot
N_{p}|}{2}.$

Hence,

$\left|e^{-\frac{y}{|n\cdot N_{p}|}}-e^{-\frac{y}{|n\cdot N_{q}|}}\right|\leq
e^{-\frac{y}{|n\cdot N_{p}|}}y\frac{4|p-q|}{R|n\cdot N_{p}|^{2}}.$

Putting together these inequalities we compute

$\begin{split}\int_{A_{3}}\int_{0}^{\infty}g_{\nu}(n)&\left|e^{-\frac{y}{|n\cdot
N_{p}|}}-e^{-\frac{y}{|n\cdot
N_{q}|}}\right|\;d\nu\;dn\leq\frac{4|p-q|}{R}\Arrowvert
g\Arrowvert_{\infty}\int_{A_{3}}dn\;e^{-\frac{y}{|n\cdot
N_{p}|}}\frac{y}{|n\cdot N_{p}|^{2}}\\\ \leq&\frac{4|p-q|}{R}\Arrowvert
g\Arrowvert_{\infty}4\pi\int_{0}^{\frac{\pi}{2}}e^{-\frac{y}{\cos(\theta)}}\frac{y\sin(\theta)}{\cos^{2}(\theta)}\;d\theta=\frac{16\pi|p-q|}{R}\Arrowvert
g\Arrowvert_{\infty}e^{-y},\end{split}$ (3.31)

where we estimated the last integral in $A_{3}$ using polar coordinates in
$\mathbb{S}^{2}$ using as reference $N_{p}$. It remains to estimate the
integral on $A_{4}$. For this term we use the inclusion

$\begin{split}A_{4}\subset&\left\\{n\in\mathbb{S}^{2}:n\cdot N_{p}<0,\;|n\cdot
N_{p}|\leq\frac{4}{R}|p-q|\right\\}\\\
\subset&\left\\{(\varphi,\theta)\in[0,2\pi]\times[0,\pi]:\theta\in\left(-\frac{\pi}{2},-\frac{\pi}{2}+C(R)|p-q|\right)\cup\left(\frac{\pi}{2}-C(R)|p-q|,\frac{\pi}{2}\right)\right\\},\end{split}$

where the last inclusion is due to the smallness of $\frac{4}{R}|p-q|<1$ and
the expansion of the arc-cosine. Moreover, $C(R)$ is a constant depending only
on $R$. Hence, as we estimated in (3.30) we have

$\int_{A_{4}}\int_{0}^{\infty}g_{\nu}(n)\left|e^{-\frac{y}{|n\cdot
N_{p}|}}-e^{-\frac{y}{|n\cdot N_{q}|}}\right|\;d\nu\;dn\leq C(R)4\pi\Arrowvert
g\Arrowvert_{\infty}|p-q|.$ (3.32)

Now, with equations (3.30),(3.31) and (3.32) we estimate the operator by

$\left|\mathcal{L}\left(\overline{u}(y,p)-\overline{u}(y,q)\right)\right|\leq
C(R)\Arrowvert g\Arrowvert_{\infty}|p-q|e^{-y},$

where $C(R)>0$ is a constant depending only on the minimal radius of curvature
$R$. Theorem 3.2 and the maximum principle imply the existence of a unique
non-negative bounded continuous function $V$ solution to the equation
$\mathcal{L}(V)(y)=e^{-y}$ for $y\geq 0$. Hence, we apply the maximum
principle of Theorem 3.1 as in Lemma 3.2 to the continuous functions
$C(R)\Arrowvert
g\Arrowvert_{\infty}|p-q|V-\left(\overline{u}(y,p)-\overline{u}(y,q)\right)$
and $C(R)\Arrowvert
g\Arrowvert_{\infty}|p-q|V-\left(\overline{u}(y,q)-\overline{u}(y,p)\right)$.
We conclude the uniformly continuity of $\overline{u}(y,p)$ in $p$ uniformly
in $y$

$\left|\overline{u}(y,p)-\overline{u}(y,q)\right|\leq C(R)\Arrowvert
g\Arrowvert_{\infty}|p-q|.$

The modulus of continuity $\omega_{1}$ is hence defined by
$\omega_{1}(r)=C(R)\Arrowvert g\Arrowvert_{\infty}r$. ∎

######  Corollary 3.3.

The limit $\overline{u}_{\infty}$ is Lipschitz continuous in
$p\in\partial\Omega$.

###### Proof.

This is a direct consequence of the previous Lemma 3.6. The modulus of
continuity of $\overline{u}_{\infty}$ is still the same $\omega_{1}$ of
$\overline{u}(y,p)$. ∎

Finally, we summarize all properties of $\overline{u}$ in the following
proposition.

######  Proposition 3.1.

Let $g_{\nu}(n)$ be as in Theorem 3.2 and $\Omega$ as in the assumption. For
every $p\in\partial\Omega$ there exists a unique non-negative bounded solution
$\overline{u}(y,p)$ to (2.8). For every $p\in\partial\Omega$ the function
$\overline{u}(\cdot,p)$ is uniformly continuous in $[0,\infty)$ and has a non-
negative limit
$\overline{u}_{\infty}(p)=\lim\limits_{y\to\infty}\overline{u}(y,p)$, which is
strictly positive if $\left\\{n\in\mathbb{S}:n\cdot N_{p}<0\text{ and
}\int_{0}^{\infty}d\nu g_{\nu}(n)\not\equiv 0\right\\}$ is not a zero measure
set. Moreover, $\overline{u}(y,p)$ is uniformly bounded in both variables and
it is Lipschitz continuous with respect to $p\in\partial\Omega$ uniformly on
$y\in\mathbb{R}_{+}$. Finally, $\overline{u}_{\infty}$ is Lipschitz continuous
and there exists a constant $C>0$ independent of $p\in\partial\Omega$ such
that $|\overline{u}(y,p)-\overline{u}_{\infty}(p)|\leq Ce^{-\frac{|y|}{2}}$.

##  4 Rigorous proof of the diffusion equilibrium approximation for constant
absorption coefficient

This section of the paper deals with the rigorous proof of the diffusion
equilibrium approximation for the constant absorption coefficient case. We
will show that the Stefan-Boltzmann law $u^{\varepsilon}(x)=4\pi\sigma
T_{\varepsilon}^{4}(x)$ for the temperature $T_{\varepsilon}$ associated to
the boundary value problem (1.5) converges pointwise as $\varepsilon\to 0$ to
$v$, the solution to the Dirichlet problem

$\begin{cases}-\Delta v=0&\text{ in }\Omega,\\\ v=\overline{u}_{\infty}&\text{
on }\partial\Omega,\end{cases}$ (4.1)

where $\overline{u}_{\infty}$ is defined as in Proposition 3.1.

###  4.1 Derivation of the equation for $u^{\varepsilon}$

Let us call $I^{\varepsilon}_{\nu}$ the solution to the initial boundary value
problem (1.5). We start with the derivation of the integral equation satisfied
by $u^{\varepsilon}=4\pi\sigma T_{\varepsilon}^{4}$. To this end we solve by
characteristics the equation

$n\cdot\nabla_{x}I_{\nu}\left(x,n\right)=\frac{1}{\varepsilon}\left(B_{\nu}\left(T\left(x\right)\right)-I_{\nu}\left(x,n\right)\right)$

Let $x\in\Omega$ and $n\in\mathbb{S}^{2}$. The convexity of $\Omega$ implies
the existence of a unique $x_{\Omega}(x,n)\in\partial\Omega$ connecting $x$ in
direction $-n$ with the boundary $\partial\Omega$. Hence,
$\frac{x-x_{\Omega}(x,n)}{\left|x-x_{\Omega}(x,n)\right|}=n$ and we define
$s(x,n)=\left|x-x_{\Omega}(x,n)\right|$. Then $x=x_{\Omega}(x,n)+s(x,n)n$.
Integrating along the characteristics equation (1.5) we get

$\begin{split}I^{\varepsilon}_{\nu}(x,n)=g_{\nu}(n)e^{-\frac{\left|x-x_{\Omega}(x,n)\right|}{\varepsilon}}+\frac{1}{\varepsilon}\int_{0}^{s(x,n)}e^{-\frac{t}{\varepsilon}}B_{\nu}\left(T\left({x-tn}\right)\right)\;dt.\end{split}$

Using the heat equation, i.e. $\nabla_{x}\cdot\mathcal{F}=0$ (see (1.5)), we
calculate

$0=\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;n\cdot\nabla_{x}I^{\varepsilon}_{\nu}(x,n)=\frac{1}{\varepsilon}\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;B_{\nu}(T_{\varepsilon})(x)-I^{\varepsilon}_{\nu}(x,n).$

We define $u^{\varepsilon}(x)=4\pi\sigma
T_{\varepsilon}^{4}(x)=\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dnB_{\nu}\left(T_{\varepsilon}(x)\right)$
according to (1.3). Hence also
$u^{\varepsilon}(x)=\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;I_{\nu}^{\varepsilon}(x,n)$.
We integrate now the expression we got for the intensity and we conclude with
the equation satisfied by $u^{\varepsilon}$ as follows

$\begin{split}u^{\varepsilon}(x)&=\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;g_{\nu}(n)e^{-\frac{\left|x-x_{\Omega}(x,n)\right|}{\varepsilon}}+\frac{1}{4\pi\varepsilon}\int_{\mathbb{S}^{2}}dn\int_{0}^{s(x,n)}e^{-\frac{t}{\varepsilon}}u^{\varepsilon}(x-tn)\;dt\\\
&=\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;g_{\nu}(n)e^{-\frac{\left|x-x_{\Omega}(x,n)\right|}{\varepsilon}}+\frac{1}{4\pi\varepsilon}\int_{\Omega}\frac{e^{-\frac{\left|x-\eta\right|}{\varepsilon}}}{\left|x-\eta\right|^{2}}u^{\varepsilon}(\eta)\;d\eta,\end{split}$

where the last equality is due to the change of variables
$\mathbb{S}^{2}\times(0,\infty)\to\Omega$ with $(n,t)\mapsto x-tn=\eta$. Hence
the sequence $u^{\varepsilon}$ of exact solutions solves

$u^{\varepsilon}(x)-\int_{\Omega}\frac{e^{-\frac{\left|x-\eta\right|}{\varepsilon}}}{4\pi\varepsilon\left|x-\eta\right|^{2}}u^{\varepsilon}(\eta)\;d\eta=\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;g_{\nu}(n)e^{-\frac{\left|x-x_{\Omega}(x,n)\right|}{\varepsilon}}.$
(4.2)

We define the kernel
$K_{\varepsilon}(x):=\frac{e^{-\frac{\left|{x}\right|}{\varepsilon}}}{4\pi\varepsilon\left|{x}\right|^{2}}$
and we notice that its integral in $\mathbb{R}^{3}$ is $1$.

######  Remark.

There exists a unique solution $u^{\varepsilon}$ continuous and bounded. We
adapt the proof in [21]. The existence and uniqueness of a solution
$u^{\varepsilon}\in L^{\infty}\left(\Omega\right)$ can be shown with the
Banach fixed-point Theorem. We define for every given $g$ and $\varepsilon>0$
the self map $A_{g}^{\varepsilon}:L^{\infty}\left(\Omega\right)\to
L^{\infty}\left(\Omega\right)$ by

$A_{g}^{\varepsilon}(u)(x)=\int_{\Omega}K_{\varepsilon}\left(\eta-x\right)u(\eta)\;d\eta+\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;g_{\nu}(n)e^{-\frac{\left|x-x_{\Omega}(x,n)\right|}{\varepsilon}}.$

Then since
$\int_{\Omega}K_{\varepsilon}(\eta-x)\;d\eta<\int_{\mathbb{R}^{3}}K_{\varepsilon}(\eta-x)\;d\eta=1$
we conclude that $A_{g}^{\varepsilon}$ is a contraction, hence there is a
unique fixed-point, which is the desired unique solution. Moreover,
$G^{\varepsilon}_{x_{\Omega}}(x):=\int_{0}^{\infty}d\nu\int_{\mathbb{S}^{2}}dn\;g_{\nu}(n)e^{-\frac{\left|x-x_{\Omega}(x,n)\right|}{\varepsilon}}$
is continuous and since $u^{\varepsilon}\in L^{\infty}\left(\Omega\right)$ and
$K_{\varepsilon}(x-\cdot)\in L^{{1}}\left(\mathbb{R}^{3}\right)$ we conclude
that the convolution
$\int_{\Omega}K_{\varepsilon}(\eta-x)u^{\varepsilon}(\eta)\;d\eta$ is
continuous and bounded. Hence, $u^{\varepsilon}$ is continuous and bounded. We
can also extend continuously $u^{\varepsilon}$ to the boundary
$\partial\Omega$ defining $|x-x_{\Omega}(x,n)|=0$ for $x\in\partial\Omega$ and
$n\cdot N_{x}\leq 0$. Then using the generalized dominated convergence theorem
we see that both integral terms in (4.2) are continuous up to the boundary.
Hence, $u^{\varepsilon}\in C\left(\overline{\Omega}\right)$. Moreover,
$u^{\varepsilon}$ is non-negative. This is because of the maximum principle as
stated in the following theorem.

######  Theorem 4.1 (Maximum Principle).

Let $v$ be bounded and continuous, $v\in C\left(\overline{\Omega}\right)$. Let
$\mathcal{L}_{\Omega}^{\varepsilon}(v)(x)=v(x)-\int_{\Omega}K_{\varepsilon}(\eta\leavevmode\nobreak\
-\leavevmode\nobreak\ x)v(\eta)\;d\eta$. Assume $v$ satisfies one of the
following properties:

  1. (i)

$\mathcal{L}_{\Omega}^{\varepsilon}(v)(x)\geq 0$ if $x\in\Omega$;

  2. (ii)

$\mathcal{L}_{\Omega}^{\varepsilon}(v)(x)\geq 0$ if $x\in O\subset\Omega$ open
and $v(x)\geq 0$ if $x\in\Omega\setminus O$.

Then, $v\geq 0$.

###### Proof.

Let $y\in\overline{\Omega}$ such that $v(y)=\min_{x\in\overline{\Omega}}v(x)$.
Assume $v(y)<0$.

Assume that property $(i)$ holds. By continuity of the the operator we have
that $\mathcal{L}^{\varepsilon}_{\Omega}(v)(x)\geq 0$ for all
$x\in\overline{\Omega}$. Then

$\begin{split}0\leq&\mathcal{L}^{\varepsilon}_{\Omega}(v)(y)=v(y)-\int_{\Omega}K_{\varepsilon}(\eta-y)v(\eta)\;d\eta\\\
=&\int_{\Omega}K_{\varepsilon}(\eta-y)\left(v(y)-v(\eta)\right)\;d\eta+v(y)\int_{\Omega^{c}}K_{\varepsilon}(\eta-y)\;d\eta<0,\end{split}$
(4.3)

where we used the normalization of the kernel $K_{\varepsilon}$. Hence, this
contradiction yields $v\geq 0$.

Assume now that $(ii)$ holds. Then in this case $y\in\overline{O}$. Then again
by the continuity of the operator we obtain exactly as in (4.3) a
contradiction. Thus the Theorem is proved. ∎

###  4.2 Uniform boundedness of $u^{\varepsilon}$

In this section we will show that the sequence $u^{\varepsilon}$ is uniformly
bounded in $\varepsilon$. We will use the maximum principle again. Indeed, we
will construct functions $\Phi^{\varepsilon}$ uniformly bounded such that
$\mathcal{L}^{\varepsilon}_{\Omega}(\Phi^{\varepsilon})(x)\geq\Arrowvert
g\Arrowvert_{1}e^{-\frac{\text{dist}(x,\partial\Omega)}{\varepsilon}}$. We
will use this to prove
$\mathcal{L}^{\varepsilon}_{\Omega}\left(\Phi^{\varepsilon}-u^{\varepsilon}\right)(x)\geq
0$ which implies using the maximum principle $0\leq
u^{\varepsilon}\leq\Phi^{\varepsilon}$. The main result of this subsection is
the following.

######  Theorem 4.2.

There exists suitable constants $0<\mu<1$, $0<\gamma(\mu)<\frac{1}{3}$,
$C_{1},\;C_{2},\;C_{3}>0$ and there exists some $\varepsilon_{0}>0$ such that
the function

$\Phi^{\varepsilon}(x)=C_{3}\Arrowvert
g\Arrowvert_{1}\left(C_{1}-\left|x\right|^{2}\right)+C_{2}\Arrowvert
g\Arrowvert_{1}\left[\left(1-\frac{\gamma}{1+\left(\frac{d(x)}{\varepsilon}\right)^{2}}\right)\wedge\left(1-\frac{\gamma}{1+\left(\frac{\mu
R}{\varepsilon}\right)^{2}}\right)\right],$

for $a\wedge b\>=\min\left(a,b\right)$, $R>0$ the minimal radius of curvature
$R=\min_{x\in\partial\Omega}R(x)$ and
$d(x):=\text{dist}\left(x,\partial\Omega\right)$, satisfies
$\mathcal{L}_{\Omega}^{\varepsilon}\left(\Phi^{\varepsilon}\right)(x)\geq\Arrowvert
g\Arrowvert_{1}e^{-\frac{d(x)}{\varepsilon}}$ in $\Omega$ uniformly for all
$\varepsilon<\varepsilon_{0}$. Moreover, the solutions $u^{\varepsilon}$ of
(4.2) are uniformly bounded in $\varepsilon$.

We split the proof of this theorem in two lemmas.

######  Lemma 4.1.

Let
$C_{1}:=2\max_{x\in\overline{\Omega}}\left|x\right|^{2}+2\;\textnormal{diam}\left(\Omega\right)^{2}+4\;\textnormal{diam}\left(\Omega\right)+4$,
let $0<\varepsilon<1$. Then

$\mathcal{L}^{\varepsilon}_{\Omega}\left(C_{1}-\left|x\right|^{2}\right)\geq
2\varepsilon^{2}.$

###### Proof.

We start computing the action of $\mathcal{L}^{\varepsilon}_{\mathbb{R}^{3}}$
on $\left|x\right|^{2}$.

$\begin{split}\mathcal{L}^{\varepsilon}_{\mathbb{R}^{3}}\left[\left|\cdot\right|^{2}\right](x)=&\left|x\right|^{2}-\int_{\mathbb{R}^{3}}K_{\varepsilon}\left(\eta-x\right)\left|\eta\right|^{2}\;d\eta\\\
=&-\int_{\mathbb{R}^{3}}K_{\varepsilon}\left(\eta-x\right)\left|\eta-x\right|^{2}\;d\eta=-2\varepsilon^{2},\end{split}$

where we expanded $|\eta|^{2}=|x+(\eta-x)|^{2}$ and we used that
$\int_{\mathbb{R}^{3}}K_{\varepsilon}=1$ and the symmetry of the kernel
$K_{\varepsilon}$.

Le $D:=\textnormal{diam}(\Omega)$ and let
$B=2\max_{x\in\overline{\Omega}}\left|x\right|^{2}$ and $\beta=2D^{2}+4D+4$.
Thus, $C_{1}=B+\beta$. Then

$\begin{split}\mathcal{L}^{\varepsilon}_{\Omega}&\left(K+\beta-\left|\cdot\right|^{2}\right)(x)=(B+\beta)\int_{\Omega^{c}}K_{\varepsilon}\left(\eta-x\right)\;d\eta\;-\mathcal{L}^{\varepsilon}_{\mathbb{R}^{3}}\left[\left|\cdot\right|^{2}\right](x)-\int_{\Omega^{c}}K_{\varepsilon}\left(\eta-x\right)\left|\eta\right|^{2}\;d\eta\\\
\geq&\left(B+\beta\right)\int_{\Omega^{c}}K_{\varepsilon}\left(\eta-x\right)\;d\eta+2\varepsilon^{2}-2\left|x\right|^{2}\int_{\Omega^{c}}K_{\varepsilon}\left(\eta-x\right)\;d\eta-2\int_{\Omega^{c}}K_{\varepsilon}\left(\eta-x\right)\left|\eta-x\right|^{2}\;d\eta,\\\
\end{split}$

where we used $\left|\eta\right|^{2}\leq
2\left|x\right|^{2}+2\left|\eta-x\right|^{2}$. Moreover using that
$B-2|x|^{2}\geq 0$ and splitting for $x\in\Omega$ the complement of the domain
as $\Omega^{c}=\left(\Omega^{c}\cap B_{D}(x)\right)\cup B_{D}^{c}(x)$ we
obtain

$\begin{split}\mathcal{L}^{\varepsilon}_{\Omega}\left(K+\beta-\left|\cdot\right|^{2}\right)(x)\geq&2\varepsilon^{2}+\int_{B^{c}_{D/\varepsilon}(0)}K_{\varepsilon}\left(\eta\right)\left(\beta-2\varepsilon^{2}\left|\eta\right|^{2}\right)\;d\eta\\\
=&2\varepsilon^{2}+\beta
e^{-\frac{D}{\varepsilon}}-e^{-\frac{D}{\varepsilon}}\left(2D^{2}+4D\varepsilon+4\varepsilon^{2}\right)\geq
2\varepsilon^{2},\end{split}$

where in the first inequality we used that $2\left|\eta-x\right|^{2}\leq
2D^{2}\leq\beta$ for $\eta,x\in B_{D}(x)$ and for the integral in
$B^{c}_{D}(x)$ we changed variables $\frac{\eta-x}{\varepsilon}\mapsto\eta$
and we computed the resulting integral using also that $\varepsilon<1$. ∎

In order to proceed further with the construction of the supersolution, we
will use repeatedly the distance function and its relation to the curvature of
the domain’s boundary. All the properties of this function can be found in the
Appendix “Boundary curvatures and distance functions” in [9]. It is well-know
that if the boundary $\partial\Omega$ is $C^{3}$, then in a neighborhood of
the boundary the distance function can be expanded by Taylor as

$d(\eta)=d(x)+\nabla
d(x)\cdot(\eta-x)+\frac{1}{2}\left(\eta-x\right)^{\top}\nabla^{2}d(x)\left(\eta-x\right)+\mathcal{O}\left(\left|\eta-x\right|^{3}\right)$
(4.4)

Moreover, the following proposition holds.

######  Proposition 4.1.

For $x\in\Omega$ in a neighborhood of the boundary the gradient of the
distance function is the inner normal, so that $\left|\nabla d(x)\right|=1$.
Moreover, denoting $R=\min_{x\in\partial\Omega}R(x)>0$ the minimal radius of
curvature and letting $\mu\in(0,1)$ we have

$\xi^{\top}\nabla^{2}d(x)\xi\leq\frac{1}{(1-\mu)R}$ (4.5)

for every $x\in\left\\{y\in\Omega:\;d(y)<R\mu\right\\}$ and
$\Arrowvert\xi\Arrowvert=1$.

###### Proof.

See 14.6, Appendix “Boundary curvatures and distance functions” ([9]). ∎

Using these properties of the distance function we can prove the next lemma.

######  Lemma 4.2.

Let
$\psi(x):=\left(1-\frac{\gamma}{1+\left(\frac{d(x)}{\varepsilon}\right)^{2}}\right)\wedge\left(1-\frac{\gamma}{1+\left(\frac{\mu
R}{\varepsilon}\right)^{2}}\right)$. Then there exists some $0<\mu<1$ small
enough, $0<\gamma(\mu)<\frac{1}{3}$, $0<\varepsilon_{1}<1$ small enough and
constants $C_{0}:=C_{0}(R,\Omega,\mu,\gamma)>0$ and $c:=c(R,\mu,\gamma)>0$
such that for all $0<\varepsilon\leq\varepsilon_{1}$

$\mathcal{L}_{\Omega}^{\varepsilon}\left(\psi\right)(x)\geq\begin{cases}C_{0}e^{-\frac{d(x)}{\varepsilon}}&0<d(x)\leq\frac{R\mu}{2}\\\
-c\varepsilon^{2}&\frac{R\mu}{2}<d(x)<R\mu\\\ 0&d(x)\geq R\mu\\\ \end{cases}$
(4.6)

###### Proof.

We start with some preliminary consideration on the distance function. We
define $\frac{d(\eta)}{\varepsilon}:=d_{\varepsilon}\left(\eta\right)$. For
every $x,\eta\in\left\\{y\in\Omega:\;d(y)<R\mu\right\\}$ we have using (4.4)

$\begin{split}d_{\varepsilon}\left(\eta\right)^{2}=&d_{\varepsilon}\left(x\right)^{2}+\frac{2d(x)\nabla
d(x)\cdot\left(\eta-x\right)}{\varepsilon^{2}}+\frac{d(x)\left(\eta-x\right)^{\top}\nabla^{2}d(x)\left(\eta-x\right)}{\varepsilon^{2}}\\\
&+\frac{\left(\nabla
d(x)\cdot(\eta-x)\right)^{2}}{\varepsilon^{2}}+\mathcal{O}\left(\frac{d(x)}{\varepsilon^{2}}\left|\eta-x\right|^{3}\right).\end{split}$
(4.7)

Then Taylor’s expansion shows

$\begin{split}\frac{1}{1+d_{\varepsilon}\left(\eta\right)^{2}}=&\frac{1}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)\left(1+\left[d_{\varepsilon}\left(\eta\right)^{2}-d_{\varepsilon}\left(x\right)^{2}\right]\frac{1}{1+d_{\varepsilon}\left(x\right)^{2}}\right)}\\\
=&Q_{\varepsilon}^{(1)}(x,\eta)+Q_{\varepsilon}^{(2)}(x,\eta)+Q_{\varepsilon}^{(3)}(x,\eta),\\\
\end{split}$ (4.8)

where we the terms $Q_{\varepsilon}^{(i)}$ are defined as follows.

$Q_{\varepsilon}^{(1)}(x,\eta)=\frac{1}{1+d_{\varepsilon}\left(x\right)^{2}}-\frac{2d(x)\nabla
d(x)\cdot\left(\eta-x\right)}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}},$
$Q_{\varepsilon}^{(2)}(x,\eta)=-\frac{d(x)\left(\eta-x\right)^{\top}\nabla^{2}d(x)\left(\eta-x\right)}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}-\frac{\left(\nabla
d(x)\cdot(\eta-x)\right)^{2}}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}+\frac{4d^{2}(x)\left(\nabla
d(x)\cdot(\eta-x)\right)^{2}}{\varepsilon^{4}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{3}},$
$Q_{\varepsilon}^{(3)}(x,\eta)=\mathcal{O}\left(\frac{d(x)}{\varepsilon^{2}}\frac{\left|\eta-x\right|^{3}}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\right)+\mathcal{O}\left(\frac{d(x)}{\varepsilon^{4}}\frac{\left|\eta-x\right|^{3}}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{3}}\right).$

We consider now the function $\psi(x)$ defined in the statement of Lemma 4.2.
We take $M=\frac{1}{\mu^{2}}$ for $0<\mu<1$ small enough and $0<\varepsilon<1$
also small enough such that $0<M\varepsilon<\frac{R\mu}{2}$, i.e.
$0<\varepsilon<\frac{R\mu^{3}}{2}$, and we decompose $\Omega$ in four disjoint
sets

$\Omega=\left\\{d(x)\geq
R\mu\right\\}\cup\left\\{d(x)<M\varepsilon\right\\}\cup\left\\{M\varepsilon\leq
d(x)\leq\frac{R\mu}{2}\right\\}\cup\left\\{\frac{R\mu}{2}<d(x)<R\mu\right\\}.$

We proceed estimating $\mathcal{L}_{\Omega}^{\varepsilon}(\psi)(x)$ for $x$ in
each of these regions of $\Omega$.

Figure 4: Decomposition of $\Omega$.

For further reference we write

$\begin{split}\mathcal{L}^{\varepsilon}_{\Omega}\left(\psi\right)(x)=&\psi(x)-\int_{\Omega\cap\\{d(\eta)<R\mu\\}}d\eta
K_{\varepsilon}(\eta-x)\left(1-\frac{\gamma}{1+d_{\varepsilon}\left(\eta\right)^{2}}\right)\\\
&-\int_{\Omega\cap\\{d(\eta)\geq R\mu\\}}d\eta
K_{\varepsilon}(\eta-x)\left(1-\frac{\gamma}{1+\left(\frac{\mu
R}{\varepsilon}\right)^{2}}\right).\end{split}$ (4.9)

In order to estimate $\mathcal{L}_{\Omega}^{\varepsilon}(\psi)(x)$ in the
region $\\{d(x)\geq R\mu\\}$ we will use the fact that the minimum of
supersolutions is again a supersolution. In the region where
$d(x)<M\varepsilon$ we will use the explicit form of the kernel to see that
the main contribution has the right sign. Finally, in the region
$\\{M\varepsilon\leq d(x)<R\mu\\}$ the idea behind the arguments we will
present is that $\mathcal{L}_{\Omega}^{\varepsilon}(\psi)(x)$ can be
approximated by $-\varepsilon^{2}\Delta\psi$ using Taylor.

Step 1: $\\{d(x)\geq R\mu\\}$

First of all we notice that if $d(x)\geq R\mu$ then
$\mathcal{L}^{\varepsilon}_{\Omega}(\psi)(x)\geq 0$. Indeed,
$\psi(\eta)\leq\psi(x)=1-\frac{\gamma}{1+\left(\frac{\mu
R}{\varepsilon}\right)^{2}}$ in the first integral of (4.9) since
$d(\eta)<R\mu$ there. Hence

$\mathcal{L}^{\varepsilon}_{\Omega}\left(\psi\right)(x)\geq\mathcal{L}^{\varepsilon}_{\Omega}\left(1-\frac{\gamma}{1+\left(\frac{\mu
R}{\varepsilon}\right)^{2}}\right)\geq 0.$ (4.10)

Step 2: $\\{d(x)<M\varepsilon\\}$

We consider now the region $\\{d(x)<M\varepsilon\\}$. After a suitable rigid
motion we can assume $0\in\partial\Omega$ and $x=(d(x),0,0)$. Hence,
$\Omega\subset\mathbb{R}_{+}\times\mathbb{R}^{2}$ and

$\begin{split}\int_{\Omega^{c}}\frac{e^{-\frac{\left|{\eta-x}\right|}{\varepsilon}}}{4\pi\varepsilon\left|{\eta-x}\right|^{2}}\;d\eta&\geq\int_{-\infty}^{-d(x)/\varepsilon}K\left(\eta\right)\;d\eta\geq\int_{-\infty}^{-M}K\left(\eta\right)\;d\eta:=\nu_{M}>0.\end{split}$

$K$ is as usual the normalized exponential integral. On the other hand, using
that $\frac{1}{1+d_{\varepsilon}\left(x\right)^{2}}\leq 1$ and choosing
$\gamma<\frac{\nu_{M}}{2}$ we can conclude

$\begin{split}\mathcal{L}^{\varepsilon}_{\Omega}\left(\psi\right)(x)=&-\frac{\gamma}{1+d_{\varepsilon}\left(x\right)^{2}}+\int_{\Omega^{c}}d\eta
K_{\varepsilon}(\eta-x)+\gamma\int_{\Omega}d\eta
K_{\varepsilon}(\eta-x)\left(\frac{1}{1+d_{\varepsilon}\left(\eta\right)^{2}}\vee\frac{1}{1+\left(\frac{\mu
R}{\varepsilon}\right)^{2}}\right)\\\
\geq&\frac{\nu_{M}}{2}\geq\frac{\nu_{M}}{2}e^{-d_{\varepsilon}\left(x\right)},\end{split}$
(4.11)

where $a\vee b=\max(a,b)$.

Step 3: $\left\\{M\varepsilon\leq d(x)\leq\frac{R\mu}{2}\right\\}$

We consider now the set $\left\\{M\varepsilon\leq
d(x)\leq\frac{R\mu}{2}\right\\}$. As first step we plug (4.8) into the right
hand side of (4.9). To this end we define three integral terms
$J_{1},\;J_{2},\;J_{3}$ as

$\begin{split}J_{1}=&1-\frac{\gamma}{1+d_{\varepsilon}\left(x\right)^{2}}-\int_{\Omega\cap\\{d(\eta)<R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\left(1-\gamma
Q_{\varepsilon}^{(1)}(x,\eta)\right)\\\ &-\int_{\Omega\cap\\{d(\eta)\geq
R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\left(1-\frac{\gamma}{1+\frac{R^{2}\mu^{2}}{\varepsilon^{2}}}\right),\end{split}$
(4.12)
$J_{2}=\int_{\Omega\cap\\{d(\eta)<R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\left(\gamma
Q_{\varepsilon}^{(2)}(x,\eta)\right),$ (4.13)
$J_{3}=\int_{\Omega\cap\\{d(\eta)<R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\left(\gamma
Q_{\varepsilon}^{(3)}(x,\eta)\right).$ (4.14)

Hence, we have

$\begin{split}\mathcal{L}_{\Omega}^{\varepsilon}(\psi)(x)=J_{1}+J_{2}+J_{3}.\end{split}$
(4.15)

The main contribution to these terms is due to $J_{2}$. Therefore we start
with this term and we show that for $0<\mu<1$ small enough there exists a
constant $\tilde{C}(\mu)>0$ independent of $\varepsilon$ such that

$\begin{split}J_{2}\geq\frac{\tilde{C}(\mu)\gamma}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}.\end{split}$
(4.16)

In order to prove this estimate we first notice that

$\begin{split}\frac{4d_{\varepsilon}\left(x\right)^{2}}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)}-1=&3-\frac{4}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)}\geq
3-\frac{4}{\left(1+M^{2}\right)}\geq 0.\end{split}$ (4.17)

Hence, multiplying this inequality by
$K_{\varepsilon}(\eta-x)\frac{\gamma\left(\nabla
d(x)\cdot\left(\eta-x\right)\right)^{2}}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}$
and integrating on $\\{d(\eta)<R\mu\\}$ we obtain

$\begin{split}\int_{\Omega\cap\\{d(\eta)<R\mu\\}}&d\eta\;K_{\varepsilon}(\eta-x)\left(-\frac{\gamma\left(\nabla
d(x)\cdot\left(\eta-x\right)\right)^{2}}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}+\frac{4\gamma
d^{2}(x)\left(\nabla
d(x)\cdot\left(\eta-x\right)\right)^{2}}{\varepsilon^{4}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{3}}\right)\\\
\geq&\frac{\gamma\left(3-\frac{4}{1+M^{2}}\right)}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\int_{B_{M\varepsilon}(x)}d\eta\;K_{\varepsilon}(\eta-x)\frac{\left(\nabla
d(x)\cdot\left(\eta-x\right)\right)^{2}}{\varepsilon^{2}}\\\
=&\frac{\gamma\left(3-\frac{4}{1+M^{2}}\right)}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\frac{1}{4\pi}\int_{0}^{2\pi}d\varphi\int_{0}^{\pi}d\theta\sin(\theta)\cos^{2}(\theta)\int_{0}^{M}dr\>e^{-r}r^{2}=\frac{\gamma
C(M)\left(3-\frac{4}{1+M^{2}}\right)}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}},\end{split}$
(4.18)

where used that $B_{M\varepsilon}(x)\subset\\{d(\eta)<R\mu\\}$ and we define
the constant
$C(M)=\frac{1}{3}\int_{0}^{M}dr\;e^{-r}r^{2}=\frac{1}{3}(2-2e^{-M}-2Me^{-M}-M^{2}e^{-M})$
which depends on $M=\frac{1}{\mu^{2}}$. Notice that $C(M)\to\frac{2}{3}$ as
$M\to\infty$ and hence for $M$ sufficiently large we have also
$C(M)\geq\frac{1}{2}$.

In order to conclude the estimate for $J_{2}$ we use the result (4.5) to
estimate the Hessian of the distance function, thus

$\frac{\gamma
d(x)\left(\eta-x\right)^{\top}\nabla^{2}d(x)\left(\eta-x\right)}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\leq\frac{\gamma\mu\left|\eta-x\right|^{2}}{\varepsilon^{2}(1-\mu)\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}$
(4.19)

and we conclude

$-\int_{\Omega\cap\\{d(\eta)<R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\frac{\gamma
d(x)\left(\eta-x\right)^{\top}\nabla^{2}d(x)\left(\eta-x\right)}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\geq-C\frac{\gamma\mu}{(1-\mu)\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}},$
(4.20)

for some constant $C>0$.

Combining (4.18) and (4.20) we obtain (4.16).

We proceed now with the term $J_{1}$ in (4.12). Using the symmetry of the
scalar product in $\mathbb{R}^{3}$ we write

$\begin{split}J_{1}=&\int_{\Omega^{c}}d\eta\;K_{\varepsilon}(\eta-x)\left(1-\gamma
Q_{\varepsilon}^{(1)}(x,\eta)\right)+\int_{\Omega\cap\\{d(\eta)\geq
R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\left(\frac{\gamma}{1+\frac{R^{2}\mu^{2}}{\varepsilon^{2}}}-\gamma
Q_{\varepsilon}^{(1)}(x,\eta)\right)\\\ =&J_{1,1}+J_{1,2}.\end{split}$ (4.21)

We proceed with the estimate for $J_{1,1}$ in (4.21). By means of a suitable
coordinate system we can assume again $0\in\partial\Omega$ and $x=(d(x),0,0)$.
We notice that if $\eta\in\left(-\infty,d(x)\right)\times\mathbb{R}^{2}$ then
$\nabla d(x)\cdot\left(\eta-x\right)=\eta_{1}-d(x)\leq 0$, while if
$\eta\in\left(d(x),\infty\right)\times\mathbb{R}^{2}$ then $\nabla
d(x)\cdot\left(\eta-x\right)\geq 0$. Hence, we obtain

$\begin{split}J_{1,1}\geq&\int_{\Omega^{c}\cap\left(-\infty,d(x)\right)\times\mathbb{R}^{2}}d\eta\;K_{\varepsilon}(\eta-x)\left(1-\gamma
Q_{\varepsilon}^{(1)}(x,\eta)\right).\\\ \end{split}$ (4.22)

We now decompose the set
$\Omega^{c}\cap\left(\left(-\infty,d(x)\right)\times\mathbb{R}^{2}\right)=\left(\left(-\infty,0\right)\times\mathbb{R}^{2}\right)\cup\left(\Omega^{c}\cap\left(\left(0,d(x)\right)\times\mathbb{R}^{2}\right)\right)$.
Using that

$\frac{d(x)}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}=\frac{1}{d(x)\left(1+d_{\varepsilon}\left(x\right)^{2}\right)}-\frac{1}{d(x)\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}$
(4.23)

and since $\gamma<\frac{1}{3}$ we have
$1-\frac{\gamma}{1+d_{\varepsilon}\left(x\right)^{2}}>0$ and therefore we
obtain

$\begin{split}\int_{\left(-\infty,0\right)\times\mathbb{R}^{2}}&d\eta\;K_{\varepsilon}(\eta-x)\left(1-\gamma
Q_{\varepsilon}^{(1)}(x,\eta)\right)\geq\int_{\left(-\infty,0\right)\times\mathbb{R}^{2}}d\eta\;K_{\varepsilon}(\eta-x)\frac{2\gamma\nabla
d(x)\cdot\left(\eta-x\right)}{d_{\varepsilon}\left(x\right)\varepsilon\left(1+d_{\varepsilon}\left(x\right)^{2}\right)}\\\
=&-\frac{2\gamma}{d_{\varepsilon}\left(x\right)\left(1+d_{\varepsilon}\left(x\right)^{2}\right)}\int_{d_{\varepsilon}\left(x\right)}^{\infty}dz\;K\left(z\right)z\geq-\frac{\gamma}{2d_{\varepsilon}\left(x\right)}\frac{1+d_{\varepsilon}\left(x\right)}{1+d_{\varepsilon}\left(x\right)^{2}}e^{-d_{\varepsilon}\left(x\right)}\\\
\geq&-\frac{\gamma
C}{M}\frac{1}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}},\end{split}$
(4.24)

where we also changed variable $(d_{\varepsilon}\left(x\right)-z)\mapsto z$,
we used the identity (2.11) for the normalized exponential integral in
Proposition 2.2, we estimated $d_{\varepsilon}\left(x\right)\geq M$ and
finally we denote by $C$ the constant such that
$\frac{(1+x^{2})^{2}}{2}e^{-|x|}\leq C$.

Concerning the integral in the set
$\Omega^{c}\cap\left(\left(0,d(x)\right)\times\mathbb{R}^{2}\right)$ we
proceed similarly using again (4.23) and also the fact that if $z>0$ then
$z-d(x)>-d(x)$. Hence, we have

$\begin{split}&\int_{\Omega^{c}\cap\left(\left(0,d(x)\right)\times\mathbb{R}^{2}\right)}d\eta\;K_{\varepsilon}(\eta-x)\left(1-\gamma
Q_{\varepsilon}^{(1)}(x,\eta)\right)\\\
\geq&\int_{\Omega^{c}\cap\left(\left(0,d(x)\right)\times\mathbb{R}^{2}\right)}d\eta\;K_{\varepsilon}(\eta-x)\left(1-\frac{\gamma}{1+d_{\varepsilon}\left(x\right)^{2}}+\frac{2\gamma\nabla
d(x)\cdot\left(\eta-x\right)}{d_{\varepsilon}\left(x\right)\varepsilon\left(1+d_{\varepsilon}\left(x\right)^{2}\right)}\right)\\\
=&\int_{\Omega^{c}\cap\left((0,d(x)\times\mathbb{R}^{2}\right)}dz\;K_{\varepsilon}(\eta-d(x)e_{1})\left(1-\frac{\gamma}{1+d_{\varepsilon}\left(x\right)^{2}}+\frac{2\gamma(\eta_{1}-d(x))}{d(x)\left(1+d_{\varepsilon}\left(x\right)^{2}\right)}\right)\geq
0\end{split}$ (4.25)

Hence, for $M\varepsilon\leq d(x)<R\mu$ and $\gamma<\frac{1}{3}$ we can
summarize

$J_{1,1}\geq-\frac{\gamma}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\frac{C}{M}.$
(4.26)

######  Remark.

Notice that the estimates (4.22)-(4.26) are valid in the whole region
$\\{M\varepsilon\leq d(x)<R\mu\\}$.

We still have to consider the integral $J_{1,2}$ in (4.21). We notice that for
all $\eta\in\Omega$ with $d(\eta)\geq R\mu$ we have on the one hand
$\left|\eta-x\right|\geq\frac{R\mu}{2}$ and on the other hand $\nabla
d(x)\cdot(\eta-x)\geq 0$ since $d(\eta)>d(x)$. We recall that
$D:=\textnormal{diam}\left(\Omega\right)$ and that $\Omega\cap\\{d(\eta)\geq
R\mu\\}\subset B_{D}(x)$. Therefore, we estimate

$\begin{split}J_{1,2}\geq&-\int_{\Omega\cap\\{d(\eta)\geq
R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\frac{\gamma}{1+d_{\varepsilon}\left(x\right)^{2}}\geq-\frac{\gamma
e^{-\frac{R\mu}{2\varepsilon}}}{1+d_{\varepsilon}\left(x\right)^{2}}\int_{B_{D}(0)}dz\;\frac{1}{4\pi\varepsilon|z|^{2}}\\\
\geq&-\gamma\frac{e^{-\frac{d_{\varepsilon}\left(x\right)}{2}}}{1+d_{\varepsilon}\left(x\right)^{2}}\frac{4D}{R\mu}\geq-\gamma
C\frac{D}{R}\frac{\mu}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\end{split}$
(4.27)

where we used the well-known estimate $xe^{-x}\leq e^{-1}$ combined with
$e^{-\frac{R\mu}{4\varepsilon}}\leq e^{-\frac{d(x)}{2\varepsilon}}$ and we
denoted by $C$ the constant such that $4x(1+x^{2})e^{-\frac{x}{2}}\leq C$ and
finally the relation $M=\frac{1}{\mu^{2}}$.

Finally we estimate the term $J_{3}$ in (4.14). Here we have to estimate the
integral term containing the error terms $Q_{\varepsilon}^{(3)}(x,\eta)$ of
the Taylor expansion (4.7). If $M\varepsilon\leq d(x)\leq\frac{R\mu}{2}$ and
if $\varepsilon<1$ we use
$\frac{x}{1+x^{2}}=\frac{1}{x}-\frac{1}{x\left(1+x^{2}\right)}$ and we
calculate

$\begin{split}\gamma\int_{\Omega\cap\\{d(\eta)<R\mu\\}}d\eta\;&K_{\varepsilon}(\eta-x)\left(\frac{d(x)}{\varepsilon^{2}}\frac{\left|\eta-x\right|^{3}}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}+\frac{d(x)}{\varepsilon^{4}}\frac{\left|\eta-x\right|^{3}}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{3}}\right)\\\
\leq&\int_{\mathbb{R}^{3}}d\eta\;\frac{\gamma
e^{-|\eta|}}{4\pi}\frac{\left|\eta\right|}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\left(d(x)\varepsilon+\frac{1}{\frac{d(x)}{\varepsilon}}-\frac{1}{\frac{d(x)}{\varepsilon}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)}\right)\\\
\leq&\frac{C\gamma}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\left(\frac{R\mu}{2}+\mu^{2}\right).\end{split}$
(4.28)

Hence, also
$J_{3}\geq-\frac{C\gamma}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\left(\frac{R\mu}{2}+\mu^{2}\right)$.

We conclude putting together estimates (4.16) (4.21), (4.26), (4.27) and
(4.28) the existence of a constant $C(\Omega)>0$ independent of
$\mu,\gamma,\varepsilon$ such that

$\mathcal{L}_{\Omega}^{\varepsilon}\left(\psi\right)(x)\geq\frac{\gamma}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\left[C(M)\left(3-\frac{4}{1+M^{2}}\right)-C(\Omega)\frac{\mu}{1-\mu}\right].$
(4.29)

Choosing $0<\mu<1$ small enough, depending only on $\Omega$, such that
$C(M)>\frac{1}{3}$ and $C(\Omega)\frac{\mu}{1-\mu}<\frac{1}{6}$ we obtain

$\mathcal{L}_{\Omega}^{\varepsilon}\left(\psi\right)(x)\geq\frac{\gamma}{6\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\geq
Ce^{-\frac{d(x)}{\varepsilon}}$ (4.30)

for $M\varepsilon\leq d(x)\leq\frac{R\mu}{2}$ and some constant $C$ depending
on $\Omega$, $R$, $\gamma$, $\mu$ but independent of $\varepsilon$.

Step 4: $\left\\{\frac{R\mu}{2}<d(x)<R\mu\right\\}$

It remains to calculate the behavior of
$\mathcal{L}_{\Omega}^{\varepsilon}(\psi)$ when $\frac{R\mu}{2}<d(x)<R\mu$.
Here, we show that there exists a constant $c(R,\mu,\gamma)$ such that
$\mathcal{L}_{\Omega}^{\varepsilon}\left(\psi\right)(x)\geq-c\varepsilon^{2}$.
We can use several results we obtained in Step 3. We decompose again the
operator $\mathcal{L}_{\Omega}^{\varepsilon}(\psi)(x)=J_{1}+J_{2}+J_{3}$
according to (4.15) using the integral terms defined in (4.12)-(4.14).

First of all (4.17) implies

$\int_{\Omega\cap\\{d(\eta)<R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\left(-\frac{\gamma\left(\nabla
d(x)\cdot\left(\eta-x\right)\right)^{2}}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}+\frac{4\gamma
d^{2}(x)\left(\nabla
d(x)\cdot\left(\eta-x\right)\right)^{2}}{\varepsilon^{4}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{3}}\right)\geq
0$

and hence we estimate $J_{2}$ using (4.19) and (4.20)

$\begin{split}J_{2}\geq&-\int_{\Omega\cap\\{d(\eta)<R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\frac{\gamma
d(x)\left(\eta-x\right)^{\top}\nabla^{2}d(x)\left(\eta-x\right)}{\varepsilon^{2}\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\\\
\geq&-C\frac{\gamma\mu}{(1-\mu)\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\geq-\frac{8\gamma
C}{(1-\mu)R^{3}}\varepsilon^{3},\end{split}$ (4.31)

where we used $1+d_{\varepsilon}\left(x\right)^{2}\geq
d_{\varepsilon}\left(x\right)^{2}\geq\left(\frac{R\mu}{2\varepsilon}\right)^{2}$
and $0<\varepsilon<\frac{R\mu^{3}}{2}$.

We now proceed to estimate $J_{1}$. To this end we use again the decomposition
(4.21). The estimate (4.26) for $J_{1,1}$ is also valid in the region
$\\{\frac{R\mu}{2}<d(x)<R\mu\\}$, as we indicated in the remark after (4.26).
Hence we have for $\varepsilon<\frac{R\mu^{3}}{2}$

$\begin{split}J_{1,1}\geq-\frac{\gamma\mu^{2}C}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}\geq-\frac{8\gamma
C}{R^{2}}\varepsilon^{3}.\end{split}$

Concerning the term $J_{1,2}$ we have to argue slightly different than in Step
3. Using now the first inequality in (4.27) and $\int_{\mathbb{R}^{3}}d\eta
K_{\varepsilon}(\eta-x)=1$ we compute

$\begin{split}J_{1,2}\geq&-\int_{\Omega\cap\\{d(\eta)\geq
R\mu\\}}d\eta\;K_{\varepsilon}(\eta-x)\frac{\gamma}{1+d_{\varepsilon}\left(x\right)^{2}}\geq-\frac{\gamma}{1+d_{\varepsilon}\left(x\right)^{2}}\geq-\frac{4\gamma}{\left(R\mu\right)^{2}}\varepsilon^{2}.\end{split}$
(4.32)

Finally, we estimate $J_{3}$ as defined in (4.14). Arguing as in (4.28) and
using $1+x^{2}\geq x^{2}$ and $0<\varepsilon<\frac{R\mu^{3}}{2}$ we compute

$\begin{split}\int_{\Omega\cap\\{d(\eta)<R\mu\\}}d\eta\;&K_{\varepsilon}(\eta-x)\left(\frac{d(x)}{\varepsilon^{2}}\frac{\left|\eta-x\right|^{3}}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{2}}+\frac{d(x)}{\varepsilon^{4}}\frac{\left|\eta-x\right|^{3}}{\left(1+d_{\varepsilon}\left(x\right)^{2}\right)^{3}}\right)\\\
\leq&\frac{\gamma\left(d(x)^{2}+1\right)}{4\pi\left(d_{\varepsilon}\left(x\right)\right)^{5}}\int_{\mathbb{R}^{3}}d\eta\;e^{-|\eta|}\left|\eta\right|\leq\frac{2\gamma
C\left({R^{2}}+2\right)}{R^{3}}\varepsilon^{2}.\end{split}$ (4.33)

Thus, also $J_{3}\geq-\frac{2\gamma
C\left({R^{2}}+2\right)}{R^{3}}\varepsilon^{2}$.

Hence, (4.31),(4.2),(4.32) and (4.33) imply the existence of a constant
$c(R,\mu,\gamma)>0$ independent of $\varepsilon$ such that

$\mathcal{L}_{\Omega}^{\varepsilon}\left(\psi\right)(x)\geq-c\varepsilon^{2}$
(4.34)

for all $\frac{R\mu}{2}<d(x)<R\mu$ .

We know summarize the results. Equations (4.10), (4.11), (4.30), (4.34) imply
the claim in (4.6). We remark that $\mu$, $\gamma$ and $\varepsilon_{1}$ are
chosen as follows. First of all $\mu$ is chosen according to Step 3 as in
(4.29), then $\gamma$ is taken according to Step 2 such that
$0<\gamma<\frac{\nu_{M}}{2}$ and finally $\varepsilon_{1}$ satisfies
$0<\varepsilon_{1}<\frac{R\mu^{3}}{2}$. This concludes the Lemma 4.2. ∎

Using Lemma 4.1 and 4.2 we can now prove Theorem 4.2.

###### (Proof of Theorem 4.2).

Let $C_{1}$ be the constant defined in Lemma 4.1 and let
$\gamma,\;\mu,\;C_{0},\;c$ be as in Lemma 4.2. We define
$C_{2}:=\frac{1}{C_{0}}$ and $C_{3}:=\frac{C_{0}+c}{2C_{0}}>\frac{1}{2}$.
Notice that all these constants are independent of $\varepsilon$. Hence, Lemma
4.1 and 4.2 imply

$\mathcal{L}_{\Omega}^{\varepsilon}\left(\Phi^{\varepsilon}\right)(x)\geq\Arrowvert
g\Arrowvert_{1}\begin{cases}e^{-\frac{d(x)}{\varepsilon}}+2C_{3}\varepsilon^{2}&0<d(x)\leq\frac{R\mu}{2},\\\
\varepsilon^{2}&\frac{R\mu}{2}<d(x)<R\mu,\\\ 2C_{3}\varepsilon^{2}&d(x)\geq
R\mu,\\\ \end{cases}\geq\Arrowvert
g\Arrowvert_{1}\begin{cases}e^{-\frac{d(x)}{\varepsilon}}&0<d(x)\leq\frac{R\mu}{2},\\\
\varepsilon^{2}&\frac{R\mu}{2}<d(x)<R\mu,\\\ \varepsilon^{2}&d(x)\geq R\mu.\\\
\end{cases}$ (4.35)

We define now $\varepsilon_{0}:=\min\left\\{1,\;a,\varepsilon_{1}\right\\}$
with $a$ such that $2a\ln(\frac{1}{a})<\frac{R\mu}{2}$ and $\varepsilon_{1}>0$
as in Lemma 4.2. Then $\varepsilon^{2}\geq e^{-\frac{R\mu}{2\varepsilon}}\geq
e^{-\frac{d(x)}{\varepsilon}}$ for all $d(x)>\frac{R\mu}{2}$.

We now apply the maximum principle in Theorem 4.1 to the function
$\Phi^{\varepsilon}-u^{\varepsilon}$. This function satisfies the continuity
and boundedness assumption. Indeed, for any $\varepsilon>0$ the function
$u^{\varepsilon}$ is continuous and bounded as we have seen at the beginning
of Section 4.1. Moreover, by construction $\Phi^{\varepsilon}$ is continuous
and it is easy to see that it is even uniformly bounded since

$0\leq\Phi^{\varepsilon}(x)\leq\Arrowvert
g\Arrowvert_{1}\left(2C_{3}C_{1}+C_{2}\right).$

We also have

$\mathcal{L}_{\Omega}^{\varepsilon}\left(\Phi^{\varepsilon}-u^{\varepsilon}\right)(x)\geq\Arrowvert
g\Arrowvert_{1}e^{-\frac{d(x)}{\varepsilon}}-\int_{0}^{\infty}d\nu\int_{n\cdot
N_{x_{\Omega}}<0}dn\;g_{\nu}(n)e^{-\frac{\left|x-x_{\Omega}(x,n)\right|}{\varepsilon}}\geq
0,$

since $\left|x-x_{\Omega}(x,n)\right|\geq d(x)$. Hence, Theorem 4.1 implies
that $\Phi^{\varepsilon}-u^{\varepsilon}\geq 0$ and thus

$0\leq u^{\varepsilon}\leq\Phi^{\varepsilon}\leq\tilde{C}<\infty$

uniformly in $\varepsilon$ and $x\in\Omega$. ∎

###  4.3 Estimates of $u^{\varepsilon}-\overline{u}$ near the boundary
$\partial\Omega$

In this subsection we will prove that for each point $p\in\partial\Omega$ the
function $\overline{u}$ defined in (2.8) is a good approximation of
$u^{\varepsilon}$ in a neighborhood of size close to
$\varepsilon^{\frac{1}{2}}$. Notice that this neighborhood is much greater
than the region of size $\varepsilon$. We will do it by means of the maximum
principle in Theorem 4.1. Now we start estimating the action of the operator
$\mathcal{L}_{\Omega}^{\varepsilon}$ on $\overline{u}-u^{\varepsilon}$.

######  Lemma 4.3.

Let $p\in\partial\Omega$ and let $\mathcal{R}_{p}$ be the isometry defined in
(1.12). Then the following holds for $x\in\Omega$, $\delta>0$ sufficiently
small and independent of $\varepsilon$ and a suitable $0<A<1$ and constant
$C>0$

$\left|\mathcal{L}_{\Omega}^{\varepsilon}\left(\overline{u}\left(\frac{\mathcal{R}_{p}(\cdot)\cdot
e_{1}}{\varepsilon},p\right)-u^{\varepsilon}\right)(x)\right|\leq
Ce^{-\frac{Ad(x)}{\varepsilon}}\begin{cases}\varepsilon^{\delta}&\text{ if
}|x-p|<\varepsilon^{\frac{1}{2}+2\delta},\\\ 1&\text{ if
}|x-p|\geq\varepsilon^{\frac{1}{2}+2\delta}.\end{cases}$ (4.36)

###### Proof.

Let us denote by $\Pi_{p}$ the half space
$\Pi_{p}:=\mathcal{R}_{p}^{-1}\left(\mathbb{R}_{+}\times\mathbb{R}^{2}\right)$.
Then the function
$\overline{U}_{\varepsilon}(x,p):=\overline{u}\left(\frac{\mathcal{R}_{p}(x)\cdot
e_{1}}{\varepsilon},p\right)$ is a continuous bounded function which maps
$\Pi_{p}\times\partial\Omega$ to $\mathbb{R}_{+}$. Notice that
$\overline{U}_{\varepsilon}(x,p)$ is the solution to the planar equation (2.8)
before rescaling and rotating. Our plan is to approximate
$\mathcal{L}_{\Omega}^{\varepsilon}\left(\overline{U}_{\varepsilon}\right)$ by
$\mathcal{L}_{\Pi_{p}}^{\varepsilon}\left(\overline{U}_{\varepsilon}\right)$.
Let $x\in\Pi_{p}$ and $p\in\Omega$. Using the definition of $\overline{u}$ in
(2.8) we can compute

$\begin{split}\int_{0}^{\infty}d\eta\;&K\left(\eta-\frac{\mathcal{R}_{p}(x)\cdot
e_{1}}{\varepsilon}\right)\overline{u}(\eta,p)=\int_{\mathbb{R}_{+}\times\mathbb{R}^{2}}d\eta\;\frac{e^{-\left|\eta-\frac{\mathcal{R}_{p}(x)}{\varepsilon}\right|}}{4\pi\left|\eta-\frac{\mathcal{R}_{p}(x)}{\varepsilon}\right|^{2}}\overline{u}\left(\eta_{1},p\right)\\\
&=\int_{\mathbb{R}_{+}\times\mathbb{R}^{2}}d\eta\;\frac{e^{-\frac{\left|\eta-\mathcal{R}_{p}(x)\right|}{\varepsilon}}}{4\pi\varepsilon\left|\eta-\mathcal{R}_{p}(x)\right|^{2}}\overline{u}\left(\frac{\eta_{1}}{\varepsilon},p\right)=\int_{\Pi_{p}}d\eta\;K_{\varepsilon}(\eta-x)\overline{U}_{\varepsilon}(\eta,p),\end{split}$

where we used in the first equality the translation invariance of the integral
with respect to the second and third variable, the definition of the planar
kernel and the definition of $y$. For the second equality we used the change
of variables $\tilde{\eta}=\varepsilon\eta$ and in the last identity the
change of variables $\tilde{\eta}=\mathcal{R}_{p}^{-1}(\eta)$ gives the
result. In order to write the value of
$\mathcal{L}_{\Pi_{p}}^{\varepsilon}\left(\overline{U}_{\varepsilon}\right)$
we use once again equation (2.8) and we define $x_{\Pi_{p}}(x,n)$ as the point
on the boundary of $\Pi_{p}$ with
$\frac{x-x_{\Pi_{p}}(x,n)}{\left|x-x_{\Pi_{p}}(x,n)\right|}=n$, i.e.
$x=x_{\Pi_{p}}(x,n)+\left|x-x_{\Pi_{p}}(x,n)\right|n$ if $n\cdot N_{p}<0$. By
construction we see that $\frac{\mathcal{R}_{p}(x)\cdot e_{1}}{|n\cdot
N_{p}|}=\left|x-x_{\Pi_{p}}(x,n)\right|$. Hence,

$\mathcal{L}_{\Pi_{p}}^{\varepsilon}\left(\overline{U}_{\varepsilon}(\cdot,p)\right)(x)=\int_{0}^{\infty}d\nu\int_{n\cdot
N_{p}<0}dn\;g_{\nu}(n)e^{-\frac{\left|x-x_{\Pi_{p}}(x,n)\right|}{\varepsilon}}.$