\documentclass[reqno]{amsart}
\usepackage{hyperref}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2016 (2016), No. 110, pp. 1--16.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2016 Texas State University.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2016/110\hfil Asymptotic expansion for shape functions]
{On the high-order topological asymptotic expansion for shape functions}

\author[M. Hassine, K. Khelifi \hfil EJDE-2016/110\hfilneg]
{Maatoug Hassine, Khalifa Khelifi}

\address{Maatoug Hassine \newline
Monastir University,
Department of Mathematics,
Faculty of Sciences, \newline
Avenue de l'Environnement 5000 ,
Monastir, Tunisia}
\email{maatoug.hassine@enit.rnu.tn}

\address{Khalifa Khelifi \newline
Monastir University,
Department of Mathematics, Faculty of Sciences, \newline
Avenue de l'Environnement 5000, Monastir, Tunisia}
\email{khalifakhelifi@hotmail.fr}

\thanks{Submitted January 3, 2016. Published April 26, 2016.}
\subjclass[2010]{35A15, 35B25, 35B40, 49K40}
\keywords{Laplace equation; calculus of variations; sensitivity analysis;
\hfill\break\indent  topological derivative; topology optimization}

\begin{abstract}
 This article concerns the topological sensitivity analysis for the Laplace
 operator with respect to the presence of a Dirichlet geometry perturbation.
 Two main results are presented in this work.
 In the first result we discuss the influence of the considered geometry
 perturbation on  the Laplace solution. In the second result we study the
 high-order topological derivatives. We derive a high-order topological
 asymptotic expansion for a large class of shape functions.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

The topological sensitivity analysis consists in studying the variation
of a shape functional with respect to the presence of a small geometry
perturbation at an arbitrary point of the domain; see
\cite{AH,BHJM,Bend,GGM,GS1,GS2,HJM,HM,SAM}.
To present the basic idea, we consider an open and bounded domain
 $\Omega\subset \mathbb{R}^3$ and a shape function $j(\Omega)=J(u_\Omega)$
to be minimized, where $u_{\Omega}$ is the solution to a given partial
differential equation defined in $\Omega$.
For $\varepsilon >0$, let
 $\Omega_{z,\varepsilon}=\Omega\setminus\overline{\omega_{z,\varepsilon}}$
be the perturbed domain obtained by removing a small part
$\omega_{z,\varepsilon}=z+\varepsilon \omega$ from the domain $\Omega$,
where $z\in\Omega$ and $\omega\subset \mathbb{R}^3$ is a given fixed and
bounded domain containing the origin. The topological sensitivity analysis
leads to an asymptotic expansion of the function $j$ in the form
$$
j(\Omega_{z,\varepsilon})=j(\Omega)+ f(\varepsilon)\delta j(z)+ o(f(\varepsilon)),
$$
where $f(\varepsilon)$ is a scalar positive function approaching zero
as $\varepsilon$ approaches zero.
The function $\delta j$ is called the topological gradient.
It gives us the best locations in $\Omega$ of the geometry perturbations
for which the shape function $j$ decrease most, i.e.
the topological gradient $\delta j$ is as negative as possible.
In fact, if $\delta j(z)<0$, we have $j(\Omega_{z,\varepsilon})<j(\Omega)$
 for small $\varepsilon$.

The topological gradient $\delta j$ has been used as a descent direction
 to solve various problems;  fluid flow optimal shape design
\cite{AH,AHM,BCD}, structural mechanics \cite{HAAS7,GGM}, geometry
inverse problems \cite{AHMas,BHJM,MHKK}, image processing \cite{BJMS},
and many other applications.

The majority of the optimization algorithms dealing with the topological
derivative are based on the first-order asymptotic expansion.
 This provides interesting optimization results in some particular
configurations like the case when the unknown domain is small and not close
to the boundary $\partial\Omega$, one can consult \emph{detection of small
cavities in Stokes flow} in BenAbda et al \cite{BHJM}.

Classically, the topological gradient $\delta j$ described by the leading
term of the first-order asymptotic expansion, dealing only with infinitesimal
geometry perturbations. However, for practical applications, we need
to detect domains of finite size. Therefore, as a natural extension of the
topological derivative concept we consider high-order terms in the asymptotic
expansion. In this context, Novotny et al. \cite{Faria2,JRAARAEC,JRFAA}
was derived a second-order topological asymptotic for the Laplace operator.
The obtained results are limited to the two dimensional case.

In this work, we consider the three dimensional case and we derive a
high-order topological asymptotic expansion for the Laplace operator with
respect to the presence of Dirichlet geometric perturbations.
The proposed approach is based on two main steps.

In the first one, we derive a high-order asymptotic expansion for the
solution of the perturbed Laplace equation with respect to $\varepsilon$.
This question has been investigated by Ammari and Kang \cite{AMK} in the
inhomogeneities case where the perturbed solution is computed in the entire
domain $\Omega$ using continuity condition on the boundary
$\partial \omega_{z,\varepsilon}$. In this work, we deal with more singular geometric
perturbation. The solution of the perturbed Laplace equation is computed in
$\Omega_{z,\varepsilon}=\Omega\setminus\overline{\omega_{z,\varepsilon}}$
with Dirichlet condition on $\partial \omega_{z,\varepsilon}$. As we will show
in Section 3, this type of perturbations leads to an asymptotic behavior
with respect to $\varepsilon$ different from that obtained in \cite{AMK}.

In the second step, we derive a high-order topological asymptotic expansion
for the Laplace operator. More precisely, we derive an asymptotic expansion
of a given shape functional $j$ in the  form
$$
j(\Omega_{z,\varepsilon})=j(\Omega)+\sum_{k=1}^{N}
f_k(\varepsilon)\delta^k j(z)+ o(f_N(\varepsilon)),
$$
where,
\begin{itemize}
\item $f_k$, $1\leq k\leq N$ are positive scalar functions satisfying
$f_{k+1}(\varepsilon)= o(f_k(\varepsilon))$ and
$ \lim_{\varepsilon\to 0 }f_k(\varepsilon)=0$.


\item $\delta^k j$ denotes the $k$-th topological derivative of the shape
function $j$.
\end{itemize}

The topological asymptotic expansion has been derived for various operators
and has been applied for many applications; one can see \cite{GS1}
for the Laplace equation, \cite{GS2,HJM} for the Stokes system, \cite{GGM,HJM}
for the elasticity problem, \cite{JPBS,SAM} for the Helmhotz equation, etc.
In all theses works, the optimization algorithms are based on the first-order
topological derivative which is only valid for small geometry perturbation size.
The use of higher-order terms in the topological asymptotic expansion of
the shape function may certainly be decisive in improving the topological
optimization algorithms without restrictions on the perturbations sizes.
The high-order topological derivative are essential when the first-order
topological derivative $\delta j$ vanishes at some critical points inside $\Omega$.

The present work can be considered as a generalization of the topological
gradient notion. The obtained results are valid for a large class of
 shape functions. The mathematic analysis is general and can be easily
adapted to other partial differential equations.

This article is organized as follows. The formulation of the problem
is presented in Section 2. In Section 3, we discuss the influence of
the geometry perturbation on the Laplace equation solution.
 We derive an asymptotic expansion for the perturbed solution with
respect to $\varepsilon$. The Section 4 is devoted to the high-order topological
derivatives. A high-order topological asymptotic expansion is derived for
a large class of shape functions. Two particular examples of shape functions
are considered in Section 5. Some concluding remarks are presented in Section 6.

\section{Formulation of the problem}

Let $\Omega$ be a bounded domain of $\mathbb{R}^3$ with smooth boundary
$\partial \Omega$. We consider the case in which $\Omega$ contains a
small geometry perturbation $\omega_{z,\varepsilon}$ that is centered at
$z\in \Omega$ and has the form  $\omega_{z,\varepsilon}=z+\varepsilon\omega$,
where $\omega \subset \mathbb{R}^3$ is a given fixed and bounded regular
domain containing the origin.

Consider now a shape function
$$
j(\Omega\setminus\overline{\omega_{z,\varepsilon}})
=J_\varepsilon(u_\varepsilon),
$$
where $J_\varepsilon$ is defined on
$H^1(\Omega\setminus\overline{\omega_{z,\varepsilon}})$
and $u_\varepsilon$ is the solution to Laplace problem in the perturbed
domain $\Omega_{z,\varepsilon}=\Omega\setminus\overline{\omega_{z,\varepsilon}}$
 with homogeneous Dirichlet condition on $\partial \omega_{z,\varepsilon}$
\begin{equation}\label{eqq1}
\begin{gathered}
-\Delta u_\varepsilon =0 \quad \text{in } \Omega_{z,\varepsilon},\\
\nabla u_\varepsilon\cdot n =\Phi_n \quad \text{on } \Gamma_n,\\
u_\varepsilon =\Phi_d \quad \text{on } \Gamma_d, \\
u_\varepsilon =0 \quad\text{on } \partial\omega_{z,\varepsilon},
\end{gathered}
\end{equation}
where $\Phi_n\in H^{-1/2}(\Gamma_n)$ and $\Phi_d\in H^{1/2}(\Gamma_d)$
are two given data, with $\Gamma_n$ and $\Gamma_d$ are two parts of
the boundary $\partial\Omega$ satisfying
$\overline{\partial\Omega}=\overline{\Gamma_n}\cup\overline{\Gamma_d}$
 and $\Gamma_d\cap\Gamma_n=\emptyset$.

Note that for $\varepsilon=0$, we have $\Omega_0=\Omega$ and $u_0$ is
the solution to
\begin{equation}\label{eqq2}
\begin{gathered}
-\Delta u_0 =0 \quad\text{in } \Omega, \\
\nabla u_0\cdot n =\Phi_n \quad \text{on } \Gamma_n, \\
u_0 =\Phi_d \quad \text{on } \Gamma_d.
\end{gathered}
\end{equation}
Using the weak formulation of \eqref{eqq1}, one can deduce that
$u_\varepsilon$ is the unique solution to the variational problem
 find  $u_\varepsilon \in  H^1(\Omega_{z,\varepsilon})$ such that
\begin{equation}\label{eqq2b}
\begin{gathered}
a_\varepsilon(u_\varepsilon,w)=l_\varepsilon(w),\quad
\forall w\in\mathcal{V}_\varepsilon,\\
u_\varepsilon = \Phi_d \quad \text{on } \Gamma_d
\end{gathered}
\end{equation}
where the function space $\mathcal{V}_\varepsilon$,
the bilinear form $a_\varepsilon$, and the linear form $l_\varepsilon$
are defined by:
\begin{gather*}
\mathcal{V}_\varepsilon = \big\{u\in H^1(\Omega_{z,\varepsilon});\,
u=0 \quad \text{on } \Gamma_d\cup\partial{\omega_{z,\varepsilon}} \big\}, \\
 a_\varepsilon(v,w) = \int_{\Omega_{z,\varepsilon}}\nabla v\cdot \nabla w\,dx,
\quad  \forall v,w\in\mathcal{V}_\varepsilon, \\
  l_\varepsilon(w) = \int_{\Gamma_n}\Phi_n w ds,\quad
\forall w\in\mathcal{V}_\varepsilon.
\end{gather*}
In the absence of any perturbation (i.e. $\varepsilon=0$), the weak formulation
of  problem \eqref{eqq2} consists in finding $ u_0\in H^1(\Omega)$ such that
\begin{gather*}
 a_0(u_0,w) = l_0(w),\quad  \forall w\in\mathcal{V}_0\\
 u_0= \Phi_d \quad \text{on } \Gamma_d.
\end{gather*}

As we have mentioned in the introduction, the aim of this work is to derive
a high-order topological asymptotic expansion for the shape function $j$
with respect to the presence of the geometry perturbation $\omega_{z,\varepsilon}$
in the domain $\Omega$. It consists in studying the variation
$j(\Omega_{z,\varepsilon})-j(\Omega)$ with respect to $\varepsilon$ and
establishing an asymptotic formula of the form
$$
j(\Omega_{z,\varepsilon})-j(\Omega)=\sum_{k=1}^{N} f_k(\varepsilon)\delta^k j(z)
+ o(f_N(\varepsilon)).
$$

To derive the expected formula, we will proceed in two steps.
Firstly, we will give a topological sensitivity analysis for the Laplace
operator in Section 3. It consists in studying the asymptotic behavior
 of the solution $u_\varepsilon$ with respect to $\varepsilon$. Secondly, we will study
the variation of a shape function $j$ with respect to the presence of a
geometry perturbation $\omega_{z,\varepsilon}$ in $\Omega$.
The general case, which is valid for a large class of shape functions,
will be discussed in Section 4. In Section 5, we will present the asymptotic
formulas for two shape functionals examples.

\section{Sensitivity analysis for the Laplace operator}\label{asymptotic-formula}

In this section, we give a sensitivity analysis for the Laplace operator
with respect to the presence of a geometry perturbation $\omega_{z,\varepsilon}$
in the domain $\Omega$. More precisely, we derive an asymptotic expansion
for the solution $u_\varepsilon$ with respect to $\varepsilon$. Our procedure is based
on the successive approximations of the variation $u_\varepsilon-u_0$.
 We start our analysis by the following estimate.

\begin{lemma}\label{Lemma0}
Let $\omega_{z,\varepsilon}=z+\varepsilon\omega$ be a topological perturbation
inside the domain $\Omega$. If $\omega_{z,\varepsilon}\subset \Omega$ is not
close to the boundary $\partial\Omega$, then the variation $u_\varepsilon-u_0$
admits the  estimate
 $$
u_\varepsilon(x)-u_0(x)=W_0((x-z)/\varepsilon)\, +\, O(\varepsilon)  \quad\text{in } \Omega_{z,\varepsilon},
$$
where the function $x \mapsto W_0((x-z)/\varepsilon)$ is the unique solution to
the Laplace exterior problem
\begin{equation}\label{w0}
\begin{gathered}
-\Delta W_0 =0 \quad \text{in } \mathbb{R}^3\setminus\overline{\omega},\\
W_0  \to  0 \quad \text{at } \infty\\
W_0  =-u_0(z) \quad \text{on } \partial\omega.
\end{gathered}
\end{equation}
\end{lemma}

\begin{proof}
 The existence of the function $W_0$ is most easily established using a single
layer potential \cite{DL}
$$
W_0(y)=\int_{\partial \omega}G(y-t)\,q_0(t)ds(t), \quad
\forall y\in \mathbb{R}^3\setminus\overline{\omega},
$$
where $G$ is the fundamental solution of the Laplace equation in $\mathbb{R}^3$,
$$
G(y)= \frac{1}{4\pi \|y\|}.
$$
The function $q_0 \in H^{-1/2}(\partial\omega)$ is the solution to the boundary
integral equation
$$
\int_{\partial \omega}G(y-t)\,q_0(t)ds(t)=-u_0(z), \forall y\in \partial \omega.
$$
Posing $R_{0,\varepsilon}(x)= u_\varepsilon(x)-u_0(x)-W_0((x-z)/\varepsilon)$.
One can easily remark that $R_{0,\varepsilon}$ is solution to the system
\begin{gather*}
-\Delta R_{0,\varepsilon} =0 \quad\text{in } \Omega_{z,\varepsilon},\\
\nabla R_{0,\varepsilon}\cdot n =-\nabla W_0((x-z)/\varepsilon)\cdot n \quad
\text{on } \Gamma_n,\\
R_{0,\varepsilon} =-W_0((x-z)/\varepsilon) \quad\text{on } \Gamma_d,\\
R_{0,\varepsilon} =-(u_0-u_0(z)) \quad\text{on }\partial \omega_{z,\varepsilon}.
\end{gather*}
Since the perturbation $\omega_{z,\varepsilon}$ is not close to the boundary
$\partial\Omega$, the function $x \mapsto W_0((x-z)/\varepsilon)$ is
regular in the neighborhood of $\Gamma_d$ and $\Gamma_n$.
It satisfies the following asymptotic behavior:  for all $x\in \Omega_{z,\varepsilon}$,
\begin{align*}
W_0((x-z)/\varepsilon)
&=\varepsilon\int_{\partial \omega}G(x-z-\varepsilon\,t)\,q_0(t)ds(t)\\
&= \varepsilon\,G(x-z)\,\int_{\partial \omega}q_0(t)ds(t) + O(\varepsilon).
\end{align*}
Similarly, the smoothness of $u_0$ near $z$ leads to $u(x)-u_0(z) = O(\varepsilon) $
on $\partial \omega_{z,\varepsilon}$. By elliptic variational inequality,
one can deduce the estimate
$$
R_{0,\varepsilon}= O(\varepsilon)\quad \text{in } \Omega_{z,\varepsilon}.
$$
Consequently, the solution $u_\varepsilon$ of the Laplace equation in the perturbed
domain admits the following asymptotic expansion
$$
u_\varepsilon (x) = u_0(x) + W_0((x-z)/\varepsilon) + O(\varepsilon) \quad\text{in }
\Omega_{z,\varepsilon}.
$$
\end{proof}

This result was proved in \cite[Proposition 3.1]{AH} for the Stokes system.
It has been used to describe the variation of the velocity field with respect
to the presence of a small obstacle.

We are now ready to present the main result of this section.
We will derive a high-order asymptotic expansion of $u_\varepsilon$ with respect to
$\varepsilon$. The obtained result is described by the following theorem.

\begin{theorem}\label{thm-asym}
Let $\omega_{z,\varepsilon}=z+\varepsilon\omega$ be a topological perturbation
inside the domain $\Omega$. If $\omega_{z,\varepsilon}\subset \Omega$ is not close
to the boundary $\partial\Omega$, then the Laplace equation solution $u_\varepsilon$
in the perturbed domain $\Omega_{z,\varepsilon}$ admits the following asymptotic expansion
$$
u_\varepsilon (x)=  \sum_{k=0}^{N} \varepsilon^k [U_k (x)+ W_k ((x-z)/\varepsilon)) ]
+ O(\varepsilon^{N+1}) \quad \text{ in } \Omega_{z,\varepsilon},
$$
where
\begin{itemize}
\item $U_k$, $0\leq k\leq N$ are smooth functions defined in $\Omega$, obtained
as the solutions to a sequence of interior Laplace problems.

\item $W_k$, $0\leq k\leq N$ are smooth functions defined in
$\mathbb{R}^3\setminus\overline{\omega}$, obtained as the solutions to a
sequence of exterior Laplace problems.
\end{itemize}
\end{theorem}

\begin{proof}
 The sequences of functions $(U_k)_{0\leq k\leq N}$ and $(W_k)_{0\leq k\leq N}$
are constructed using an iterative process with $U_0=u_0$ and $W_0$ is
the solution to \eqref{w0}. As we will show later, for all $1\leq k \leq N$:

$\bullet$ The term $U_k$ will be defined as the solution of the Laplace
equation in  $\Omega$ with boundaries conditions depending on the function
$x\mapsto W_{l}((x-z)/\varepsilon)$, $0\leq l\leq k-1$.

$\bullet$ The term $W_k$, will be defined as the solution of the Laplace
equation in  $\mathbb{R}^3\setminus\overline{\omega}$ with a boundary
condition depending on the functions $U_l$, $0\leq l\leq k$.

 Using a single layer potential \cite{DL}, the functions
$W_k$, $ 0\leq k\leq N$ can be written on the following general form
$$
W_k(y)=\int_{\partial \omega}G(y-t)\,q_k(t)ds(t), \quad \forall
y\in \mathbb{R}^3\setminus\overline{\omega},
$$
where $q_k$ is the solution to a boundary integral equation defined on
$\partial \omega$.

To present our construction process, we start our analysis by studying
the variation of the function $x\mapsto W_k((x-z)/\varepsilon)$ with respect
to $\varepsilon$. For each $x\in \mathbb{R}^3\setminus\overline{\omega_{z,\varepsilon}}$
we have
\begin{align*}
W_k((x-z)/\varepsilon)
&= \int_{\partial \omega}G((x-z)/\varepsilon-t)\,q_k(t)ds(t) \\
&= \varepsilon \int_{\partial \omega}G((x-z)-\varepsilon \,t)\,q_k(t)ds(t).
\end{align*}
Using the fact that the perturbation $\omega_{z,\varepsilon}$ is not close to
the boundary $\partial\Omega$, one can remark that for all $t\in \partial\omega$
and for all $x\in \Omega_{z,\varepsilon}$, the function
$\varphi_{x-z,t}:\varepsilon \mapsto \varphi_{x-z,t}(\varepsilon)=\varepsilon\,G((x-z)-\varepsilon\,t)$
is smooth with respect to $\varepsilon$ and satisfies
$$
 \varphi_{x-z,t}(\varepsilon)=\sum_{p=1}^{N}\frac{\varepsilon^p}{p!}\varphi^{(p)}_{x-z,t}(0)
+ O(\varepsilon^{N+1}),
$$
where $\varphi^{(p)}_{x-z,t}(0)$ is the $p$-th derivative of
$\varphi_{x-z,t}$ at $\varepsilon=0$. It depends on the $p$-th derivative
of the function $G$ at the point $x-z$.

 Consequently, the function $x\mapsto W_k((x-z)/\varepsilon)$ admits the
 asymptotic expansion
\begin{equation}\label{fun-wkp0}
W_{k}((x-z)/\varepsilon)=\sum_{p=1}^{N}\varepsilon^p\,W_k^{(p)}(x-z)
+O(\varepsilon^{N+1}),
\end{equation}
where $W_k^{(p)}$ is the smooth function defined in $\mathbb{R}^3\setminus\{z\}$ by
\begin{equation}\label{fun-wkp}
W_k^{(p)}(x-z)=\frac{1}{ p!}
\int_{\partial\omega}\varphi^{(p)}_{x-z,t}(0)q_k(t)ds(t),\quad
\forall x\in \mathbb{R}^3\setminus\{z\}.
\end{equation}
We are now ready to present the main steps of our construction procedure.
\smallskip

\noindent \textbf{First order term:}
It is described by the function $x\mapsto U_1(x) + W_1((x-z)/\varepsilon)$,
$x\in \Omega_{z,\varepsilon}$ which is constructed as follows:
\begin{itemize}
\item The term $U_1$ depends on $W_0$ and solves the interior problem
\begin{equation}
\begin{gathered}
-\Delta U_1 =0 \quad \text{in } \Omega,\\
\nabla U_1\cdot n =-\nabla W_0^{(1)}(x-z)\cdot n \quad\text{on } \Gamma_n,\\
U_1 =-W_0^{(1)}(x-z) \quad \text{on } \Gamma_d,
\end{gathered}
\end{equation}
with $W_0^{(1)}$ is defined by \eqref{fun-wkp} in the particular case
$k=0$ and $p=1$. One can easily check that
$$
W_0^{(1)}(x-z)= G(x-z)\,\int_{\partial \omega}q_0(t)ds(t),
$$
where $q_0$ is the density associated to $W_0$.

\item The term $W_1$ depends on $U_0$ and $U_1$, and solves the following
exterior problem
\begin{equation}\label{eqq6}
\begin{gathered}
-\Delta W_{1} =0 \quad\text{in } \mathbb{R}^3\setminus\overline{\omega},\\
W_1  \to  0\quad  \text{at }  \infty\\
W_1  =-U_1(z)-DU_0(z)(y)  \quad \text{on } \partial\omega.
\end{gathered}
\end{equation}
\end{itemize}

\noindent\textbf{Higher-order terms:}
Let us assume that we have already calculated the first $k-1$ terms.
The $k$-th order term is  described by the function
$x\mapsto U_k(x) + W_k((x-z)/\varepsilon)$, $x\in \Omega_{z,\varepsilon}$
which is defined as follows:
\begin{itemize}
\item The term $U_k$ depends on $W_{j}$, $0\leq j \leq k-1$ and solves the
 interior problem
\begin{equation}
\begin{gathered}
-\Delta U_k =0 \quad\text{in } \Omega,\\
\nabla U_k\cdot n =-\sum_{p=1}^{k}\nabla W_{k-p}^{(p)}(x-z)\cdot n \quad
\text{on } \Gamma_n,\\
U_k = -\sum_{p=1}^{k}W_{k-p}^{(p)}(x-z)\quad \text{on } \Gamma_d,
\end{gathered}
\end{equation}
with $W_j^{(p)}$ is defined by \eqref{fun-wkp}.

\item The term $W_k$ depends on $U_j$, $0\leq j \leq k$ and solves the
 exterior problem
\begin{equation}\label{eqq6b}
\begin{gathered}
-\Delta W_{k} =0 \quad\text{in } \mathbb{R}^3\setminus\overline{\omega},\\
W_k  \to  0\quad \text{at }  \infty\\
W_k  =-U_k(z)- \sum_{p=1}^{k}\frac{1}{p!}D^p{U}_{k-p}(z)(y^p) \quad\text{on }
\partial\omega,
\end{gathered}
\end{equation}
where $D^p{U}_{k-p}(z)$ is the $p$-th derivative of the harmonic function
${U}_{k-p}$ at the point $z\in\Omega$ and $y^p=(y,\dots, y)\in (\mathbb{R}^3)^p$.
\end{itemize}

To prove the desired estimate, we introduce the function
$R_{N,\varepsilon}$ defined in $\Omega_{z,\varepsilon}$ by
\begin{align*}
R_{N,\varepsilon}(x)
&= U_0 (x)+ W_0 ((x-z)/\varepsilon)+\varepsilon\,(U_1 (x)+ W_1
((x-z)/\varepsilon))+\dots \\
&\quad + \varepsilon^N (U_N (x)+ W_N ((x-z)/\varepsilon)-u_\varepsilon (x).
\end{align*}
It is easy to see that $R_{N,\varepsilon}$ is harmonic in $\Omega_{z,\varepsilon}$
and satisfies the following boundary conditions:

On $\partial\omega_{z,\varepsilon}$:
\begin{equation}
\begin{aligned}
R_{N,\varepsilon}(x)
&=  U_0 (x)+ W_0 ((x-z)/\varepsilon) +\sum_{k=1}^{N}\varepsilon^k [U_k (x)+ W_k ((x-z)/\varepsilon)]\\
&=  \sum_{k=0}^{N}\varepsilon^k U_k (x) - \sum_{k=0}^{N} \varepsilon^k
 \Big[ \sum_{p=0}^{k}\frac{1}{p!}D^p{U}_{k-p}(z)(((x-z)/\varepsilon)^p)\Big].
\end{aligned}
\end{equation}
Using the multi-linearity of $D^p{U}_{k-p}(z)$, it follows
\begin{align*}
 \sum_{k=1}^{N} \varepsilon^k\Big[
\sum_{p=0}^{k}\frac{1}{p!}D^p{U}_{k-p}(z)(((x-z)/\varepsilon)^p)\Big]
&= \sum_{k=0}^{N}\sum_{p=0}^{k}\frac{\varepsilon^{k-p}}{p!}D^p{U}_{k-p}(z)((x-z)^p)\\
&= \sum_{k=0}^{N}\varepsilon^k \sum_{p=0}^{N-k}\frac{1}{p!}D^p{U}_{k}(z)((x-z)^p).
\end{align*}
Then, one can deduce
\[
 R_{N,\varepsilon}(x)= \sum_{k=0}^{N}\varepsilon^k
\Big[U_k(x)-\sum_{p=0}^{N-k}\frac{1}{p!}D^p{U}_{k}(z)((x-z)^p)\Big].
\]
Due to Taylor's Theorem and the fact that $\|x-z\| = O(\varepsilon)$ on
$\partial \omega_{z,\varepsilon}$, we obtain
\[
R_{N,\varepsilon}(x)=  O(\varepsilon^{N+1}) \quad \text{on } \partial \omega_{z,\varepsilon}.
\]


On $\Gamma_d$:
\begin{align*}
 R_{N,\varepsilon}(x)
&= \sum_{k=0}^{N}\varepsilon^k W_k((x-z)/\varepsilon)
- \sum_{k=1}^{N}\varepsilon^k \big[\sum_{p=1}^{k}W_{k-p}^{(p)}(x-z)\big]\\
&= \sum_{k=0}^{N}\varepsilon^k W_k((x-z)/\varepsilon) - \sum_{k=0}^{N-1}\varepsilon^k
\big[\sum_{p=1}^{N-k}\varepsilon^p W_{k}^{(p)}(x-z)\big].
\end{align*}
The last equality can be rewritten as
\[
 R_{N,\varepsilon}(x)= \varepsilon^N W_N((x-z)/\varepsilon) +
\sum_{k=0}^{N-1}\varepsilon^k \big[W_k((x-z)/\varepsilon) - \sum_{p=1}^{N-k}
\varepsilon^p W_{k}^{(p)}(x-z)\big].
\]
Then, using the asymptotic expansion \eqref{fun-wkp0} we obtain
$$
R_{N,\varepsilon}(x)= O(\varepsilon^{N+1}) \quad \text{on } \Gamma_d.
$$

On $\Gamma_n$: using the same analysis, one can derive
\[
\nabla R_{N,\varepsilon}\cdot n= O(\varepsilon^{N+1}) \quad \text{on } \Gamma_n.
\]
\end{proof}

\section{High-order topological asymptotic expansion}\label{asymp-expansion}

This section is focused on high-order topological derivatives.
It consists in studying the variation of a shape function $j$ with
respect to the topology perturbation of the domain.
The topology perturbation is described by the hole $\omega_{z,\varepsilon}$
created at an arbitrary point $z\in \Omega$ and having the form
$\omega_{z,\varepsilon}=z+\varepsilon\omega$. We derive a high-order
topological asymptotic expansion for a large class of shape functions.
More precisely, the obtained results are valid for all shape function
 $j$ having the form
$$
j(\Omega_{z,\varepsilon})= J_\varepsilon(u_\varepsilon),
$$
with $J_\varepsilon$ is a scalar function defined on
$H^1(\Omega_{z,\varepsilon})$ and satisfying the assumptions:

\begin{itemize}
\item[(A1)] The function  $J_0$ is differentiable with respect to $u$.
\item[(A2)] There exist real numbers  $\delta^1 J(z),\dots,\delta^N J(z) $,
such that for all $\varepsilon>0$,
$$
J(u_\varepsilon)-J_0(u_0)=DJ_0(u_0)(u_\varepsilon-u_0)
+\sum_{k=1}^{N}\varepsilon^k\delta^k J(z) + o(\varepsilon^N).
$$
\end{itemize}
In the last equality, the solution $u_\varepsilon$ is extended by zero
inside the domain $\omega_{z,\varepsilon}$.
Its extension will be denoted by $u_\varepsilon$ throughout the rest of the paper.

Under the considered assumptions, the variation of the shape function $j$ reads
$$
j(\Omega_{z,\varepsilon})-j(\Omega)
=a_0(u_0-u_\varepsilon,v_0)+\sum_{k=1}^{N} \varepsilon^k
\delta^k J(z) +o(\varepsilon^N),
$$
where $v_0\in\mathcal{V}_0$ is the solution to the  adjoint
problem
\begin{equation}\label{adj-pb}
a_0(w,v_0)=-DJ_0(u_0)(w),\quad \forall w\in\mathcal{V}_0.
\end{equation}
Next, we will derive an asymptotic expansion of the term
$a_0(u_0-u_\varepsilon,v_0)$ which can be written as
\begin{align*}
 a_0(u_0-u_\varepsilon,v_0)
&= \int_{\Omega} (\nabla u_0-\nabla u_\varepsilon)\cdot\nabla v_0dx   \\
&= \int_{\omega_{z,\varepsilon}}\nabla u_0\cdot\nabla
v_0dx+\int_{\Omega_{z,\varepsilon}}(\nabla u_0-\nabla
u_\varepsilon)\cdot\nabla v_0dx.
\end{align*}
Using Green formula, it follows that
\begin{equation}\label{e77}
a_0(u_0-u_\varepsilon,v_0)
= \int_{\omega_{z,\varepsilon}}\nabla
u_0\cdot\nabla v_0dx+\int_{\partial \omega_{z,\varepsilon}}\nabla
(u_0-u_\varepsilon)\cdot nv_0ds.
\end{equation}
By Theorem \ref{thm-asym}, we have
\begin{align*}
 \int_{\partial\omega_{z,\varepsilon}}\nabla(u_0-u_\varepsilon)\cdot nv_0ds
&=  -\sum_{k=1}^{N} \varepsilon^k \int_{\partial\omega_{z,\varepsilon}}
\nabla U_k(x)\cdot n(x)\,v_0(x)ds \\
&\quad -   \sum_{k=0}^{N} \varepsilon^k
\int_{\partial\omega_{z,\varepsilon}}\nabla_x W_k ((x-z)/\varepsilon))\cdot
n\,v_0ds + O(\varepsilon^{N+1}).
\end{align*}
Consequently, the term $a_0(u_0-u_\varepsilon,v_0)$ can be decomposed as
\begin{equation}\label{decomp} %\label{e77}
\begin{aligned}
&a_0(u_0-u_\varepsilon,v_0) \\
&= \int_{\omega_{z,\varepsilon}}\nabla u_0\cdot\nabla v_0dx
- \sum_{k=0}^{N} \varepsilon^k \int_{\partial\omega_{z,\varepsilon}}
 \nabla_x W_k ((x-z)/\varepsilon))\cdot n\,v_0ds\\
&\quad - \sum_{k=1}^{N} \varepsilon^k
\int_{\partial\omega_{z,\varepsilon}} \nabla U_k(x)\cdot
n(x)\,v_0(x)ds +  O(\varepsilon^{N+1}).
\end{aligned}
\end{equation}
In the next section, we will derive an estimate for each term on the
right-hand-side of the equality \eqref{decomp}.

\subsection{Preliminary estimates}\label{prel-estim}

The following lemma gives an estimate for the first term.

\begin{lemma}\label{Lemma1}
 The first term on the right-hand-side of the equality \eqref{decomp} admits
the  asymptotic expansion
$$
\int_{\omega_{z,\varepsilon}}\nabla u_0\cdot\nabla v_0dx
=\sum_{k=3}^{N}\varepsilon^k \,\mathcal T_{u_0,v_0}^{1,k-3}(z)+ O(\varepsilon^{N+1}),
$$
where the functions $z\mapsto \mathcal T_{u_0,v_0}^{1,k}(z)$,
$0\leq k\leq N$ are defined in $\Omega$ by
\begin{equation}\nonumber
\mathcal T_{u_0,v_0}^{1,k}(z)=
\sum_{p=0}^{k}\frac{1}{p!(k-p)!} \int_{\omega}  \nabla^{(p+1)} u_0
(z)(y^p)\cdot  \nabla^{(k-p+1)} v_0 (z)(y^{k-p})dy,
\end{equation}
with $y^k=(y,\dots,y)\in (\mathbb{R}^3)^k$ and $\nabla^{(k)}w(z)$
denotes the $k$-th derivative of the function $w$ at the point $z$.
\end{lemma}

\begin{proof}
The proof of this lemma is based on the well known Taylor-Young formula.
Since $u_0$ and $v_0$ are sufficiently regular in $\omega_{z,\varepsilon}$,
we have
\begin{gather*}
 \nabla u_0(z+\varepsilon y)=\nabla u_0(z)
 + \sum_{k=1}^{N-1}\frac{\varepsilon^k}{k!} \nabla^{(k+1)} u_0 (z)(y^k)+O(\varepsilon^{N})\\
 \nabla v_0(z+\varepsilon y)=\nabla v_0(z)
 + \sum_{k=1}^{N-1}\frac{\varepsilon^k}{k!} \nabla^{(k+1)} v_0
(z)(y^k)+O(\varepsilon^{N}).
\end{gather*}
Using the change of variable $x=z+\varepsilon y$, we derive
\begin{align*}
&\int_{\omega_{z,\varepsilon}}\nabla u_0\cdot\nabla v_0dx \\
&= \varepsilon^3 \int_{\omega}\nabla u_0(z+\varepsilon y)\cdot
 \nabla v_0(z+\varepsilon y)dy\\
&=\varepsilon^3 \int_{\omega}
\Big[\sum_{k=0}^{N-1}\frac{\varepsilon^k}{k!} \nabla^{(k+1)} u_0 (z)(y^k)\Big]
\Big[\sum_{k=0}^{N-1}\frac{\varepsilon^k}{k!} \nabla^{(k+1)} v_0 (z)(y^k)\Big]dy
+O(\varepsilon^{N+1}).
\end{align*}
Using the Cauchy product formula, we obtain the desired result
\begin{align*}
&\int_{\omega_{z,\varepsilon}}\nabla u_0\cdot\nabla v_0dx \\
&= \sum_{k=0}^{N-3}\varepsilon^{k+3}
\Big(\sum_{p=0}^{k}\frac{1}{p!(k-p)!} \int_{\omega}  \nabla^{(p+1)} u_0
(z)(y^p)\cdot  \nabla^{(k-p+1)} v_0 (z)(y^{k-p})dy \Big)\\
&\quad +O(\varepsilon^{N+1}).
\end{align*}
\end{proof}


\begin{lemma}\label{Lemma2}
The second term on the right-hand-side of the equality \eqref{decomp}
admits the  asymptotic expansion
\[
 \sum_{k=0}^{N} \varepsilon^k
\int_{\partial\omega_{z,\varepsilon}}\nabla_x W_k
((x-z)/\varepsilon))\cdot n\,v_0ds
= -\sum_{k=1}^{N} \varepsilon^{k}  \mathcal
T_{W,v_0}^{2,k-1}(z) +  O(\varepsilon^{N+1}),
\]
where the functions $z\mapsto \mathcal T_{W,v_0}^{2,k}(z)$,
$0\leq k\leq N$ are defined in $\Omega$ by
\begin{equation}\nonumber
\mathcal T_{W,v_0}^{2,k}(z)=
-\sum_{p=0}^{k}\frac{1}{p!}
\int_{\partial\omega}\nabla_y W_{k-p}(y)\cdot
n(y)[\nabla^{(p)} v_0 (z)(y^{p})]ds(y).
\end{equation}
\end{lemma}

\begin{proof} 
Using the change of variable $x=z+\varepsilon y$, we obtain
\begin{equation}\label{w_k-v0}
\int_{\partial\omega_{z,\varepsilon}} \nabla_x W_k ((x-z)/\varepsilon))\cdot n(x)v_0(x)ds
=\varepsilon \int_{\partial\omega} \nabla_y W_k(y)\cdot n(y)\,v_0(z+\varepsilon y)ds(y).
\end{equation}
Using the fact that $v_0$ is smooth in  a neighborhood of $z$, one can derive
\begin{align*}
 v_0(z+\varepsilon y)
&=  v_0(z) + \sum_{p=1}^{N-1}\frac{\varepsilon^p}{p!} 
 \nabla^{(p)} v_0 (z)(y^p)+O(\varepsilon^{N})\\
&=  \sum_{p=0}^{N-1}\frac{\varepsilon^p}{p!} \nabla^{(p)} v_0 (z)(y^p)
+O(\varepsilon^{N}).
\end{align*}
It leads to the  asymptotic expansion of the term \eqref{w_k-v0},
\begin{align*}
&\int_{\partial\omega_{z,\varepsilon}}
\nabla_x W_k ((x-z)/\varepsilon))\cdot n(x)v_0(x)ds  \\
&=\sum_{p=0}^{N-1}\frac{\varepsilon^{p+1}}{p!}
\int_{\partial\omega}\nabla_y W_k(y)\cdot n(y)[\nabla^{(p)} v_0
(z)(y^{p})]ds(y)+O(\varepsilon^{N+1}).
\end{align*}
Consequently, 
\begin{align*}
&\sum_{k=0}^{N} \varepsilon^k \int_{\partial\omega_{z,\varepsilon}}
\nabla_x W_k ((x-z)/\varepsilon))\cdot n\,v_0\,ds  \\ 
&=\sum_{k=0}^{N} \varepsilon^k\sum_{p=0}^{N-1}\frac{\varepsilon^{p+1}}{p!} 
 \int_{\partial\omega} \nabla_y W_k(y)\cdot n(y)[\nabla^{(p)} v_0 (z)(y^{p})]ds(y)
 +O(\varepsilon^{N+1})\\ 
&=\sum_{k=1}^{N} \varepsilon^{k}\sum_{p=0}^{k-1}\frac{1}{p!}
\int_{\partial\omega}\nabla_y W_{k-p-1}(y)\cdot
n(y)[\nabla^{(p)} v_0 (z)(y^{p})]ds(y)+O(\varepsilon^{N+1}).
\end{align*}
\end{proof}

\begin{lemma}\label{Lemma3} 
The third term on the right-hand-side of the equality \eqref{decomp} 
admits the following  expansion
\[
\sum_{k=1}^{N} \varepsilon^k
\int_{\partial\omega_{z,\varepsilon}} \nabla
U_k(x)\cdot n(x) v_0(x)ds 
=-\sum_{k=3}^{N} \varepsilon^{k} \mathcal
T_{U,v_0}^{3,k-3}(z)  +  O(\varepsilon^{N+1}).
\]
where the functions $z\mapsto \mathcal T_{U,v_0}^{3,k}(z)$,
$0\leq k\leq N$ are defined in $\Omega$ by
\begin{align*}
&\mathcal T_{U,v_0}^{3,k}(z)\\
&=- \sum_{p=0}^{k} \sum_{q=0}^{p}\frac{1}{q!(p-q)!}
\int_{\partial\omega}  [ \nabla^{(q+1)} U_{k-p+1}
(z)(y^q)] \cdot n(y) [\nabla^{(p-q)} v_0 (z)(y^{p-q})]ds(y).
\end{align*}
\end{lemma}

\begin{proof} Using the change of variable $x=z+\varepsilon y$, we obtain
\begin{equation}\label{u_k-v0}
\int_{\partial\omega_{z,\varepsilon}} 
 \nabla U_k(x)\cdot n(x)v_0(x)ds 
=\varepsilon^2 \int_{\partial\omega} \nabla U_k(z+\varepsilon y)\cdot n(z+\varepsilon y) 
v_0(z+\varepsilon y)ds(y).
\end{equation}
From the fact that $v_0$ is smooth in  a neighborhood of $z$, one can derive
\begin{align*}
 v_0(z+\varepsilon y)
&=  v_0(z) + \sum_{p=1}^{N-1}\frac{\varepsilon^p}{p!} \nabla^{(p)} v_0 (z)(y^p)
 +O(\varepsilon^{N})\\
&=  \sum_{p=0}^{N-1}\frac{\varepsilon^p}{p !} \nabla^{(p)} v_0 (z)(y^p)
 +O(\varepsilon^{N}).
\end{align*}
Similarly, $U_k$ is smooth in  a neighborhood of $z$, it  can be estimated as
\[
\nabla U_k(z+\varepsilon y)=\sum_{q=0}^{N-1}\frac{\varepsilon^q}{q!}
 \nabla^{(q+1)} U_k (z)(y^q)+O(\varepsilon^{N}).
\]
Then, it follows that
\begin{align*}
&\int_{\partial\omega_{z,\varepsilon}} \nabla U_k(x)\cdot n(x)v_0(x)ds \\ 
&=\varepsilon^2  \int_{\partial\omega} [\sum_{q=0}^{N-1}\frac{\varepsilon^q}{q!} 
\nabla^{(q+1)} U_k (z)(y^q)] \cdot n(y) 
\big[\sum_{p=0}^{N-1}\frac{\varepsilon^p}{p!} \nabla^{(p)} v_0 (z)
(y^p)\big]ds(y)+O(\varepsilon^{N+1}).
\end{align*}
Using the Cauchy product formula, one can check the following asymptotic 
expansion of the term \eqref{u_k-v0},
\begin{align*}
&\int_{\partial\omega_{z,\varepsilon}}
\nabla U_k(x)\cdot n(x)v_0(x)ds \\ 
& =\sum_{p=0}^{N-2}\varepsilon^{p+2}
\sum_{q=0}^{p}\frac{1}{q!(p-q)!} \\
&\quad\times \int_{\partial\omega} [
\nabla^{(q+1)} U_k (z)(y^q)] \cdot n(y) [\nabla^{(p-q)} v_0
(z)(y^{p-q})]ds(y)+O(\varepsilon^{N+1}).
\end{align*}
Consequently, 
\begin{align*}
& \sum_{k=1}^{N} \varepsilon^k \int_{\partial\omega_{z,\varepsilon}}  
\nabla U_k(x)\cdot n(x)\,v_0(x)ds \\
& =\sum_{k=1}^{N}\sum_{p=0}^{N-2}\varepsilon^{k+p+2}
 \sum_{q=0}^{p}\frac{1}{q!(p-q)!}\\
&\quad\times  \int_{\partial\omega}  
[ \nabla^{(q+1)} U_{k} (z)(y^q)] \cdot n(y)
[\nabla^{(p-q)} v_0 (z)(y^{p-q})]ds(y)+O(\varepsilon^{N+1})\\
&=\sum_{k=3}^{N} \varepsilon^{k} \sum_{p=0}^{k-3}
\sum_{q=0}^{p}\frac{1}{q!(p-q)!} \\
&\quad\times \int_{\partial\omega}  [ \nabla^{(q+1)} U_{k-p-2}
(z)(y^q)] \cdot n(y) [\nabla^{p-q} v_0
(z)(y^{(p-q)})]ds(y)+O(\varepsilon^{N+1}).
\end{align*}
\end{proof}


\subsection{Asymptotic expansion}\label{asymptotic}
We are now ready to present the main results of this section. 
Based on the previous estimates, we derive a high-order topological 
asymptotic expansion for all shape function satisfying the assumptions
(A1) and (A2).


\begin{theorem}\label{thm-top}
  Let $\omega_{z,\varepsilon}=z+\varepsilon\omega$ be a small topological 
perturbation in $\Omega$ and $j$ a shape function of the form 
$j(\Omega_{z,\varepsilon})=J_{\varepsilon}(u_\varepsilon)$. 
If $J_\varepsilon$ satisfies the assumptions {\rm (A1) and (A2)}, then $j$ 
 admits the  asymptotic expansion
$$
j(\Omega_{z,\varepsilon})-j(\Omega)=\sum_{k=1}^{N}\varepsilon^k\delta^k j(z)
+o(\varepsilon^N),
$$
where $\delta^k j$ is the $k$-th topological derivative defined in 
$\Omega$ by
\begin{equation}\nonumber
\delta^k j(z)=\begin{cases}
\mathcal T_{W,v_0}^{2,k-1}(z)+\delta^k J(z) &\text{if } k=1,2\\
\mathcal T_{u_0,v_0}^{1,k-3}(z)+\mathcal T_{W,v_0}^{2,k-1}(z)
 +\mathcal T_{U,v_0}^{3,k-3}(z)+\delta^k J(z) &\text{if } 3\leq k\leq N.
\end{cases}
\end{equation}
\end{theorem}

\begin{proof}
 Using the fact that $j$ satisfies assumptions (A1) and (A2), we have
$$
J_\varepsilon(u_\varepsilon)-J_0(u_0)
=DJ_0(u_0)(u_\varepsilon-u_0)+\sum_{k=1}^{N}\varepsilon^k\delta^k
J(z) + o(\varepsilon^N).
$$ 
Using \eqref{adj-pb}, we derive
$$
DJ_0(u_0)(u_\varepsilon-u_0)=a_0(u_0-u_\varepsilon,v_0),
$$
Using the decomposition \eqref{decomp} and according to  Lemmas
\ref{Lemma1}, \ref{Lemma2} and \ref{Lemma3}, we derive
\begin{align*}
DJ_0(u_0)(u_\varepsilon-u_0)
&= \sum_{k=3}^{N}\varepsilon^k
 \mathcal T_{u_0,v_0}^{1,k-3}(z)+\sum_{k=1}^{N} \varepsilon^{k}
\mathcal T_{W,v_0}^{2,k-1}(z)\\
 &\quad +\,\sum_{k=3}^{N} \varepsilon^{k} \mathcal T_{U,v_0}^{3,k-3}(z)
 + O(\varepsilon^{N+1}).
\end{align*}
By combining the above equalities we obtain the desired result.
\end{proof}

\section{Shape function examples}\label{funct-exples} 

We now discuss the assumptions (A1) and (A2). We present two examples 
of shape functions satisfying the considered assumptions and we calculate
 their variations $\delta^1J$, $\delta^2J$, \dots, and $\delta^N J$.


\subsection{First example} We consider  the linear function
\begin{equation}
 J_{\varepsilon}(u)=\int_{\Omega_{z,\varepsilon}}g\,udx, \quad
\forall u\in H^1(\Omega_{z,\varepsilon}),
\end{equation}
with $g\in H^1(\Omega)$ is a given function.

\begin{proposition}\label{prop1} 
The function $J_\varepsilon$ satisfies the assumptions {\rm (A1)} and {\rm(A2)} with
\[
 D J_0(w)=\int_{\Omega}gwdx,\quad \forall w\in
\mathcal{V}_0, \text{ and for any } 1\leq k\leq N,\quad 
\delta^k J(z)=0 \quad \text{in } \Omega.
\]
Then the associated shape function 
\[
 j(\Omega_{z,\varepsilon})=\int_{\Omega_{z,\varepsilon}}g\,u_\varepsilon dx
\]
 admits the  high-order asymptotic expansion
$$
j(\Omega_{z,\varepsilon})-j(\Omega)=\sum_{k=1}^{N}
\varepsilon^k\delta^k j(z)+o(\varepsilon^N),
$$
where $\delta^k j$ is the $k$-th topological derivative of $j$
defined in $\Omega$ by
\begin{equation}
\delta^k j(z)=\begin{cases}
\mathcal T_{W,v_0}^{2,k-1}(z) &\text{if } k=1,2\\
\mathcal T_{u_0,v_0}^{1,k-3}(z)+\mathcal T_{W,v_0}^{2,k-1}(z)
 +\mathcal T_{U,v_0}^{3,k-3}(z) &\text{if } 3\leq k\leq N.
\end{cases}
\end{equation}
\end{proposition}

\begin{proof} 
The function $J_0$ is differentiable and we have 
$$ 
D J_0(w)=\int_{\Omega}gw\,dx,\quad \forall w\in \mathcal{V}_0.
$$ 
The variation of $j$ is given by
\[
j(\Omega_{z,\varepsilon})-j(\Omega)
= \int_{\Omega_{z,\varepsilon}}gu_\varepsilon \,dx-\int_{\Omega}gu_0\,dx 
= DJ_0(u_0)(u_\varepsilon-u_0).\nonumber
\]
Hence the function $J_\varepsilon$ satisfies the assumptions (A1) and (A2) with
\begin{gather*}
D J_0(w)=\int_{\Omega}gw\,dx\quad \forall w\in \mathcal{V}_0,\\
\delta^k J(z)=0 \quad  \text{for each $1\leq k\leq N$ and all } z\in\Omega.
\end{gather*}
The asymptotic expansion of $j$ follows immediately from Theorem \ref{thm-top}.
\end{proof} 


\subsection{Second example} 
We consider the semi-norm function associated to the $H^1$ Sobolev space
\begin{equation}
J_{\varepsilon}(u)=\int_{\Omega_{z,\varepsilon}}|\nabla
u-\nabla U_d|^2dx, \quad \forall u\in H^1(\Omega_{z,\varepsilon})
\end{equation}
with $U_d\in H^1(\Omega)$ is a given desired (objective) state, 
smooth in a neighborhood of $z$.

\begin{proposition}\label{prop2}
The function $J_{\varepsilon}$ satisfies the assumptions
{\rm (A1)} and {\rm (A2)} with
$$ 
D J_0(w)= 2\int_{\Omega}\nabla(u_0-U_d)\cdot\nabla
wdx,\quad \forall w\in \mathcal{V}_0,
$$ 
where
\[
 \delta^k J(z)= \begin{cases}
\mathcal T_{W,u_0}^{2,k-1}(z) &\text{if } k=1,2\\
\mathcal T_{W,u_0}^{2,k-1}(z)+
\mathcal T_{u_0,u_0}^{1,k-3}(z)+\mathcal T_{U_d,U_d}^{1,k-3}(z)
+\mathcal T_{U,u_0}^{3,k-3}(z)
&\text{if } 3\leq k\leq N.
\end{cases}
\]
\end{proposition}

\begin{proof}
 The function $J_0$ is differentiable and we have
$$ 
DJ_0(u_0)(w)=2\int_{\Omega}[\nabla u_0-\nabla U_d]\cdot\nabla w dx,
$$
and
\begin{align*}
j(\Omega_{z,\varepsilon})-j(\Omega)
&= \int_{\Omega_{z,\varepsilon}}|\nabla u_\varepsilon-\nabla U_d|^2dx
 -\int_{\Omega}|\nabla u_0-\nabla U_d|^2dx   \\
&= DJ_0(u_0)(u_\varepsilon-u_0)+\int_{\omega_{z,\varepsilon}}
  |\nabla u_0|^2dx \\
&\quad -\int_{\omega_{z,\varepsilon}}|\nabla
U_d|^2dx+\int_{\Omega_{z,\varepsilon}}|\nabla u_0-\nabla
u_\varepsilon|^2dx.\nonumber
\end{align*}
Thanks to the regularity of $u_0$ and $U_d$ in $\omega_{z,\varepsilon}$, one obtains
\begin{gather*}
\int_{\omega_{z,\varepsilon}}|\nabla u_0|^2dx 
=\sum_{k=3}^{N}\varepsilon^k \mathcal
T_{u_0,u_0}^{1,k-3}(z)+ O(\varepsilon^{N+1}) , \\
\int_{\omega_{z,\varepsilon}}|\nabla U_d|^2dx 
=\sum_{k=3}^{N}\varepsilon^k \mathcal T_{U_d,U_d}^{1,k-3}(z)+ O(\varepsilon^{N+1}).
\end{gather*}
By the Green formula, it follows that
$$
\int_{\Omega_{z,\varepsilon}}|\nabla u_0-\nabla u_\varepsilon|^2dx
=-\int_{\partial \omega_{z,\varepsilon}}\nabla (u_0-u_\varepsilon)\cdot nu_0ds.
$$
Applying the technique developed in Section \ref{asymp-expansion}, one can derive
$$
\int_{\Omega_{z,\varepsilon}}|\nabla u_0-\nabla u_\varepsilon|^2dx
=\sum_{k=1}^{N} \varepsilon^{k}  \mathcal T_{W,u_0}^{2,k-1}(z)
 +\sum_{k=3}^{N} \varepsilon^{k} \mathcal T_{U,u_0}^{3,k-3}(z)  +  O(\varepsilon^{N+1}).
$$
By combining the above equalities we obtain the desired result.
\end{proof}

\subsection*{Concluding remarks}
Two main results are presented in this paper.

The first result is devoted to a high-order  asymptotic expansion for the 
Laplace equation solution with respect to the presence of a Dirichlet 
geometry perturbation. This question has been investigated by Ammari 
and Kang \cite{AMK} in the inhomogeneities case. Here, we extend this
 result for a more singular case described by a Dirichlet perturbation.

The second result deals with the high-order topological derivatives. 
A high-order topological asymptotic expansion is derived for a large 
class of shape functions. The use of higher-order terms in the topological 
asymptotic expansion of the shape function may certainly be decisive in 
improving the topological optimization algorithms without restrictions 
on the perturbations sizes. The high-order topological derivative are 
essential when the first-order topological derivative $\delta j$ vanishes 
at some critical points inside $\Omega$.

The present work can be considered as a generalization of the topological 
gradient notion. The mathematic analysis is general and can be easily 
adapted to other partial differential equations.


\begin{thebibliography}{00}
\bibitem{AH} M. Abdelwahed, M. Hassine;
\emph{Topological optimization method for a geometric control problem in Stokes flow.}
 Appl. Numer. Math. 59(8), 1823-1838, 2009.

\bibitem{AHM} M. Abdelwahed , M. Hassine, M. Masmoudi;
\emph{Optimal shape design for fluid flow using topological perturbation technique}, 
J. Math. Anal. and Applic., 356, 548-563, 2009.

\bibitem{AMK} H. Ammari, H. Kang;
\emph{High-order terms in the asymptotic expansions of the steady-state voltage 
potentials in the presence of inhomogeneities of small diameter.} SIAM J Math Analy
34(5):115-166.

\bibitem{SA} S. Amstutz;
 \emph{Aspects th\'eoriques et num\'eriques en optimisation de forme topologique.} 
Th\`ese, Institut National des Siences Appliqu\'ees de Toulouse, 2003.

\bibitem{AHMas} S. Amstutz, I. Horchani, M. Masmoudi;
\emph{Crack detection by the topological gradient method}. Control and Cybernetics, 
Vol. 34, No. 1, pp 81-101, 2005.

\bibitem{BCD} M. Badra, F.Caubet, M. Dambrine;
\emph{Detecting an obstacle immersed in a fluid by shape optimization methods.} 
M3AS, 21 (10), 2069-2101, 2011.

\bibitem{BHJM} A. Ben Abda, M. Hassine , M. Jaoua, M. Masmoudi, 
\emph{Topological sensitivity analysis for the location of small cavities 
in Stokes flow}, SIAM J. Contr. Optim, 48 (5), 2871--2900, 2009.

\bibitem{BJMS} L. Belaid, M.Jaoua, M. Masmoudi, L. Siala;
\emph{Image restoration and edge detection by topological asymptotic expansion}, 
C.R.Acad.Sci. Ser.I 342, p. 313-318, 2006.

\bibitem{Bend} M. Bendsoe;
 \emph{Optimal topology design of continuum structure: an introduction}. 
Technical report, Departement of mathematics, Technical University of 
Denmark, DK2800 Lyngby, Denmark, september 1996.

\bibitem{JRFAA} J. Rocha de Faria, A.A. Novotny;
 \emph{On the second order topologial asymptotic expansion.} Structural 
and Multidis-ciplinary Optimization, 39(6): 547-555, 2009.

\bibitem{Faria2} J. Rocha de Faria, A. A. Novotny, R. A. Feijoo, E. Taroco;
 \emph{First- and second-order topological sensitivity analysis for inclusions} 
J. Inverse Problems in Science and Engineering, vol. 17, no. 5,
pp. 665-679, 2009.

\bibitem{JRAARAEC} J. Rocha de Faria, A. A. Novotny, R. A. Feijao, E. Taroco, C.
Padra;
 \emph{Second order topological sensitivity analysis} International Journal 
of Solids and Structures, vol. 44, no. 14, pp. 4958-4977, 2007

\bibitem{DL} R. Dautray, J. Lions;
 \emph{Analyse math\'emathique et calcul num\'erique pour les sciences et 
les techniques.} Collection CEA, Masson, 1987.

\bibitem{HAAS7} H. A. Eschenauer and A. Schumacher, Topology and 
shape optimization procedures using hole positioning criteria �theory
and applications, in Topology optimization in structural mechanics,
CISM Courses and Lectures 374, Springer, Vienna (1997) 13-96.

\bibitem{GGM} S. Garreau, Ph. Guillaume, M. Masmoudi;
\emph{The topological asymptotic for PDE systems: The elastics case},
 SIAM J. contr. Optim. 39(6), 1756-1778, 2001.

\bibitem{GS1} Ph. Guillaume, K. Sid Idris;
 \emph{The topological asymptotic expansion for the Dirichlet Problem.} 
SIAM J. Control. Optim. 41,  1052-1072, 2002.

\bibitem{GS2} Ph. Guillaume, K. Sid Idris;
\emph{Topological sensitivity and shape optimization for the Stokes equations},
 SIAM J. Control Optim. 43 (1) 1-31, 2004.

\bibitem{GH} Ph. Guillaume, M. Hassine;
\emph{Removing holes in topological shape optimization}, ESAIM, Control, 
Optimisation and Calculus of Variations, 14 (1), 2008, 160-191.

\bibitem{HJM} M. Hassine, S. Jan, M. Masmoudi;
 \emph{From differential calculus to $0-1$ topological optimization}, 
SIAM J. Cont. Optim, 45 (6), 1965-1987, 2007.

\bibitem{MHKK} M. Hassine, Kh. Khelifi;
 \emph{Topology optimization method using the Kohn-Vogelius formulation 
and the topological sensitivity analysis}, Proceedings of the Magrebian 
conference of Applied Mathematics (TAMTAM'13), 188-195, 2013.

\bibitem{HM} M. Hassine, M. Masmoudi;
\emph{The topological asymptotic expansion for the quasi-Stokes problem},
 ESAIM COCV J. 10(4) 478-504, 2004.

\bibitem{mas} M. Masmoudi;
 \emph{The topological asymptotic, Computational method for control applications}, 
Ed.H. Kawarada and J.P\'eriaux, International Series GAKUTO, 2002.

\bibitem{JPBS} J. Pommier, B. Samet;
 \emph{The topological asymptotic for the Helmholtz equation with Dirichlet
 condition on the boundary of an arbitrarily shaped hole}, SIAM J. Control
Optim. 43(3), pp. 899-921, 2004.

\bibitem{SAM} B. Samet, S. Amstutz, M. Masmoudi;
 \emph{The topological asymptotic for the Helmholtz equation}, SIAM
J. Control Optim. Vol. 42(5), pp. 1523-1544, 2003.

\bibitem{SZ} J. Sokolowski., A. Zochowski;
 \emph{On the topological derivative in shape optimization,} 
SIAM J. Control Optim. 37(4) 1251-1272, 1999.

\end{thebibliography}

\end{document}