# On the new coupled complex boundary method in shape optimization framework
for solving stationary free boundary problems111Some of the results in this
paper were presented at the SAB-ATAN 2021: International Conference on
Mathematical Sciences and Applications held virtually on November 24-26, 2021
and hosted by the University of the Philippines Baguio.

Julius Fergy T. Rabago

(Faculty of Mathematics and Physics
Institute of Science and Engineering
Kanazawa University, Kanazawa 920-1192, Japan
<EMAIL_ADDRESS><EMAIL_ADDRESS>
)

###### Abstract

We expose here a novel application of the so-called coupled complex boundary
method – first put forward by Cheng et al. (2014) to deal with inverse source
problems – in the framework of shape optimization for solving the exterior
Bernoulli problem, a prototypical model of stationary free boundary problems.
The idea of the method is to transform the overdetermined problem to a complex
boundary value problem with a complex Robin boundary condition coupling the
Dirichlet and Neumann boundary conditions on the free boundary. Then, we
optimize the cost function constructed by the imaginary part of the solution
in the whole domain in order to identify the free boundary. We also prove the
existence of the shape derivative of the complex state with respect to the
domain. Afterwards, we compute the shape gradient of the cost functional, and
characterize its shape Hessian at the optimal domain under a strong, and then
a mild regularity assumption on the domain. We then prove the ill-posedness of
the proposed shape problem by showing that the latter expression is compact.
Also, we devise an iterative algorithm based on a Sobolev gradient scheme via
finite element method to solve the minimization problem. Finally, we
illustrate the applicability of the method through several numerical examples,
both in two and three spatial dimensions.

Keywords coupled complex boundary method, stationary free boundary problem,
Bernoulli problem, shape calculus, shape derivatives, adjoint method

##  1 Introduction

Let $A$ and $B$ be two bounded and simply connected domains in
$\mathbb{R}^{d}$, $d\in\\{2,3\\}$, with respective boundaries
$\Gamma:=\partial A$ and $\Sigma:=\partial B$, such that
$B\supset\overline{A}$. Denote the annular domain
$\Omega=B\setminus\overline{A}$ (which is assumed, unless stated otherwise, to
be non-empty bounded open Lipschitz subset of $\mathbb{R}^{d}$ throughout the
paper), then an exterior free boundary problem may be given as follows: for
given functions $f$, $g$, $h$, one tries to find $\Omega$ with the associated
function $u:=u(\Omega)$ such that the overdetermined problem

$-\Delta u=f\quad\text{in}\ \Omega,\qquad u=g\quad\text{on}\ \Gamma,\qquad
u=0\quad\text{and}\quad\partial_{\bm{\mathbf{n}}}{u}=h\quad\text{on}\ \Sigma,$
(1)

is satisfied, where $\partial_{\bm{\mathbf{n}}}{u}:=\nabla
u\cdot\bm{\mathbf{n}}$ is the outward normal derivative of $u$.

In this work, we are primarily interested in solving the free boundary problem
(FBP) (1) through the novel application of the so-called coupled complex
boundary method or CCBM in solving stationary FBPs through the context of
shape optimization. For simplicity of discussion, we will consider the
prototypical case of (1) in two spatial dimensions popularly known as the
exterior Bernoulli problem wherein $f\equiv 0$, $g\equiv 1$, and $h=\lambda$,
where $\lambda<0$ is a fixed constant in (1). That is, we consider the problem
of finding a pair $(\Omega,u):=(\Omega,u(\Omega))$ that solves the system

$-\Delta u=0\ \text{in}\ \Omega,\qquad u=1\ \text{on}\ \Gamma,\qquad
u=0\quad\text{and}\quad\partial_{\bm{\mathbf{n}}}{u}=\lambda\ \text{on}\
\Sigma.$ (2)

The system of partial differential equations (PDEs) admits a classical
solution for simply connected bounded domain $\Omega$ for any $\lambda<0$, and
uniqueness can be guaranteed for bounded convex domains $A$ [29]. Moreover, in
the said case, it was shown in [41, Thm 1.1] that the free boundary is
$\mathcal{C}^{2,\alpha}$ regular. For additional qualitative properties of
solutions, as well as their numerical treatment, we refer the readers, for
instance, to [29]. The problem in consideration is also known in the
literature as the Alt-Caffarelli problem [4]. It originates from the
description of free surfaces for ideal fluids [30], but copious industrial
applications leading to similar formulations to (2) arises in many other
related contexts, see [9, 28, 29]. We mention that, in this paper, we do not
tackle the question of existence of optimal shape solutions for the proposed
shape optimization problem. Instead, we will tacitly assume the existence of
optimal domains which is sufficiently regular to carry out a second-order
shape calculus. Even so, we mention that existence proofs and tools developed
in this direction are issued in [6, 36, 37, 38].

The methods of shape optimization is a well-established tool to solve free
boundary problems, and in the case of (2), the method can be applied in
several ways. The usual strategy is to choose one of the boundary conditions
on the free boundary to obtain a well-posed state equation and then track the
remaining boundary data in $L^{2}(\Sigma)$ (see [24, 25, 35, 37, 42, 54, 55])
or utilize the Dirichlet energy functional as a shape functional (see [22,
34]). Another option is to consider an energy-gap type cost function which
consists of two auxiliary states; one that is a solution of a pure Dirichlet
problem and one that satisfies a mixed Dirichlet-Neumann problem (see [1, 5,
26]). The latter formulation is precisely given as the minimization problem

$J_{KV}(\Omega):=\frac{1}{2}\int_{\Omega}{\left|\nabla\left(u_{\text{N}}-u_{\text{D}}\right)\right|^{2}}{\,{\operatorname{\mathnormal{d}}}x}\
\to\ \inf,$ (3)

where the state functions $u_{\text{N}}$ and $u_{\text{D}}$ respectively
satisfy the well-posed systems

$\displaystyle-\Delta u_{\text{N}}=0\ \text{in}\ \Omega,$ $\displaystyle
u_{\text{N}}=1\ \text{on}\ \Gamma,\quad\qquad$
$\displaystyle\partial_{\bm{\mathbf{n}}}{u_{\text{N}}}=\lambda\ \text{on}\
\Sigma;$ (4) $\displaystyle-\Delta u_{\text{D}}=0\ \text{in}\ \Omega,$
$\displaystyle u_{\text{D}}=1\ \text{on}\ \Gamma,\quad\qquad$ $\displaystyle
u_{\text{D}}=0\ \text{on}\ \Sigma.$ (5)

The equivalence between the above shape optimization formulation and the
exterior Bernoulli problem (2) issues from the following statement. If
$(u,\Omega)$ is a solution of (2), then $u_{\text{N}}=u_{\text{D}}=u$;
therefore, $J_{KV}(\Omega)\equiv 0$. Conversely, if $J_{KV}(\Omega)=0$, then
$u=u_{\text{N}}=u_{\text{D}}$ is a solution of problem (2) because
$u_{\text{D}}-u_{\text{N}}\in H^{1}_{\Gamma,0}:=\\{\varphi\in
H^{1}(\Omega)\mid\varphi=0\ \text{on}\ \Gamma\\}$,
$J_{KV}(\Omega)=\frac{1}{2}\big{|}u_{\text{D}}-u_{\text{N}}\big{|}_{H^{1}(\Omega)}$
is a norm on $H^{1}_{\Gamma,0}(\Omega)$. Formulation (3) is better known in
the literature as the Kohn-Vogelius approach (see, e.g. [5, 26, KV87]).

Recently, some modifications of the above mentioned existing methods were
offered and examined in [51, 52], including the so-called Dirichlet-data-gap
cost functional approach – firstly proposed in [53]. These approaches make use
of a Robin problem as one of the state problem. Additionally, we mention that
there is another numerical technique – closely related to shape optimization
method and was originally conceived to solve moving boundary problems via
finite element method – called the comoving mesh method or CMM developed in
[60] which can be applied to solve the free boundary problem (2); see [60,
Section 3] for details.

In this paper, we want to offer yet another shape optimization approach to
solve (2). The idea is somewhat similar to [51, 62], but applies the concept
of complex PDEs. More exactly, we showcase here a new application of the so-
called coupled complex boundary method for solving stationary free boundary
problems. The idea of the method is simple: we couple the Dirichlet and
Neumann data in a Robin boundary condition in such a way that the Dirichlet
data and the Neumann data are the respective real and imaginary parts of the
Robin boundary condition. As a result, the conditions that have to be
satisfied on the free boundary are transformed into one condition that needs
to be satisfied on the domain. With the new method, as in problem (3), the
objective function can then be defined as a volume integral, see (10).

The CCBM was first put forward by Cheng et al. in [8] for solving an inverse
source problem (see also [7]) and was then used to solve the Cauchy problem
stably in [7]. It is later on applied to an inverse conductivity problem with
one measurement in [32] and also to parameter identification in elliptic
problems in [64]. Much more recently, CCBM was also applied in solving inverse
obstacle problems by Afraites in [2]. To the best of our knowledge – as the
method has not been applied yet to solving FBPs in previous investigations –
this is the first time that CCBM will be explored to deal with stationary
FBPs, particularly as a numerical resolution to the exterior Bernoulli
problem. This, in turn, provides new directions in treating related free
boundary/surface problems in the context of shape optimization.

The remainder of the paper is as follows. In Section 2, we describe how CCBM
actually applies to solving problem (2) in the framework of shape
optimization. We also give further motivations about the present study –
providing merits to considering CCBM as a way to solve FBPs. Section 3 is
devoted to proving that the map $t\mapsto u^{t}$ is $\mathcal{C}^{1}$ in a
neighborhood of zero, and to the characterization of its derivative. Analogous
results for the map $t\mapsto J(\Omega_{t})$ are also exhibited therein,
including especially the first-order shape derivative of the cost. We also
derive its second-order shape derivative in subsection 4.2, looking
particularly on its structure at the solution of (2). With the latter
expression, we examine the instability of the shape optimization problem in
subsection 4.3 by proving the compactness of the shape Hessian of $J$ at a
critical shape. In Section 5, we devise a numerical algorithm (subsection 5.1)
based on Sobolev gradient method to solve the shape optimization problem in
consideration. This is followed by an intermediate subsection (subsection
5.2), where we state and prove a small result concerning a stationary point of
an evolving domain that evolves according to a given pseudo flow field related
to Sobolev gradient method used in this investigation. Then, we illustrate the
feasibility of the proposed method through several numerical examples both in
two and three spatial dimensions (subsection 5.3). Lastly, in Appendix A.1, a
rigorous proof of the existence of the material derivative of the state, as
well as the characterization of its corresponding shape derivative, are
provided.

##  2 CCBM in shape optimization setting

We present here the proposed CCBM formulation of (2) and give the motivation
of the method for solving the free boundary problem.

###  2.1 Formulation and Notations

In what follows we present the CCBM formulation of (2) and its shape
optimization reformulation by introducing the least-squares fitting for the
imaginary part of the complex PDE solution. The discussion will be issued in
the case of two dimensions, but the results easily extend in three dimensions.
The main point of departure of the method is to recast (2) into the complex
PDE system

$-\Delta u=0\ \text{in}\ \Omega,\qquad\quad u=1\ \text{on}\
\Gamma,\qquad\quad\partial_{\bm{\mathbf{n}}}{u}+iu=\lambda\ \text{on}\
\Sigma,$ (6)

where $i=\sqrt{-1}$ stands for the unit imaginary number. Letting
$u=u_{1}+iu_{2}$ denote the solution of (6), it can be verified that the real-
valued functions $u_{1}$ and $u_{2}$ respectively satisfy the real PDE
systems:

$\displaystyle-\Delta u_{1}=0\ \text{in}\ \Omega,\qquad\qquad u_{1}=1\
\text{on}\ \Gamma,$
$\displaystyle\partial_{\bm{\mathbf{n}}}{u_{1}}-u_{2}=\lambda\ \text{on}\
\Sigma;$ (7) $\displaystyle-\Delta u_{2}=0\ \text{in}\ \Omega,\qquad\qquad
u_{2}=0\ \text{on}\ \Gamma,$ $\displaystyle
u_{1}+\partial_{\bm{\mathbf{n}}}{u_{2}}=0\ \text{on}\ \Sigma.$ (8)

As mentioned earlier, a Robin-type boundary condition has already been used in
[51, 52, 53, 62]. The formulation issued here, however, is not covered by any
of the aforementioned study. This is because the Robin coefficient appearing
in (6) is a complex number while the one used in [51, 52, 53] is a positive
real number. Moreover, the Robin coefficient in [62] is assumed to be a non-
negative real number.

######  Remark 1.

Let us note that if $u_{2}=0$ in $\Omega$, then we have
$u_{2}=\partial_{\bm{\mathbf{n}}}{u_{2}}=0$ on $\Sigma$ and $u_{1}=0$ on
$\Sigma$. In view of (7) and (8), we see that the pair $(\Omega,u_{1})$ solves
the original free boundary problem (2). Conversely, if $(\Omega,u)$ is the
solution to the free boundary problem (2), then clearly $u_{1}$ and $u_{2}$
satisfy (7) and (8).

We infer from the previous remark that the original free boundary problem (2)
can be recast into an equivalent shape problem given as follows.

######  Problem 1.

Given a fixed interior boundary $\Gamma$ and a real number $\lambda<0$, find
an annular domain $\Omega$, with the exterior boundary
$\Sigma:=\partial\Omega\setminus\Gamma$, and a function $u:=u(\Omega)$ such
that $u_{2}=0$ in $\Omega$ and $u=u_{1}+iu_{2}$ solves the PDE system (6).

Notations. The notations for the function spaces used in the paper are as
follows. We denote by $W^{m,p}(\Omega)$ the standard real Sobolev space with
the norm $\|\cdot\|_{W^{m,p}(\Omega)}$. Let $W^{0,p}(\Omega)=L^{p}(\Omega)$,
and particularly, $H^{m}(\Omega)$ represents $W^{m,2}(\Omega)$ with the
corresponding inner product $(\cdot,\cdot)_{m,\Omega}$ and norm
$\|\cdot\|_{H^{m}(\Omega)}$. We let $\bm{\mathbf{\mathsf{H}}}^{m}(\Omega)$ be
the complex version of $H^{m}(\Omega)$ with the inner product
$(\\!(\cdot,\cdot)\\!)_{m,\Omega}$ and norm
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\,\cdot\,\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{H}}}^{m}(\Omega)}$
defined as follows: for all $u,v\in\bm{\mathbf{\mathsf{H}}}^{m}(\Omega)$,
$(\\!(u,v)\\!)_{m,\Omega}=(u,\overline{v})_{m,\Omega}$ and
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{H}}}^{m}(\Omega)}^{2}=(\\!(v,v)\\!)_{m,\Omega}$.
Similarly, we denote $V(\Omega):=H^{1}_{\Gamma,0}(\Omega)$,
$\bm{\mathbf{\mathsf{V}}}(\Omega):=\bm{\mathbf{\mathsf{H}}}^{1}_{\Gamma,0}(\Omega)$,
$Q=L^{2}(\Omega)$,
$\bm{\mathbf{\mathsf{Q}}}=\bm{\mathbf{\mathsf{L}}}^{2}(\Omega)$,
$S=L^{2}(\Sigma)$,
$\bm{\mathbf{\mathsf{S}}}=\bm{\mathbf{\mathsf{L}}}^{2}(\Sigma)$. Lastly, we
define
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{V}}}(\Omega)}^{2}:=\int_{\Omega}{\nabla{v}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+\int_{\Sigma}{{v}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$.

Before we go further, we remark that the complex PDE system (6) is well-posed.
Indeed, with the sesquilinear form $\bm{\mathbf{\mathsf{a}}}$ defined on
$\bm{\mathbf{\mathsf{V}}}(\Omega)\times\bm{\mathbf{\mathsf{V}}}(\Omega)$ by

$\bm{\mathbf{\mathsf{a}}}({u},{v})=\int_{\Omega}{\nabla{u}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{{u}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma},\quad\forall{u},{v}\in\bm{\mathbf{\mathsf{V}}}(\Omega),$

and the linear form
$l(v)=\lambda\int_{\Sigma}{{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$,
we may state (6) in variational form as follows:

$\text{find ${u}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)$, \ $u=1$ on $\Gamma$,
\ \ such that\ \ $\bm{\mathbf{\mathsf{a}}}({u},{v})=l({v})$, \ for all
$v\in\bm{\mathbf{\mathsf{V}}}(\Omega)$}.$ (9)

The existence and uniqueness of solution of the problem follows from the
complex version of the Lax-Milgram lemma [13, p. 376] (see also [56, Lem.
2.1.51, p. 40]).

Now, to solve Problem 1, we introduce the cost functional

$J(\Omega)=\frac{1}{2}\|u_{2}\|^{2}_{L^{2}(\Omega)}=\frac{1}{2}\int_{\Omega}{|u_{2}|^{2}}{\,{\operatorname{\mathnormal{d}}}x},$
(10)

where $u_{2}$ is, of course, subject to the state problem (8). The shape
optimization problem that we consider then is the problem of minimizing
$J(\Omega)$ over a set of admissible domains $\mathcal{O}_{\text{ad}}$. Here,
$\mathcal{O}_{\text{ad}}$ is essentially the set of $\mathcal{C}^{1,1}$
annular domains $\Omega$ with (fixed) interior boundary $\Gamma$ and (free)
exterior boundary $\Sigma$. It actually suffices to consider $\Gamma$ to be
only Lipschitz regular in deriving the (first-order) shape derivative of the
functions and shape functionals involved, but for simplicity, we also assume
it to be $\mathcal{C}^{1,1}$ regular. For second-order shape sensitivity
analysis, we will require the domain $\Omega$ be of class $\mathcal{C}^{2,1}$
– at least when applying the chain rule approach.

To numerically solve the optimization problem $J(\Omega)\to\inf$, we will
apply a shape-gradient-based descent method based on finite element method
(FEM). We will not, however, employ any kind of adaptive mesh refinement in
our numerical scheme as opposed to [51, 52, 53]. In this way, we can further
assess the stability of the new method in comparison with the Kohn-Vogelius
approach. Our method is of course different from [24, 25, 26] which make use
of the boundary element method to numerically solve the shape optimization
problem studied in the said papers. The expression for the shape derivative of
the cost will be exhibited in the next section using shape calculus [17, 40,
46, 57, 59]. We point out though that as opposed to [24, 25, 26] where the
chain rule approach is used to compute the shape derivative of the cost, our
approach does not use the strong form of the shape (Eulerian) derivative of
the state, but instead only make use of the weak form of the material
(Lagrangian) derivative of the state problem (see Remark 5). This technique
also applies to computing the second-order shape derivative of the cost which
requires the domain to be only of class $\mathcal{C}^{1,1}$.

Throughout the rest of the paper, we shall refer to the Kohn-Vogelius method
as KVM, or just KV, and to our proposed shape optimization method simply by
CCBM.

###  2.2 Motivations

To further motivate the study of the proposed shape optimization formulation
$J(\Omega)\to\inf$, we discuss here about some key advantages of CCBM over the
KVM (3) (and to other traditional approaches) while giving some additional
observations about the two formulations.

From the theoretical point of view, the computation of the first-order shape
derivative of the cost function $J$ is simpler compared to $J_{KV}$ when
computed via a Lagrange formulation or through the minimax formulation [14].
This is because the corresponding Lagrangian expression for CCBM is composed
of only one equality constraint (this corresponds to problem (6)) in addition
to the objective function $J$. In the case of KVM, the Lagrangian consists of
the cost function $J_{KV}$ and two equality constraints corresponding to a
mixed Dirichlet-Neumann problem and a pure Dirichlet problem. Moreover, in
order to take into account in the Lagrangian functional the Dirichlet data on
the free boundary for KVM, one needs to introduce a suitable Lagrangian
multiplier in constructing the functional. Such additional requirement is not
needed in the case of the CCBM formulation. In addition, when applying instead
the chain rule approach to get the shape derivative of the cost $J$ through
the shape derivative of the states, one only needs to compute and justify the
existence of the shape derivative of a single complex PDE in contrast to the
case of KVM where one has to deal with the shape derivative of two systems of
PDEs. Certainly, the latter would be more involved, especially when
considering more general and complicated boundary value problems. Also, we add
that the formulation – under suitable assumptions – offers more regularity for
the solution of the corresponding adjoint problem (see Remark 3) than in the
case of formulating (2) into a shape optimization setting that utilizes a
boundary-data-tracking-type cost functional. Meanwhile, from the numerical
viewpoint, CCBM appears to have nearly the same computational complexity with
respect to KVM. Nonetheless, it features some notable advantages in terms of
overall performance over KVM, but also presents some downsides. On the one
hand, CCBM only needs to solve one complex state problem to evaluate its
corresponding cost function. On the other hand, it also requires to solve
another PDE system (this corresponds to the adjoint problem (14) associated
with the formulation) in order to calculate its shape gradient. This is as
opposed to KVM which does not need the introduction of an adjoint problem.
Also, we stress here that, when applying FEM for instance, one does not need
to solve the coupled problems (7) and (8) simultaneously since it is enough to
solve the variational problem corresponding to (6). Even so, the time needed
to compute the solution to (6) is practically the same with the time required
to solve the two state problems associated with the KV formulation via a
finite element scheme. These observations may not seem to provide any
advantages to CCBM over KVM, but computational results – specifically obtained
under large deformations of the domains – provided in Examples 2–3 of
subsection 5.3 show otherwise. In fact, in many of these performed
experiments, CCBM seems to require less overall computational-time-per-
iteration of the approximate shape solution to (2) compared to KVM when solve
via a gradient-based scheme, not to mention that, in some instances, KVM
converges prematurely and sometimes tends to overshoot the optimal shape not
like CCBM. In general, KVM converges faster than CCBM and that the two methods
complete the same number of iterations at nearly the same time. Moreover, as
we observed in our numerical experiments, CCBM seems to have some sort of a
similar smoothing effect that has been observed in [51] when replacing the
Dirichlet boundary condition with a Robin boundary condition in a Neumann-
data-tracking approach. In addition, we become aware – after conducting
numerous numerical experiments – that the cost function $J_{KV}$ is less
sensitive compared to $J$. In fact, we notice that $J_{KV}(\Omega)$ is
insensitive to large deformations of $\Omega$. In some sense, this means that,
after some iterations, the $k$th approximant $\Omega^{k}$ of the exact shape
solution $\Omega^{\ast}$ gets stuck to some certain geometry and so will not
get closer to the exact solution even after an additional number of
iterations. The insensitivity of $J_{KV}$ to large variations limits the
algorithm to take smaller step sizes at each iterations – so to avoid being
stuck at certain point which may not yet be considered optimal. Surely, in
some situations, smaller step sizes may result to slower convergence behavior.
Be that as it may, the two methods provide almost the same optimal solution to
the given problem. These claims are made evident in subsection 5.3 where
numerical examples are provided. Last but not least, we point out that there
is another domain-integral-type penalization that may be considered. More
exactly, one may opt to penalize, instead of $J_{KV}$, by the cost functional

$J_{L^{2}}(\Omega):=\frac{1}{2}\int_{\Omega}{(u_{\text{N}}-u_{\text{D}})^{2}}{\,{\operatorname{\mathnormal{d}}}x}.$

Compared with $J_{KV}$, however, its shape gradient is more complex in
structure and would be numerically expensive to evaluate. In fact, it can be
checked that in this case one would need to solve five real BVPs: two state
equations, two adjoint equations (cf. [44]), and one to approximate the mean
curvature of the free boundary implicitly.

The discussion given above, as well as the numerical findings issued in
subsection 5.3, clearly warrants the proposal of the new method as a way to
solve stationary free boundary problems such as (1) and (2).

##  3 Shape Derivatives

This section is devoted to proving that the map $t\mapsto u^{t}$ is
$\mathcal{C}^{1}$ in a neighborhood of zero (Proposition 3.1), and to the
characterization of its derivative. Similar results for the map $t\mapsto
J(\Omega_{t})$ (Theorem 3.5) is also presented herein. On this purpose, we
need the notions of shape derivatives of functions and shape functionals, and
those we give in the first subsection. Meanwhile, a rigorous proof of
existence of the Lagrangian derivative of the state – presented specifically
in its weak form – will be given in subsection 3.2. Finally, the shape
derivative of the cost functional is exhibited in subsection 3.3.

###  3.1 Some elements of shape calculus

To accomplish our main tasks, we apply the concept of velocity or speed method
(see, e.g., [17, Chap. 4] or [16]). We let
$\mathcal{D}^{k}(\mathbb{R}^{d};\mathbb{R}^{d})$ be the space of $k$-times
continuously differentiable functions with compact support contained in
$\mathbb{R}^{d}$. Consider
$\bm{\mathbf{V}}\in\mathcal{E}^{k}:=\mathcal{C}([0,\varepsilon);\mathcal{D}^{k}(\mathbb{R}^{d};\mathbb{R}^{d}))$,
$\mathbb{N}\ni k\geqslant 2$, and let $\varepsilon>0$ be a small real number.
The field $\bm{\mathbf{V}}(t)(x)=\bm{\mathbf{V}}(t,x)$, $x\in\mathbb{R}^{d}$,
generates the transformations $T_{t}(\bm{\mathbf{V}})(X):=T_{t}(X)=x(t;X)$,
$t\geqslant 0$, $X\in\mathbb{R}^{d}$, through the differential equation
$\dot{x}(t;x_{0})=\bm{\mathbf{V}}(t,x(t;x_{0}))$, $x(0;x_{0})=x_{0}$, with the
initial value $x_{0}$ specified. We denote the “transformed domain”
$T_{t}(\bm{\mathbf{V}})(\Omega)$ at $t\geqslant 0$ by
$\Omega_{t}=:T_{t}(\Omega)$, and assumed that all admissible deformations of
$\Omega$ is contained in a larger open, bounded, and connected set
$U\subset\mathbb{R}^{d}$ (e.g., a ball in $\mathbb{R}^{d}$) of class
$\mathcal{C}^{k,1}$. In this work, the evolutions of the reference domain
$\Omega$ is described using time-independent velocity fields $\bm{\mathbf{V}}$
such that

$\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{k}:=\\{\bm{\mathbf{V}}\in\mathcal{C}^{k,1}(\overline{\Omega})^{d}\mid\bm{\mathbf{V}}=0\
\text{on}\ \Gamma\cup\partial U\\}.$ (11)

Here, we are fixing the interior boundary $\Gamma$ by taking
$\bm{\mathbf{V}}=0$ on $\Gamma$. We note that it is enough to consider
transformations $T_{t}$ that only change the position of the free part
$\Sigma$ of $\partial\Omega$, but do not rotate it. In other words, we may
just consider vector fields that have zero tangential part along $\Sigma$;
that is,
$\bm{\mathbf{V}}|_{\Sigma}=(\bm{\mathbf{V}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}=V_{n}\bm{\mathbf{n}}$.

We say that the function $u$ has a material derivative $\dot{u}$ and a shape
derivative $u^{\prime}$ at $0$ in the direction of the vector field
$\bm{\mathbf{V}}$ if the limits $\dot{u}=\lim_{t\searrow
0}\frac{1}{t}\left(u(\Omega_{t})\circ T_{t}-u(\Omega)\right)$ and
$u^{\prime}=\lim_{t\searrow
0}\frac{1}{t}\left(u(\Omega_{t})-u(\Omega)\right)$, exist, respectively, where
$(u(\Omega_{t})\circ T_{t})(x)=u(\Omega_{t})(T_{t}(x))$. These expressions are
related by $u^{\prime}=\dot{u}-(\nabla u\cdot\bm{\mathbf{V}})$ provided that
$\nabla u\cdot\bm{\mathbf{V}}$ exists in some appropriate function space [17,
59].

Given a shape functional $\mathsf{j}:\Omega\to\mathbb{R}$, we say that it has
a directional Eulerian derivative at $\Omega$ in the direction of
$\bm{\mathbf{V}}$ if the limit $\lim_{t\searrow
0}\frac{1}{s}\left(\mathsf{j}(\Omega_{t})-\mathsf{j}(\Omega)\right)=:\operatorname{\mathnormal{d}}\mathsf{j}(\Omega)[\bm{\mathbf{V}}]$
exists (cf. [17, Eq. (3.6), p. 172]). If the map
$\bm{\mathbf{V}}\mapsto\operatorname{\mathnormal{d}}\mathsf{j}(\Omega)[\bm{\mathbf{V}}]$
is linear and continuous, then $j$ is shape differentiable at $\Omega$, and
the map is referred to as the shape gradient of $j$. Similarly, the second-
order Eulerian derivative of $j$ at $\Omega$ along the two fields
$\bm{\mathbf{V}}$ and $\bm{\mathbf{W}}$ is given by $\lim_{s\searrow
0}\frac{1}{s}\left(\operatorname{\mathnormal{d}}\mathsf{j}(\Omega_{s}(\bm{\mathbf{W}}))[\bm{\mathbf{V}}]-\operatorname{\mathnormal{d}}\mathsf{j}(\Omega)[\bm{\mathbf{V}}]\right)=:{\operatorname{\mathnormal{d}}}^{2}\mathsf{j}(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]$
if the limit exists [15, Def. 2.3]. In addition, $\mathsf{j}$ is said to be
twice shape differentiable if, for all $\bm{\mathbf{V}}$ and
$\bm{\mathbf{W}}$,
${\operatorname{\mathnormal{d}}}^{2}\mathsf{j}(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]$
exists, and is bilinear and continuous with respect to
$\bm{\mathbf{V}},\bm{\mathbf{W}}$. Accordingly, we call the expression as the
shape Hessian of $\mathsf{j}$.

To complete our preparation, we also mention some properties of
$T_{t}:=T_{t}(\bm{\mathbf{V}})$, $\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{1}$,
that are essential to our investigation. For sufficiently small $t>0$, $T_{t}$
and its inverse $T_{t}^{-1}$ are both in
$\mathcal{D}^{1}(\mathbb{R}^{d};\mathbb{R}^{d})$ [17, Chap. 4] where $d$ is
the space dimension. Moreover, for $t\in[0,\varepsilon)$, $\varepsilon>0$ can
be chosen sufficiently small so that $I_{t}:=\det\,{\mathnormal{D}}T_{t}>0$.
Throughout the paper, we use the notations
$A_{t}:=I_{t}({\mathnormal{D}}T_{t}^{-1})({\mathnormal{D}}T_{t})^{-\top}$ and
$B_{t}=I_{t}|({\mathnormal{D}}T_{t})^{-\top}\bm{\mathbf{n}}|$, and assume
that, for all $t\in[0,\varepsilon)$, $I_{t}>0$ and there exist constants
$\Lambda_{0},\Lambda_{1}$ ($0<\Lambda_{0}<\Lambda_{1}$) such that
$\Lambda_{0}\leqslant B_{t}\leqslant\Lambda_{1}$ and
$\Lambda_{0}|\xi|^{2}\leqslant A_{t}\xi\cdot\xi\leqslant\Lambda_{1}|\xi|^{2}$
for all $\xi\in\mathbb{R}^{d}$. We also note the following regularities:
$[t\mapsto
I_{t}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}(\overline{\Omega}))$,
$[t\mapsto
A_{t}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}(\overline{\Omega})^{d\times
d})$, and $[t\mapsto
B_{t}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}(\Sigma))$, on a neighborhood
$\mathcal{I}$ of $0$ such that, for $t\in\mathcal{I}$, $T_{t}$ is a
diffeomorphism of $\Omega$ onto $\Omega_{t}$ (see, e.g., [42]). Lastly, we
note the following derivatives:
$\dot{I}_{t}\big{|}_{t=0}:=(d/dt)I_{t}\big{|}_{t=0}={\operatorname{div}}\,\bm{\mathbf{V}}$,
$\dot{A}_{t}\big{|}_{t=0}=A:=({\operatorname{div}}\bm{\mathbf{V}})\bm{\mathbf{I}}-\operatorname{\mathnormal{D}}\\!\bm{\mathbf{V}}-(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{V}})^{\top}$,
and
$\dot{B}_{t}\big{|}_{t=0}={\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}={\operatorname{div}}\bm{\mathbf{V}}\big{|}_{\Sigma}-(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{V}}\bm{\mathbf{n}})\cdot\bm{\mathbf{n}}$.

###  3.2 Lagrangian derivatives of the state

The main point of this section is to prove the following proposition
concerning the material derivative of the state. Throughout the section,
$\Omega$ is assumed to be of class $\mathcal{C}^{1,1}$ and
$\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{1}$. This implies that, for
sufficiently small $t>0$, $T_{t}$ is a $\mathcal{C}^{1,1}$ diffeomorphism of
$\Omega$ onto $\Omega_{t}$.

######  Proposition 3.1.

Let $\Omega$ be of class $\mathcal{C}^{1,1}$ and
$\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{1}$. The map $t\mapsto
u^{t}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)$ is $\mathcal{C}^{1}$ in a
neighborhood of $0$. Its Lagrangian derivative at $0$, denoted by $\dot{u}$,
belongs to $\bm{\mathbf{\mathsf{V}}}(\Omega)$ and satisfies

$\displaystyle\int_{\Omega}{\nabla{\dot{u}}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{\dot{u}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle=-\int_{\Omega}{A\nabla{u}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{({\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}){u}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\qquad+\lambda\int_{\Sigma}{({\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}){\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma},\quad\forall{v}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$

To prove the above result, we need the following lemmas.

######  Lemma 3.2.

The sesquilinear form $\bm{\mathbf{\mathsf{a}}}_{t}$ defined on
$\bm{\mathbf{\mathsf{V}}}(\Omega)\times\bm{\mathbf{\mathsf{V}}}(\Omega)$ by

$\bm{\mathbf{\mathsf{a}}}_{t}({u},{v})=\int_{\Omega}{A_{t}\nabla{u}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{B_{t}{u}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma},\quad\forall{u},{v}\in\bm{\mathbf{\mathsf{V}}}(\Omega),$

is bounded and coercive on
$\bm{\mathbf{\mathsf{V}}}(\Omega)\times\bm{\mathbf{\mathsf{V}}}(\Omega)$ for
sufficiently small $t>0$.

###### Proof.

For ${u}$, ${v}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$ and sufficiently small
$t>0$, we have, by Cauchy-Schwarz and trace inequalities, the estimate
$|\bm{\mathbf{\mathsf{a}}}_{t}({u},{v})|\leqslant
c_{t}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{V}}}(\Omega)}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{V}}}(\Omega)}$
where $c_{t}=\max\left(|A_{t}|_{\infty},|B_{t}|_{\infty}\right)$. This shows
that $\bm{\mathbf{\mathsf{a}}}_{t}$ is bounded. Also, we have that222Here,
$\bm{\mathbf{I}}$ stands for the identity matrix in $d=2$ dimensions.

$\bm{\mathbf{\mathsf{a}}}_{t}({u},{u})=\int_{\Omega}{\left\\{\nabla
u\cdot\nabla\overline{u}+(A_{t}-\bm{\mathbf{I}})\nabla{u}\cdot\nabla{\overline{u}}\right\\}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{\left\\{u\overline{u}+(B_{t}-1){u}\overline{u}\right\\}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$

Because the maps $t\mapsto A_{t}$ and $t\mapsto B_{t}$ are continuous at $0$,
then – for sufficiently small $t$ –
$\max\left(|A_{t}-\bm{\mathbf{I}}|_{\infty},|B_{t}-1|_{\infty}\right)<1$.
Consequently,
$\Re\\{\bm{\mathbf{\mathsf{a}}}_{t}({u},{u})\\}\gtrsim\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\nabla{u}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\bm{\mathbf{\mathsf{Q}}}}+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{\bm{\mathbf{\mathsf{S}}}}\right)$.
Hence, $\bm{\mathbf{\mathsf{a}}}_{t}$ is coercive for small enough $t$. ∎

######  Lemma 3.3.

The function $u^{t}=u_{1}^{t}+iu_{2}^{t}$ uniquely solves in
$\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)$ the equation

$\displaystyle\int_{\Omega}{A_{t}\nabla{u^{t}}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{B_{t}{u^{t}}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle=\lambda\int_{\Sigma}{B_{t}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma},\
\ \forall{v}\in\bm{\mathbf{\mathsf{V}}}(\Omega),$ (12) $\displaystyle u^{t}$
$\displaystyle=1\ \text{on}\ \Gamma.$

###### Proof.

The function ${u_{t}}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega_{t})$ satisfies
${u_{t}}=1$ on $\Gamma$, and solves the variational problem

$\int_{\Omega_{t}}{\nabla{u_{t}}\cdot\nabla\overline{\varphi}_{t}}{\,{\operatorname{\mathnormal{d}}}x_{t}}+i\int_{\Sigma_{t}}{{u_{t}}\overline{\varphi}_{t}}{\,{\operatorname{\mathnormal{d}}}\sigma_{t}}=\lambda\int_{\Sigma_{t}}{\overline{\varphi}_{t}}{\,{\operatorname{\mathnormal{d}}}\sigma_{t}},\quad\forall\varphi_{t}\in\bm{\mathbf{\mathsf{V}}}(\Omega_{t}),$

where
$\bm{\mathbf{\mathsf{V}}}(\Omega_{t})=\\{\varphi_{t}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega_{t})\mid\varphi_{t}=0\
\text{on}\ \Gamma\\}$. Using the relation $u^{t}={u_{t}}\circ
T_{t}:={u_{1t}}\circ T_{t}+i{u_{2t}}\circ T_{t}$, the identity
$(\nabla\varphi_{t})\circ T_{t}=DT_{t}^{-\top}\nabla\varphi^{t}$ which holds
for any $\varphi_{t}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega_{t})$ and
$\varphi^{t}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)$, and the change of
variables (cf. [17, subsec. 9.4.2–9.4.3, pp. 482–484]), we get ${u^{t}}=1$ on
$\Gamma$ and the variational equation above transforms to

$\int_{\Omega}{A_{t}\nabla{u^{t}}\cdot\nabla\overline{\varphi}^{t}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{B_{t}{u^{t}}\overline{\varphi}^{t}}{\,{\operatorname{\mathnormal{d}}}\sigma}=\lambda\int_{\Sigma}{B_{t}\overline{\varphi}^{t}}{\,{\operatorname{\mathnormal{d}}}\sigma},\quad\forall\varphi^{t}\in\bm{\mathbf{\mathsf{V}}}(\Omega).$

We immediately get (12) by taking $v=\varphi^{t}$ above. Now, from Lemma 3.2,
$\bm{\mathbf{\mathsf{a}}}_{t}(\cdot,\cdot):[\bm{\mathbf{\mathsf{V}}}(\Omega)]^{2}\to\mathbb{R}$
is bounded and coercive. We let $u_{0}\in\bm{\mathbf{\mathsf{H}}}^{1}(U)$ be a
fixed function such that $u_{0}=1$ on $\Gamma$. Then,
$u^{t}-u_{0}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$, and by (12), we have

$\displaystyle\int_{\Omega}{A_{t}\nabla{(u^{t}-u_{0})}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{B_{t}{(u^{t}-u_{0})}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\qquad=-\int_{\Omega}{A_{t}\nabla
u_{0}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{B_{t}u_{0}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}+\lambda\int_{\Sigma}{B_{t}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma},\quad\forall{v}\in\bm{\mathbf{\mathsf{V}}}(\Omega).$

For $t>0$ small enough, the following estimates hold

$\displaystyle\left|\int_{\Sigma}{B_{t}u_{0}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}\right|$
$\displaystyle\lesssim|B_{t}|_{\infty}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u_{0}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{H}}}^{1}(U)}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)},$
$\displaystyle\left|\lambda\int_{\Sigma}{B_{t}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}\right|$
$\displaystyle\lesssim|B_{t}|_{\infty}|\Sigma|^{1/2}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)},$
$\displaystyle\left|\int_{\Omega}{A_{t}\nabla
u_{0}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}\right|$
$\displaystyle\leqslant|A_{t}|_{\infty}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u_{0}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{H}}}^{1}(U)}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)}.$

By Lax-Milgram lemma, $z^{t}=u^{t}-u_{0}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$
is the unique solution in $\bm{\mathbf{\mathsf{V}}}(\Omega)$ of

$\displaystyle\int_{\Omega}{A_{t}\nabla{z^{t}}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{B_{t}{z^{t}}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\qquad=-\int_{\Omega}{A_{t}\nabla
u_{0}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{B_{t}u_{0}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}+\lambda\int_{\Sigma}{B_{t}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma},\quad\forall{v}\in\bm{\mathbf{\mathsf{V}}}(\Omega).$

Now, let $u^{t}=z^{t}+u_{0}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)$. Then, we
have

$\displaystyle\int_{\Omega}{A_{t}\nabla
u^{t}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{B_{t}u^{t}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\quad\qquad\qquad=\int_{\Omega}{A_{t}\nabla(z^{t}+u_{0})\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{B_{t}(z^{t}+u_{0}){\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\quad\qquad\qquad=\lambda\int_{\Sigma}{B_{t}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma},\quad\forall{v}\in\bm{\mathbf{\mathsf{V}}}(\Omega).$

Clearly, $u^{t}=z^{t}+1=1$ on $\Gamma$ because
$z^{t}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$. Uniqueness of $u^{t}$ follows from
the uniqueness of $z^{t}$. That is, $u^{t}$ uniquely solves (12) in
$\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)$. ∎

######  Lemma 3.4.

The map $t\mapsto{u^{t}}$ is $\mathcal{C}^{1}$ in a neighborhood of $0$.

###### Proof.

We prove the statement using the implicit function theorem (IFT). In view of
(12), we see that ${u^{t}}-{u}$ is the unique element in
$\bm{\mathbf{\mathsf{V}}}(\Omega)$ that satisfies

$\displaystyle\int_{\Omega}{A_{t}\nabla({u^{t}}-{u})\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{B_{t}({u^{t}}-{u}){\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\qquad\qquad=-\int_{\Omega}{A_{t}\nabla{u}\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{B_{t}{u}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}+\lambda\int_{\Sigma}{B_{t}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma},\quad\forall{v}\in\bm{\mathbf{\mathsf{V}}}(\Omega).$

Denoting by $\langle\cdot,\cdot\rangle$ the duality pairing between
$\bm{\mathbf{\mathsf{V}}}(\Omega)$ and its dual space
$\bm{\mathbf{\mathsf{V}}}^{\prime}(\Omega)$, we consider the function
${\mathcal{F}}:\mathcal{I}\times\bm{\mathbf{\mathsf{V}}}(\Omega)\to\bm{\mathbf{\mathsf{V}}}^{\prime}(\Omega)$
defined by

$\langle{\mathcal{F}}(t,w),v\rangle:=\int_{\Omega}{A_{t}\nabla(w+{u})\cdot\nabla{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}x}+{i}\int_{\Sigma}{B_{t}(w+{u}){\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma}-\lambda\int_{\Sigma}{B_{t}{\overline{v}}}{\,{\operatorname{\mathnormal{d}}}\sigma},$

for all $v$, $w\in\bm{\mathbf{\mathsf{V}}}(\Omega)$. Clearly, ${\mathcal{F}}$
is $\mathcal{C}^{1}$ because $t\mapsto A_{t}$ and $t\mapsto B_{t}$ are
$\mathcal{C}^{1}$ in a neighborhood of $0$. We observe that ${u^{t}}-{u}$
uniquely solves ${\mathcal{F}}(t,{u^{t}}-{u})=0$ in
$\bm{\mathbf{\mathsf{V}}}(\Omega)$. In addition,
$\langle\operatorname{\mathnormal{D}}\\!_{w}{\mathcal{F}}(0,0)w,v\rangle=\bm{\mathbf{\mathsf{a}}}_{t}(w,v)$.
Using Lemma 3.2, we deduce via the complex version of Lax-Milgram lemma that
$\operatorname{\mathnormal{d}}\\!_{w}{\mathcal{F}}(0,0)$ is an isomorphism
from $\bm{\mathbf{\mathsf{V}}}(\Omega)$ to
$\bm{\mathbf{\mathsf{V}}}^{\prime}(\Omega)$. By the IFT, we conclude that the
map $t\mapsto{u^{t}}-{u}$ is $\mathcal{C}^{1}$ in a neighborhood of $0$. Now
let $\dot{u}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$ be its derivative at $t=0$.
Differentiating ${\mathcal{F}}(t,{u^{t}}-{u})=0$ with respect to $t$, we
obtain
$\left\langle\operatorname{\mathnormal{D}}\\!_{w}{\mathcal{F}}(0,0)\dot{u},v\right\rangle+\left\langle\frac{\partial}{\partial
t}{\mathcal{F}}(0,0),v\right\rangle=0$, for all
$v\in\bm{\mathbf{\mathsf{V}}}(\Omega)$, leading to (3.1). ∎

###  3.3 Computation of the shape gradient

We now prove the differentiability of the map $t\mapsto J(\Omega_{t})$ and
characterize its derivative.

######  Theorem 3.5 (Shape gradient of $J$).

Let $\Omega$ be of class $\mathcal{C}^{1,1}$ and
$\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{1}$. The map $t\mapsto J(\Omega_{t})$
is $\mathcal{C}^{1}$ in a neighborhood of $0$, and its derivative at $0$ is
given by
${\operatorname{\mathnormal{d}}}J(\Omega)[\bm{\mathbf{V}}]=\int_{\Sigma}{G\nu\cdot\bm{\mathbf{V}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
where

$\displaystyle G$
$\displaystyle=\frac{1}{2}|u_{2}|^{2}-\Big{[}\nabla_{\Sigma}p_{1}\cdot\nabla_{\Sigma}u_{2}-\nabla_{\Sigma}p_{2}\cdot\nabla_{\Sigma}u_{1}$
(13)
$\displaystyle\qquad\qquad\qquad+p_{1}\left(\partial_{\bm{\mathbf{n}}}{u_{1}}+\kappa
u_{1}\right)+p_{2}\left(\partial_{\bm{\mathbf{n}}}{u_{2}}+\kappa
u_{2}\right)+\lambda\kappa p_{2}\Big{]},$

$\kappa$ stands for the mean curvature of the free boundary $\Sigma$,
$u=u_{1}+iu_{2}$ is the unique solution to (6), and $p=p_{1}+ip_{2}$ uniquely
solves the adjoint system

$-\Delta p=u_{2}\ \text{in $\Omega$},\qquad\quad p=0\ \text{on
$\Gamma$},\qquad\quad\partial_{\bm{\mathbf{n}}}{p}-ip=0\ \text{on $\Sigma$}.$
(14)

######  Remark 2.

The weak formulation of (14) reads as follows:

$\text{ find $p\in\bm{\mathbf{\mathsf{V}}}(\Omega)$ such
that}\int_{\Omega}{\nabla
p\cdot\nabla{\overline{\varphi}}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{p{\overline{\varphi}}}{\,{\operatorname{\mathnormal{d}}}\sigma}=\int_{\Omega}{u_{2}{\overline{\varphi}}}{\,{\operatorname{\mathnormal{d}}}x},\
\forall\varphi\in\bm{\mathbf{\mathsf{V}}}(\Omega).$ (15)

The existence and uniqueness of solution to the variational problem (15) is
again a consequence of the Lax-Milgram lemma.

######  Remark 3.

In general, for $\Omega$ of class $\mathcal{C}^{k+1,1}$, $k$ a non-negative
integer, it can be shown that the weak solution
$u\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)$ to the variational problem (9) is
also $\bm{\mathbf{\mathsf{H}}}^{k+2}(\Omega)$. In particular, $u_{2}\in
H^{k+2}(\Omega)$. Consequently, we find that the weak solution $p$ of problem
(15) is not only in $\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)$, but is also an
element of $\bm{\mathbf{\mathsf{H}}}^{k+4}(\Omega)$.

Theorem 3.5 relies on following lemma. The
$\bm{\mathbf{\mathsf{H}}}^{2}(\Omega)$ regularity of ${u}$, which holds true
since $\Omega$ is assumed to be $\mathcal{C}^{1,1}$ and
$\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{1}$, will be used subsequently
without further notice.

######  Lemma 3.6.

Let $\Omega$ be of class $\mathcal{C}^{1,1}$ and
$\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{1}$. The solution ${u}$ of the state
problem (6) and the adjoint variable $p$ which is the solution to (14) satisfy
the equation

$\displaystyle\int_{\Omega}{A\nabla{u}\cdot\nabla\overline{p}}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle=-\int_{\Omega}{u_{2}\bm{\mathbf{V}}\cdot\nabla{u}}{\,{\operatorname{\mathnormal{d}}}x}+\int_{\Sigma}{(\nabla\overline{{p}}\cdot\nabla{u})V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}+i\int_{\Sigma}{{u}(\bm{\mathbf{V}}\cdot\nabla\overline{{p}})}{\,{\operatorname{\mathnormal{d}}}\sigma}$
(16)
$\displaystyle\qquad\quad+i\int_{\Sigma}{\overline{p}(\bm{\mathbf{V}}\cdot\nabla{u})}{\,{\operatorname{\mathnormal{d}}}\sigma}-\lambda\int_{\Sigma}{(\bm{\mathbf{V}}\cdot\nabla\overline{{p}})}{\,{\operatorname{\mathnormal{d}}}\sigma}.$

###### Proof.

The result follows from the formula

$\displaystyle\int_{\Omega}{A\nabla\varphi\cdot\nabla\overline{\psi}}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle=\int_{\Omega}{(\Delta\varphi)\bm{\mathbf{V}}\cdot\nabla\overline{\psi}}{\,{\operatorname{\mathnormal{d}}}x}+\int_{\Omega}{(\Delta\overline{\psi})\bm{\mathbf{V}}\cdot\nabla\varphi}{\,{\operatorname{\mathnormal{d}}}x}-\int_{\Sigma}{\partial_{\bm{\mathbf{n}}}{\varphi}(\bm{\mathbf{V}}\cdot\nabla\overline{\psi})}{\,{\operatorname{\mathnormal{d}}}\sigma}$
(17)
$\displaystyle\quad\qquad-\int_{\Sigma}{\partial_{\bm{\mathbf{n}}}{\overline{\psi}}(\bm{\mathbf{V}}\cdot\nabla\varphi)}{\,{\operatorname{\mathnormal{d}}}\sigma}+\int_{\Sigma}{(\nabla\overline{\psi}\cdot\nabla\varphi)V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma},$

which holds for all functions $\varphi$,
$\psi\in\bm{\mathbf{\mathsf{V}}}(\Omega)\cap\bm{\mathbf{\mathsf{H}}}^{2}(\Omega)$
and $\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{1}$, where
$A=({\operatorname{div}}\bm{\mathbf{V}})\bm{\mathbf{I}}-\operatorname{\mathnormal{D}}\\!\bm{\mathbf{V}}-(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{V}})^{\top}$
(see, e.g., [17]). Letting $\varphi={u}$ and $\psi=p$, and noting that
$\partial_{\bm{\mathbf{n}}}{\overline{p}}=-i\overline{p}$ and
$\partial_{\bm{\mathbf{n}}}{u}=-i{u}+\lambda$ on $\Sigma$ and that
$-\Delta\overline{p}=u_{2}$ in $\Omega$, we obtain

$\displaystyle\int_{\Omega}{A\nabla{u}\cdot\nabla\overline{p}}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle=-\int_{\Omega}{u_{2}\bm{\mathbf{V}}\cdot\nabla{u}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{{u}(\bm{\mathbf{V}}\cdot\nabla\overline{{p}})}{\,{\operatorname{\mathnormal{d}}}\sigma}+i\int_{\Sigma}{\overline{p}(\bm{\mathbf{V}}\cdot\nabla{u})}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\qquad-\lambda\int_{\Sigma}{(\bm{\mathbf{V}}\cdot\nabla\overline{{p}})}{\,{\operatorname{\mathnormal{d}}}\sigma}+\int_{\Sigma}{(\nabla\overline{{p}}\cdot\nabla{u})V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma},$

as desired. This proves the lemma. ∎

###### Proof of Theorem 3.5.

By change of variables, we write
$J(\Omega_{t})=\frac{1}{2}\int_{\Omega}{I_{t}|u_{2}^{t}|^{2}}{\,{\operatorname{\mathnormal{d}}}x}$.
Because $[t\mapsto
I_{t}]\in\mathcal{C}^{1}(\mathcal{I},\mathcal{C}(\overline{\Omega}))$ and
$[t\mapsto
u^{t}]\in\mathcal{C}^{1}(\mathcal{I},\bm{\mathbf{\mathsf{H}}}^{1}(\Omega))$,
then the mapping $t\mapsto J(\Omega_{t})$ is also $\mathcal{C}^{1}$ in the
neighborhood $\mathcal{I}$ of $0$, and we have

${\operatorname{\mathnormal{d}}}J(\Omega)[\bm{\mathbf{V}}]=\frac{1}{2}\int_{\Omega}{({\operatorname{div}}\,\bm{\mathbf{V}})|u_{2}|^{2}}{\,{\operatorname{\mathnormal{d}}}x}+\int_{\Omega}{u_{2}\dot{u}{}_{2}}{\,{\operatorname{\mathnormal{d}}}x}.$
(18)

On the one hand, using the identity
$-\int_{\Omega}{({\operatorname{div}}\,\bm{\mathbf{v}})\phi}{\,{\operatorname{\mathnormal{d}}}x}=\int_{\Omega}{\bm{\mathbf{v}}\cdot\nabla\phi}{\,{\operatorname{\mathnormal{d}}}x}-\int_{\Sigma}{\phi(\bm{\mathbf{v}}\cdot\nu)}{\,{\operatorname{\mathnormal{d}}}\sigma}$,
which holds (on bounded Lipschitz domain $\Omega$) for any vector field
$\bm{\mathbf{v}}\in\mathcal{C}^{1}(\overline{\Omega})^{2}$ and scalar function
$\phi\in W^{1,1}(\Omega)$, the first integral can be equivalently written as

$\frac{1}{2}\int_{\Omega}{({\operatorname{div}}\,\bm{\mathbf{V}})|u_{2}|^{2}}{\,{\operatorname{\mathnormal{d}}}x}=-\int_{\Omega}{u_{2}\bm{\mathbf{V}}\cdot\nabla
u_{2}}{\,{\operatorname{\mathnormal{d}}}x}+\frac{1}{2}\int_{\Sigma}{|u_{2}|^{2}V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$
(19)

On the other hand, the second integral is actually not useful for practical
applications, especially in the numerical realization of the present shape
minimization problem because it requires the solution of (3.1) for each
velocity field $\bm{\mathbf{V}}$. A way to resolve this issue is to rewrite
the integrand in terms of another variable – getting rid of the term
$\dot{u}{}$ – through the adjoint method. To this end, we utilize the
variational problem (15) with the test function
$\varphi=\dot{u}{}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$ to obtain

$\int_{\Omega}{u_{2}\overline{\dot{u}{}}}{\,{\operatorname{\mathnormal{d}}}x}=\int_{\Omega}{\nabla
p\cdot\nabla\overline{\dot{u}{}}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{p\overline{\dot{u}{}}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$
(20)

Now, choosing $v=p\in\bm{\mathbf{\mathsf{V}}}(\Omega)$ as the test function in
(3.1) gives us

$\displaystyle\int_{\Omega}{\nabla{\dot{u}}\cdot\nabla{\overline{p}}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma}{\dot{u}{\overline{p}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle=-\int_{\Omega}{A\nabla{u}\cdot\nabla{\overline{p}}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{({\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}){u}{\overline{p}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
(21)
$\displaystyle\qquad+\lambda\int_{\Sigma}{({\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}){\overline{p}}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$

Hence, taking the complex conjugate of both sides of (20) and then comparing
the resulting equation to (21) leads us to the identity

$\displaystyle\int_{\Omega}{u_{2}\dot{u}{}}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle=-\int_{\Omega}{A\nabla{u}\cdot\nabla{\overline{p}}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{({\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}){u}{\overline{p}}}{\,{\operatorname{\mathnormal{d}}}\sigma}+\lambda\int_{\Sigma}{({\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}){\overline{p}}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$

Applying (16) in Lemma 3.6, we can further write

$\displaystyle\int_{\Omega}{u_{2}\dot{u}{}}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle=\int_{\Omega}{u_{2}\bm{\mathbf{V}}\cdot\nabla{u}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{\overline{p}(\bm{\mathbf{V}}\cdot\nabla{u})}{\,{\operatorname{\mathnormal{d}}}\sigma}-i\int_{\Sigma}{{u}(\bm{\mathbf{V}}\cdot\nabla\overline{{p}})}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\qquad-i\int_{\Sigma}{{u}{\overline{p}}\,{\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}}{\,{\operatorname{\mathnormal{d}}}\sigma}-\int_{\Sigma}{(\nabla\overline{{p}}\cdot\nabla{u})V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\qquad+\lambda\int_{\Sigma}{\nabla\overline{{p}}\cdot\bm{\mathbf{V}}}{\,{\operatorname{\mathnormal{d}}}\sigma}+\lambda\int_{\Sigma}{{\overline{p}}\,{\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$

At this point we apply the following version of the tangential Green’s
formula333A proof of this formula can be found in [46]., which is valid, for
instance, when $\Sigma$ is $\mathcal{C}^{1,1}$,

$\int_{\Sigma}{\left(\nabla\phi\cdot\bm{\mathbf{V}}+\phi\,{\operatorname{div}}_{\Sigma}\bm{\mathbf{V}}\right)}{\,{\operatorname{\mathnormal{d}}}\sigma}=\int_{\Sigma}{\left(\partial_{\bm{\mathbf{n}}}{\phi}+\phi\,{\operatorname{div}}_{\Sigma}\bm{\mathbf{n}}\right){V_{n}}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$
(22)

Here, the function $\phi$ is supposed to be $W^{2,1}(U)$ regular. We also note
the fact that $\kappa={\operatorname{div}}_{\Sigma}\bm{\mathbf{n}}$, where
$\kappa$ is the mean curvature of $\Sigma$, and utilize the identity
$\nabla\overline{{p}}\cdot\nabla{u}=\nabla_{\Sigma}\overline{{p}}\cdot\nabla_{\Sigma}{u}+\partial_{\bm{\mathbf{n}}}{\overline{p}}\partial_{\bm{\mathbf{n}}}{u}$444Here
$\nabla_{\Sigma}$ is the tangential gradient operator on $\Sigma$. The
intrinsic definition of the operator is given in [17, Chap. 5., Sec. 5.1, p.
492]. More discussion on tangential calculus can be found in the same
referenced text. on $\Sigma$ to obtain

$\displaystyle\int_{\Omega}{u_{2}\dot{u}{}}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle=\int_{\Omega}{u_{2}\bm{\mathbf{V}}\cdot\nabla{u}}{\,{\operatorname{\mathnormal{d}}}x}-i\int_{\Sigma}{\overline{p}\left(\partial_{\bm{\mathbf{n}}}{u}+\kappa
u\right)V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle\qquad-\int_{\Sigma}{\nabla_{\Sigma}\overline{{p}}\cdot\nabla_{\Sigma}{u}V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}+\lambda\int_{\Sigma}{\kappa\overline{p}V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$

The desired expression immediately follows by comparing the respective complex
parts on both sides of the above equation, and together with (19). ∎

The next conclusion can be drawn easily from (13), (14), and Remark 1.

######  Corollary 3.7 (Necessary condition).

Let the domain $\Omega^{\ast}$ be such that the state $u=u(\Omega^{\ast})$
satisfies the overdetermined boundary-value problem (2), i.e., there holds

$u_{2}=0\quad\text{on $\Omega^{\ast}$}.$ (23)

Then, the domain $\Omega^{\ast}$ is stationary for the shape problem
$\frac{1}{2}\int_{\Omega}{|u_{2}|^{2}}{\,{\operatorname{\mathnormal{d}}}x}\to\inf$,
where $u_{2}$ is subject to (8). That is, it fulfills the necessary optimality
condition

${\operatorname{\mathnormal{d}}}J(\Omega^{\ast})[\bm{\mathbf{V}}]=0,\quad\text{for
all $\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{2}$}.$ (24)

###### Proof.

By the assumption that $u_{2}=0$ on $\Omega^{\ast}$, one finds that
$p=p(\Omega^{\ast})$ on $\Omega^{\ast}$. Thus, it follows that $G\equiv 0$ on
$\Sigma^{\ast}$ which implies that
${\operatorname{\mathnormal{d}}}J(\Omega^{\ast})[\bm{\mathbf{V}}]=0$, for any
$\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{1}$. ∎

######  Remark 4.

In connection with the previous result, we remark that solutions of the
necessary condition (24) might exist such that the state does not satisfy
equation (23). However, only in the case of exact matching of boundary data a
stationary domain $\Omega^{\ast}$ is a global minimum because
$J(\Omega^{\ast})=0$.

######  Remark 5.

To obtain the form
${\operatorname{\mathnormal{d}}}J(\Omega)[\bm{\mathbf{V}}]=\int_{\Sigma}{G\nu\cdot\bm{\mathbf{V}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$
– in accordance with Hadamard–Zolésio’s structure theorem – via the chain rule
approach, the authors in [24, 25, 26, 52, 53] utilize the strong form of the
equation as well as the boundary conditions satisfied by the shape derivative
of $u$. The existence of the said derivative, however, requires more
regularity on $u$, in addition to assuming that the domain is at least
$\mathcal{C}^{2,\alpha}$ ($\alpha\in(0,1]$). Allowing this higher regularity
assumptions and taking the shape differentiability of the state $u$ for
granted, an alternative and more direct argument to derive the shape
derivative of $J$ could be given, see Appendix A.2. The derivation of the
shape gradient of the cost issued above is based on the Lagrangian derivative
$\dot{u}{}$ which is rigorously accounted for in the proof of Theorem 3.1. We
emphasize that only the weak form of the Lagrangian derivatives is utilized.
This bypasses the need to use the strong form of the shape derivative of $u$.
Alternatively, Theorem 3.5 could be shown via the strong form of the equation,
and of the Eulerian derivatives of the state, as well as the boundary
conditions satisfied by the shape derivative of ${u}$ coupled with Hadamard’s
domain differentiation formula (see, e.g., [17, Thm. 4.2, p. 483]), [40, eq.
(5.12), Thm. 5.2.2, p. 194] or [59, eq. (2.168), p. 113]):

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

##  4 Instability analysis of the critical shape

In this section, we investigate the question of stability of the proposed
shape optimization formulation of (2) at a critical shape $\Omega^{\ast}$. To
do this, we need the expression for the shape Hessian of $J$ at
$\Omega^{\ast}$ which we exhibit in subsection 4.2. Afterwards, we deal with
the stability condition of the present shape optimization problem in
subsection 4.3. A quick review of the sufficient second-order optimality
condition for the present shape optimization problem is given first in the
following subsection.

###  4.1 Sufficient Conditions

Before we proceed to the computation of the shape Hessian, we briefly provide
here a rough discussion about the regular local optimality of second-order of
the critical shape $\Omega^{\ast}$ (the optimal solution of the proposed shape
problem) on the basis of [27, Thm. 3.3]. Hereinafter, whenever we state the
phrase ‘$\mathcal{C}^{\circ}$ regular/smooth’ we mean that the object it
describes is either $\mathcal{C}^{0,1}\cap\mathcal{C}^{2,\alpha}$ regular,
$\alpha>0$, or simply $\mathcal{C}^{2,1}$ regular. Similar definition is given
for $\mathcal{C}^{\circ}(\Omega)$ and $\mathcal{C}^{\circ}(\Sigma)$. Moreover,
$\bm{\mathbf{\mathcal{C}}}^{\circ}(\cdot):=[\mathcal{C}^{\circ}(\cdot)]^{2}$
(e.g.,
$\bm{\mathbf{\mathcal{C}}}^{2,\alpha}(\Omega):=[\mathcal{C}^{2,\alpha}(\Omega)]^{2}$).
In the sequel, we say that $\Omega\in\mathcal{C}^{\circ}$ is an admissible
perturbation of $\widehat{\Omega}\in\mathcal{C}^{\circ}$ if for some fixed
small number $\delta\in(0,1)$ the following inequality condition holds:

$\|\phi-\hat{\phi}\|_{\mathcal{C}^{\circ}(\mathbb{R};\mathbb{R}^{2})}<\delta,$

where $\phi,\hat{\phi}\in\mathcal{C}^{\circ}(\mathbb{R};\mathbb{R}^{2})$ are
respectively a periodic parametrization of the free boundary of $\Omega$ and
$\widehat{\Omega}$. We denote the collection of such perturbations of
$\widehat{\Omega}$ by $U_{\delta}({\hat{\phi}})$. Note that, essentially,
$U_{\delta}({\hat{\phi}})$ is the set of all perturbations $\Omega$ that are
‘sufficiently near’ – based on some appropriate metric – to
$\widehat{\Omega}$. In fact, in the next proposition, the result only holds
for all perturbations $\Omega$ that lies in some specific vicinity of the
stationary shape $\Omega^{\ast}$ (i.e., domains whose free boundaries are
contained within some certain tubular neighborhood of $\Sigma^{\ast}$). We let
$U:=U_{1}(\cdot)$ be another open bounded set containing all sets
$U_{\delta}(\cdot)$, $\delta\in(0,1)$.

######  Proposition 4.1 (sufficient second-order optimality condition).

Let the necessary condition

${\operatorname{\mathnormal{d}}}J(\Omega^{\ast})[\bm{\mathbf{V}}]=0,\quad\text{for
all
$\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{2}\cap\bm{\mathbf{\mathcal{C}}}^{\circ}(\overline{U})$},$
(A1)

holds for a certain critical shape $\Omega^{\ast}\in\mathcal{C}^{\circ}$. For
all admissible perturbation $\Omega=\Omega(\phi)\in U_{\delta}(\phi^{\ast})$
of $\Omega^{\ast}$, we suppose that there is a constant $c_{b}>0$ which
depends continuously on $\Omega$ via the parametric function $\phi$ such that
the bilinear form imposed by the shape Hessian satisfies the following
inequality

$|{\operatorname{\mathnormal{d}}}^{2}J(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]|\leqslant
c_{b}(\Omega)\|\bm{\mathbf{V}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)}\|\bm{\mathbf{W}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)},$
(A2)

if $\Omega\in\overline{U_{\delta}(\phi^{\ast})}$, and the remainder estimate

$|{\operatorname{\mathnormal{d}}}^{2}J(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]-{\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]|\leqslant\eta\left(\|\phi-\phi^{\ast}\|_{\mathcal{C}^{\circ}(\mathbb{R};\mathbb{R}^{2})}\right)\|\bm{\mathbf{V}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)}\|\bm{\mathbf{W}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)},$
(A3)

for all
$\bm{\mathbf{V}},\bm{\mathbf{W}}\in\bm{\mathbf{\Theta}}^{2}\cap\bm{\mathbf{\mathcal{C}}}^{\circ}(\overline{U})$,
where $\eta:\mathbb{R}^{+}_{0}\to\mathbb{R}^{+}_{0}$ is a decreasing function
that satisfies the condition $\eta(s)\to 0$ as $s\to 0$. Then, the domain
$\Omega^{\ast}$ is a strong regular local optimum of second-order with respect
to a specific constant $\hat{c}_{e}>0$,

$J(\Omega)-J(\Omega^{\ast})\leqslant\hat{c}_{e}\|\phi-\phi^{\ast}\|_{\bm{\mathbf{H}}^{1}(\Sigma^{\ast})},\quad\text{for
all $\Omega\in\overline{U_{\hat{\delta}}(\phi^{\ast})}$}.$ (26)

if and only if the shape Hessian satisfies the strong coercivity estimate

${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{V}}]\geqslant
c_{e}\|\bm{\mathbf{V}}\|^{2}_{\bm{\mathbf{H}}^{1}(\Sigma^{\ast})},\quad\text{for
all
$\bm{\mathbf{V}}\in\bm{\mathbf{\Theta}}^{2}\cap\bm{\mathbf{\mathcal{C}}}^{\circ}(\overline{U})$},$
(A4)

for some constant $c_{e}>0$.

As remarked in many references (see, e.g., [21, 22, 27, 23, 24, 25, 26]), it
is generally impossible to realize coercivity with respect, for instance, to
the space ${\mathcal{C}^{2,\alpha}}$ where the involved objects are defined
and differentiation is undertaken (hence Assumption (A4)). In subsection 4.3,
we will show that the energy space of the bilinear form imposed by the shape
Hessian ${\operatorname{\mathnormal{d}}}^{2}J(\Omega)$ is the Sobolev space
$\bm{\mathbf{H}}^{1}(\Sigma)$, and that coercivity of the shape Hessian
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})$ can only be achieved in
the weaker space $\bm{\mathbf{H}}^{1/2}(\Sigma^{\ast})$. The latter result
implies that the shape problem under consideration is mildly ill-posed.

###  4.2 Shape Hessian of the cost function at the critical shape

Our goal here is to get the structure of the shape Hessian at a critical shape
(see Proposition 4.3). To do this we shall assume that $\Omega$ is of class
$\mathcal{C}^{2,1}$ and take the deformation fields from the set
$\bm{\mathbf{\Theta}}^{2}$. Given this assumption, the existence of the shape
derivative of the state is guaranteed. This implies that we only need to apply
the chain rule to get the desired expression. By this approach, however, we
first need to exhibit the shape derivative of the states which we give in the
next lemma. A revision of Proposition 4.3 with the mild $\mathcal{C}^{1,1}$
regularity assumption on the domain, however, is given in Proposition 4.4.

######  Lemma 4.2.

Let $\Omega\in\mathcal{C}^{2,1}$ and
$\bm{\mathbf{V}}\in\in\bm{\mathbf{\Theta}}^{2}$. Then,
$u\in\bm{\mathbf{\mathsf{H}}}^{3}(\Omega)$ is shape differentiable with
respect to $\Omega$ in the direction of $\bm{\mathbf{V}}$, and its shape
derivative $u^{\prime}\in\bm{\mathbf{\mathsf{H}}}^{1}(\Omega)$ uniquely
satisfies the boundary value problem

$-\Delta{u}^{\prime}=0\ \text{in $\Omega$},\quad{u}^{\prime}=0\ \text{on
$\Gamma$},\quad\partial_{\bm{\mathbf{n}}}{{u}^{\prime}}+i{u}^{\prime}=\Upsilon(u)[V_{n}]\
\text{on $\Sigma$},$ (27)

where
$\Upsilon(u)[V_{n}]={\operatorname{div}}_{\Sigma}(V_{n}\nabla_{\Sigma}{u})-i(\partial_{\bm{\mathbf{n}}}{{u}}+\kappa{u})V_{n}+\lambda\kappa
V_{n}$.

Because the above results is new and is in fact not yet available in the
literature, we provide the proof of the lemma in Appendix A.1. Now with the
lemma at hand, the first of the two main results of this subsection is now in
order.

######  Proposition 4.3.

Let $\Omega\in\mathcal{C}^{2,1}$ and
$\bm{\mathbf{V}},\bm{\mathbf{W}}\in\bm{\mathbf{\Theta}}^{2}$. Then, the shape
Hessian of $J$ at the solution $\Omega^{\ast}$ of the Bernoulli problem (2)
has the following structure:

${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]=-\int_{\Sigma^{\ast}}{\lambda(p_{1,W}^{\prime}+\kappa
p_{2,W}^{\prime}){V_{n}}}{\,{\operatorname{\mathnormal{d}}}\sigma},$ (28)

where $p^{\prime}_{W}=p_{1,W}^{\prime}+ip_{2,W}^{\prime}$ satisfies the
complex PDE system

$-\Delta{p}^{\prime}_{W}=u^{\prime}_{2,W}\ \text{in
$\Omega^{\ast}$},\qquad{p}^{\prime}_{W}=0\ \text{on
$\Gamma$},\qquad\partial_{\bm{\mathbf{n}}}{{p}^{\prime}_{W}}-i{p}^{\prime}_{W}=0\
\text{on $\Sigma^{\ast}$}.$ (29)

###### Proof.

The higher regularity assumption on $\Omega$ and on the deformation fields
allows us to apply formula (25) twice. Indeed, we may write
${\operatorname{\mathnormal{d}}}J(\Omega)[\bm{\mathbf{V}}]=\int_{\Omega}{u_{2}u_{2}^{\prime}}{\,{\operatorname{\mathnormal{d}}}x}+\int_{\Sigma}{\operatorname{div}\left(\frac{1}{2}|u_{2}|^{2}\bm{\mathbf{V}}\right)}{\,{\operatorname{\mathnormal{d}}}\sigma}$
by Stokes’ theorem. Then, by the same approach, we get

$\displaystyle{\operatorname{\mathnormal{d}}}^{2}J(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]$
$\displaystyle=\int_{\Omega}{u_{2}u^{\prime\prime}_{2,V,W}}{\,{\operatorname{\mathnormal{d}}}x}+\int_{\Omega}{u^{\prime}_{2,V}u^{\prime}_{2,W}}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle\qquad+\int_{\Omega}{\operatorname{div}\left[u_{2}u^{\prime}_{2,W}\bm{\mathbf{V}}+\operatorname{div}\left(\frac{1}{2}u_{2}^{2}\bm{\mathbf{V}}\right)\bm{\mathbf{W}}\right]}{\,{\operatorname{\mathnormal{d}}}x}.$

Clearly, at $\Omega=\Omega^{\ast}$, the above expression for the shape Hessian
reduces to
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]=\int_{\Omega^{\ast}}{u^{\prime}_{2,V}u^{\prime}_{2,W}}{\,{\operatorname{\mathnormal{d}}}x}$.
To obtain (28) from this equation, we utilize the adjoint problem (29). Note
that the shape derivative of the state and of the adjoint variable at the
stationary solution $\Omega^{\ast}$ of (2) are given by

$-\Delta{u}^{\prime}_{V}=0\ \text{in $\Omega^{\ast}$},\quad{u}^{\prime}_{V}=0\
\text{on
$\Gamma$},\quad\partial_{\bm{\mathbf{n}}}{{u}^{\prime}_{V}}+i{u}^{\prime}_{V}=\lambda(\kappa-i)V_{n}\
\text{on $\Sigma^{\ast}$},$ (30)

and by system (29), respectively. Multiplying each of these systems by
$\overline{p}^{\prime}_{W}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$ and
${u}^{\prime}_{V}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$, respectively, and then
applying integration by parts will eventually lead us to the identity
$\int_{\Omega^{\ast}}{u^{\prime}_{2,W}{u}^{\prime}_{V}}{\,{\operatorname{\mathnormal{d}}}x}=\int_{\Sigma}{\lambda\left(\kappa-i\right)\overline{p}^{\prime}_{W}V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}$.
Comparing the real and imaginary parts on both sides of this equation in the
end will give us (28). ∎

######  Remark 6.

Observing from problem (30), we see that $u^{\prime}_{2,V}\neq 0$ provided
$\bm{\mathbf{V}}\not\equiv\bm{\mathbf{0}}$ on $\Sigma^{\ast}$. Therefore, from
the proof of Proposition 4.3, we have
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{V}}]=\int_{\Omega^{\ast}}{|u^{\prime}_{2,V}|^{2}}{\,{\operatorname{\mathnormal{d}}}x}=\|u^{\prime}_{2,V}\|_{L^{2}(\Omega^{\ast})}^{2}>0$
for any $\bm{\mathbf{V}}\not\equiv\bm{\mathbf{0}}$ on $\Sigma^{\ast}$. This
inequality, however, does not mean that the shape optimization problem is
already well-posed. The matter will be discussed further in the next
subsection.

We remark that the second-order shape derivative
${\operatorname{\mathnormal{d}}}^{2}J(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]$
can actually be expressed using the material derivatives of $u$ and $p$
instead of their shape derivatives. This allows for a more general formula
since one does not need additional regularity for $u$ and $p$. Howbeit, the
order two analysis was only carried out in this work in order to examine the
ill-posedness of the proposed shape optimization problem (see also Remark 11).
Hence, we only need the expression for the shape Hessian at a critical shape
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]$,
and since we already assume that $\Omega$ is of class $\mathcal{C}^{2,1}$, the
second-order shape derivative
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]$
can be obtained with less effort and computations as the shape derivative for
$u$ and $p$ are already available, and that we do not need to find a
simplified form of
$\int_{\Omega}{u_{2}u^{\prime\prime}_{2,V,W}}{\,{\operatorname{\mathnormal{d}}}x}$
in terms of $u_{2}$ because the function already vanishes at
$\Omega=\Omega^{\ast}$. Still, we stress that it is possible to obtain the
second-order shape derivative
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]$
without using $u^{\prime\prime}_{2,V,W}$. This can be done by writing
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]$
using the material derivative $\dot{u}_{2}$ which then also allow one to
weaken the regularity assumption on the domain. In this case, we only need
$\Omega$ be of class $\mathcal{C}^{1,1}$. Indeed, from (18), we have
${\operatorname{\mathnormal{d}}}J(\Omega_{s})[\bm{\mathbf{V}}]=\frac{1}{2}\int_{\Omega}{({\operatorname{div}}\,\bm{\mathbf{V}})I_{s}|u_{2}^{s}|^{2}}{\,{\operatorname{\mathnormal{d}}}x}+\int_{\Omega}{u_{2}^{s}\dot{u}{}_{2}^{s}}{\,{\operatorname{\mathnormal{d}}}x}$.
Note here that the deformation field $\bm{\mathbf{V}}$ is independent of
$s>0$. Therefore, differentiating the integral expression with respect to $s$,
and then evaluating at $s=0$, we get

$\displaystyle{\operatorname{\mathnormal{d}^{2}}}J(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]$
$\displaystyle=\frac{1}{2}\int_{\Omega}{\left[({\operatorname{div}}\,\bm{\mathbf{V}})^{2}|u_{2}|^{2}+2({\operatorname{div}}\,\bm{\mathbf{V}})u_{2}\dot{u}_{2,W}\right]}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle\qquad+\int_{\Omega}{\left(\dot{u}_{2,V}\dot{u}_{2,W}+u_{2}\ddot{u}_{2,V,W}\right)}{\,{\operatorname{\mathnormal{d}}}x},$

for any given vector fields
$\bm{\mathbf{V}},\bm{\mathbf{W}}\in\bm{\mathbf{\Theta}}^{1}$. At
$\overline{\Omega}=\overline{\Omega}^{\ast}$, $u_{2}\equiv 0$, and so we
obtain
${\operatorname{\mathnormal{d}^{2}}}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]=\int_{\Omega}{\dot{u}_{2,V}\dot{u}_{2,W}}{\,{\operatorname{\mathnormal{d}}}x}$.
This integral can of course be expressed in terms of the material derivative
of $p$ using the same technique used in the proof of Theorem 3.5. In fact,
however, we can just simply introduce an adjoint problem in order to get rid
of the first-order derivative terms. To this end, let us consider the PDE
system

$-\Delta q_{W}=\dot{u}_{2,W}\ \text{in $\Omega$},\qquad q_{W}=0\ \text{on
$\Gamma$},\qquad\partial_{\bm{\mathbf{n}}}{}q_{W}-iq_{W}=0\ \text{on
$\Sigma$},$ (31)

at $\overline{\Omega}=\overline{\Omega}^{\ast}$ which can be shown to be well-
posed, and that $q_{W}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$ because of
Proposition 3.1. Here, of course, $q_{W}=q_{1,W}+iq_{2,W}$, where $q_{1,W}$
and $q_{2,W}$ denote the real and imaginary parts of $q_{W}$, respectively.
Multiplying the above equation by $\dot{u}_{V}$, and then applying integration
by parts, yields
$\int_{\Omega^{\ast}}{\nabla\bar{q}_{W}\cdot\nabla\dot{u}_{V}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma^{\ast}}{\bar{q}_{W}\dot{u}_{V}}{\,{\operatorname{\mathnormal{d}}}\sigma}=\int_{\Omega^{\ast}}{\dot{u}_{2,W}\dot{u}_{V}}{\,{\operatorname{\mathnormal{d}}}x}$.
Let us also consider equation (3.1) with
${v}=q_{W}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$,
$\overline{\Omega}=\overline{\Omega}^{\ast}$, and $\dot{u}$ replaced by
$\dot{u}_{V}$. Then,
$\int_{\Omega^{\ast}}{\nabla\dot{u}_{V}\cdot\nabla\bar{q}_{W}}{\,{\operatorname{\mathnormal{d}}}x}+i\int_{\Sigma^{\ast}}{\dot{u}_{V}\bar{q}_{W}}{\,{\operatorname{\mathnormal{d}}}\sigma}=-\int_{\Omega^{\ast}}{A\nabla{u_{1}}\cdot\nabla\bar{q}_{W}}{\,{\operatorname{\mathnormal{d}}}x}+\lambda\int_{\Sigma^{\ast}}{({\operatorname{div}}_{\Sigma}\bm{\mathbf{V}})\bar{q}_{W}}{\,{\operatorname{\mathnormal{d}}}\sigma}$.
Comparing this equation with the previous one, we get the equation
$\int_{\Omega^{\ast}}{\dot{u}_{2,W}\dot{u}_{V}}{\,{\operatorname{\mathnormal{d}}}x}=-\int_{\Omega^{\ast}}{A\nabla{u_{1}}\cdot\nabla\bar{q}_{W}}{\,{\operatorname{\mathnormal{d}}}x}+\lambda\int_{\Sigma^{\ast}}{({\operatorname{div}}_{\Sigma}\bm{\mathbf{V}})\bar{q}_{W}}{\,{\operatorname{\mathnormal{d}}}\sigma}$.
Meanwhile, it can be verified that a similar identity to (16) holds for
$\dot{u}_{2}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$ and
$\bar{q}_{W}\in\bm{\mathbf{\mathsf{V}}}(\Omega)$ at
$\overline{\Omega}=\overline{\Omega}^{\ast}$ that is given by
$-\int_{\Omega^{\ast}}{A\nabla{u_{1}}\cdot\nabla\bar{q}_{W}}{\,{\operatorname{\mathnormal{d}}}x}=-\int_{\Sigma^{\ast}}{(\nabla\bar{q}_{W}\cdot\nabla{u_{1}})V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}-i\int_{\Sigma^{\ast}}{\bar{q}_{W}(\bm{\mathbf{V}}\cdot\nabla{u_{1}})}{\,{\operatorname{\mathnormal{d}}}\sigma}+\lambda\int_{\Sigma^{\ast}}{(\bm{\mathbf{V}}\cdot\nabla\bar{q}_{W})}{\,{\operatorname{\mathnormal{d}}}\sigma}$.
Putting together the last two equations, and then appealing to formula (22),
we get

$\displaystyle\int_{\Omega^{\ast}}{\dot{u}_{2,W}\dot{u}_{V}}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle=-i\int_{\Sigma^{\ast}}{\lambda\bar{q}_{W}V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}+\lambda\int_{\Sigma^{\ast}}{\kappa\bar{q}_{W}V_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$

Finally, comparing the real and imaginary parts on both sides of the above
equation, we arrive at the final expression for the shape Hessian of $J$ at a
critical shape $\Omega^{\ast}$:

${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]=-\int_{\Sigma^{\ast}}{\lambda(q_{1,W}+\kappa
q_{2,W}){V_{n}}}{\,{\operatorname{\mathnormal{d}}}\sigma}.$

The above integral expression is clearly identical to (28) in terms of
structure.

To close this subsection, we formally summarize the previous result into a
proposition given below which is essentially a revision of Proposition 4.3
with weaker regularity assumptions on the domain $\Omega$ and on the
deformation fields $\bm{\mathbf{V}}$ and $\bm{\mathbf{W}}$, and combined with
a different adjoint problem.

######  Proposition 4.4.

Let $\Omega\in\mathcal{C}^{1,1}$ and
$\bm{\mathbf{V}},\bm{\mathbf{W}}\in\bm{\mathbf{\Theta}}^{1}$. Then, the shape
Hessian of $J$ at the solution $\Omega^{\ast}$ of the Bernoulli problem (2)
has the following structure:

${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]=-\int_{\Sigma^{\ast}}{\lambda(q_{1,W}+\kappa
q_{2,W}){V_{n}}}{\,{\operatorname{\mathnormal{d}}}\sigma},$ (32)

where $q_{W}=q_{1,W}+iq_{2,W}$ satisfies the complex PDE system

$-\Delta{q}_{W}=\dot{u}_{2,W}\ \text{in $\Omega^{\ast}$},\qquad{q}_{W}=0\
\text{on $\Gamma$},\qquad\partial_{\bm{\mathbf{n}}}{{q}_{W}}-i{q}_{W}=0\
\text{on $\Sigma^{\ast}$}.$ (33)

######  Remark 7.

As was shown in the computation of expression (32) issued above, the shape
derivative of the state is not needed to exhibit the structure of the shape
Hessian at a critical shape. It therefore goes without saying that the
expression for the shape Hessian of the cost function, for general domains of
class $\mathcal{C}^{1,1}$, can also be obtained without using the second-order
shape derivative of the state by using the same technique applied to show
(32). Here, we omit the computation of the said expression since we are only
interested in the structure of
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})$ (see Remark 11).

###  4.3 Compactness of the Hessian at the optimal domain

Having computed the form of the shape Hessian at a stationary domain
$\Omega^{\ast}$, we now examine the stability or instability of the presently
proposed shape optimization formulation. Before that, we first make a few
remarks regarding the structure of the shape Hessian.

It can be verified that the exact form of the shape Hessian
${\operatorname{\mathnormal{d}}}^{2}J(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]$
depends on the shape derivative $\kappa^{\prime}$ of the mean curvature
$\kappa$. In fact, with the first-order shape derivative of $J$ given as
${\operatorname{\mathnormal{d}}}J(\Omega)[\bm{\mathbf{V}}]=\int_{\Sigma}{G\bm{\mathbf{n}}\cdot\bm{\mathbf{V}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$,
where $G$ is given by (13), the shape Hessian can be shown to have the
structure

$\displaystyle{\operatorname{\mathnormal{d}}}^{2}J(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]=\int_{\Sigma}{\left[{G}^{\prime}_{W}V_{n}+\left(\partial_{\bm{\mathbf{n}}}{G}+\kappa{G}\right)V_{n}W_{\text{n}}-{G}K+{G}(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{V}})W_{\text{n}}\right]}{\,{\operatorname{\mathnormal{d}}}\sigma},$

where
$K=\bm{\mathbf{v}}_{\Sigma}\cdot(\operatorname{\mathnormal{D}}\\!_{\Sigma}\bm{\mathbf{n}})\bm{\mathbf{w}}_{\Sigma}+\bm{\mathbf{n}}\cdot(\operatorname{\mathnormal{D}}\\!_{\Sigma}\bm{\mathbf{v}})\bm{\mathbf{w}}_{\Sigma}+\bm{\mathbf{n}}\cdot(\operatorname{\mathnormal{D}}\\!_{\Sigma}\bm{\mathbf{w}})\bm{\mathbf{v}}_{\Sigma}$,
$\bm{\mathbf{v}}=\bm{\mathbf{V}}|_{\Sigma}$,
$\bm{\mathbf{v}}=\bm{\mathbf{v}}_{\Sigma}+v_{\text{n}}\bm{\mathbf{n}}:=(\bm{\mathbf{v}}\cdot\tau)\tau+(\bm{\mathbf{v}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}$
and $\operatorname{\mathnormal{D}}\\!_{\Sigma}$ denotes the tangential
differential operator called the _tangential Jacobian matrix_ given as
$\operatorname{\mathnormal{D}}\\!_{\Sigma}\bm{\mathbf{v}}=\operatorname{\mathnormal{D}}\\!\bm{\mathbf{V}}|_{\Sigma}-(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{V}}\bm{\mathbf{n}})\bm{\mathbf{n}}^{\top}$
(see, e.g., [17, Eq. (5.2), p. 495]). Obviously, from (13), the term
${G}^{\prime}$ appearing in the shape Hessian is composed of the derivative
$\bm{\mathbf{n}}^{\prime}$ of the normal vector and and $\kappa^{\prime}$ of
the mean-curvature. The form of $\bm{\mathbf{n}}^{\prime}_{W}$ and
$\kappa^{\prime}_{W}$ obtained along the deformation field
$\bm{\mathbf{W}}\in\bm{\mathbf{\Theta}}^{2}$ are respectively given by
$(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{W}}\bm{\mathbf{n}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}-(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{W}})^{\top}\bm{\mathbf{n}}-(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{n}})\bm{\mathbf{W}}$
and
$\text{trace}\left\\{\operatorname{\mathnormal{D}}\\!\left[(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{W}}\bm{\mathbf{n}}\cdot\bm{\mathbf{n}})\bm{\mathbf{n}}-(\operatorname{\mathnormal{D}}\\!\bm{\mathbf{W}})^{\top}\bm{\mathbf{n}}\right]-\operatorname{\mathnormal{D}}\\!\bm{\mathbf{n}}\operatorname{\mathnormal{D}}\\!\bm{\mathbf{W}}\right\\}-\nabla\kappa\cdot\bm{\mathbf{W}}$
(see [17, 59]). Clearly, the latter expression consists of a second-order
tangential derivative of the velocity field $\bm{\mathbf{W}}$, and this
derivative actually exists due to our assumption that $\Omega$ is of class
$\mathcal{C}^{2,1}$ [17, 59]. Hence, we can actually decompose the shape
Hessian into three parts: $h_{1}(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]$,
$h_{2}(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]$, and
$h_{3}(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]$, where $h_{2}$ is composed of
the term $\bm{\mathbf{n}}^{\prime}_{W}$ while $h_{3}$ consists of the
expression $\kappa^{\prime}_{W}$, and such that we have the estimates
$|h_{1}(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]|\lesssim\|\bm{\mathbf{V}}\|_{\bm{\mathbf{L}}^{2}(\Sigma)}\|\bm{\mathbf{W}}\|_{\bm{\mathbf{L}}^{2}(\Sigma)}$,
$|h_{2}(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]|\lesssim\|\bm{\mathbf{V}}\|_{\bm{\mathbf{H}}^{1/2}(\Sigma)}\|\bm{\mathbf{W}}\|_{\bm{\mathbf{H}}^{1/2}(\Sigma)}$,
and
$|h_{3}(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]|\lesssim\|\bm{\mathbf{V}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)}\|\bm{\mathbf{W}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)}$
(cf. [26, proof of Thm. 3]). From these, we can infer that the shape Hessian
defines a continuous bilinear form
${\operatorname{\mathnormal{d}}}^{2}J(\Omega):\bm{\mathbf{H}}^{1}(\Sigma)\times\bm{\mathbf{H}}^{1}(\Sigma)\to\mathbb{R}$;
that is,
$|{\operatorname{\mathnormal{d}}}^{2}J(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]|\lesssim\|\bm{\mathbf{V}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)}\|\bm{\mathbf{W}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)}$.555The
notation $\bm{\mathbf{H}}^{1}(\cdot)$ stands for the Sobolev space
$\bm{\mathbf{H}}^{1}(\cdot):=\\{\bm{\mathbf{v}}:=(v_{1},v_{2}):v_{1},v_{2}\in
H^{1}(\cdot)\\}$ and is equipped with the norm
$\|\bm{\mathbf{v}}\|^{2}_{\bm{\mathbf{H}}^{1}(\cdot)}=\|v_{1}\|^{2}_{H^{1}(\cdot)}+\|v_{2}\|^{2}_{H^{1}(\cdot)}$.
Similar definition is also given to the
$\bm{\mathbf{H}}^{1}_{\Gamma,0}(\cdot)$-space. In addition, one can also
obtained a remainder estimate given by
$|{\operatorname{\mathnormal{d}}}^{2}J(\Omega)[\bm{\mathbf{V}},\bm{\mathbf{W}}]-{\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{W}}]|\leqslant\eta\left(\mathcal{M}(\Omega,\Omega^{\ast})\right)\|\bm{\mathbf{V}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)}\|\bm{\mathbf{W}}\|_{\bm{\mathbf{H}}^{1}(\Sigma)}$,
where $\mathcal{M}(B_{1},B_{2})$ is some appropriate metric measuring the
distance between the two sets $B_{1}$ and $B_{2}$ in $\mathbb{R}^{d}$ while
$\eta:\mathbb{R}^{+}_{0}\to\mathbb{R}^{+}_{0}$ is a decreasing function that
satisfies $\eta(s)\to 0$ as $s\to 0$ (cf. Assumption A3).

In connection with the above discussion, it is natural to ask whether we also
have the estimate
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{V}}]\gtrsim\|\bm{\mathbf{V}}\|^{2}_{\bm{\mathbf{H}}^{1}(\Sigma^{\ast})}$.
This inequality condition actually has something to do with the stability of a
local minimizer $\Omega^{\ast}$ of $J$. A result regarding sufficient second-
order conditions from [10, 11, 20]666Independently of [10, 11], also Eppler
derived second-order sufficient optimality conditions in [20]; see
particularly Section 4 of the said paper. in fact states that a local
minimizer $\Omega^{\ast}$ is stable if and only if the shape Hessian
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})$ is strictly coercive in
its corresponding energy space which is the
$\bm{\mathbf{H}}^{1}(\Sigma^{\ast})$ space in the present case. Such strict
coercivity, however, cannot be established for the shape Hessian (28). The
aforesaid lack of coercivity is known, especially in the shape optimization
literature, as the two-norm discrepancy, see [10, 19, 27] for more details.
Meanwhile, for a more recent study concerning the question of stability in the
field of shape optimization – focusing especially on the strategy using
second-order shape derivatives – we refer the readers to [12].

To analyze the shape Hessian (28), we will write it into an equivalent
expression (see [21, 22, 24, 25, 26]) and adapt the method already used in [2,
3]. We first introduce the operators $\mathcal{L}:\
\bm{\mathbf{H}}^{1/2}(\Sigma^{\ast})\to H^{1/2}(\Sigma^{\ast})$ and
$\mathcal{K}:\ H^{1/2}(\Sigma^{\ast})\to H^{1/2}(\Sigma^{\ast})$ (see [52,
Sec. 3.4]), defined respectively as $\mathcal{L}\bm{\mathbf{V}}:=\lambda
V_{n}$ and $\mathcal{K}v:=\kappa v$, which we shall utilize in our
argumentation. The continuity of the operators $\mathcal{L}$ and $\mathcal{K}$
follow from the next lemma (cf. [21, Lem. 3.3] and see also [24, 25, 26]).

######  Lemma 4.5.

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded Lipschitz domain with boundary
$\Gamma:=\partial\Omega$. Then, the map $v\mapsto\phi v$ is continuous in
$H^{1/2}(\Gamma)$ for any $v\in H^{1/2}(\Gamma)$ and
$\phi\in\mathcal{C}^{0,1}(\Gamma)$.

###### Proof.

By McShane-Whitney extension theorem, there is some function
$\tilde{\phi}\in\mathcal{C}^{0,1}(\overline{\Omega})$ such that
$\tilde{\phi}|_{\Gamma}=\phi$. Also, from [45, Thm. 3.37, p. 102], there is a
bounded linear extension operator $\tilde{E}:H^{1/2}(\Gamma)\to
H^{1}(\Omega)$. Hence, in view of [33, Thm. 1.4.1.1, p. 21] and by trace
theorem, we see that the operator $\mathcal{K}v=\phi v$ is continuous in
$H^{1/2}(\Gamma)$ since
$\mathcal{K}v={\operatorname{trace}}(\tilde{\phi}\tilde{E}v)$ is a composition
of bounded operators. ∎

######  Remark 8.

Let us note that for $\mathcal{C}^{1,1}$ smooth boundary $\Sigma$, $\kappa$ is
well defined almost everywhere on the boundary and actually belongs to
$L^{\infty}$ on account of Rademacher’s theorem (see, e.g., [56, Thm. 2.7.1,
p. 67]). Moreover, the mean curvature is continuous for boundaries of class
$\mathcal{C}^{2,1}$. Due to our smoothness assumption on the free boundary,
the operator $\mathcal{K}$ is actually a continuous map from
$H^{s}(\Sigma^{\ast})$ to $H^{s}(\Sigma^{\ast})$ for all $s\in[0,1]$, and the
same is also true for the (bijective) map $\mathcal{L}$.

In addition to $\mathcal{L}$ and $\mathcal{K}$, we also introduce the operator
$\mathcal{P}_{i}$, $i\in\\{1,2\\}$, by
$\mathcal{P}_{i}:\bm{\mathbf{H}}^{1/2}(\Sigma^{\ast})\to
H^{-1/2}(\Sigma^{\ast})$, $\bm{\mathbf{V}}\mapsto p_{i,V}^{\prime}$, for all
$i\in\\{1,2\\}$. Accordingly, we may write the shape Hessian at the optimal
domain $\Omega^{\ast}$ as

$\displaystyle{\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{V}}]$
$\displaystyle=-\int_{\Sigma^{\ast}}{(\lambda{V_{n}}p_{1,V}^{\prime}+\kappa\lambda{V_{n}}p_{2,V}^{\prime})}{\,{\operatorname{\mathnormal{d}}}\sigma}$
$\displaystyle=-\langle\mathcal{L}\bm{\mathbf{V}},\mathcal{P}_{1}\bm{\mathbf{V}}\rangle-\langle\mathcal{K}\mathcal{L}\bm{\mathbf{V}},\mathcal{P}_{2}\bm{\mathbf{V}}\rangle,$

where $\langle\,\cdot\,,\,\cdot\,\rangle$ is the duality product between
$H^{1/2}(\Sigma^{\ast})$ and $H^{-1/2}(\Sigma^{\ast})$.

Now, in relation to Remark 8, we state our final result which claims the ill-
posedness of the proposed CCBM formulation of (2).

######  Proposition 4.6.

Let $\Omega^{\ast}$ be the stationary solution to (2), then the Riesz operator
associated to
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{V}}]:\bm{\mathbf{H}}^{1/2}(\Sigma^{\ast})\to
H^{-1/2}(\Sigma^{\ast})$ is compact.

###### Proof.

The idea of the proof is to express the shape Hessian as a composition of
linear continuous operators and a compact one (the compactness being obtained
using the compactness of the imbedding between two Sobolev spaces). As shown
above, the operators $\mathcal{L}$ and $\mathcal{K}$ are continuous, but the
operators $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ are compact. We verify this
claim by decomposing these operators as composition of continuous and compact
ones. For the product
$-\langle\mathcal{L}\bm{\mathbf{V}},\mathcal{P}_{1}\bm{\mathbf{V}}\rangle$, we
decompose the map $\mathcal{P}_{1}$ by first expressing it as the composition
$\mathcal{P}_{1}=\mathcal{Q}_{2}\circ\mathcal{Q}_{1}$ where
$\mathcal{Q}_{1}:\bm{\mathbf{H}}^{1/2}(\Sigma^{\ast})\to
H^{1}(\Omega^{\ast})$, $\bm{\mathbf{V}}\mapsto u_{2}^{\prime}$, and
$\mathcal{Q}_{2}:H^{1}(\Omega^{\ast})\to H^{-1/2}(\Sigma^{\ast})$,
$\phi\mapsto{w}$. Here, $u_{2}^{\prime}$ solves (27) and
${w}={w}_{1}+i{w}_{2}$ satisfies

$-\Delta{w}=\phi\ \text{in $\Omega^{\ast}$},\qquad{w}=0\ \text{on
$\Gamma$},\qquad\partial_{\bm{\mathbf{n}}}{{w}}-i{w}=0\ \text{on
$\Sigma^{\ast}$}.$ (34)

Clearly, $\mathcal{Q}_{1}$ is continuous. Next, we further write
$\mathcal{Q}_{2}$ as
$\mathcal{Q}_{2}:=\mathcal{R}_{3}\circ\mathcal{R}_{2}\circ\mathcal{R}_{1}$
where

$\displaystyle\mathcal{R}_{1}:H^{1}(\Omega^{\ast})\to
H^{3}(\Omega^{\ast}),\quad\phi\mapsto{w}_{1},$
$\displaystyle\qquad\mathcal{R}_{2}:H^{3}(\Omega^{\ast})\to
H^{5/2}(\Sigma^{\ast}),\quad{w}_{1}\mapsto{w}_{1},$
$\displaystyle\qquad\text{and}\qquad\mathcal{R}_{3}:H^{5/2}(\Sigma^{\ast})$
$\displaystyle\to H^{-1/2}(\Sigma^{\ast}),\quad{w}_{1}\mapsto{w}_{1}.$

The operators $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ are continuous while
$\mathcal{R}_{3}$ is the compact embedding of $H^{5/2}(\Sigma^{\ast})$ into
$H^{-1/2}(\Sigma^{\ast})$.777An embedding result for fractional Sobolev spaces
$H^{s}(\Omega)$ under bounded $\mathcal{C}^{k,\alpha}$-domains,
$k=0,1,\ldots$, $\alpha\in[0,1]$, can be found in [63, Thm. 7.9, p. 119]. See
also [56, Thm. 2.5.5, p. 61], but for Lipschitz and $\mathcal{C}^{k}$ domains.

For the product
$-\langle\mathcal{K}\mathcal{L}\bm{\mathbf{V}},\mathcal{P}_{2}\bm{\mathbf{V}}\rangle$,
a similar decomposition is applied, but with an additional decomposition to
get the compactness result. Indeed, let us decompose the map $\mathcal{P}_{2}$
as the composition $\mathcal{P}_{2}=\mathcal{Q}_{2}\circ\mathcal{Q}_{1}$ where
$\mathcal{Q}_{1}:\bm{\mathbf{H}}^{1/2}(\Sigma^{\ast})\to
H^{1}(\Omega^{\ast})$, $\bm{\mathbf{V}}\mapsto u_{2}^{\prime}$, and
$\mathcal{Q}_{2}:H^{1}(\Omega^{\ast})\to H^{-1/2}(\Sigma^{\ast})$,
$\phi\mapsto{w}$. Here, again, $u_{2}^{\prime}$ and ${w}$ respectively solves
(27) and (34). Next, we write $\mathcal{Q}_{2}$ as
$\mathcal{Q}_{2}:=\mathcal{S}_{4}\circ\mathcal{S}_{3}\circ\mathcal{S}_{2}\circ\mathcal{S}_{1}$
where

$\displaystyle\mathcal{S}_{1}:H^{1}(\Omega^{\ast})\to H^{3}(\Omega^{\ast}),\ \
\phi\mapsto{w}_{1};\qquad\mathcal{S}_{2}:H^{3}(\Omega^{\ast})\to
H^{3/2}(\Sigma^{\ast}),\ \ {w}_{1}\mapsto\partial_{\bm{\mathbf{n}}}{{w}_{1}};$
$\displaystyle\mathcal{S}_{3}:H^{3/2}(\Sigma^{\ast})\to
H^{5/2}(\Sigma^{\ast}),\ \
\partial_{\bm{\mathbf{n}}}{{w}_{1}}\mapsto{w}_{2};\qquad\mathcal{S}_{4}:H^{5/2}(\Sigma^{\ast})\to
H^{-1/2}(\Sigma^{\ast}),\ \ {w}_{2}\mapsto{w}_{2}.$

The first three operators $\mathcal{S}_{1}$, $\mathcal{S}_{2}$, and
$\mathcal{S}_{3}$ are continuous while $\mathcal{S}_{4}$ is the compact
embedding of $H^{5/2}(\Sigma^{\ast})$ into $H^{-1/2}(\Sigma^{\ast})$. This
proves the compactness result. ∎

Related results, specifically given for $L^{2}$-tracking type and compact
gradient tracking functionals, are issued in [21, Prop. 3.1] and [23, Prop.
2.10]. It is worth to remark that the compactness of the shape Hessian stated
in Proposition 4.6 opposes the strict coercivity
${\operatorname{\mathnormal{d}}}^{2}J(\Omega^{\ast})[\bm{\mathbf{V}},\bm{\mathbf{V}}]\gtrsim\|\bm{\mathbf{V}}\|^{2}_{\bm{\mathbf{H}}^{1}(\Sigma^{\ast})}$
that corresponds to the well-posedness of the proposed shape optimization
formulation of the free boundary problem in consideration (on a related note,
see [21, Rem. 3.2–3.3] and [23, Rem. 2.11]).

##  5 Numerical approximation

The numerical resolution to our proposed shape optimization approach to (2) is
carried out using a Sobolev gradient-based method. The implementation is
realized in line with the author’s previous work using the finite element
method, see [51, 52, 53], but with some changes which are crucial for
assessing the numerical performance of the new method over the classical Kohn-
Vogelius cost functional approach. The proposed shape optimization
reformulation can of course be solved numerically using other methods such as
the level-set method (see [50]) – an Eulerian-like type numerical scheme –
employed, for instance, in [1, 35, 42], or via a boundary element method
through the concept of boundary integral equations used in [22, 24, 25, 26,
34].

###  5.1 Numerical algorithm

For completeness, we give below the important details of our algorithm.

Choice of descent direction. The choice
$\bm{\mathbf{V}}=V_{n}\bm{\mathbf{n}}=-G\bm{\mathbf{n}}$, $G\in
L^{2}(\Sigma)$, $G\not\equiv 0$, provides a descent direction for the cost
function $J(\Omega)$ given in (10). Indeed, in general, the inequality
condition
$J(\Omega_{t})=J(\Omega)+t\left.\frac{\operatorname{\mathnormal{d}}}{\operatorname{\mathnormal{d}}\varepsilon}J(\Omega_{\varepsilon})\right|_{\varepsilon=0}+O(t^{2})=J(\Omega)+t\int_{\Sigma}{GV_{n}}{\,{\operatorname{\mathnormal{d}}}\sigma}+O(t^{2})=J(\Omega)-t\int_{\Sigma}{|G|^{2}}{\,{\operatorname{\mathnormal{d}}}\sigma}+O(t^{2})<J(\Omega)$,
holds for sufficiently small real number $t>0$. However, as alluded above, we
make use of the Riesz representation of the shape gradient. More exactly, we
apply an extension-regularization technique by taking the descent direction
$\bm{\mathbf{V}}$ as the solution in
$\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega)$ to the variational
problem
$a(\bm{\mathbf{V}},\bm{\mathbf{\varphi}})=-\int_{\Sigma}{G\bm{\mathbf{n}}\cdot\bm{\mathbf{\varphi}}}{\,{\operatorname{\mathnormal{d}}}\sigma}$,
for all
$\bm{\mathbf{\varphi}}\in\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega)$,
where $a$ is the ${\bm{\mathbf{H}}^{1}(\Omega)}(:=H^{1}(\Omega)^{d}$)-inner
product in $d$-dimension, $d\in\\{2,3\\}$ (cf. eq. (36) in next subsection).
In this sense, the Sobolev gradient $\bm{\mathbf{V}}$ [47] becomes a smoothed
preconditioned extension of $-G\bm{\mathbf{n}}$ over the entire domain
$\Omega$. For more discussion about discrete gradient flows for shape
optimization, we refer the readers to [18].

The main algorithm. The main steps of the iterative algorithm, computing the
$k$th domain $\Omega^{k}$, is summarized as follows:

1\. Initilization

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

2\. Iteration

For $k=0,1,2,\ldots$

  1. 2.1

Solve the state and adjoint state systems on the current domain $\Omega^{k}$.

  2. 2.2

Choose $t^{k}>0$, and compute the Sobolev gradient $\bm{\mathbf{V}}^{k}$ in
$\Omega^{k}$.

  3. 2.3

Update the current domain by setting
$\Omega^{k+1}=(\operatorname{id}+t^{k}\bm{\mathbf{V}}^{k})\Omega^{k}$.

3\. Stop Test

Repeat the Iteration until convergence.

######  Remark 9 (Step-size computation).

The step size $t^{k}$ is computed via a backtracking line search procedure
using the formula $t^{k}=\mu
J(\Omega^{k})/|\bm{\mathbf{V}}^{k}|^{2}_{\mathbf{H}^{1}(\Omega^{k})}$ at each
iteration, where $\mu>0$ is a given real number. This choice of the step size
is based on an Armijo-Goldstein-like condition for the shape optimization
method, see, for example, [53, p. 281]. In our application, however, the step
size parameter $\mu$ is not limited to the interval $(0,1)$ in contrast to
[51, 52, 53]. The less restrictive choice of the value for $\mu$ allows us to
numerically evaluate the sensitivity of the cost functions $J(\Omega)$ and
$J_{KV}(\Omega)$ through large deformations of the domain $\Omega$.

######  Remark 10 (Stopping conditions).

The algorithm is stopped as soon as $\Omega^{k}$, $\Sigma^{k}$, and
$\bm{\mathbf{V}}^{k}$ satisfy the inequality condition

$\max\left(\sqrt{a(\bm{\mathbf{V}}^{k},\bm{\mathbf{V}}^{k})},\|\bm{\mathbf{V}}^{k}\|_{\mathcal{C}(\Sigma^{k})^{d}},J(\Omega^{k})\right)<\texttt{Tol},$

for some fixed small value $\texttt{Tol}>0$. We also terminate the algorithm
as soon as the absolute difference $|J(\Omega^{k})-J(\Omega^{k-1})|$ is small
enough, or after a finite number of iterations.

######  Remark 11.

The convergence behavior of a gradient-based iterative scheme can be improved
by incorporating the Hessian information in the numerical procedure. The
drawback, however, of a second-order method is that, typically, it demands
additional computational burden and time to carry out the calculation,
especially when the Hessian is complicated [48, 58]. Here, we will not employ
a second-order method to numerically solve the optimization problem. The order
two analysis was performed here only to carry out the stability analysis for
the proposed optimization problem.

###  5.2 The extension-regularization technique preserves the critical shape

In this intermediate subsection, we issue a small result concerning the
stationary point $\Omega^{\ast}$ of the evolving boundary $\Omega(t)$, where
$t\geqslant 0$ (interpreted here as a “pseudo-time” step), that evolves from
an initial geometric profile $\Omega(0)$ under the pseudo flow field
$V_{n}\bm{\mathbf{n}}$. This flow field is related to the extension-
regularization technique used to compute the descent field at the beginning of
the previous section. Here, there will be a slight abuse of notation. In
previous discussions $\bm{\mathbf{V}}$ stands for the deformation field that
deforms the reference domain $\Omega$ using the application of the operator
$T_{t}$. In the arguments given below, we will be using the same notation to
represent an extended-regularized (normal) flow field for the evolving
boundary $\Sigma(t)$ which, in some sense, related to the extension-
regularization technique presented in 5.1 for the computation of the
associated descent direction (i.e., in other words, we somehow view the
evolution of the free boundary $\Sigma$ generated by the approximation process
as an evolving boundary problem). We emphasize that the discussion given below
does not attempt to prove the existence and/or convergence of approximate
shape solutions concerning the present shape optimization problem. Indeed, a
careful analysis along the lines of arguments used in [27, Sec. 3] (see also
[38]) is needed to address such delicate issue and actually goes beyond the
scope of the present study. Besides, some key assumptions on the admissible
set and on the functional $J$ have to be imposed (see [27, Sec. 3]). Moreover,
as in many optimization problems, it is important to define the notion of
convergence of $\Omega(t)$ to the stationary point $\Omega^{\ast}$ if one
wants to view the sequence $\\{\Omega(t_{k})\\}$, $k=0,1,\ldots$, as a
sequence of approximations of $\Omega^{\ast}$ in the discrete setting; see
[38, Chap. 2]. For a closely related topic concerning a finite element
approximation for shape optimization problems with mixed boundary conditions,
we refer the readers to [61]. In the sequel, the phrases ‘moving boundary’,
‘evolving boundary’, and ‘free boundary’ are used interchangeably.

To proceed, let us consider the following abstract autonomous evolving
boundary problem.

######  Problem 2.

Let $\Omega$ be a bounded annular domain with $\mathcal{C}^{\circ}$ regular
boundaries $\Gamma$ and $\Sigma$ such that $\Sigma$ is exterior to the fixed
(non-moving) boundary $\Gamma$. Given an initial geometry $\Sigma^{0}$ that is
$\mathcal{C}^{\circ}$ regular and a real-valued function $\Phi$ defined on
$\Sigma$ (i.e., $\Phi(\,\cdot\,;\Sigma):\Sigma\to\mathbb{R}$), find a moving
boundary/surface $\Sigma(t)$, $t\geqslant 0$, with the normal speed $V_{n}$
which satisfies

$V_{n}(x,t)=\Phi(x;\Sigma(t)),\quad x\in\Sigma(t),\quad t\geqslant
0,\qquad\Sigma(0)=\Sigma^{0}.$ (35)

In (35), to keep the regularity of the initial domain – which we assume to be
at least $\mathcal{C}^{\circ}$ regular – during evolution, one needs the
function $\Phi(x;\Sigma(t))$ to also be at least $\mathcal{C}^{\circ}$ for all
$x\in\Sigma(t)$, for $t>0$ (unless specified, this will be assumed in the rest
of the discussion). Given this assumption on $\Phi$, for sufficiently small
$\varepsilon>0$, it can be shown that the moving boundary $\Sigma(t)$ actually
maintains the $\mathcal{C}^{\circ}$ regularity throughout the short time
interval $[0,\varepsilon)$.

We next define a stationary solution to Problem 2 as follows.

######  Definition 1.

A domain $\Omega^{\ast}$ is said to be a stationary solution to Problem 2 if
$\Sigma^{\ast}=\partial\Omega^{\ast}\setminus\Gamma$, and
$\Phi(x;\Sigma^{\ast})=0$, for all $x\in\Sigma^{\ast}$.

Also, in accordance with the discussion issued in the previous subsection
about the extension-regularization technique used to compute for the
perturbation field, we associate with Problem 2 the extended-regularized
evolving boundary problem stated as follows.

######  Problem 3.

Given an annular domain $\Omega^{0}\in\mathcal{C}^{\circ}$ with boundary
$\Gamma\cup\Sigma^{0}$ ($\Sigma^{0}$ is exterior to $\Gamma$) and a real-
valued function $\Phi(\,\cdot\,;\Sigma(t))\in
L^{2}(\Sigma(t))\cap\mathcal{C}^{1,\alpha}(\Sigma(t))$, $\alpha>0$, $t>0$, we
seek to find an evolving boundary/surface $\Sigma(t)$ with initial profile
$\Sigma(0)=\Sigma^{0}$ that solves the problem

$\left\\{\begin{array}[]{rcll}-\Delta\bm{\mathbf{V}}+\bm{\mathbf{V}}&=&\bm{\mathbf{0}}&\text{in
$\Omega(t)$},\\\ {\bm{\mathbf{V}}}&=&\bm{\mathbf{0}}&\text{on $\Gamma$},\\\
\nabla\bm{\mathbf{V}}\cdot\bm{\mathbf{n}}&=&\Phi(\,\cdot\,;\Sigma(t)){\bm{\mathbf{n}}}&\text{on
$\Sigma(t)$},\\\ V_{n}&=&\bm{\mathbf{V}}\cdot\bm{\mathbf{n}}&\text{on
$\Sigma(t)$,\quad$t\geqslant 0$},\\\
\Omega(0)&=&\Omega^{0}.\end{array}\right.$ (36)

For a domain $\Omega$ that is $\mathcal{C}^{2,\alpha}$ regular, $\alpha>0$,
the (outward unit) normal vector $\bm{\mathbf{n}}$ is
$\mathcal{C}^{1,\alpha}(\partial\Omega)$ smooth. Hence, for fixed $t\geqslant
0$ and $\Omega(t)\in\mathcal{C}^{2,\alpha}$, it can be proved (using, for
example, the results from [31, 43]) that the first three equations in (36)
admits a unique (classical) solution
$\bm{\mathbf{V}}\in\bm{\mathbf{\mathcal{C}}}^{2,\alpha}(\Omega(t))$, for any
given $\Phi(\,\cdot\,;\Sigma(t))\in\mathcal{C}^{1,\alpha}(\Sigma(t))$,
$\alpha>0$. Notice here that, for fixed $t$, $\Phi(\,\cdot\,;\Sigma(t))$ only
needs to be $\mathcal{C}^{1,\alpha}(\Sigma(t))$ regular and
$\Omega(t)\in\mathcal{C}^{2,\alpha}$ for $\bm{\mathbf{V}}$ to be
$\bm{\mathbf{\mathcal{C}}}^{2,\alpha}(\Omega(t))$ smooth. We also remark that
the evolution of $\Omega(t)$ is essentially assumed here as perturbations of
$\Omega^{0}$ that can also be obtained via a diffeomorphic map which is close
to the identity.

We called equation (36) an extended-regularized evolving boundary problem
because, originally in Problem 2, $V_{n}$ is only supported on $\Sigma(t)$,
and, in (36), we want not only to extend the vector $V_{n}\bm{\mathbf{n}}$ in
the entirety of $\Omega(t)$, for $t>0$ via equation (36)1, but also to add
more regularity on the normal flow field through equation (36)3. Meanwhile, we
have equation (36)2 on $\Gamma$ since we want the interior boundary to remain
fixed during evolution.

Finally, in relation to Problem 3, a stationary solution $\Omega^{\ast}$ is
define next.

######  Definition 2.

A domain $\Omega^{\ast}\in\mathcal{C}^{\circ}$ is said to be a stationary
solution to Problem 3 if $\Sigma^{\ast}=\partial\Omega^{\ast}\setminus\Gamma$
and
${\bm{\mathbf{V}}}\in\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega^{\ast})\cap\bm{\mathbf{\mathcal{C}}}^{\circ}(\Omega^{\ast})$888Here
we define
$\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega^{\ast}):=H_{\Gamma,\bm{\mathbf{0}}}^{1}(\Omega^{\ast};\mathbb{R}^{d})$,
$d\in\\{2,3\\}$. satisfies the equation

$\displaystyle\
\int_{\Omega^{\ast}}{\left(\nabla{\bm{\mathbf{V}}}:\nabla\bm{\mathbf{\varphi}}+{\bm{\mathbf{V}}}\cdot\bm{\mathbf{\varphi}}\right)}{\,{\operatorname{\mathnormal{d}}}x}$
$\displaystyle=\int_{\Sigma^{\ast}}{\Phi(\cdot;\Sigma^{\ast}){\bm{\mathbf{n}}}\cdot\bm{\mathbf{\varphi}}}{\,{\operatorname{\mathnormal{d}}}\sigma},\quad\forall\bm{\mathbf{\varphi}}\in\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega^{\ast}),$
(37) $\displaystyle\text{and}\qquad\bm{\mathbf{V}}\cdot\bm{\mathbf{n}}$
$\displaystyle=0\quad\text{on $\Sigma^{\ast}$}.$

For a domain $\Omega^{\ast}$ of class $\mathcal{C}^{0,1}$ and function
$\Phi(\,\cdot\,;\Sigma)\in L^{2}(\Sigma)$, the variational problem (37) can be
shown to have a weak solution
${\bm{\mathbf{V}}}\in\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega^{\ast})$
via Lax-Milgram lemma. With the definitions given above, we will now issue the
main point of this subsection which is given in the next proposition. Here, we
will tacitly assume – for the sake of argument – that $\Phi$ vanishes within a
short time interval and that we have the convergence of the evolving domains
to a stationary point $\Omega^{\ast}$ at that time interval without referring
to a formal mathematical notion of convergence of sets. On a related note, we
point out that small and smooth perturbations of a regular domain may be
“uniquely” described by normal deformations of the boundary of the domain, see
[49].

######  Proposition 5.1.

Let $\Omega^{\ast}\in\mathcal{C}^{\circ}$ and $\Phi(\,\cdot\,;\Sigma)\in
L^{2}(\Sigma)\cap\mathcal{C}^{\circ}(\Omega)$. Then, $\Omega^{\ast}$ is a
stationary solution to Problem 2 if and only if $\Omega^{\ast}$ is a
stationary solution to Problem 3.

###### Proof.

Consider equation (36) over the stationary shape $\overline{\Omega^{\ast}}$
with Lipschitz boundary $\partial\Omega^{\ast}=\Gamma\cup\Sigma^{\ast}$. For
the necessity part, we assume that
$L^{2}(\Sigma)\cap\mathcal{C}^{\circ}(\Omega)\ni\Phi(\,\cdot\,;\Sigma)=0$, and
we need to verify that $\bm{\mathbf{V}}\cdot\bm{\mathbf{n}}=0$ on
$\Sigma^{\ast}$. To this end, we multiply the first equation in (37) by
$\bm{\mathbf{V}}\in\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega^{\ast})\cap\bm{\mathbf{\mathcal{C}}}^{\circ}(\Omega^{\ast})$
and then apply integration by parts – noting that
$\bm{\mathbf{V}}=\bm{\mathbf{0}}$ on $\Gamma$ – to obtain
$0\leqslant\int_{\Omega^{\ast}}{|\bm{\mathbf{V}}|^{2}}{\,{\operatorname{\mathnormal{d}}}x}=-\int_{\Sigma^{\ast}}{\partial_{\bm{\mathbf{n}}}{\bm{\mathbf{V}}}\cdot\bm{\mathbf{V}}}{\,{\operatorname{\mathnormal{d}}}\sigma}=-\int_{\Omega^{\ast}}{|\nabla\bm{\mathbf{V}}|^{2}}{\,{\operatorname{\mathnormal{d}}}x}\leqslant
0$. Clearly, $\bm{\mathbf{V}}=\bm{\mathbf{0}}$ in $\Omega^{\ast}$. Moreover,
because $\bm{\mathbf{V}}\big{|}_{\Gamma}=\bm{\mathbf{0}}$, then, by the
maximum principle, $\bm{\mathbf{V}}\equiv\bm{\mathbf{0}}$ on
$\overline{\Omega}^{\ast}$, In particular, we have
$\bm{\mathbf{V}}\cdot{\bm{\mathbf{n}}}=0$ on $\Sigma^{\ast}$.

For the sufficiency part, we need to prove that if
$\bm{\mathbf{V}}\cdot\bm{\mathbf{n}}=0$ on $\Sigma^{\ast}$, where
$\bm{\mathbf{V}}$ solves problem (36) on $\overline{\Omega^{\ast}}$, then
$\Phi=0$ on $\Sigma^{\ast}$. We take
$\bm{\mathbf{\varphi}}=\bm{\mathbf{V}}\in\bm{\mathbf{H}}_{\Gamma,\bm{\mathbf{0}}}^{1}(\Omega^{\ast})$
in (37) from which we get
$\int_{\Omega^{\ast}}{\left(\nabla{\bm{\mathbf{V}}}:\nabla{\bm{\mathbf{V}}}+{\bm{\mathbf{V}}}\cdot{\bm{\mathbf{V}}}\right)}{\,{\operatorname{\mathnormal{d}}}x}=\int_{\Sigma^{\ast}}{\Phi(\,\cdot\,;\Sigma^{\ast}){\bm{\mathbf{n}}}\cdot\bm{\mathbf{V}}}{\,{\operatorname{\mathnormal{d}}}\sigma}=0$.
Clearly, $\bm{\mathbf{V}}\equiv\bm{\mathbf{0}}$ on $\overline{\Omega}^{\ast}$.
Going back to (37), we obtain
$\int_{\Sigma^{\ast}}{\Phi(\,\cdot\,;\Sigma^{\ast}){\bm{\mathbf{n}}}\cdot\bm{\mathbf{\varphi}}}{\,{\operatorname{\mathnormal{d}}}\sigma}=0$,
for all
$\bm{\mathbf{\varphi}}\in\bm{\mathbf{H}}_{\Gamma,\bm{\mathbf{0}}}^{1}(\Omega^{\ast})$,
from which we conclude that $\Phi=0$ on $\Sigma^{\ast}$. This proves the
assertion. ∎

Let us look at the situation when $\Phi=-G$ in (35). In this case, $\Phi$ not
only depends on some geometric quantities on the free boundary, but also to
some functions which are solutions to specific equality constraints (the state
and adjoint state equations to be exact). For a domain $\Omega$ that is
$\mathcal{C}^{k,\alpha}$ regular, $k\in\mathbb{N}$, $k\geqslant 2$, it can be
shown that both the state and the adjoint state systems (6) and (14),
respectively, admit a unique (classical) solution in the space
$\mathcal{C}^{k,\alpha}(\overline{\Omega};\mathbb{C}^{d})$. In particular, the
imaginary part $u_{2}$ of the state variable $u$ is an element of the set
$\mathcal{C}^{k,\alpha}(\overline{\Omega};\mathbb{C}^{d})$. This implies that
the adjoint state $p$ is even more regular, and is in fact
$\mathcal{C}^{k+2,\alpha}$ smooth in $\overline{\Omega}$. Meanwhile, the
normal vector $\bm{\mathbf{n}}$ to $\Omega\in\mathcal{C}^{3,\alpha}$ is only
$\bm{\mathbf{\mathcal{C}}}^{2,\alpha}$ regular, implying that the mean
curvature $\kappa$ is $\mathcal{C}^{1,\alpha}(\partial\Omega)$ smooth. Since
$G$ consists of $\kappa$, it appears that one even needs the initial domain
$\Omega^{0}$ be at least $\mathcal{C}^{4,\alpha}$ regular for the evolving
domain $\Omega(t)$ to be $\mathcal{C}^{2,\alpha}$ smooth for some short time
interval $[0,\varepsilon)$. On the other hand, considering Problem 3 with
$\Phi=-G$, it appears that it is enough to assume that
$\Omega^{0}\in\mathcal{C}^{3,\alpha}$ for the moving domain $\Omega(t)$ to be
$\mathcal{C}^{2,\alpha}$ smooth in some short time interval
$[0,\varepsilon)\ni t$. Indeed, in this case,
$\Phi\bm{\mathbf{n}}=-G\bm{\mathbf{n}}\in\bm{\mathbf{\mathcal{C}}}^{1,\alpha}(\Sigma(t))$,
and so $\bm{\mathbf{V}}\in\bm{\mathbf{\mathcal{C}}}^{2,\alpha}(\Omega(t))$,
for $t\in[0,\varepsilon)$. Hence, it can actually be shown that the weak
solution
${\bm{\mathbf{V}}}\in\bm{\mathbf{H}}^{1}_{\Gamma,\bm{\mathbf{0}}}(\Omega(t))$
of (36) is also a classical solution of the problem over the small time
interval $[0,\varepsilon)$. Therefore, $\Omega(t)$ in (36) remains
$\mathcal{C}^{2,\alpha}$ smooth throughout the short time interval
$[0,\varepsilon)$.

Although not so important, we provide additional comments about the evolving
boundaries. The evolution of the moving boundary $\Sigma(t)$ due to (35) and
that of (36) through time are not the same. Indeed, only in the case that
Proposition 5.1 is true we are sure that the two evolving boundaries coincide
(except of course with the initial shape). Nevertheless, a more accurate
extended-regularized version of (35) can be formulated by replacing the
Neumann condition in (36) by the Dirichlet condition
$\bm{\mathbf{V}}=\Phi(\,\cdot\,;\Sigma(t)){\bm{\mathbf{n}}}$ on $\Sigma(t)$,
for $t>0$. In this case, however, the regularity of the moving boundary
$\Sigma(t)$ is not preserved when the evolving boundary evolves according to
the normal speed $V_{n}=\bm{\mathbf{V}}\cdot\bm{\mathbf{n}}$. Even so, one can
address the issue to some extent by considering an approximation of the
Dirichlet condition using a Robin condition. That is, we can add more
regularity to the vector $V_{n}\bm{\mathbf{n}}$ by setting
$\beta\nabla\bm{\mathbf{V}}+\bm{\mathbf{V}}=\Phi(\,\cdot\,;\Sigma(t)){\bm{\mathbf{n}}}$
on $\Sigma(t)$, for $t>0$, where $\beta>0$. Still, with this alternative
formulation, there is a trade-off (controlled by $\beta$) between the accuracy
of the evolution and the regularity of the moving boundary. Here, we do not
bother about the accuracy of the extended-regularized evolving boundary with
respect to the corresponding original one as we are only interested in
extending in $\Omega(t)$ the normal velocity vector $\Phi\bm{\mathbf{n}}$ of
the evolving boundary $\Sigma(t)$ while adding more regularity to it.

###  5.3 Numerical examples

We now illustrate the feasibility and applicability of the new method in
solving concrete examples of the free boundary problem (2). We first test our
method in two dimensions, and carry out a comparison with KVM (Examples 1–3).
Also, since most of the previous studies only dealt with problems in two
dimensions (except in [24] and also in [34] which applies the Newton scheme to
the Dirichlet energy functional), we also put our attention on testing CCBM to
three dimensional cases (Examples 4–8). In the case of three spatial
dimensions, we let $\lambda=-10$, unless specified.

######  Remark 12 (Details of the computational setup and environment).

The numerical simulations conducted here are all implemented in the
programming software FreeFem++, see [39]. Every variational problem involved
in the procedure is solved using $\mathbb{P}1$ finite element discretization
and are solved in FreeFEM++ via the command line problem with the default
solver (sparsesolver or LU if any direct sparse solver is available999see,
e.g., p. 380 of FreeFem Documentation, Release 4.8, May 16, 2022,
www.freefem.org.). Moreover, all mesh deformations are carried out without any
kind of adaptive mesh refinement as opposed to what has been usually done in
earlier works, see [51, 52, 53]. We emphasize that, in this way, we can
further assess the stability of CCBM in comparison with KVM. The computations
are all performed on a MacBook Pro with Apple M1 chip computer having 16GB RAM
processors.

######  Remark 13.

As pointed out in Remark 12, we avoid the generation of a new triangulation of
the domain at every iterative step. Obviously, this is achieved by moving not
only the boundary, but also the internal nodes of the mesh triangulation at
every iteration (as mentioned earlier). By doing so, the mesh only needs to be
generated at initial iteration. In order to move the boundary and internal
nodes simultaneously, we solve the discretized version of (36) and then move
the domain in the direction of the resulting vector field scaled with the
pseudo-time step size $t^{k}$ (see Remark 9 for the computation of $t^{k}$).
That is, we find
${{\bm{\mathbf{V}}_{h}^{k}}}\in\mathbb{P}_{1}(\Omega_{h}^{k})^{d}$ such that
it solves the equation

$-\Delta{\bm{\mathbf{V}}_{h}^{k}}+{\bm{\mathbf{V}}_{h}^{k}}=\bm{\mathbf{0}}\ \
\text{in $\Omega_{h}^{k}$},\qquad{\bm{\mathbf{V}}_{h}^{k}}=\bm{\mathbf{0}}\ \
\text{on
$\Gamma^{h}$},\qquad{\nabla{\bm{\mathbf{V}}_{h}^{k}}}\cdot{\bm{\mathbf{n}}_{h}^{k}}=-G^{k}{\bm{\mathbf{n}}_{h}^{k}}\
\ \text{on $\Sigma_{h}^{k}$},$

where we suppose a polygonal domain $\overline{\Omega_{h}^{k}}$ and its
triangulation
$\mathcal{T}_{h}(\overline{\Omega_{h}^{k}})=\\{K^{k}_{l}\\}^{N_{e}}_{l=1}$
($K^{k}_{l}$ is a closed triangle for $d=2$, or a closed tetrahedron for
$d=3$) are given, and $\mathbb{P}_{1}(\Omega_{h}^{k})^{d}$ denotes the
$\mathbb{R}^{d}$-valued piecewise linear function space on
$\mathcal{T}_{h}(\overline{\Omega_{h}^{k}})$. Then, we update the domain or,
equivalently, move the nodes of the mesh by defining $\Omega^{k+1}_{h}$ and
$\mathcal{T}_{h}(\overline{\Omega^{k+1}_{h}})=\\{K^{k+1}_{l}\\}^{N_{e}}_{l=1}$
respectively as
$\overline{\Omega^{k+1}_{h}}:=\left\\{x+t^{k}{\bm{\mathbf{V}}_{h}^{k}}(x)\
\middle|\ x\in\overline{\Omega_{h}^{k}}\right\\}$ and
$K^{k+1}_{l}:=\left\\{x+t^{k}{\bm{\mathbf{V}}_{h}^{k}}(x)\ \middle|\ x\in
K^{k}_{l}\right\\}$, for all $k=0,1,\ldots$.

We are now ready to give our first numerical example.

######  Example 1 (Testing the accuracy of the gradient).

Consider two concentric circles centered at the origin $\mathbf{0}$ with
radius $r>0$ and $R>r$ given by $C(\bm{0},r)$ and $C(\bm{0},R)$, respectively.
Then, problem (5) can be expressed as the PDE system

$-\partial_{\rho\rho}^{2}u-\rho^{-1}\partial_{\rho}u=0\ \ \text{for
$r<\rho<R$},\qquad u(r)=1,\qquad\text{and}\qquad u(R)=0,$

whose exact solution can be computed as $u(\rho)=\log{(\rho/R)}/\log{(r/R)}$.
Moreover, in this case, $\partial_{\bm{\mathbf{n}}}{u(R)}=1/[R\log{(r/R)}]$.
Therefore, the exterior Bernoulli FBP (2) with
$\Gamma=\\{x\in\mathbb{R}^{2}:|x|=r\\}$ and
$\lambda=1/\left[R\log{\left(r/R\right)}\right]$, $0<r<R$, has the unique
exact free boundary solution $\Sigma^{\ast}=C(\bm{0},R)$. With this in mind,
we let $r=0.5$ (i.e., $\Gamma=\Gamma_{\texttt{C}}:=C(\bm{0},0.5)$) and
$R^{\ast}=0.7$, which gives us $\lambda=-4.24573$, and set
$\Sigma^{0}=C(\bm{0},1.25)$ as our initial guess. In this experiment, we
examine the sensitivity of the cost functions $J$ and $J_{KV}$ by testing the
methods with large variations. To this end, we consider three test cases for
KV by varying the step size parameter $\mu$: (KV)1: $\mu=2.0$; (KV)2:
$\mu=1.0$; and (KV)3: $\mu=0.5$, while keeping $\mu=2.0$ for CCBM. Recall that
$\mu$ dictates how large the magnitude of $t^{k}$ can be at every iteration,
see Remark 9. Moreover, we discretize the initial domain with uniform mesh
sizes and look at the effect of accuracy of the methods under different mesh
widths. Finally, we stop our algorithm as soon as the absolute difference
between consecutive cost values is less than $10^{-6}$. Figure 1 shows the
histories of cost values for KVM and CCBM, as well as the histories of
Hausdorff distances $d_{\text{H}}(\Sigma^{k},\Sigma^{\ast})$. A summary of
Hausdorff distances, cpu-time, and cpu-time-per-iteration against mesh sizes
$h=0.2,0.1,0.05,0.025,0.0125$ in the form of plots are plotted in Figure 2.
Based from these results, we draw the following observations:

  * •

$J_{KV}$ is less sensitive to $J$ in terms of large variations;

  * •

KVM tends to converge prematurely for coarse meshes unlike CCBM;

  * •

KVM and CCBM nearly have the same convergence behavior for finer meshes;

  * •

KVM and CCBM almost have the same accuracy for finer meshes, but CCBM is more
accurate in general (primarily because of the second point);

  * •

KVM and CCBM complete the iteration procedure at almost the same time,

  * •

however, CCBM requires less computing-time-per-iteration than KVM,

  * •

and the latter converges in less number of iterations (hence, the previous
point).

Based from these observations, it seems that, in the case of two dimensional
problems, CCBM features some merits over the KVM in terms of computational
aspects. Of course, we expect that these advantages can be exploited
especially when dealing with three dimensional problems as we can instead
resort to coarse meshes without concerning much ourselves with the accuracy of
the approximation process – at least when dealing with axisymmetric cases.
Also, it appears that we have more freedom to take large step sizes in the
case of CCBM than when applying KVM as the former is less prone to premature
convergence. However, as we already mentioned, KVM requires less number of
iterations than CCBM, and that, under small step sizes, the two methods
converges to almost identical optimal shape solutions.

Figure 1: Histories of cost values (first row) and Hausdorff distances (second
row) for Example 1

Figure 2: Hausdorff distances, cpu-time, and cpu-time-per-iteration against
the mesh size $h$

In the next two examples, we carry out three different experiments for KVM to
further highlight the sensitivity of $J$ in comparison with $J_{KV}$ in the
case of a bit more complex geometries for the fixed boundary with concave
regions. We again test the effect of the step size parameter $\mu$ in the
approximation process. For CCBM, we fix $\mu$ to $1.0$. The test experiments
are as follows: test (KV)a: $\mu=1.0$; test (KV)b: $\mu=0.5$; and test (KV)c:
$\mu=0.25$, and we look at the problems where $\lambda=-10,-9,\ldots,-1$.
Also, we choose $C(\bm{0},1.25)$ as the initial profile $\Sigma^{0}$ of the
free boundary, and we stop the algorithm after completing $100$ iterations.

######  Example 2 (An L-shape fixed boundary).

Let us consider an L-shape geometry for $\Gamma$ given by the boundary of the
domain $D=(-0.25,0.25)^{2}\setminus[0,0.25]^{2}$, i.e., we set
$\Gamma=\Gamma_{\texttt{L}}:=\partial D$. A comparison between the computed
free boundaries using the KVM and CCBM for each test experiments are shown in
Figure 3 (first row). In those figures, the outermost boundary corresponds to
$\lambda=-1$ while the innermost (exterior) boundary corresponds to
$\lambda=-10$. The histories of $d_{\text{H}}(\Sigma^{k},\Sigma_{100})$-values
are shown in Figure 4 while the histories of cost values are plotted in Figure
5. It appears based on the results that CCBM almost has the same convergence
behavior with KVM in the case of larger values for $\lambda$, and a bit faster
for the case of smaller values of $\lambda$. It must be evident, however, that
KVM tends to converge to a stationary shape in a fewer number of iterations
than CCBM, but to a less accurate geometry for the free boundary. In fact, KVM
converges prematurely in some instances; see, for example, Figure 8 where the
evolutions of the shapes (plotted at every ten iterates with $\mu=2.0$) are
illustrated in the case $\lambda=-5$. Nevertheless, the two methods with small
step sizes provide almost identical optimal solutions to the minimization
problem. Furthermore, we notice from the histories of Hausdorff distances and
cost values that $J$ is more sensitive to large variations than $J_{KV}$.
Based on these observations, we can say that CCBM has some advantages over KVM
in terms of computational performance, specifically with respect to overall
computing-time-per-iteration (see left plot in Figure 9) and accuracy when
taking large deformations of the domain. However, as evident in Figures 4–5,
KVM converges in fewer iterations than CCBM.

######  Example 3 (A ribbon shape fixed boundary).

Next, we consider a ribbon shape fixed boundary similar to the one examined in
[22] which is parametrized as
$\Gamma_{\texttt{R}}:=\\{(0.45\cos\theta,0.3\sin\theta(1.25+\cos
2\theta))^{\top}\mid 0\leqslant\theta\leqslant 2\pi\\}$. The results of the
test experiments are depicted in Figure 3 (second row). Meanwhile, the
histories of the Hausdorff distance values
$d_{\text{H}}(\Sigma^{k},\Sigma_{100})$ and of the cost values are shown in
Figure 6 and Figure 7, respectively. In these computational results, we
observe similar findings to the previous example in that $J_{KV}$ is less
sensitive to large variations than $J$. Moreover, in some cases, KVM tends to
overshoot the approximate optimal shape unlike CCBM. Even so, as in the case
of the previous test experiments, the present algorithm with the KVM uses
fewer number of iterations to converge compared to CCBM.

Figure 3: Computational results for Example 2 (first row) and Example 3
(second row)

Figure 4: iteration $k$ vs distance $d_{H}(\Sigma^{k},\Sigma_{100})$ (legend:
— (KV)a — (KV)b — (KV)c — (CCBM))

Figure 5: iteration $k$ vs cost values (legend: — (KV)a — (KV)b — (KV)c —
(CCBM))

Figure 6: iteration $k$ vs distance $d_{H}(\Sigma^{k},\Sigma_{100})$ (legend:
— (KV)a — (KV)b — (KV)c — (CCBM))

Figure 7: iteration $k$ vs cost values (legend: — (KV)a — (KV)b — (KV)c —
(CCBM))

Figure 8: Evolution of shapes (plotted at every ten iterates and were obtained
with step size parameter value $\mu=2.0$) for Example 2 when $\lambda=-7$
using KVM (left-most plot) and CCBM (middle plot), and a direct comparison
with the computed optimal shapes (right-most plot)

Figure 9: $\lambda$ vs computational-time-per-iteration for Example 2 (left)
and Example 3 (right)

The rest of our examples will focus on three dimensional cases using CCBM.

######  Example 4 (Axisymmetric case in 3D).

Let us first test the method to a simple axisymmetric 3D-case with an
analytical solution. On this purpose, we consider the spheres
$\mathcal{S}(\bm{\mathbf{0}},r):=\\{\xi\in\mathbb{R}^{3}:|\xi|=r\\}$ and
$\mathcal{S}(\bm{\mathbf{0}},R):=\\{\xi\in\mathbb{R}^{3}:|\xi|=R\\}$ centered
at $\bm{\mathbf{0}}$ with radius $r>0$ and $R>r$, respectively. With $u(r)=1$
and $u(R)=0$, the solution to the Dirichlet problem (5) is exactly given by
$u(\rho)=r(R-\rho)/[\rho(R-r)]$, $\rho\in(r,R)$ with normal derivative
$\partial_{\rho}u(\rho)=-Rr/[\rho^{2}(R-r)]$. So, on the exterior surface, we
have $\partial_{\rho}u(R)=-r/[R(R-r)]=:\lambda$. Thus, problem (2) with
$\Gamma=\mathcal{S}(\bm{\mathbf{0}},r)$ and $\lambda=-r/[R(R-r)]$, $0<r<R$,
has the unique exact free boundary solution
$\Sigma^{\ast}=\mathcal{S}(\bm{\mathbf{0}},R^{\ast})$. For a concrete example,
we let $r=0.3$ and $R^{\ast}=0.5$, giving us $\lambda=-3$, and take
$\Sigma^{0}=\mathcal{S}(\bm{\mathbf{0}},0.6)$ as the initial guess. With
$\texttt{Tol}=10^{-6}$, initial maximum mesh size $h_{\max}=0.1$ on the
surfaces, and maximum volume $0.001$ for the tetrahedra, the procedure is
completed after $53$ sec. The nodes on $\Sigma_{f}$ have mean radii of
$\bar{R}=0.4787$ which gives a $4.26\%$ error with respect to the exact radius
$R^{\ast}=0.5$.

In the next three examples, the algorithm is terminated after 600 iterations
in addition to setting $\texttt{Tol}=10^{-8}$ in the stopping condition (see
Remark 10).

######  Example 5 (Perturbed sphere).

Let us define $\Gamma$ as the surface of a perturbed sphere having strict
concave regions. With $\Sigma^{0}=\mathcal{S}(\bm{\mathbf{0}},1.5)$ as the
initial guess, the results of the computation are summarized in Figure 10
(first row), including the plots for histories of cost values and Sobolev
gradient norms.

######  Example 6 (Torus).

We also look at the case where the fixed surface $\Gamma$ is given by a torus
and set $\Sigma^{0}=\mathcal{S}(\bm{\mathbf{0}},0.8)$ as our initial guess.
The computational results for this example are shown in Figure 10 (second
row).

######  Example 7 (Four disjoint spheres).

Let us also consider the case where $\Gamma$ is the union of four disjoints
spheres having the same exact radius $r=0.25$, and let
$\Sigma^{0}=\mathcal{S}(\bm{\mathbf{0}},0.9)$. The results for this test case
are summarized in Figure 10 (third row).

A summary of mesh details (number of boundary elements, triangles, and
vertices) used in the last three examples are tabulated in Table 1. The
computed cost at final iterate and over-all cpu times for the three examples
are also shown in the table.

Figure 10: Computational results for Examples 5–7 | Example 5 | Example 6 | Example 7
---|---|---|---
number of boundary elements | 30,762 | 9,232 | 20,838
number of triangles | 65,272 | 21,549 | 49,900
number of vertices | 18,879 | 5,905 | 13,629
cost | 2.92e-10 | 5.26e-09 | 4.84e-06
cpu time | 14,436 s | 2,212 s | 9,680 s
Table 1: Mesh details and additional computational results for Examples 5–7

######  Example 8 (L-block figure).

Finally, we consider $\Sigma$ as the surface of an L-block figure, and
consider the values $\lambda=-1,-7,-10$. The results of the computations are
shown in Figure 11. The over-all cpu times for the case $\lambda=-7$ and $-10$
are less than $2,200$ sec, while the iteration process was completed after
$6,200$ sec for $\lambda=-1$.

Figure 11: Computational results for Example 8

##  6 Conclusions and Future Works

We have proposed here a complex coupled boundary method in shape optimization
framework as a numerical resolution to the exterior Bernoulli problem. The
shape gradient of the cost corresponding to the formulation is computed under
a mild regularity assumption on the domain, and by using only the weak form of
the equation satisfied by the material derivative of the state problem whose
existence is shown in a rigorous manner. The shape Hessian at a critical shape
is also characterized through the chain rule approach under enough smoothness
assumption on the domain. Also, by the same technique used to derive the shape
gradient, the same expression is recomputed, but now with a weaker assumption
on the regularity of the domain. The aforesaid expression is then analyzed in
order to study the algebraic ill-posedness of the proposed method which is
done by showing that the Riesz operator associated to the quadratic shape