# Detecting immersed obstacle in Stokes fluid flow using the coupled complex
boundary method

J. F. T. Rabago  L. Afraites  H. Notsu

###### Abstract

A non-conventional shape optimization approach is introduced to address the
identification of an obstacle immersed in a fluid described by the Stokes
equation within a larger bounded domain, relying on boundary measurements on
the accessible surface. The approach employs tools from shape optimization,
utilizing the coupled complex boundary method to transform the over-specified
problem into a complex boundary value problem by incorporating a complex Robin
boundary condition. This condition is derived by coupling the Dirichlet and
Neumann boundary conditions along the accessible boundary. The identification
of the obstacle involves optimizing a cost function constructed based on the
imaginary part of the solution across the entire domain. The subsequent
calculation of the shape gradient of this cost function, rigorously performed
via the rearrangement method, enables the iterative solution of the
optimization problem using a Sobolev gradient descent algorithm. The
feasibility of the method is illustrated through numerical experiments in both
two and three spatial dimensions, demonstrating its effectiveness in
reconstructing obstacles with pronounced concavities under high-level noise-
contaminated data, all without perimeter or volume functional penalization.

Keywords: coupled complex boundary method, obstacle problem, shape inverse
problem, Stokes equation, shape optimization, rearrangement method, and
adjoint method

###### Contents

  1. 1 Introduction
    1. 1.1 Problem setting
    2. 1.2 Shape optimization techniques for shape recovery
    3. 1.3 The CCBM and contribution to the literature
  2. 2 Preliminaries
    1. 2.1 The CCBM formulation
    2. 2.2 Well-posedness of the CCBM formulation
    3. 2.3 The shape optimization problem
  3. 3 First-order shape sensitivity analysis
    1. 3.1 Some concepts from shape calculus
    2. 3.2 The shape gradient of $J$
    3. 3.3 Computation of the shape gradient: Proof of Theorem 3.1
  4. 4 Numerical experiments
    1. 4.1 Numerical algorithm
    2. 4.2 Forward problem
    3. 4.3 Inversion procedure
    4. 4.4 Numerical tests in 2D
      1. 4.4.1 Large convex obstacles
      2. 4.4.2 Large non-convex obstacle
      3. 4.4.3 Effect of the choice of initial guess
      4. 4.4.4 Noisy data
      5. 4.4.5 Influence of the size of the obstacle
      6. 4.4.6 An obstacle with reentrant corner
      7. 4.4.7 Detecting more than one object
    5. 4.5 Numerical tests in 3D
      1. 4.5.1 Large convex obstacle
      2. 4.5.2 Large non-convex obstacle
      3. 4.5.3 Two small cubic-shaped obstacles
  5. 5 Conclusion
  6. A Lemmata Proofs
    1. A.1 Proof of Lemma 3.4
    2. A.2 Proof of Lemma 3.5
    3. A.3 Proof of Lemma 3.6

##  1 Introduction

In this study, we are interested in the problem of reconstructing the shape of
an obstacle or inclusion $\omega$ immersed in a fluid living within a larger
bounded domain $D$ from measurements taken on the accessible boundary. Our
main objective is to numerically solve this shape inverse problem through a
non-conventional shape optimization approach, which consists of a system of
complex partial differential equations (PDEs).

###  1.1 Problem setting

Let $D\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$, be a given bounded Lipschitz
domain, i.e., $D$ is a bounded, open, and connected set. We fix a small number
$\delta>0$ and denote by $\mathcal{O}_{\delta}$ the set of all subdomain
$\omega$ with $\mathcal{C}^{1,1}$ boundary $\Gamma:=\partial\omega$ such that
$d(x,\Sigma)>\delta$ for all $x\in\omega$, where $\Sigma:=\partial{D}$, and
$\Omega:=D\setminus\overline{\omega}$ is connected. We call
$\mathcal{O}_{\delta}$ the set of admissible geometries (or inclusions).
Hereinafter, we call $\Omega$ an admissible domain if
$\Omega=D\setminus\overline{\omega}$ for some $\omega\in\mathcal{O}_{\delta}$.

For a given Neumann boundary data $\bm{\mathbf{g}}\in
H^{-1/2}(\Sigma)^{d}:=(H^{1/2}(\Sigma;\mathbb{R}^{d}))^{\prime}$,
$\bm{\mathbf{g}}\not\equiv\bm{\mathbf{0}}$, and a corresponding Dirichlet
boundary measurement $\bm{\mathbf{f}}\in
H^{1/2}(\Sigma)^{d}:=H^{1/2}(\Sigma;\mathbb{R}^{d})$ (or vice versa), the
inclusion $\omega$, the vector function
$\bm{\mathbf{u}}=\bm{\mathbf{u}}(\Omega)$ representing the velocity of fluid,
and a scalar function $p=p(\Omega)$ representing the pressure, satisfy the
over-specified boundary value problem (BVP):

$-\nabla\cdot\sigma(\bm{\mathbf{u}},p)=\bm{\mathbf{0}}\ \text{in
$\Omega$},\quad\operatorname{div}{\bm{\mathbf{u}}}=0\ \text{in
$\Omega$},\quad\bm{\mathbf{u}}=\bm{\mathbf{f}}\ \ \text{and}\ \
\sigma(\bm{\mathbf{u}},p)\bm{\mathbf{n}}=\bm{\mathbf{g}}\ \text{on
$\Sigma$},\quad\bm{\mathbf{u}}=\bm{\mathbf{0}}\ \text{on $\Gamma$}.$ (1)

For technical reasons, we assume – without further notice – that the boundary
data $\bm{\mathbf{f}}\in H^{1/2}(\Sigma)^{d}$,
$\bm{\mathbf{f}}\not\equiv\bm{\mathbf{0}}$, satisfies the standard flux
compatibility condition

$\int_{\Sigma}{\bm{\mathbf{u}}\cdot\bm{\mathbf{n}}}{\,{d}\sigma}=\int_{\Sigma}{\bm{\mathbf{f}}\cdot\bm{\mathbf{n}}}{\,{d}\sigma}=\langle\bm{\mathbf{f}},\bm{\mathbf{n}}\rangle_{\Sigma}=0.$
(2)

In (1), $\sigma$ denotes the stress tensor defined by
$\sigma(\bm{\mathbf{u}},p):={\alpha}\big{(}\nabla\bm{\mathbf{u}}+\nabla\bm{\mathbf{u}}^{\top}\big{)}-p\,{id}$,
and ${\alpha}$ is the coefficient of kinematic viscosity while $id$ stands for
the $d$-dimensional identity matrix111Here, $id$ also denotes the identity
operator. If there is no confusion, this abuse of notation is used throughout
the paper.. The notation $\nabla^{\top}\bm{\mathbf{u}}$ stands for the
transposed of the matrix $\nabla{\bm{\mathbf{u}}}$. Because
$\operatorname{div}{\bm{\mathbf{u}}}=0$ in $\Omega$, then we can actually
write
$-\operatorname{div}{\sigma(\bm{\mathbf{u}},p)}=-{\alpha}\Delta\bm{\mathbf{u}}+\nabla
p$ in $\Omega$. In this case, we simply express
$\sigma(\bm{\mathbf{u}},p)\bm{\mathbf{n}}={\alpha}\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}-p\bm{\mathbf{n}}$
on $\Sigma$, where $\partial_{n}\bm{\mathbf{u}}$ denotes the outward normal
derivative (with $\bm{\mathbf{n}}$ being the outward unit normal vector to
$\Omega$). Here, we will use both of these expressions interchangeably.

Henceforth, we assume here that there exists $\omega\in\mathcal{O}_{\delta}$
such that (1) admits a solution. This means that the measurement
$\bm{\mathbf{f}}$ is perfect, that is to say without error. Therefore, we can
consider the shape inverse problem formulated as follows:

######  Problem 1.1.

Find $w\in\mathcal{O}_{\delta}$ and the pair $(\bm{\mathbf{u}},p)$ that
satisfies the over-specified system (1).

If $(\omega,\bm{\mathbf{u}},p)$ solves Problem 1.1, then we call the triplet a
solution of (1).

As mentioned earlier, we will tackle the inverse problem of reconstructing
$\omega$, employing tools from shape optimization. In this direction, we first
need to discuss the identifiability of the inclusion from boundary
measurements. This fundamental question in the Dirichlet case was addressed by
Alvarez-Conca-Friz-Kavian-Ortega in [ACF+05]:

######  Theorem 1.2 ([ACF+05]).

Let ${D}\subset\mathbb{R}^{d}$, $d\in\\{2,3\\}$, be a bounded
$\mathcal{C}^{1,1}$ domain, and $\Sigma^{0}$ be a non-empty open subset of
$\Sigma:=\partial{D}$. Moreover, let

$\omega^{1},\omega^{2}\in\\{\omega\Subset{D}\mid\ \text{$\omega$ is open,
Lipschitz, and $\Omega:=D\setminus\overline{\omega}$ is connected}\\},$

and $\bm{\mathbf{f}}\in H^{1/2}(\Sigma)^{d}$ with
$\bm{\mathbf{f}}\neq\bm{\mathbf{0}}$, satisfying the flux condition
$\displaystyle\int_{\Sigma}{\bm{\mathbf{f}}\cdot\bm{\mathbf{n}}}{\,{d}\sigma}=0$.

Let $(\bm{\mathbf{u}}^{j},p^{j})$, for $j=1,2$, be a solution of

$\left\\{\begin{array}[]{rcll}-\operatorname{div}{\sigma(\bm{\mathbf{u}}^{j},p^{j})}&=&\bm{\mathbf{0}}&\
\text{in $\Omega^{j}$},\\\ \operatorname{div}{\bm{\mathbf{u}}^{j}}&=&0&\
\text{in $\Omega^{j}$},\\\ \bm{\mathbf{u}}^{j}&=&\bm{\mathbf{f}}&\ \text{on
$\Sigma:=\partial{D}$},\\\ \bm{\mathbf{u}}^{j}&=&\bm{\mathbf{0}}&\ \text{on
$\Gamma^{j}:=\partial\omega^{j}$},\end{array}\right.$

If $(\bm{\mathbf{u}}^{1},p^{1})$ and $(\bm{\mathbf{u}}^{2},p^{2})$ are such
that

$\sigma(\bm{\mathbf{u}}^{1},p^{1})\bm{\mathbf{n}}=\sigma(\bm{\mathbf{u}}^{2},p^{2})\bm{\mathbf{n}}\quad\text{on
$\Sigma^{0}$},$

then

$\omega^{1}\equiv\omega^{2}.$

A similar result can also be proven for the case when the working problem is
given with purely the Neumann boundary condition, and with
${\alpha}\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}-p\bm{\mathbf{n}}$ is
given instead of $\sigma(\bm{\mathbf{u}},p)\bm{\mathbf{n}}$, $\Omega$ a
Lipschitz regular domain, and $\bm{\mathbf{f}}\in H^{1/2}(\Sigma)^{d}$, see
[BCD11, Thm. 2.2]. This aforementioned result states that given a fixed
$\bm{\mathbf{f}}$, two different geometries $\omega^{1}$ and $\omega^{2}$ in
$\mathcal{O}_{\delta}$ yield two different measures ${\bm{\mathbf{g}}}^{1}$
and ${\bm{\mathbf{g}}}^{2}$.

###  1.2 Shape optimization techniques for shape recovery

To recover the shape of the obstacle $\omega$, a common strategy involves
minimizing a cost function. Various types of cost functions may be considered
for this purpose; see [BCD11, Cau13, CDKT13] for examples of existing cost
functions. A conventional approach is to consider the Neumann-least-squares-
fitting cost function given by [BCD11, Cau13]

$J_{N}(\omega)=\cfrac{1}{2}\int_{\Sigma}|\sigma(\bm{\mathbf{u}}_{D}(\omega),p_{D}(\omega))\bm{\mathbf{n}}-{\bm{\mathbf{g}}}|^{2},$
(3)

where ${\bm{\mathbf{g}}}$ is the boundary measurement and
$(\bm{\mathbf{u}}_{D},p_{D})$ is the solution of the pure Dirichlet BVP:

$-{\alpha}\Delta\bm{\mathbf{u}}_{D}+\nabla p_{D}=\bm{\mathbf{0}}\ \text{ in
}\Omega,\qquad\operatorname{div}{\bm{\mathbf{u}}_{D}}=0\ \text{ in
}\Omega,\qquad\bm{\mathbf{u}}_{D}=\bm{\mathbf{0}}\ \text{ on
}\Gamma,\qquad\bm{\mathbf{u}}_{D}=\bm{\mathbf{f}}\ \text{ on }\Sigma.$ (4)

Given the compatibility condition (2), problem (4) admits a unique solution
once a normalization condition on the pressure $p$ is imposed (see, e.g.,
[Tem01, Chap. 1]). In the minimization setting ‘$J_{N}(\omega)\to\inf$’, one
can choose the normalization
$\langle\sigma(\bm{\mathbf{u}}_{D},p_{D})\bm{\mathbf{n}},\bm{\mathbf{n}}\rangle_{\Sigma}=\langle\bm{\mathbf{g}},\bm{\mathbf{n}}\rangle_{\Sigma}$.
To ensure that the cost function $J_{N}$ is well-defined, it is necessary to
assume that the Dirichlet data input is such that $\bm{\mathbf{f}}\in
H^{3/2}(\Sigma)^{d}$. To avoid this higher regularity assumption, it is
usually preferable to fit Dirichlet data rather than Neumann data. In this
case, one evaluates the misfit of the Dirichlet data using the least-squares
method:

$J_{D}(\omega)=\cfrac{1}{2}\int_{\Sigma}|\bm{\mathbf{u}}_{N}(\omega)-\bm{\mathbf{f}}|^{2},$
(5)

where $(\bm{\mathbf{u}}_{N},p_{N})$ satisfies the mixed Dirichlet-Neumann BVP:

$-{\alpha}\Delta\bm{\mathbf{u}}_{N}+\nabla p_{N}=\bm{\mathbf{0}}\ \text{ in
}\Omega,\quad\operatorname{div}{\bm{\mathbf{u}}_{N}}=0\ \text{ in
}\Omega,\quad\bm{\mathbf{u}}_{N}=\bm{\mathbf{0}}\ \text{ on
}\Gamma,\quad\sigma(\bm{\mathbf{u}}_{N},p_{N})\bm{\mathbf{n}}={\bm{\mathbf{g}}}\
\text{ on }\Sigma.$ (6)

It can be shown that system (6) admits a unique solution once a normalization
condition on $\bm{\mathbf{u}}_{N}$ is imposed.

Another cost function one may take into consideration is the called Kohn-
Vogelius functional defined as follows [CDK13]:

$J_{KV}(\omega)=\cfrac{1}{2}\int_{\Omega}{{\alpha}|{\nabla(\bm{\mathbf{u}}_{D}(\omega)-\bm{\mathbf{u}}_{N}(\omega))}|^{2}}{\,{d}x},$
(7)

where $\bm{\mathbf{u}}_{N}$ is the solution of Neumman problem (6) and
$\bm{\mathbf{u}}_{D}$ is the solution of the Dirichlet one (4). The cost
function $J_{KV}$ measures the energy-gap between the solutions of (4) and
(6). This function vanishes only if $\bm{\mathbf{u}}_{D}=\bm{\mathbf{u}}_{N}$
which is the case when the shape $\omega$ fits the exact inclusion.

The optimization problems and shape functionals mentioned above were studied
in the context of closely related shape inverse problems for the Stokes
equations, utilizing tools from shape optimization, as discussed in [CDK13,
CDKT13]. This involves employing the shape derivative of the associated cost
functional with respect to the shape of the region to be reconstructed. In the
three works mentioned above, the authors characterized the shape gradient of
the same functionals to facilitate numerical resolution. Inspired by these
works, we will do the same here, but through a proposal of a new shape
optimization formulation.

######  Remark 1.3 (Existence, uniqueness, and regularity).

For detailed discussions on the existence, uniqueness, and regularity of
solutions to the Stokes problem with Dirichlet boundary conditions, we
recommend referring to [BF06] or [Gal94]. Meanwhile, for the existence,
uniqueness, and local regularity results for the Stokes’ problem with mixed
conditions (6), one may consult Theorem A.1 and Theorem A.2 in [CDKT13, Appx.
A].

In this study, we propose the use of the Coupled Complex Boundary Method
(CCBM) in a shape optimization setting to numerically approximate a solution
to the inverse problem. We first emphasize that, in the classical approach (3)
and (5), the required data for fitting is carried out on the boundary.
However, this condition often renders the numerical solution unstable. In
contrast, in the case of CCBM, the data needed for fitting is transferred from
the boundary to the interior of the domain – similar to (7). This is because
CCBM couples the Neumann and Dirichlet data to derive a Robin condition.
Consequently, the Neumann data and the Dirichlet data are represented as the
real and imaginary parts, respectively, of the Robin boundary condition. As a
result, this method is more favorable in terms of numerical realization.
However, we note that the new boundary value problem is complex, leading to an
increase in the dimension of the corresponding discrete system. Nevertheless,
the method proves to be more robust, and effective numerical methods could be
explored based on the new framework.

###  1.3 The CCBM and contribution to the literature

Initially introduced by Cheng et al. in [CGHZ14] for an inverse source
problem, CCBM has found applications in various problmes. Gong et al. utilized
it for an inverse conductivity problem in [GCH17], and Cheng et al. applied it
to parameter identification in elliptic problems in [ZCG20]. Afraites extended
the method to shape inverse problems for the Laplace equation in [Afr22] and
geometric inverse source problems in [AMN22]. Rabago pioneered its use in free
boundary problems, specifically solving the exterior case of the Bernoulli
problem in [Rab22]. Rabago and Notsu further extended the use of CCBM to free
surface problems for the Stokes equations in [RN24]. Notably, there has been
no prior application of CCBM to obstacle problems except in [OCNN22], where it
was used to solve an inverse Cauchy Stokes problem. In this study, the authors
employed CCBM as a numerical strategy to determine fluid velocity and flux
over a specific boundary portion, reconstructing an obstacle immersed in a
fluid described by the Stokes equation through numerical shape optimization
techniques. The novelty of our work lies in applying CCBM to shape inverse
obstacle problems in fluid flows. Emphasizing the motivation, it is worth
noting that, in the context of free boundary problems [Rab22] and free surface
problems [RN24], CCBM outperforms the classical least-squares approach in
terms of stability and accuracy in obtaining the expected optimal shape
solution. Additionally, in terms of computational cost, CCBM is roughly
comparable to the Kohn-Vogelius method, see [Rab22].

We highlight the contribution of this examination to the literature. From a
theoretical point of view, we emphasize the following.

  * •

This study introduces a novel cost functional (14) to address the shape
inverse obstacle Problem 1.1. The formulation is derived from the CCBM
formulation (8). Notably, this functional has not been explored in the current
literature within the context of the present problem.

  * •

The study is centered on the computation of the first-order shape derivative
of the cost functional associated with CCBM (Theorem 3.1). This computation is
achieved only through the Hölder continuity of state variables (Lemma 3.6),
following the methodology outlined in [RN24]. Our approach deviates from the
conventional chain rule method (see, e.g., [BCD11, CDKT13]) and the
traditional Lagrangian formulation technique (see [KO22]), eliminating the
need to calculate material or shape derivatives of states.

Meanwhile, from a numerical standpoint, we depart innovatively from the
methodology applied in [CDKT13]. Our approach distinguishes itself in three
key aspects.

  * •

In the previously mentioned work, the authors addressed the ill-posedness of
the problem by employing parametric regularization, incorporating a parametric
model of shape variations. This involved imposing restrictions on domains,
confining them to star-shaped configurations, and using polar coordinates for
parametrization. This approach led to a unique method for calculating
deformation fields. In contrast, our approach relies on a straightforward
Sobolev gradient-based method inspired by [DMNV07]. Unlike the parametric
method, we avoid the complexities associated with determining the number of
parameters needed for a satisfactory approximation of the target shape. The
parametric approach, as presented in [CDKT13], faces challenges in accurately
determining this number. Insufficient parameters may hinder the detection of
non-trivial shapes, while an excess may result in functional degeneracy.
Consequently, this approach, in our view, poses implementation difficulties,
especially in three dimensions (3D).

  * •

In our numerical simulations, we use lower-degree finite elements, in contrast
to [CDKT13], which employed higher-degree elements (refer to Remark 4.3). It
is crucial to highlight that utilizing higher-degree finite elements, like
$P4/P3$ finite elements, for computing velocity-pressure variables in the
Stokes equations, especially in 3D cases, is computationally expensive.

  * •

The work in [KO22] solely computes the shape derivative of a Lagrangian
functional without including any numerical experiments. Additionally, to the
best of our knowledge, prior numerical investigations on shape inverse
obstacle problems have exclusively concentrated on two-dimensional (2D)
settings (refer to, e.g., [CDKT13]). In contrast, our current study extends
this scope by evaluating the effectiveness of our proposed methodology in 3D
cases. This further warrants the execution of the present work.

The outline of this paper is organized as follows. In Section 2, we present
the coupled complex boundary method appropriate for our inverse obstacle
problem. We also reformulate it into a shape optimization setting, introducing
the least-squares fitting for the imaginary part of the complex PDE’s
solution. Then, in Section 3, we rigorously derive the shape derivative of
this cost functional. After that, in Section 4, we provide an algorithm based
on the Sobolev gradient method with the finite element scheme. We then offer
various numerical examples in two and three spatial dimensions to illustrate
the feasibility and effectiveness of the proposed method. We conclude the
paper with a brief summary in Section 5. Finally, Appendix A is dedicated to
proving the key lemmas used to derive the shape gradient via the rearrangement
method.

##  2 Preliminaries

Our immediate goal is to rephrase the inverse problem (5) using CCBM-based
shape optimization. To start, we express (1) using the concept of CCBM and
analyze the well-posedness of the resulting PDE system. To save space, we
introduce notations aligned with those in [RN24].

Notations. Let $\partial_{i}:=\partial/\partial x_{i}$ and
$\nabla:=\begin{pmatrix}\partial_{1},\ldots,\partial_{d}\end{pmatrix}^{\top}$.
We define
$\partial_{\bm{\mathbf{n}}}{}:=\partial/\partial\bm{\mathbf{n}}=\sum_{i=1}^{d}n_{i}\partial_{i}$,
where $\bm{\mathbf{n}}=(n_{1},\ldots,n_{d})^{\top}\in\mathbb{R}^{d}$ is the
outward unit normal vector to $\Omega$. The inner product of column vectors
$\bm{\mathbf{a}},\bm{\mathbf{b}}\in\mathbb{R}^{d}$ is denoted by
$\bm{\mathbf{a}}\cdot\bm{\mathbf{b}}:=\bm{\mathbf{a}}^{\top}\bm{\mathbf{b}}\equiv\langle\bm{\mathbf{a}},\bm{\mathbf{b}}\rangle_{\mathbb{R}^{d}}=\langle\bm{\mathbf{a}},\bm{\mathbf{b}}\rangle$,
the latter being used when there is no confusion. For simplicity, we
occasionally use $\varphi_{n}$ instead of
$\bm{\mathbf{\varphi}}\cdot\bm{\mathbf{n}}$, where $\bm{\mathbf{\varphi}}$ is
a vector-valued function. For a vector-valued function
$\bm{\mathbf{u}}:=(u_{1},\ldots,u_{d})^{\top}:\Omega\to\mathbb{R}^{d}$,
$(\nabla\bm{\mathbf{u}})_{ij}:=(\partial_{i}u_{j})_{i,j=1,\ldots,d}$
represents its gradient, while its Jacobian is denoted by
$(D\bm{\mathbf{u}})_{ij}=(\partial_{j}u_{i})_{i,j=1,\ldots,d}\equiv\nabla^{\top}\bm{\mathbf{u}}$.
Accordingly, we use
$\partial_{\bm{\mathbf{n}}}\bm{\mathbf{u}}:=(D\bm{\mathbf{u}})\bm{\mathbf{n}}$.

Let $1\leqslant p\leqslant\infty$ and $m\in\mathbb{N}\cup\\{0\\}$. The
function spaces $W^{m,p}(\Omega)$, $W^{0,p}(\Omega)=L^{p}(\Omega)$, and
$H^{m}(\Omega)=W^{m,2}(\Omega)$ denote the standard real Sobolev spaces
equipped with their natural norms, also expressed in usual notations [DL98,
Chap. IV] (e.g., $\|\cdot\|_{H^{1}(\Omega)}$ is the standard
$H^{1}(\Omega)$-norm). We define the space
$H^{m}(\Omega)^{d}:=\\{\bm{\mathbf{u}}:\Omega\to\mathbb{R}^{d}\mid u_{i}\in
H^{m}(\Omega),\text{for $i=1,\ldots,d$}\\}$ with norm
$\|\bm{\mathbf{u}}\|_{H^{m}(\Omega)^{d}}:=(\sum_{i=1}^{d}|u_{i}|^{2}_{H^{m}(\Omega)})^{1/2}$.
A similar definition is given when $\Omega$ is replaced by $\partial\Omega$,
with $m$ as a rational number.

In this paper, we let $\bm{\mathbf{\mathsf{H}}}^{m}(\Omega)^{d}$ be the
complex version of $H^{m}(\Omega)^{d}$ with the inner product
$(\\!(\cdot,\cdot)\\!)_{m,\Omega,d}$ and norm
${|\kern-1.07639pt|\kern-1.07639pt|\,\cdot\,|\kern-1.07639pt|\kern-1.07639pt|}_{m,\Omega,d}$
defined respectively as follows: for all
$\bm{\mathbf{u}},\bm{\mathbf{v}}\in\bm{\mathbf{\mathsf{H}}}^{m}(\Omega)^{d}$,
$(\\!(\bm{\mathbf{u}},\bm{\mathbf{v}})\\!)_{m,\Omega,d}=\sum_{j=1}^{d}(\bm{\mathbf{u}},\overline{\bm{\mathbf{v}}})_{m,\Omega}$
and
${|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{v}}|\kern-1.07639pt|\kern-1.07639pt|}_{m,\Omega,d}=\sqrt{(\\!(\bm{\mathbf{v}},\bm{\mathbf{v}})\\!)_{m,\Omega,d}}$.
Note that, for real valued functions,
${|\kern-1.07639pt|\kern-1.07639pt|\cdot|\kern-1.07639pt|\kern-1.07639pt|}_{m,\Omega,d}=\|\cdot\|_{W^{m,p}(\Omega)}$.
Here, the operation ‘$:$’ stands for the Frobenius inner product and is
defined as
$\nabla\bm{\mathbf{v}}:\nabla\overline{\bm{\mathbf{w}}}=\sum_{i,j=1}^{d}{\partial_{i}{v}_{j}}{\partial_{i}{\overline{w}}_{j}}\in\mathbb{R}$.
With these definitions, we sometimes write
${|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{v}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}={|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{v}}|\kern-1.07639pt|\kern-1.07639pt|}_{1,\Omega,d}$
and
${|\kern-1.07639pt|\kern-1.07639pt|q|\kern-1.07639pt|\kern-1.07639pt|}_{Q}={|\kern-1.07639pt|\kern-1.07639pt|q|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega,d}$,
dropping $d$ for convenience. For ease of writing, we use the function space
notations $X:=\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)^{d}$,
${V_{\Gamma}}:=\bm{\mathbf{\mathsf{H}}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega)^{d}=\\{\bm{\mathbf{\varphi}}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)^{d}\mid\bm{\mathbf{\varphi}}=\bm{\mathbf{0}}\
\text{on $\Gamma$}\\}$, $Q:=\bm{\mathbf{\mathsf{L}}}^{2}(\Omega)$, and
${Q}_{0}=\left\\{q\in Q\mid\int_{\Omega}{q}{\,{d}x}=0\right\\}$.

We use the symbol $c$ to represent a positive constant, and its specific value
may vary in different contexts. The symbol ‘$\lesssim$’ denotes that if
$x\lesssim y$, then there exists a constant $c>0$ such that $x\leqslant cy$.
Similarly, $y\gtrsim x$ is defined as $x\lesssim y$. Lastly, we adhere to
standard notations, emphasizing them only for clarity.

###  2.1 The CCBM formulation

The CCBM formulation of (1) is given by the complex PDE system

$-{\alpha}\Delta{\bm{\mathbf{u}}}+\nabla p={\bm{\mathbf{0}}}\ \text{ in
}\Omega,\quad\operatorname{div}{\bm{\mathbf{u}}}={\bm{\mathbf{0}}}\ \text{ in
}\Omega,\quad\sigma(\bm{\mathbf{u}},p)\bm{\mathbf{n}}+i\bm{\mathbf{u}}={{\bm{\mathbf{g}}}+i\bm{\mathbf{f}}}\
\text{ on }\Sigma,\quad\bm{\mathbf{u}}={\bm{\mathbf{0}}}\ \text{ on }\Gamma,$
(8)

where $i=\sqrt{-1}$ is the imaginary unit. We call the solution of (8) as
states or state solutions.

Let
$(\bm{\mathbf{u}},p)=({\bm{\mathbf{u}}_{r}}+i{\bm{\mathbf{u}}_{i}},p_{r}+i{p_{i}})$222This
notation is used hereon without further notice. be a solution pair to (8).
Then, the real-valued functions ${\bm{\mathbf{u}}_{r}}$,
${\bm{\mathbf{u}}_{i}}$, $p_{r}$, and $p_{i}$ satisfy the following systems of
real PDEs:

$-{\alpha}\Delta{\bm{\mathbf{u}}_{r}}+\nabla p_{r}=\bm{\mathbf{0}}\ \text{ in
}\Omega,\quad\operatorname{div}{{\bm{\mathbf{u}}_{r}}}=0\ \text{ in
}\Omega,\quad\sigma({\bm{\mathbf{u}}_{r}},p_{r})\bm{\mathbf{n}}={\bm{\mathbf{g}}}+{\bm{\mathbf{u}}_{i}}\
\text{ on }\Sigma,\quad{\bm{\mathbf{u}}_{r}}=\bm{\mathbf{0}}\ \text{ on
}\Gamma;$ (9)
$-{\alpha}\Delta{\bm{\mathbf{u}}_{i}}+\nabla{p_{i}}=\bm{\mathbf{0}}\ \text{ in
}\Omega,\quad\operatorname{div}{{\bm{\mathbf{u}}_{i}}}=0\ \text{ in
}\Omega,\quad\sigma({\bm{\mathbf{u}}_{i}},{p_{i}})\bm{\mathbf{n}}=\bm{\mathbf{f}}-{\bm{\mathbf{u}}_{r}}\
\text{ on }\Sigma,\quad{\bm{\mathbf{u}}_{i}}=\bm{\mathbf{0}}\ \text{ on
}\Gamma.$ (10)

######  Remark 2.1.

If ${\bm{\mathbf{u}}_{i}}=\bm{\mathbf{0}}$ and ${p_{i}}=0$ in $\Omega$, then
${\bm{\mathbf{u}}_{i}}=\partial_{n}{\bm{\mathbf{u}}_{i}}=\bm{\mathbf{0}}$ on
$\Sigma$. As a result, from the PDE systems (9) and (10),
$(\omega,{\bm{\mathbf{u}}_{r}},p_{r})$ is a solution to the original problem
(1). Conversely, if $(\omega,\bm{\mathbf{u}},p)$ is a solution of (1), then
clearly it satisfies (8).

Remark 2.1 implies that Problem 1.1 can equivalently be formulated as follows:

######  Problem 2.2.

Find $\omega\in\mathcal{O}_{\delta}$ such that
${\bm{\mathbf{u}}_{i}}=\bm{\mathbf{0}}$ and ${p_{i}}=0$ in $\Omega$, where
$(\bm{\mathbf{u}},p)$ satisfy (8).

###  2.2 Well-posedness of the CCBM formulation

Here, we briefly examine the well-posedness of the complex PDE system (8). To
do so, we introduce the following forms:

$\left\\{\begin{aligned}
a({\bm{\mathbf{\varphi}}},{\bm{\mathbf{\psi}}})&=\int_{\Omega}{{\alpha}\nabla{\bm{\mathbf{\varphi}}}:\nabla{\overline{\bm{\mathbf{\psi}}}}}{\,{d}x}+i\int_{\Sigma}{{\bm{\mathbf{\varphi}}}\cdot{\overline{\bm{\mathbf{\psi}}}}}{\,{d}\sigma},\quad\text{where}\
{\bm{\mathbf{\varphi}}},{\bm{\mathbf{\psi}}}\in{V_{\Gamma}},\\\
b(\bm{\mathbf{\varphi}},\lambda)&=-\int_{\Omega}{\lambda(\nabla\cdot\overline{\bm{\mathbf{\varphi}}})}{\,{d}x},\quad\text{where}\
{\bm{\mathbf{\varphi}}}\in{V_{\Gamma}},\ \lambda\in{Q},\\\
F(\bm{\mathbf{\psi}})&=\int_{\Sigma}{(\bm{\mathbf{g}}+i\bm{\mathbf{f}})\cdot\overline{\bm{\mathbf{\psi}}}}{\,{d}\sigma},\quad\text{where}\
{\bm{\mathbf{\psi}}}\in{V_{\Gamma}}.\end{aligned}\right.$ (11)

With the above forms, we can state the weak formulation of (8) as follows:

$\text{Find\quad$({\bm{\mathbf{u}}},p)\in V_{\Gamma}\times Q$\quad such that
\quad}\left\\{\begin{array}[]{rcll}a({\bm{\mathbf{u}}},{\bm{\mathbf{\varphi}}})+b(\bm{\mathbf{\varphi}},p)&=&F(\bm{\mathbf{\varphi}}),&\
\forall{\bm{\mathbf{\varphi}}}\in{V_{\Gamma}},\\\
b(\bm{\mathbf{u}},\lambda)&=&0,&\ \forall\lambda\in{Q}.\\\ \end{array}\right.$
(12)

The well-posedness of the above variational problem is issued in the following
lemma.

######  Lemma 2.3.

Let $(\bm{\mathbf{f}},{\bm{\mathbf{g}}})\in H^{1/2}(\Sigma)^{d}\times
H^{-1/2}(\Sigma)^{d}$ be given. Then, problem (12) admits a (unique) solution
$(\bm{\mathbf{u}},p)\in{V_{\Gamma}}\times Q$. This solution depends
continuously on the data and we have

${|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}|\kern-1.07639pt|\kern-1.07639pt|}_{X},\
{|\kern-1.07639pt|\kern-1.07639pt|p|\kern-1.07639pt|\kern-1.07639pt|}_{Q}\lesssim{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{f}}|\kern-1.07639pt|\kern-1.07639pt|}_{H^{1/2}(\Sigma)^{d}}+{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{g}}|\kern-1.07639pt|\kern-1.07639pt|}_{H^{-1/2}(\Sigma)^{d}}=:\left\|{(\bm{\mathbf{f}},\bm{\mathbf{g}})}\right\|_{1/2,\Sigma}.$

The above claim can be validated through standard steps (cf. [RN24, Appx. A]).
In fact, it follows from the continuity of the sesquilinear form
$a(\cdot,\cdot)$ on $V_{\Gamma}\times V_{\Gamma}$, its coercivity on
$V_{\Gamma}\subset X$ (i.e., there is a constant $c_{a}>0$ such that the
relaxed $V_{\Gamma}$-ellipticity assumption
$\Re(a(\bm{\mathbf{\varphi}},\bm{\mathbf{\varphi}}))\geqslant
c_{a}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}$,
for all $\bm{\mathbf{\varphi}}\in{V_{\Gamma}}$), the boundedness of $F$, and
on an inf-sup condition; that is, there exists a constant $\beta_{0}>0$ such
that

$\inf_{\begin{subarray}{c}{\lambda}\in{Q}\\\ {\lambda}\neq
0\end{subarray}}\sup_{\begin{subarray}{c}\bm{\mathbf{\varphi}}\in{V_{\Gamma}}\\\
\bm{\mathbf{\varphi}}\neq\bm{\mathbf{0}}\end{subarray}}\frac{b(\bm{\mathbf{\varphi}},{\lambda})}{{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}{|\kern-1.07639pt|\kern-1.07639pt|{\lambda}|\kern-1.07639pt|\kern-1.07639pt|}_{Q}}\geqslant\beta_{0}.$
(13)

###  2.3 The shape optimization problem

We now present the solution method we use to resolve the shape inverse problem
2.2. To this end, we introduce the cost functional

$J(\Omega):=J(D\setminus\overline{\omega})=\frac{1}{2}\Big{\\{}\big{\|}{\bm{\mathbf{u}}_{i}}\big{\|}_{L^{2}(\Omega)^{d}}^{2}+\big{\|}{p_{i}}\big{\|}_{L^{2}(\Omega)}^{2}\Big{\\}}$
(14)

and consider the following shape optimization problem:

######  Problem 2.4.

Find $\omega^{\ast}\in\mathcal{O}_{\delta}$ such that
$\omega^{\ast}=\displaystyle\inf_{\omega\in\mathcal{O}_{\delta}}J(\Omega)$.

######  Remark 2.5.

CCBM enables us to define the cost function $J$ across the entire domain
$\Omega$, offering advantages in robustness for reconstruction, such as the
Kohn-Vogelius method, when compared to the cost functionals of boundary-
integral-least-squares fitting, namely $J_{N}$ and $J_{D}$.

To find the shape derivative of the cost function (14) with respect to
$\omega$ (or equivalently, to $\Omega$), we can use the classical chain rule
approach (see, e.g., [DZ11, HP18]). This method relies on having the
expression for the shape derivative of the state problem (8). However,
obtaining the shape derivative of the states requires higher regularity
assumptions on the unknown variables, typically derived from the regularities
imposed on the domain and on the boundary data. Our goal here is to avoid such
unnecessary requirements for the first-order shape sensitivity of $J$.
Consequently, we establish the shape gradient of $J$ – following the Hadamard-
Zolésio structure theorem [DZ11, p. 479] – using only a mild
$\mathcal{C}^{1,1}$ regularity assumption on the domain, combined with the
assumption that $(\bm{\mathbf{f}},\bm{\mathbf{g}})\in
H^{3/2}(\Sigma)^{d}\times H^{1/2}(\Sigma)^{d}$. This task is achieved through
the application of an alternative technique in computing shape derivatives of
shape functionals known as the rearrangement method [IKP06, IKP08]. This
approach bypasses the computation of the material and shape derivative of the
states and utilizes only the Hölder continuity of the state variables to
obtain the expression for the shape derivative of the cost. It is worth noting
that, despite its application in computing shape derivatives for various
problems, such as shape optimization reformulations of free boundary problems
[IKP06, HIK+09, BP13, RN24], shape optimal problems [IKP08, KK11, KK14], and
shape design problems [SN22b, SN22a], there is, to the best of our knowledge,
no documented demonstration of the method applied to shape inverse problems.
This observation serves as a compelling motivation for applying this method to
our present problem.

##  3 First-order shape sensitivity analysis

To iteratively solve Problem 2.4 through a numerical procedure, our goal is to
characterize the shape gradient of the cost function $J$, as enunciated in
Theorem 3.1. To facilitate this, we will make necessary preparations by
recalling some relevant concepts from shape calculus.

###  3.1 Some concepts from shape calculus

Let us define $T_{t}$ as the perturbation of the identity $id$ given by the
map

$T_{t}=T_{t}({\bm{\mathbf{\theta}}})=id+t\bm{\mathbf{\theta}},\qquad(T_{t}:D\longmapsto
D),$ (15)

where $\bm{\mathbf{\theta}}$ is a $t$-independent deformation field belonging
to the admissible space

$\bm{\mathbf{\Theta}}:=\\{\bm{\mathbf{\theta}}=(\theta_{1},\ldots,\theta_{d})^{\top}\in\mathcal{C}^{1,1}(\overline{D})^{d}\mid\operatorname{supp}{\bm{\mathbf{\theta}}}\subset{\overline{D}}_{\delta}\\},$
(16)

where $\\{x\in{D}\mid
d(x,\partial{D})>\delta/2\\}\subset{D}_{\delta}\subset\\{x\in{D}\mid
d(x,\partial{D})>\delta/3\\}$. Clearly, by definition, we have that
$\Sigma_{t}:=T_{t}(\Sigma)\equiv\Sigma$ and $\Gamma_{t}:=T_{t}(\Gamma)$. In
addition, $\Omega_{0}=\Omega$ and $\Sigma_{0}=\Sigma$.

Hereinafter, we assume that $t>0$ is sufficiently small such that $T_{t}$ is a
diffeomorphism from $\Omega\in\mathcal{C}^{1,1}$ onto its image. Specifically,
we let $\varepsilon\in\mathbb{R}^{+}$ be small enough so that $[t\mapsto
T_{t}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}^{1,1}(\overline{D})^{d})$,
where $\mathcal{I}:=[0,\varepsilon]$. Accordingly, we define the set of all
admissible perturbations of $\Omega$ denoted by $\mathcal{O}_{ad}$ as follows:

$\mathcal{O}_{ad}=\left\\{T_{t}({\bm{\mathbf{\theta}}})(\overline{\Omega})\subset
D\mid\Omega=D\setminus{\overline{\omega}}\in\mathcal{C}^{1,1},\omega\in\mathcal{O}_{\delta},t\in\mathcal{I},\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}\right\\}.$
(17)

Note that we do not necessarily require the fixed boundary
$\Sigma=\partial\Omega$ to be $\mathcal{C}^{1,1}$ regular (it is sufficient
for it to be Lipschitz), but for the sake of simplifying the discussion, we
assume the stated regularity.

The functional $J:\mathcal{O}_{ad}\to\mathbb{R}$ has a directional first-order
Eulerian derivative at $\Omega$ in the direction of the field
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}$ if the limit $\lim_{t\searrow
0}\frac{J(\Omega_{t})-J(\Omega)}{t}=:{d}J(\Omega)[\bm{\mathbf{\theta}}]$
exists (see, e.g., [DZ11, Sec. 4.3.2, Eq. (3.6), p. 172]). The function $J$ is
considered shape differentiable at $\Omega$ in the direction of
$\bm{\mathbf{\theta}}$ if the mapping
$\bm{\mathbf{\theta}}\mapsto{d}J(\Omega)[\bm{\mathbf{\theta}}]$ is both linear
and continuous. In such cases, we call it the shape gradient of $J$.

We introduce a few more notations for ease of writing. The Jacobian matrix of
$T_{t}$, its inverse, and its inverse transpose are denoted respectively by
$DT_{t}$, $(DT_{t})^{-1}$, and $(DT_{t})^{-\top}:=((DT_{t})^{\top})^{-1}$. For
convenience, we define

$I_{t}:=\det\,DT_{t},\quad
A_{t}:=I_{t}(DT_{t}^{-1})(DT_{t})^{-\top},\quad\text{and}\quad
B_{t}:=I_{t}|{M_{t}\bm{\mathbf{n}}}|,\quad M_{t}:=(DT_{t})^{-\top},$

and we assume that, for all $t\in\mathcal{I}$, we have

$0\leqslant\Lambda_{1}\leqslant
I_{t}\leqslant\Lambda_{2}\qquad\text{and}\qquad
0<\Lambda_{3}|{\xi}|^{2}\leqslant
A_{t}\xi\cdot\xi\leqslant\Lambda_{4}|{\xi}|^{2},$ (18)

for all $\xi\in\mathbb{R}^{d}$, for some constants $\Lambda_{1}$,
$\Lambda_{2}$, $\Lambda_{3}$, and $\Lambda_{4}$ ($\Lambda_{1}<\Lambda_{2}$,
$\Lambda_{3}<\Lambda_{4}$).

At $t=0$, $I_{0}=1$, $DT_{0}^{-1}=id$, $(DT_{0})^{-\top}=id$, $A_{0}=id$, and
$B_{0}=1$. The maps $t\mapsto(DT_{t})^{\pm 1}$, $t\mapsto I_{t}$, and
$t\mapsto A_{t}$ are differentiable (see (A1)), and their respective
derivative are given as follows:

$\displaystyle\frac{d}{dt}(DT_{t})^{\pm 1}\big{\rvert}_{t=0}$
$\displaystyle=\lim_{t\to 0}\dfrac{(DT_{t})^{\pm 1}-id}{t}=\pm
D\bm{\mathbf{\theta}},\qquad\qquad\frac{d}{dt}I_{t}\big{\rvert}_{t=0}=\lim_{t\to
0}\frac{I_{t}-1}{t}=\operatorname{div}\bm{\mathbf{\theta}},$ (19)
$\displaystyle\quad\frac{d}{dt}A_{t}\big{\rvert}_{t=0}$
$\displaystyle=\lim_{t\to
0}\frac{A_{t}-id}{t}=(\operatorname{div}\bm{\mathbf{\theta}}){id}-{D}\bm{\mathbf{\theta}}-({D}\bm{\mathbf{\theta}})^{\top}=:A.$

We also introduce the following notations for economy of space:

$\mathfrak{i}_{t}:=I_{t}-1,\qquad\mathfrak{m}_{t}:=M_{t}^{\top}-id,\qquad\mathfrak{a}_{t}:=A_{t}-id,\qquad\text{and}\quad\theta_{n}:=\bm{\mathbf{\theta}}\cdot\bm{\mathbf{n}}.$
(20)

Since we intend to map a function $\bm{\mathbf{\varphi}}_{t}:\Omega_{t}\to
D\subset\mathbb{R}^{d}$ to a reference domain using the transformation
$T_{t}$, we will utilize the following notation:

$\bm{\mathbf{\varphi}}^{t}:=\bm{\mathbf{\varphi}}_{t}\circ T_{t}:\Omega\to
D\subset\mathbb{R}^{d}.$

To conclude this review and preparation section, we present some results
regarding the limits of functions composed with the map $T_{t}$.

For a vector-valued function
$\bm{\mathbf{\varphi}}=(\varphi_{1},\ldots,\varphi_{d})^{\top}\in
H^{2}(D)^{d}$, we have

$\lim_{t\to 0}\|\bm{\mathbf{\varphi}}\circ
T_{t}-\bm{\mathbf{\varphi}}\|_{H^{1}(D)^{d}}=\lim_{t\to
0}\left(\sum_{i=1}^{d}\|\varphi_{i}\circ
T_{t}-\varphi_{i}\|_{H^{1}(D)}\right)^{1/2}=0.$ (21)

Moreover, the maps $t\mapsto\bm{\mathbf{\varphi}}\circ T_{t}$ from
$\mathcal{I}\to H^{1}(\Omega)^{d}$ and $t\to I_{t}\bm{\mathbf{\varphi}}\circ
T_{t}$ from $\mathcal{I}$ to $L^{2}(\Omega)^{d}$ are differentiable at $t=0$,
and we have

$\lim_{t\to 0}\frac{1}{t}\left(\bm{\mathbf{\varphi}}\circ
T_{t}-\bm{\mathbf{\varphi}}\right)=D\bm{\mathbf{\varphi}}\bm{\mathbf{\theta}}\quad\text{and}\quad\lim_{t\to
0}\frac{1}{t}\left(I_{t}\bm{\mathbf{\varphi}}\circ
T_{t}-\bm{\mathbf{\varphi}}\right)=\bm{\mathbf{\nabla}}\cdot(\bm{\mathbf{\varphi}}\otimes\bm{\mathbf{\theta}}),$
(22)

where
$\bm{\mathbf{\nabla}}\cdot(\bm{\mathbf{\varphi}}\otimes\bm{\mathbf{\theta}})=(\operatorname{div}(\varphi_{1}\bm{\mathbf{\theta}}),\operatorname{div}(\varphi_{2}\bm{\mathbf{\theta}}),\ldots,\operatorname{div}(\varphi_{d}\bm{\mathbf{\theta}}))^{\top}$,
and $\otimes$ denotes the outer product; i.e.,
$\bm{\mathbf{\varphi}}\otimes\bm{\mathbf{\theta}}=\bm{\mathbf{\varphi}}\bm{\mathbf{\theta}}^{\top}=(\varphi_{j}\theta_{k})_{jk}$
where $j,k=1,\ldots,d$.

The proofs of the above results are omitted since they are standard and follow
similar arguments used for the case of scalar functions, see, e.g., [DZ11,
IKP06].

###  3.2 The shape gradient of $J$

Our main result in this section is the expression for the shape gradient of
$J$, which we present in the following theorem.

######  Theorem 3.1.

Let $(\bm{\mathbf{f}},\bm{\mathbf{g}})\in H^{3/2}(\Sigma)^{d}\times
H^{1/2}(\Sigma)^{d}$ be a given Cauchy pair, $\Omega$ is an admissible domain,
and $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}$. Then, $J$ is shape
differentiable with respect to $\Omega$ in the direction of the deformation
field $\bm{\mathbf{\theta}}$. Its shape derivative is given by

$\displaystyle{d}J(\Omega)[\bm{\mathbf{\theta}}]$
$\displaystyle=\left\langle{{G}\bm{\mathbf{n}}},{\bm{\mathbf{\theta}}}\right\rangle_{\Gamma}=\int_{\Gamma}{\left(\Im\left\\{\overline{\sigma(\bm{\mathbf{v}},q)\bm{\mathbf{n}}}\cdot\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\right\\}+\frac{1}{2}|p_{i}|^{2}\right)\bm{\mathbf{n}}\cdot\bm{\mathbf{\theta}}}{\,{d}\sigma},$
(23)

where $\Im$ denotes the imaginary part, $(\bm{\mathbf{u}},p)$ is the unique
pair of solution to (8), while the pair of adjoints
$(\bm{\mathbf{v}},q)=(\bm{\mathbf{v}}_{r}+i\bm{\mathbf{v}}_{i},q_{r}+iq_{i})$
uniquely solves the adjoint system

$-{\alpha}\Delta\bm{\mathbf{v}}+\nabla q=\bm{\mathbf{u}}_{i}\ \text{in
$\Omega$},\quad-\nabla\cdot\bm{\mathbf{v}}=p_{i}\ \text{in
$\Omega$},\quad\bm{\mathbf{v}}=\bm{\mathbf{0}}\ \text{on
$\Gamma$},\quad\sigma(\bm{\mathbf{v}},q)\bm{\mathbf{n}}-i\bm{\mathbf{v}}=\bm{\mathbf{0}}\
\text{on $\Sigma$},$ (24)

with compatibility condition
$\langle\bm{\mathbf{v}},\bm{\mathbf{n}}\rangle_{\Sigma}=-(p_{i},1)_{\Omega}$.

We emphasize here that if, instead, $\Omega$ is of class $\mathcal{C}^{2,1}$,
$(\mathbf{f},\bm{\mathbf{g}})\in H^{5/2}(\Sigma)^{d}\times
H^{3/2}(\Sigma)^{d}$, and $\bm{\mathbf{\theta}}$ is a $\mathcal{C}^{2,1}$
smooth deformation field, then the existence of the material and shape
derivatives of the states is guaranteed. Because then we will have
$\bm{\mathbf{u}}_{r}$, $\bm{\mathbf{u}}_{i}\in H^{3}(\Omega)^{d}$ and $p_{r}$,
$p_{i}\in H^{2}(\Omega)$. With this regularity, the shape gradient of the cost
can be computed with ease using Hadamard’s domain differentiation formula
(see, e.g., [DZ11, Thm. 4.2, p. 483]), [HP18, eq. (5.12), Thm. 5.2.2, p. 194]
or [SZ92, eq. (2.168), p. 113]):

$\displaystyle\left.\left\\{\frac{{d}}{{d}t}\int_{\Omega_{t}}{f(t,x)}{\,{d}x_{t}}\right\\}\right|_{t=0}=\int_{\Omega}{\frac{\partial}{\partial
t}f(0,x)}{\,{d}x}+\int_{\partial\Omega}{f(0,\sigma)\theta_{n}}{\,{d}\sigma}$
(25)

Before we provide a rigorous proof of Theorem 3.1 , let us exhibit shortly the
shape gradient of $J$ via the chain rule approach. Let us assume, for now,
that $\Omega$ is an admissible domain of class $\mathcal{C}^{2,1}$, and
$\bm{\mathbf{\theta}}\in\mathcal{C}^{2,1}(\overline{D})^{d}$ such that
$\operatorname{supp}{\bm{\mathbf{\theta}}}\subset{\overline{D}}_{\delta}$, and
that $(\bm{\mathbf{f}},\bm{\mathbf{g}})\in H^{5/2}(\Sigma)^{d}\times
H^{3/2}(\Sigma)^{d}$. Then, we can apply formula (25) to obtain – noting that
$\bm{\mathbf{\theta}}=\bm{\mathbf{0}}$ on $\Sigma$ and
$\bm{\mathbf{u}}=\bm{\mathbf{0}}$ on $\Gamma$ – the derivative

$\displaystyle{{d}}J(\Omega)[\bm{\mathbf{\theta}}]=\int_{\Omega}{\left(\bm{\mathbf{u}}_{i}\cdot\bm{\mathbf{u}}_{i}^{\prime}+p_{i}p_{i}^{\prime}\right)}{\,{d}x}+\frac{1}{2}\int_{\Gamma}{|p_{i}|^{2}\theta_{n}}{\,{d}\sigma}=:L_{1}+L_{2}.$
(26)

Here, the shape derivative $(\bm{\mathbf{u}}^{\prime},p^{\prime})$ of the
states $(\bm{\mathbf{u}},p)$ is characterized by the PDE system

$-{\alpha}\Delta\bm{\mathbf{u}}^{\prime}+\nabla p^{\prime}=\bm{\mathbf{0}}\
\text{in $\Omega$},\quad\nabla\cdot\bm{\mathbf{u}}^{\prime}=\bm{\mathbf{0}}\
\text{in
$\Omega$},\quad\sigma(\bm{\mathbf{u}}^{\prime},p^{\prime})+i\bm{\mathbf{u}}^{\prime}=\bm{\mathbf{0}}\
\text{on
$\Sigma$},\quad\bm{\mathbf{u}}^{\prime}=-\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\theta_{n}\
\text{on $\Gamma$}.$ (27)

We can confirm this claim using standard techniques from shape calculus [DZ11,
HP18], with arguments similar to the proof in [RN24, Thm. C.1], so we omit it.

We focus on eliminating the derivatives in the first integral $L_{1}$. To do
this, we make use of the adjoint method, utilizing the adjoint system (24).

Let us multiply by $\overline{\bm{\mathbf{u}}}^{\prime}$ equation (24), apply
integration by parts (IBP), and then use the boundary conditions in (24) to
obtain

$\int_{\Omega}{{\alpha}\nabla{\bm{\mathbf{v}}}:\nabla{\overline{\bm{\mathbf{u}}}^{\prime}}}{\,{d}x}-i\int_{\Sigma}{\bm{\mathbf{v}}\cdot\overline{\bm{\mathbf{u}}}^{\prime}}{\,{d}\sigma}-\int_{\Gamma}{({\alpha}\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{v}}}-q\bm{\mathbf{n}})\cdot\overline{\bm{\mathbf{u}}}^{\prime}}{\,{d}\sigma}=\int_{\Omega}{\bm{\mathbf{u}}_{i}\cdot\overline{\bm{\mathbf{u}}}^{\prime}}{\,{d}x}.$

Taking the complex conjugate on both sides of the above equation, we get
(using the sesquilinear form in (11))
$a(\bm{\mathbf{u}}^{\prime},\bm{\mathbf{v}})-\int_{\Gamma}{\overline{({\alpha}\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{v}}}-q\bm{\mathbf{n}})}\cdot\bm{\mathbf{u}}^{\prime}}{\,{d}\sigma}=\int_{\Omega}{\bm{\mathbf{u}}_{i}\cdot\bm{\mathbf{u}}^{\prime}}{\,{d}x}.$
Meanwhile, multiplying equation (27) by $\overline{\bm{\mathbf{v}}}$ and then
employing IBP, while taking into account the boundary conditions in (27), we
get
$a(\bm{\mathbf{u}}^{\prime},\bm{\mathbf{v}})=-\int_{\Omega}{p_{i}p^{\prime}}{\,{d}x}.$
Subtracting this equation from the previous one, we see that

$\displaystyle
L_{1}=\Im\left\\{\int_{\Omega}{\left(\bm{\mathbf{u}}\cdot\bm{\mathbf{u}}_{i}^{\prime}+pp_{i}^{\prime}\right)}{\,{d}x}\right\\}$
$\displaystyle=\Im\left\\{\int_{\Gamma}{\overline{\sigma(\bm{\mathbf{v}},q)\bm{\mathbf{n}}}\cdot\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\theta_{n}}{\,{d}\sigma}\right\\},$

from which, after combining with $L_{2}$, we recover the expression in (23).

Our next objective is to demonstrate that applying the rearrangement method
yields the same expression, with the only requirement being mild
$\mathcal{C}^{1,1}$ regularity assumptions on $\Omega$, with perturbations in
$\mathcal{O}_{ad}$, and regularity conditions $(\mathbf{f},\bm{\mathbf{g}})\in
H^{3/2}(\Sigma)^{d}\times H^{1/2}(\Sigma)^{d}$.

######  Remark 3.2.

The conclusion in Theorem 3.1 holds even when $(\mathbf{f},\bm{\mathbf{g}})\in
H^{1/2}(\Sigma)^{d}\times H^{-1/2}(\Sigma)^{d}$. In this case, the expression
of $\left\langle{{G}},{\theta_{n}}\right\rangle_{\Gamma}$ has to be seen as a
duality product $H^{1/2}\times H^{-1/2}$.

###  3.3 Computation of the shape gradient: Proof of Theorem 3.1

We are now in the position to compute the shape derivative of $J$ rigorously
via rearrangement method given that $(\bm{\mathbf{f}},\bm{\mathbf{g}})\in
H^{3/2}(\Sigma)^{d}\times H^{1/2}(\Sigma)^{d}$,
$\Omega=D\setminus\overline{\omega}$ of class $\mathcal{C}^{1,1}$ such that
$\omega\in\mathcal{O}_{\delta}$, and
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}$. The discussion then lays out
the proof of Theorem 3.1, closely following [RN24].

We start with the weak formulation of (24) which is given as follows:

$\text{find \quad$({\bm{\mathbf{v}}},q)\in V_{\Gamma}\times Q$\quad such that
\quad}\left\\{\begin{array}[]{rcll}\tilde{a}({\bm{\mathbf{v}}},{\bm{\mathbf{\psi}}})+b(\bm{\mathbf{\psi}},q)&=&\tilde{F}(\bm{\mathbf{\psi}}),&\
\forall{\bm{\mathbf{\psi}}}\in{V_{\Gamma}},\\\
b(\bm{\mathbf{v}},\mu)&=&(\mu,p_{i}),&\ \forall\mu\in{Q},\end{array}\right.$
(28)

where

$\left\\{\begin{aligned}
\tilde{a}({\bm{\mathbf{\varphi}}},{\bm{\mathbf{\psi}}})&=\int_{\Omega}{{\alpha}\nabla{\bm{\mathbf{\varphi}}}:\nabla{\overline{\bm{\mathbf{\psi}}}}}{\,{d}x}-i\int_{\Sigma}{{\bm{\mathbf{\varphi}}}\cdot{\overline{\bm{\mathbf{\psi}}}}}{\,{d}\sigma},\quad\text{with}\
{\bm{\mathbf{\varphi}}},{\bm{\mathbf{\psi}}}\in{V_{\Gamma}},\\\
\tilde{F}(\bm{\mathbf{\psi}})&=\int_{\Omega}{\bm{\mathbf{u}}_{i}\cdot\overline{\bm{\mathbf{\psi}}}}{\,{d}x},\quad\text{with
$\bm{\mathbf{u}}_{i}\in H^{1}(\Omega)^{d}$ and
${\bm{\mathbf{\psi}}}\in{V_{\Gamma}}$}.\end{aligned}\right.$ (29)

The demonstration of the well-posedness of (28) follows a standard procedure.
It mirrors the proof outlined in Proposition 2.3, anchoring on the distinct
solvability of (12) within the space $V_{\Gamma}\times Q$. Notably, it can be
established that the variational problem satisfies an inf-sup condition.

At this juncture, we note that for $\Omega$ of class $\mathcal{C}^{k+1,1}$,
where $k\in\mathbb{N}_{0}$, and assuming $(\bm{\mathbf{f}},\bm{\mathbf{g}})\in
H^{k+3/2}(\Sigma)^{d}\times H^{k+1/2}(\Sigma)^{d}$, it can be demonstrated
that the weak solution $(\bm{\mathbf{u}},p)\in V_{\Gamma}\times Q$ of the
variational problem (12) also belongs to
$\bm{\mathbf{\mathsf{H}}}^{k+2}(\Omega)^{d}\times\bm{\mathbf{\mathsf{H}}}^{k+1}(\Omega)$
– specifically, $\bm{\mathbf{u}}_{r},\bm{\mathbf{u}}_{i}\in
H^{k+2}(\Omega)^{d}\cap H^{1}_{\Gamma,0}(\Omega)^{d}$ and $p_{r},p_{i}\in
H^{k+1}(\Omega)$. Consequently, the weak solution $(\bm{\mathbf{v}},q)$ of
problem (28) is not only in $V_{\Gamma}\times Q$ but also lies within
$\bm{\mathbf{\mathsf{H}}}^{k+2}(\Omega)^{d}\times\bm{\mathbf{\mathsf{H}}}^{k+1}(\Omega)$.

Next, we introduce the sesquilinear forms $a^{t},\tilde{a}^{t}\in
V_{\Gamma}\times V_{\Gamma}\to\mathbb{R}$ and linear forms,
$b^{t}:V_{\Gamma}\times Q\to\mathbb{R}$ and
$F^{t},\tilde{F}^{t}:V_{\Gamma}\to\mathbb{R}$:

$\left\\{\begin{aligned}
a^{t}({\bm{\mathbf{\varphi}}},{\bm{\mathbf{\psi}}})&=\int_{\Omega}{{\alpha}A_{t}\nabla{\bm{\mathbf{\varphi}}}:\nabla{\overline{\bm{\mathbf{\psi}}}}}{\,{d}x}+i\int_{\Sigma}{{\bm{\mathbf{\varphi}}}\cdot{\overline{\bm{\mathbf{\psi}}}}}{\,{d}\sigma},\qquad
b^{t}(\bm{\mathbf{\varphi}},\lambda)=-\int_{\Omega}{I_{t}\bar{\lambda}(M_{t}^{\top}:\nabla\bm{\mathbf{\varphi}})}{\,{d}x},\\\
F^{t}(\bm{\mathbf{\varphi}})&=\int_{\Sigma}{B_{t}(\bm{\mathbf{g}}^{t}+i\bm{\mathbf{f}}^{t})\cdot\overline{\bm{\mathbf{\varphi}}}}{\,{d}\sigma}\equiv
F(\bm{\mathbf{\varphi}})=\int_{\Sigma}{(\bm{\mathbf{g}}+i\bm{\mathbf{f}})\cdot\overline{\bm{\mathbf{\varphi}}}}{\,{d}\sigma},\\\
\tilde{a}^{t}({\bm{\mathbf{\varphi}}},{\bm{\mathbf{\psi}}})&=\int_{\Omega}{{\alpha}A_{t}\nabla{\bm{\mathbf{\varphi}}}:\nabla{\overline{\bm{\mathbf{\psi}}}}}{\,{d}x}-i\int_{\Sigma}{{\bm{\mathbf{\varphi}}}\cdot{\overline{\bm{\mathbf{\psi}}}}}{\,{d}\sigma},\qquad\tilde{F}^{t}(\bm{\mathbf{\varphi}})=\int_{\Omega}{I_{t}\bm{\mathbf{u}}_{i}\cdot\overline{\bm{\mathbf{\varphi}}}}{\,{d}x},\end{aligned}\right.$
(30)

where $\bm{\mathbf{u}}_{i}$ is the imaginary part of the velocity solution to
the state system (8).

In the above, notice that
$F^{t}(\bm{\mathbf{\varphi}})=F(\bm{\mathbf{\varphi}})$ holds true for any
${\bm{\mathbf{\varphi}}}\in{V_{\Gamma}}$ because $\Sigma_{t}\equiv\Sigma$ for
all $t>0$, and thus $T_{t}(x)=x$ for all $x\in\Sigma_{t}$, for all $t>0$. This
is clear from the fact that $B_{t}=1$ on $\Sigma$. As a consequence, we have
$\bm{\mathbf{f}}^{t}=\bm{\mathbf{f}}$ and
$\bm{\mathbf{g}}^{t}=\bm{\mathbf{g}}$ on $\Sigma$. Further below, we will use
the given notations to express the respective variational formulations of the
transformed versions of the state and adjoint state systems (8) and (24).

Before deriving the shape gradient of $J$, we need to establish some key
lemmas.

######  Lemma 3.3.

For $t\in\mathcal{I}$, the sesquilinear forms $a^{t}$ and $\tilde{a}^{t}$
defined on $X\times X$ are bounded and coercive on $X\times X$. In addition,
the linear forms $F^{t}(\bm{\mathbf{\varphi}})$ and
$\tilde{F}^{t}(\bm{\mathbf{\varphi}})$ are also bounded. Moreover, the
bilinear form $b$ satisfies the condition that there is a constant
$\beta_{1}>0$ such that

$\inf_{\begin{subarray}{c}{\lambda}\in{Q}\\\ {\lambda}\neq
0\end{subarray}}\sup_{\begin{subarray}{c}\bm{\mathbf{\varphi}}\in{V_{\Gamma}}\\\
\bm{\mathbf{\varphi}}\neq\bm{\mathbf{0}}\end{subarray}}\frac{b^{t}(\bm{\mathbf{\varphi}},{\lambda})}{{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}{|\kern-1.07639pt|\kern-1.07639pt|{\lambda}|\kern-1.07639pt|\kern-1.07639pt|}_{Q}}\geqslant\beta_{1}.$
(31)

###### Proof.

The proof follows a parallel line of reasoning as presented in [RN24, Appx.
A], so we omit it. However, it is important to note that one must consider the
properties of $A_{t}$ and $B_{t}$ given in (A1), along with the bounds
specified in (18), to derive the desired inequality. ∎

The following two lemmas address the well-posedness of the transported
perturbed version of (12) and the (uniform) boundedness of its solution,
respectively. Regarding the result, it is important to note that the
compatibility condition
$\int_{\Sigma_{t}}{\bm{\mathbf{u}}_{t}\cdot\bm{\mathbf{n}}_{t}}{\,{d}\sigma_{t}}=0$,
assumed hereafter, is equivalent to (2) because $\Sigma$ is invariant.

######  Lemma 3.4.

The pair
$(\bm{\mathbf{u}}^{t},p^{t})=(\bm{\mathbf{u}}_{r}^{t}+i\bm{\mathbf{u}}_{i}^{t},p_{r}^{t}+ip_{i}^{t})$
uniquely solves in $V_{\Gamma}\times Q$ the variational problem

$\left\\{\begin{array}[]{rcll}a^{t}({\bm{\mathbf{u}}^{t}},{\bm{\mathbf{\varphi}}})+b^{t}(\bm{\mathbf{\varphi}},p^{t})&=&F^{t}(\bm{\mathbf{\varphi}}),&\
\forall{\bm{\mathbf{\varphi}}}\in{V_{\Gamma}},\\\
b^{t}(\bm{\mathbf{u}}^{t},\lambda)&=&0,&\
\forall\lambda\in{Q}.\end{array}\right.$ (32)

######  Lemma 3.5.

For $t\in\mathcal{I}$, the solution pair $(\bm{\mathbf{u}}^{t},p^{t})$ of (32)
are uniformly bounded in $X\times Q$. More precisely, for all
$t\in\mathcal{I}$, we have
${|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X},\
{|\kern-1.07639pt|\kern-1.07639pt|p^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{Q}\lesssim\left\|{(\bm{\mathbf{f}},\bm{\mathbf{g}})}\right\|_{1/2,\Sigma}$.

To complete our preparations, we next issue the Hölder continuity of the state
variables $(\bm{\mathbf{u}}^{t},p^{t})$ with respect to $t$.

######  Lemma 3.6.

Let $(\bm{\mathbf{u}},p)\in V_{\Gamma}\times Q$ be the solution of (12). Then,
the following limit holds

$\lim_{t\to
0^{+}}\frac{1}{\sqrt{t}}\left({|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}+{|\kern-1.07639pt|\kern-1.07639pt|p^{t}-p|\kern-1.07639pt|\kern-1.07639pt|}_{Q}\right)=0,$

where $({\bm{\mathbf{u}}^{t}},p^{t})\in V_{\Gamma}\times Q$ solves (32), for
$t\in\mathcal{I}$.

The proofs of Lemma 3.4, Lemma 3.5, and Lemma 3.6, which closely follow the
techniques in [RN24], are deferred to Appendix A.

Additionally, by the higher regularity property of $\bm{\mathbf{u}}^{t}$,
$\bm{\mathbf{u}}$, $p^{t}$, and $p$, we easily obtain the following corollary
of Lemma 3.6.

######  Corollary 3.7.

The state $(\bm{\mathbf{u}},p)$ belongs to
$(\bm{\mathbf{\mathsf{H}}}^{2}(\Omega)^{d}\times\bm{\mathbf{\mathsf{H}}}^{1}(\Omega))\cap(V_{\Gamma}\cap
Q)$. Moreover, the following limit holds

$\lim_{t\to
0^{+}}\frac{1}{\sqrt{t}}\left({|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}|\kern-1.07639pt|\kern-1.07639pt|}_{\bm{\mathbf{\mathsf{H}}}^{2}(\Omega)^{d}}+{|\kern-1.07639pt|\kern-1.07639pt|p^{t}-p|\kern-1.07639pt|\kern-1.07639pt|}_{\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)}\right)=0,$

where
$({\bm{\mathbf{u}}^{t}},p^{t})\in(\bm{\mathbf{\mathsf{H}}}^{2}(\Omega)^{d}\times\bm{\mathbf{\mathsf{H}}}^{1}(\Omega))\cap(V_{\Gamma}\cap
Q)$ solves (32), for $t\in\mathcal{I}$.

We are almost ready to present the proof of Theorem 3.1 without relying on the
shape derivative of the state. We only need to establish additional groundwork
by introducing a set of auxiliary results, which were previously established
in [DZ11] and [RN24]. The lemmas derived from these results encompass various
identities crucial for our purpose and serve as key components in deriving the
boundary expression for the shape derivative of $J$, aligning with the
Hadamard-Zolésio structure theorem [DZ11, p. 479].

######  Lemma 3.8 ([DZ11]).

For $\varphi$, $\psi\in H^{2}(\Omega)$ such that $\varphi=0$ and $\psi=0$ on
$\Sigma$ and $\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}$, we have

$\displaystyle-\int_{\Omega}{A\nabla\varphi\cdot\nabla\overline{\psi}}{\,{d}x}$
$\displaystyle=-\int_{\Omega}{(\Delta\varphi)\bm{\mathbf{\theta}}\cdot\nabla\overline{\psi}}{\,{d}x}-\int_{\Omega}{(\Delta\overline{\psi})\bm{\mathbf{\theta}}\cdot\nabla\varphi}{\,{d}x}+\int_{\Gamma}{\partial_{\bm{\mathbf{n}}}{\varphi}(\bm{\mathbf{\theta}}\cdot\nabla\overline{\psi})}{\,{d}\sigma}$
$\displaystyle\
\qquad+\int_{\Gamma}{\partial_{\bm{\mathbf{n}}}{\overline{\psi}}(\bm{\mathbf{\theta}}\cdot\nabla\varphi)}{\,{d}\sigma}-\int_{\Gamma}{(\nabla\overline{\psi}\cdot\nabla\varphi)\theta_{n}}{\,{d}\sigma}.$

######  Lemma 3.9 ([RN24, Lem. 3.11]).

The state solution $\bm{\mathbf{u}}\in X$ of (12) satisfies the equation

$\lim_{t\to 0}b(\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}},\lambda)=-\lim_{t\to
0}\int_{\Omega}{\overline{\lambda}\nabla\cdot\left(\frac{\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}}{t}\right)}{\,{d}x}=-\int_{\Omega}{\overline{\lambda}(D\bm{\mathbf{\theta}}:\nabla\bm{\mathbf{u}})}{\,{d}x},\quad\forall\lambda\in{Q}.$

######  Lemma 3.10 ([RN24, Lem. 3.13]).

For sufficiently smooth $\bm{\mathbf{\varphi}}$ and $\bm{\mathbf{\theta}}$, it
holds that

$\operatorname{div}(\nabla\bm{\mathbf{\varphi}}^{\top}\bm{\mathbf{\theta}})=(\bm{\mathbf{\theta}}\cdot\nabla)(\nabla\cdot\bm{\mathbf{\varphi}})+D\bm{\mathbf{\theta}}:\nabla\bm{\mathbf{\varphi}}=(\bm{\mathbf{\theta}}\cdot\nabla)(\nabla\cdot\bm{\mathbf{\varphi}})+\nabla\bm{\mathbf{\varphi}}^{\top}:D\bm{\mathbf{\theta}}^{\top}.$
(33)

The next lemma is derived from the arguments presented in the proof of Theorem
3.4 in [RN24], in conjunction with the fact that
$\bm{\mathbf{v}}=\bm{\mathbf{0}}$ on $\Gamma$.

######  Lemma 3.11.

Let $p$ be the pressure solution of the state system (8) and $\bm{\mathbf{v}}$
be the velocity solution of the adjoint problem (24). Then, for any
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}$, the following identity holds

$\int_{\Omega}{\left[p(\nabla\cdot\bm{\mathbf{\theta}})(\nabla\cdot\overline{\bm{\mathbf{v}}})+p(\bm{\mathbf{\theta}}\cdot\nabla)(\nabla\cdot\overline{\bm{\mathbf{v}}})\right]}{\,{d}x}=\int_{\Gamma}{\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}\cdot
p\bm{\mathbf{n}}\theta_{n}}}{\,{d}\sigma}+\int_{\Omega}{(\bm{\mathbf{\theta}}\cdot\nabla
p)p_{i}}{\,{d}x}.$

We now provide the proof of our main result, Theorem 3.1.

###### Proof of Theorem 3.1 .

The proof essentially proceeds in two steps. First, we evaluate the limit
$\lim_{t\to 0}[J(\Omega_{t})-J(\Omega)]/t$. Then, using the regularity of the
domain as well as the state and adjoint variables, we characterized the
boundary integral expression for the computed limit.

We begin by applying a change-of-variables formula (see [DZ11, eq. (4.2), p.
482]):

$\int_{\Omega_{t}}{\varphi_{t}}{\,{d}x_{t}}=\int_{\Omega}{\varphi_{t}\circ
T_{t}I_{t}}{\,{d}x}=\int_{\Omega}{\varphi^{t}I_{t}}{\,{d}x},$ (34)

which holds for $\varphi_{t}\in L^{1}(\Omega_{t})$ combined with the identity
$\eta^{2}-\zeta^{2}=(\eta-\zeta)^{2}+2\zeta(\eta-\zeta)$ to obtain the
following sequence of calculations:

$\displaystyle J(\Omega_{t})-J(\Omega)$
$\displaystyle=\frac{1}{2}\int_{\Omega}{{\mathfrak{i}}_{t}(|{\bm{\mathbf{u}}_{i}^{t}}|^{2}-|{\bm{\mathbf{u}}_{i}}|^{2})}{\,{d}x}+\frac{1}{2}\int_{\Omega}{{\mathfrak{i}}_{t}|{\bm{\mathbf{u}}_{i}}|^{2}}{\,{d}x}+\frac{1}{2}\int_{\Omega}{|{\bm{\mathbf{y}}_{i}^{t}}|^{2}}{\,{d}x}+\int_{\Omega}{\bm{\mathbf{y}}_{i}^{t}\cdot\bm{\mathbf{u}}_{i}}{\,{d}x}$
$\displaystyle\quad+\frac{1}{2}\int_{\Omega}{{\mathfrak{i}}_{t}(|{p_{i}^{t}}|^{2}-|{p_{i}}|^{2})}{\,{d}x}+\frac{1}{2}\int_{\Omega}{{\mathfrak{i}}_{t}|{p_{i}}|^{2}}{\,{d}x}+\frac{1}{2}\int_{\Omega}{|{r_{i}^{t}}|^{2}}{\,{d}x}+\int_{\Omega}{r_{i}^{t}p_{i}}{\,{d}x}$
$\displaystyle=:\sum_{i=1}^{8}{J}_{i}(t),$

where $\bm{\mathbf{y}}_{i}^{t}=\bm{\mathbf{u}}_{i}^{t}-\bm{\mathbf{u}}_{i}$
and $r_{i}^{t}=p_{i}^{t}-p_{i}$. In view of (19)1 and using Lemma 3.6, we
easily infer that
$\dot{J}_{1}(0)=\dot{J}_{3}(0)=\dot{J}_{5}(0)=\dot{J}_{7}(0)=0$. The limits in
(19)1 also reveal that

$\dot{J}_{2}(0)+\dot{J}_{6}(0)=\frac{1}{2}\int_{\Omega}{\operatorname{div}\bm{\mathbf{\theta}}(|{\bm{\mathbf{u}}_{i}}|^{2}+|{p_{i}}|^{2})}{\,{d}x}.$

We expand the above integral expression through the identity
$\operatorname{div}(\psi\bm{\mathbf{\varphi}})=\psi\operatorname{div}\bm{\mathbf{\varphi}}+\bm{\mathbf{\varphi}}\cdot\nabla\psi$,
and then apply Green’s theorem to obtain

$\displaystyle\dot{J}_{2}(0)+\dot{J}_{6}(0)$
$\displaystyle=-\int_{\Omega}{(\bm{\mathbf{u}}_{i}\cdot\nabla\bm{\mathbf{u}}_{i}^{\top}\bm{\mathbf{\theta}}+p_{i}\bm{\mathbf{\theta}}\cdot\nabla
p_{i})}{\,{d}x}+\frac{1}{2}\int_{\Gamma}{|{p_{i}}|^{2}\theta_{n}}{\,{d}\sigma}.$
(35)

The calculations for the remaining two expressions, $\dot{J}{4}(0)$ and
$\dot{J}{8}(0)$, demand more effort and necessitate the application of the
adjoint system (28). Considering that
$\bm{\mathbf{y}}^{t}=\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}\in V_{\Gamma}$ and
$r^{t}=p^{t}-p\in Q$, we can express

$\displaystyle{J}_{4}(t)+{J}_{8}(t)$
$\displaystyle=\Im\left\\{\int_{\Omega}{{\alpha}\nabla\overline{\bm{\mathbf{v}}}\cdot\nabla\bm{\mathbf{y}}^{t}}{\,{d}x}+i\int_{\Sigma}{\overline{\bm{\mathbf{v}}}\cdot\bm{\mathbf{y}}^{t}}{\,{d}\sigma}-\int_{\Omega}{\overline{q}\operatorname{div}\bm{\mathbf{y}}^{t}}{\,{d}x}-\int_{\Omega}{r^{t}\operatorname{div}\bm{\mathbf{v}}}{\,{d}x}\right\\}$
$\displaystyle\equiv\Im\left\\{a(\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}},\bm{\mathbf{v}})+b(\bm{\mathbf{v}},p_{t}-p)+b(\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}},q)\right\\}$
$\displaystyle=\Im\\{\Phi^{t}(\bm{\mathbf{v}})+b(\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}},q)\\},$

where the latter equation follows from (A4) and $\Phi^{t}(\bm{\mathbf{v}})$ is
given by (see (20) for notations)

$\Phi^{t}(\bm{\mathbf{\varphi}})=-\int_{\Omega}{{\alpha}\mathfrak{a}_{t}\nabla{\bm{\mathbf{u}}^{t}}:\nabla{\overline{\bm{\mathbf{\varphi}}}}}{\,{d}x}+\int_{\Omega}{{\mathfrak{i}}_{t}p^{t}(M_{t}^{\top}:\nabla\overline{\bm{\mathbf{\varphi}}})}{\,{d}x}+\int_{\Omega}{p^{t}(\mathfrak{m}_{t}:\nabla\overline{\bm{\mathbf{\varphi}}})}{\,{d}x},$
(36)

with $\bm{\mathbf{\varphi}}=\bm{\mathbf{v}}\in V_{\Gamma}$. Using (19), (22),
Lemma 3.9, and identity (33), we obtain the following limit

$\displaystyle\dot{J}_{4}(0)+\dot{J}_{8}(0)$ $\displaystyle=\lim_{t\to
0}\frac{1}{t}\left[\Im\\{\Phi^{t}(\bm{\mathbf{v}})+b(\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}},q)\\}\right]$
(37)
$\displaystyle=\Im\Big{\\{}-\int_{\Omega}{{\alpha}A\nabla{\bm{\mathbf{u}}}:\nabla{\overline{\bm{\mathbf{v}}}}}{\,{d}x}+\int_{\Omega}{(\nabla\cdot\bm{\mathbf{\theta}})p(\nabla\cdot\overline{\bm{\mathbf{v}}})}{\,{d}x}+\int_{\Omega}{p[(-D\bm{\mathbf{\theta}}):\nabla\overline{\bm{\mathbf{v}}}]}{\,{d}x}$
$\displaystyle\quad\qquad-\int_{\Omega}{\overline{q}(D\bm{\mathbf{\theta}}:\nabla\bm{\mathbf{u}})}{\,{d}x}\Big{\\}}$
$\displaystyle\quad=:\Im\left\\{K_{1}+K_{2}+K_{3}+K_{4}\right\\}.$

To write this expression as a boundary integral, we invoke the regularity
assumptions on $\bm{\mathbf{f}}$, $\bm{\mathbf{g}}$, and $\Omega$. Let us note
that $\bm{\mathbf{u}}=\bm{\mathbf{v}}=\bm{\mathbf{0}}$ on $\Gamma$ which
implies that
$\nabla{\bm{\mathbf{u}}}=\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\bm{\mathbf{n}}$
and
$\nabla{\overline{\bm{\mathbf{v}}}}=\partial_{\bm{\mathbf{n}}}{\overline{\bm{\mathbf{v}}}}\bm{\mathbf{n}}$
on $\Gamma$. Thus, through Lemma 3.8, the integral $K_{1}$ can be expanded as
follows

$\displaystyle-\int_{\Omega}{{\alpha}A\nabla{\bm{\mathbf{u}}}:\nabla{\overline{\bm{\mathbf{v}}}}}{\,{d}x}$
$\displaystyle\quad=\sum_{j=1}^{d}\Bigg{\\{}\int_{\Omega}{(-{\alpha}\Delta{u}_{j})\bm{\mathbf{\theta}}\cdot\nabla\overline{{v}}_{j}}{\,{d}x}-\int_{\Omega}{({\alpha}\Delta\overline{{v}}_{j})\bm{\mathbf{\theta}}\cdot\nabla{u}_{j}}{\,{d}x}+\int_{\Gamma}{\mu\partial_{\bm{\mathbf{n}}}{{u}_{j}}(\bm{\mathbf{\theta}}\cdot\nabla\overline{{v}}_{j})}{\,{d}\sigma}$
$\displaystyle\quad\quad\qquad+\int_{\Gamma}{\mu\partial_{\bm{\mathbf{n}}}{\overline{{v}}_{j}}(\bm{\mathbf{\theta}}\cdot\nabla{u}_{j})}{\,{d}\sigma}-\int_{\Gamma}{\mu(\nabla\overline{{v}}_{j}\cdot\nabla{u}_{j})\theta_{n}}{\,{d}\sigma}\Bigg{\\}}$
$\displaystyle\quad=\sum_{j=1}^{d}\Bigg{\\{}\int_{\Omega}{(-{\alpha}\Delta{u}_{j})\bm{\mathbf{\theta}}\cdot\nabla\overline{{v}}_{j}}{\,{d}x}-\int_{\Omega}{({\alpha}\Delta\overline{{v}}_{j})\bm{\mathbf{\theta}}\cdot\nabla{u}_{j}}{\,{d}x}+\int_{\Gamma}{\mu\partial_{\bm{\mathbf{n}}}{\overline{{v}}_{j}}(\bm{\mathbf{\theta}}\cdot\nabla{u}_{j})}{\,{d}\sigma}\Bigg{\\}}$
$\displaystyle\quad=\int_{\Omega}{p\operatorname{div}{(\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}})}}{\,{d}x}-\int_{\Gamma}{p\bm{\mathbf{n}}\cdot\nabla\overline{\bm{\mathbf{v}}}^{\top}\bm{\mathbf{\theta}}}{\,{d}\sigma}+\int_{\Omega}{\bm{\mathbf{u}}_{i}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}}{\,{d}x}+\int_{\Omega}{\overline{q}\operatorname{div}{(\nabla{\bm{\mathbf{u}}}^{\top}\bm{\mathbf{\theta}})}}{\,{d}x}$
$\displaystyle\qquad\qquad+\int_{\Gamma}{\overline{\sigma(\bm{\mathbf{v}},q)}\bm{\mathbf{n}}\cdot\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\theta_{n}}{\,{d}\sigma}$
$\displaystyle\quad=:K_{11}+K_{12}+K_{13}+K_{14}+K_{15}.$

In light of Lemma 3.10 (considering $\nabla\cdot{\bm{\mathbf{u}}}=0$ in
$\Omega$) and subsequent application of Lemma 3.11, we determine that
$K_{14}+K_{4}=0$ and
$K_{11}+K_{12}+K_{2}+K_{3}=\int_{\Omega}{(\bm{\mathbf{\theta}}\cdot\nabla{p})p_{i}}{\,{d}x}$.
Thus, from equation (37), we have the limit

$\dot{J}_{4}(0)+\dot{J}_{8}(0)=\Im\left\\{\int_{\Omega}{\left[\bm{\mathbf{u}}_{i}\cdot\nabla\bm{\mathbf{u}}^{\top}\bm{\mathbf{\theta}}+(\bm{\mathbf{\theta}}\cdot\nabla{p})p_{i}\right]}{\,{d}x}\right\\}+\Im\left\\{\int_{\Gamma}{\overline{\sigma(\bm{\mathbf{v}},q)}\bm{\mathbf{n}}\cdot\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}\theta_{n}}{\,{d}\sigma}\right\\}.$

Combining this expression with $\dot{J}_{2}(0)+\dot{J}_{6}(0)$ from (35)
yields the desired expression (23), concluding the proof of the theorem. ∎

######  Remark 3.12.

In the proof of Theorem 3.1, we only require $\Omega$ to be
$\mathcal{C}^{1,1}$ regular, in addition to the condition that
$(\bm{\mathbf{f}},\bm{\mathbf{g}})\in H^{3/2}(\Sigma)^{d}\times
H^{1/2}(\Sigma)^{d}$. These assumptions enable us to attain higher regularity
for the state and adjoint variables, facilitating the application of Lemma 3.8
to express the shape derivative of $J$ in a boundary integral form. It
actually suffices to assume that $\Sigma$ is only Lipschitz smooth to obtain
the shape derivative of $J$. However, for the sake of simplicity in our
arguments, we have adopted the stated regularity assumptions.

The following result can be drawn easily from (23), (24), and the equivalence
between the shape optimization problem and the original inverse problem.

######  Corollary 3.13 (Necessary optimality condition).

Let the domain $\Omega^{\ast}$ be such that
$(\bm{\mathbf{u}},p)=(\bm{\mathbf{u}}(\Omega^{\ast}),p(\Omega^{\ast}))$
satisfies (1), i.e., we have
$-p\bm{\mathbf{n}}+{\alpha}\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{u}}}=\bm{\mathbf{g}}$
and $\bm{\mathbf{u}}=\bm{\mathbf{f}}$ on $\Sigma^{\ast}$, or equivalently,

$\bm{\mathbf{u}}_{i}=\bm{\mathbf{0}}\quad\text{and}\quad p_{i}=0\quad\text{on
$\Omega^{\ast}$},$ (38)

with $(\bm{\mathbf{u}},p)$ satisfying (8). Then, the domain $\Omega^{\ast}$ is
a stationary solution for the shape optimization problem $J(\Omega)\to\inf$
which is subject to the state problem (8); i.e., it fulfills the necessary
optimality condition

${d}J(\Omega^{\ast})[\bm{\mathbf{\theta}}]=0,\quad\text{for all
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}$}.$ (39)

###### Proof.

Using the assumption that $\bm{\mathbf{u}}_{i}=\bm{\mathbf{0}}$ and $p_{i}=0$
on $\Omega^{\ast}$, one obtains – in view of (24) – that
$\bm{\mathbf{v}}_{i}=\bm{\mathbf{0}}$ and $q_{i}=0$ on $\Omega^{\ast}$. Thus,
it follows that ${G}=0$ on $\Gamma^{\ast}$ which therefore gives us the
conclusion ${d}J(\Omega^{\ast})[\bm{\mathbf{\theta}}]=0$, for any
$\bm{\mathbf{\theta}}\in\bm{\mathbf{\Theta}}$. ∎

With Theorem 3.1 and Corollary 3.13 at hand, we can devise a numerical
procedure to solve Problem 2.4. We will present this procedure in the next
section, along with various numerical examples to test the proposed method.

##  4 Numerical experiments

To numerically solve the optimization problem $J(\Omega)\to\inf$, we will
apply a preconditioned shape-gradient-based descent method based on FEM.

######  Remark 4.1.

Newton-type methods, as shown in [EH05, ADEK07, RA20], could be developed for
the numerical procedure using shape Hessian information. However, we refrain
from applying the method for two main reasons. Firstly, computing second-order
shape derivatives is challenging, involving substantial effort, preparations,
and technical intricacies; see, e.g., [ADK08]. In the computation of the shape
Hessian of $J$, two additional systems of complex PDEs, which correspond to
the new adjoints, must be introduced.The second reason is primarily numerical.
In general, the implementation and computational effort required for computing
the Hessian are not justified by a significant reduction in the number of
iterations. Indeed, it is important to note that evaluating the shape Hessian
requires the solution of (at least) four variational problems (two for the
states and another two for the adjoints), making it prohibitively costly.
Despite these drawbacks, we defer the examination of the shape Hessian of $J$
for future work.

###  4.1 Numerical algorithm

Our main algorithm follows a standard approximation procedure (see, e.g.,
[CDKT13, RN24]), the important details of which we provide as follows.

Choice of descent direction. We will employ a Riesz representation of the
shape gradient $G$ to avoid undesired oscillations on the unknown boundary,
which could result in instabilities in the approximation. To do this, we fix
$\eta\in(0,1]$ and generate a smooth descent direction for the cost function
$J$ by seeking a vector $\bm{\mathbf{\theta}}\in H_{\Sigma,0}^{1}(\Omega)^{d}$
that solves the variational equation

$\eta\int_{\Omega}{\nabla\bm{\mathbf{\theta}}:\nabla\bm{\mathbf{\varphi}}}{\,{d}x}+(1-\eta)\int_{\Gamma}{\nabla_{\Gamma}\bm{\mathbf{\theta}}:\nabla_{\Gamma}\bm{\mathbf{\varphi}}}{\,{d}\sigma}=-\int_{\Gamma}{{{G}}\bm{\mathbf{n}}\cdot\bm{\mathbf{\varphi}}}{\,{d}\sigma},\
\text{for all $\bm{\mathbf{\varphi}}\in H_{\Sigma,0}^{1}(\Omega)^{d}$}.$ (40)

Here, we choose $\eta=0.5$. In this way we obtain a Sobolev gradient [Neu97]
representation $\bm{\mathbf{\theta}}$ of $G$. More importantly, this approach
produces a smoothed preconditioned extension of $-{G}\bm{\mathbf{n}}$ over the
entire domain $\Omega$. Such extensions allow us to deform the discretized
computational domain by moving the (movable) nodes of the computational mesh
and not only the boundary nodes that are subject for perturbation. The
boundary integral appearing on the right side of (40) corresponds to the
Laplace-Beltrami operator on $\Gamma$. Its purpose is to ensure an additional
level of smoothness for the descent direction $\bm{\mathbf{\theta}}$ within
$\Gamma$. Further discussion about discrete gradient flows for shape
optimization problems are issued in [DMNV07].

######  Remark 4.2.

We reiterate here that our approach differs completely from that applied in
[CDKT13]. In the cited work, the authors employed parametric regularization to
address the ill-posedness of the problem, utilizing a parametric model of
shape variations for numerical realization. This process involved restricting
the domains to star-shaped configurations and using polar coordinates for
parametrization, resulting in an explicit computation of
$\bm{\mathbf{\theta}}$.

To compute the $k$th domain $\Omega^{k}$, we carry out the following
procedures:

1\. Initilization

Choose an initial guess $\Omega^{0}$.

2\. Iteration

For $k=0,1,2,\ldots$, do the following:

  1. 2.1

Solve the state’s and adjoint’s variational equations on the current domain
$\Omega^{k}$.

  2. 2.2

Choose $t^{k}>0$, and compute the descent vector $\bm{\mathbf{\theta}}^{k}$
using (40) in $\Omega^{k}$.

  3. 2.3

Update the current domain by setting
$\Omega^{k+1}:=\\{x+t^{k}\bm{\mathbf{\theta}}^{k}(x)\in\mathbb{R}^{d}\mid
x\in\Omega^{k}\\}$.

3\. Stop Test

Repeat Iteration until convergence.

Step-size computation and stopping condition. In Step 2.2, $t^{k}$ is computed
using a backtracking line search procedure inspired by [RA20, p. 281] with the
formula $t^{k}=\mu
J(\Omega^{k})/|\bm{\mathbf{\theta}}^{k}|^{2}_{\mathbf{H}^{1}(\Omega^{k})}$ at
each iteration, where $\mu>0$ is a scaling factor. This step size is further
reduced to avoid reversed triangles within the mesh after the update.
Meanwhile, we terminate the algorithm as soon as the condition that
$|{J^{k-1}-J^{k}}|<\varepsilon_{J}$ or $t^{k}<\varepsilon_{T}$ is satisfied
where $\varepsilon_{J},\varepsilon_{T}>0$ are sufficiently small real numbers,
or the algorithm reached the maximum allowable number of iterations. In all
experiments, we set $\mu=0.5$, $\varepsilon_{J}=\varepsilon_{T}=10^{-12}$, and
the maximum number of iterations to $300$.

###  4.2 Forward problem

We let $\alpha=1.0$ and consider the unit circle as the shape of the medium.
For the input data, we will consider
$\bm{\mathbf{g}}=(\sin\theta,-\cos\theta)^{\top}$, where $\theta\in[0,2\pi)$
represents the prescribed flux. To obtain additional boundary measurements on
the accessible part $\Sigma$, we will use synthetic data. Specifically, given
a prescribed Neumann flux $\bm{\mathbf{g}}$, we generate the Dirichlet data
$\bm{\mathbf{f}}$ on the accessible boundary $\Sigma$ by numerically solving
the well-posed PDE problem (6) using the finite element method. To avoid
‘inverse crimes’ (as discussed by Colton and Kress [CK19, p. 179]) in
producing the measurements, we generate the synthetic data with a different
numerical scheme. This involves using a larger number of discretization points
(with finer meshes of uniform mesh width $h^{\ast}=1/256$) and applying
${P}{3}/{P}{2}$ finite element basis functions in the FreeFem++ [Hec12] code
than in the inversion process. For experiments with noise-contaminated data
with a specified noise level $\delta$, we perturb the data with Gaussian
noise.

###  4.3 Inversion procedure

In the inversion procedure, all variational problems are solved using Hood-
Taylor (${P}{2}/{P}{1}$) finite elements, and we discretize the domain with
coarse meshes. Specifically, in all experiments, the initial domain is
discretized with uniform mesh size having $70$ nodes on the interior boundary
and $100$ nodes on the exterior boundary.

######  Remark 4.3.

We emphasize here that we are, in fact, using lower-degree finite elements in
our numerical simulations, in contrast to the computational setup in [CDKT13],
where the authors used ${P}{4}/{P}{3}$ finite elements for the forward problem
while employing ${P}{3}/{P}{2}$ finite elements in the inversion process.

We will now test our proposed method by solving concrete 2D and 3D examples of
the inverse obstacle problem (1.1).

###  4.4 Numerical tests in 2D

In all plotted figures illustrating the computed shapes, the surface of the
medium is represented by the thicker black solid lines, while the initial
guesses are depicted with black dotted lines. The exact obstacles are
illustrated by red solid lines, whereas the intermediate and final approximate
shapes are represented by black dashed lines (or grey dashed lines), usually
accompanied by markers. We initiate the study with numerical experiments using
exact data before progressing to cases involving noisy data.

Before we proceed with our numerical examples, we recall that, by assumption,
$\Omega$ is of class $\mathcal{C}^{1,1}$. However, in our experiments, we will
also consider cases of unknown obstacles with Lipschitz boundaries (e.g., a
square or an L-shape obstacle). Although these test examples violate the
regularity assumption, we will see below that the algorithm was still able to
produce fair reconstructions of the considered obstacle types.

####  4.4.1 Large convex obstacles

The method performs exceptionally well with large convex obstacles exhibiting
smooth shapes (see the left plot in Figure 1) and functions quite effectively
even with large convex obstacles that have corners (see the right plot in
Figure 1) – as expected.

Figure 1: Shape reconstructions with large convex obstacles

####  4.4.2 Large non-convex obstacle

Meanwhile, when dealing with more complicated shapes featuring pronounced
concavities, the method performs only fairly. For instance, for obstacles
parametrized as follows

$\displaystyle\Gamma^{\ast}_{1}$
$\displaystyle=\left\\{\begin{pmatrix}0.195+0.4(\cos{t}+0.65\cos{2t}),0.55\sin{t}\end{pmatrix}^{\top},\
\forall t\in[0,2\pi)\right\\},$ (C.I) $\displaystyle\Gamma_{2}^{\ast}$
$\displaystyle=\left\\{\begin{pmatrix}0.64\cos{t},0.48\sin{t}(1.8+\cos{(2t)})\end{pmatrix}^{\top},\
\forall t\in[0,2\pi)\right\\},$ (C.II) $\displaystyle\Gamma^{\ast}_{3}$
$\displaystyle=\left\\{\begin{pmatrix}-0.25+\displaystyle\frac{0.6+0.54\cos{t}+0.06\sin{2t}}{1+0.75\cos{t}}\cos{t}\\\\[7.5pt]
0.05+\displaystyle\frac{0.6+0.54\cos{t}+0.06\sin{2t}}{1+0.75\cos{t}}\sin{t}\end{pmatrix},\
\forall t\in[0,2\pi)\right\\},$ (C.III) $\displaystyle\Gamma^{\ast}_{4}$
$\displaystyle=\left\\{0.4(1+0.75\cos(5t+\pi))\begin{pmatrix}\cos{t},\sin{t}\end{pmatrix}^{\top},\
\forall t\in[0,2\pi)\right\\},$ (C.IV)

the reconstructed shapes are only reasonably accurate, as depicted in Figure
2. It seems that obstacles with deep concavities are more challenging to
reconstruct with good accuracy – as expected (see the rightmost plot in Figure
2). Nevertheless, the proposed method accurately predicts the positions of
concavities within the obstacles, as evident in Figure 2.

Figure 2: Shape reconstructions with complicated-shaped obstacles

####  4.4.3 Effect of the choice of initial guess

We also tested various initial guesses to assess their impact on the
reconstruction. Based on our findings, the effect on the reconstruction is
negligible (see Figure 3), unless the initial guess is excessively large or
small. In such cases, it becomes imperative to implement remeshing, involving
the use of coarser meshes, after a certain number of iterations. This
necessity arises because, with continued iterations, the mesh quality
deteriorates, leading to the flattening of some triangles, especially when
employing large step sizes. Addressing this issue is feasible by opting for
smaller step sizes in each iteration.

Figure 3: Effect of the choice of initial guess

####  4.4.4 Noisy data

We present reconstructions using noisy data and emphasize that, in the
presence of noise, reconstructions typically involve regularization. One
common approach is perimeter penalization (see, e.g., [CDKT13, RA18]). The
problem is inherently ill-posed, necessitating the application of
regularization methods for numerical solutions. Introducing a perimeter
functional term, for instance, results in well-posed problems (see, e.g.,
[Dam02, BB05]). However, in our numerical investigation with noisy data, we
refrain from adding any additional regularization term to our method. This
decision is based on our observations from conducted tests, where the proposed
CCBM formulation exhibited a regularizing effect similar to the Kohn-Vogelius
method; see, e.g., [AR22]. The results of our numerical experiments, depicted
in Figures 4 (also referring to Figures 5 and 6), demonstrate reasonable
reconstructions even in the presence of high noise levels. Additionally, the
method generally exhibited the capability to accurately detect concavities
within obstacles. It is worth noting, however, that the presence of noise
compels the approximation to converge towards larger cost values, indicating a
greater deviation of the reconstructed shape from the exact shape, as
expected.

Figure 4: Shape reconstructions with noisy data

####  4.4.5 Influence of the size of the obstacle

We next look at the influence of the obstacle’s size on the reconstruction. We
observed that when attempting to reconstruct a smaller object, the results are
worse, as evident in Figure 5 – as expected. This phenomenon, related to the
obstacle’s size, could be further explained using the explicit calculus of the
shape Hessian, as done, for example, in [ADEK07, CDKT13]. We defer this aspect
for future examination, as the first-order method already yields satisfactory
reconstructions of shapes, even in conditions characterized by a high level of
noise ($\delta=15\%$). See also Remark 4.1.

Figure 5: Influence of the size of the obstacle

####  4.4.6 An obstacle with reentrant corner

We also tested our method for reconstructing an obstacle with a reentrant
corner, such as an L-block (domain see Figure 6). As anticipated, the
reconstruction is much worse compared to obstacles with smooth shapes and
concavities. However, it is noteworthy that the reconstruction accuracy can be
improved by choosing a good initial guess, as shown in Figure 6. We emphasize
that, although the test examples violate the regularity assumption, it appears
that the algorithm was still able to generate fair reconstruction of the type
of inclusion considered.

Figure 6: L-shape obstacle with noisy data ($\delta=15\%$)

####  4.4.7 Detecting more than one object

Let us note here that the main assumption about the domain is that
$\Omega=D\setminus\overline{\omega}$ is connected. This assumption does not
exclude the case of two or more inclusions in $D$. Thus, we aim to numerically
detect two objects: a square and a circle with relatively small sizes. Here,
we apply remeshing every ten iterations to further avoid instabilities during
mesh deformations. As evident in Figure 7, the reconstructions with data
contaminated with a $30\%$ noise level are reasonable, and we point out that
the computational times for the iteration processes in detecting multiple
obstacles are comparable to the case of just one inclusion. For instance, in
the case of a small square centered at the origin shown in Figure 5, it only
requires 291 seconds to complete 300 iterations, while in the present test
case, the average computational time is 300 seconds.

######  Remark 4.4.

In [CDKT13], the authors noted longer computational times for detecting
multiple obstacles compared to a single object – a contrast to our method.
Through remeshing at selected iterates, we achieve fast and fairly accurate
reconstructions of obstacles, even in high-noise conditions, thanks to our use
of coarser meshes and lower-degree finite element methods.

Figure 7: Detecting more than one object with noisy data ($\delta=30\%$)

######  Remark 4.5.

Obviously, the Lagrangian method adopted here to detect objects does not
permit the modification of the topology of the object $\overline{\omega}$.
Thus, in order to detect two or more inclusions, we need to know how many
objects are included in $D$. One solution could be to initialize the algorithm
using the notion of topological gradient (see, for example, [AHJM10, GGM01,
SZ99]), which could provide us with the number of inclusions and their rough
locations, offering initial shapes for our optimization method.

###  4.5 Numerical tests in 3D

In our investigation, we extend the application of our proposed method to 3D
cases. The computational setup closely mirrors that employed in 2D settings.
Specifically, the forward problems are solved using tetrahedra characterized
by smaller mesh widths and volumes, while larger ones are used in the
inversion process. Furthermore, the finite elements utilized differ between
the forward and inverse problems. Specifically, $P3/P2$ finite elements are
employed in the context of the forward problems, whereas $P2/P1$ finite
elements are used for the inverse problems. In the subsequent numerical
experiments, we prescribe the input data as $\bm{\mathbf{g}}=[1,1,1]$ and set
$\mu=1$ during the step size computation. The maximum number of iterations is
set to $30$ for the first two examples, while it is set to $100$ for the third
example. We do not employ adaptive remeshing to reduce the computational cost
associated with the approximation.

####  4.5.1 Large convex obstacle

For large convex shapes, such as an elliptical obstacle shown in Figure 8, the
method works perfectly even in the presence of a high level of noise (e.g.,
$\delta=30\%$). As an illustration, the plot shown in Figure 9 is a
reconstruction of the exact shape in Figure 8 with an initial guess given by a
sphere of radius $0.65$.

Figure 8: Exact ellipsoidal shape obstacle

Figure 9: Identified shape with noisy data ($\delta=30\%$)

####  4.5.2 Large non-convex obstacle

In the case of non-convex obstacles, such as a star-shaped obstacle depicted
in Figure 10, achieving accurate reconstructions appears to be more
challenging, as expected. For instance, in Figure 11, we can only achieve a
reasonably accurate reconstruction of the precise geometry of the obstacle
depicted in Figure 10 at high noise levels with $\delta=30\%$ (and an initial
guess given by a sphere of radius 0.65.). Nevertheless, we are able to
identify some concavities within the shape. Based on our experience, it seems
that regions with pronounced concavities are easier to reconstruct, while
sharper or less pronounced concavities, especially those distant from the
measurement region, prove difficult to detect.

Figure 10: Exact star-shape obstacle

Figure 11: Identified shape with noisy data ($\delta=30\%$)

####  4.5.3 Two small cubic-shaped obstacles

Finally, we evaluate our algorithm using two cubic-shaped obstacles, as
illustrated in Figure 12. Despite violating the regularity assumption on the
unknown interior boundary, this test example demonstrates the algorithm’s
ability to generate a reasonable reconstruction of the inclusion type, even in
the presence of high noise levels ($\delta=30\%$), as depicted in Figure 13.
As expected, it is difficult to accurately reconstruct the edges of the cubic
obstacles, see Figure 14. Nevertheless, it correctly predicts the location of
the two obstacles, reaching convergence in the cost function after
approximately 60 iterations.

Figure 12: Exact geometries

Figure 13: Identified shapes with noisy data ($\delta=30\%$)

Figure 14: Cross-comparisons of exact (blackest color), initial (whitest
color), and identified (cyan color) shapes corresponding to Figures 12 and 13,
along with histories of cost and gradient norms (rightmost plot)

To conclude, we comment that in the first two examples, the algorithm took
approximately $390$ seconds for $30$ iterations with remeshing every five
steps. In the third example, it took around $10,000$ seconds to complete $100$
iterations with remeshing every two steps.

In general, our experiments reveal that smaller obstacles and parts distant
from the measurement region present greater challenges for reconstruction.
Nevertheless, the algorithm effectively detects the obstacle’s location and
provides a reasonable reconstruction of its geometry, even in the presence of
noise in the measurements.

##  5 Conclusion

The coupled complex boundary method is employed to address a shape-inverse
obstacle problem in Stokes fluid flow. The study reformulates the inverse
problem as a shape optimization setup, utilizing least-squares fitting of the
imaginary part of the complex PDEs’ solution in the flow domain. The unknown
inclusion is treated as a nonparametric shape. The shape gradient of the
associated cost functional is rigorously characterized using the rearrangement
method, enabling effective numerical resolution through a Sobolev gradient
descent method in a finite element setting. Numerical results in both two and
three-dimensional cases demonstrate the method’s effectiveness in detecting
non-trivial obstacles, even under high-level noise-contaminated data.

## Acknowledgements

The work of JFTR is supported by the JSPS Grant-in-Aid for Early-Career
Scientists under Japan Grant Number JP23K13012. HN is partially supported by
JSPS Grants-in-Aid for Scientific Research under Grant Numbers JP20KK0058,
JP21H04431, and JP20H01823. JFTR and HN are also partially supported by the
JST CREST Grant Number JPMJCR2014.

## References

  * [ACF+05] C. Alvarez, C. Conca, L. Friz, O. Kavian, and J. H. Ortega.  Identification of immersed obstacles via boundary measurements.  Inverse Problems, 21:1531–1552, 2005.
  * [ADEK07] L. Afraites, M. Dambrine, K. Eppler, and D. Kateb.  Detecting perfectly insulated obstacles by shape optimization techniques of order two.  Discrete Contin. Dyn. Syst. Ser. B, 8(2), 2007.
  * [ADK08] L. Afraites, M. Dambrine, and D. Kateb.  On second order shape optimization methods for electrical impedance tomography.  SIAM J. Control Optim., 47(3):1556–1590, 2008.
  * [Afr22] L. Afraites.  A new coupled complex boundary method (CCBM) for an inverse obstacle problem.  Discrete Contin. Dyn. Syst. Ser. S, 15(1):23–40, 2022.
  * [AHJM10] A. Ben Abda, M. Hassine, M. Jaoua, and M. Masmoudi.  Topological sensitivity analysis for the location of small cavities in stokes flow.  SIAM J. Control Optim., 48:2871–2900, 2009/10.
  * [AMN22] L. Afraites, C. Masnaoui, and M. Nachaoui.  Shape optimization method for an inverse geometric source problem and stability at critical shape.  Discrete Contin. Dyn. Syst. Ser. S, 15(1):1–21, 2022.
  * [AR22] L. Afraites and J. F. T. Rabago.  Shape optimization methods for detecting an unknown boundary with the Robin condition by a single measurement.  Discrete Contin. Dyn. Syst. Ser. S, 2022.
  * [BB05] D. Bucur and G. Buttazzo.  Variational Methods in Shape Optimization Problems, volume 65 of Progr. Nonlinear Differential Equations Appl. Birkhàuser-Verlag, Basel, 2005.
  * [BCD11] M. Badra, F. Caubet, and M. Dambrine.  Detecting an obstacle immersed in a fluid by shape optimization methods.  Math. Meth. Appl. Sci., 21(10):2069–2101, 2011.
  * [BF06] F. Boyer and P. Fabrie.  Éléments d’Analyse Pour L’étude de Quelques Modèles D’écoulements de Fluides Visqueux Incompressibles.  Mathématiques & Applications. Springer-Verlag, Berlin, 2006.
  * [BP13] J. B. Bacani and G. H. Peichl.  On the first-order shape derivative of the Kohn-Vogelius cost functional of the Bernoulli problem.  Abstr. Appl. Anal., 2013:19 pp. Article ID 384320, 2013.
  * [Cau13] F. Caubet.  Instability of an inverse problem for the stationary Navier-Stokes equations.  SIAM J. Control Optim., 51(4):2949–2975, 2013.
  * [CDK13] F. Caubet, M. Dambrine, and D. Kateb.  Shape optimization methods for the inverse obstacle problem with generalized impedance boundary conditions.  Inverse Problems, 29:Art. 115011 (26pp), 2013.
  * [CDKT13] F. Caubet, M. Dambrine, D. Kateb, and C. Z. Timimoun.  A Kohn-Vogelius formulation to detect an obstacle immersed in a fluid.  Inverse Prob. Imaging, 7(1):123–157, 2013.
  * [CGHZ14] X. L. Cheng, R. F. Gong, W. Han, and X. Zheng.  A novel coupled complex boundary method for solving inverse source problems.  Inverse Problems, page Article 055002, 2014.
  * [CK19] D. Colton and R. Kress.  Inverse Acoustic and Electromagnetic Scattering Theory.  Springer-Verlag, New York, 4th edition, 2019.
  * [Dam02] M. Dambrine.  On variations of the shape Hessian and sufficient conditions for the stability of critical shapes.  Rev. R. Acad. Cienc. Exactas Fis. Nat. Ser. A. Mat., 96:95–121, 2002.
  * [DL98] R. Dautray and J.-L. Lions.  Mathematical Analysis and Numerical Methods for Science and Technology, volume 2.  Springer, 1998.
  * [DMNV07] G. Doǧan, P. Morin, R.H. Nochetto, and M. Verani.  Discrete gradient flows for shape optimization and applications.  Comput. Methods Appl. Mech. Engrg., 196:3898–3914, 2007.
  * [DZ11] M. C. Delfour and J.-P. Zolésio.  Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization, volume 22 of Adv. Des. Control.  SIAM, Philadelphia, 2nd edition, 2011.
  * [EH05] K. Eppler and H. Harbrecht.  A regularized Newton method in electrical impedance tomography using shape hessian information.  Control Cybern., 34:203–225, 2005.
  * [Gal94] G. P. Galdi.  An introduction to the mathematical theory of the Navier-Stokes equations.  Springer Tracts in Natural Philosophy. Springer-Verlag, New York, 1994\.
  * [GCH17] R. Gong, X. Cheng, and W. Han.  A coupled complex boundary method for an inverse conductivity problem with one measurement.  Appl. Anal., 96(5):869–885, 2017.
  * [GGM01] S. Garreau, P. Guillaume, and M. Masmoudi.  The topological asymptotic for pde systems: the elasticity case.  SIAM J. Control Optim., 39:1756–1778, 2001.
  * [Hec12] F. Hecht.  New development in FreeFem++.  J. Numer. Math., 20:251–265, 2012.
  * [HIK+09] J. Haslinger, K. Ito, T. Kozubek, K. Kunisch, and G. H. Peichl.  On the shape derivative for problems of Bernoulli type.  Interfaces Free Bound., 11:317–330, 2009.
  * [HP18] A. Henrot and M. Pierre.  Shape Variation and Optimization: A Geometrical Analysis, volume 28 of Tracts in Mathematics.  European Mathematical Society, Zürich, 2018.
  * [IKP06] K. Ito, K. Kunisch, and G. H. Peichl.  Variational approach to shape derivative for a class of Bernoulli problem.  J. Math. Anal. Appl., 314(2):126–149, 2006.
  * [IKP08] K. Ito, K. Kunisch, and G. H. Peichl.  Variational approach to shape derivatives.  ESAIM Control Optim. Calc. Var., 14:517–539, 2008.
  * [KK11] H. Kasumba and K. Kunisch.  On shape sensitivity analysis of the cost functional without shape sensitivity of the state variable.  Control. Cybern., 40(4):989–1017, 2011.
  * [KK14] H. Kasumba and K. Kunisch.  On computation of the shape hessian of the cost functional without shape sensitivity of the state variable.  J. Optim. Theory Appl., 162:779–804, 2014.
  * [KO22] V. A. Kovtunenko and H. Ohtsuka.  Inverse problem of shape identification from boundary measurement for Stokes equations: shape differentiability of Lagrangian.  J. Inverse Ill-Posed Probl. 2022, 30(4):461–474, 2022.
  * [Neu97] J. W. Neuberger.  Sobolev Gradients and Differential Equations.  Springer-Verlag, Berlin, 1997.
  * [OCNN22] H. Ouaissa, A. Chakib, A. Nachaoui, and M. Nachaoui.  On numerical approaches for solving an inverse Cauchy Stokes problem.  Appl. Math. Optim., 85(Art. 3):37 pp., 2022.
  * [RA18] J. F. T. Rabago and H. Azegami.  Shape optimization approach to defect-shape identification with convective boundary condition via partial boundary measurement.  Japan J. Indust. Appl. Math., 31(1):131–176, 2018.
  * [RA20] J. F. T. Rabago and H. Azegami.  A second-order shape optimization algorithm for solving the exterior Bernoulli free boundary problem using a new boundary cost functional.  Comput. Optim. Appl., 77(1):251–305, 2020.
  * [Rab22] J. F. T. Rabago.  On the new coupled complex boundary method in shape optimization framework for solving stationary free boundary problems.  Math. Control Relat. Fields, 2022.
  * [RN24] J. F. T. Rabago and H. Notsu.  Numerical solution to a free boundary problem for the stokes equation using the coupled complex boundary method in shape optimization setting.  Appl. Math. Optim., 89:Art. 2, 2024.
  * [SN22a] J. S. H. Simon and H. Notsu.  A shape design problem for the Navier-Stokes flow with a convective boundary condition.  Comput. Appl. Math., 41:Art. 167, 2022.
  * [SN22b] J. S. H. Simon and H. Notsu.  A shape optimization problem constrained with the Stokes equations to address maximization of vortices.  Evol. Equ. Control Theory, 11:1873–1902, 2022.
  * [SZ92] J. Sokołowski and J.-P. Zolésio.  Introduction to Shape Optimization: Shape Sensitivity Analysis.  Springer Series in Computational Mathematics. Springer-Verlag, Berlin, Heidelberg, 1992.
  * [SZ99] J. Sokołowski and A. Żochowski.  On the topological derivative in shape optimization.  SIAM J. Control Optim., 37:1251–1272, 1999.
  * [Tem01] R. Temam.  Navier-Stokes Equations, Theory and Numerical Analysis.  Amer. Math. Soc., reprint of the 1984 edition edition, 2001.
  * [ZCG20] X. Zheng, X. Cheng, and R. Gong.  A coupled complex boundary method for parameter identification in elliptic problems.  Int. J. Comput. Math., 97(5):998–1015, 2020.

##  Appendix A Lemmata Proofs

In this appendix, we present the proofs of the key lemmas employed in the
previous section. Let us first note that the following regularities hold (see,
e.g., [IKP06, IKP08]):

$\left\\{\begin{aligned}
&\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}^{1,1}(\overline{D})^{d}),&&[t\mapsto
T_{t}^{-1}]\in\mathcal{C}(\mathcal{I},\mathcal{C}^{1}(\overline{D})^{d}),\\\
[t\mapsto
DT_{t}]&\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}^{0,1}(\overline{D})^{d\times
d}),&&[t\mapsto(DT_{t})^{-\top}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}(\overline{D})^{d\times
d}),\\\ [t\mapsto
I_{t}]&\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}(\overline{\Omega})),&&\
[t\mapsto
I_{t}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}^{0,1}(\overline{D})),\\\
[t\mapsto
A_{t}]&\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}(\overline{\Omega})^{d\times
d}),&&[t\mapsto
A_{t}]\in\mathcal{C}(\mathcal{I},\mathcal{C}(\overline{D})^{d\times
d}).\end{aligned}\right.$ (A1)

###  A.1 Proof of Lemma 3.4

The pair
$({\bm{\mathbf{u}}_{t}},{p_{t}})\in{\bm{\mathbf{\mathsf{H}}}_{\Gamma,\bm{\mathbf{0}}}^{1}(\Omega_{t})}^{d}\times\bm{\mathbf{\mathsf{L}}}^{2}(\Omega_{t})$
solves the system

$\left\\{\begin{aligned}
\int_{\Omega_{t}}{{\alpha}\nabla{\bm{\mathbf{u}}_{t}}:\nabla{\overline{\bm{\mathbf{\varphi}}}}}{\,{d}x_{t}}+i\int_{\Sigma_{t}}{{\bm{\mathbf{u}}_{t}}\cdot{\overline{\bm{\mathbf{\varphi}}}}}{\,{d}\sigma_{t}}&-\int_{\Omega_{t}}{p_{t}(\nabla\cdot\overline{\bm{\mathbf{\varphi}}})}{\,{d}x_{t}}\\\
&=\int_{\Sigma_{t}}{(\bm{\mathbf{g}}_{t}+i\bm{\mathbf{f}}_{t})\cdot\overline{\bm{\mathbf{\varphi}}}}{\,{d}\sigma_{t}},\quad\forall{\bm{\mathbf{\varphi}}}\in{\bm{\mathbf{\mathsf{H}}}_{\Gamma,\bm{\mathbf{0}}}^{1}(\Omega_{t})}^{d},\\\
-\int_{\Omega_{t}}{\bar{\lambda}(\nabla\cdot\bm{\mathbf{u}}_{t})}{\,{d}x_{t}}&=0,\quad\forall\lambda\in{\bm{\mathbf{\mathsf{L}}}^{2}(\Omega_{t})}.\end{aligned}\right.$

By change of variables (cf. [DZ11, pp. 482–484]), the fact that B${}_{t}=1$ on
$\Sigma$, for all $t>0$, $\bm{\mathbf{u}}^{t}={\bm{\mathbf{u}}_{t}}\circ
T_{t}=({\bm{\mathbf{u}}_{rt}}+i{\bm{\mathbf{u}}_{it}})\circ T_{t}$, and
together with the identities

$(\nabla\bm{\mathbf{\varphi}}_{t})\circ
T_{t}=DT_{t}^{-\top}\nabla\bm{\mathbf{\varphi}}^{t}\quad\text{and}\quad(\nabla\cdot\bm{\mathbf{\varphi}}_{t})\circ
T_{t}=(DT_{t})^{-1}:\nabla\bm{\mathbf{\varphi}}^{t},$

where
$\bm{\mathbf{\varphi}}_{t}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega_{t})^{d}$ and
$\bm{\mathbf{\varphi}}^{t}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)^{d}$, we get

$\left\\{\begin{aligned}
\int_{\Omega}{{\alpha}A_{t}\nabla{\bm{\mathbf{u}}^{t}}:\nabla{\overline{\bm{\mathbf{\varphi}}}}}{\,{d}x}+i\int_{\Sigma}{{\bm{\mathbf{u}}^{t}}\cdot{\overline{\bm{\mathbf{\varphi}}}}}{\,{d}\sigma}-\int_{\Omega}{I_{t}p^{t}(M_{t}^{\top}:\nabla\overline{\bm{\mathbf{\varphi}}})}{\,{d}x}&=\int_{\Sigma}{(\bm{\mathbf{g}}^{t}+i\bm{\mathbf{f}}^{t})\cdot\overline{\bm{\mathbf{\varphi}}}}{\,{d}\sigma},\quad\forall{\bm{\mathbf{\varphi}}}\in{V_{\Gamma}},\\\
-\int_{\Omega}{I_{t}\bar{\lambda}(M_{t}^{\top}:\nabla{\bm{\mathbf{u}}^{t}})}{\,{d}x}&=0,\quad\forall\lambda\in{Q}.\end{aligned}\right.$

Using the notations in (30), we recover (32).

The rest of the proof is again a routine step. The well-posedness is achieved
using standard arguments similar to the proof of Proposition 2.3, which
follows the same lines of argumentation for the real case. More precisely, the
argumentations make use of Lemma 3.3 which exploits the properties of $I_{t}$,
$A_{t}$, and $B_{t}$ issued in (18) and (A1). Concerning uniqueness of
solution, the proof (with a compatibility condition imposed) is also standard,
so we omit it.

###  A.2 Proof of Lemma 3.5

We first obtain a uniform bound for $\bm{\mathbf{u}}^{t}$ in $X$ for all
$t\in\mathcal{I}$. To do this, we take
$(\bm{\mathbf{\varphi}},\lambda)=(\bm{\mathbf{u}}^{t},p^{t})\in
V_{\Gamma}\times Q$ in (32), utilize the properties of $A_{t}$ and $B_{t}$
given in (18) and (A1). as well as the coercivity of $a^{t}$ on $X\times X$
(Lemma 3.3). From (32), we have
$a^{t}({\bm{\mathbf{u}}^{t}},{\bm{\mathbf{u}}^{t}})=F^{t}({\bm{\mathbf{u}}^{t}})$.
Using the relaxed $V_{\Gamma}$-ellipticity property of $a^{t}$ (cf. Lemma
3.3), we get the following sequence of inequalities

${|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{2}\lesssim\Re\left\\{\int_{\Omega}{{\alpha}A_{t}|{\nabla{\bm{\mathbf{u}}^{t}}}|^{2}}{\,{d}x}+i\int_{\Sigma}{|{\bm{\mathbf{u}}^{t}}|^{2}}{\,{d}\sigma}\right\\}\\\
\lesssim\left\|{(\bm{\mathbf{f}},\bm{\mathbf{g}})}\right\|_{1/2,\Sigma}{|\kern-1.07639pt|\kern-1.07639pt|\overline{\bm{\mathbf{u}}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Sigma}.$

Applying the trace theorem, and the fact that
$\bm{\mathbf{g}}^{t}=\bm{\mathbf{g}}$ and
$\bm{\mathbf{f}}^{t}=\bm{\mathbf{f}}$ on $\Sigma$, we get the estimate

${|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}\lesssim\left\|{(\bm{\mathbf{f}},\bm{\mathbf{g}})}\right\|_{1/2,\Sigma},$
(A2)

which shows that $\bm{\mathbf{u}}^{t}$ is bounded in $X$ for all
$t\in\mathcal{I}$.

For the boundedness of $p^{t}$ in $Q$ for all $t\in\mathcal{I}$, we employ
(31), which is equivalent to the following inequality condition

$\sup_{\begin{subarray}{c}\bm{\mathbf{\varphi}}\in{V_{\Gamma}}\\\
\bm{\mathbf{\varphi}}\neq\bm{\mathbf{0}}\end{subarray}}\frac{b^{t}(\bm{\mathbf{\varphi}},{\lambda})}{{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}}\geqslant\beta_{1}{|\kern-1.07639pt|\kern-1.07639pt|{\lambda}|\kern-1.07639pt|\kern-1.07639pt|}_{Q},\quad\forall{\lambda}\in
Q.$

We let $\lambda=p^{t}\in Q$. Then, from (3.4) and (18) we obtain the following
sequence of inequalities

$\displaystyle\beta_{1}{|\kern-1.07639pt|\kern-1.07639pt|{p^{t}}|\kern-1.07639pt|\kern-1.07639pt|}_{Q}$
$\displaystyle\leqslant\sup_{\begin{subarray}{c}\bm{\mathbf{\varphi}}\in{V_{\Gamma}}\\\
\bm{\mathbf{\varphi}}\neq\bm{\mathbf{0}}\end{subarray}}\frac{b^{t}(\bm{\mathbf{\varphi}},{p^{t}})}{{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}}=\sup_{\begin{subarray}{c}\bm{\mathbf{\varphi}}\in{V_{\Gamma}}\\\
\bm{\mathbf{\varphi}}\neq\bm{\mathbf{0}}\end{subarray}}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{-1}\Bigg{\\{}F^{t}(\bm{\mathbf{\varphi}})-a^{t}({\bm{\mathbf{u}}^{t}},{\bm{\mathbf{\varphi}}})\Bigg{\\}}$
$\displaystyle\leqslant\sup_{\begin{subarray}{c}\bm{\mathbf{\varphi}}\in{V_{\Gamma}}\\\
\bm{\mathbf{\varphi}}\neq\bm{\mathbf{0}}\end{subarray}}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{-1}\Bigg{\\{}\int_{\Sigma}{(\bm{\mathbf{g}}+i\bm{\mathbf{f}})\cdot\overline{\bm{\mathbf{\varphi}}}}{\,{d}\sigma}-\int_{\Omega}{{\alpha}A_{t}\nabla{\bm{\mathbf{u}}^{t}}:\nabla{\overline{\bm{\mathbf{\varphi}}}}}{\,{d}x}-\int_{\Sigma}{{\bm{\mathbf{u}}^{t}}\cdot{\overline{\bm{\mathbf{\varphi}}}}}{\,{d}\sigma}\Bigg{\\}}$
$\displaystyle\leqslant
c\sup_{\begin{subarray}{c}\bm{\mathbf{\varphi}}\in{V_{\Gamma}}\\\
\bm{\mathbf{\varphi}}\neq\bm{\mathbf{0}}\end{subarray}}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{-1}\Bigg{\\{}\left\|{(\bm{\mathbf{f}},\bm{\mathbf{g}})}\right\|_{1/2,\Sigma}{|\kern-1.07639pt|\kern-1.07639pt|\overline{\bm{\mathbf{\varphi}}}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega}+{|\kern-1.07639pt|\kern-1.07639pt|\nabla{\bm{\mathbf{u}}^{t}}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega}{|\kern-1.07639pt|\kern-1.07639pt|\nabla{\overline{\bm{\mathbf{\varphi}}}}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega}+{|\kern-1.07639pt|\kern-1.07639pt|{\bm{\mathbf{u}}^{t}}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Sigma}{|\kern-1.07639pt|\kern-1.07639pt|{\overline{\bm{\mathbf{\varphi}}}}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Sigma}\Bigg{\\}}$
$\displaystyle\leqslant
c\left(\left\|{(\bm{\mathbf{f}},\bm{\mathbf{g}})}\right\|_{1/2,\Sigma}+{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}\right).$

Therefore, using the estimate for
${|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}$
given in (A2), we also have the estimate
${|\kern-1.07639pt|\kern-1.07639pt|{p^{t}}|\kern-1.07639pt|\kern-1.07639pt|}_{Q}\lesssim\left\|{(\bm{\mathbf{f}},\bm{\mathbf{g}})}\right\|_{1/2,\Sigma}$,
completing the proof of the lemma.

###  A.3 Proof of Lemma 3.6

The proof of the lemma mainly consists of two key steps: the first step is to
show that $\lim_{t\to 0}\bm{\mathbf{u}}^{t}=\bm{\mathbf{u}}$ in $X$ and
$\lim_{t\to 0}p^{t}=p$ in $Q$, and the second one is to prove that $\lim_{t\to
0^{+}}\frac{1}{\sqrt{t}}({|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}+{|\kern-1.07639pt|\kern-1.07639pt|p^{t}-p|\kern-1.07639pt|\kern-1.07639pt|}_{Q})=0$.

Step 1. We consider the difference between the variational equations (32) and
(12), and define $\bm{\mathbf{y}}^{t}:=\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}$
and $r^{t}:=p^{t}-p$. By making $\varepsilon>0$ smaller if necessary, it can
be shown that the following expansions hold for all $t\in\mathcal{I}$,

$\left\\{\begin{aligned}
I_{t}&=1+t\operatorname{div}\bm{\mathbf{\theta}}+t^{2}\tilde{\rho}(t,\bm{\mathbf{\theta}}),\quad\text{where}\
\tilde{\rho}\in\mathcal{C}(\mathbb{R},\mathcal{C}^{0,1}(D)),\\\
M_{t}^{\top}&=(DT_{t})^{-1}=(id+tD\bm{\mathbf{\theta}})^{-1}=id-
tD\bm{\mathbf{\theta}}+O{(t^{2})}id,\quad\text{where}\
0\leqslant{{t}}\leqslant\varepsilon<|\lambda_{max}|^{-1},\end{aligned}\right.$

where $\lambda_{max}$ is the maximum singular value of
$D\bm{\mathbf{\theta}}$. Moreover, we denote

$\rho(t):=t^{2}\tilde{\rho}(t,\bm{\mathbf{\theta}}),\quad
R(t):=O{(t^{2})}id,\quad\rho_{1}(t):=t\tilde{\rho}(t,\bm{\mathbf{\theta}}),\quad\text{and}\quad
R_{1}(t):=O{(t)}id.$

Let us consider the variational equation

$b^{t}(\bm{\mathbf{u}}^{t},\lambda)-b(\bm{\mathbf{u}},\lambda)=-\int_{\Omega}{I_{t}\bar{\lambda}(M_{t}^{\top}:\nabla{\bm{\mathbf{u}}^{t}})}{\,{d}x}-\left(-\int_{\Omega}{\overline{\lambda}\nabla\cdot\bm{\mathbf{u}}}{\,{d}x}\right)=0,\quad\forall\lambda\in
Q.$

By applying the expansion for $M_{t}^{\top}$ and taking
$\lambda=r^{t}=p^{t}-p\in Q$ in the above equation, it can be shown that

$\displaystyle
b(\bm{\mathbf{y}}^{t},r^{t})=\int_{\Omega}{r^{t}\nabla\cdot\overline{\bm{\mathbf{y}}}^{t}}{\,{d}x}$
$\displaystyle=-\int_{\Omega}{(t\nabla\cdot\bm{\mathbf{\theta}}+\rho(t))r^{t}(M_{t}^{\top}:\nabla\overline{\bm{\mathbf{u}}}^{t})}{\,{d}x}$
(A3)
$\displaystyle\qquad-\int_{\Omega}{r^{t}[(-tD\bm{\mathbf{\theta}}+R(t)):\nabla\overline{\bm{\mathbf{u}}}^{t}]}{\,{d}x}.$

On the other hand, we also have the equation

$a(\bm{\mathbf{y}}^{t},\bm{\mathbf{\varphi}})+b(\bm{\mathbf{\varphi}},r^{t})=\Phi^{t}(\bm{\mathbf{\varphi}}),\qquad\forall{\bm{\mathbf{\varphi}}}\in
V_{\Gamma},$ (A4)

where $\Phi^{t}(\bm{\mathbf{\varphi}})$ is given by (36), and
$a:V_{\Gamma}\times V_{\Gamma}\to\mathbb{R}$ and $b:V_{\Gamma}\times
Q\to\mathbb{R}$ are given in (11). In above, we have used the fact that
$B_{t}=1$ on $\Sigma$, for all $t>0$, and that
$\bm{\mathbf{\theta}}=\bm{\mathbf{0}}$ on $\Sigma$. Choosing
$\bm{\mathbf{\varphi}}=\bm{\mathbf{y}}^{t}\in V_{\Gamma}$ in (A4) and
utilizing identity (A3), it follows that

$\displaystyle
c_{a}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{2}$
$\displaystyle\
\leqslant\left|\Phi^{t}(\bm{\mathbf{y}}^{t})-b(\bm{\mathbf{y}}^{t},p^{t}-p)\right|$
$\displaystyle\
\leqslant{{\alpha}}\left|{\mathfrak{a}_{t}}\right|_{\infty}{|\kern-1.07639pt|\kern-1.07639pt|\nabla{\bm{\mathbf{u}}^{t}}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega}{|\kern-1.07639pt|\kern-1.07639pt|\nabla{\overline{\bm{\mathbf{y}}}^{t}}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega}+\left(\left|{\mathfrak{i}}_{t}\right|_{\infty}|{M_{t}^{\top}}|_{\infty}+\left|{\mathfrak{m}_{t}}\right|_{\infty}\right){|\kern-1.07639pt|\kern-1.07639pt|p^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega}{|\kern-1.07639pt|\kern-1.07639pt|\nabla{\overline{\bm{\mathbf{y}}}^{t}}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega}$
$\displaystyle\
\quad+{{t}}\left[\left(|{\nabla\cdot\bm{\mathbf{\theta}}}|_{\infty}+|{\rho_{1}(t)}|_{\infty}\right)|{M_{t}^{\top}}|_{\infty}+|{D\bm{\mathbf{\theta}}}|_{\infty}+|{R_{1}{(t)}}|_{\infty}\right]{|\kern-1.07639pt|\kern-1.07639pt|r^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega}{|\kern-1.07639pt|\kern-1.07639pt|\nabla\overline{\bm{\mathbf{u}}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{0,\Omega}.$

Therefore, we have the estimate

$c_{a}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{2}\leqslant{m}_{t}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}+t{\Xi}^{t}{|\kern-1.07639pt|\kern-1.07639pt|r^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{Q},$
(A5)

where

$\left\\{\begin{aligned}
{m}_{t}&:={{\alpha}}\left|{\mathfrak{a}_{t}}\right|_{\infty}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}+\left(\left|{I}_{t}\right|_{\infty}|{M_{t}^{\top}}|_{\infty}+\left|{\mathfrak{m}_{t}}\right|_{\infty}\right){|\kern-1.07639pt|\kern-1.07639pt|p^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{Q},\\\
{\Xi}^{t}&:=\left[\left(|{\nabla\cdot\bm{\mathbf{\theta}}}|_{\infty}+|{\rho_{1}(t)}|_{\infty}\right)|{M_{t}^{\top}}|_{\infty}+|{D\bm{\mathbf{\theta}}}|_{\infty}+|{R_{1}{(t)}}|_{\infty}\right]{|\kern-1.07639pt|\kern-1.07639pt|{\bm{\mathbf{u}}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}.\end{aligned}\right.$
(A6)

Note that the quantities ${m}_{t}$ and ${\Xi}^{t}$ are finite for all
$t\in\mathcal{I}$ by Lemma 3.5.

Now, to proceed further, we take $\lambda=r^{t}=p^{t}-p\in Q$ in the inf-sup
condition (13) and consider equation (A4) to obtain

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|{r^{t}}|\kern-1.07639pt|\kern-1.07639pt|}_{Q}\leqslant\beta_{0}^{-1}\sup_{\begin{subarray}{c}\bm{\mathbf{\varphi}}\in{V_{\Gamma}}\\\
\bm{\mathbf{\varphi}}\neq\bm{\mathbf{0}}\end{subarray}}\frac{b(\bm{\mathbf{\varphi}},{r^{t}})}{{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}}$
$\displaystyle=\beta_{0}^{-1}\sup_{\begin{subarray}{c}\bm{\mathbf{\varphi}}\in{V_{\Gamma}}\\\
\bm{\mathbf{\varphi}}\neq\bm{\mathbf{0}}\end{subarray}}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{\varphi}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{-1}\Big{\\{}\Phi^{t}(\bm{\mathbf{\varphi}})-a(\bm{\mathbf{y}}^{t},\bm{\mathbf{\varphi}})\Big{\\}}$
(A7)
$\displaystyle\leqslant\beta_{0}^{-1}\left({m}_{t}+\max\\{{{\mu}},1\\}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}\right).$

Going back to (A5) and utilizing the above estimate, we get

$\displaystyle
c_{a}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{2}$
$\displaystyle\leqslant
t{\Xi}^{t}\beta_{0}^{-1}{m}_{t}+\left({m}_{t}+t{\Xi}^{t}\beta_{0}^{-1}\max\\{{{\mu}},1\\}\right){|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}.$
(A8)

We apply Peter-Paul inequality to the second summand in above inequaltiy to
obtain the following estimate

$\displaystyle\left({m}_{t}+t{\Xi}^{t}\beta_{0}^{-1}\max\\{{{\mu}},1\\}\right){|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}$
$\displaystyle\leqslant\frac{\left({m}_{t}+t{\Xi}^{t}\beta_{0}^{-1}\max\\{{{\mu}},1\\}\right)^{2}}{2\varepsilon_{1}}+\frac{\varepsilon_{1}}{2}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{2}.$

for some constant $\varepsilon_{1}>0$. We choose and fixed $\varepsilon_{1}$
such that $\bar{c}:=c(c_{a},\varepsilon_{1}):=2c_{a}-\varepsilon_{1}>0$ so
that, from our first estimate (A8), we get

${|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}\leqslant\bar{c}^{-\frac{1}{2}}\left(t{\Xi}^{t}\beta_{0}^{-1}{m}_{t}+\frac{\left({m}_{t}+t{\Xi}^{t}\beta_{0}^{-1}\max\\{{{\mu}},1\\}\right)^{2}}{2\varepsilon_{1}}\right)^{1/2}.$
(A9)

Because $I_{t}|M_{t}\bm{\mathbf{n}}|^{-1}=(1)(|id\bm{\mathbf{n}}|^{-1})=1$ at
$t=0$, then, in view of (A6) together with Lemma 3.5, (21), and (22), we see
that ${m}_{t}\to 0$ as $t\to 0$. Moreover, since ${\Xi}^{t}$ is uniformly
bounded for all $t\in\mathcal{I}$ by Lemma 3.5, we deduce – by Lebesgue’s
dominated convergence theorem – that

$\lim_{t\to
0}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}=0\qquad\Longleftrightarrow\qquad\lim_{t\to
0}\bm{\mathbf{u}}^{t}=\bm{\mathbf{u}}\quad\text{in $X$}.$ (A10)

Similarly, based from the above discussion, we know that the quantities on the
right side of the inequality (A7) vanish as $t\to 0$. Therefore, we also have
the limit

$\lim_{t\to
0}{|\kern-1.07639pt|\kern-1.07639pt|r^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{Q}=0\qquad\Longleftrightarrow\qquad\lim_{t\to
0}p^{t}=p\quad\text{in $Q$}.$ (A11)

Step 2. For the second part of the proof, we start by noting that, for
sufficiently small $t>0$, $\frac{1}{t}\bm{\mathbf{y}}^{t}\in V_{\Gamma}$ and
$\frac{1}{t}r^{t}\in Q$. In addition, we recall that the derivatives of
$I_{t}$ and $M_{t}$ with respect to $t$ exists in $L^{\infty}(\Omega)$ and
$L^{\infty}(\Omega)^{d\times d}$, respectively. Now, to finish the proof, we
go back to inequality (A9) to obtain, after dividing by $t>0$, the following
estimate

$\displaystyle\frac{1}{t}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}^{2}\leqslant\bar{c}^{-1}\left({\Xi}^{t}\beta_{0}^{-1}{m}_{t}+\frac{1}{2\varepsilon_{1}}\left(\frac{1}{\sqrt{t}}{m}_{t}+\sqrt{t}{\Xi}^{t}\beta_{0}^{-1}\max\\{{{\mu}},1\\}\right)^{2}\right).$

Observe that we have

$\frac{1}{\sqrt{t}}{m}_{t}=\sqrt{t}{{\mu}}\left|\frac{\mathfrak{a}_{t}}{t}\right|_{\infty}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}\
+\sqrt{t}\left(\left|\frac{{\mathfrak{i}}_{t}}{t}\right|_{\infty}|{M_{t}^{\top}}|_{\infty}+\left|\frac{\mathfrak{m}_{t}}{t}\right|_{\infty}\right){|\kern-1.07639pt|\kern-1.07639pt|p^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{Q}.$

Thus, using (A1), (19), (22), Lemma 3.5, (A10), and (A11), we deduce that the
following limit holds

$\lim_{t\to
0^{+}}\frac{1}{\sqrt{t}}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{u}}^{t}-\bm{\mathbf{u}}|\kern-1.07639pt|\kern-1.07639pt|}_{X}=0.$
(A12)

Similarly, from estimate (A7), we know that (after dividing by $\sqrt{t}>0$)
the following inequality holds

$\displaystyle\frac{1}{\sqrt{t}}{|\kern-1.07639pt|\kern-1.07639pt|{r^{t}}|\kern-1.07639pt|\kern-1.07639pt|}_{Q}\leqslant\beta_{0}^{-1}\left[\frac{1}{\sqrt{t}}{m}_{t}+\max\\{{{\mu}},1\\}\left(\frac{1}{\sqrt{t}}{|\kern-1.07639pt|\kern-1.07639pt|\bm{\mathbf{y}}^{t}|\kern-1.07639pt|\kern-1.07639pt|}_{X}\right)\right].$

Again, applying (A1), (19), (22), and Lemma 3.5, but now combined with (A12),
we infer that

$\lim_{t\to
0^{+}}\frac{1}{\sqrt{t}}{|\kern-1.07639pt|\kern-1.07639pt|p^{t}-p|\kern-1.07639pt|\kern-1.07639pt|}_{Q}=0,$
(A13)

proving the lemma.

J. .F. .T. Rabago, Faculty of Mathematics and Physics, Institute of Science
and Engineering, Kanazawa University, Kakumamachi, Kanazawa 920-1192, Japan

E-mail address, J. .F. .T. Rabago<EMAIL_ADDRESS>
L. Afraites, EMI, FST, Université Sultan Moulay Slimane, Béni-Mellal 23000,
Morocco

E-mail address, L. Afraites<EMAIL_ADDRESS>
H. Notsu, Faculty of Mathematics and Physics, Institute of Science and
Engineering, Kanazawa University, Kakumamachi, Kanazawa 920-1192, Japan

E-mail address, H. Notsu<EMAIL_ADDRESS>