On second order shape optimization methods for electrical impedance tomography

TL;DR

This paper introduces a second order method for shape optimization in electrical impedance tomography, analyzing derivatives, stability, and ill-posedness to improve shape recovery accuracy.

Contribution

It is the first to develop a second order approach for shape optimization in electrical impedance tomography, including derivative analysis and stability assessment.

Findings

Derived second order derivatives of the state with respect to shape perturbations.

Proposed a Newton method for shape reconstruction.

Proved the ill-posedness via Hessian compactness at the solution.

Abstract

This paper is devoted to the analysis of a second order method for recovering the \emph{a priori} unknown shape of an inclusion $\omega$ inside a body $\Omega$ from boundary measurement. This inverse problem - known as electrical impedance tomography - has many important practical applications and hence has focussed much attention during the last years. However, to our best knowledge, no work has yet considered a second order approach for this problem. This paper aims to fill that void: we investigate the existence of second order derivative of the state $u$ with respect to perturbations of the shape of the interface $\partial\omega$, then we choose a cost function in order to recover the geometry of $\partial \omega$ and derive the expression of the derivatives needed to implement the corresponding Newton method. We then investigate the stability of the process and explain why this inverse problem is severely ill-posed by proving the compactness of the Hessian at the global minimizer.

Equations

$$-\text{mbox}\mathrm{div}\left(\sigma (x)\nabla u\right)=0\text{ in }\Omega ,$$

$$u=f\text{ on }\partial \Omega .$$

$${\Lambda }_{\sigma }:f\mapsto \sigma {\left({\partial }_{\mathbf{n}}u\right)}_{|\partial \Omega },$$

$$\sigma ={\sigma }_1{\chi }_{\Omega \setminus \omega }+{\sigma }_2{\chi }_{\omega }.$$

$${\sigma }_1{\partial }_{\mathbf{n}}u=g\text{ on }\partial \Omega ,$$

$$J_{KV}(\omega )={\int }_{\Omega }\sigma |\nabla (u_d-u_n){|}^2$$

$$\mathcal{H}=\{\mathbf{h}\in {\mathcal{C}}^{4,\alpha }({\mathbb{R}}^d,{\mathbb{R}}^d),Supp(\mathbf{h})\subset {\Omega }_{d_0}\}.$$

$$\left\{\begin{array}[]{rcl}\Delta u_{d}^{\prime}&=&0\text{ in }\Omega\setminus\overline{\omega}\text{ and in }\omega,\\ \left[u_{d}^{\prime}\right]&=&h_{n}~{}\cfrac{[\sigma]}{\sigma_{1}}~{}\partial_{\mathbf{n}}u_{d}^{-}~{}\text{ on }\partial\omega,\\ \left[\sigma\partial_{n}u_{d}^{\prime}\right]&=&[\sigma]\mbox{\rm div}_{\tau}\left(h_{n}\nabla_{\tau}u_{d}\right)\text{ on }\partial\omega,\\ u_{d}^{\prime}&=&0\text{ on }\partial\Omega.\end{array}\right.$$

$$\left\{\begin{array}[]{rcl}\Delta u_{n}^{\prime}&=&0\text{ in }\Omega\setminus\overline{\omega}\text{ and in }\omega,\\ \left[u_{n}^{\prime}\right]&=&h_{n}~{}\cfrac{[\sigma]}{\sigma_{1}}~{}\partial_{\mathbf{n}}u_{n}^{-}~{}\text{ on }\partial\omega,\\ \left[\sigma\partial_{n}u_{n}^{\prime}\right]&=&[\sigma]\mbox{\rm div}_{\tau}\left(h_{n}\nabla_{\tau}u_{n}\right)\text{ on }\partial\omega,\\ \partial u_{n}^{\prime}&=&0\text{ on }\partial\Omega.\end{array}\right.$$

$$\left\{\begin{array}[]{rcl}\Delta u_{d}^{\prime\prime}&=&0\text{ in }\Omega\setminus\overline{\omega}\text{ and in }\omega,\\ \left[u_{d}^{\prime\prime}\right]&=&\left(h_{1,n}h_{2,n}H-\mathbf{{{\bf{h}}}_{1}}_{\tau}.(D\mathbf{n}\,\mathbf{{{\bf{h}}}_{2}}_{\tau})\right)[\partial_{\mathbf{n}}u_{d}]-\left(h_{1,n}[\partial_{\mathbf{n}}(u_{d})^{\prime}_{2}]+h_{2,n}[\partial_{\mathbf{n}}(u_{d})^{\prime}_{1}]\right)\\ &&~{}~{}~{}~{}~{}~{}+\left(\mathbf{{{\bf{h}}}_{1}}_{\tau}.\nabla h_{2,n}+\mathbf{{{\bf{h}}}_{2}}_{\tau}.\nabla h_{1,n}\right)[\partial_{\mathbf{n}}u_{d}]\text{ on }\partial\omega,\\ \left[\sigma\partial_{n}u_{d}^{\prime\prime}\right]&=&\mbox{\rm div}_{\tau}\left(h_{2,n}\left[\sigma\nabla_{\tau}(u_{d})^{\prime}_{1}\right]+h_{1,n}\left[\sigma\nabla_{\tau}(u_{d})^{\prime}_{2}\right]+\mathbf{h_{1}}_{\tau}.(D\mathbf{n}\,\mathbf{h_{2}}_{\tau})[\sigma\nabla_{\tau}u_{d}]\right)\\ &&~{}~{}~{}~{}~{}~{}-\mbox{\rm div}_{\tau}\left((\mathbf{{{\bf{h}}}_{1}}_{\tau}.\nabla_{\tau}h_{2,n}+\nabla_{\tau}h_{1,n}.\mathbf{h_{2}}_{\tau})\left[\sigma\nabla_{\tau}u_{d}\right]\right)\\ &&~{}~{}~{}~{}~{}~{}+\mbox{\rm div}_{\tau}\left(h_{2,n}h_{1,n}(2D\mathbf{n}-HI)\left[\sigma\nabla_{\tau}u_{d}\right]\right)\text{ on }\partial\omega,\\ u_{d}^{\prime\prime}&=&0\text{ on }\partial\Omega.\end{array}\right.$$

$$D{J}_{KV}(\omega )\mathbf{h}=[\sigma ]{\int }_{\partial \omega }\left[\frac{{\sigma }_1}{{\sigma }_2}\left(|{\partial }_{\mathbf{n}}{u}_d^{+}{|}^2-|{\partial }_{\mathbf{n}}{u}_n^{+}{|}^2\right)+|{\nabla }_{\tau }{u}_d{|}^2-|{\nabla }_{\tau }{u}_n{|}^2\right]h_n.$$

$$\begin{split}D^{2}J_{KV}(\omega)({\bf{h}}_{1},{\bf{h}}_{2})&=\int_{\partial\omega}\left[\sigma|\nabla v|^{2}\right]\left(\mathbf{{\bf{h}}_{1}}_{\tau}.\nabla_{\tau}(h_{2,n})+\mathbf{{\bf{h}}_{2}}_{\tau}.\nabla_{\tau}(h_{1,n})-\mathbf{h_{2}}_{\tau}.(D\mathbf{n}\,\mathbf{{\bf{h}}_{1}}_{\tau})\right)\\ &-\int_{\partial\omega}\partial_{\mathbf{n}}\left(\left[\sigma|\nabla v|^{2}\right]\right)h_{1,n}h_{2,n}+2\big{[}\sigma\nabla v.(h_{1,n}\nabla v^{\prime}_{2}+h_{2,n}\nabla v^{\prime}_{1})\big{]}\\ &+\int_{\partial\omega}\left[\sigma\left(\partial_{\mathbf{n}}(u_{n})^{\prime}_{1}v^{\prime}_{2}+\partial_{\mathbf{n}}(u_{n})^{\prime}_{2}v^{\prime}_{1}-\partial_{\mathbf{n}}v^{\prime}_{1}(u_{d})^{\prime}_{2}-\partial_{\mathbf{n}}v^{\prime}_{2}(u_{d})^{\prime}_{1}\right)\right]\\ &+2\int_{\partial\omega}v\left[\sigma\partial_{\mathbf{n}}(u_{n})^{\prime\prime}_{1,2}\right]-\sigma_{1}\partial_{\mathbf{n}}v^{+}\left[(u_{d})^{\prime\prime}_{1,2}\right]\end{split}$$

$$D{J}_{KV}({\omega }^{\ast })(\mathbf{h})=0,$$

$$D^2{J}_{KV}({\omega }^{\ast })(\mathbf{h},\mathbf{h})={\int }_{\Omega }\sigma |\nabla v^{\prime }{|}^2.$$

$$\dot{y}(\Omega ;\mathbf{h})=\underset{t\to 0}{\lim }\frac{y({\Omega }_t)\circ T_t-y(\Omega )}{t},$$

$$y^{\prime }=\dot{y}-\mathbf{h}.\nabla y.$$

$$\text{mbox}{\mathrm{div}}_{\tau }\left(\mathbf{V}\right)=\text{mbox}\mathrm{div}\left(\mathbf{V}\right)-D\mathbf{V}.\mathbf{n}.\mathbf{n},$$

$$\text{mbox}{\mathrm{div}}_{\tau }\left(\mathbf{V}\right)=\text{mbox}\mathrm{div}\left(\tilde{\mathbf{V}}\right)-(D\tilde{\mathbf{V}}.\mathbf{n}).\mathbf{n},$$

$${\nabla }_{\tau }f=\nabla \tilde{f}{|}_{\partial \Omega }-\nabla \tilde{f}.\mathbf{n}\mathbf{n}\text{ on }\partial \Omega .$$

$${\int }_{\partial \Omega }{\nabla }_{\tau }f.\mathbf{F}=-{\int }_{\partial \Omega }f\text{mbox}{\mathrm{div}}_{\tau }\left(\mathbf{F}\right),$$

$$\text{mbox}\mathrm{div}\left(\mathcal{N}\right)=H\text{ on }\partial \Omega .$$

$${\int }_{\partial \Omega }\nabla f.\mathbf{F}+f\text{mbox}{\mathrm{div}}_{\tau }\left(\mathbf{F}\right)={\int }_{\partial \Omega }\left(\nabla f.\mathbf{n}+Hf\right)\mathbf{F}.\mathbf{n}.$$

$${\Delta }_{\tau }f=\text{mbox}{\mathrm{div}}_{\tau }\left({\nabla }_{\tau }f\right).$$

$$\Delta f={\Delta }_{\tau }f+H{\partial }_{\mathbf{n}}f+{\partial }_{nn}^{2}f,\text{ on }\partial \Omega .$$

$$\dot{\mathbf{n}}=-{\nabla }_{\tau }(\mathbf{h}.\mathbf{n})+D\mathbf{n}{\mathbf{h}}_{\tau },$$

$${\mathbf{n}}^{\prime }=-{\nabla }_{\tau }(\mathbf{h}.\mathbf{n}).$$

$$\dot{\overline{{\nabla }_{\tau }f}}=\nabla \dot{f}+(D^{2}f\mathbf{h}{)}_{\tau }-\nabla f.\mathbf{n}\dot{n}-\nabla f.\dot{\mathbf{n}}\mathbf{n}$$

$$\dot{\overline{\nabla f}}=\nabla f^{\prime }+D^{2}f\mathbf{h}$$

$$\dot{\overline{\nabla f.\mathbf{n}\mathbf{n}}}=\nabla f.\dot{\mathbf{n}}\mathbf{n}+\nabla f.\mathbf{n}\dot{\mathbf{n}}+\nabla f^{\prime }.\mathbf{n}\mathbf{n}+(D^{2}f\mathbf{h}).\mathbf{n}\mathbf{n}.$$

$$\int_{\partial\Omega}\left[\left(\Delta_{\tau,t}f\right)\circ T_{t}~{}\gamma_{\tau}(t)\right]\phi~{}=-\int_{\partial\Omega}\left[C(t)~{}\Big{(}\nabla(f\circ T_{t})-(B(t)~{}\mathbf{n}).\nabla(f\circ T_{t})\Big{)}\right].\nabla\phi,\,\,\forall\phi\in\mathcal{D}(\mathbb{R}^{d}).$$

Full text

**Keywords: **inverse problems, identification of inhomogenities, shape calculus, order two methods.

1 Introduction and statement of the results.

Let $\Omega$ be a bounded open set with smooth boundary in $\mathbb{R}^{2}$ or $\mathbb{R}^{3}$. Consider a $L^{\infty}$ function $\sigma$ such that there exists a real $c$ with $\sigma(x)\geq c>0$. Consider the elliptic equation

[TABLE]

with the Dirichlet boundary condition

[TABLE]

Define the Dirichlet-to-Neumann map as

[TABLE]

where $u$ solves (1),(2) and $\mathbf{n}$ is the outer unit normal vector to $\partial\Omega$. The inverse conductivity problem of Calderón is to determine $\sigma$ from $\Lambda_{\sigma}$. Electrical impedance tomography aims to form an image of the conductivity distribution $\sigma$ from the knowledge of $\Lambda_{\sigma}$. When $\sigma$ is smooth enough, one can reconstruct $\sigma$ from $\Lambda_{\sigma}$ (see the works of Sylvester and Uhlmann [21], Nachmann [15, 16] and Novikov [17]). When the conductivity distribution is only ${\rm L}^{\infty}$, Astala and Päivärinta have recently shown in [3] that, in dimension two, the map $\Lambda_{\sigma}$ determines $\sigma\in{\rm L}^{\infty}(\Omega)$.

We are interested in a particular case of that problem: when a body is inserted inside a given object with a distinct conductivity, the question of determining its shape from boundary measurement arises in many fields of modern technology. In the context of the inverse problem of conductivity of Calderón, we restrict the range of admissible conductivity distributions to the family of piecewise constant functions which take only two distinct values $\sigma_{1},\sigma_{2}>0$ which are assumed to be known. The conductivity distribution is then defined by an open subset $\omega$ as

[TABLE]

Here, the only unknown of the problem is $\omega$ a subdomain of $\Omega$ with a smooth boundary $\partial\omega$; its outer unit normal vector is denoted by $\mathbf{n}$. The notation $\chi_{\omega}$ (respectively. $\chi_{\Omega\setminus\omega}$) denotes the characteristic function of $\omega$ (respectively. ${\Omega\setminus\omega}$). The second main difference arises from practical considerations: it is unrealistic from the point of view of applications to know the full graph of Dirichlet-to-Neumann. Therefore, we will assume that one has access to a single point in that graph. This non destructive testing problem is usually written from a numerical point of view as the minimization of a cost function: typically a least-square matching criterion. Many authors have investigated the steepest descent method for this problem [13, 7, 10, 18, 1] with the methods of shape optimization since the unknown parameter is a geometrical domain.

This work is devoted to the study of second order methods for this problem that has only be considered before for simplified models in [5, 2]. By introducing second order methods, one aims to reach two distinct objectives.

• •

On one hand, we provide all the needed material to design a Newton algorithm. We will give differentiability results for the state function and for the objective that we have chosen to study in this work. Nevertheless, we point out that the discretization of a Newton method for this problem turns out to be very delicate; this is why, in the present paper, we will neither discuss about this problem nor present numerical examples. This topic is actually the main objective of a work in progress.

• •

On the other hand, we analyze rigorously the well-posedness of the optimization method. This is justified by the huge numerical literature devoted to the numerical study of this question in the field of inverse problems; the numerical experiments insist on the ill-posedness of this problem. We will explain the instability in the continuous settings in terms of shape optimization. We show that the shape Hessian is not coercive -in fact its Riesz operator is compact – and this explains the unstability of the minimization process.

Let us describe the precise problem under consideration and the notations. We consider a bounded domain $\Omega\subset\mathbb{R}^{d}$ ($d=2$ or $3$) with a $\mathcal{C}^{2}$ boundary. It is filled with a material whose conductivity is $\sigma_{1}$ and with an unknown inclusion $\omega$ in $\Omega$ of conductivity $\sigma_{2}\neq\sigma_{1}$. We search to reconstruct the shape of $\omega$ by measuring on $\partial\Omega$, the input voltage and the corresponding output current. In the sequel, we fix $d_{0}>0$ and consider inclusions $\omega$ such that $\omega\subset\subset\Omega_{d_{0}}=\{x\in\omega,~{}d(x,\partial\Omega)>d_{0}\}$. We also assume that the boundary $\partial\omega$ is of class $\mathcal{C}^{4,\alpha}$. The inverse problem arises when one has access to the normal vector derivative of the potential $u$ that solves (1)-(2) when the conductivity distribution is defined by (3) . Assume that ones knows

[TABLE]

then the problem (1)-(2)-(4) is overdetermined. The electrical impedance tomography problem we consider is to recover the shape of $\omega$ from the knowledge of the single Cauchy pair $(f,g)$.

In order to recover the shape of the inclusion $\omega$, an usual strategy is to minimize a cost function. Many choices are possible; however it turns out that a Kohn and Vogelius type objective leads to a minimization problem with nicer properties than the least squares fitting approaches (we refer to [1] for a comparison of different objectives with order one methods and to [2] for the case of a perfectly insulated inclusion). Therefore, we study such a cost function in this work.

Let us define this criterion. Its distinctive feature is to involve two state functions $u_{d}$ and $u_{n}$: the state $u_{d}$ solves (1)-(2) while $u_{n}$ solves (1)-(4). The Kohn -Vogelius objective $J_{KV}$ is then defined as:

[TABLE]

Let us sum up the results of this paper concerning the minimization of this objective. We first prove differentiability results for the state $u_{d}$. In the sequel, we use the convention that a bold character denotes a vector. If ${\bf{h}}$ denotes a deformation field, it can be written as ${\bf{h}}=\mathbf{h}_{\tau}+h_{n}\mathbf{n}$ on $\partial\omega$. Note also that in the following lines, $\mathbf{n}$ denotes the outer normal field to $\partial\omega$ pointing into $\Omega\setminus\overline{\omega}$. Hence, for $x\in\partial\omega$, we define, when the limit exists, $u^{\pm}(x)$ (resp. $(\partial_{n}u)^{\pm}(x)$) as the limit of $u(x\pm t\mathbf{n}(x))$ (resp. $\langle\nabla u(x\pm t\mathbf{n}(x),\mathbf{n}(x))$) when $t>0$ tends to [math]. Note that $\mathbf{h}_{\tau}$ is a vector while $h_{n}$ is a scalar quantity.

The admissible deformation fields have to preserve $\partial\Omega$ and the regularity of the boundaries: therefore the space of admissible fields is

[TABLE]

The following result concerns the first order derivative of the state functions $u_{d}$ and $u_{n}$. It was derived in [7, 18, 1].

Theorem 1

Let $\Omega$ be an open smooth subset of $\mathbb{R}^{d}$ ($d=2$ or $3$) and let $\omega$ be an element of $\Omega_{d_{0}}$ with a boundary of class $\mathcal{C}^{4,\alpha}$. Then the state functions $u_{d}$ and $u_{n}$ are shape differentiable; furthermore their shape derivative $u_{d}^{\prime}$ and $u_{n}^{\prime}$ belongs to ${\rm H}^{1}(\Omega\setminus\overline{\omega})\cup{\rm H}^{1}(\omega)$ and satisfy

[TABLE]

[TABLE]

The main result of this work concerns the second order derivative. It is given is the following theorem.

Theorem 2

Let $\Omega$ be an open smooth subset of $\mathbb{R}^{d}$ ($d=2$ or $3$) and let $\omega$ be an element of $\Omega_{d_{0}}$ with a $\mathcal{C}^{4,\alpha}$ boundary. Let ${{\bf{h}}}_{1}$ and ${{\bf{h}}}_{2}$ be two deformation fields in $\mathcal{H}$. Then the state $u_{d}$ has a second order shape derivative $u_{d}^{\prime\prime}\in{\rm H}^{1}(\Omega\setminus\overline{\omega})\cup{\rm H}^{1}(\omega)$ that solves

[TABLE]

Here, $(u_{d})^{\prime}_{i}$ denotes the first order derivative of $u$ in the direction of $h_{i}$ as given in (6), $D\mathbf{n}$ stands for the second fundamental form of the manifold $\partial\omega$ and $H$ stands for the mean curvature of $\partial\omega$. The twin result concerning $u_{n}$ is an easy adaption of Theorem 2. Once the differentiability of the state function has been established, one can consider the objectives. In [1], we have shown the first order result.

Theorem 3

Let $\Omega$ be an open smooth subset of $\mathbb{R}^{d}$ ($d=2$ or $3$) and let $\omega$ be an element of $\Omega_{d_{0}}$ with a $\mathcal{C}^{4,\alpha}$ boundary. Let ${{\bf{h}}}_{1}$ and ${{\bf{h}}}_{2}$ be two deformation fields in $\mathcal{H}$. The Kohn-Vogelius objective is differentiable with respect to the shape and its derivative in the direction of a deformation field ${\bf{h}}$ is given by:

[TABLE]

We now give the second-order derivative of the Kohn and Vogelius criterion.

Theorem 4

Let $\Omega$ be an open smooth subset of $\mathbb{R}^{d}$ ($d=2$ or $3$) and $\omega$ be an element of $\Omega_{d_{0}}$ with a $\mathcal{C}^{4,\alpha}$ boundary. Let ${{\bf{h}}}_{1}$ and ${{\bf{h}}}_{2}$ be two deformation fields in $\mathcal{H}$. The Kohn-Vogelius objective is twice differentiable with respect to the shape and its second derivative in the directions ${\bf{h}}_{1}$ and ${\bf{h}}_{2}$ is given by:

[TABLE]

where we have set $v=u_{d}-u_{n}$.

To investigate the properties of stability of this cost function, we are led to consider an admissible inclusion $\omega^{*}$ to solve both (1)-(2) and (1)-(4) in order to obtain the corresponding measurements $f^{*}$ and $g^{*}$. It is obvious that the domain $\omega^{*}$ realizes the absolute minimum of the criterion $J_{KV}$ since, by construction, we can write $u_{d}=u_{n}$ in $\Omega$ and hence $J_{KV}(\omega^{*})=0$. We will check that the Euler equation

[TABLE]

holds. We will also prove that

[TABLE]

Moreover, if $h_{n}\neq 0$, then $D^{2}J_{KV}(\omega^{*})({\bf{h}},{\bf{h}})>0$ holds. Nevertheless, (11) does not mean that the minimization problem is well-posed. In fact, it is the following theorem that explains the instability of standard minimization algorithms.

Theorem 5

Assume that $\omega^{*}$ is a critical shape of $J_{KV}$ for which the additional condition $u_{n}=u_{d}$ holds. Then the Riesz operator corresponding to $D^{2}J_{KV}(\omega^{*})$ defined from ${\rm H}^{1/2}(\partial\omega^{*})$ with values in ${\rm H}^{-1/2}(\partial\omega^{*})$ is compact. Moreover, the minimization problem is severely ill-posed in the following sense: if the target domain is $\mathcal{C}^{\infty}$ and if $\lambda_{n}$ denotes the $n^{th}$ eigenvalue of $D^{2}J_{KV}(\omega^{*})$, then $\lambda_{n}=o(n^{-s})$ for all $s>0$.

Theorem 5 has two main consequences. First, the shape Hessian at the global minimizer is not coercive. This means that this minimizer may not be a local strict minimum of the criterion. Moreover, the criterion provides no control of the distance between the parameter $\omega$ and the target $\omega^{*}$. The second consequence concerns any numerical scheme used to obtain this optimal domain $\omega^{*}$. One has to face this difficulty and this explains why frozen Newton or Levenberg-Marquard schemes have been used to solve numerically this problem [7, 1].

The paper is organized as follows. In a first section, we state some preliminary results. Some are well known facts in shape optimization and will be recalled without proof for the sake of readability. Some of them (e.g the derivatives of a Laplace-Beltrami operator and the tangential regularity of the solution to (1)-(2) along the discontinuity of the conductivity distribution) are less known and will be proved thanks to potential layer methods. Hence we will tackle the computations in Section 3 that we consider as the core of this work : it is essentially devoted to prove Theorem 2. After a first part where we prove the existence of a second order derivative for the state, we propose two distinct methods to find the boundary value problem solved by this second order derivative. The first method (subsection 3.3) follows the lines of classical proofs of shape differentiability by differentiating the weak formulation of problem (1)-(2) and interpreting the result in terms of differential operator and boundary conditions. The alternative method (subsection 3.4) consists in a direct differentiation of the boundary conditions. Finally, Section 4 is devoted to the analysis of the criterion, we establish Theorem 4 and Theorem 5. We will present their consequences on the stability of critical shapes.

2 Preliminary results.

2.1 Elements of shape calculus

Before entering the proof of Theorem 2, we recall without proof some basic facts from shape optimization (see [6] for references). Let ${\bf{h}}$ be a deformation field in $\mathcal{C}^{2}(\Omega,\mathbb{R}^{d})$ with $\|{\bf{h}}\|_{\mathcal{C}^{2}}<1$. We set $T_{t}(h,.)=Id+t{\bf{h}}$ and denote by $\Omega_{t}$ the transported domain $\Omega_{t}=T_{t}(\Omega)$. To avoid heavy notations, we will misuse the notation $T_{t}$ instead of $T_{t}(h,.).$

Material and shape derivatives. Classically, in mechanics of continuous media, the material derivative is defined as being a positive limit. In our context, for any vector field ${\bf{h}}\in\mathcal{H}$, we define the material derivative of the domain functional $y=y(\Omega)$ at $\Omega$ in an admissible direction ${\bf{h}}$ as the limit

[TABLE]

Similarly, one can define the material derivative $\dot{y}(\partial\Omega,{\bf{h}})$ for any domain functional $y=y(\partial\Omega)$ which depends on $\partial\Omega$. Another kind of derivative occurs : it is called the shape derivative of $y(\Omega,{\bf{h}})$. It is viewed as a first local variation. Its definition is given by the following

Definition 1

The shape derivative $y^{\prime}=y^{\prime}(\Omega;{\bf{h}})$ of a functional $y(\Omega)$ at $\Omega$ in the direction of a vector field ${\bf{h}}$ is given by

[TABLE]

For more details on these derivations, the reader can consult [20, 6].

Elements of tangential derivatives. We will need in the sequel to manipulate the tangential differential operators on a manifold. For the reader’s convenience, we recall from [4, 6] some definitions and also some useful rules of calculus.

Definition 2

The tangential divergence of a vector field $\mathbf{V}\in C^{1}(\mathbb{R}^{d},\mathbb{R}^{d})$ is given by

[TABLE]

where the notation $D\mathbf{V}$ denotes the Jacobian matrix of $\mathbf{V}$. When the vector $\mathbf{V}\in C^{1}(\partial\Omega,\mathbb{R}^{d})$ is defined on $\partial\Omega$, then the following notation is used to define the tangential divergence

[TABLE]

where $\tilde{\mathbf{V}}$ stands for an arbitrary $C^{1}$ extension of $\mathbf{V}$ on an open neighborhood of $\partial\Omega$.

We introduce now, the notion of tangential gradient $\nabla_{\tau}$ of any smooth scalar function $f$ in $\mathcal{C}^{1}(\partial\Omega,\mathbb{R}^{d})$.

Definition 3

Let an element $f\in\mathcal{C}^{1}(\partial\Omega,\mathbb{R}^{d})$ be given and let $\tilde{f}$ be an extension of $f$ in the sense that $\tilde{f}\in\mathcal{C}^{1}(U)$ and $\tilde{f}|_{\partial\Omega}=f$ and where $U$ is an open neighborhood of $\partial\Omega$. Then the following notation is used to defined the tangential gradient

[TABLE]

The details for the existence of such an extension can be found in [4]. Let us remark that these definitions do not depend on the choice of the extension. Furthermore, one can show the important relation

[TABLE]

for all elements $f\in\mathcal{C}^{1}(\partial\Omega)$ and all vector fields $F\in\mathcal{C}^{1}(\partial\Omega,\mathbb{R}^{d})$ satisfying $F_{n}=\langle F,n\rangle=0$.

Integration by parts on $\partial\Omega$. In general, the condition above $F_{n}=0$ is not always satisfied. We are then led to find another formula to extend the formula in the general case. The extension of this integration by parts formula to fields with a normal vector component involves curvature.

First, we point out that the curvature is connected to the normal vector via the tangential divergence operator. Recall that the mean curvature of $\partial\Omega$ is defined as $H=\text{div}_{\tau}(\mathbf{n})$. Making use of the form of $\text{div}_{\tau}(\mathbf{n})$ on the boundary, one shows straightforwardly the following statement.

Proposition 1

Let $\Omega$ be an open subset of $\mathbb{R}^{3}$ with a $\mathcal{C}^{2}$ boundary. For any unitary extension ${\cal N}$ of $n$ on a neighborhood of $\partial\Omega$, one has

[TABLE]

Assume that the manifold $\partial\Omega$ has no borders. If $\mathbf{F}\in{\rm H}^{2}(\partial\Omega)^{3}$ and $f\in{\rm H}^{2}(\partial\Omega)$, then we have

[TABLE]

We assume now that the domain $\Omega$ has a $\mathcal{C}^{3}$ boundary. The simplest second-order derivative is the Laplace Beltrami operator; it is defined as follows (see [20, 4, 6]) thanks to the following usual chain rule.

Definition 4

Let $f\in H^{2}(\partial\Omega)$. The Laplace-Beltrami $\Delta_{\tau}$ of $f$ is defined as follows

[TABLE]

There is a relation connecting the Laplace operator and the Laplace-Beltrami operator. Let us denote by $\partial^{2}_{nn}f=(D^{2}f.\mathbf{n}).\mathbf{n}$ where $D^{2}f$ stands for the Hessian of $f$.

Proposition 2

Let $\Omega$ be a domain with a boundary $\partial\Omega$ of class $\mathcal{C}^{3}$. For all functions $f\in H^{3}(\Omega)$, it holds

[TABLE]

We need to compute shape and material derivative of special vector fields: the outer unit normal vector $\mathbf{n}$, the tangential gradient and the Laplace-Beltrami operator applied to a function. While the derivative of the normal vector is obtained by a straightforward calculus, we have to transport from $\partial\Omega_{t}$ to $\partial\Omega$ the Laplace-Beltrami operator and the tangential gradient in order to compute the other derivatives.

Derivatives of the normal vector. We describe the material and shape derivatives of the normal vector. We will denote by $\mathbf{n}$ the gradient of the signed distance to $\partial\Omega$. This is an unitary extension of the unitary normal vector $\mathbf{n}$ at $\partial\Omega$ which is smooth in the vicinity of $\partial\Omega$. This extension furnishes a symmetric Jacobian $D\mathbf{n}$ that satisfies $D\mathbf{n}\,\mathbf{n}=0$ on $\partial\Omega$. The direction ${\bf{h}}$ will be supposed to be in $\mathcal{C}^{2}(\mathbb{R}^{d},\mathbb{R}^{d})$ or in $\mathcal{C}^{2}(\partial\Omega,\mathbb{R}^{d})$.

Proposition 3

The material derivative $\dot{\mathbf{n}}$ of the normal vector $\mathbf{n}$ at $\Omega$ in the direction of a vector field ${\bf{h}}\in\mathcal{C}^{1}(\mathbb{R}^{d},\mathbb{R}^{d})$ is given by

[TABLE]

where $\mathbf{{\bf{h}}}_{\tau}={\bf{h}}-{\bf{h}}.\mathbf{n}~{}\mathbf{n}$.

Concerning its shape derivative defined as $\mathbf{n}^{\prime}=(\partial_{t}\mathbf{n}_{t})|_{t=0}$ where $\mathbf{n}_{t}$ is any smooth unitary extension of $\mathbf{n}$ to $\partial\Omega_{t}$, we obtain.

Proposition 4

The shape boundary $\mathbf{n}^{\prime}$ in the direction of ${\bf{h}}$ is given by

[TABLE]

Derivative of the tangential gradient. For $f\in{\rm H}^{3}(\partial\Omega)$, we compute the material derivative of $\nabla_{\tau}f$. We first compute the difference $\dot{\overline{\nabla_{\tau}f}}-\nabla\dot{f}$.

Proposition 5

For all functions $f\in\mathcal{C}^{2}(\mathbb{R}^{3})$ and directions ${\bf{h}}\in C^{2}(\partial\Omega,\mathbb{R}^{3})$, one has

[TABLE]

• **Proof of Proposition 5. ** We differentiate $\nabla f$ and $\nabla f.\mathbf{n}~{}\mathbf{n}$ and obtain

[TABLE]

while

[TABLE]

The two former equations give the desired result. $\blacksquare$

Derivative of the Laplace-Beltrami operator. Now, we want to compute the material derivative $\dot{\overline{\Delta_{\tau}f}}$. We begin to study how to transport the Laplace-Beltrami operator when one works on $\partial\Omega_{t}$. Let $\Delta_{\tau,t}$ denote the Laplace-Beltrami operator on the manifold $\partial\Omega_{t}$. To compute the derivative of a Laplace-Beltrami operator, we need the following proposition that we quote from [20].

Proposition 6

Let $f\in{\rm H}^{5/2}(\mathbb{R}^{d})$, then

[TABLE]

In the former proposition, we set

[TABLE]

A straightforward computation gives

[TABLE]

Theorem 6

Let $f\in{\cal{D}}(\mathbb{R}^{d})$. The material derivative of $\Delta_{\tau}f$ in the direction ${\bf{h}}$ is given by

[TABLE]

• **Proof of Theorem 6 ** : Formula (24) is shown in a weak sense. For each test function $\phi\in\mathcal{C}^{\infty}(\partial\Omega)$, there exists an extension $\tilde{\phi}\in\mathcal{D}(\mathbb{R}^{d})$ such that $\partial_{\mathbf{n}}\tilde{\phi}=0$; this can be done by extending $\phi$ as a constant along the orbits of the gradient of the signed distance function to $\partial\Omega$ and the use of a cut-off function. For $f\in\mathcal{D}(\mathbb{R}^{d})$, we set

[TABLE]

After an integration by parts on $\partial\Omega$, we obtain:

[TABLE]

Since $\partial_{\mathbf{n}}\tilde{\phi}=0$ and $C(0)=I$, we get

[TABLE]

When $t\rightarrow 0$, it then comes

[TABLE]

Expanding the double divergence term, we obtain:

[TABLE]

In order to explicit these derivatives, we let appear the curvatures of $\partial\Omega$ by means of

[TABLE]

and this ends the proof of the theorem (24). $\blacksquare$

3 Existence of the second order derivative of the state.

Proof of Theorem 2.

The section is devoted to prove Theorem 2. We follow the usual strategy to derive existence in shape optimization. In section 3.2, we will write the weak formulation of the problem, then transport it on the reference domain, pass to the limit and obtain existence of the material derivative. In a second time, we will seek a boundary value problem solved by the material derivative. This will provide a characterization of the second order shape derivative. Two strategies, that we will detail, are possible: the first one explored in section 3.3 consists in working on the variational formulation while the second one uses the tangential differential calculus by differentiating the boundary conditions. This last approach will be presented in section 3.4. The computations that will be made in subsections 3.3 and 3.4 require some regularity of the traces of the state $u_{d}$ on the interface of discontinuity $\partial\omega$. For the sake of readability, we postponed in subsection 3.5 all the needed justifications.

3.1 Preliminary results.

In the sequel, we will use some technical formulae. To preserve the readability of the proof of the main result, we state them in this paragraph. The tools needed for proving these results can be found in [20]. Given a smooth vector field ${\bf{h}}$, we denote

[TABLE]

We begin with the following formula.

Lemma 1

It holds:

[TABLE]

Given two smooth vector fields ${\bf{h}}_{1}$ and ${\bf{h}}_{2}$, we set

[TABLE]

and

[TABLE]

Here, the notation $(A_{{{\bf{h}}}_{1}})^{\prime}({{\bf{h}}}_{2})$ stands for the matrix defined by its elements

[TABLE]

Lemma 2

One has:

[TABLE]

We need the following crucial result

Lemma 3

If $u$ is harmonic then

[TABLE]

• **Proof of Lemma 3 **

For any harmonic function $u$ in $\Omega$ and for every test function $\phi\in\mathcal{D}(\Omega)$, we can write

[TABLE]

then

[TABLE]

Since $\dot{u}=u^{\prime}+{\bf{h}}.\nabla u$ and since $u^{\prime}$ is harmonic in $\Omega$, we obtain the result. $\blacksquare$

3.2 Proof of existence of the second order derivative.

We follow Hettlich and Rundell [8] and Simon [19] to define the second order derivative of an operator with respect to a domain. We compute the second derivative by considering two admissible deformations ${\bf{h}}_{1},{\bf{h}}_{2}\in\mathcal{H}$ that will describe the small variations of $\partial\omega$. Simon shows that the second derivative $F^{\prime\prime}(\partial\omega;{\bf{h}}_{1},{\bf{h}}_{2})$ of $F(\partial\omega)$ is defined as a bounded bilinear operator satisfying

[TABLE]

For more details, the reader can consult the appendix in page 613 of [8].

Let us begin the proof. Let ${{{\bf{h}}}_{1}},{{{\bf{h}}}_{2}}\in\mathcal{H}$ be two vector fields. The direction ${{\bf{h}}}_{1}$ being fixed, we consider $\dot{u}_{1,{{\bf{h}}}_{2}}$ the variation of ${\dot{u_{1}}}$ with respect to the direction ${{\bf{h}}}_{2}$. We recall from [1] that the material derivative ${\dot{u_{1}}}$ of $u$ in the direction ${{\bf{h}}}_{1}$ satisfies

[TABLE]

Let $\phi_{2}:\Omega\mapsto\Omega$ be the diffeomorphism defined by $\phi_{2}(x)=x+{{\bf{h}}}_{2}(x)$ and we set $\psi_{2}=\phi_{2}^{-1}$. Setting $\omega_{{{\bf{h}}}_{2}}=\left\{x+{{\bf{h}}}_{2}(x),~{}x\in\omega\right\}$, $\Omega_{{{\bf{h}}}_{2}}=\left\{x+{{\bf{h}}}_{2}(x),~{}x\in\Omega\right\}=\Omega$ and $\sigma_{{\bf{h}}_{2}}=\sigma\circ\phi_{2}$, we get

[TABLE]

where $u_{{{\bf{h}}}_{2}}$ is the solution of the original problem with $\omega_{{\bf{h}}_{2}}$ instead of $\omega$. Making the change of variables $x=\phi_{2}(X)$, we get the integral identity on the fixed domain $\Omega$ :

[TABLE]

with the notations $\tilde{u}=u\circ\phi_{2}$ and $\widetilde{A_{{{\bf{h}}}_{1}}}=A_{{{\bf{h}}}_{1}}\circ\phi_{2}$. Since the material derivative $\dot{u}_{1}$ of $u$ with respect to the direction ${{\bf{h}}}_{1}$ satisfies

[TABLE]

the difference of (30) and (31) gives

[TABLE]

We quote from [13] and [8] the following asymptotic formulae

[TABLE]

Making the adequate substitutions, we easily check that the material derivative of $\dot{u}_{1}$ with respect to ${{\bf{h}}}_{2}$ exists. This derivative, denoted by $\ddot{u}_{1}$, satisfies

[TABLE]

where $\mathfrak{A}$ is defined in (26).

3.3 Derivation of (8) from the weak formulation.

We want to make explicit the problem solved by $(u^{\prime})^{\prime}$. To achieve this, we should write the right hand side

[TABLE]

as the sum of an integral with $\nabla v$ in factor and an integral of a divergence to identify the jump conditions on $\partial\omega$. To that end, we will use algebraic identities that involve second order derivatives of $u,\dot{u}_{i}$ and of the test function $v\in\mathcal{D}(\Omega)$. Using Lemma 1, we obtain:

[TABLE]

[TABLE]

Concerning the remaining terms, we use Lemma 2 to get

[TABLE]

We apply Lemma 3 and gather the expressions obtained for $F$.

[TABLE]

Using (25), we remove the dependency on $A_{{{\bf{h}}}_{1}}\nabla v$:

[TABLE]

Therefore, we write $F=F_{1}+F_{2}$ where

[TABLE]

The connection between second order material and shape derivatives is given by:

[TABLE]

incorporating this expression in (34), we rewrite (32) as:

[TABLE]

Testing it against $v\in\mathcal{D}(\Omega\setminus\partial\omega)$, we get $\Delta(u^{\prime}_{1})^{\prime}_{2}=0$ in $\Omega\setminus\overline{\omega}$ and in $\omega$. We now deduce the jump conditions for $(u_{1}^{\prime})^{\prime}_{2}$. To obtain the jump of the potential, we simply write that $\ddot{u}_{1}\in{\rm H}^{1}_{0}(\Omega)$, hence $[\ddot{u}_{1}]=0$ on $\partial\omega$ and then

[TABLE]

To express the jump of the flux, we then apply the Gauss formula in (35) to get

[TABLE]

The second term $F_{2}$ contains all the jumps of the flux on the interface $\partial\omega$.

A simplified expression of $F_{2}$. To get a simplified formula for $F_{2}$ under a boundary integral, some lengthy but straightforward calculations are needed. We summarize the result by means of the following lemma

Lemma 4

One has:

[TABLE]

• **Proof of lemma ** First, write :

[TABLE]

Note that the normal vector is oriented from $\omega$ to $\Omega\setminus\overline{\omega}$. In the same spirit, we write

[TABLE]

By a argument of symmetry, we then can write:

[TABLE]

To drop the dependency in $A_{{{\bf{h}}}_{1}}$, we use (25) and get after expansion:

[TABLE]

[TABLE]

After integrating by parts, we conclude thanks to the state equation and obtain

[TABLE]

We substitute the shape derivative $u^{\prime}$ to the material one $\dot{u}$:

[TABLE]

First, we use the continuity of the flux on $\partial\omega$, then we integrate by parts on $\partial\omega$ and finally we incorporate the expressions of the jumps of the shape derivatives $u^{\prime}$ to obtain

[TABLE]

This leads to a simplified expression for $F_{2}$:

[TABLE]

Let us study each term of this sum. Using Gauss formula and integrating by parts on the manifold $\partial\omega$, we obtain

[TABLE]

By symmetry, we also get:

[TABLE]

We now turn to the term with a double divergence. We first write it as a boundary integral thanks to Gauss formula as

[TABLE]

then, we use (14) to introduce the tangential operators

[TABLE]

We study each of these terms. We start with the one involving tangential derivatives: we expand the tangential divergence to incorporate the jump relation for the state $u$.

[TABLE]

Then, the first term becomes:

[TABLE]

We use the integration by parts formula (18) to get:

[TABLE]

Expanding

[TABLE]

we obtain the new expression:

[TABLE]

Now, we consider the term involving normal components. We have

[TABLE]

Then, we get

[TABLE]

A straightforward calculus leads to

[TABLE]

where $D^{2}u$ is the Hessian matrix of $u$. From (20) and from the jump conditions for the state $u$, we deduce that

[TABLE]

When one differentiates the relation expressing the continuity of the flux for the state along the tangential direction $\nabla_{\tau}v$, one gets ([6], p 235):

[TABLE]

In the same spirit, it comes that

[TABLE]

Since $D\mathbf{n}$ is a symmetric matrix and $D\mathbf{n}\,\mathbf{n}=0$, one checks $\nabla v.(D\mathbf{n}[\sigma\nabla_{\tau}u])=[\sigma\nabla u].(D\mathbf{n}\,\nabla_{\tau}v)$. Then

[TABLE]

We integrate this expression on $\partial\omega$ and obtain after some integration by parts:

[TABLE]

Hence

[TABLE]

[TABLE]

Gathering all the terms, we write $F_{2}$ as:

[TABLE]

We end the proof after expanding the tangential divergence of the last term of $F_{2}$. $\blacksquare$

Let us return to the weak formulation (36) of the derivative. By identification, we get

[TABLE]

It remains to compute the jump of the flux for the second order derivative. Since

[TABLE]

where $u^{\prime}_{D{\bf{h}}_{1}\,{\bf{h}}_{2}}$ is the first shape derivative of $u$ in the direction of the vector field $D{\bf{h}}_{1}\,{\bf{h}}_{2}$. Thanks to (6), we can write the jump under the form

[TABLE]

Let us split the field ${\bf{h}}_{2}$ in two parts: $D{\bf{h}}_{1}\,{\bf{h}}_{2}.\mathbf{n}=h_{2,n}\mathbf{n}.D{\bf{h}}_{1}\,\mathbf{n}+D{\bf{h}}_{1}\,\mathbf{h_{2}}_{\tau}.\mathbf{n}$. In the spirit of (40), we obtain

[TABLE]

Thanks to (39), the jump $[\sigma\partial_{\mathbf{n}}u^{\prime}_{D{\bf{h}}_{1}\,{\bf{h}}_{2}}]$ then can be written under the form

[TABLE]

Gathering all the terms, simplifications occur and we get:

[TABLE]

To get the jumps of the potential, we use (41) and obtain

[TABLE]

Thanks to the jump of the potential for the first order shape derivative given in (6), it comes that

[TABLE]

and then:

[TABLE]

Computing the other jumps that appeared in the former expression, we get

[TABLE]

With the help of formula (43), we obtain:

[TABLE]

[TABLE]

[TABLE]

Finally, we gather the results of these computations to write

[TABLE]

3.4 How to recover (8) by formal differentiation of the boundary conditions.

The aim of this section is to retrieve the expression of the flux jump $[\sigma\partial_{\mathbf{n}}u^{\prime\prime}]$ by computing the normal derivatives of each of the expressions $\dot{\overline{[\sigma\nabla u^{\prime}].\mathbf{n}}}$ and $\dot{\overline{\mbox{\rm div}_{\tau}\left(h_{1,n}[\sigma\nabla_{\tau}u]\right)}}$. Since

[TABLE]

then, we get

[TABLE]

In order to avoid lengthy computations, we shall concentrate on each normal derivative appearing in the above formula. Some of the results are straightforward and their proof will be left to the reader. Combining propositions (3) and (5), we conclude that

[TABLE]

In the same manner, we also get

[TABLE]

Hence, we can write

[TABLE]

It remains to simplify the terms $A=(D^{2}u.{\bf{h}}_{2})_{\tau}.\nabla_{\tau}h_{1,n}$ and $B=[\sigma\nabla_{\tau}u].(D^{2}h_{1,n}.{\bf{h}}_{2})_{\tau}$. We obtain:

[TABLE]

We tackle the computation of $(\partial_{\mathbf{n}}u^{\prime})^{\prime}$. For the sake of clearness, we subdivide the work in several steps.

First step. We compute $\dot{\overline{\mbox{\rm div}_{\tau}\left(h_{1,n}[\sigma\nabla_{\tau}u]\right)}}$. We expand:

[TABLE]

Hence, after substitution, one gets

[TABLE]

Second step. We compute $\dot{\overline{[\sigma\partial_{\mathbf{n}}u^{\prime}_{1}]}}$. From the expression of $\dot{\mathbf{n}}$, we get after some straightforward computations:

[TABLE]

Third step. We compute $\sigma\partial_{\mathbf{n}}(u^{\prime}_{1})^{\prime}_{2}$. From the jump condition on the flux of the derivative (6) and (47) and (48), we obtain:

[TABLE]

Taking account of the following calculation,

[TABLE]

it comes

[TABLE]

This formula remains hard to handle. To get a more convenient one, we decide to derive tangentially to the direction ${\bf{h}}_{2}$ the boundary identity

[TABLE]

This leads to:

[TABLE]

From (24) and subtracting (50) from (49), we can write

[TABLE]

From (24), we obtain

[TABLE]

and using the relation between the material and shape derivative, we get

[TABLE]

Injecting these relations in (51) and applying them for $\mathbf{{\bf{h}}_{2}}_{\tau}$, we get

[TABLE]

This last fact allows us to conclude.

3.5 Justification of the formal computations.

We have to justify rigorously that the right-hand sides of (6),(7),(8) make sense. They involve tangential derivatives of $u_{n}$ and $u_{d}$ along the interface $\partial\omega$ up to the order three. The existence of these derivatives is not clear a priori since the gradient of the solution has a discontinuity along this interface. Our first aim is to precise the tangential regularity along the interface $\partial\omega$ of the solution $u$ of (1) with either Dirichlet or Neumann boundary conditions.

We should access to the trace of $u$ on the interface $\partial\omega$. Any numerical discretization needs also to compute the state, its derivatives with respect to the shape and the normal derivatives along the interface $\partial\omega$. To that end, we introduce for any $\alpha\in H^{1/2}(\partial\omega)$ and $\beta\in H^{-1/2}(\partial\omega)$ the following boundary value problems

[TABLE]

where $(f_{1},g_{1})\in H^{1/2}(\partial\Omega)\times H^{-1/2}(\partial\Omega)$. Note that for $\alpha=0$, $\beta=0$ and $(f_{1},g_{1})=(f,g)$ then $(u_{d})$ and $u_{n}$ solve respectively (D) and (N); furthermore the choice of

[TABLE]

leads to (6) and (7) when we take $(f_{1},g_{)}=(0,0).$

Existence of solutions to (D) and (N). To study these problems, we use the integral representation in terms of layer potentials. In a first step, we recall some definitions. The Newtonian potential $\Gamma$ is defined as:

[TABLE]

The integral equations applying to direct problem will be obtained from a study of the classical single- and double-layer potentials. We begin to introduce the following operators

[TABLE]

Note that all these operators have a smooth kernel since the boundaries $\partial\omega$ and $\partial\Omega$ are assumed to have no common point. We also denote

[TABLE]

[TABLE]

We now obtain some systems of integral equations to compute the state function and their shape derivatives. Since $v$ is harmonic in $\Omega\setminus\overline{\omega}$ and for all $x\in\partial\Omega\cup\partial\omega$, it has the classical boundary representation:

[TABLE]

Similarly since $v$ harmonic in $\omega$, for all $x\in\partial\omega$ we can write

[TABLE]

Let us denote by $v_{d}$ the solution of the boundary values problem (D) in (52). Let us show how to compute their restrictions and also their normal vector derivatives on the boundaries. Incorporating the jump conditions, a straightforward computation leads to the following boundary integral equations

[TABLE]

where $\mu=[\sigma]/(\sigma_{1}+\sigma_{2})$. Thanks to (55), the quantity $(\partial_{\mathbf{n}}v_{d})^{+}$ is then given by

[TABLE]

Concerning $v_{n}$, the solution of the Neumann problem (N) in (52), the same kind of computations gives

[TABLE]

Finally, the computation of $(\partial_{\mathbf{n}}v_{n})^{+}_{|\partial\omega}$ is given by

[TABLE]

Concerning the well-posedness of (56), we can state the following result.

Theorem 7

The linear system of integral equation (56) has an unique solution in ${\rm H}^{1/2}(\partial\omega)\times{\rm H}^{-1/2}(\partial\Omega)$.

• **Proof of Theorem 7 **

Let $A$ be the matricial operator defined on ${\rm H}^{{1/2}}(\partial\omega)\times{\rm H}^{-1/2}(\partial\Omega)$ as

[TABLE]

The main argument of the proof is based on the Fredholm alternative. In a first step, we have to show that the adjoint operator $A^{*}$ is injective. Since the boundaries are bounded, the adjoint operator $A^{*}$ defined on ${\rm H}^{{-1/2}}(\partial\omega)\times{\rm H}^{1/2}(\partial\Omega)$ can be written under the form

[TABLE]

Let $(u,v)\in{\rm H}^{{-1/2}}(\partial\omega)\times{\rm H}^{{1/2}}(\partial\Omega)$ be in the kernel of $A^{*}$. Consider the potential $W$ defined for each $x\in\mathbb{R}^{d}$ by

[TABLE]

In a first step, we show that $W=0$. The function $W$ satisfies $\Delta W=0$ in $\mathbb{R}^{d}\setminus(\partial{\omega}\cup\partial\Omega)$ by construction. We check that $W|_{\partial\Omega}=0$ from the equation corresponding to the second line of $A^{*}$. By the properties of the single layer potential, $[W]=0$ on $\partial\omega$. Furthermore, it holds $[\sigma{\partial_{\mathbf{n}}W}]=0$ on $\partial\omega$. Indeed, we can have ([11])

[TABLE]

and

[TABLE]

hence,

[TABLE]

This corresponds to the first line of $A^{*}(u,v)$. Then, $W$ solves the Laplace equation (1) with homogeneous Dirichlet boundary conditions. By the uniqueness of the solution, we get $W=0$ in $\Omega$.

In a second step, we deduce that $u=v=0$. Since $W=0$ in $\Omega$, we see that $[\partial_{\mathbf{n}}W]=0$ on $\partial\omega$. Since $[\partial_{\mathbf{n}}W]=\sigma_{1}u/(\sigma_{1}+\sigma_{2})$ on $\partial\omega$ , we deduce $u=0$. From the second line of $A^{*}(u,v)=0$, we see that $S_{\Omega}v=0$ on $\partial\Omega$. Since the single layer potential operator $S_{\Omega}:{\rm H}^{-1/2}(\partial\Omega)\mapsto{\rm H}^{1/2}(\partial\Omega)$ is an isomorphism, $v=0$ holds. The injectivity of $A^{*}$ is proved. Since $2A=I+C$ where $C$ is a compact operator, we conclude that $A$ has a continuous inverse thanks to the Fredholm alternative. $\blacksquare$

In a similar way, the problem (57) is well-posed under some additional assumptions. We define the adequate space

[TABLE]

We can state the following result.

Theorem 8

If we impose the normalizing condition

[TABLE]

then there exists one unique couple $((v_{n})_{|\partial\omega},(v_{n})_{|\partial\Omega})\in{\rm H}^{1/2}(\partial\omega)\times{\rm H}^{1/2}_{\diamondsuit}(\partial\Omega)$ solution of (57) .

• **Proof of Theorem 8 ** Set

[TABLE]

the operator defined on ${\rm H}^{{1/2}}(\partial\omega)\times{\rm H}_{\diamondsuit}^{1/2}(\partial\Omega)$. The adjoint $B^{*}$ can be written under the form

[TABLE]

In a first step, we begin to show that $B^{*}$ is injective. Let $(u,v)\in{\rm H}^{{1/2}}(\partial\omega)\times{\rm H}^{{1/2}}(\partial\Omega)$ be in the kernel of $B^{*}$. We introduce the potential

[TABLE]

We can see that $Z$ is a harmonic function in ${\mathbb{R}}^{d}\backslash(\partial\omega\cup\partial\Omega)$, satisfying $\partial_{\mathbf{n}}Z~{}~{}=0$ on $\partial\Omega$. By the properties of the single layer potential, $[Z]=0$ Furthermore, a straightforward calculation shows that $[\sigma\partial_{\mathbf{n}}Z]=0$ on $\partial\omega$. Hence, $Z$ solves the boundary value problem

[TABLE]

The function is therefore constant in $\Omega$. Writing $[\partial_{\mathbf{n}}Z]=0$ on $\partial\omega$, we get easily $u=0$ and then $(-\frac{1}{2}+K^{*}_{\Omega})v=0$. Since the operator $\lambda I-K^{*}_{\Omega}$ is one to one on ${\rm H}_{\diamondsuit}^{1/2}(\partial\Omega)$, we deduce that $v=0$. We conclude the proof thanks to the Fredholm alternative. $\blacksquare$

Tangential regularity results. Let us consider now the particular case where both $\alpha$ and $\beta$ are the zero function and $(f_{1},g_{1})=(f,g)$ where $f$ and $g$ are respectively the Dirichlet and Neumann boundary data. To recover the tangential regularity of the solution $u$ along $\partial\omega$, we look at the first line of (56) to deduce that

[TABLE]

and

[TABLE]

It is easy to deduce that $(u_{d})_{|\partial\omega}\in\mathcal{C}^{3,\alpha}(\partial\omega)$. Indeed, from (63) that we consider as an equation in $(u_{d})_{|\partial\omega}$ with data $f$ and $(\partial_{\mathbf{n}}u_{d})_{|\partial\Omega}=g$, we see that $(f,(\partial_{\mathbf{n}}v_{d})_{|\partial\Omega})$ belongs to ${\rm H}^{1/2}(\partial\Omega)\times{\rm H}^{-1/2}(\partial\Omega)$, thanks to Theorem 7.

In order to give a sense to the jump conditions arising in (6),(7),(8), we need to work in space of functions of higher regularity. We choose the framework of Hölder spaces. We quote [12] to precise the behavior of the layer potentials on these spaces.

Theorem 9

(Kirsch [12])

If $\partial\omega$ is of class $C^{2,\alpha},~{}0<\alpha<1$ then the operators $S_{\omega}$ and $K_{\omega}$ map $C^{\beta}(\partial\omega)$ continuously into $C^{1,\beta}$ for all $0<\beta\leq\alpha$. 2. 2.

Let $k\in\mathbb{N}$ with $k\neq 0$. If $\partial\omega$ is of class $C^{k+1,\alpha}$ with $0<\alpha<1$, then the operators $S_{\omega}$ and $K_{\omega}$ map $C^{k,\beta}(\partial\omega)$ continuously into $C^{k+1,\beta}(\partial\omega)$ for all $0<\beta\leq\alpha$. 3. 3.

Let $k$ be an integer. If $\partial\omega$ is of class $C^{k+2,\alpha}$, then $K^{*}_{\omega}$ maps $C^{k,\beta}$ continuously into $C^{k+1,\beta}(\partial\omega)$ for all $0<\beta\leq\alpha$.

We go back to the proof. Since the two boundaries have no intersection point and since $\partial\omega$ is of class $C^{4,\alpha}$, it follows that the right hand side of the former equation is of class $\mathcal{C}^{3,\alpha}(\partial\omega)$. We then conclude the solution of (63) will be of class $\mathcal{C}^{3,\alpha}$ since the operator $1/2I+\mu K_{\omega}$ is an isomorphism from $\mathcal{C}^{3,\alpha}(\partial\omega)$ into itself. With the same arguments, we show straightforwardly that $(\partial_{\mathbf{n}}u_{n})^{+}_{|\partial\omega}\in\mathcal{C}^{2,\alpha}$.

About the regularity of the jumps of the second derivative.  The equations giving the jump conditions $[u_{d}^{\prime}]$ and $[\partial u_{d}^{\prime}]$ show obviously that $[u_{d}^{\prime}]$ and $[\partial_{n}u_{d}^{\prime}]$ belong respectively to $C^{2,\alpha}(\partial\omega)$ and $C^{1,\alpha}(\partial\omega)$. Hence, it comes straightforwardly that $[u_{d}^{\prime\prime}]\in\mathcal{C}^{1,\alpha}$. With the same arguments, we show that $[\partial_{n}u_{d}^{\prime\prime}]\in\mathcal{C}^{0,\alpha}$(see [20] for more details) and then that all the formal computations to get the equations describing the second derivative have a sense.

Remark 1

In a view of a numerical discretization of the state equation, one has to emphasize that the choice of a finite elements method seems inappropriate: one should extract tangential derivative of high order on the interface $\partial\omega$. The obtained numerical accuracy is not sufficient to incorporate the results in an optimization scheme. On the converse, the systems of boundary integral equations (56) and (57) are well-suited for this kind of computation. Nevertheless, a discussion of adapted schemes should be precise and is out of the scope of this manuscript.

3.6 Case of Neumann boundary conditions.

Since the admissible deformation fields have a support with no intersection points with the outer boundary, it is a straightforward application of the preceding computations to show that $u_{n}$ solution to (1)-(4) is twice differentiable with respect to the shape. Furthermore, its second order derivative $u_{n}^{\prime\prime}$ belongs to ${\rm H}^{1}(\Omega\setminus\overline{\omega})\cup{\rm H}^{1}(\omega)$ and solves

[TABLE]

where we use the notations of Theorem 2.

4 Second order derivatives for the criterion.

4.1 Proof of Theorem 4.

The differentiability of the objective is a direct application of Theorem 2. The computation we make here is based on the relation

[TABLE]

To obtain (10), we compute in a first step the shape gradient in the direction ${\bf{h}}_{1}$. Then, in a second step, we differentiate the obtained expression in the direction of ${\bf{h}}_{2}$. In the sequel, we adopt the notation $v=u_{d}-u_{n}$ to obtain concise expressions.

[TABLE]

where

[TABLE]

Now, we use the classical formulae to differentiate a domain integral to get

[TABLE]

The terms $DB_{i},~{}i=1,2$ require more precisions. First, we write

[TABLE]

Note that we used the Green formula twice to keep the symmetry in ${\bf{h}}_{1}$ and ${\bf{h}}_{2}$. We also use the fact that the derivatives $(u_{d})^{\prime}_{i}$ are harmonic in $\Omega\setminus\overline{\omega}$ to transform the boundary integral on the exterior boundary into an integral on the moving boundary. We obtain

[TABLE]

By the same methods, we get

[TABLE]

We regroup the different terms and after some straightforward computations, we obtain:

[TABLE]

In order to compute $D^{2}J_{KV}(\omega)({\bf{h}}_{1},{\bf{h}}_{2})$, the first order derivative of the Kohn-Vogelius objective is needed. It can be written as follows:

[TABLE]

Gathering (66),(41) and (42), we write the second derivative of the Kohn-Vogelius criterion as:

[TABLE]

Let us give a more simplified version for the first term. We decompose the field ${\bf{h}}_{2}$ into normal vector and tangential parts and we use (43). After some elementary computations, we obtain

[TABLE]

Finally, the second order derivative of the Kohn-Vogelius objective becomes:

[TABLE]

4.2 Analysis of stability. Proof of Theorem 5

Now, we specify the domain $\omega$ that is assumed to be a critical shape for $J_{KV}$. Moreover, we assume that the additional condition $u_{d}=u_{n}$ holds. To emphasize that we deal with such a special domain, we will denote it $\omega^{*}$. The assumptions mean that the measurements are compatible and that $\omega^{*}$ is a global minimum of the criterion. From the necessary condition of order two at a minimum, the shape Hessian is positive at such a point.

Let us notice that only the normal component of ${\bf{h}}$ appears. Let us also emphasize that there is no hope to get ${\bf{h}}=0$ from the structure theorem for second order shape derivative ([6]). The deformation field ${\bf{h}}$ appears in $D^{2}J_{KV}(\omega^{*})({\bf{h}},{\bf{h}})$ only thought its normal component $h_{n}$ since $\omega^{*}$ is a critical point for $J_{KV}$. This remark explains why we consider in the statement of Theorem 5 the scalar Sobolev space corresponding to the normal components of the deformation field.

We now prove Theorem 5. From (67), we deduce

[TABLE]

where $\langle,\rangle$ denotes the duality between $H^{1/2}(\partial\omega^{*})\times H^{-1/2}(\partial\omega^{*})$ . Let us introduce the operators

[TABLE]

The Hessian can then be written under the form :

[TABLE]

From the classical results of Maz’ya and Shaposhnikova on multipliers ([14], [22]), we get easily that $T_{1}$ and $T_{2}$ are continuous operators. In fact, the compactness of the Hessian is a consequence of the fact that both operators $M_{1}$ and $M_{2}$ are compact. We use a regularity argument : we remark that $M_{1}$ is the composition of the operators:

[TABLE]

where $\psi$ is the trace on $\partial\omega^{*}$ of $\Psi$ solution of

[TABLE]

While $R_{1}$ is a continuous operator, we prove that $R_{2}$ is compact. Let us express $u_{|\partial\omega^{*}}=\psi$. We use the integral formula of $u$ to obtain:

[TABLE]

The matricial operator arising in this equation appeared also in $\eqref{systemederivee:ud}$. It has a continuous inverse thanks to Theorem 7. Let us express $u|_{\partial\omega^{*}}=\psi$:

[TABLE]

Since the operators $K_{\partial\Omega\partial\omega^{*}}$ and $S_{\partial\Omega\partial\omega^{*}}$ are compact, the operator $R_{2}$ is compact, hence $M_{1}$ is compact. The proof of compactness of $M_{2}$ is similar. Let us mention that a similar strategy of proof can be found in [5].

The natural question is then to quantify how is this optimization problem ill-posed. This question is directly in related to the rate at which the singular values of the Hessian operator are decreasing. Equation (69) shows that this rate is the one of the operators $K_{\partial\Omega\partial\omega^{*}}$ and $S_{\partial\Omega\partial\omega^{*}}$. Now, since for every $u\in{\rm H}^{1/2}(\partial\Omega)$, the functions $K_{\partial\Omega\partial\omega^{*}}u$ and $S_{\partial\Omega\partial\omega^{*}}u$ are harmonic outside of $\partial\Omega$ and therefore in $\Omega$, their restrictions on $\partial\omega^{*}$ are as smooth as $\partial\omega^{*}$. We conclude that if $\partial\omega^{*}$ is $\mathcal{C}^{\infty}$ then the restriction belongs to each ${\rm H}^{s}(\partial\omega^{*})$ for $s>1/2$ then that if $\lambda_{n}$ denotes the $n^{th}$ eigenvalue of $D^{2}J_{KV}(\omega^{*})$, then $\lambda_{n}=o(n^{-s})$ for all $s>0$.