# A Nonlocal Laplacian-Based Model for Bituminous Surfacing Crack Recovery and its MPI Implementation

Noémie Debroux, Carole Le Guyader, Luminita A Vese

► To cite this version:

Noémie Debroux, Carole Le Guyader, Luminita A Vese. A Nonlocal Laplacian-Based Model for Bituminous Surfacing Crack Recovery and its MPI Implementation. Journal of Mathematical Imaging and Vision, 2020, 62 (6-7), pp.1007-1033. 10.1007/s10851-020-00968-3. hal-03034349

HAL Id: hal-03034349

https://uca.hal.science/hal-03034349

Submitted on 1 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Nonlocal Laplacian-Based Model for Bituminous Surfacing Crack Recovery and its MPI Implementation

Noémie Debroux · Carole Le Guyader · Luminita A. Vese

Received: date / Accepted: date

**Abstract** This paper is devoted to the challenging problem of fine structure detection with applications to bituminous surfacing crack recovery. Drogoul (2014) shows that such structures can be suitably modelled by a sequence of smooth functions whose Hessian matrices blow up in the perpendicular direction to the crack, while their gradient is null. This observation serves as the basis of the introduced model that also handles the natural dense and highly oscillatory texture exhibited by the images: we propose weighting $\left|\frac{\partial^2 u}{\partial x_1^2}\right|^2 + \left|\frac{\partial^2 u}{\partial x_2^2}\right|^2, u$ denoting the reconstructed image, by a variable that annihilates great expansion of this quantity, making then a connection with the elliptic approximation of the Blake-Zisserman functional. Extending then the ideas developed in the case of first order nonlocal regularization to higher order derivatives, we derive and analyze a nonlocal version of the model, and provide several theoretical results among which a $\Gamma$-convergence result as well as a detailed algorithmic approach and an MPI implementation based on a natural domain decomposition approach.

This article is dedicated to our wonderful friend and researcher, Mila Nikolova.

Noémie Debroux

Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, UK

E-mail: nd448@cam.ac.uk

Carole Le Guyader

Normandie Univ, Institut National des Sciences Appliquées de Rouen, Laboratory of Mathematics, 76000 Rouen, France
E-mail: carole.le-guyader@insa-rouen.fr

Luminita A. Vese

Department of Mathematics, University of California Los Angeles, Hilgard Avenue, Los Angeles, CA 90095-1555, USA

E-mail: lvese@math.ucla.edu

**Keywords** $\Gamma$-convergence · viscosity solutions · nonlocal second order operators · Message Passing Interface (MPI) parallelization

**Mathematics Subject Classification (2010)** 35D40 · 45E · 49J45 · 65D18 · 68W10 · 68U10 · 90C06

## 1 Introduction

Segmentation is a cornerstone step in image processing that attempts to reproduce the ability of human beings to track down significant patterns and automatically gather them into relevant and identified structures. More precisely, image segmentation consists in identifying meaningful constituents of a given image (e.g., homogeneous regions, edges, textures, etc.) for quantitative analysis or visualization purposes.

As emphasized by Zhu et al. ([46]), this task is challenging and ill-posed since the definition of an ‘object’ encompasses various acceptations: it can be something material —a thing—or a periodic pattern, this heterogeneity entailing the design of suitable methodologies for each specific application. For a relevant overview of the existing literature on this topic, we refer the reader to [46] in which an exhaustive classification of segmentation methods into three main categories (*fully supervised methods*, *semi/weakly supervised methods*, *unsupervised methods*) is provided, combined with a description / analysis of each methodology. In this paper, we focus on unsupervised methods in a variational formulation and phrased in the continuous domain.

While a contour is classically viewed as the boundary of a non zero volume object and is thus defined as the set of intensity discontinuities with jump (detected using the spatial gradient information carried by the image), fine structures exhibit discontinuities without jump. To
exemplify this latter statement, we can think of a curve
$\gamma \subset \Omega$, $\Omega$ being an open bounded subset of $\mathbb{R}^2$, and
a related image $I$ defined by $I = 1$ on $\gamma$ and $I = 0$
on $\Omega \setminus \bar{\gamma}$. The discontinuity related to the object of zero
Lebesgue measure $\{I = 1\}$ is a discontinuity without
jump (see [20,21] for further details). These thin/fila-
ment like structures cannot be captured efficiently with
differential operators of order 1. In [20,21], Drogoul
provides a heuristic illustration of this observation as
follows. Let $f$ be the function defined by $f(x) = 0$ if
$x \neq 0$ and $f(x) = 1$ if $x = 0$. This function can be
approximated by the sequence of functions ($f_\varepsilon$) such
that $f_\varepsilon(x) = \begin{cases} 0 & \text{if } x \ge \varepsilon \\ \frac{2}{\varepsilon^3} |x|^3 - \frac{3}{\varepsilon^2} |x|^2 + 1 & \text{if } x < \varepsilon \end{cases}$. On the one
hand, $f'_\varepsilon(0) = 0$, showing that the differential operator
of order 1 does not recover the singularity at 0. On the
other hand, $f''_\varepsilon(0) = -\frac{6}{\varepsilon^2}$, meaning that $f''_\varepsilon$ becomes
singular at 0. This observation is then mathematically
formalized in the two dimensional case through the fol-
lowing lemma taken from [20,21] that we recall for the
sake of completeness. Let $\mathcal{D}(\mathbb{R}^2)$ be the space of $C^\infty$
functions with compact support in $\mathbb{R}^2$ and let $\mathcal{D}'(\mathbb{R}^2)$
be the space of distributions on $\mathbb{R}^2$.

**Lemma 1 (Lemma 2.1 of [20])** Let $\varphi : \mathbb{R}^2 \to \mathbb{R}$ be a Lipschitz continuous function, and let $(g_h)_{h>0}$ be a sequence of functions defined by $g_h(x) = \frac{1}{\theta_1(h)} e^{-\frac{\varphi^2(x)}{\theta_2(h)}}$, where $\theta_i : \mathbb{R}^+ \to \mathbb{R}^+$ and $\lim_{h \to 0} \theta_i(h) = 0$.

Let $\Gamma$ be a smooth closed curve or a smooth infinite
curve of $\mathbb{R}^2$ delimiting two subdomains $\mathbb{R}_\Gamma^{2-}$ and $\mathbb{R}_\Gamma^{2+}$
forming a partition of $\mathbb{R}^2$. Let then take $\varphi$ as the signed
distance function to $\Gamma$ defined by $\varphi(x) = \text{dist}(x, \mathbb{R}_\Gamma^{2-}) -
\text{dist}(x, \mathbb{R}_\Gamma^{2+})$. $\mathbb{R}_\Gamma^{2+}$ (resp. $\mathbb{R}_\Gamma^{2-}$) identifies to the sub-domain
$\{\varphi > 0\}$ (resp. $\{\varphi < 0\}$). Taking the following scalings
$\theta_1(h) = \sqrt{\pi h}$ and $\theta_2(h) = h$, one has $g_h \xrightarrow{h \to 0} \delta_\Gamma$

in $\mathcal{D}'(\mathbb{R}^2)$. Moreover for all $x \in \Gamma$, $\nabla g_h(x) = [0,0]^T$
and $\text{spec}(\nabla^2 g_h(x)) = \{-\frac{2}{\sqrt{\pi h^{3/2}}}, 0\}$, with $\text{spec}(M)$
denoting the eigenvalues of the matrix $M$. The associated eigenvectors on $\Gamma$ are $(\nabla\varphi(x), \nabla\varphi(x)^\perp)$, where
$\nabla\varphi(x) = \mathbf{n}$ is the unit normal of $\Gamma$.

Since filaments/cracks can be modelled by a Dirac
distribution supported by a smooth curve $\Gamma$, it results
from this lemma that any such structure can be ap-
proximated by a sequence of smooth functions whose
Hessian matrices blow up in the perpendicular direc-
tion to $\Gamma$ while their gradient are null. This observation
serves as the core of the introduced material that falls
within second order variational models, includes an ad-
ditional variable encoding the singularities of the quan-
tity $\left|\frac{\partial^2 u}{\partial x_1^2}\right|^2 + \left|\frac{\partial^2 u}{\partial x_2^2}\right|^2$ (object we aim to recover) —the
cracks being mainly oriented in the vertical or horizon-

tal direction—, and that also handles the natural dense
and highly oscillatory texture exhibited by the images.
More precisely, we propose weighting $\left|\frac{\partial^2 u}{\partial x_1^2}\right|^2 + \left|\frac{\partial^2 u}{\partial x_2^2}\right|^2$, $u$
denoting the reconstructed image function, by a varia-
ble $v_2$ with values in $[0,1]$ annihilating great expansion
of this quantity so that the set of cracks is now mod-
elled by $\{v_2 = 0\}$. This leads to a variant of the elliptic
approximation of the Blake-Zisserman functional ([10])
given by Ambrosio, Faina and March ([2]).

Before introducing and analyzing in depth the proposed
model, we review some prior works related to thin struc-
ture detection and enhancement. Recently, the question
of differentiating objects within an image by contrast,
size or structure —the size being the feature of interest
in our case —has been addressed in [27]. It relies on the
spectral total-variation transform that yields a continu-
ous, multi-scale, fully edge-preserving, local descriptor,
referred to as spectral total variation local scale sig-
natures. These signatures demonstrate nice properties
such as sensitivity to size, local contrast and compo-
sition of structures, as well as invariance to rotation,
translation, flip and linear illumination change. In [24],
Frangi *et al.* tackle the issue of vessel enhancement with
the eventual goal of vessel segmentation in a multiscale
approach. It is based on the construction of a vessel
likelihood measure that relies on the computation of the
eigenvalues of the Hessian. In [6,20], Aubert and Dro-
goul introduce a topological gradient-based method for
the detection of fine structures in 2D, the major dif-
ference with our model lying in the introduction of
an auxiliary variable encoding the crack-type singulari-
ties. Still in a variational setting, Bergounioux and Vi-
cente ([9]) propose a model derived from the Mumford-
Shah functional and constrained by geometrical pri-
ors to perform the segmentation of tube-like structures
with small diameter. A limitation of this method rests
on its inability to handle junctions of tubes. Other va-
riational approaches have been investigated, dedicated
to specific applications, such as automatic road network
extraction in [36] in which nonlocal regularizers enfor-
cing straightness are applied, or detection/completion
of thin filaments in noisy blurred 2D or 3D images in [8],
relying on the Ginzburg-Landau energy. In [4], Aubert
et al. propose detecting image singularities of codimen-
sion greater or equal to 2, inspired again by Ginzburg-
Landau models.

Obviously, the range of the methods addressing this is-
sue of thin structure detection does not limit to the
variational framework: morphological approaches are
developed in [45] for automatic detection of vessel-like
patterns, though sensitive to the noise type and time-
consuming. In [41], stochastic methods in which a thin
network is simulated by a point process penalizing dis-connected segments and favoring aligned pieces are elabo-rated.

We emphasize that this paper is an extension of [18] not reduced to a simple supplement, containing no overlap with [18] and including significant theoretical and computational aspects. It completes [18] since it includes substantial results, such as a Γ-convergence result (main result) relating the local basis functional with the derived nonlocal one, and an exhaustive and substantive PDE analysis, in the viscosity solution theory framework, of the resulting evolution equation satisfied by the variable Q introduced hereafter. It enriches [18] with the derivation of an MPI parallelization of the designed algorithm inspired by the Nonlocal Means ([13]), this algorithmic part being entirely new and numerically analyzed (in particular, the choice of the weights is discussed).

The next section (Section 2) is dedicated to the mathematical material. The original local functional intertwining decomposition and fine structure detection is first introduced. Inspired then by the works of Bourgain et al. ([11]), Dávila ([17]) or Ponce ([35]) devoted to the design of nonlocal counterparts of Sobolev and BV semi-norms using radial mollifiers and integral operators, we extend their results to the case of second order operators (also inspired by the work of Lellmann et al. ([30])) and rephrase our original model in the obtained nonlocal setting. The main theoretical result —Γ-convergence result—connecting the nonlocal minimization problem to the local one is provided. Section 3 deals with the nonlocal algorithm that is devised, motivated by the Nonlocal Means approach ([13]), to include additional information related to the image content via some weights. Optimality conditions are obtained, and a thorough analysis, in the viscosity solution framework, of the PDE resulting from the evolution equation in variable Q, including existence/uniqueness and regularity in space and time is carried out.

Some limitations regarding the theoretical expression of the nonlocal component are then highlighted, which leads us to introduce a numerical compromise removing in particular the assumption of radiality of the kernel. The resulting algorithm is depicted and to improve the computation efficiency, an MPI parallelization of the code is proposed. This parallelization is driven by the natural geometry of the problem, making the partition of the image domain into subdomains supporting simultaneous local computations relevant. The paper concludes with Section 4 that intends to qualitatively assess the proposed algorithm on several experiments in comparison to prior works.

## 2 Mathematical modelling of the problem

### 2.1 Mathematical preliminaries

Let us first introduce some notations that will be useful throughout the paper. The image domain is denoted by $\Omega$ which is assumed to be a connected, bounded, open subset of $\mathbb{R}^2$ of class $C^1$. Let $f: \bar{\Omega} \to \mathbb{R}$ be the two dimensional image. For theoretical purposes, we suppose that it is in $L^\infty(\Omega)$, which is not restrictive in practice. Let $(e_1, e_2)$ be the canonical basis of $\mathbb{R}^2$. We use $dx$ ($x = (x_1, x_2)$) for integration with respect to the Lebesgue measure on $\mathbb{R}^2$ and $dt, ds, dh$ for various integrations with respect to the Lebesgue measure on $\mathbb{R}$. The differentiation indices will be a pair $\alpha = (\alpha_1, \alpha_2)$, where $\alpha_i$ is the order of the partial derivative in the variable $x_i$, and the total order of the derivative is denoted by $|\alpha| = \alpha_1 + \alpha_2$. We will use the shortened notation $D^\alpha u = \frac{\partial^{|\alpha|} u}{\partial x_1^{\alpha_1} \partial x_2^{\alpha_2}}$. Given an integer $j \ge 0$, we define the family of spaces $C_b^j(\mathbb{R}^2)$ ([19, Definition 2.2.1, p. 69]) by setting

$$C_b^j(\mathbb{R}^2) = \{u \in C^j(\mathbb{R}^2) \mid \forall \alpha \in \mathbb{N}^2, |\alpha| \le j, \exists K_\alpha, \|D^\alpha u\|_{L^\infty(\mathbb{R}^2)} \le K_\alpha\}$$

For a positive real number $0 < \lambda < 1$, the subspace $C_b^{j,\lambda}(\mathbb{R}^2)$ consists of the functions in $C_b^j(\mathbb{R}^2)$ such that if $|\alpha| \le j$, then

$$\exists C_{\alpha, \lambda}, \forall x, y \in \mathbb{R}^2, |D^\alpha u(x) - D^\alpha u(y)| \le C_{\alpha, \lambda} |x - y|^\lambda.$$

We also recall the definition of Sobolev spaces for $k \in \mathbb{N}$, and $1 \le p \le \infty$, $W^{k,p}(\Omega) = \{u \in L^p(\Omega) \mid D^\alpha u \in L^p(\Omega), \forall |\alpha| \le k\}$, and the subspaces $W_0^{k,p}(\Omega) = \{u \in W^{k,p}(\Omega) \mid \gamma_0(u) = 0\}$, with $\gamma_0$ the trace operator. We then introduce the Hilbert space $H(div)$ defined as $H(div) = \{\sigma \in (L^2(\Omega))^2 \mid \mathrm{div}\sigma \in L^2(\Omega)\}$ endowed with the scalar product $\langle\sigma_1, \sigma_2\rangle_{H(div)} = \langle\sigma_1, \sigma_2\rangle_{L^2(\Omega)}^2 + \langle\mathrm{div}\sigma_1, \mathrm{div}\sigma_2\rangle_{L^2(\Omega)}$, $\forall (\sigma_1, \sigma_2) \in H(div)^2$. Finally, we need the operator called normalized infinity Laplacian and defined as $\Delta_\infty(u) = \langle \frac{Du}{|Du|}, D^2 u \frac{Du}{|Du|} \rangle$. We now turn to the motivation and the depiction in depth of our model.

### 2.2 Original local basis model

As emphasized in the introduction, the bituminous surfacing images we consider in this work exhibit a natural dense and highly oscillatory texture (corresponding to white and black spots) that appears to be unnecessary and superfluous for the crack detection task. We thus want to discriminate and separate the texture component $v$ that catches redundant and oscillatory patterns from the simplified version of the image
denoted by $u$ in which useless information has been removed, while keeping the crack. We therefore introduce a joint decomposition/segmentation model in which the thin-structure recognition task operates on the component $u$. Texture modelling consists in finding the best functional space to represent the oscillatory patterns and has been extensively studied, either on its own or as a close counterpart of image denoising (see [25], [26], [29], [31], [32], [33], [34], [37], [39], [38] or [44]). In this work, to model the texture, we propose to use the space $G(\mathbb{R}^2)$ introduced by Meyer ([32]), space of distributions $v$ that can be written as $v = \text{divg}$ where $g = (g_1, g_2) \in (L^\infty(\mathbb{R}^2))^2$ and endowed with the norm defined by

$$ \|v\|_{G(\mathbb{R}^2)} = \inf \left\{ \|\sqrt{g_1^2 + g_2^2}\|_{L^\infty(\mathbb{R}^2)} |v = \text{divg}\right\}. $$

This space is rather large to capture the texture exhibited by the bituminous surfacing images, while its norm is quite easily approximated by a numerical efficient expression. A further justification to use this functional space is its link with the notion of scale presented by Strong et al. in [43]: if $v \in G(\mathbb{R}^2)$, then $\|v\|_G = \sup_{E \subset \Omega} \frac{|\int_E v|}{P(E, \Omega)}$ and scale($\mathbf{x}$) = $\frac{|E|}{|\partial E|}$ for $\mathbf{x} \in E \subset \Omega$, with $P(E, \Omega)$ denoting the perimeter of $E$ in $\Omega$. Therefore, the stronger the penalization of $\|v\|_G$ is, the smaller the scale of the features kept in $v$ is. Based on the following observation: by modelling the crack as a rectangle of $k \times 1$ pixels with $k \gg 1$, its scale behaves like $1/(2n)$ for an image of size $n \times n$ while small features of a pixel size have a scale behaving like $1/(4n)$, one can discriminate these small details from the crack by choosing adequately the penalization of the $G$ norm. To approximate this norm, we rely on the work of Osher et al. [34] and introduce an auxiliary variable naturally stemming from Helmholtz-Hodge decomposition as follows: $g = \nabla Q + P$, with $P$ a divergence-free vector regarded afterwards. An $L^2$ penalization term ensures the closeness between $g$ and $\nabla Q$ and the minimization of the $L^\infty$ norm is applied to $\nabla Q$, yielding a problem related to the absolutely minimizing Lipschitz extension and to the infinity Laplacian.

Equipped with this material, we propose combining the decomposition and Laplacian-based extraction tasks into a single variational framework as follows:

$$ \begin{aligned} \inf \bar{\mathcal{F}}_{\varepsilon}(u, g, Q, v_1, v_2) = & \|f - u - \text{div} g\|_{L^2(\Omega)}^2 \\ & + \mu \|\nabla Q\|_{L^\infty(\Omega)} + \frac{\gamma}{2} \|g - \nabla Q\|_{L^2(\Omega)}^2 \\ & + \rho \int_{\Omega} (v_2^2 + \kappa_{\varepsilon}) (|D^{(2,0)}u|^2 + |D^{(0,2)}u|^2) dx \quad (1) \\ & + \xi_{\varepsilon} \int_{\Omega} (v_1^2 + \zeta_{\varepsilon}) |\nabla u|^2 dx + (\alpha - \beta) \mathcal{G}_{\varepsilon}(v_1) + \beta \mathcal{G}_{\varepsilon}(v_2), \end{aligned} $$

with $\mathcal{G}_{\varepsilon}(l) = \int_{\Omega} \left[ \varepsilon |\nabla l|^2 + \frac{(l-1)^2}{4\varepsilon} \right] dx$, $\kappa_{\varepsilon}, \xi_{\varepsilon}, \zeta_{\varepsilon}$ suitable infinitesimals, and $\varepsilon$ a given parameter. The first three terms cope with the decomposition of $f$ into its texture component $v$ such that $v = \text{divg}$, $g$ being close to $\nabla Q$, and a simplified version of it, $u$, including the crack, while the useless small oscillatory patterns have been removed. The remaining terms are devoted to the crack detection process on $u$. $v_2$ is close to 1 almost everywhere yielding a regularization of $u$, except where $|D^{(2,0)}u|^2 + |D^{(0,2)}u|^2$ is large, $v_2$ being close to 0 in this case. This term thus encodes the fine structures we aim to recover, that is to say the set of cracks is now modelled as $\{v_2 = 0\}$. Similarly, $v_1$ is close to 1 almost everywhere, yielding a regularization of $u$ except in a neighborhood of the jump set of $u$ where $|\nabla u|^2$ is large, $v_1$ being close to 0 in this case to preserve it. It encodes the discontinuity set of $u$, that is to say, the edges with jumps by $\{v_1 = 0\}$. The regularization on $v_1$ and $v_2$ ensures their smoothness and that they are close to 1 almost everywhere, while the terms combining $v_2^2$ and $|D^{(2,0)}u|^2 + |D^{(0,2)}u|^2$ (resp. $v_1^2$ and $|\nabla u|^2$) assure that they are close to 0 on thin structures (resp. on the jump sets). The fidelity term guarantees similarity between $u$ and the observation $f$, while the component $Q$ regulates the size of the features captured by divg. The minimizers are searched in the respective following functional spaces:

- $g$ in $H(\text{div})$.

- $Q$ in $\{Q \in W^{1,\infty}(\Omega) | \int_{\Omega} Q dx = 0\}$ —we are interested in $\nabla Q$ rather than in $Q$ itself, so this condition is not restrictive —.

- $v_1$ in $W^{1,2}(\Omega, [0, 1]) = \{f \in L^2(\Omega, [0, 1]), |\forall \alpha$ such that $|\alpha| \le 1, D^\alpha f \in L^2(\Omega)\}$.

- $u$ in $W_0^{1,2}(\Omega) \cap W^{2,2}(\Omega)$ enabling $u$ to be extended by 0 outside the domain $\Omega$ without being too restrictive.

- $v_2$ in $W^{1,2}(\Omega, [0, 1])$.

The reasons for such requirements will be made clearer in the following. Nevertheless, these assumptions are reasonable and not restrictive if we assume for instance that the observed image $f$ is with compact support. Existence of minimizers is guaranteed as stated in [18, Theorem 11].

## 2.3 Towards a nonlocal related model

Inspired by the strength and robustness of nonlocal methods exemplified in many image processing tasks in [28] (color image denoising, color image deblurring in the presence of Gaussian or impulse noise, color image
inpainting, color image super-resolution, and color fil-
ter demosaicing), and also motivated by their mathe-
matical aspects, we introduce a nonlocal version of our
functional. Based on prior related works by Bourgain,
Brezis and Mironescu [11] introducing a new nonlocal
characterization of the Sobolev spaces W¹, p, 1 < p <
+∞, Aubert and Kornprobst [5] questioning the uti-
lity of this characterization to actually solve variational
problems, Dàvila [17] and Ponce [35] dedicated to the
design of nonlocal counterparts of first order Sobolev
and BV seminorms, we introduce a sequence of radial
mollifiers {ρn}n∈N satisfying: ∀n ∈ N, ∀x ∈ R, ρn(x) =
ρn(|x|); ∀n ∈ N, ρn ≥ 0; ∀n ∈ N, ∫R ρn(x) dx = 1;
∀δ > 0, lim n→+∞ ∫+∞ δ ρn(r) dr = 0, and extend the pre-
vious results to second order operators. Our nonlocal
extension is also deeply inspired by the formulation of
a nonlocal Hessian and the novel characterization of
higher Sobolev and BV-spaces introduced in [30]. This
new operator is derivative free and only requires the
considered function u to be in an Lp space as in our
model. The main difference relies in our choice to treat
independently the directional derivatives and to neglect
the cross derivatives based on the empirical observa-
tion that the crack is most likely horizontal, making it
not too restrictive to consider only the canonical direc-
tions. Besides, thanks to the theory of tempered dis-
tributions, one can show that if u, ∂²u/∂x₁², ∂²u/∂x₂² ∈ L²(R²),
then u ∈ W²,²(R²). We thus introduce a sequence of
functionals Fε,n depending on n and such that the com-
ponent ∫Ω (v²₂ + κε) (|D(²,⁰)u|² + |D(⁰,²)u|²) dx is ap-
proximated by an integral operator involving a differ-
ential quotient and the radial mollifiers. The resulting
model inherits fine analytical properties, has the ad-
vantage of being numerically more tractable than [30]
and of being straightforwardly connected to our imag-
ing problem. We would like to point out that the qual-
ifying term “nonlocal” might sound inadequate in the
sense that the “nonlocal counterpart” includes a param-
eter n, via the mollifier ρn, that is destined to tend to
infinity, concentrating the measure around the point of
interest and removing in some way the nonlocal nature
of the component. Nevertheless, it takes on mathema-
tical interests since, to the best of our knowledge, this
kind of approximation with an independent treatment
of the directional derivatives for second order operators
has not been investigated yet except in [18] where it
has been initiated. It also represents a good compro-
mise, shown by the new Γ-convergence result, between
our local model and the actual numerical model we im-
plement which is introduced and motivated later, and
falls within the “true nonlocal algorithms”. In practice,
for “true nonlocal methods”, the computations are re-

stricted to a small area around the point of interest
(like the NL-means algorithm [13]) and it makes these
methods quite similar to our model from our point of
view.

We now introduce the related nonlocal problem whose
construction relies on a derivative free nonlocal formu-
lation of the L²-norms ∫R² |D(²,⁰)u|² dx and
∫R² |D(⁰,²)u|² dx respectively, and more precisely on the
following theorem ([18, Theorem 10]).

**Theorem 1** Let $u \in W^{2,2}(\mathbb{R}^2)$. Then

$$
\int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{|u(x + he_i) - 2u(x) + u(x - he_i)|^2}{|h|^4} \rho_n(h) dh dx \\
\underset{n \to +\infty}{\longrightarrow} \begin{cases} \|D^{(2,0)}u\|_{L^2(\mathbb{R}^2)}^2 & \text{if } i = 1, \\ \|D^{(0,2)}u\|_{L^2(\mathbb{R}^2)}^2 & \text{if } i = 2. \end{cases}
$$

Equipped with this theoretical result and characteriza-
tion, we are henceforth able to derive a nonlocal coun-
terpart of our local problem (1), phrased in terms of
the functional $\overline{\mathcal{F}}_{\varepsilon,n}$ defined below and for which an ex-
istence result is provided (see [18, Theorem 12]).

**Theorem 2** Let $\Omega$ be a regular bounded open subset of $\mathbb{R}^2$ with boundary of class $C^2$. Let us assume that the functions $t \mapsto \rho_n(t)$, $t \mapsto t^q \rho_n(t)$ are non-increasing for $t \ge 0$ and $q \in [0, 1[$. (Such a function $\rho_n$ exists: for instance, with $q \in ]0, 1[$, $\rho(t) = \frac{e^{-|t|}}{|t|^q}$ and $\rho_n(t) = C n \rho(nt)$ with $C = \frac{1}{\int_{\mathbb{R}} \rho(t) dt}$). With $\kappa_\epsilon, \xi_\epsilon, \zeta_\epsilon > 0$, for any $n \in \mathbb{N}^*$,

$$
\begin{align}
\inf \overline{\mathcal{F}}_{\varepsilon,n}(u, \mathbf{g}, Q, v_1, v_2) &= \|f - u - \operatorname{div} \mathbf{g}\|_{L^2(\Omega)}^2 \nonumber \\
&\quad + \mu \|\nabla Q\|_{L^\infty(\Omega)} + \frac{\gamma}{2} \|\mathbf{g} - \nabla Q\|_{L^2(\Omega)}^2 \nonumber \\
&\quad + \rho \int_{\mathbb{R}^2} (v_{2,e}^2(x) + \kappa_\varepsilon) \tag{2} \\
&\qquad + \sum_{i=1}^2 \int_{\mathbb{R}} \frac{|u_e(x + he_i) - 2u_e(x) + u_e(x - he_i)|^2}{|h|^4} \rho_n(h) dh dx \nonumber \\
&\qquad + \xi_\varepsilon \int_{\Omega} (v_1^2 + \zeta_\varepsilon) |\nabla u|^2 dx + (\alpha - \beta) \mathcal{G}_\varepsilon(v_1) + \beta \mathcal{G}_\varepsilon(v_2), \nonumber
\end{align}
$$

where $v_{2,e}$ and $u_e$ are respectively the extensions of $v_2$
according to [12, Theorem IX.7, p. 158] —by construc-
tion, $0 \le v_{2,e} \le 1$ a.e. — and of $u$ on $\mathbb{R}^2$ by 0 (—with
the regularity assumed on $\Omega$, $v_{2,e}$ and $u_e$ are in $W^{1,2}(\mathbb{R}^2)$
—), admits minimizers $(\overline{u}_n = \overline{u}_{\varepsilon,n}, \overline{g_n} = \overline{g}_{\varepsilon,n}, \overline{Q}_n =$
$\overline{Q}_{\varepsilon,n}, \overline{v}_{1,n} = \overline{v}_{1,\varepsilon,n}, \overline{v}_{2,n} = \overline{v}_{2,\varepsilon,n})$ on $W_0^{1,2}(\Omega) \cap W^{s,2}(\Omega) \times$
$H(div) \times \{Q \in W^{1,\infty}(\Omega) | \int_{\Omega} Q dx = 0\} \times W^{1,2}(\Omega, [0, 1])$
$\times W^{1,2}(\Omega, [0, 1])$, with $s \in [\frac{3}{2}, 2[$.
## 2.4 Main theoretical result

Given the above material, we can now state the main theoretical result, a Γ-convergence one, relating the non-local minimization problem to the local one as n tends to +∞.

**Theorem 3 (Γ-convergence)**

Let $(u_n, \mathbf{g}_n, Q_n, v_{1,n}, v_{2,n}) \in H_0^1(\Omega) \cap W^{s,2}(\Omega) \times H(\mathrm{div}) \times \{Q \in W^{1,\infty}(\Omega) | \int_{\Omega} Q dx = 0\} \times W^{1,2}(\Omega, [0,1]) \times W^{1,2}(\Omega, [0,1])$ with $s \in [\frac{3}{2}, 2[$ be a sequence of minimizers of (2) for each $n \in \mathbb{N}^*$. Let us assume additionally for technical purposes, that

$v_{2,n} \in \{v_2 \in W^{1,2}(\Omega, [0,1])\} \cap W^{1,\infty}(\Omega)$ with

$\sup_{n \in \mathbb{N}^*} \|\nabla v_{2,n}\|_{L^\infty(\Omega)} \le C_1 < \infty$. Then there exist a sub-sequence still denoted by $(u_n, \mathbf{g}_n, Q_n, v_{1,n}, v_{2,n})$ and a minimizer $(\bar{u}, \bar{\mathbf{g}}, \bar{Q}, \bar{v}_1, \bar{v}_2) \in W^{2,2}(\Omega) \cap H_0^1(\Omega) \times H(\mathrm{div})$

$\times \{Q \in W^{1,\infty}(\Omega) | \int_{\Omega} Q dx = 0\} \times W^{1,2}(\Omega, [0,1]) \times W^{1,2}(\Omega, [0,1])$ of (1) such that $u_n \underset{n \to +\infty}{\to} \bar{u}$ in $W^{\frac{3}{2},2}(\Omega)$,

$\mathbf{g}_n \underset{n \to +\infty}{\to} \bar{\mathbf{g}}$ in $H(\mathrm{div})$, $Q_n \underset{n \to +\infty}{\overset{*}{\to}} \bar{Q}$ in $W^{1,\infty}(\Omega)$, $v_{1,n} \underset{n \to +\infty}{\to} \bar{v}_1$ in $W^{1,2}(\Omega)$, $v_{2,n} \underset{n \to +\infty}{\to} \bar{v}_2$ in $W^{1,2}(\Omega)$ and $\bar{F}_{n,\varepsilon}(u_n, \mathbf{g}_n, Q_n, v_{1,n}, v_{2,n}) \underset{n \to +\infty}{\to} \bar{F}_{\varepsilon}(\bar{u}, \bar{\mathbf{g}}, \bar{Q}, \bar{v}_1, \bar{v}_2)$.

*Proof* The proof is subdivided into three main steps that we sketch hereafter. The details can be found in Appendix A.

Step 1. We have proved in what precedes that for any $n \in \mathbb{N}^*$, there exists a solution to (2). We consider a sequence of such minimizers $(u_n, \mathbf{g}_n, Q_n, v_{1,n}, v_{2,n}) \in H^s(\Omega) \cap H_0^1(\Omega) \times H(\mathrm{div}) \times \{Q \in W^{1,\infty}(\Omega) | \int_{\Omega} Q dx = 0\} \times W^{1,2}(\Omega, [0,1]) \times W^{1,2}(\Omega, [0,1])$, with $s \in [\frac{3}{2}, 2[$, and prove successively that up to a subsequence, $(Q_n)$ weakly-* converges to $\bar{Q}$ in $W^{1,\infty}(\Omega)$ with $\int_{\Omega} \bar{Q} dx = 0$, $(v_{1,n})$ weakly converges to $\bar{v}_1$ in $W^{1,2}(\Omega)$ with $\bar{v}_1 \in W^{1,2}(\Omega, [0,1])$, $(v_{2,n})$ weakly converges to $\bar{v}_2$ in $W^{1,2}(\Omega)$ with $\bar{v}_2 \in W^{1,2}(\Omega, [0,1])$, $(\mathbf{g}_n)$ weakly converges to $\bar{\mathbf{g}}$ in $H(\mathrm{div})$, $(u_n)$ weakly converges in $W^{1,2}(\Omega)$ to $\bar{u}$ with $\bar{u} \in W_0^{1,2}(\Omega)$, and $(u_{n,e})$ weakly converges in $H^{\frac{3}{2}}(\mathbb{R}^2)$ to $\bar{u}_e$. Some additional regularity is then obtained showing that $\bar{u} \in W^{2,2}(\Omega) \cap W_0^{1,2}(\Omega)$.

Step 2. In a second stage, we demonstrate that

$$ \limsup_{n \to +\infty} \bar{F}_{n,\varepsilon}(u_n, \mathbf{g}_n, Q_n, v_{1,n}, v_{2,n}) \\ \le \bar{F}_{\varepsilon}(\bar{u}, \bar{\mathbf{g}}, \bar{Q}, \bar{v}_1, \bar{v}_2). $$

Step 3. The last step consists in showing that

$$ \bar{F}_{\varepsilon}(\bar{u}, \bar{\mathbf{g}}, \bar{Q}, \bar{v}_1, \bar{v}_2) \\ \le \liminf_{n \to +\infty} \bar{F}_{n,\varepsilon}(u_n, \mathbf{g}_n, Q_n, v_{1,n}, v_{2,n}), \\ (\bar{u}, \bar{\mathbf{g}}, \bar{Q}, \bar{v}_1, \bar{v}_2) \text{ being a minimizer of (1).} $$

# 3 Numerical implementation

## 3.1 Optimality conditions

We first derive the Euler-Lagrange equations according to each unknown (—for the sake of simplicity, we did not include the interpolation constraints—), with $x = (x_1, x_2)$ and $\mathbf{n} = (n_{x_1}, n_{x_2})$, the unit outward normal to the boundary. Making use of the absolutely minimizing Lipschitz extensions ([3]) for the equation in $Q$, we get:

$$ \left\{ \begin{array}{l} v_1 = \frac{\alpha - \beta}{2\varepsilon} + 2(\alpha - \beta)\varepsilon\Delta v_1, \\ v_2 = \frac{\alpha - \beta}{2\varepsilon} + \frac{2\beta\varepsilon\Delta v_2}{2\rho \sum_{i=1}^2 \int_{\mathbb{R}} \frac{|u_e(x+he_i)-2u_e(x)+u_e(x-he_i)|^2}{|h|^4} \rho_n(h) dh + \frac{\beta}{2\varepsilon}}, \\ g_1 = \partial_{x_1} Q - \frac{2}{\gamma} \partial_{x_1} (f - u - \mathrm{div} \, \mathbf{g}), \\ g_2 = \partial_{x_2} Q - \frac{2}{\gamma} \partial_{x_2} (f - u - \mathrm{div} \, \mathbf{g}), \\ -\mu\Delta_{\infty} Q - \gamma\Delta Q + \gamma\mathrm{div} \, \mathbf{g} = 0, \end{array} \right. $$

combined with suitable boundary conditions. Focusing now on the variable $u$, denoting by $\eta$ a test function, by $\epsilon$ a real number, and setting

$$ J(\epsilon) = \|f - u - \epsilon\eta - \mathrm{div}\,\mathbf{g}\|_{L^2(\Omega)}^2 + \rho \int_{\mathbb{R}^2} (v_{2,e}^2(x) + \kappa_\epsilon) \\ \sum_{i=1}^2 \int_{\mathbb{R}} \frac{|u_e(x+he_i)+\epsilon\eta(x+he_i)-2u_e(x)-2\epsilon\eta(x)+u_e(x-he_i)+\epsilon\eta(x-he_i)|^2}{|h|^4} \\ \rho_n(h) dh \, dx + \xi_\epsilon \int_{\Omega} (v_1^2 + \zeta_\epsilon) |\nabla u + \epsilon\nabla\eta|^2 \, dx, $$

we obtain, with suitable boundary conditions:

$$ \begin{aligned} J'(\epsilon) = & -2 \int_{\Omega} (f-u-\epsilon\eta-\mathrm{div}\mathbf{g})\eta \,dx + \rho \int_{\mathbb{R}^2} (v_{2,e}^2(x) + \kappa_\epsilon) \\ & \sum_{i=1}^2 \int_{\mathbb{R}} \frac{2(u_e(x+he_i)+\epsilon\eta(x+he_i)-2u_e(x)-2\epsilon\eta(x)+u_e(x-he_i)+\epsilon\eta(x-he_i))}{|h|^4} \\ & (\eta(x+he_i)-2\eta(x)+\eta(x-he_i))\rho_n(h) \,dh \,dx \\ & + 2\xi_\epsilon \int_{\Omega} (v_1^2 + \zeta_\epsilon) \langle \nabla u + \epsilon\nabla\eta, \nabla\eta \rangle \,dx, \\ J'(0) = & - \int_{\Omega} 2(f-u-\mathrm{div}\mathbf{g})\eta \,dx + \rho \int_{\mathbb{R}^2} (v_{2,e}^2(x) + \kappa_\epsilon) \\ & \sum_{i=1}^2 \int_{\mathbb{R}} \frac{2(u_e(x+he_i)-2u_e(x)+u_e(x-he_i))(\eta(x+he_i)-2\eta(x)+\eta(x-he_i))}{|h|^4} \\ & \rho_n(h) dh \,dx + 2\xi_\epsilon \int_{\Omega} (v_1^2 + \zeta_\epsilon) \langle \nabla u, \nabla\eta \rangle \,dx = 0, \\ J'(0) = & - \int_{\Omega} (f-u-\mathrm{div}\mathbf{g})\eta \,dx - \xi_\epsilon \int_{\Omega} \mathrm{div}((v_1^2+\zeta_\epsilon)\nabla u)\eta \,dx \\ & + \rho \int_{\mathbb{R}^2} \sum_{i=1}^2 \int_{\mathbb{R}} \left( \frac{-2(v_{2,e}^2(x)+\kappa_\epsilon)(u_e(x+he_i)-2u_e(x)+u_e(x-he_i))}{|h|^4} \right. \\ & \left. + \frac{(v_{2,e}^2(x-he_i)+\kappa_\epsilon)(u_e(x)-2u_e(x-he_i)+u_e(x-2he_i))}{|h|^4} \right) \\ & + \frac{(v_{2,e}^2(x+he_i)+\kappa_\epsilon)(u_e(x+2he_i)-2u_e(x+he_i)+u_e(x))}{|h|^4} \eta\rho_n(h) \,dh \,dx \\ & = 0, \end{aligned} $$
which yields the desired equation in $u$. The equation in $Q$ encodes the most technical elements (nonlinearity, singularity at $\nabla Q = 0$). The next subsection is thus devoted to the derivation and analysis of the evolution equation satisfied by $Q$ and involving the infinity Laplacian, in the viscosity solution theory framework. This theory applies to certain partial differential equations and allows merely continuous functions to be solutions of fully nonlinear equations. We refer the reader to [15] for a general introduction including the main material and definitions.

## 3.2 Theoretical results: existence, uniqueness and regularity of the solution of the evolution equation in $Q$

In that purpose, let us introduce an artificial time and embed the equation in $Q$ in a time-dependent setting. Let $T > 0$ be given. The parabolic equation in the unknown $Q$ (which is the one we solve in practice) is thus given by

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial Q}{\partial t} = \mu \Delta_{\infty} Q + \gamma \Delta Q - \gamma \operatorname{div} g \quad \text{on } \mathbb{R}^2 \times (0, T), \\ Q(x, 0) = Q_0(x) \quad \text{on } \mathbb{R}^2, \end{array} \right. \qquad (EE) $$

with $Q_0 \in W^{1,\infty}(\mathbb{R}^2)$ and $B_0$ its Lipschitz constant. (To remove the problem of boundary conditions, we work on $\mathbb{R}^2$ for the spatial domain). We now give a thorough analysis of the PDE in $Q$ in the viscosity solution theory framework. Let us first state the additional following assumption:

$\operatorname{div} g$ is bounded and is Lipschitz continuous uniformly in time with $\kappa_g$ its Lipschitz constant independent of time.

For the sake of conciseness and using the normalized version of the infinity Laplacian, the evolution equation is now written in the form

$$ \frac{\partial Q}{\partial t} + G(x, t, \nabla Q, \nabla^2 Q) = 0, $$

with $G: \mathbb{R}^2 \times [0, T) \times \mathbb{R}^2 \times S^2$ ($S^2$ being the set of symmetric $2 \times 2$ matrices equipped with its Loewner partial order) defined by

$$ \begin{align*} G(x, t, \mathbf{p}, X) &= -\gamma \operatorname{trace}(X) - \mu \left\langle \frac{\mathbf{p}}{|\mathbf{p}|}, X \frac{\mathbf{p}}{|\mathbf{p}|} \right\rangle + \gamma \operatorname{div} g \\ &= -\gamma \operatorname{trace}(X) - \mu \operatorname{trace} \left( \frac{\mathbf{p} \otimes \mathbf{p}}{|\mathbf{p}|^2} X \right) + \gamma \operatorname{div} g, \\ &= E(X) + F(\mathbf{p}, X) + \gamma \operatorname{div} g, \end{align*} $$

⊗ representing the Kronecker product, and with the following properties:

1. The operators $G, E: X \mapsto -\gamma \operatorname{trace}(X)$ and $F: (\mathbf{p}, X) \mapsto -\mu \operatorname{trace} \left( \frac{\mathbf{p} \otimes \mathbf{p}}{|\mathbf{p}|^2} X \right)$ are independent of $Q$ and are elliptic, i.e., $\forall X, Y \in S^2, \forall \mathbf{p} \in \mathbb{R}^2 \setminus \{0\}_{\mathbb{R}^2}$, if $X \preceq Y$ then $F(\mathbf{p}, X) \ge F(\mathbf{p}, Y)$ since $F(\mathbf{p}, X) - F(\mathbf{p}, Y) = -\mu \operatorname{trace} \left( \frac{\mathbf{p} \otimes \mathbf{p}}{|\mathbf{p}|^2} X \right) + \mu \operatorname{trace} \left( \frac{\mathbf{p} \otimes \mathbf{p}}{|\mathbf{p}|^2} Y \right) = -\mu \operatorname{trace} \left( \frac{\mathbf{p} \otimes \mathbf{p}}{|\mathbf{p}|^2} (X - Y) \right) = -\mu \left\langle \frac{\mathbf{p}}{|\mathbf{p}|}, (X - Y) \frac{\mathbf{p}}{|\mathbf{p}|} \right\rangle \ge 0$ as $X \preceq Y$.

2. $F$ is locally bounded on $\mathbb{R}^2 \times S^2$, continuous on $\mathbb{R}^2 \setminus \{0\}_{\mathbb{R}^2} \times S^2$, and $F^*(0, 0) = F_*(0, 0) = 0$, where $F^*(\text{resp. } F_*)$ is the upper semicontinuous (usc) envelope (resp. lower semicontinuous (lsc) envelope) of $F$. Indeed, using Rayleigh quotient and its properties, it is not difficult to see that for nonzero vector $\mathbf{p}$, $\lambda_{min}(X) \le \operatorname{trace} \left( \frac{\mathbf{p} \otimes \mathbf{p}}{|\mathbf{p}|^2} X \right) = \frac{\langle \mathbf{p}, X \mathbf{p} \rangle}{\langle \mathbf{p}, \mathbf{p} \rangle} = R(X, \mathbf{p}) \le \lambda_{max}(X)$, $\lambda_{min}$ (resp. $\lambda_{max}$) denoting the smallest (resp. biggest) eigenvalue of $X$.

We begin our analysis by establishing a comparison principle stating that if a sub-solution and a super-solution are ordered at the initial time, then they are ordered at any time.

**Theorem 4 (Comparison principle, adapted from [23])** Assume $(H)$ and let $q^{-} : \mathbb{R}^2 \times [0, T) \to \mathbb{R}$ be a bounded upper semicontinuous sub-solution and $q^{+} : \mathbb{R}^2 \times [0, T) \to \mathbb{R}$ be a bounded lower semicontinuous super-solution of (EE). Assume that $q^{-}(x, 0) \le Q_0(x) \le q^{+}(x, 0)$ in $\mathbb{R}^2$, then $q^{-} \le q^{+}$ in $\mathbb{R}^2 \times [0, T)$.

*Proof* We refer the reader to Appendix B.

We now construct barriers in order to use Perron's method.

**Theorem 5 (Construction of barriers, adapted from [23])** Assume $(H)$ and let $Q_0 \in W^{1,\infty}(\mathbb{R}^2)$. Then $q^+ = \sup_{x \in \mathbb{R}^2} Q_0(x) + \gamma \| \operatorname{div} g \|_{L^\infty(\mathbb{R}^2 \times [0, T))} t$ and $q^- = \inf_{x \in \mathbb{R}^2} Q_0(x) - \gamma \| \operatorname{div} g \|_{L^\infty(\mathbb{R}^2 \times [0, T))} t$ are respectively super- and sub-solution of (EE).

*Proof* The proof of this theorem is straightforward. Indeed, since $Q_0 \in W^{1,\infty}(\mathbb{R}^2)$, then $q^-$ and $q^+$ are twice differentiable in space and once differentiable in time and are bounded. Besides $q^-(x, 0) = \inf_{x \in \mathbb{R}^2} Q_0(x) \le Q_0(x)$, $q^+(x, 0) = \sup_{x \in \mathbb{R}^2} Q_0(x) \ge Q_0(x)$, $\forall x \in \mathbb{R}^2$. Furthermore, $\forall (x, t) \in \mathbb{R}^2 \times [0, T)$, $q_t^-(x, t) = -\gamma \| \operatorname{div} g \|_{L^\infty(\mathbb{R}^2 \times [0, T))}$, $Dq^-(x, t) = 0$ and $D^2q^-(x, t) = 0$;
$q_t^+(x,t) = \gamma \|\text{div } \mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0,T))}$, $Dq^+(x,t) = 0$ and $D^2q^+(x,t) = 0$.

It follows then that $\forall (x,t) \in \mathbb{R}^2 \times [0,T)$,

$$q_t^-(x,t) + G_*(t,x, q^-(x,t), Dq^-(x,t), D^2q^-(x,t)) = - \gamma \|\text{div } \mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0,T))} + \gamma \text{divg}(x,t) \le 0,$$

and

$$q_t^+(x,t) + G^*(t,x, q^+(x,t), Dq^+(x,t), D^2q^+(x,t)) = \gamma \|\text{div } \mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0,T))} + \gamma \text{divg}(x,t) \ge 0.$$

Hence $q^-$ is a sub-solution and $q^+$ a super-solution of (EE).

We are now able to use Perron's method and to give the following existence result.

**Theorem 6 (Existence and uniqueness of a solution)** Assume (H) and $Q_0 \in W^{1,\infty}(\mathbb{R}^2)$. Then there exists a unique bounded continuous solution of (EE) in $\mathbb{R}^2 \times [0,T)$.

*Proof* We follow Perron's method ([15, Theorem 4.1]). Indeed, we have constructed a bounded sub-solution and a bounded super-solution which fall within the comparison principle.

Then function $Q = \sup\{w, w$ sub-solution of (EE) and $\forall (x,t) \in \mathbb{R}^2 \times [0,T)$, $q^-(x,t) \le w(x,t) \le q^+(x,t)\}$ is a solution of (EE) potentially discontinuous. Clearly, $Q$ is bounded since $q^+$ and $q^-$ are bounded and since $Q$ is a solution of (EE), then $Q_*$ is a super-solution and $Q^*$ is a sub-solution. Applying the comparison principle yields $Q^* \le Q_*$. But by definition, $Q_* \le Q^*$, so we finally get $Q = Q^* = Q_*$ meaning that $Q$ is continuous on $\mathbb{R}^2 \times [0,T)$.

We now analyze the regularity of the solution. Let us first consider the regularity in space.

**Theorem 7 (Lipschitz regularity in space, adapted from [23])** Assume (H) and that $\|\nabla Q_0\|_{L^\infty(\mathbb{R}^2)} \le B_0$ with $B_0 > 0$. Then the solution of (EE) is Lipschitz continuous in space and satisfies $\|\nabla Q(\cdot, t)\|_{L^\infty(\mathbb{R}^2)} \le B(t)$, with $B(t) = \gamma \kappa_g t + B_0$.

*Proof* We refer the reader to Appendix C.

Besides, by analyzing the smoothness in time of the solution, we can show that this solution is uniformly continuous in time.

**Theorem 8 (adapted from [23])** Assume (H), and that divg is uniformly continuous in time with $\omega_{\text{divg}}$ its modulus of continuity. Then the solution of (EE) is uniformly continuous in time.

*Proof* We refer the reader to Appendix D.

This concludes the analysis of the evolution equation in $Q$ in the viscosity solution theory setting.

### 3.3 Handling of the nonlocal component and resulting algorithm:

Some limitations can be identified regarding the theoretical expression of the nonlocal component, which leads us to introduce a numerical compromise. In particular, the properties of the kernels $\rho_n$ are not all well-adapted for imaging problems. We recall that they satisfy $\forall n \in \mathbb{N}, \forall x \in \mathbb{R}, \rho_n(x) = \rho_n(|x|); \forall n \in \mathbb{N}, \rho_n \ge 0; \forall n \in \mathbb{N}, \int_{\mathbb{R}} \rho_n(x) dx = 1; \forall \delta > 0, \lim_{n \to +\infty} \int_{\delta}^{+\infty} \rho_n(r) dr = 0$. While the last three properties are consistent with the discrete setting as in practice the weights are non-negative, normalized and concentrated near the current window center, the radiality property is not relevant for imaging problems. In fact, we would like to define a nonlocal version of our second order operator at a given point $x$ that grants more weight to points that belong to the same region as $x$ not only based on the difference between intensities but for instance, on the difference between patches around those points and thus not favoring spatial proximity. This first observation led us to reconsider the definition of our nonlocal weights. We would like also to stress that, to the best of our knowledge, this assumption of radiality cannot be removed in the theoretical analysis (see [11], [17] and [35]). We thus redefine the weights $\frac{\rho_n(h)}{|h|^4}$ involved in the theoretical expression, inspired by the NL-means algorithm ([13]). Integrating additional information related to the image content $f$ is pertinent here since it strengthens/favors similar behaviors. More precisely, we put more weight to neighboring pixels that have similar edges/creases identified by the gradient / Hessian behaviour as seen in Lemma 1 and that are geographically close. In that purpose, we consider the following nonlocal weights

$$\frac{\rho_n(h)}{|h|^4} \approx w_{f,i,x}(h) = \exp\left(-\frac{d_{f,x,i}(h)}{\alpha^2}\right) \text{ where } d_{f,x,i}(h) = \int_{\mathbb{R}^2} G_a(t) \|\nabla f(x+t) - \nabla f(x+t+he_i)\|^2 dt \text{ or } d_{f,x,i}(h) = \int_{\mathbb{R}^2} G_a(t) \|\nabla^2 f(x+t) - \nabla^2 f(x+t+he_i)\|^2 dt, G_a \text{ being a Gaussian kernel with standard deviation } \alpha \text{ controlling the patch size, and } \alpha \text{ is the filtering parameter. The data structure summed-area table also known as integral image is used to efficiently compute the involved sums. Later in the text is a pseudo-code associated with the computation of the nonlocal weights (Algorithm 1). The resulting nonlocal algorithm relies on an alternating strategy in which we solve the Euler-Lagrange equations related to each unknown using specifically, a time-dependent scheme for } Q = Q(x_1, x_2, t) \text{ (nonlinear overrelaxation method, see [14, Section 4])}, \text{ and a stationary semi-implicit fixed-point scheme in } v_1 = v_1(x_1, x_2), v_2 = v_2(x_1, x_2) \text{ and } g = g(x_1, x_2). \text{ Neumann boundary conditions for } u \text{ and } Q, \text{ and Dirichlet boundary conditions for } v_1 - 1, v_2 - 1 \text{ and } g \text{ are applied. We denote by}$$
$\Delta x_1 = \Delta x_2 = h = 1$ the space steps in each direction, by $\Delta t$ the time step, and by $f_{i,j}$, $u_{i,j}^n$, $Q_{i,j}^n$, $v_{1,i,j}^n$, $v_{2,i,j}^n$, $g_{i,j}^n = (g_{1,i,j}^n, g_{2,i,j}^n)^T$ the discrete versions of $f$, $u$, $Q$, $v_1$, $v_2$ and $g$ at iteration $n \ge 0$, for $1 \le i \le M$, $1 \le j \le N$.

**Algorithm 1** Computation of the nonlocal weights inspired by the NL-means algorithm.

**Input:** Initial image $f$
**Output:** Weights in the first direction $w_{f,x,1}$, the shifts of the selected neighbors in the first direction $\text{indice}_1$, weights in the second direction $w_{f,x,2}$, the shifts of the selected neighbors in the second direction $\text{indice}_2$.

1. Define $w := \text{window size}$, $p := \text{patch size}$, $h := 0.25$, NbNeigh := number of actual required neighbors including the closest one.

2. Compute the extended image by symmetry.

**for all pixels $x = (x_1, x_2)$ do**

3. Compute the distances (we do not make use in practice of the Gaussian kernel $G_a$, yielding equal weights for the contributions)

$$d_{f,x,1}(y) = \sum_{i=-\frac{p-1}{2}}^{\frac{p-1}{2}} \sum_{j=-\frac{p-1}{2}}^{\frac{p-1}{2}} \left( \frac{\|\nabla f(x_1+i, x_2+j) - \nabla f(x_1+y+i, x_2+j)\|^2}{2} + \frac{\|\nabla f(x_1+i, x_2+j) - \nabla f(x_1+i-y, x_2+j)\|^2}{2} \right)$$

or

$$d_{f,x,1}(y) = \sum_{i=-\frac{p-1}{2}}^{\frac{p-1}{2}} \sum_{j=-\frac{p-1}{2}}^{\frac{p-1}{2}} \left( \frac{\|\nabla^2 f(x_1+i, x_2+j) - \nabla^2 f(x_1+y+i, x_2+j)\|^2}{2} + \frac{\|\nabla^2 f(x_1+i, x_2+j) - \nabla^2 f(x_1+i-y, x_2+j)\|^2}{2} \right)$$

between all patches centered at $x + ye_1$ of size $p$ with $0 \le y \le \frac{iw-1}{2}$ and the patch centered at the current $x$. Compute the distance

$$d_{f,x,2}(y) = \sum_{i=-\frac{p-1}{2}}^{\frac{p-1}{2}} \sum_{j=-\frac{p-1}{2}}^{\frac{p-1}{2}} \left( \frac{\|\nabla f(x_1+i, x_2+j) - \nabla f(x_1+i, x_2+y+j)\|^2}{2} + \frac{\|\nabla f(x_1+i, x_2+j) - \nabla f(x_1+i, x_2+j-y)\|^2}{2} \right).$$

or

$$d_{f,x,2}(y) = \sum_{i=-\frac{p-1}{2}}^{\frac{p-1}{2}} \sum_{j=-\frac{p-1}{2}}^{\frac{p-1}{2}} \left( \frac{\|\nabla^2 f(x_1+i, x_2+j) - \nabla^2 f(x_1+i, x_2+y+j)\|^2}{2} + \frac{\|\nabla^2 f(x_1+i, x_2+j) - \nabla^2 f(x_1+i, x_2+j-y)\|^2}{2} \right).$$

between all patches centered at $x + ye_2$ of size $p$ with $0 \le y \le \frac{iw-1}{2}$ and the patch centered at the current $x$.

**end for**
**for all pixels $x$ do**

4.1. Sort the previous distances in ascending order for each direction and keep only the lowest NbNeigh-1 values with the corresponding shift ($y$) for the first direction in $\text{indice}_1$ and for the second direction in $\text{indice}_2$.

4.2. Add the closest neighbor in geographical sense and in each direction to make the weights more similar to the theoretical ones.

**end for**
**for all pixels $x$ do**

5.1. Compute $w_{f,x,1}$ by the following formula: $w_{f,x,1}(y) = 0$ if $y$ does not belong to the previous list of neighbors, $w_{f,x,1}(y) = \exp\left\{-\frac{d_{f,x,1}(y)}{h^2}\right\}$ otherwise. Compute $w_{f,x,2}$ by the following formula: $w_{f,x,2}(y) = 0$ if $y$ does not belong to the previous list of neighbors, $w_{f,x,2}(y) = \exp\left\{-\frac{d_{f,x,2}(y)}{h^2}\right\}$ otherwise.

**end for**
**Return** $w_{f,x,1}$, $w_{f,x,2}$, $\text{indice}_1$, $\text{indice}_2$.

**Algorithm 2** Nonlocal alternating algorithm for crack recovery.

1. [Initialization step]:
   $u^0 = f$, $g^0 = 0$, $Q^0 = 0$, $v_1^0 = 1$ and $v_2^0 = 1$.

2. [Main step]:
   **for** $n=1$ to $n=500$ **do**

**if** $n\%100==0$ **then**

2.1 Compute the nonlocal weights associated to $u^n$:
$[w_{f,x,1}, w_{f,x,2}, \text{indice}_1, \text{indice}_2] = \text{compute\_weights}(u^n)$.

**end if**
2.2 Compute for all pixels $(i, j)$:

$$|\nabla u^n|^2_{i,j} = (u_{i+1,j}^n - u_{i,j}^n)^2 + (u_{i,j+1}^n - u_{i,j}^n)^2,$$

$$v_{1,i,j}^{n+1} = \frac{\frac{(\alpha-\beta)}{2\epsilon} + 2(\alpha-\beta)\epsilon (v_{1,i+1,j}^n + v_{1,i-1,j}^n + v_{1,i,j+1}^n + v_{1,i,j-1}^n - 4v_{1,i,j}^n)}{2\epsilon |\nabla u^n|^2_{i,j} + \frac{(\alpha-\beta)}{2\epsilon}},$$

$$\text{nlnormu}_{i,j}^2 = \sum_{h \in \text{indice}_1} (u^n(i+h, j) - 2u^n(i, j) + u^n(i-h, j))^2 w_{f, i, j, 1}(h) + \sum_{h \in \text{indice}_2} (u^n(i, j+h) - 2u^n(i, j) + u^n(i, j-h))^2 w_{f, i, j, 2}(h),$$

$$v_{2,i,j}^{n+1} = \frac{\frac{\beta}{2\epsilon} + 2\beta\epsilon (v_{2,i+1,j}^n + v_{2,i-1,j}^n + v_{2,i,j+1}^n + v_{2,i,j-1}^n - 4v_{2,i,j}^n)}{2\rho \text{nlnormu}_{i,j}^2 + \frac{\beta}{2\epsilon}},$$

$$\frac{u_{i,j}^{n+1} - u_{i,j}^n}{\Delta t} = (f_{i,j} - u_{i,j}^n - \frac{g_{1,i+1,j}^n - g_{1,i-1,j}^n}{2} - \frac{g_{2,i,j+1}^n - g_{2,i,j-1}^n}{2} + \xi_\epsilon \left[ (v_{1,i,j}^n)^2 (u_{i,j+1}^n - u_{i,j}^n) - (v_{1,i,j-1}^n)^2 (u_{i,j}^n - u_{i,j-1}^n) \right] + \xi_\epsilon \left[ (v_{1,i,j}^n)^2 (u_{i+1,j}^n - u_{i,j}^n) - (v_{1,i-1,j}^n)^2 (u_{i,j}^n - u_{i-1,j}^n) \right] - \sum_{h \in \text{indice}_1} \rho w_{f,i,j,1}(h) \left[ (v_{2,i+h,j}^n)^2 (u_{i+2h,j}^n - 2u_{i+h,j}^n + u_{i,j}^n) - 2(v_{2,i,j}^n)^2 (u_{i+h,j}^n - 2u_{i,j+h}^n + u_{i,j-2h,j}^n) \right] + \xi_\epsilon \left[ (v_{2,i,j}^n)^2 (u_{i,j-2h,j}^n - 2u_{i,j+h}^n + u_{i,j+2h}^n) + (v_{2,i,j+h}^n)^2 (u_{i,j+2h}^n - 2u_{i,j+h}^n + u_{i,j-2h}^n) \right],$$

using symmetry if it does not belong to the image domain, and equations derived in the same way for $g_1^n$, $g_2^n$ and $Q^n$.
**end for**
*Remark 1* Provided that we solve exactly each subproblem, a result of convergence of the alternating strategy to existing partial minimisers can be proved.

In order to improve the computation efficiency, we propose an MPI parallelization of our code, which motivates our choice of a rather simple alternating minimization method for which a decomposition domain approach is well-suited.

## 3.4 MPI parallelization

The parallelization of the C code is motivated by the natural geometry of the problem—an image is defined on a rectangle domain $\Omega$—making the partition of the image domain into subdomains supporting simultaneous local computations relevant. Note also that the computational complexity increases with the image size (in practice we have worked with some images of size $2248 \times 4000$), requiring more memory to store the data and the results, this fact being particularly marked in the nonlocal model that involves the resolution of a nonlocal partial differential equation.

The meshing is made of $ntx$ interior points in the row direction (we removed the first and last rows of points) and $nty$ interior points in the column direction (we removed the leftmost and rightmost columns of points). The implementation revolves around the following steps:

(i) we generate a Cartesian topology (see Figure 1 for an example), each subdomain comprising $w-1$ ($w$ is the window size $-$) rows of ghost cells above, $w-1$ rows of ghost cells below, and similarly for the columns, in order to store the data exchanged with neighboring subdomains. Either the developer selects the number of nodes in each direction, or it is left to the MPI library. Some latitude is also given to the user in terms of periodicity ($-$ a periodicity can be applied on the grid in each direction if required thanks to the array **periods** $-$) and reorganization ($-$ if the user wants the processes to keep the same rank as in the original communicator $-$).

(ii) For each subdomain, we recover the bounds with respect to the original image reference frame of the indices $i$ and $j$ that are then stored in the 1d array **tabBounds**: $tabBounds[0]=sx$, $tabBounds[1]=ex$, $tabBounds[2]=sy$ and $tabBounds[3]=ey$. The numbering of the original image reference frame starts at 0 and the origin is the top left corner. The indices $(sx-w+1, ..., sx-1, ex+1, ..., ex+w-1, sy-w+1, ..., sy-1, ey+1, ..., ey+w-1)$

<figure><img src="image_2.png" /><figcaption>Fig. 1: Cartesian topology in the local case with 4 processes on a (10 × 10) image before reversing the row axis:</figcaption></figure>

Process 0: $(sx = 2, ex = 4, sy = 5, ey = 7)$;

Process 1: $(sx = 5, ex = 7, sy = 5, ey = 7)$;

Process 2: $(sx = 2, ex = 4, sy = 2, ey = 4)$;

Process 3: $(sx = 5, ex = 7, sy = 5, ey = 7)$.

are used to store the data sent by the 8-connected neighboring subdomains.

*Created function: void domaine(MPI.Comm comm2d, int rang, int * coordinates, int ntx, int nty, int * dimensions, int * tabBounds, int * periods, int reorganization, int nb_procs)*

(iii) For each subdomain, the neighboring subdomain ranks are returned. This is achieved thanks to the routine **voisinage** and these ranks are stored respectively in the 1d array voisin (for the 4-connected blocks) and voisin_diagonale (for the diagonally connected blocks). The necessity of storing the ranks of diagonally connected blocks arises from the numerical schemes used to discretize the partial differential equations satisfied by $u, Q, g_1$ and $g_2$ that involve for instance components like $Q_{i-1,j+1}$.

*Created function: void voisinage(MPI.Comm comm2d, int * voisin, int * voisin_diagonale, int * coordinates, int * dimensions)*

(iv) Once the Cartesian topology is created, one needs to distribute the data file (image data) to each subdomain in parallel (more precisely, the portion of the data file that must be visible for the related process). The general file manipulation function **MPI.File_open** is called. We then create a dataset type **MPI.Datatype mysubarray** describing a two-dimensional subarray (the portion of the image related to the current subdomain) of a bigger two-dimensional array (the image here) (routine **MPI_Type_create_subarray**). The **MPI_File_set_view** routine allows to change the process view of the data in the file: the beginning of the data accessible in the file through that view is set to 0, the type of data is set to MPI_DOUBLE, and the
distribution of data to processes is set to mysubar-ray. Then the **MPI_File_read** routine enables us to read the file starting at the specified location.

(v) Derived datatypes are created, describing the rows, columns and 2×2 diagonal arrays involved in the MPI communications: **MPI_Datatype_type**(w-1)colonnes, type_ligne, type_colonne_mono, type**(w-1)**lignes among others. In that purpose, the routines **MPI_Type_contiguous** —creating a contiguous datatype, here a single row of data—and **MPI_Type_vector** —general constructor that allows replication of a datatype into locations that consist of equally spaced blocks, here a single column, a group of (w-1) adjacent columns and a group of (w-1) adjacent rows—are used. For the communications with diagonally connected subdomains, a datatype describing a two-dimensional subarray of size 2×2 of a bigger two-dimensional array is created for each spatial configuration: top left corner, top right corner, bottom left corner and bottom right corner.

(vi) The communications are then handled with the routine **MPI_Sendrecv**. This send-receive operation combines in one call the sending of a message to one destination and the receiving of another message from another process. As an illustration (see also Figure 2), **u_local)_mat** denoting the local array describing **u** :

(a) **u_local mat**( *sx* : *sx* + *w* − 2, :) is sent to the northern neighbor that receives data in **u_local mat**( *ex* + 1 : *ex* + *w* − 1, ;);

(b) **u_local mat**( *ex* − *w* + 2 : *ex*, :) is sent to the southern neighbor that receives data in **u_local mat**( *sx* − *w* + 1 : *sx* − 1, ;);

(c) **u_local mat**(, *sy* : *sy* + *w* − 2) is sent to the western neighbor that receives data in **u_local mat**(, *ey* + 1 : *ey* + *w* − 1);

(d) **u_local mat**(, *ey* − *w* + 2 : *ey*) is sent to the eastern neighbor that receives data in **u_local mat**(, *sy* − *w* + 1 : *sy* − 1);

(e) **u_local mat**( *sx* : *sx* + 1, *ey* − 1 : *ey*) is sent to the northeast neighbor that receives data in **u_local mat**( *ex* + 1 : *ex* + 2, *sy* − 2 : *sy* − 1);

(f) **u_local mat**( *ex* − 1 : *ex*, *ey* − 1 : *ey*) is sent to the southeast neighbor that receives data in **u_local mat**( *sx* − 2 : *sx* − 1, *sy* − 2 : *sy* − 1);

(g) **u_local mat**( *ex* − 1 : *ex*, *sy* : *sy* + 1) is sent to the southwest neighbor that receives data in **u_local mat**( *sx* − 2 : *sx* − 1, *ey* + 1 : *ey* + 2);

(h) **u_local mat**( *sx* : *sx* + 1, *sy* : *sy* + 1) is sent to the northwest neighbor that receives data in **u_local mat**( *ex* + 1 : *ex* + 2, *ey* + 1 : *ey* + 2).

For subdomains with at least one edge included in
the image domain boundary, diagonal communica-
tions do not occur and are replaced by communica-
tions with the ghost cells of the involved contiguous
subdomain.

Created function : void communication(double
** u, double ** v1_local mat, double **
v2_local mat, double ** g1_local mat, dou-
ble ** g2_local mat, double ** Q_local mat,
int * tab bounds, MPI.Comm comm2d,int *
voisin, int * voisin _diagonale)

<figure><img src="image_3.png" /><figcaption>Fig. 2: Phantom cells and communications, example on pro-
cess 0 for the local case.</figcaption></figure>

(vii) The values of $u^{n+1}$, $Q^{n+1}$, $v_1^{n+1}$, $v_2^{n+1}$, $g_1^{n+1}$ and $g_2^{n+1}$ are computed using the above mentioned finite difference schemes ($v_1$ and $v_2$ have been initialized to 1, $g_1$, $g_2$ and $Q$ to 0, while $u$ has been set to the values of the original image) and the question of boundary conditions is addressed. For the sake of simplicity, we have assumed homogeneous Neumann boundary conditions for $u$ and $Q$, simulating reflection of the array through its boundaries, and Dirichlet boundary conditions for $v_1-1$, $v_2-1$ and $g$. We thus identify subdomains with at least one edge included in the image domain boundary. As an illustration, for the variable $u$, we replicate the row of index $iw-1$ (resp. $ex-sx+iw-1$) of the local matrix in the row of index $iw-2$ (resp. $ex-sx+iw$), similarly for the columns. The diagonal components are
processed using point symmetry. The newly computed values are then transmitted to neighboring subdomains and a new time step is achieved.

(viii) We finally write the result (component $v_2$ encoding the fine structure we aim to recover) in a file using again MPI I/O. The general file manipulation function **MPI_File_open** is called. We then create a new datatype describing the two-dimensional sub-array extracted from the original local array $v_2$ local mat when removing the ghost cells. The routine **MPI_File_set_view** changes the process view of the data in the file, while the routine **MPI_File_write_all** writes the file, starting at the locations specified by individual file pointers.

Created function : **void ecrire_mpi(double * v2_localimox,int ntx,int nty,int * tab� bounds, MPI_Cmm comm2d)**

(ix) The computation of the nonlocal weights is done by every process on its associated subdomain.

The computations have been made with the supercomputer Myria operated by the CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie, http://www.criann.fr/). Myria is an ATOS BULL solution with 11144 computing cores, with a power of 403 TFlops Xeon, 170 TFlops GPU and 27 TFlops Xeon Phi KNL. Myria also has a fast storage space of about 2.5 Po. Submission of the work is done through the SLURM software.

The local algorithm has been applied to a bituminous surfacing image of size 2248 × 4000, requiring 500 time step iterations both in the case of the sequential algorithm, execution time: 289.6 seconds, and parallelized algorithm with 224 tasks, execution time: 1.7 seconds. The efficiency is of 74% from 1 to 224 tasks and of 100% from 28 to 224 tasks. The nonlocal algorithm has been applied to the same bituminous surfacing image of size 2248 × 4000, requiring 360 time step iterations both in the case of the sequential algorithm, execution time: 596.582511 seconds, and parallelized algorithm with 224 tasks, execution time: 4.409619 seconds.

The numerical applications presented below have been obtained in this MPI setting, with the same orders of reduction factor when passing from a single to 224 tasks. We no longer comment on the optimization of the computation code but now focus on the quantitative and qualitative analysis of the results.

## 4 Numerical experiments

Before showing some numerical simulations on real datasets resulting from the application of the above algorithm (with different weights) and discussing its ac-

currency, we would like to give some insight on how to choose the parameters involved in our model. The more sensitive ones are $α, β, ρ$. Smaller values for these parameters induce less regularization and consequently more information (edges/creases) in the $v_1$ and $v_2$ components. Parameter $μ$ weights the $L^∞$-norm of $|∇Q|$ and thus influences the size of the texture captured in $v = \text{divg}$. A higher parameter $μ$ leads to smaller scale features caught by $v = \text{divg}$. $ε$ plays on the thickness of the contours in $v_1$ and $v_2$: larger values of $ε$ imply thicker contours.

We now apply the proposed algorithm to crack detection on concrete walls, both on Figure 3 (size 501 × 501) and Figure 4 (size 285 × 429), courtesy of A. Drogoul. We depict the three main components of the decomposition/segmentation process, i.e., $u$ local, $v = \text{div } g$ local, $v_2$ local, for the local model (related to problem (1) and for the sake of comparison with the nonlocal model), and $u$ nonlocal and $v_2$ nonlocal, for the three versions of our nonlocal algorithm (weights based on gradient stands for the algorithm in which the weights are computed thanks to the distance between image gradients, weights based on Hessian is for the algorithm in which the weights are given by comparing the Hessians of the image, and only 4 closest neighbors refers to our nonlocal implementation in which we consider only the four closest neighbors), as well as the results obtained with Aubert and Drogoul's topological gradient method ([6, 20]). The cracks are correctly enhanced while the oscillatory patterns are well captured by the $v = \text{div } g$ component. We also observe that the crack is smoother with Aubert and Drogoul's method and our method seems to detect more accurately the center of the crack. The contrast is the highest for the nonlocal method with weights based on Hessian information and the connectivity of the crack appears clearer too. Again, the role of the decomposition part of the algorithm is highlighted (Figure 3) by depicting the obtained $v_2$ component when decomposition is turned off in our local implementation (spurious details are visible on the top of the image). Also, the piecewise linear nature of the component $u$ in Figure 4 is properly returned.

We conclude the paper with two applications dedicated to crack detection on bituminous surfacing Figure 5 (size 231 × 650) and 6 (size 201 × 640), courtesy of CEREMA (Centre d'Études et d'Expertise sur les Risques, l'Environnement et l'Aménagement), France, that motivated our work. The two considered slices of bitumen, in addition to long and thin cracks, exhibit high oscillatory patterns and white spots of varying sizes, which makes the straight application of our algorithm difficult. Indeed, in terms of scale, the crack and some of these spots could be comparable and could
<figure><img src="image_4.png" /><figcaption>Fig. 3: Crack detection: $\mu = 0.001$, $\xi_e = 3.5$, $\alpha = 0.14$, $\beta = 0.07$, $\rho = 3.5$, $\varepsilon = 1$, $\gamma = 0.5$, 50 iterations; nonlocal (weights based on Hessian): $\mu = 0.0001$, $\xi_e = 3.5$, $\alpha = 0.14$, $\beta = 0.07$, $\rho = 3.5$, $\varepsilon = 1.5$, $\gamma = 1.0$, $w = 15$, $p = 7$, 550 iterations; nonlocal (only 4 closest neighbors): $\mu = 0.0001$, $\xi_e = 3.5$, $\alpha = 0.06$, $\beta = 0.03$, $\rho = 3.5$, $\varepsilon = 1.5$, $\gamma = 1.0$, 550 iterations; nonlocal (weights based on gradient): $\mu = 0.0001$, $\xi_e = 3.5$, $\alpha = 0.06$, $\beta = 0.03$, $\rho = 3.5$, $\varepsilon = 1.0$, $\gamma = 1.0$, $w = 9$, $p = 5$, 1000 iterations.</figcaption></figure>

pre-processing step is applied. It consists in apprehen-
ding the problem first as an inpainting one ([1]), and
by considering these white spots as missing parts of the
image that need to be filled. This is achieved with the
MATLAB® function **imfill** (https://fr.mathworks.
com/help/images/ref/imfill.html) —to fill holes in
a grayscale image) applied to the inverse image, yiel-
ding an image that serves as input of our algorithm.
In both cases, the cracks are well recovered in the $v_2$
component which does not include superfluous infor-
mation. The edge detector $v_1$ also recovers parts of the
crack but contains spurious information regarding the
problem we address, such as asphalt defect boundaries.
It thus justifies the use of a second order method.
Besides, Figure 5-(g)-(h)-(i) and 6-(g)-(h)-(i) are the
results obtained without considering the decomposition
part. Thanks to Figure 5-(j)-(k)-(l) and 6-(j)-(k)-(l) show-
ing the absolute difference between both results, we ob-
serve that $u$ is less noisy with our method, $v_1$ and $v_2$
also exhibit better contrast with less superfluous infor-
mation. The results obtained with the three versions of
our nonlocal algorithm are comparable with the ones
obtained with our local implementation.

not be properly discriminated, resulting in superfluous
information in the $v_2$ component. Think for instance
of a white spot assimilated to a ball of radius 2 pixels
—if the image domain is the $n \times n$ discretized unit
square, then the scale behaves like $\frac{1}{n}$—), and of a long
thin crack of width 2 pixels and length $k$ pixels ($k \gg 1$)
leading to a similar scale. To circumvent this issue, a
<figure><img src="image_5.png" /><figcaption>(a) Original image</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) $u$ local</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) $v = \operatorname{div} g$ local</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) $u$ nonlocal (weights based on gradient)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(f) $v_2$ nonlocal (weights based on gradient)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(g) $u$ nonlocal (weights based on Hessian)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(h) $v_2$ nonlocal (weights based on Hessian)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(i) $u$ nonlocal (only 4 closest neighbors)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(j) $v_2$ nonlocal (only 4 closest neighbors)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(k) Structures detected by the topological indicator</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(a) Original image</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Filtered image: input of our algorithm</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) $u$ local</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) $v = \operatorname{div} g$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(e) $v_2$ local</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(f) $v_1$ local</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(g) $u_{wd}$ local without decomposition</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(h) $v_{2,wd}$ local without decomposition</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(i) $v_{1,wd}$ local without decomposition</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(j) $|u - u_{wd}|$ local</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(k) $|v_1 - v_{1,wd}|$ local</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(l) $|v_2 - v_{2,wd}|$ local</figcaption></figure>

## 5 Conclusion

This work was devoted to crack detection on bituminous surfacing images. To overcome the limitation of classical first order methods, namely the inability to capture singularities without jumps, and motivated by Drogoul et al.'s observation, a second order approach has been proposed based on a modified elliptic approximation of the Blake-Zisserman functional with a variable encoding the fine structures. A decomposition model has been added to handle the dense and highly
<figure><img src="image_27.png" /><figcaption>Fig. 5: Crack detection: $\mu = 0.0001$, $\xi_{\epsilon} = 3.5$, $\alpha = 0.14$, $\beta = 0.07$, $\rho = 3.5$, $\epsilon = 1$, $\gamma = 0.9$, 270 iterations; nonlocal (weights based on Hessian): $\mu = 0.0001$, $\xi_{\epsilon} = 3.5$, $\alpha = 0.14$, $\beta = 0.07$, $\rho = 3.5$, $\epsilon = 1.0$, $\gamma = 1.0$, $w = 15$, $p = 7$, 500 iterations; nonlocal (only 4 closest neighbors): $\mu = 0.0001$, $\xi_{\epsilon} = 3.5$, $\alpha = 0.14$, $\beta = 0.07$, $\rho = 3.5$, $\epsilon = 1.0$, $\gamma = 1.0$, 500 iterations; nonlocal (weights based on gradient): $\mu = 0.0001$, $\xi_{\epsilon} = 3.5$, $\alpha = 0.14$, $\beta = 0.07$, $\rho = 3.5$, $\epsilon = 1.0$, $\gamma = 1.0$, $w = 9$, $p = 3$, 500 iterations.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 6: Crack detection: $\mu = 0.001$, $\xi_{\epsilon} = 2.5$, $\alpha = 0.1$, $\beta = 0.05$, $\rho = 2.5$, $\epsilon = 1$, $\gamma = 0.9$, 270 iterations; nonlocal (weights based on Hessian): $\mu = 0.0001$, $\xi_{\epsilon} = 3.5$, $\alpha = 0.14$, $\beta = 0.07$, $\rho = 3.5$, $\epsilon = 1.0$, $\gamma = 1.0$, $w = 9$, $p = 3$, 500 iterations; nonlocal (only 4 closest neighbors): $\mu = 0.0001$, $\xi_{\epsilon} = 3.5$, $\alpha = 0.14$, $\beta = 0.07$, $\rho = 3.5$, $\epsilon = 1.0$, $\gamma = 1.0$, 500 iterations; nonlocal (weights based on gradient): $\mu = 0.0001$, $\xi_{\epsilon} = 4.5$, $\alpha = 0.14$, $\beta = 0.07$, $\rho = 4.5$, $\epsilon = 1.0$, $\gamma = 1.0$, $w = 15$, $p = 7$, 2000 iterations.</figcaption></figure>
oscillatory texture of the images and to discriminate
it from the crack. To further improve the quality of
the detection, a nonlocal version of the model has been
theoretically studied with a new Γ-convergence result,
and numerically examined by providing a detailed algo-
rithmic approach and an efficient MPI implementation.
The method has been tested on real images given by
the CEREMA and shows good performance with still
room for improvement. A first lead would consist in dif-
fusing along the obtained cracks in order to highlight
the contrast and fill the missing data.

Appendix A - Proof of Theorem 3: Γ-convergence result

**Step 1.**

We have proved in what precedes that for any $n \in \mathbb{N}^*$, there exists a solution to (2). Let us consider a sequence of such minimizers $(u_n, g_n, Q_n, v_{1,n}, v_{2,n}) \in H^s(\Omega) \cap H_0^1(\Omega) \times H(\text{div}) \times \{Q \in W^{1,\infty}(\Omega) \mid \int_{\Omega} Q dx = 0\} \times W^{1,2}(\Omega, [0, 1]) \times W^{1,2}(\Omega, [0, 1])$, with $s \in [\frac{3}{2}, 2]$.

Let $u \in W^{2,2}(\Omega) \cap W_0^{1,2}(\Omega)$, $g \in H(\text{div})$, $Q \in \{Q \in W^{1,\infty}(\Omega) \mid \int_{\Omega} Q dx = 0\}$, $v_1 \in W^{1,2}(\Omega, [0, 1])$, $v_2 \in W^{1,2}(\Omega, [0, 1])$, then $\forall n \in \mathbb{N}^*$, $\bar{F}_{n,e}(u_n, g_n, Q_n, v_{1,n}, v_{2,n}) \le \bar{F}_{n,e}(u, g, Q, v_1, v_2) \le \rho(1+\kappa_e) \|u\|_{W^{2,2}(\Omega)}^2 + \xi_e \int_{\Omega} (v_1^2(x) + \zeta_e) |\nabla u(x)|^2 dx + \|f - u - \text{div}g\|_{L^2(\Omega)}^2 + \mu \|\nabla Q\|_{L^\infty(\Omega)}^2 + \frac{\gamma}{2} \|g - \nabla Q\|_{L^2(\Omega)}^2 + (\alpha - \beta) \int_{\Omega} \frac{(v_1(x)-1)^2}{4\epsilon} + \epsilon |\nabla v_1(x)|^2 dx + \beta \int_{\Omega} \frac{(v_2(x)-1)^2}{4\epsilon} + \epsilon |\nabla v_2(x)|^2 dx = C < +\infty$ using [18, Lemma 2].

We deduce that:

- $\|\nabla Q_n\|_{L^\infty(\Omega)}$ is uniformly bounded with respect to $n$ and since $\int_\Omega Q_n \,dx = 0$ for any $n \in \mathbb{N}^*$, then using Poincaré-Wirtinger inequality, we get that $(Q_n)$ is uniformly bounded in $W^{1,\infty}(\Omega)$. Thus we can extract a subsequence still denoted by $(Q_n)$ weakly-* converging to $\bar{Q}$ in $W^{1,\infty}(\Omega)$. Besides, the compact embedding $W^{1,\infty}(\Omega) \subset C_0(\bar{\Omega})$ leads to $\int_{\Omega} \bar{Q} \,dx = 0$.

- $\|\nabla v_{1,n}\|_{L^2(\Omega)}$ is uniformly bounded with respect to $n$ and since $v_{1,n} \in W^{1,2}(\Omega, [0, 1])$ for any $n \in \mathbb{N}^*$ and $\Omega$ being bounded, then by using Poincaré-Wirtinger inequality, we get that $(v_{1,n})$ is uniformly bounded in $W^{1,2}(\Omega)$. We can therefore extract a subsequence still denoted by $(v_{1,n})$ weakly converging to $\bar{v}_1$ in $W^{1,2}(\Omega)$ and thanks to the compact embedding $W^{1,2}(\Omega) \subset L^2(\Omega)$, the convergence is also pointwise almost everywhere up to a subsequence, leading to $\bar{v}_1 \in W^{1,2}(\Omega, [0, 1])$.

- Using the same arguments, we get that $v_{2,n} \xrightarrow{n \to +\infty} \bar{v}_2$ in $W^{1,2}(\Omega)$ with $\bar{v}_2 \in W^{1,2}(\Omega, [0, 1])$.

- $\|\nabla u_n\|_{L^2(\Omega)}$ is uniformly bounded with respect to $n$ and since $u_n \in W_0^{1,2}(\Omega)$ for any $n \in \mathbb{N}^*$, then using Poincaré inequality, we can extract a subsequence still denoted by $(u_n)$ weakly converging in $W^{1,2}(\Omega)$ to $\bar{u}$. By continuity of the trace operator, we have $\bar{u} \in W_0^{1,2}(\Omega)$.

- $C \ge \frac{1}{4} \|\text{div} g_n\|_{L^2(\Omega)}^2 - \frac{1}{2} \|u_n\|_{L^2(\Omega)}^2 - \|f\|_{L^2(\Omega)}^2 + \frac{\gamma}{4} \|g_n\|_{L^2(\Omega)}^2 - \frac{\gamma}{2} \|\nabla Q_n\|_{L^2(\Omega)}^2$. So, $(g_n)$ is uniformly bounded in $H(\text{div})$ and we can extract a subsequence still denoted by $(g_n)$ weakly converging to $\bar{g}$ in $H(\text{div})$.

Now, let us show that $(u_{n,e})$ weakly converges in $H_2^{\frac{3}{2}}(\mathbb{R}^2)$ to $\bar{u}_e$.

We set $E_n(h) = \int_{\mathbb{R}^2} |u_{n,e}(x+he_1) - 2u_{n,e}(x) + u_{n,e}(x- he_1)|^2 dx$. It satisfies $E_n(2h) \le 16E_n(h)$, $\forall h \in \mathbb{R}$. By using Fubini-Tonelli theorem, we have

$$
\begin{align*}
& \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{|u_{n,e}(x+he_1)-2u_{n,e}(x)+u_{n,e}(x- he_1)|^2}{|h|^4} \rho_n(h) \, dh \, dx \\
&= \int_{\mathbb{R}} \frac{E_n(h)}{|h|^4} \rho_n(h) \, dh = 2 \int_0^\infty \frac{E_n(h)}{|h|^4} \rho_n(h) \, dh \le C. \text{ We then} \\
& \text{apply [5, Lemma 3.1] with } M = \delta = 1, g(t) = \frac{E_n(t)}{t^{q+1}}, k(t) = \\
& \qquad t^{q-3} \rho_n(t) \text{ and we get}
\end{align*}
$$

$C(1) \int_0^1 t^{q-3} \rho_n(t) dt \le \int_0^1 \frac{E_n(t)}{t^{q+1}} dt \le \int_0^1 \frac{E_n(t)}{t^q} \rho_n(t) dt$. (We will see further that the condition of monotonicity on $k$ is fulfilled). We now need $g$ to verify the assumption of this lemma, that is to say, $g(\frac{t}{2}) \ge g(t)$. We know that $g(\frac{t}{2}) = \frac{E_n(\frac{t}{2})}{t^{q+1}} 2^{q+1} \ge 2^{q-3}g(t)$. Thus if $q \ge 3$, this condition is fulfilled. By using the properties of $\rho_n$, we deduce first that $\int_0^1 \frac{E_n(t)}{t^{q+1}} dt \le C$ with $C$ independent of $n$ for $q = 3$ and $n$ large enough since then $\int_0^1 t^{q-3} \rho_n(t) dt = \int_0^1 \rho_n(t) dt = 1 - \int_1^\infty \rho_n(t) dt$, $\forall n \in \mathbb{N}^*$ with $\lim_{n \to +\infty} \int_1^\infty \rho_n(t) dt = 0$ and so for $n$ large

$$
\begin{align*}
& \text{enough, } \int_0^1 \rho_n(t) dt \in [\frac{1}{2}, 1]. \text{ Then } \int_1^\infty \frac{E_n(t)}{t^{q+1}} dt \le \\
& C' \|u_{n,e}\|_{L^2(\mathbb{R}^2)}^2 \int_1^\infty \frac{1}{t^{q+1}} dt, \text{ C' being a constant and the last} \\
& \text{integral being convergent since } q=3, \text{ resulting in the uniform} \\
& \text{boundedness of } \int_0^\infty \frac{E_n(t)}{t^{q+1}} dt. \text{ Besides, } \int_{\mathbb{R}} \frac{E_n(t)}{|t|^{q+1}} dt = \\
& \qquad \int_{\mathbb{R}} \frac{1}{|t|^{q+1}} \int_{\mathbb{R}^2} |u_{n,e}(x+he_1)-2u_{n,e}(x)+u_{n,e}(x- he_1)|^2 \rho_n(h) \, dx \, dh \\
&= \int_{\mathbb{R}} \frac{1}{|h|^{q+1}} \|\tau_{he_1} u_{n,e} - 2u_{n,e} + \tau_{-he_1} u_{n,e}\|_{L^2(\mathbb{R}^2)}^2 \, dh = \\
& \qquad \int_{\mathbb{R}} \frac{1}{|h|^{q+1}} \int_{\mathbb{R}^2} |e^{2i\pi h \xi_1} - 2 + e^{-2i\pi h \xi_1}|^2 |\mathcal{F}(u_{n,e})|^2(\xi) \, d\xi \, dh \text{ by} \\
& \text{Plancherel's theorem } (\tau \text{ denoting the usual translation ope-} \\
& \text{rator and } \xi = (\xi_1, \xi_2)). \text{ Then one can prove that } \int_{\mathbb{R}} \frac{E_n(t)}{|t|^{q+1}} dt = \\
& C'' \int_{\mathbb{R}^2} |\mathcal{F}(u_{n,e})(\xi)|^2 |\xi_1|^q \int_{\mathbb{R}} \frac{\sin^4(u)}{|u|^{q+1}} du \, d\xi \le C'' \text{ (the constant} \\
& \text{C'' hereafter may change line to line). The generalized inte-} \\
& \text{gral in } u \text{ converges if and only if } q \in [3, 4 \text{ so for } q = 3.} \\
& \text{By using the same arguments in the other direction } (e_2), \text{ we} \\
& \text{get that } |\cdot|^{\frac{3}{2}} \mathcal{F}(u_{n,e})(\cdot) \in L^2(\mathbb{R}^2) \text{ and so } u_{n,e} \in H_2^{\frac{3}{2}}(\mathbb{R}^2) \\
& \text{(being a Hilbert space) and is uniformly bounded for the as-} \\
& \text{sociated norm for } n \text{ large enough. There exists a subsequence} \\
& \text{still denoted by } (u_{n,e}) \text{ weakly converging to } \bar{u} \text{ in } H_2^{\frac{3}{2}}(\mathbb{R}^2).} \\
& \text{Besides, we know that } u_{n,e} = u_n \text{ on } \Omega \text{ and } D^{(1,0)} u_{n,e} = \\
& \qquad (D^{(1,0)} u_n)_e = D^{(1,0)} u_n \text{ on } \Omega, \text{ and } D^{(0,1)} u_{n,e} = (D^{(0,1)} u_n)_e = \\
& \qquad D^{(0,1)} u_n \text{ on } \Omega. \text{ Thus, } \|u_n\|_{W^{\frac{3}{2},2}(\Omega)}^2 = \|u_n\|_{W^{1,2}(\Omega)}^2 + \\
& \qquad \int_{\Omega} \int_{\Omega} \frac{|\nabla u_n(y) - \nabla u_n(x)|^2}{|x-y|^3} \, dx \, dy \le \\
& C + \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(y) - \nabla u_{n,e}(x)|^2}{|x-y|^3} \, dx \, dy, \text{ with } C \text{ independent of} \\
& \text{n. From [19, Lemma 4.33, p. 200], we know that} \\
& \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(x) - \nabla u_{n,e}(y)|^2}{|x-y|^3} \, dx \, dy < \infty \\
& \Leftrightarrow \int_{\mathbb{R}} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(x) - \nabla u_{n,e}(x+he_1)|^2}{|h|^2} \, dx \, dh < \infty \text{ and} \\
& \int_{\mathbb{R}} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(x) - \nabla u_{n,e}(x+he_2)|^2}{|h|^2} \, dx \, dh < \infty. \text{ Let us now}
\end{align*}
$$

$$
\begin{align*}
& \int_{\Omega} \int_{\Omega} \frac{|\nabla u_n(y) - \nabla u_n(x)|^2}{|x-y|^3} \, dx \, dy \le \\
& C + \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(y) - \nabla u_{n,e}(x)|^2}{|x-y|^3} \, dx \, dy, \text{ with } C \text{ independent of} \\
& n. \text{ From [19, Lemma 4.33, p. 200], we know that} \\
& \int_{\mathbb{R}^2} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(x) - \nabla u_{n,e}(y)|^2}{|x-y|^3} \, dx \, dy < \infty \\
& \Leftrightarrow \int_{\mathbb{R}} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(x) - \nabla u_{n,e}(x+he_1)|^2}{|h|^2} \, dx \, dh < \infty \text{ and} \\
& \int_{\mathbb{R}} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(x) - \nabla u_{n,e}(x+he_2)|^2}{|h|^2} \, dx \, dh < \infty. \text{ Let us now}
\end{align*}
$$

prove that

$$
\int_{\mathbb{R}} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(x) - \nabla u_{n,e}(x + he_1)|^2}{|h|^2} \, dx \, dh < \infty \text{ and} \\
\int_{\mathbb{R}} \int_{\mathbb{R}^2} \frac{|\nabla u_{n,e}(x) - \nabla u_{n,e}(x + he_2)|^2}{|h|^2} \, dx \, dh < \infty \text{ independently of } n.
$$

$$
\begin{align*}
& \text{We have } \int_{\mathbb{R}} \frac{1}{|h|^2} \|\tau_{he_1}(\nabla u_{n,e}) - \nabla u_{n,e}\|_{L^2(\mathbb{R}^2)}^2 \, dh = \\
& \qquad \bar{C} \int_{\mathbb{R}} \frac{\sin^2(u)}{|u|^2} \int_{\mathbb{R}^2} |\xi_1| (|\xi_1|^2 + |\xi_2|^2) |\mathcal{F}(u_{n,e})(\xi)|^2 \, d\xi \, du \\
& \le C \|u_{n,e}\|_{H_2^{\frac{3}{2}}(\mathbb{R}^2)}^2 \text{ by using Plancherel's theorem and with} \\
& \text{C independent of } n. \text{ By doing the same computations in}
\end{align*}
$$
the other direction, we prove that $||u_n||_{W^{\frac{3}{2},2}(\Omega)}$ is uniformly bounded and so up to a subsequence, $u_n \xrightarrow{n \to +\infty} \bar{u}$ in $W^{\frac{3}{2},2}(\Omega)$, since $W^{\frac{3}{2},2}(\Omega) \subset W^{1,2}(\Omega)$ and the mapping $T: W^{\frac{3}{2},2}(\Omega) \to W^{1,2}(\Omega)$, $u \mapsto Tu = u$ is linear and continuous for the strong topology, so for the weak topology ([19, Theorem III.9, p. 39]). As $W^{\frac{3}{2},2}(\Omega) \subset C_b^{0,\lambda}(\Omega)$, with $\lambda < \frac{1}{2}$, then $(u_n)$ strongly converges to $\bar{u}$ in $C_b^{0,\lambda}(\Omega)$ and so pointwise everywhere on $\Omega$. Then $\bar{u} = \bar{u}$ on $\Omega$ and $\bar{u} = 0$ on $\partial\Omega$, by uniqueness of the weak limit. Now, $H_2^{\frac{3}{2}}(\mathbb{R}^2) \subset L^2(\mathbb{R}^2) \subset S'(\mathbb{R}^2) \subset D'(\mathbb{R}^2)$ with continuous embeddings. $\forall \varphi \in D(\mathbb{R}^2)$, $\underbrace{\int_{\mathbb{R}^2} (u_{n,e} - \bar{u})\varphi dx}_{\xrightarrow{n \to +\infty} 0} = \underbrace{\int_{\Omega} (u_n - \bar{u})\varphi dx}_{\xrightarrow{n \to +\infty} 0} + \underbrace{\int_{\mathbb{R}^2 \setminus \Omega} (u_{n,e} - \bar{u})\varphi dx}_{\xrightarrow{n \to +\infty} 0}$

Consequently, $\forall \varphi \in D(\mathbb{R}^2)$, $\int_{\mathbb{R}^2 \setminus \Omega} \bar{u}\varphi dx = \int_{\mathbb{R}^2 \setminus \Omega} \bar{u}\varphi dx = 0$ since $\bar{u} \in H_2^{\frac{3}{2}}(\mathbb{R}^2) \circ C(\mathbb{R}^2)$. In particular, $\forall \varphi \in D(\mathbb{R}^2 \setminus \bar{\Omega})$, $\int_{\mathbb{R}^2 \setminus \bar{\Omega}} \bar{u}\varphi dx = 0$, meaning that $\bar{u} = 0$ on $\mathbb{R}^2 \setminus \bar{\Omega}$ in the sense of distributions. Due to the continuity of $\bar{u}$, we deduce that $\bar{u} = 0$ everywhere on $\mathbb{R}^2 \setminus \Omega$ and so $\bar{u} = \bar{u}_e \in H_2^{\frac{3}{2}}(\mathbb{R}^2)$. Now, let us prove that $\bar{u} \in H^2(\Omega) \cap H_0^1(\Omega)$ and $\bar{u}_e \in H^2(\mathbb{R}^2)$. Since $u_n \xrightarrow{n \to +\infty} \bar{u}$ in $L^2(\Omega)$ then $u_{n,e} \xrightarrow{n \to +\infty} \bar{u}_e$ in $L^2(\mathbb{R}^2)$. Let $\varphi \in C_0^\infty(\mathbb{R}^2)$. We denote by $\bar{u}_e$ and by $u_{n,e}$ the extensions by 0 of $\bar{u}$ and $u_n$ for any $n \in \mathbb{N}^*$ on $\mathbb{R}^2$. Since $\Omega$ is of class $C^2$ then $\bar{u}_e, u_{n,e} \in W^{1,2}(\mathbb{R}^2)$, for any $n \in \mathbb{N}^*$. We have

$$ \begin{aligned} & \left| \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{\varphi(x + he_1) - 2\varphi(x) + \varphi(x - he_1)}{h^2} \rho_n(h) dh(u_{n,e}(x) + (\bar{u}_e(x) - u_{n,e}(x))) dx \right|, \\ & \le \left| \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)}{h^2} \rho_n(h) dh(u_{n,e}(x)) dx \right| \\ & \quad + \left| \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)}{h^2} \rho_n(h) dh(\bar{u}_e(x) - u_{n,e}(x)) dx \right|, \\ & \le \left( \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{(\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1))^2}{h^4} \rho_n(h) dh dx \right)^{\frac{1}{2}}. \\ & \quad \left( \int_{\mathbb{R}^2} \int_{\mathbb{R}} |u_{n,e}(x) - \bar{u}_e(x)|^2 \rho_n(h) dh dx \right)^{\frac{1}{2}} \\ & \quad + \left| \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{u_{n,e}(x+he_1)-2u_{n,e}(x)+u_{n,e}(x-he_1)}{h^2} \rho_n(h) dh(\varphi(x)) dx \right|, \end{aligned} $$

owing to Hölder's inequality with respect to the measure $\rho_n(h) dh \, dx$. Since $\varphi \in C_0^\infty(\mathbb{R}^2)$ (justification of the integrability and the possibility to change the order of integration is postponed - in fact, it is not difficult to check that Fubini-Lebesgue theorem applies),

$$ \begin{aligned} & \le ||D^{(2,0)}\varphi||_{L^2(\mathbb{R}^2)} ||u_{n,e} - \bar{u}_e||_{L^2(\mathbb{R}^2)} \\ & \quad + \left( \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{(u_{n,e}(x+he_1)-2u_{n,e}(x)+u_{n,e}(x-he_1))^2}{h^4} \rho_n(h) dh dx \right)^{\frac{1}{2}} \\ & \quad \left( \int_{\mathbb{R}^2} \int_{\mathbb{R}} \varphi^2(x) \rho_n(h) dh dx \right)^{\frac{1}{2}}, \\ & \le ||D^{(2,0)}\varphi||_{L^2(\mathbb{R}^2)} ||u_{n,e} - \bar{u}_e||_{L^2(\mathbb{R}^2)} \\ & \quad + \left( \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{(u_{n,e}(x+he_1)-2u_{n,e}(x)+u_{n,e}(x-he_1))^2}{h^4} \rho_n(h) dh dx \right)^{\frac{1}{2}} ||\varphi||_{L^2(\mathbb{R}^2)}, \\ & \le ||D^{(2,0)}\varphi||_{L^2(\mathbb{R}^2)} ||u_{n,e} - \bar{u}_e||_{L^2(\mathbb{R}^2)} + \sqrt{C} ||\varphi||_{L^2(\mathbb{R}^2)}. \end{aligned} $$

Since $\varphi \in C_0^\infty(\mathbb{R}^2)$, then for almost every $x \in \mathbb{R}^2$, $\int_{\mathbb{R}} \frac{\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)}{h^2} \rho_n(h) dh \bar{u}_e(x) \xrightarrow{n \to +\infty} D^{(2,0)}\varphi(x)\bar{u}_e(x)$ thanks to [18, Theorem 9]. Besides, $\forall n \in \mathbb{N}^*$ and almost every $x \in \mathbb{R}^2$,

$$ \begin{aligned} & \left| \int_{\mathbb{R}} \frac{\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)}{h^2} \rho_n(h) dh \bar{u}_e(x) \right| \le |\bar{u}_e(x)|. \\ & \int_{\mathbb{R}} \frac{|\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)|}{|h|^2} \rho_n(h) dh \end{aligned} $$

by applying Jensen's inequality with respect to the measure $\rho_n(h) dh$ (see [40, Theorem 139, p. 56]). But $\bar{u}_e \in L^2(\mathbb{R}^2)$ with $||\bar{u}_e||_{L^2(\mathbb{R}^2)} = ||\bar{u}||_{L^2(\Omega)} \ge c||\bar{u}||_{L^1(\Omega)} = c||\bar{u}_e||_{L^1(\mathbb{R}^2)}$, since $\bar{u}$ is extended by 0 outside $\Omega$. As $\varphi \in C_0^\infty(\mathbb{R}^2)$, then

$$ \begin{aligned} & \int_{\mathbb{R}} \frac{|\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)|}{|h|^2} \rho_n(h) dh \le L^1(\mathbb{R}^2) \\ & \int_{\mathbb{R}} \frac{|\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)|}{|h|^2} \rho_n(h) dh \le \end{aligned} $$

$$ \begin{aligned} & \int_{\mathbb{R}} \int_0^1 \int_0^1 |D^{(2,0)}\varphi(x + (t + s - 1)he_1)| \rho_n(h) ds dt dh \\ & \le ||D^{(2,0)}\varphi||_{L^\infty(\mathbb{R}^2)}. \end{aligned} $$

$$ \begin{aligned} & \int_{\mathbb{R}^2} \left( \int_{\mathbb{R}} \frac{\varphi(x + he_1) - 2\varphi(x) + \varphi(x - he_1)}{h^2} \rho_n(h) dh \right)^2 dx \\ & \le \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{(\varphi(x + he_1) - 2\varphi(x) + \varphi(x - he_1))^2}{h^4} \rho_n(h) dh dx \\ & \le ||D^{(2,0)}\varphi||_{L^2(\mathbb{R}^2)}^2 < +\infty. \end{aligned} $$

Thus $\int_{\mathbb{R}} \frac{\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)}{h^2} \rho_n(h) dh \in L^2(\mathbb{R}^2)$, and $\bar{u}_e$ being $L^2(\mathbb{R}^2)$, it leads to $\int_{\mathbb{R}} \frac{\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)}{h^2} \rho_n(h) dh \bar{u}_e(x) \in L^1(\mathbb{R}^2)$ with $\left| \int_{\mathbb{R}} \frac{\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)}{h^2} \rho_n(h) dh \bar{u}_e(x) \right| \le |\bar{u}_e(x)||D^{(2,0)}\varphi||_{L^\infty(\mathbb{R}^2)} \in L^1(\mathbb{R}^2)$. From the dominated convergence theorem, we can conclude that

$$ \begin{aligned} & \lim_{n \to +\infty} \int_{\mathbb{R}} \bar{u}_e(x) \int_{\mathbb{R}} \frac{\varphi(x+he_1)-2\varphi(x)+\varphi(x-he_1)}{h^2} \rho_n(h) dh dx = \\ & \int_{\mathbb{R}^2} \bar{u}_e(x) D^{(2,0)}\varphi(x) dx. \end{aligned} $$
By letting $n$ tend to infinity in the previous inequality, we get $\int_{\mathbb{R}^2} \bar{u}_e(x) D^{(2,0)}\varphi(x) dx \le \sqrt{C} ||\varphi||_{L^2(\mathbb{R}^2)}$. Eventually, $D^{(2,0)}\bar{u}_e \in L^2(\mathbb{R}^2)$ using [12, Proposition VIII.3]. By applying the same reasoning, we also get $D^{(0,2)}\bar{u}_e \in L^2(\mathbb{R}^2)$. Since $\bar{u}_e \in L^2(\mathbb{R}^2)$ then it is a tempered distribution and so are its successive derivatives. We can now take the Fourier transform $\hat{D^{(1,1)}\bar{u}_e}(\xi) = -\xi_1\xi_24\pi^2\hat{u}_e(\xi)$, $\hat{D^{(0,2)}\bar{u}_e}(\xi) = -\xi_2\xi_24\pi^2\hat{u}_e(\xi)$, $\hat{D^{(2,0)}\bar{u}_e}(\xi) = -\xi_1\xi_14\pi^2\hat{u}_e(\xi)$. Since $D^{(0,2)}\bar{u}_e \in L^2(\mathbb{R}^2)$ and $D^{(2,0)}\bar{u}_e \in L^2(\mathbb{R}^2)$, then $-\xi_2\xi_2\hat{u}_e(\xi) \in L^2(\mathbb{R}^2)$ and $-\xi_1\xi_1\hat{u}_e(\xi) \in L^2(\mathbb{R}^2)$. But

$$ \int_{\mathbb{R}^2} 4\pi^2 \xi_1 \xi_2 \hat{u}_e^2(\xi) \, d\xi \le 2\pi^2 \int_{\mathbb{R}^2} \xi_1^4 \hat{u}_e^2(\xi) \, d\xi \\ + 2\pi^2 \int_{\mathbb{R}^2} \xi_2^4 \hat{u}_e^2(\xi) \, d\xi < +\infty, $$

which means that $\hat{D^{(1,1)}\bar{u}_e}(\xi) \in L^2(\mathbb{R}^2)$ and by Plancherel's theorem we can conclude that $\hat{D^{(1,1)}\bar{u}_e} \in L^2(\mathbb{R}^2)$. This proves that $\bar{u}_e \in W^{2,2}(\mathbb{R}^2)$ since $\bar{u}_e \in W^{1,2}(\mathbb{R}^2)$ by construction. As $\Omega \in C^2$ and $\bar{u}_e$ is the extension of $\bar{u}$ by 0 outside $\Omega$, then $\bar{u} \in W^{2,2}(\Omega) \cap W_0^{1,2}(\Omega)$.

### Step 2.

By definition of the sequence $(u_n, g_n, Q_n, v_{1,n}, v_{2,n})$, we have $\forall n \in \mathbb{N}^*$, $\bar{F}_{n,\varepsilon}(u_n, g_n, Q_n, v_{1,n}, v_{2,n}) \le \bar{F}_{n,\varepsilon}(\bar{u}, \bar{g}, \bar{Q}, \bar{v}_1, \bar{v}_2)$ and by taking the lim sup when $n \to +\infty$, it yields

$$ \limsup_{n \to +\infty} \bar{F}_{n,\varepsilon}(u_n, g_n, Q_n, v_{1,n}, v_{2,n}) \le \bar{F}_{\varepsilon}(\bar{u}, \bar{g}, \bar{Q}, \bar{v}_1, \bar{v}_2). $$

Indeed, thanks to [18, Theorem 9], we know that $(\bar{v}_2, e)(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+he_1)-2u(x)+u(x-he_1)|^2}{|h|^4} \rho_n(h) dh \xrightarrow{n \to +\infty} (\bar{v}_2, e)(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2(x)$ everywhere in $\mathbb{R}^2$ and for all $u \in C_c^4(\mathbb{R}^2)$. Without loss of generality, we assume that $\text{supp}(u) \subset$
$B(0, R)$ with $R > 0$. We now aim to prove that $\forall \varepsilon > 0, \exists L = L(\varepsilon) > 1$ (we believe that the confusion with the $\varepsilon$ from the elliptic approximation is not possible) such that

$$ \sup_n \int_{B(0,LR)^c} \left| (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+he_1)-2u(x)+u(x-he_1)|^2}{|h|^4} \rho_n(h) dh \right| dx \le \varepsilon. \text{ We have that} $$

$$ \begin{aligned} & \int_{B(0,LR)^c} \left| (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+he_1)-2u(x)+u(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh \right| dx \\ & \le \frac{(1+\kappa_\varepsilon)4}{(L-1)^4 R^4} \|u\|_{L^2(\mathbb{R}^2)}^2 \text{ and the conclusion follows. As} \end{aligned} $$

$$ \begin{aligned} & \int_{\mathbb{R}^2} \left| (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+he_1)-2u(x)+u(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh - (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u(x)|^2 \right| dx \\ & = \int_{B(0,LR)} \left| (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+\bar{h}e_1)-2u(x)+u(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh - (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u(x)|^2 \right| dx \\ & \qquad + \int_{B(0,LR)^c} \left| (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+\bar{h}e_1)-2u(x)+u(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh - (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u(x)|^2 \right| dx, \\ & \le \int_{B(0,LR)} \left| (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+\bar{h}e_1)-2u(x)+u(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh - (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u(x)|^2 \right| dx + \varepsilon. \end{aligned} $$

We know that $(\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+\bar{h}e_1)-2u(x)+u(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh \xrightarrow{n \to +\infty} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2(x)$ everywhere in $B(0, LR)$ and
$$ \left| (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+\bar{h}e_1)-2u(x)+u(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh - (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2(x) \right| \le 2(1+\kappa_\varepsilon) \|D^{(2,0)}u\|_{L^\infty(\mathbb{R}^2)}^2 \in L^1(B(0, LR)). $$ Using the dominated convergence theorem, we get that

$$ \begin{aligned} & \int_{B(0,LR)} \left| (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+\bar{h}e_1)-2u(x)+u(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh - (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u(x)|^2 \right| dx \xrightarrow{n \to +\infty} 0 \text{ and by letting } \varepsilon \\ & \text{tend to 0, we conclude that} \end{aligned} $$

$$ \begin{aligned} & \int_{\mathbb{R}^2} \left| (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) \int_{\mathbb{R}} \frac{|u(x+\bar{h}e_1)-2u(x)+u(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh - (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u(x)|^2 \right| dx \xrightarrow{n \to +\infty} 0. \text{ We now extend} \\ & \text{this result to } u \in H^2(\mathbb{R}^2). \text{ Let } \varepsilon > 0. \text{ By density, there exists } v_\varepsilon \in C_0^\infty(\mathbb{R}^2) \text{ such that } \|D^{(2,0)}u - D^{(2,0)}v_\varepsilon\|_{L^2(\mathbb{R}^2)} \le \varepsilon. \text{ Let us set } u_n(x,h) = (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon)^{\frac{1}{2}} \frac{u(x+\bar{h}e_1)-2u(x)+u(x-\bar{h}e_1)}{\bar{h}^2} \rho_n^{\frac{1}{2}}(h) \\ & \in L^2(\mathbb{R}^2 \times \mathbb{R}) \text{ and } v_{n,\varepsilon}(x,h) = (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon)^{\frac{1}{2}} \frac{v_\varepsilon(x+\bar{h}e_1)-2v_\varepsilon(x)+v_\varepsilon(x+\bar{h}e_1)}{\bar{h}^2} \rho_n^{\frac{1}{2}}(h) \in L^2(\mathbb{R}^2 \times \mathbb{R}). \text{ We have} \end{aligned} $$

$$ \begin{aligned} & \|u_n - v_{n,\varepsilon}\|_{L^2(\mathbb{R}^2 \times \mathbb{R})} \le (1 + \kappa_\varepsilon)^{\frac{1}{2}} \|D^{(2,0)}u - D^{(2,0)}v_\varepsilon\|_{L^2(\mathbb{R}^2)} \\ & \le (1 + \kappa_\varepsilon)^{\frac{1}{2}} \varepsilon, \end{aligned} $$

yielding, thanks to the second triangle inequality,

$$ \|u_n\|_{L^2(\mathbb{R}^2 \times \mathbb{R})} - \|v_{n,\varepsilon}\|_{L^2(\mathbb{R}^2 \times \mathbb{R})} \le (1 + \kappa_\varepsilon)^{\frac{1}{2}} \varepsilon. $$

Now,

$$ \begin{aligned} & \left| \|u_n\|_{L^2(\mathbb{R}^2 \times \mathbb{R})} - \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2 dx \right)^{1/2} \right| \\ & \le \|u_n\|_{L^2(\mathbb{R}^2 \times \mathbb{R})} - \|v_{n,\varepsilon}\|_{L^2(\mathbb{R}^2 \times \mathbb{R})} \quad (A.1) \\ & + \left| \|v_{n,\varepsilon}\|_{L^2(\mathbb{R}^2 \times \mathbb{R})} - \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon|^2 dx \right)^{1/2} \right| \\ & + \left| \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon|^2 dx \right)^{1/2} - \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2 dx \right)^{1/2} \right|. \end{aligned} $$

Also,

$$ \begin{aligned} & \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon + D^{(2,0)}u|^2 dx = \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon|^2 dx \\ & + \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2 dx + 2 \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) D^{(2,0)}u D^{(2,0)}v_\varepsilon dx, \\ & \le \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon|^2 dx + \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2 dx + \\ & + 2 \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2 dx)^{1/2} \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon|^2 dx \right)^{1/2}, \end{aligned} $$

from Cauchy-Schwarz inequality. It follows that

$$ \begin{aligned} & \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon + D^{(2,0)}u|^2 dx \right)^{1/2} \\ & \le \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2 dx \right)^{1/2} \\ & + \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon|^2 dx \right)^{1/2}. \end{aligned} $$

It results from this that:

$$ \begin{aligned} & \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon|^2 dx \right)^{1/2} \\ & = \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon - D^{(2,0)}u + D^{(2,0)}u|^2 dx \right)^{1/2}, \\ & \le \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon - D^{(2,0)}u|^2 dx \right)^{1/2} \\ & + \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2 dx \right)^{1/2}. \end{aligned} $$

Exchanging the role of $v_\varepsilon$ and $u$, we obtain:

$$ \begin{aligned} & \left| \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon|^2 dx \right)^{1/2} - \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}u|^2 dx \right)^{1/2} \right| \\ & \le \left( \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_\varepsilon) |D^{(2,0)}v_\varepsilon - D^{(2,0)}u|^2 dx \right)^{1/2}, \\ & \le (1 + \kappa_\varepsilon)^{\frac{1}{2}} \varepsilon. \end{aligned} $$

Combining these different elements in (A.1) yields the desired result. Since $\Omega \in \mathcal{C}^2$, then $\bar{u}_\varepsilon$ is the extension by 0 of $\bar{u}$ from $W^{2,2}(\Omega)$ to $W^{2,2}(\mathbb{R}^2)$. So, $D^{(2,0)}\bar{u}_\varepsilon = (D^{(2,0)}\bar{u})_e$ and $D^{(2,0)}\bar{u}_\varepsilon(x) = 0$ almost everywhere on $\mathbb{R}^2 \setminus \Omega$. We finally have $\int_{\mathbb{R}^2} \int_{\mathbb{R}} (\bar{v}_{2,e}^2(x)^2 + \kappa_\varepsilon) |\bar{u}_\varepsilon, n(x,h)|^2 dh dx \xrightarrow{n \to +\infty} \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x)^2 + \kappa_\varepsilon) |D^{(2,0)}\bar{u}_\varepsilon(x)|^2 dx \le \int_{\Omega} (\bar{v}_{2,e}^2(x)^2 + \kappa_\varepsilon) |D^{(2,0)}\bar{u}_\varepsilon(x)|^2 dx.$

### Step 3

It remains to prove that $\bar{F}_\varepsilon(\bar{u}, \bar{g}, \bar{Q}, \bar{v}_1, \bar{v}_2) \le \liminf_{n \to +\infty} \bar{F}_{n,\varepsilon}(u_n, g_n, Q_n, v_1, v_2)$.

$$ \int_{\Omega} (\bar{v}_2(x)^2 + \kappa_\varepsilon) |D^{(2,0)}\bar{u}(x)|^2 dx $$
$$
\le \liminf_{n \to +\infty} \int_{\mathbb{R}^2} (v_2^2, n, e(x) + \kappa_\epsilon) \\
\int_{\mathbb{R}} \frac{|u_{n,e}(x+he_1)-2u_{n,e}(x)+u_{n,e}(x-\bar{h}e_1)|^2}{h^4} \rho_n(h) \, dh \, dx \text{ and} \\
\int_{\Omega} (\bar{v}_2^2(x) + \kappa_\epsilon) |D^{(0,2)}\bar{u}(x)|^2 \, dx \le \liminf_{n \to +\infty} \int_{\mathbb{R}^2} (v_2^2, n, e(x) + \kappa_\epsilon) \\
\int_{\mathbb{R}} \frac{|u_{n,e}(x+he_2)-2u_{n,e}(x)+u_{n,e}(x-\bar{h}e_2)|^2}{h^4} \rho_n(h) \, dh \, dx. \text{ Let } \eta \in C_0^\infty(\mathbb{R}^2) \text{ be a non-negative radial function satisfying } \int_{\mathbb{R}^2} \eta \, dx = 1, \text{ supp}(\eta) \subset B(0,1). \text{ We then define a regularized function associated with } f \text{ by } f_\delta(x) = \frac{1}{\delta^2} \int_{\mathbb{R}^2} f(y)\eta(\frac{x-y}{\delta}) \, dy, \forall x \in \mathbb{R}^2. \text{ We focus on the direction } e_1 \text{ and compute with } \delta > 0 \text{ and } \delta' > 0:
$$

$$
\int_{\mathbb{R}^2} (v_2^2, n, e, \delta'(x) + \kappa_\epsilon) \\
\int_{\mathbb{R}} \frac{|u_{n,e,\delta}(x+he_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-\bar{h}e_1)|^2}{|h|^4} \rho_n(h) \, dh \, dx = \\
\int_{\mathbb{R}^2} (v_2^2, n, e, \delta'(x) + \kappa_\epsilon) \\
\int_{\mathbb{R}} \frac{|\int_{B(0,\delta)} \frac{u_{n,e}(x+\bar{h}e_1-z)-2u_{n,e}(x-z)+u_{n,e}(x-\bar{h}e_1-z)}{\delta^2|\bar{h}|^4} \eta(\frac{z}{\delta}) dz|^2}{|\bar{h}|^4} \rho_n(h) \, dh \, dx.
$$

We apply Jensen's inequality with respect to the measure $\frac{\eta(\frac{z}{\delta})}{\delta^2} dz$, yielding

$$
\le \int_{\mathbb{R}^2} (v_2^2, n, e, \delta'(x) + \kappa_\epsilon) \\
\int_{\mathbb{R}} \int_{B(0,\delta)} \frac{|u_{n,e}(x+\bar{h}e_1-z)-2u_{n,e}(x-z)+u_{n,e}(x-\bar{h}e_1-z)|^2}{\delta^2|\bar{h}|^4} \eta(\frac{z}{\delta}) \rho_n(h) \, dz \, dh \, dx.
$$

We make the following change of variable: $v = x - z$,

$$
\le \int_{\mathbb{R}^2} \int_{\mathbb{R}} \int_{B(0,\delta)} (v_2^2, n, e, \delta'(v+z) + \kappa_\epsilon) \\
\frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{\delta^2|\bar{h}|^4} \eta(\frac{z}{\delta}) \rho_n(h) \, dz \, dh \, dv.
$$

Since $v_{2,n,e,\delta'} \in C^\infty(\mathbb{R}^2)$, then $v_{2,n,e,\delta'}^2(v+z) = v_{2,n,e,\delta'}^2(v) + \int_0^1 \langle \nabla v_{2,n,e,\delta'}^2(v+sz), z \rangle ds$ and by introducing it in the previous inequality, we get:

$$
\le \int_{\mathbb{R}^2} \int_{\mathbb{R}} \int_{B(0,\delta)} (v_2^2, n, e, \delta'(v) + \kappa_\epsilon) \\
\frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{\delta^2|\bar{h}|^4} \eta(\frac{z}{\delta}) \rho_n(h) \, dz \, dh \, dv \\
+ \int_{\mathbb{R}^2} \int_{\mathbb{R}} \int_{B(0,\delta)} \int_0^1 ((\nabla v_{2,n,e,\delta'}^2(v+sz), z)) \, ds \\
\frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{\delta^2|\bar{h}|^4} \eta(\frac{z}{\delta}) \rho_n(h) \, dz \, dh \, dv.
$$

We now integrate with respect to $z$ in the first integral and use Cauchy-Schwarz inequality, and the change of variable $u = \frac{z}{\delta}$ in the second one, after bounding the component above by $||\nabla v_{2,n,e,\delta'}^2||_{L^\infty(\mathbb{R}^2)}$:

$$
\le \int_{\mathbb{R}^2} \int_{\mathbb{R}} (v_2^2, n, e, \delta'(v) + \kappa_\epsilon) \\
\frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) \, dh \, dv \\
+ \int_{\mathbb{R}^2} \int_{\mathbb{R}} \int_{B(0,1)} (\|\nabla v_2^2, n, e, \delta'\|_{L^\infty(\mathbb{R}^2)} \delta|u|) \\
\frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \eta(u) \rho_n(h) \, du \, dh \, dv, \\
\le \int_{\mathbb{R}^2} \int_{\mathbb{R}} (v_2^2, n, e, \delta'(v) + \kappa_\epsilon) \\
\frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) \, dh \, dv \\
+ \|\nabla v_2^2, n, e, \delta'\|_{L^\infty(\mathbb{R}^2)} \\
\delta \int_{\mathbb{R}^2} \int_{\mathbb{R}} \frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) \, dh \, dv, \\
\le \int_{\mathbb{R}^2} \int_{\mathbb{R}} (v_2^2, n, e, \delta'(v) + \kappa_\epsilon) \\
\frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) \, dh \, dv \\
+ \|\nabla v_2^2, n, e, \delta'\|_{L^\infty(\mathbb{R}^2)} \delta C, \text{ from the coercivity inequality.}
$$

Since we assumed that $v_{2,n} \in W^{1,\infty}(\Omega, [0,1])$ with $\sup_{n \in \mathbb{N}^*} \|v_{2,n,e}\|_{W^{1,\infty}(\mathbb{R}^2)} \le \sup_{n \in \mathbb{N}^*} \|v_{2,n}\|_{W^{1,\infty}(\Omega)} \le C_1 < \infty$ and so $\sup_{n \in \mathbb{N}^*} \|\nabla v_{2,n,e}\|_{L^\infty(\mathbb{R}^2)} \le C_1 + 1$ where $v_{2,n,e}$ is the extension of $v_{2,n}$ from $W^{1,2}(\Omega, [0,1]) \cap W^{1,\infty}(\Omega)$ to $W^{1,2}(\mathbb{R}^2, [0,1]) \cap W^{1,\infty}(\mathbb{R}^2)$ and $\bar{v}_{2,e}$ is the extension of $\bar{v}_2$ from $W^{1,2}(\Omega, [0,1])$ to $W^{1,2}(\mathbb{R}^2, [0,1])$, then we have thanks to [22, Theorem 1 p.123] that $\forall x \in \mathbb{R}^2$, $D^{(1,0)}v_{2,n,e,\delta'}(x) = \eta_{\delta'} * D^{(1,0)}v_{2,n,e}(x) = \int_{\mathbb{R}^2} \frac{1}{\delta^{\frac{1}{2}}} D^{(1,0)}v_{2,n,e}(y)\eta(\frac{x-y}{\delta'}) \, dy \le \|D^{(1,0)}v_{2,n,e}\|_{L^\infty(\mathbb{R}^2)}$, $D^{(0,1)}v_{2,n,e,\delta'}(x) = \eta_{\delta'} * D^{(0,1)}v_{2,n,e}(x) = \int_{\mathbb{R}^2} \frac{1}{\delta^2} D^{(0,1)}v_{2,n,e}(y)\eta(\frac{x-y}{\delta'}) \, dy \le \|D^{(0,1)}v_{2,n,e}\|_{L^\infty(\mathbb{R}^2)}$ and so $\|\nabla v_{2,n,e,\delta'}\|_{L^\infty(\mathbb{R}^2)} \le \|\nabla v_{2,n,e}\|_{L^\infty(\mathbb{R}^2)} \le C_1 + 1$. Eventually, we get:

$$
\int_{\mathbb{R}^2} (v_2, n, e, \delta'(x)^2 + \kappa_\epsilon) \\
\int_{\mathbb{R}} \frac{|u_{n,e,\delta}(x+\bar{h}e_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) \, dh \, dx \\
\le \int_{\mathbb{R}^2} \int_{\mathbb{R}} (v_2^2, n, e, \delta'(v) + \kappa_\epsilon) \\
\frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) \, dh \, dv \\
+ 2(C_1 + 1)\delta C.
$$

As $W^{1,\infty}(\mathbb{R}^2) \subset C(\mathbb{R}^2)$, then $v_{2,n,e} \in C(\mathbb{R}^2)$ and [22, Theorem 1 p.123] gives us that $v_{2,n,e,\delta'} \xrightarrow{\delta' \to 0} v_{2,n,e}$ uniformly on compact subsets of $\mathbb{R}^2$ and so pointwise almost everywhere on $\mathbb{R}^2$ as $\mathbb{R}^2$ is locally compact. We thus have $(v_{2,n,e,\delta'}^2(v) + \kappa_\epsilon) \frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) \xrightarrow{\delta' \to 0} (v_2^2, n, e(v) + \kappa_\epsilon) \frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h)$ and $(v_{2,n,e,\delta'}^2(v) + \kappa_\epsilon) \frac{|u_{n,e,\delta}(v+\bar{h}e_1)-2u_{n,e,\delta}(v)+u_{n,e,\delta}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) \xrightarrow{\delta' \to 0} (v_2^2, n, e(v) + \kappa_\epsilon) \frac{|u_{n,e,\delta}(v+\bar{h}e_1)-2u_{n,e,\delta}(v)+u_{n,e,\delta}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h)$ almost
everywhere on $\mathbb{R}^2$ and everywhere on $\mathbb{R}$. Besides,

$$ \int_{\mathbb{R}} (v_{2,n,e,\delta'}^2(v) + \kappa_{\varepsilon}) \frac{|u_{n,e}(v+he_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh \\ \leq (1+\kappa_{\varepsilon}) \int_{\mathbb{R}} \frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh \in L^1(\mathbb{R}^2) $$

from the coercivity inequality, and

$$ \int_{\mathbb{R}} (v_{2,n,e,\delta'}^2(v) + \kappa_{\varepsilon}) \frac{|u_{n,e,\delta}(v+\bar{h}e_1)-2u_{n,e,\delta}(v)+u_{n,e,\delta}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh \\ \leq (1+\kappa_{\varepsilon}) \int_{\mathbb{R}} \frac{|u_{n,e,\delta}(v+\bar{h}e_1)-2u_{n,e,\delta}(v)+u_{n,e,\delta}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh $$

$\in L^1(\mathbb{R}^2)$ from [18, Theorem 9]. We can thus apply the dominated convergence theorem and get

$$ \lim_{\delta' \to 0} \int_{\mathbb{R}^2} (v_{2,n,e,\delta'}(x)^2 + \kappa_{\varepsilon}) \int_{\mathbb{R}} \frac{|u_{n,e,\delta}(x+\bar{h}e_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh dx \\ \leq \lim_{\delta' \to 0} \left( \int_{\mathbb{R}^2} \int_{\mathbb{R}} (v_{2,n,e,\delta'}^2(v) + \kappa_{\varepsilon}) \frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh dv \right. \\ \left. + 2(C_1+1)\delta C \right), $$

$$ \int_{\mathbb{R}^2} (v_{2,n,e}(x)^2 + \kappa_{\varepsilon}) \int_{\mathbb{R}} \frac{|u_{n,e,\delta}(x+\bar{h}e_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh dx, \\ \leq \int_{\mathbb{R}^2} \int_{\mathbb{R}} (v_{2,n,e}(v) + \kappa_{\varepsilon}) \frac{|u_{n,e}(v+\bar{h}e_1)-2u_{n,e}(v)+u_{n,e}(v-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh dv \\ + 2(C_1+1)\delta C. $$

Let us first show that $\lim_{n \to +\infty} \int_{\mathbb{R}^2} \int_{\mathbb{R}} (v_{2,n,e}^2(x) + \kappa_{\varepsilon}) \frac{|u_{n,e,\delta}(x+\bar{h}e_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-\bar{h}e_1)|^2}{|\bar{h}|^4} \rho_n(h) dh dx = \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x)+\kappa_{\varepsilon})|D^{(2,0)}\bar{u}_{e,\delta}(x)|^2 dx = \lim_{n \to +\infty} \int_{\mathbb{R}^2} (v_{2,n,e}^2(x)+\kappa_{\varepsilon})|D^{(2,0)}\bar{u}_{e,\delta}(x)|^2 dx$, this last equality being obvious. We have

$$ \int_{\mathbb{R}^2} \int_{\mathbb{R}} (v_{2,n,e}^2(x) + \kappa_{\varepsilon}) \frac{|u_{n,e,\delta}(x+\bar{h}e_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-\bar{h}e_1)-\bar{h}^2 D^{(2,0)} \bar{u}_{e,\delta}(x)|^2}{|\bar{h}|^4} \rho_n(h) dh dx \\ \leq 2 \int_{\mathbb{R}^2} \int_0^R \left| \int_0^1 D^{(2,0)} u_{n,e,\delta}(x+(t+s-1)\bar{h}e_1) - D^{(2,0)} \bar{u}_{e,\delta}(x) ds dt \right|^2 \\ \rho_n(h) dh dx \\ + 2 \int_{\mathbb{R}^2} \int_R^{+\infty} (v_{2,n,e}^2(x) + \kappa_{\varepsilon}) \frac{|u_{n,e,\delta}(x+\bar{h}e_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-\bar{h}e_1)-\bar{h}^2 D^{(2,0)} \bar{u}_{e,\delta}(x)|^2}{|\bar{h}|^4} \rho_n(h) dh dx. $$

Using Jensen's inequality, we get

$$ \int_{\mathbb{R}^2} \int_0^R \left| \int_0^1 \int_0^1 D^{(2,0)} u_{n,e,\delta}(x+(t+s-1)\bar{h}e_1) - D^{(2,0)} \bar{u}_{e,\delta}(x) ds dt \right|^2 \\ \rho_n(h) dh dx \\ \leq \int_{\mathbb{R}^2} \int_0^R \int_0^1 |D^{(2,0)} u_{n,e,\delta}(x+(t+s-1)\bar{h}e_1) - D^{(2,0)} \bar{u}_{e,\delta}(x)|^2 ds dt \\ \rho_n(h) dh dx, \\ \leq 2 \int_{\mathbb{R}^2} \int_0^R \int_0^1 \int_0^1 |D^{(2,0)} u_{n,e,\delta}(x+(t+s-1)\bar{h}e_1) \\ - D^{(2,0)} \bar{u}_{e,\delta}(x+(t+s-1)\bar{h}e_1)|^2 ds dt \rho_n(h) dh dx \\ + 2 \int_{\mathbb{R}^2} \int_0^R \int_0^1 |D^{(2,0)} \bar{u}_{e,\delta}(x+(t+s-1)\bar{h}e_1) - D^{(2,0)} \bar{u}_{e,\delta}(x)|^2 ds dt \\ \rho_n(h) dh dx. $$

We then use Fubini's theorem, Jensen's inequality, and a Taylor's development, and get

$$ \leq 2 \int_0^R \int_0^1 \int_{\mathbb{R}^2} |D^{(2,0)} u_{n,e,\delta}(x+(t+s-1)\bar{h}e_1) - D^{(2,0)} \bar{u}_{e,\delta}(x+(t+s-1)\bar{h}e_1)|^2 ds dt \rho_n(h) dh dx \\ + 2 \int_{\mathbb{R}^2} \int_0^R \int_0^1 \int_0^1 |D^{(3,0)} \bar{u}_{e,\delta}(x+k(t+s-1)\bar{h}e_1)|^2 \\ |\bar{h}|^2 |t+s-1|^2 dk ds dt \rho_n(h) dh dx, \\ \leq 2 \|D^{(2,0)} u_{n,e,\delta} - D^{(2,0)} \bar{u}_{e,\delta}\|_{L^2(\mathbb{R}^2)}^2 \\ + \frac{1}{6} \|D^{(3,0)} \bar{u}_{e,\delta}\|_{L^2(\mathbb{R}^2)}^2 \int_0^R h^2 \rho_n(h) dh, \\ \leq 2 \|D^{(2,0)} u_{n,e,\delta} - D^{(2,0)} \bar{u}_{e,\delta}\|_{L^2(\mathbb{R}^2)}^2 \\ + \frac{R^2}{6} \|D^{(3,0)} \bar{u}_{e,\delta}\|_{L^2(\mathbb{R}^2)}^2. $$

As $u_{n,e} \xrightarrow{n \to +\infty} \bar{u}_e$ in $L^2(\mathbb{R}^2)$ and $D^{(2,0)}u_{n,e,\delta} = D^{(2,0)}\eta_\delta * u_{n,e}$, then $\|D_{n,e,\delta}^{(2,0)} - D^{(2,0)}\bar{u}_{e,\delta}\|_{L^2(\mathbb{R}^2)} \xrightarrow{n \to +\infty} 0$. Also by letting $R$ tend to 0, $\lim_{n \to +\infty} \int_{\mathbb{R}^2} \int_0^R \left| \int_0^1 \int_0^1 D^{(2,0)} u_{n,e,\delta}(x+(t+s-1)\bar{h}e_1) - D^{(2,0)} \bar{u}_{e,\delta}(x) ds dt \right|^2 \rho_n(h) dh dx = 0$. Fur-
thermore,

$$
\begin{align*}
& \int_{\mathbb{R}^2} \int_{R}^{+\infty} (v_{2,n,e}^2(x) + \kappa_{\varepsilon}) \\
& \qquad \frac{|u_{n,e,\delta}(x+he_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-he_1)-h^2D^{(2,0)}\bar{u}_{e,\delta}(x)|^2}{|h|^4} \rho_n(h) \, dh \, dx \\
& \le 2 \int_{\mathbb{R}^2} \int_{R}^{+\infty} (1 + \kappa_{\varepsilon}) \\
& \qquad \frac{|u_{n,e,\delta}(x+he_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-he_1)|^2}{|h|^4} \rho_n(h) \, dh \, dx \\
& \quad + 2 \int_{\mathbb{R}^2} \int_{R}^{+\infty} (1 + \kappa_{\varepsilon}) |D^{(2,0)}\bar{u}_{e,\delta}(x)|^2 \rho_n(h) \, dh \, dx, \\
& \le 2 \int_{\mathbb{R}^2} \int_{R}^{+\infty} (1 + \kappa_{\varepsilon}) \frac{1}{|h|^4} \\
& \qquad \left( |u_{n,e,\delta}(x+he_1) - \bar{u}_{e,\delta}(x+he_1) - 2u_{n,e,\delta}(x) + 2\bar{u}_{e,\delta}(x) + u_{n,e,\delta}(x-he_1) - \bar{u}_{e,\delta}(x-he_1) + \bar{u}_{e,\delta}(x+he_1) - 2\bar{u}_{e,\delta}(x) + \bar{u}_{e,\delta}(x-he_1)|^2 \right) \rho_n(h) \, dh \, dx \\
& \quad + 2(1 + \kappa_{\varepsilon}) \|D^{(2,0)}\bar{u}_{e,\delta}(x)\|_{L^2(\mathbb{R}^2)}^2 \int_{R}^{+\infty} \rho_n(h) \, dh, \\
& \le 4 \int_{\mathbb{R}^2} \int_{R}^{+\infty} (1 + \kappa_{\varepsilon}) \frac{1}{|h|^4} \\
& \qquad \left( |u_{n,e,\delta}(x+he_1) - \bar{u}_{e,\delta}(x+he_1)|^2 + 4|u_{n,e,\delta}(x) - \bar{u}_{e,\delta}(x)|^2 + 4|\bar{u}_{e,\delta}(x+he_1) - 2\bar{u}_{e,\delta}(x) + \bar{u}_{e,\delta}(x-he_1)|^2 + 4|u_{n,e,\delta}(x-he_1) - \bar{u}_{e,\delta}(x-he_1)|^2 \right) \rho_n(h) \, dh \, dx \\
& \quad + 2(1 + \kappa_{\varepsilon}) \|D^{(2,0)}\bar{u}_{e,\delta}(x)\|_{L^2(\mathbb{R}^2)}^2 \int_{R}^{+\infty} \rho_n(h) \, dh, \\
& \le \left( \frac{36(1 + \kappa_{\varepsilon})}{R^4} \|u_{n,e,\delta} - \bar{u}_{e,\delta}\|_{L^\infty(\mathbb{R}^2)}^2 + 16(1 + \kappa_{\varepsilon}) \|D^{(2,0)}\bar{u}_{e,\delta}\|_{L^2(\mathbb{R}^2)}^2 + 2(1 + \kappa_{\varepsilon}) \|D^{(2,0)}\bar{u}_{e,\delta}(x)\|_{L^2(\mathbb{R}^2)}^2 \right) \int_{R}^{+\infty} \rho_n(h) \, dh.
\end{align*}
$$

Since $\int_{R}^{+\infty} \rho_n(h) dh \xrightarrow{n \to +\infty} 0$, then

$$
\int_{\mathbb{R}^2} \int_{R}^{+\infty} (v_{2,n,e}^2(x) + \kappa_{\varepsilon}) \frac{|u_{n,e,\delta}(x+he_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-he_1)-h^2D^{(2,0)}\bar{u}_{e,\delta}(x)|^2}{|h|^4} \rho_n(h) \, dh \, dx \xrightarrow{n \to +\infty} 0.
$$

This result gives that

$$
\lim_{n \to +\infty} \left( \int_{\mathbb{R}^2} \int_{\mathbb{R}} (v_{2,n,e}^2(x) + \kappa_{\varepsilon}) \frac{|u_{n,e,\delta}(x+he_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-he_1)|^2}{|h|^4} \rho_n(h) \, dh \, dx \right)^{1/2} \\
= \lim_{n \to +\infty} \left( \int_{\mathbb{R}^2} (v_{2,n,e}^2(x) + \kappa_{\varepsilon}) |D^{(2,0)} \bar{u}_{e,\delta}(x)|^2 \, dx \right)^{1/2},
$$

or equivalently, the result claimed. Eventually, we have

$$
\begin{align*}
& \int_{\mathbb{R}^2} \int_{\mathbb{R}} (\bar{v}_{2,e}^2(x) + \kappa_{\varepsilon}) |D^{(2,0)} \bar{u}_{e,\delta}(x)|^2 \rho_n(h) \, dh \, dx \\
&= \liminf_{n \to +\infty} \int_{\mathbb{R}^2} \int_{\mathbb{R}} ((v_{2,n,e}^2(x) + \kappa_{\varepsilon}) \\
& \qquad \frac{|u_{n,e,\delta}(x+he_1)-2u_{n,e,\delta}(x)+u_{n,e,\delta}(x-he_1)|^2}{|h|^4} \rho_n(h) \, dh \, dx, \\
& \le \liminf_{n \to +\infty} \int_{\mathbb{R}^2} \int_{\mathbb{R}} ((v_{2,n,e}^2(x) + \kappa_{\varepsilon}) \\
& \qquad \frac{|u_{n,e}(x+he_1)-2u_{n,e}(x)+u_{n,e}(x-he_1)|^2}{|h|^4} \rho_n(h) \, dh \, dx \\
& \quad + 2C(1+C_1)\delta.
\end{align*}
$$

Since $\bar{u}_e \in W^{2,2}(\mathbb{R}^2)$ then $D^{(2,0)}\bar{u}_{e,\delta} \xrightarrow{\delta \to 0} D^{(2,0)}\bar{u}_e$ in $L^2(\mathbb{R}^2)$ and we deduce by letting $\delta$ tend to 0 that

$$
\begin{align*}
& \int_{\mathbb{R}^2} (\bar{v}_{2,e}^2(x) + \kappa_{\varepsilon}) |D^{(2,0)} \bar{u}_e(x)|^2 \, dx \\
&\le \liminf_{n \to +\infty} \int_{\mathbb{R}^2} \int_{\mathbb{R}} ((v_{2,n,e}^2(x) + \kappa_{\varepsilon}) \\
&\qquad \frac{|u_{n,e}(x+he_1)-2u_{n,e}(x)+u_{n,e}(x-he_1)|^2}{|h|^4} \rho_n(h) \, dh \, dx.
\end{align*}
$$

We finally get $\lim_{n \to +\infty} \bar{F}_{n,\varepsilon}(u_n, g_n, Q_n, v_{1,n}, v_{2,n})$
$$= \bar{F}_{\varepsilon}(\bar{u}, \bar{g}, \bar{Q}, \bar{v}_1, \bar{v}_2) \text{ and } \forall (u, g, Q, v_1, v_2) \in W^{2,2}(\Omega)$$
$$\cap W_0^{1,2}(\Omega) \times H(\text{div}) \times \{Q \in W^{1,\infty}(\Omega) \mid \int_{\Omega} Q \, dx = 0\}$$
$$\times W^{1,2}(\Omega, [0, 1]) \times W^{1,2}(\Omega, [0, 1]),$$
$$\bar{F}_{n,\varepsilon}(u_n, g_n, Q_n, v_{1,n}, v_{2,n}) \le \bar{F}_{\bar{n},\varepsilon}(u, g, Q, v_1, v_2), \text{ and by}$$
$$\text{letting } n \text{ tend to infinity } \bar{F}_{\varepsilon}(\bar{u}, \bar{g}, \bar{Q}, \bar{v}_1, \bar{v}_2) \le \bar{F}_{\varepsilon}(u, g, Q, v_1, v_2)$$
$$\text{thanks to Theorem 1. This concludes the proof.}$$

Appendix B - Proof of Theorem 4: Comparison principle

We follow the arguments of [15] for parabolic equations and
[23]. We first remark that for any $\lambda > 0$, $\tilde{q}^{-} = q^{-} - \frac{\lambda}{T-t}$ is
also a sub-solution of (EE) and satisfies:

$$
\tilde{q}_t^{-} + G_{*}(x, t, D\tilde{q}^{-}, D^2\tilde{q}^{-}) \leq - \frac{\lambda}{(T-t)^2} \leq - \frac{\lambda}{T^2},
$$

since $\tilde{q}_t^- = q_t^- - \frac{\lambda}{(T-t)^2}$, $D\tilde{q}^- = Dq^-$ and $D^2\tilde{q}^- = D^2q^-$.
As $q^- \le q^+$ comes from $\tilde{q}^- \le q^+$ in the limit when $\lambda$ tends to
0, it is sufficient to prove the comparison under the additional
assumption:

$$
\left\{
\begin{array}{ll}
(i) & q_t^- + G_*(x, t, Dq^-, D^2q^-) \leq -\frac{\lambda}{T^2} \\
(ii) & \lim_{t \to T} q^-(x, t) = -\infty
\end{array}
\right.
$$

Let us set $M = \sup_{(x,t) \in \mathbb{R}^2 \times [0,T]} q^-(x,t) - q^+(x,t)$. We want to prove that $M \le 0$. To do so, we argue by contradiction and assume that $M > 0$ so there exists $(x^*, t^*) \in \mathbb{R}^2 \times [0,T)$ such that $q^-(x^*, t^*) - q^+(x^*, t^*) > 0$.

We introduce $M_0$ defined by $M_0 = \sup_{(x,y,t) \in \mathbb{R}^2 \times \mathbb{R}^2 \times [0,T]} \{f(x,y,t) = q^-(x,t) - q^+(y,t) - \frac{1}{4\epsilon}|x-y|^4 - \frac{\alpha}{2}(|x|^2 + |y|^2)\}, \epsilon > 0$. We notice that $M_0 \ge q^-(x^*, t^*) - q^+(x^*, t^*) - \alpha|x^*|^2 > 0$ and so $M_0 > 0$ for $\alpha$ small enough. This supremum is reached owing to the term $\frac{\alpha}{2}(|x|^2 + |y|^2)$, $M_{q^-}$ the bound above of $q^-$ and $-m_{q^+}$ the bound above of $-q^+$, the upper semicontinuity of $f$, the fact that $f$ is proper and $\lim_{t \to T, |x| \to +\infty, |y| \to +\infty} f(x, y, t) = -\infty$ since $f(x, y, t) \le$
$M_{q^{-}} - m_{q^{+}} - \frac{\alpha}{2}(|x|^{2} + |y|^{2})$. We denote by $(x_0, y_0, t_0) \in \mathbb{R}^{2} \times \mathbb{R}^{2} \times [0, T)$ a maximum. Consequently, we get the following lemma:

**Lemma 2 (adapted from [23])**

Let $M' = \lim_{h \to 0} \sup_{|x-y| \le h, t \in [0,T]} (q^{-}(x,t) - q^{-}(y,t))$. Then

$$1. \lim_{\alpha \to 0} \alpha x_0 = \lim_{\alpha \to 0} \alpha y_0 = 0.$$

$$2. \lim_{\epsilon \to 0} |x_0 - y_0|^4 = 0.$$

$$3. \lim_{\epsilon \to 0} \lim_{\alpha \to 0} M_0 = M'.$$

$$4. \lim_{\epsilon \to 0} \lim_{\alpha \to 0} \frac{1}{\epsilon} |x_0 - y_0|^4 = 0.$$

$$5. \lim_{\epsilon \to 0} \lim_{\alpha \to 0} \alpha (|x_0|^2 + |y_0|^2) = 0.$$

*Proof* By boundedness of $q^{-}$ and $q^{+}$, the function $(x, y) \mapsto q^{-}(x, t) - q^{+}(y, t)$ is bounded for any $t \in [0, T)$. Besides, we assume that $M_0 > 0$ ($\alpha$ is assumed small enough). We also have $\frac{|x_0 - y_0|^4}{4\epsilon} + \frac{\alpha}{2}(|x_0|^2 + |y_0|^2) \le C$, $C$ depending on the bounds of $q^{-}$ and $q^{+}$. We deduce that $\frac{\alpha}{2}(|x_0|^2 + |y_0|^2) \le C$ and $|x_0 - y_0|^4 \le 4C\epsilon$ leading to $\lim_{\epsilon \to 0} |x_0 - y_0| = 0$ and $\lim_{\alpha \to 0} \alpha x_0 = \lim_{\alpha \to 0} \alpha y_0 = 0$.

Let us now set $M_h = \sup_{|x-y| \le h, t \in [0,T)} (q^{-}(x,t) - q^{-}(y,t))$.

Let $(x_n^h, y_n^h, t_n^h) \in \mathbb{R}^2 \times \mathbb{R}^2 \times [0, T)$ be a sequence such that $\forall n \in \mathbb{N}^*, q^{-}(x_n^h, t_n^h) - q^{+}(y_n^h, t_n^h) \ge M_h - \frac{1}{n}$ with $|x_n^h - y_n^h| \le h$. This sequence is independent of $\alpha$. We then get $M_h - \frac{1}{n} - \frac{h^4}{4\epsilon} - \frac{\alpha}{2}(|x_n^h|^2 + |y_n^h|^2) \le q^{-}(x_n^h, t_n^h) - v(y_n^h, t_n^h) - \frac{|x_n^h - y_n^h|^4}{4\epsilon} - \frac{\alpha}{2}(|x_n^h|^2 + |y_n^h|^2) \le M_0 \le q^{-}(x_0, t_0) - q^{+}(y_0, t_0)$.

Let $\alpha$ tend to 0 in what precedes, then $M_h - \frac{1}{n} - \frac{h^4}{4\epsilon} \le \liminf_{\alpha \to 0} (q^{-}(x_0, t_0) - v(y_0, t_0)) \le \limsup_{\alpha \to 0} (q^{-}(x_0, t_0) - v(y_0, t_0))$.

Now let $h$ tend to 0, we get $M' - \frac{1}{n} \le \liminf_{\alpha \to 0} (q^{-}(x_0, t_0) - q^{+}(y_0, t_0))$

$\le \limsup_{\alpha \to 0} (q^{-}(x_0, t_0) - q^{+}(y_0, t_0))$. Finally, we let $\epsilon$ tend to 0

and obtain:

$$M' - \frac{1}{n} \le \liminf_{\epsilon \to 0} \liminf_{\alpha \to 0} (q^{-}(x_0, t_0) - q^{+}(y_0, t_0))$$

$$\le \liminf_{\epsilon \to 0} \limsup_{\alpha \to 0} (q^{-}(x_0, t_0) - q^{+}(y_0, t_0))$$

$$\le \limsup_{\epsilon \to 0} \limsup_{\alpha \to 0} (q^{-}(x_0, t_0) - q^{+}(y_0, t_0))$$

$$\le \limsup_{\epsilon \to 0} \limsup_{\alpha \to 0} \sup_{|x-y| \le C \frac{1}{4}, t \in [0,T)} (q^{-}(x,t) - q^{+}(y,t))$$

$$\le \limsup_{h \to 0} \sup_{|x-y| \le h, t \in [0,T)} (q^{-}(x,t) - q^{+}(y,t)) = M'$$

So, $\lim_{\epsilon \to 0} \lim_{\alpha \to 0} (q^{-}(x_0, t_0) - q^{+}(y_0, t_0)) = M'$. In the same way, we get $\lim_{\epsilon \to 0} \lim_{\alpha \to 0} M_0 = M'$.

Thus $\lim_{\epsilon \to 0} \lim_{\alpha \to 0} (\frac{|x_0 - y_0|^4}{4\epsilon} + \frac{\alpha}{2}(|x_0|^2 + |y_0|^2)) = \lim_{\epsilon \to 0} \lim_{\alpha \to 0} (q^{-}(x_0, t_0) - q^{+}(y_0, t_0) - M_0) = 0$.

$$\text{So} \left\{ \begin{array}{l} \lim_{\epsilon \to 0} \lim_{\alpha \to 0} \frac{|x_0 - y_0|^4}{4\epsilon} = 0 \\ \lim_{\epsilon \to 0} \lim_{\alpha \to 0} \frac{\alpha}{2} (|x_0|^2 + |y_0|^2) = 0 \end{array} \right.$$

We then distinguish two cases.

1. $\forall \epsilon > 0, \exists \alpha \in (0, \epsilon)$ such that $t_0 = 0$. Then there exists $\epsilon_n \underset{n \to +\infty}{\longrightarrow} 0$ and $\alpha_n \underset{n \to +\infty}{\longrightarrow} 0$ with $t_n = 0$ and $0 < M_0 \le q^{-}(x_0, n, 0) - q^{+}(y_0, n, 0) \le Q_0(x_0, n) - Q_0(y_0, n) \le B_0 |x_0 - y_0, n|$ where $B_0$ is the Lipschitz constant of $Q_0$. We obtain a contradiction since $\lim_{n \to +\infty} |x_0, n - y_0, n| = 0$.

2. $\exists \epsilon > 0$ such that $\forall \alpha \in (0, \epsilon), t_0 > 0$. We can choose $\epsilon$ small enough (otherwise we use the first argument), i.e., such that $\frac{|x_0 - y_0|^4}{\epsilon} \le \frac{\lambda}{2T^2 \gamma \kappa_g}$. We consider $\tilde{q}^{-}(x,t) = q^{-}(x,t) - \frac{\alpha}{2}|x|^2$ and $\tilde{q}^{+}(x,t) = q^{+}(x,t) + \frac{\alpha}{2}|x|^2$ and so $M_0 = \sup_{(x,y,t) \in \mathbb{R}^2 \times \mathbb{R}^2 \times [0,T)} \{\tilde{q}^{-}(x,t) - \tilde{q}^{+}(y,t) - \frac{|x-y|^4}{4\epsilon}\}$.

Let us take the test function $\psi(x, y, t) = \frac{|x-y|^4}{4\epsilon}$ and set $p_0 = x_0 - y_0$.

By using the parabolic version of Ishii's lemma ([15, Theorem 8.3.]) with the same notations, it comes that $\tau = 0$, $p_1 = \frac{|p_0|^2 p_0}{\epsilon} = p_2$, $A = \frac{2}{\epsilon} |p_0|^2 \begin{pmatrix} Z & -Z \\ -Z & Z \end{pmatrix}$ with $Z = \frac{1}{2} + \frac{p_0 \otimes p_0}{|p_0|^2}$ and for each $\beta > 0$ such that $\beta A < I$, there exist $(X, Y) \in (S^2)^2$ and two reals $\tau_1$ and $\tau_2$ such that

$$\tau_1 - \tau_2 = 0,$$

$$(\tau_1, \frac{|p_0|^2 p_0}{\epsilon}, X) \in \bar{P}^{+} \tilde{u}(x_0, t_0),$$

$$(\tau_2, \frac{|p_0|^2 p_0}{\epsilon}, Y) \in \bar{P}^{-} \tilde{v}(y_0, t_0),$$

$$\frac{-1}{\beta} \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix} \le \begin{pmatrix} X & 0 \\ 0 & -Y \end{pmatrix} \le (I - \beta A)^{-1} A.$$

The last inequality implies that $X \preceq Y$. We then use the following lemma.

**Lemma 3 (adapted from [23])** We have the following estimations on the matrix A:

$$\frac{1}{\|A\|} A \preceq I,$$

$$\|A\| \le \frac{6|p_0|^2}{\epsilon},$$

$$\text{If } \delta = \frac{1}{2\|A\|}, \text{ then } (I - \delta A)^{-1} A \preceq 2 \|A\| I \le \frac{12}{\epsilon} |p_0|^2 I,$$

$$\text{with } \|A\| = \sup_{\xi \in \mathbb{R}^4} \frac{|A\xi \cdot \xi|}{\xi \cdot \xi}, \text{ the denoting the inner product in } \mathbb{R}^4.$$

*Proof* By definition of $\|A\|$, we have $\frac{A\xi \cdot \xi}{\xi \cdot \xi} \le \|A\|$ for any $\xi \in \mathbb{R}^4 \setminus \{0, \mathbb{R}^4\}$ and so $\frac{A\xi \cdot \xi}{\|A\|} \preceq I\xi \cdot \xi$, which gives the first result of the lemma.

Let $\xi = \begin{pmatrix} \xi_1 \\ \xi_2 \end{pmatrix} \in \mathbb{R}^4$, $\xi_i \in \mathbb{R}^2$, $i = \{1, 2\}$. We have:

$$A\xi \cdot \xi = \frac{2}{\epsilon} |p_0|^2 Z(\xi_1 - \xi_2) \cdot (\xi_1 - \xi_2),$$

$$\le \frac{4}{\epsilon} |p_0|^2 \|Z\| \|\xi\|^2.$$

So it suffices to show that $\|Z\| \le \frac{3}{2}$.

$$Z\xi_1 \cdot \xi_1 = \frac{|\xi_1|^2}{2} + \frac{p_0 \otimes p_0}{|p_0|^2} \xi_1 \cdot \xi_1$$

$$\le \frac{|\xi_1|^2}{2} + |\xi_1|^2,$$

$$\le \frac{3}{2} |\xi_1|^2.$$

For the last point, it suffices to notice that if $B \ge 0$ and $C \ge 0$ with $B, C \in S^2$ such that $BC = CB$, then $CB \ge 0$. Indeed, $CB\xi \cdot \xi = CB \frac{1}{2}\xi \cdot B \frac{1}{2}\xi \ge 0$ since $C$ is non-negative and $B$ is symmetric non-negative. So, if $B \succeq C$ and $D \succeq 0$ with $D(B - C) = (B - C)D$ then $DB \succeq DC$. It thus suffices to prove that $A \preceq 2\|A\|(I - \delta A) = 2\|A\|I - A$, i.e., $A \preceq \|A\|I$ which is true by definition of the norm.
By using this lemma with $\beta = \frac{1}{2\|A\|}$, the inequality becomes
$$ \frac{-12|p_0|^2}{\epsilon} \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix} \preceq \begin{pmatrix} X & 0 \\ 0 & -Y \end{pmatrix} \preceq \frac{12|p_0|^2}{\epsilon} \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix}. $$
Thus $X$ and $Y$ are bounded independently of $\alpha$ (since $p_0$ is bounded independently of $\alpha$ having $\frac{|x_0-y_0|^4}{\epsilon} \le \frac{\lambda}{2T^2\gamma\kappa_g}$). In particular, we have $\frac{-12|p_0|^2}{\epsilon} I \preceq X \preceq \frac{12|p_0|^2}{\epsilon} I$. So there exists a sequence $(\alpha_n)$ such that $\alpha_n \xrightarrow{n \to +\infty} 0$, $t_0 \xrightarrow{n \to +\infty} 0$, $p_0 \xrightarrow{n \to +\infty} p_\infty$, $X \xrightarrow{n \to +\infty} X_\infty$ and $Y \xrightarrow{n \to +\infty} Y_\infty$. Furthermore, since $q^-$ is a sub-solution, $q^+$ is a supersolution and using the additional assumption, we get:

$$ \begin{aligned} \tau_1 + G_*(x_0, t_0, \frac{|p_0|^2 p_0}{\epsilon} + \alpha_n x_0, X + \alpha_n I) &\le \frac{-\lambda}{T^2}, \\ \tau_2 + G^*(y_0, t_0, \frac{|p_0|^2 p_0}{\epsilon} - \alpha_n y_0, Y - \alpha_n I) &\ge 0. \end{aligned} $$

Then using the matrix inequality $X \preceq Y$ and the ellipticity of the equation, we obtain $\tau_2 + G^*(y_0, t_0, \frac{|p_0|^2 p_0}{\epsilon} - \alpha_n y_0, X - \alpha_n I) \ge 0$. Subtracting these two inequalities yields successively

$$ \begin{aligned} & \frac{\lambda}{T^2} + G_* \left( x_0, t_0, \frac{|p_0|^2 p_0}{\epsilon} + \alpha_n x_0, X + \alpha_n I \right) \\ & \le G^* \left( y_0, t_0, \frac{|p_0|^2 p_0}{\epsilon} - \alpha_n y_0, X - \alpha_n I \right), \\ & \frac{\lambda}{T^2} \le -\gamma \operatorname{divg}(x_0, t_0) + \gamma \operatorname{divg}(y_0, t_0) - E(X + \alpha_n I) \\ & - F_* \left( \frac{|p_0|^2 p_0}{\epsilon} + \alpha_n x_0, X + \alpha_n I \right) + E(X - \alpha_n I) \\ & + F^* \left( \frac{|p_0|^2 p_0}{\epsilon} - \alpha_n y_0, X - \alpha_n I \right). \end{aligned} $$

Let $n$ tend to infinity. As $\operatorname{divg}$ is Lipschitz continuous uniformly in time with $\kappa_g$ the Lipschitz constant, we get:

$$ \begin{aligned} & \gamma \kappa_g |p_\infty| - F_* \left( \frac{|p_\infty|^2 p_\infty}{\epsilon}, X_\infty \right) + F^* \left( \frac{|p_\infty|^2 p_\infty}{\epsilon}, X_\infty \right) \\ & \ge \frac{\lambda}{T^2}. \end{aligned} $$

We distinguish here two cases:

- First case: $p_\infty = 0$. From the matrix inequality
$$ \frac{-12|p_0|^2}{\epsilon} \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix} \preceq \begin{pmatrix} X & 0 \\ 0 & -Y \end{pmatrix} \preceq \frac{12|p_0|^2}{\epsilon} \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix}, $$
it yields $X_\infty = Y_\infty = 0$. Having $F^*(0,0) = F^*(0,0) = 0$ leads to $0 \ge \frac{\lambda}{T^2}$, which is absurd.

- Second case: $p_\infty \ne 0$. Since $F$ is continuous on $\mathbb{R}^2 \setminus 0_{\mathbb{R}^2} \times S^2$, then
$$ F_* \left( \frac{|p_\infty|^2 p_\infty}{\epsilon}, X_\infty \right) = F^* \left( \frac{|p_\infty|^2 p_\infty}{\epsilon}, X_\infty \right) = F \left( \frac{|p_\infty|^2 p_\infty}{\epsilon}, X_\infty \right). $$
Now, $|p_0| \le \left( \frac{\epsilon \lambda}{2T^2 \gamma \kappa_g} \right)^{1/4}$ and so $|p_\infty| \le \left( \frac{\epsilon \lambda}{2T^2 \gamma \kappa_g} \right)^{1/4}$. Letting $\epsilon$ tend to 0 raises again a contradiction.

## Appendix C - Proof of Theorem 7: Lipschitz regularity in space

We follow the arguments of [23, Lemma 4.15]. We have proved that $Q$ is bounded and continuous on $\mathbb{R}^2 \times [0, T)$. Let us set $\Phi^\epsilon(x, y, t) = B(t)(|x-y|^2 + \epsilon^2)^{\frac{1}{2}}$. We aim to show that $Q(x, t) - Q(y, t) \le \Phi^\epsilon(x, y, t)$ for any $(x, y, t) \in \mathbb{R}^2 \times \mathbb{R}^2 \times [0, T)$.

We introduce $M = \sup_{(x,y) \in \mathbb{R}^2 \times \mathbb{R}^2, t \in [0, T)} (Q(x,t) - Q(y,t) - \Phi^\epsilon(x, y, t))$.
We assume that $M > 0$. We denote by $\bar{M} = \sup_{(x,y) \in \mathbb{R}^2 \times \mathbb{R}^2, t \in [0, T)} \{Q(x,t) - Q(y,t) - \frac{\alpha}{2}(|x|^2 + |y|^2) - \frac{\lambda}{T-t} - \Phi^\epsilon(x, y, t)\}$. For $\lambda$ and $\alpha$ small enough, we have $\bar{M} > 0$. As $Q$ is bounded and continuous with $|x| \to +\infty, |y| \to +\infty, t \to T$, $Q(x,t) - Q(y,t) - \frac{\lambda}{T-t} - \Phi^\epsilon(x, y, t) = -\infty$, the maximum is attained at $(\bar{x}, \bar{y}, \bar{t})$ with $\bar{x} \ne \bar{y}$ and $\frac{\alpha}{2}(|\bar{x}|^2 + |\bar{y}|^2) \le C$ owing to $\bar{M} > 0$ and $Q$ bounded. Therefore, $\lim_{\alpha \to 0} \alpha|\bar{x}| = \lim_{\alpha \to 0} \alpha|\bar{y}| = 0$.

We now prove that $\bar{t} > 0$ by contradiction. Let us assume that $\bar{t} = 0$. Then $Q(\bar{x}, 0) - Q(\bar{y}, 0) - \Phi^\epsilon(\bar{x}, \bar{y}, 0) > 0$, that is to say, $Q_0(\bar{x}) - Q_0(\bar{y}) > B(0)(|\bar{x} - \bar{y}|^2 + \epsilon^2)^{\frac{1}{2}} > B_0|\bar{x} - \bar{y}|$, which is absurd since $\|\nabla Q_0\|_{L^\infty(\mathbb{R}^2)} \le B_0$.

We set $\bar{p} = D_x \Phi^\epsilon(\bar{x}, \bar{y}, \bar{t}) = -D_y \Phi^\epsilon(\bar{x}, \bar{y}, \bar{t}) = B(\bar{t})(\bar{x} - \bar{y})(|\bar{x} - \bar{y}|^2 + \epsilon^2)^{-\frac{1}{2}} \ne 0$, $Z = D_x^2 \Phi^\epsilon(\bar{x}, \bar{y}, \bar{t}) = D_y^2 \Phi^\epsilon(\bar{x}, \bar{y}, \bar{t}) = B(\bar{t})\left(|(\bar{x} - \bar{y}|^2 + \epsilon^2)^{-\frac{1}{2}} I - \left(|\bar{x} - \bar{y}|^2 + \epsilon^2\right)^{-\frac{3}{2}} (\bar{x} - \bar{y}) \otimes (\bar{x} - \bar{y})\right)$ and $A = \begin{pmatrix} Z & -Z \\ -Z & Z \end{pmatrix}$. Then by the parabolic version of Ishii's lemma ([15, Theorem 8.3.]) applied to $\tilde{u}(x,t) = Q(x,t) - \frac{\alpha}{2}|x|^2$, $\tilde{v}(x,t) = Q(x,t) + \frac{\alpha}{2}|x|^2$, and $\Phi(x, y, t) = \Phi^\epsilon(x, y, t) + \frac{\lambda}{T-t}$, for every $\beta$ such that $\beta A < I$, there exist $\tau_1 \in \mathbb{R}$, $\tau_2 \in \mathbb{R}$ and $X \in S^2$, $Y \in S^2$ such that

$$ \begin{aligned} \tau_1 - \tau_2 &= \frac{\lambda}{(T-t)^2} + B'(\bar{t})(|\bar{x} - \bar{y}|^2 + \epsilon^2)^{\frac{1}{2}} = \frac{\lambda}{(T-t)^2} \\ &\quad + \gamma \kappa_g (|\bar{x} - \bar{y}|^2 + \epsilon^2)^{\frac{1}{2}}, \\ (\tau_1, \bar{p} + \alpha \bar{x}, X + \alpha I) &\in \bar{P}^+ Q(\bar{x}, \bar{t}), \\ (\tau_2, \bar{p} - \alpha \bar{y}, Y - \alpha I) &\in \bar{P}^- Q(\bar{y}, \bar{t}), \\ \frac{-1}{\beta} \begin{pmatrix} I & 0 \\ 0 & I \end{pmatrix} &\le \begin{pmatrix} X & 0 \\ 0 & -Y \end{pmatrix} \le (I - \beta A)^{-1} A. \end{aligned} $$

Thus we have

$$ \begin{aligned} \tau_1 + \operatorname{divg}(\bar{x}, \bar{t}) + E(X + \alpha I) + F_* (\bar{p} + \alpha \bar{x}, X + \alpha I) &\le 0, \\ \tau_2 + \operatorname{divg}(\bar{y}, \bar{t}) + E(Y - \alpha I) + F^* (\bar{p} - \alpha \bar{y}, Y - \alpha I) &\ge 0. \end{aligned} $$

As the matrix inequality implies $X \preceq Y$ and using the ellipticity of the equation, we deduce that

$$ \tau_2 + \operatorname{divg}(\bar{y}, \bar{t}) + E(X - \alpha I) + F^*(\bar{p} - \alpha \bar{y}, X - \alpha I) \ge 0. $$

Subtracting both inequalities and using again the Lipschitz continuity of $\operatorname{divg}$, we get successively

$$ \begin{aligned} & \frac{\lambda}{(T-t)^2} + \gamma \kappa_g (|\bar{x} - \bar{y}|^2 + \epsilon^2)^{\frac{1}{2}} - \gamma \kappa_g |\bar{x} - \bar{y}| \\ & + E(X + \alpha I) - E(X - \alpha I) \\ & + F_* (\bar{p} + \alpha \bar{x}, X + \alpha I) - F^* (\bar{p} - \alpha \bar{y}, X - \alpha I) \le 0, \\ & \frac{\lambda}{T^2} + E(X + \alpha I) - E(X - \alpha I) \\ & + F_* (\bar{p} + \alpha \bar{x}, X + \alpha I) - F^* (\bar{p} - \alpha \bar{y}, X - \alpha I) \le 0. \end{aligned} $$

Letting $\alpha$ tend to 0 ($\bar{p}$ and $X$ are bounded independently of $\alpha$ so we can extract convergent subsequences) and using similar arguments to those in [23] to handle $F_*(\bar{p} + \alpha \bar{x}, X + \alpha I) - F^*(\bar{p} - \alpha \bar{y}, X - \alpha I)$ yields $\frac{\lambda}{T^2} \le 0$, which is absurd.

Then $Q(x, t) - Q(y, t) \le \Phi^\epsilon(x, y, t)$. Letting $\epsilon$ tend to 0, we
get $Q(x,t) - Q(y,t) \le B(t)|x-y|$. By exchanging $x$ and $y$ in what precedes, we get $|Q(x,t) - Q(y,t)| \le B(t)|x-y|$, which concludes the proof.

## Appendix D - Proof of Theorem 8: Uniform continuity in time

We follow again the arguments of [23, Lemma 4.15]. Let us set $\delta > 0$. For any $(x, t) \in \mathbb{R}^2 \times (0, T)$ such that $t + \delta \le T$, we set $q_1(x, t) = Q(x, t + \delta)$. Then $q_1$ is a sub-solution of $w_t - \omega_{\text{div}\mathbf{g}}(\delta) + \text{div}\mathbf{g}(x, t) + E(D^2w) + F(Dw, D^2w) = 0$ on $\mathbb{R}^2 \times (0, T - \delta)$ with $w(x, 0) = Q(x, \delta)$ on $\mathbb{R}^2$ in the sense of viscosity solutions theory. Indeed, we have $q_{1,t}(x, t) = Q_t(x, t + \delta)$, $q_{1,t} + \text{div}\mathbf{g}(x, t + \delta) + E(D^2q_1) + F(Dq_1, D^2q_1) = 0$ and $\text{div}\mathbf{g}(x, t + \delta) \ge -\omega_{\text{div}\mathbf{g}}(\delta) + \text{div}\mathbf{g}(x, t)$ leading to $q_{1,t} - \omega_{\text{div}\mathbf{g}}(\delta) + \text{div}\mathbf{g}(x, t) + E(D^2q_1) + F(Dq_1, D^2q_1) \le 0$. Besides, $\tilde{q} = Q + \sup_{x \in \mathbb{R}^2} ((Q(x, \delta) - Q_0(x))^+) + \omega_{\text{div}\mathbf{g}}(\delta)t$ ($-\omega_+ = \max\{0, \alpha\}$) is a super-solution with $q_1(x, 0) \le \tilde{q}(x, 0)$ for any $x \in \mathbb{R}^2$. Using the comparison principle, we get, $\forall (x, t) \in \mathbb{R}^2 \times (0, T - \delta)$, $Q(x, t + \delta) - Q(x, t) \le \sup_{x \in \mathbb{R}^2} ((Q(x, \delta) - Q_0(x))^+) + \omega_{\text{div}\mathbf{g}}(\delta)t$.

We will now focus on an auxiliary problem. Let us first assume that $Q_0 \in C_b^2(\mathbb{R}^2)$. We set $q^+/ = Q_0(x) + / - C_1t$ with $C_1 = \inf_{x \in \mathbb{R}^2} \{-E(D^2Q_0) - F^*(DQ_0, D^2Q_0), E(D^2Q_0) + F_*(DQ_0, D^2Q_0)\}$. We easily check that $q^+$ is a super-solution and $q^-$ a sub-solution of (SP): $u_t + E(D^2u) + F(Du, D^2u) = 0$, $u(x, 0) = Q_0(x)$. Then there exists a unique solution $q$ of this problem and by the comparison principle, the following holds: $\forall t \in [0, T)$, $\forall x \in \mathbb{R}^2$, $|q(x, t) - Q_0(x)| \le C_1t$. We then set $q_1(x, t) = q(x, t + h)$. So $q_1$ is also a solution to the problem (SP) with $q_1(x, 0) = q(x, h)$, $\forall x \in \mathbb{R}^2$ and by the comparison principle, we get $q(x, t + h) - q(x, t) \le \sup_{x \in \mathbb{R}^2} \{q(x, h) - Q_0(x)\} \le C_1h$. Similarly, we have $q(x, t) - q(x, t + h) \le \sup_{x \in \mathbb{R}^2} \{q(x, h) - Q_0(x)\} \le C_1h$ and so $|q(x, t + h) - q(x, t)| \le \sup_{x \in \mathbb{R}^2} \{q(x, h) - Q_0(x)\} \le C_1h$. We now assume that

$Q_0 \in W^{1, \infty}(\mathbb{R}^2)$. We set $Q_\epsilon^0 = Q_0 * \rho_\epsilon$, where $\rho_\epsilon$ is a regularizing sequence satisfying $\rho_\epsilon = \frac{1}{\epsilon^2}\rho(\cdot, \epsilon)$, $\rho \in C_c^\infty(\mathbb{R}^2, \mathbb{R})$, $\rho \ge 0$, $\text{supp}(\rho) \subset \bar{B}(0, 1)$, and $\int_{\mathbb{R}^2} \rho(x) dx = 1$. Then $Q_\epsilon^0 \in C_c^2(\mathbb{R}^2)$ and $\|DQ_\epsilon^0\|_{L^\infty(\mathbb{R}^2)} \le B_0$, $\|D^2Q_\epsilon^0\|_{L^\infty(\mathbb{R}^2)} \le \frac{B_0C_2}{\epsilon}$. Indeed using properties stated in [22, Theorem 1, p. 123]⁶, we have

$$ 
\begin{aligned}
|DQ_\epsilon^0(x)| &= |DQ_0 * \rho_\epsilon(x)| = \left| \int_{\mathbb{R}^2} DQ_0(x-y)\rho_\epsilon(y) dy \right|, \\
&\le \int_{\mathbb{R}^2} |DQ_0(x-y)\rho_\epsilon(y)| dy \le B_0 \int_{\mathbb{R}^2} \rho_\epsilon(y) dy, \\
&\le B_0,
\end{aligned}
$$

and

$$ 
\begin{aligned}
|D^2 Q_\epsilon^0(x)| &= |DQ_0 * D\rho_\epsilon(x)| = \left| \int_{\mathbb{R}^2} DQ_0(x-y) D\rho_\epsilon(y) dy \right|, \\
&\le \int_{\mathbb{R}^2} |DQ_0(x-y) D\rho_\epsilon(y)| dy, \\
&\le \frac{B_0}{\epsilon} \int_{\mathbb{R}^2} \frac{1}{\epsilon^2} \left| D\rho\left(\frac{y}{\epsilon}\right) \right| dy \le \frac{B_0}{\epsilon} \|D\rho\|_{L^1(\mathbb{R}^2)}.
\end{aligned}
$$

Moreover, $\|Q_0 - Q_\epsilon^0\|_{L^\infty(\mathbb{R}^2)} \le B_0\epsilon$ since

$$ 
\begin{aligned}
|Q_0(x) - Q_\epsilon^0(x)| &\le \int_{\mathbb{R}^2} |Q_0(x) - Q_0(x-y)| \rho_\epsilon(y) dy, \\
&\le B_0 \int_{\bar{B}(0, \epsilon)} |y| \rho_\epsilon(y) dy, \\
&\le B_0 \epsilon \int_{\bar{B}(0, \epsilon)} \rho_\epsilon(y) dy = B_0 \epsilon.
\end{aligned}
$$

We denote by $q_\epsilon$ the solution with initial condition $Q_\epsilon^0$. Then, by the comparison principle, $\|q_\epsilon'(., t) - q_\epsilon(., t)\|_{L^\infty(\mathbb{R}^2)} \le \|Q_\epsilon^0 - Q_\epsilon^0\|_{L^\infty(\mathbb{R}^2)}$, and so $(u_\epsilon)$ uniformly converges to $q$ the solution of (SP) with initial condition $Q_0$ since $(Q_\epsilon^0)$ converges uniformly to $Q_0$. We then have by the comparison principle that $\|q_\epsilon(., t) - q(., t)\|_{L^\infty(\mathbb{R}^2)} \le \|Q_\epsilon^0 - Q_0\|_{L^\infty(\mathbb{R}^2)}$ and deduce that

$$ 
\begin{aligned}
\|q(., t+h) - q(., t)\|_{L^\infty(\mathbb{R}^2)} &\le 2\|Q_\epsilon^0 - Q_0\|_{L^\infty(\mathbb{R}^2)} \\
&\quad + \|q_\epsilon(., t+h) - q(., t)\|_{L^\infty(\mathbb{R}^2)}, \\
&\le 2B_0\epsilon + C_1 \left( B_0, \frac{B_0C_2}{\epsilon} \right) h.
\end{aligned}
$$

By taking the minimum on $\epsilon$, we obtain the modulus of continuity of $q$, called $\omega_F$ which depends only on $B_0$. And so $q$ satisfies $Q_0(x) - \omega_F(t) \le q(x, t) \le Q_0(x) + \omega_F(t)$. We now prove that $\|Dq\|_{L^\infty(\mathbb{R}^2 \times (0, T))} \le \|DQ_0\|_{L^\infty(\mathbb{R}^2)}$. Indeed, let us consider the function $q^h(x, t) = q(x+h, t) + \|DQ_0\|_{L^\infty(\mathbb{R}^2)}|h|$. Then $q^h$ is still a solution of (SP) with initial condition $q^h(x, 0) = Q_0(x+h) + \|DQ_0\|_{L^\infty(\mathbb{R}^2)}|h| \ge Q_0(x)$. So, by the comparison principle, we have $q^h \ge q$. We deduce that for every $h \in \mathbb{R}^2$, $q(x, t) - q(x+h, t) \le \|DQ_0\|_{L^\infty(\mathbb{R}^2)}|h|$, and so $\|Dq\|_{L^\infty(\mathbb{R}^2 \times (0, T))} \le \|DQ_0\|_{L^\infty(\mathbb{R}^2)}$. We then set $Q^+(x, t) = q(x, t) + \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}t$ and come back to the initial problem. Then

$\|DQ^+\|_{L^\infty(\mathbb{R}^2 \times (0, T))} \le \|Dq\|_{L^\infty(\mathbb{R}^2 \times (0, T))} \le B_0$ and $Q^+$ is a super-solution of (EE) since $Q^+(x, 0) = q(x, 0) = Q_0(x)$ and $Q_t^+(x, t) + G^*(x, t, Q^+, D^2Q^+) = q_t(x, t) + \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))} + \|\text{div}\mathbf{g}(x, t) + E(D^2q(x, t)) + F^*(Dq(x, t), D^2q(x, t))\| = \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))} + \|\text{div}\mathbf{g}(x, t) \ge 0$ and satisfies $Q^+(x, t) \le Q_0(x) + \omega_F(t) + \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}t$. Similarly, we construct a sub-solution $Q^-$ such that $Q^-(x, t) \ge Q_0(x) - \omega_F(t) - \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}t$ by setting $Q^-(x, t) = q(x, t) - \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}t$. By applying the comparison principle, we get that

$$ 
\begin{aligned}
Q_0(x) - \omega_F(t) - \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}t &\le Q^-(x, t) \le Q(x, t) \\
&\le Q^+(x, t) \le Q_0(x) + \omega_F(t) + \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}t.
\end{aligned}
$$

Thus sup$_{x \in \mathbb{R}^2}$ $((Q(x, \delta) - Q_0(x))^+) \le \omega_F(\delta) + \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}\delta$ and so $Q(x, t + \delta) - Q(x, t) \le \omega_F(\delta) + \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}\delta + \omega_{\text{div}\mathbf{g}}(\delta)T$. Similarly, $v$ is a super-solution of $w_t + \omega_{\text{div}\mathbf{g}}(\delta) + \text{div}\mathbf{g}(x, t) + E(D^2w) + F(Dw, D^2w) = 0$ with $w(x, 0) = Q(x, \delta)$ on $\mathbb{R}^2$ and $\tilde{u} = Q(x, t) - \sup_{x \in \mathbb{R}^2} ((Q(x, \delta) - Q_0(x))^-) - \omega_{\text{div}\mathbf{g}}(\delta)t$

($\omega_-$) is a sub-solution. So, by the comparison principle, we have

$$ 
\begin{aligned}
Q(x, t) - Q(x, t + \delta) &\le \sup_{x \in \mathbb{R}^2} ((Q(x, \delta) - Q_0(x))^-) + \omega_{\text{div}\mathbf{g}}(\delta)t, \\
&\le \omega_F(\delta) + \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}\delta \\
&\quad + \omega_{\text{div}\mathbf{g}}(\delta)T.
\end{aligned}
$$

And so,

$$ 
\begin{aligned}
|Q(x, t) - Q(x, t + \delta)| &\le \omega_F(\delta) + \|\text{div}\mathbf{g}\|_{L^\infty(\mathbb{R}^2 \times [0, T))}\delta \\
&\quad + \omega_{\text{div}\mathbf{g}}(\delta)T,
\end{aligned}
$$
which achieves the proof.

## References

1. Adam, R. D., Peter, P., Weickert, J., *Denoising by in-painting*. In Lauze, F., Dong, Y., Dahl, A. B. (eds) Scale Space and Variational Methods in Computer Vision: 6th International Conference, SSVM 2017, Kolding, Denmark, June 4-8, 2017, Proceedings, pp. 121-132. Springer International Publishing (2017).

2. Ambrosio, L., Faina, L., March, R., Variational Approximation of a Second Order Free Discontinuity Problem in Computer Vision. *SIAM J. Math. Anal.* **32**(6) (2001), 1171-1197.

3. Aronsson, G., Crandall, M., Juutinen, P., A Tour of the Theory of Absolutely Minimizing Functions. *B. Am. Math. Soc.* **41** (2004), 439-505.

4. Aubert, G., Aujol, J.-F., Blanc-Féraud, L., Detecting Codimension-Two Objects in an Image with Ginzburg-Landau Models. *Int. J. Comput. Vision* **65**(1) (2005), 29-42.

5. Aubert, G., Kornprobst, P., Can the nonlocal characterization of Sobolev spaces by Bourgain et al. be useful for solving variational problems? *SIAM J. Numer. Anal.* **47**(2) (2009), 844-860.

6. Aubert, G., Drogoul, A., Topological Gradient for a Fourth Order Operator Used in Image Analysis. *ESAIM: COCV* **21**(4) (2015), 1120-1149.

7. Aujol, J.-F., Chambolle, A., Dual norms and image decomposition models. *International Journal of Computer Vision* **63**(1) (2005), 85-104.

8. Baudour, A., Aubert, G., Blanc-Féraud, L., *Detection and Completion of Filaments: A Vector Field and PDE Approach*. In: Sgallari, F., Murli, A., Paragios, N. (eds.) SSVM 2007. LNCS, vol. 4485, pp. 451-460. Springer, Heidelberg (2007).

9. Bergounioux, M., Vicente, D., Parameter Selection in a Mumford-Shah Geometrical Model for the Detection of Thin Structures. *Acta Appl. Math.* **141**(1) (2016), 17-48.

10. Blake, A., Zisserman, A., *Visual Reconstruction*. MIT Press, Cambridge. (1989).

11. Bourgain, J., Brezis, H., Mironescu, P., *Another look at Sobolev spaces*. In: Menaldi, J. L., Rofman, E., Sulem, A. (eds) Optimal Control and Partial Differential Equations, in honour of Professor Alain Bensoussan's 60th Birthday, pp. 439-455 (2001).

12. Brezis, H., *Analyse fonctionnelle*. Dunod Paris (2005).

13. Buades, A., Coll, B., Morel, J.-M., A Review of Image Denoising Algorithms, with a New One. *Multiscale Model. Simul.* **4**(2) (2005), 490-530.

14. Caselles, V., Morel, J.-M., Sbert, C., An Axiomatic Approach to Image interpolation. *IEEE T. Image Process.* **7**(3) (1998), 376-386.

15. Crandall, M., Ishii, H., Lions P.-L., User's guide to viscosity solutions of second order partial differential equations. *Bulletin of the American Mathematical Society* **27**(1) (1992), 1-67.

16. Crow, F.C., *Summed-Area Tables for Texture Mapping*. In: Proceedings of the 11th Annual Conference on Computer Graphics and Interactive Techniques. SIGGRAPH '84. ACM (1984), 207-212.

17. Dávila, J., On an open question about functions of bounded variation. *Calc. Var. Partial Differential Equations* **15**(4) (2002), 519-527.

18. Debroux, N., Le Guyader, C., Vese, L.A., A second order free discontinuity model for bituminous surfacing crack recovery and analysis of a nonlocal version of it. *Annals of Mathematical Sciences and Applications* **3**(1) (2018), 49-88.

19. Demengel, F., Demengel, G., *Functional Spaces for the Theory of Elliptic Partial Differential Equations*. Springer London (2012).

20. Drogoul, A., Numerical Analysis of the Topological Gradient Method for Fourth Order Models and Applications to the Detection of Fine Structures in Imaging. *SIAM J. Imaging Sci.* **7**(4) (2014), 2700-2731.

21. Drogoul, A., Méthode du gradient topologique pour la détection de contours et de structures fines en imagerie. PhD thesis, Université de Nice-Sophia Antipolis, France, 2014.

22. Evans, L. C., Gariepy, R. F., *Measure Theory and Fine Properties of Functions*. CRC Press (1992).

23. Forcadel, N., Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness. *Differ. Integral Equ.* **21**(3-4) (2008), 285-304.

24. Frangi A.F., Niessen W.J., Vincken K.L., Viergever M.A., *Multiscale vessel enhancement filtering*. In: Wells W.M., Colchester A., Delp S. (eds) Medical Image Computing and Computer-Assisted Intervention MICCAI'98. MICCAI 1998. Lecture Notes in Computer Science, vol. 1496. Springer, Berlin, Heidelberg (1998).

25. Gao, Y., Bredies, K., Infimal convolution of oscillation total generalized variation for the recovery of images with structured texture. *SIAM J. Imaging Sci.* **11**(3) (2018), 2021-2063.

26. Gilboa, G., Osher, S., Nonlocal operators with applications to image processing. *Multiscale Model. Simul.* **7**(3) (2008), 1005-1028.

27. Hait, E., Gilboa, G., Spectral Total-Variation Local Scale Signatures for Image Manipulation and Fusion. *IEEE T. Image Process.* **28**(2) (2019), 880-895.

28. Jung, M., Bresson, X., Chan, T., Vese, L., Nonlocal Mumford-Shah regularizers for color image restoration. *IEEE transactions on image processing* **20**(6) (2011), 1583-1598.

29. Kim, Y., Vese, L., Image recovery using functions of bounded variation and Sobolev spaces of negative differentiability. *Inverse Problems and Imaging* **3**(1) (2009), 43-68.

30. Lellman, J., Papafitsoros, K., Schönlieb, C., Spector, D., Analysis and application of a non-local Hessian. *SIAM J. Imaging Sci.* **8**(4) (2015), 2161-2202.

31. Lieu, L., Vese, L., Image restoration and decomposition via bounded total variation and negative Hilbert-Sobolev spaces. *Applied Mathematics & Optimization* **58**(2) (2008), 167-193.

32. Meyer, Y., Oscillating Patterns in Image Processing and Nonlinear Evolution Equations. *Amer. Math. Soc.* **22** (2002).

33. Ono, S., Miyata, T., Yamada, I., Cartoon-Texture Image Decomposition Using Blockwise Low-Rank Texture Characterization. *IEEE T. Image Process.* **23**(3) (2014), 1128-1142.

34. Osher, S., Solé, A., Vese, L., Image decomposition and restoration using total variation minimization and the $H^{-1}$ norm. *Multiscale Modeling & Simulation* **1**(3) (2003), 349-370.

35. Ponce, A. C., A new approach to Sobolev spaces and connections to $\Gamma$-convergence. *Calc. Var. Partial Differential Equations* **19**(3) (2004), 229-255.
36. Rochery, M., Jermyn, I.H., Zerubia, J., Higher Order Active Contours. *Int. J. Comput. Vision* **69**(1) (2006), 27–42.

37. Rudin, L., Osher, S., Fatemi, E., Nonlinear Total Variation Based Noise Removal Algorithms. *Phys. D* **60**(1-4) (1992), 259–298.

38. Schaeffer, H., Osher, S., A Low Patch-Rank Interpretation of Texture. *SIAM J. Imaging Sci.* **6**(1) (2013), 226–262.

39. Schaeffer, H., Garnett, J., Vese, L. A., A texture model based on a concentration of measure. *Inverse Problems & Imaging* **7**(3) (2013), 927–946.

40. Spector, D., Characterization of Sobolev and BV Spaces. PhD thesis, Carnegie Mellon University (2011).

41. Stoica, R., Descombes, X., Zerubia, J., A Gibbs Point Process for Road Extraction from Remotely Sensed Images. *Int. J. Comput. Vision* **57**(2) (2004), 121–136.

42. Strang, G., *L¹ and L∞ Approximation of Vector Fields in the Plane*. In: Nonlinear Partial Differential Equations in Applied Science; Proceedings of The U.S.-Japan Seminar, 273–288. North-Holland (1982).

43. Strong, D. M. and Aujol, J.-F. and Chan, T. F., Scale Recognition, Regularization Parameter Selection, and Meyer's *G* Norm in Total Variation Regularization. *Multiscale Modeling & Simulation* **5**(1) (2006), 273–303.

44. Vese, L., Osher, S., Modeling Textures With Total Variation Minimization And Oscillating Patterns In Image Processing. *J. Sci. Comput.*, 19, 553-572 (2003)

45. Zana, F., Klein, J.-C., Segmentation of Vesell-Like Patterns Using Mathematical Morphology and Curvature Evaluation. *IEEE T. Image Process.* **10**(7) (2001), 1010–1019.

46. Zhu, H., Meng, F., Cai, J. and Lu, S., Beyond pixels: A comprehensive survey from bottom-up to semantic image segmentation and cosegmentation. *J. Vis. Commun. Image R.* **34** (2016), 12–27.

<figure><img src="image_29.png" /></figure>

Noémie Debroux received the French engineering degree (equivalent to M.Sc. degree) in mathematical engineering and the Postgraduate degree in engineering sciences from the National Institute of Applied Sciences (NIAS), Rouen, France, in 2015, and the Ph.D. degree in applied mathematics from NIAS, Rouen, in 2018. Her thesis was dedicated to joint models for segmentation and registration studied from both theoretical and applied point of views. She held an Associate Research position at the University of Cambridge from April 2018 to August 2019. She is currently an Assistant Professor at the University of Clermont Auvergne, France. Her research interests include joint variational methods in image processing, medical image analysis (liver and heart MRI), spatiotemporal modelling and segmentation under topology constraints.

<figure><img src="image_30.png" /></figure>

Carole Le Guyader received the French engineering degree (equivalent to M.Sc. degree) in mathematical engineering and the Postgraduate degree in engineering sciences from the National Institute of Applied Sciences (NIAS), Rouen, France, in 2002, and the Ph.D. degree in applied mathematics from NIAS, Rouen, in 2004. Her thesis was dedicated to segmentation under geometrical constraints from a theoretical and applied point of view. She held visiting research positions at the University of California, Los Angeles (Summer 2005 and 2007) and then a position of Assistant Professor at the NIAS, Rennes. She is currently a Full Professor at the NIAS, Rouen, France. Her research interests include variational methods in image processing, hybridization of variational and machine-learning based methods, curve and surface evolution using PDEs, viscosity solutions, and approximation of surfaces.

<figure><img src="image_31.png" /></figure>

Luminita Vese is Professor of Mathematics at the University of California, Los Angeles. She obtained the B.S. and M.S. in mathematics from the West University of Timisoara, Romania (1992), and the M.S. and Ph.D. from the University of Nice Sophia Antipolis, France (1996). During 1996 and 2002 she held postdoctoral, teaching and research positions at the University of Nice, University of Paris IX-Dauphine, and UCLA. Her research interests are in calculus of variations and partial differential equations with applications to image processing and analysis. She is co-author with Carole Le Guyader of the book *Variational Methods in Image Processing*.