\begin{document}

\title{Boundary perturbations due to the presence of small linear cracks in an elastic body\thanks{\footnotesize This work was
supported by the ERC Advanced Grant Project MULTIMOD--267184 and
NRF grants No. 2009-0090250, 2010-0004091, and 2010-0017532.}}

\author{Habib Ammari\thanks{\footnotesize Department of Mathematics and Applications, Ecole Normale Sup\'erieure,
45 Rue d'Ulm, 75005 Paris, France (habib.ammari@ens.fr).} \and
Hyeonbae Kang\thanks{Department of Mathematics, Inha University, Incheon
402-751, Korea (hbkang@inha.ac.kr, hdlee@inha.ac.kr, jisunlim@inha.ac.kr).}  \and
Hyundae Lee\footnotemark[3]  \and Jisun Lim\footnotemark[3]}

\maketitle

\begin{abstract}
In this paper, Neumann cracks in elastic bodies are
considered. We establish a rigorous asymptotic expansion for the boundary
perturbations of the displacement (and traction) vectors that are due to the presence of a small elastic linear
crack. The formula reveals that the leading order term is $\varepsilon^2$ where $\varepsilon$ is the length of the crack, and the $\varepsilon^3$-term vanishes. We obtain an asymptotic expansion of the elastic potential energy as an immediate consequence of the boundary perturbation formula. The derivation is based on layer potential
techniques. It is expected that the formula would lead to very
effective direct approaches for locating a collection of small
elastic cracks and estimating their sizes and orientations.
\end{abstract}

\noindent {\footnotesize {\bf Mathematics subject classification
(MSC2000)}: 35B30, 74B05}

\noindent {\footnotesize {\bf Keywords}:  elastic crack, expansion
formula, boundary perturbations}

\section{Introduction}

The displacement (or traction) vector can be perturbed due to the presence of a small crack in an elastic
medium. The aim of this paper is to derive an asymptotic formula for
the boundary perturbations of the displacement as the length of the crack tends to zero. The focus is on cracks with homogeneous Neumann boundary
conditions, {\it i.e.},  perfectly insulating cracks. We
consider the linear isotropic elasticity system in two dimensions and assume that the crack is a
line segment of small size. The derivation of the asymptotic formula is based on layer potential techniques.

The paper extends recent asymptotic results that have been used
for an efficient imaging of small defects. In \cite{AKLP} an
electrostatic model, where the crack is perfectly conducting, was considered and an asymptotic expansion of
the boundary perturbations that are due to the presence of a small
linear crack was derived. The asymptotic formula leads us to efficient algorithms to detect cracks using
 boundary measurements \cite{AGKPS, AKLP}. Their resolution and stability of the algorithms with respect to medium
and measurement noises were investigated in \cite{AGJK, AGS}.
There were some work on boundary perturbation due to the presence of
small inclusions in linear elasticity; the effect of small
inclusions on boundary measurements has been studied in
\cite{AK2004, AKNT}. The effect of thin elastic inclusions on
boundary measurements was quantified in \cite{BF,BFV}. Direct
reconstruction algorithms for locating small or thin elastic
defects were developed in \cite{ACI, BFKL, KKL1, KKL2}. We
emphasize that the results of this paper (on cracks) can {\it not}
be obtained as a limiting case of thin inclusions.

The results of this paper reveals that the leading order term of the boundary perturbation
is $\varepsilon^2$ where $\varepsilon$ is the length of the crack, and its intensity is given by the traction
force of the background solution on the crack (see Theorem \ref{thm-neuman}.) We also prove that
the $\varepsilon^3$-order term vanishes. By integrating the boundary perturbation
formula against the given traction, we are able to derive an asymptotic expansion for
the perturbation of the elastic potential energy, which is an improvement over already existing
results \cite{GN, NDTP} (see the discussion at the end of Section 5).

The boundary perturbation formula derived in this paper carries information about the location,
size, and orientation of the crack, and we expect,  as in the electrostatic case, that the formula will provide a
powerful tool to solve the inverse problem of identifying the
cracks in terms of boundary measurements. The implementation of imaging
algorithms based on the present expansion and the analysis of
their resolution and stability will be the subject of a
forthcoming paper.

The paper is organized as follows. In section \ref{sect2}, a
representation formula for the solution of the problem in the
presence of a Neumann crack is derived. Section \ref{sect3} is
devoted to making explicit the hyper-singular character involved in the
representation formula. Using analytical results for the finite
Hilbert transform, we derive in section \ref{sect4} an asymptotic
expansion of the effect of a small Neumann crack on the boundary
values of the solution. Section \ref{sect5} aims to derive the
topological derivative of the elastic potential energy functional.
Appendix contains technical calculation of the double layer potential.

\section{A representation formula} \langlebel{sect2}

Let $\Omega \subset \mathbb{R}^2$ be an open bounded domain, whose
boundary $\partial \Omega$ is of class $\mathcal{C}^{1, \alpha}$ for some $\alpha>0$. We assume that $\Omega$
is a homogeneous isotropic elastic body so that its elasticity
tensor $\mathbb{C}=(C_{ijkl})$ is given by
\begin{equation}
C_{ijkl}= \langlembda \delta_{ij}\delta_{kl} + \mu (\delta_{ik}\delta_{jl}+ \delta_{il}\delta_{jk})
\end{equation}
with the Lam\'{e} coefficients $\langlembda$ and $\mu$  satisfying
$\mu >0$ and $\langlembda + \mu >0$. Let ${\gamma}_{\varepsilon} \subset
\Omega$ be a small straight crack with  size $\varepsilon$, located at
some fixed distance $d_0$ from $\partialartial \Omega$,  {\it i.e.},
$$\mbox{dist}({\gamma}_{\varepsilon},
\partialartial \Omega) \ge d_0.$$
We denote by ${\bf e}^\partialerp$ a unit normal to $\gamma_\varepsilon$.

Let
$$\Psi := \bigg\{ \partialsi: \partialartial_i \partialsi_j + \partialartial_j \partialsi_i = 0, \quad
1\leq i,j\leq 2 \bigg\},$$ or equivalently,
$$\Psi = \mbox{span}\bigg\{\begin{bmatrix} 1 \\
0
\end{bmatrix}, \begin{bmatrix} 0
\\ 1 \end{bmatrix}, \begin{bmatrix} y \\ -x \end{bmatrix}
\bigg\}.$$ Introduce the space
$$
L^2_\Psi(\partialartial \Omega) := \bigg\{ \mathbf{f} \in L^2(\partialartial
\Omega):  \int_{\partialartial \Omega} \mathbf{f} \cdot {\partialsi} \,
d\sigma = 0 \quad \mbox{for all } \partialsi \in \Psi \bigg\}.
$$
Let ${\bf u}_\varepsilon$ be the displacement vector caused by the traction
${\bf g} \in L^2_\Psi(\partialartial \Omega)$ applied on the boundary $\partial\Omega$ in the presence of $\gamma_\varepsilon$.
Then ${\bf u}_\varepsilon$ is the solution to
\begin{equation}\langlebel{prob1}
\begin{cases}
\nabla \cdot \sigma({\bf u}_\varepsilon) = 0 \quad  & \mbox{in}\; \Omega \setminus \bar{\gamma}_{\varepsilon}, \\
\sigma({\bf u}_\varepsilon)\,{\bf n} = {\bf g}   & \mbox{on} \; \partialartial \Omega,\\
\sigma({\bf u}_\varepsilon)\,{\bf e}^\partialerp = 0   & \mbox{on} \; \gamma_{\varepsilon},
\end{cases}
\end{equation}
where ${\bf n}$ is the outward unit  normal to $\partial\Omega$ and
$\sigma({\bf u}_\varepsilon)$ is the stress defined by
\begin{equation}
\sigma({\bf u}_\varepsilon) = \mathbb{C} \nabla^s {\bf u}_\varepsilon := \frac{1}{2} \mathbb{C} (\nabla {\bf u}_\varepsilon +\nabla {\bf u}_\varepsilon^T ) .
\end{equation}
Here, $\nabla^s {\bf u}_\varepsilon=\frac{1}{2}  (\nabla {\bf u}_\varepsilon +\nabla {\bf u}_\varepsilon^T )$ is the
strain tensor and the superscript $T$ denotes the transpose of a
matrix. Note that the functions in $\Psi$ are solutions to the
homogeneous problem (\ref{prob1}) with ${\bf g}=0$. So we impose
orthogonality condition on
 ${\bf u}_\varepsilon$ to guarantee the uniqueness of a solution to \eqnref{prob1}:
\begin{equation}\langlebel{nor_cond}
\int_{\partial\Omega} {\bf u}_\varepsilon \cdot \partialsi \, d\sigma=0 \quad
\mbox{for all } \partialsi \in \Psi .
\end{equation}
Let  ${\bf u}_0$ be the solution in the absence of the crack, {\it
i.e.}, the solution to
\begin{equation}\langlebel{modelu0}
\begin{cases}
\nabla \cdot \sigma({\bf u}_0) = 0 \quad  & \mbox{in}\; \Omega , \\
\sigma({\bf u}_0)\,{\bf n} = {\bf g}   & \mbox{on} \; \partialartial \Omega,
\end{cases}
\end{equation}
with the orthogonality condition: ${\bf u}_0 |_{\partialartial \Omega} \in
L^2_\Psi(\partialartial \Omega)$ (or equivalently, \eqnref{nor_cond}
with ${\bf u}_\varepsilon$ replaced with ${\bf u}_0$).

It is well-known that the solution ${\bf u}_\varepsilon$ to \eqnref{prob1}
belongs to $H^1(\Omega \setminus \gamma_\varepsilon)$. In fact, we have
 \begin{equation}\langlebel{bound}
 \| {\bf u}_\varepsilon \|_{H^1(\Omega \setminus \gamma_\varepsilon)} \le C
 \end{equation}
for some $C$ independent of $\varepsilon$. To see this, we introduce the potential energy functional
 \begin{equation}\langlebel{energy}
 J_\varepsilon[{\bf u}] := -\frac{1}{2} \int_{\Omega \setminus \gamma_\varepsilon} \sigma({\bf u}) : \nabla^s
 {\bf u}.
 \end{equation}
The solution ${\bf u}_\varepsilon$ of (\ref{prob1}) is the maximizer
of $J_\varepsilon$, {\it i.e.},
 \begin{equation}
 J_\varepsilon[{\bf u}_\varepsilon] = \max J_\varepsilon[{\bf u}],
 \end{equation}
where the maximum is taken over all ${\bf u} \in H^1(\Omega \setminus
\gamma_\varepsilon)$ satisfying $\sigma({\bf u}) \,{\bf n} = {\bf g}$ on $\partial\Omega$ and
$\sigma({\bf u})\,{\bf e}^\partialerp = 0$ on $\gamma_\varepsilon$.  Let ${\bf v}$ be a smooth
function with a compact support in $\Omega$ such that
$\sigma({\bf v})\,{\bf e}^\partialerp = - \sigma({\bf u}_0)\,{\bf e}^\partialerp$ on $\gamma_\varepsilon$. We may
choose  ${\bf v}$ so that $J_\varepsilon[{\bf v}]$ is independent of $\varepsilon$.
Since $0 \ge J_\varepsilon[{\bf u}_\varepsilon] \ge J_\varepsilon[{\bf u}_0 + {\bf v}]$, we have
 $$
 \| \nabla^s {\bf u}_\varepsilon \|_{L^2(\Omega \setminus \gamma_\varepsilon)} \le C.
 $$
We then have from the Korn's inequality that there is a constant $C$ independent of $\varepsilon$ such that
 \begin{equation}\langlebel{korn}
 \| {\bf u}_\varepsilon - {\bf u}_0 \|_{H^1(\Omega \setminus \gamma_\varepsilon)} \le C (\| \nabla^s ({\bf u}_\varepsilon - {\bf u}_0) \|_{L^2(\Omega\setminus \gamma_\varepsilon)} + \| {\bf u}_\varepsilon - {\bf u}_0 \|_{H^{1/2}(\partial\Omega)}).
 \end{equation}
Since $\| {\bf u}_\varepsilon - {\bf u}_0 \|_{H^{1/2}(\partial\Omega)}$ is bounded regardless of $\varepsilon$ as we shall show later (Theorem \ref{thm-neuman}), we obtain \eqnref{bound}.

Let
\begin{equation}\langlebel{phi}
\varphi_\varepsilon({\bf x}):={\bf u}_\varepsilon|_{+}({\bf x}) - {\bf u}_\varepsilon|_{-}({\bf x}), \qquad {\bf x} \in \gamma_{\varepsilon},
\end{equation}
where $+\,(\mbox{resp. } -)$ indicates the limit on the crack
$\gamma_{\varepsilon}$ from the given normal direction ${\bf e}^\partialerp$ (resp.
opposite direction), {\it i.e.},
$$ {\bf u}_{\partialm}(x):= \lim_{t\to 0} {\bf u}({\bf x}\partialm t {\bf e}^\partialerp).$$
We sometimes denote $\sigma({\bf u})\,{\bf n}$, the traction on $\partial\Omega$ (or on $\gamma_\varepsilon$),  by $\partial {\bf u} /\partial \nu$, {\it i.e.},
\begin{equation}\langlebel{conormal1}
\frac{\partial {\bf u}}{\partial \nu} := \langlembda (\nabla \cdot {\bf u}) {\bf n} + \mu( \nabla {\bf u} + \nabla {\bf u}^T ) {\bf n} \quad \mbox{on} \; \partial \Omega.
\end{equation}
If ${\Phi}=(\Phi_{ij})_{2\times 2}$ is the Kelvin matrix of the fundamental solutions of Lam\'{e} system, {\it i.e.},
\begin{equation}\langlebel{lamefund}
   \Phi_{ij}({\bf x}):=\frac{A}{2\partiali}\delta_{ij}\log|{\bf x}|- \frac{B}{2\partiali}\frac{x_ix_j}{|{\bf x}|^2},\qquad {\bf x} \neq 0 \in \mathbb{R}^2,
\end{equation}
where
\begin{equation}
A= \frac{\langlembda+3\mu}{2\mu(\langlembda+2\mu)} \quad\mbox{and}\quad
B= \frac{\langlembda+\mu}{2\mu(\langlembda+2\mu)},
\end{equation}
then the solution ${\bf u}_\varepsilon$ to \eqref{prob1} is represented as
\begin{align}\langlebel{uep_integral}
{\bf u}_\varepsilon({\bf x})&=
 \int_{\partial \Omega} \frac{\partial  \Phi}  {\partial \nu_{\bf y}} ({\bf x}-{\bf y}) {\bf u}_\varepsilon({\bf y})\,d\sigma({\bf y}) -\int_{\partial \Omega}  \Phi({\bf x}-{\bf y}) \frac{\partial {\bf u}_\varepsilon} {\partial \nu}({\bf y}) \,d\sigma({\bf y})\nonumber\\
&\quad  -   \int_{\gamma_{\varepsilon}} \frac{\partial \Phi } {\partial \nu_{\bf y}}({\bf x}-{\bf y}) \varphi_\varepsilon({\bf y}) \,d\sigma({\bf y}), \qquad {\bf x} \in \Omega\setminus \gamma_\varepsilon.
\end{align}
The solution ${\bf u}_0$ to \eqref{modelu0} is represented as
 $$
 {\bf u}_0({\bf x})=
 \int_{\partial \Omega} \frac{\partial  \Phi}  {\partial \nu_{\bf y}} ({\bf x}-{\bf y}) {\bf u}_0({\bf y})\,d\sigma({\bf y}) -\int_{\partial \Omega}  \Phi({\bf x}-{\bf y}) \frac{\partial {\bf u}_0} {\partial \nu}({\bf y}) \,d\sigma({\bf y}).
 $$

Let
 \begin{equation}
 {\bf w}_\varepsilon:= {\bf u}_\varepsilon - {\bf u}_0.
 \end{equation}
Since $\frac{\partial {\bf u}_\varepsilon} {\partial \nu}= \frac{\partial {\bf u}_0} {\partial \nu}$ on $\partial\Omega$, by subtracting above two identities, we have
 \begin{equation}\langlebel{intgralform}
 {\bf w}_\varepsilon ({\bf x}) - \mathcal{D}_\Omega[{\bf w}_\varepsilon]({\bf x}) = - \mathcal{D}_\varepsilon [\varphi_\varepsilon]({\bf x}), \qquad {\bf x} \in \Omega,
 \end{equation}
where the double layer potentials $\mathcal{D}_\Omega$ and $\mathcal{D}_\varepsilon$ are
defined by
 \begin{equation}
 \mathcal{D}_\Omega[{\bf w}_\varepsilon]({\bf x}):=  \int_{\partial \Omega} \frac{\partial  \Phi}
  {\partial \nu_{\bf y}} ({\bf x}-{\bf y}) {\bf w}_\varepsilon({\bf y})\,d\sigma({\bf y}), \quad {\bf x} \in \Omega
 \end{equation}
and
 \begin{equation}
 \mathcal{D}_\varepsilon [\varphi_\varepsilon]({\bf x}):= \int_{\gamma_{\varepsilon}} \frac{\partial \Phi } {\partial \nu_{\bf y}}({\bf x}-{\bf y}) \varphi_\varepsilon({\bf y}) \,d\sigma({\bf y}), \quad {\bf x} \in \Omega.
 \end{equation}

Let ${\bf n}_{\bf x}$ denote the outward normal to $\partial\Omega$ at ${\bf x} \in \partial\Omega$ and let
 $$
 \mathcal{D}_\Omega[\varphi] \big|_{-}({\bf x}) = \lim_{t \to 0+} \mathcal{D}_\Omega[\varphi]({\bf x}+t {\bf n}_{\bf x}).
 $$
Then, it is well known (see, for example, \cite{book2}) that
 $$
 \mathcal{D}_\Omega[\varphi] \big|_{-}({\bf x}) = \left(\frac{1}{2} I + \mathcal{K}_\Omega \right) [\varphi]({\bf x}), \quad\mbox{a.e. } {\bf x} \in \partial\Omega,
 $$
where $\mathcal{K}_\Omega$ is the boundary integral operator defined by
 $$
 \mathcal{K}_\Omega[\varphi]({\bf x}) := \mbox{p.v.} \int_{\partial \Omega} \frac{\partial  \Phi}  {\partial \nu_{\bf y}} ({\bf x}-{\bf y}) {\bf w}_\varepsilon({\bf y})\,d\sigma({\bf y}), \quad {\bf x} \in \partial\Omega,
 $$
and $I$ is the identity operator. Here, p.v. stands for the Cauchy principal value. It then follows from \eqnref{intgralform} that
 \begin{equation}
 \left(-\frac{1}{2} I + \mathcal{K}_\Omega \right) [{\bf w}_\varepsilon] ({\bf x}) = \mathcal{D}_\varepsilon [\varphi_\varepsilon]({\bf x}), \quad {\bf x} \in \partial\Omega.
 \end{equation}
Since $-\frac{1}{2} I + \mathcal{K}_\Omega$ is invertible on
$L^2_\Psi(\partial\Omega)$ (see, for instance, \cite{DKV88}), we have
 \begin{equation}\langlebel{14}
 {\bf w}_\varepsilon({\bf x})= \int_{\gamma_{\varepsilon}} \frac{\partial} {\partial \nu_{\bf y}} (-\frac{1}{2} I + \mathcal{K}_\Omega)^{-1}
 \left[\Phi(\cdot-{\bf y})\right]({\bf x}) \varphi_\varepsilon({\bf y}) \,d\sigma({\bf y}), \quad {\bf x} \in \partial \Omega.
 \end{equation}
Note that
 \begin{equation}\langlebel{fundNeumann}
 (-\frac{1}{2} I + \mathcal{K}_\Omega)^{-1} \left[\Phi(\cdot-{\bf y})\right]({\bf x}) = {\bf N}({\bf x}, {\bf y}), \quad {\bf x} \in \partial\Omega, \ \ {\bf y} \in \Omega,
 \end{equation}
modulo a function in $\Psi$, where
${\bf N}({\bf x}, {\bf y})$ is the Neumann function for the Lam\'{e} system on
$\Omega$, namely, for ${\bf y} \in \Omega$, ${\bf N}({\bf x}, {\bf y})$ is the solution
to
\begin{equation}
\begin{cases}
 \nabla \cdot \sigma({\bf N}(\cdot,{\bf y})) = - \delta_{{\bf y}} {\bf I} \quad & \mbox{in}\; \Omega, \\
 \noalign{
}
 \displaystyle \sigma({\bf N}(\cdot, {\bf y}))\,{\bf n} = - \frac{1}{|\partial \Omega|} {\bf I}    \quad & \mbox{on} \; \partial \Omega,
\end{cases}  
\langlebel{neumann_def}
\end{equation}
subject to the orthogonality condition:
$$\int_{\partial \Omega} {\bf N}({\bf x}, {\bf y}) \cdot \partialsi({\bf x}) \, d\sigma({\bf x}) = 0 \quad \mbox{for all }
\partialsi \in \Psi. $$
Here, ${\bf I}$ is the $2 \time 2$ identity matrix. See \cite{book2, AKL2009,Kup} for properties of the Neumann
function and a proof of (\ref{fundNeumann}). Thus we obtain from
\eqnref{14} that
 \begin{equation}\langlebel{Neumann_repr}
 {\bf u}_\varepsilon ({\bf x}) = {\bf u}_0({\bf x})+ \int_{\gamma_{\varepsilon}} \frac{\partial} {\partial \nu_{\bf y}} {\bf N}({\bf x}, {\bf y}) \varphi_\varepsilon({\bf y}) \,d\sigma({\bf y}), \quad {\bf x} \in \partial \Omega.
 \end{equation}

We now describe the scheme to derive an asymptotic expansion of ${\bf u}_\epsilon-{\bf u}_0$ on $\partial\Omega$. Since $\partiald{{\bf u}_\epsilon}{\nu} =\sigma({\bf u}_\epsilon){\bf e}^\partialerp=0$ on $\gamma_\varepsilon$, we use \eqnref{intgralform} to obtain
 \begin{equation}\langlebel{inteqn1}
 \partiald{{\bf u}_0}{\nu} + \partiald{}{\nu}\mathcal{D}_\Omega[{\bf w}_\varepsilon] = \partiald{}{\nu} \mathcal{D}_\varepsilon [\varphi_\varepsilon] \quad\mbox{on } \gamma_\epsilon.
 \end{equation}
We solve this integral equation for $\varphi_\varepsilon$ and then
substitute it into \eqnref{Neumann_repr} to derive an asymptotic
expansion of ${\bf u}_\varepsilon$ as $\varepsilon \to 0$.

\section{Derivation of an explicit integral equation} \langlebel{sect3}

In view of \eqnref{inteqn1}, we need to compute $\frac{\partial  }{\partial \nu_{\bf x}}\left(\frac{\partial \Phi} {\partial \nu_{\bf y}}({\bf x}-{\bf y}) \right)$ on $\gamma_\varepsilon$. As before ${\bf e}^\partialerp$ is the unit normal to $\gamma_\varepsilon$, and denoted by ${\bf e}^\partialerp=(n_1, n_2)$. It is worth mentioning that $n_1$ and $n_2$ are constant since $\gamma_\varepsilon$ is a line segment. It is convenient to use the following expression of the conormal derivative:
 \begin{equation}
 \frac{\partial {\bf u}}{\partial \nu} =  T(\partial){\bf u} ,
 \end{equation}
where the operator $T(\partial) =T(\partial_1,\partial_2)$, where $\partial_j=\partiald{}{x_j}$, is defined by
\begin{equation}\langlebel{op-T}
T(\xi_1,\xi_2):=\begin{bmatrix}
 (\langlembda+2\mu) n_1 \xi_1 + \mu n_2 \xi_2  &
 \mu n_2 \xi_1 + \langlembda n_1 \xi_2   
 \\
   \langlembda n_2 \xi_1 + \mu n_1 \xi_2 &
    \mu n_1 \xi_1 + (\langlembda+2\mu) n_2 \xi_2
\end{bmatrix}.
\end{equation}

We first obtain the following formula whose derivation will be
given in Appendix \ref{appendix}. For ${\bf x} \neq {\bf y}$, we have
\begin{align}\langlebel{repDij}
\displaystyle {\bf i}g(\frac{\partial \Phi}{\partial \nu_{\bf y}}({\bf x}-{\bf y}){\bf i}g)_{ij}=&
 \left[ a\,\delta_{ij} +b\frac{(x_i-y_i)(x_j-y_j)}{|{\bf x}-{\bf y}|^2}\right]
 \sum_{l=1}^2 \frac{n_l(x_l-y_l)}{|{\bf x}-{\bf y}|^2} \nonumber \\
&\qquad   -\,a \frac{n_j(x_i-y_i)-n_i(x_j-y_j)}{|{\bf x}-{\bf y}|^2},
  \end{align}
where
 \begin{equation}
 a= - \frac{\mu}{2\partiali( \langlembda+2\mu)}, \quad
 b =  - \frac{( \langlembda+\mu)}{\partiali( \langlembda+2\mu)} .
 \end{equation}
Let $v_{ij}=(\frac{\partial \Phi}{\partial\nu_{\bf y}}({\bf x}-{\bf y}))_{ij}$ for convenience and let
\begin{equation}\langlebel{W-fn}
W({\bf x}-{\bf y}):=\frac{\partial  }{\partial \nu_{\bf x}}\left(\frac{\partial \Phi} {\partial \nu_{\bf y}}({\bf x}-{\bf y})\right).
\end{equation}
Then one can use \eqnref{op-T} to derive
\begin{eqnarray*}
W({\bf x}-{\bf y})_{11}
&=& n_1[ (\langlembda +2\mu)  \partial_1v_{11}
+\langlembda {\partial_2}v_{21}] + n_2 ( \mu{\partial_2}v_{11}+  \mu{\partial_1}v_{21}),\\
W({\bf x}-{\bf y})_{12}
&=&n_1[(\langlembda+2\mu) \partial_1v_{12} + \langlembda{\partial_ {x_2}}v_{22} ]+n_2 ( \mu\partial_2v_{12} +  \mu{\partial_1}v_{22}),\\
W({\bf x}-{\bf y})_{21}
&=& n_1 (\mu\partial_2v_{11}+\mu\partial_1v_{21})
+n_2[ \langlembda \partial_1v_{11}+(\langlembda +2\mu)\partial_2v_{21}],\\
W({\bf x}-{\bf y})_{22}
&=&n_1(\mu {\partial_2}v_{12}+
\mu {\partial_1}v_{22}) +n_2[\langlembda {\partial_1}v_{12}
+(\langlembda+2\mu) {\partial_2}v_{22}].
\end{eqnarray*}

Since  the crack which we consider is a line segment with length $\varepsilon$ in the domain $\Omega \subset \mathbb{R}^2$,
we may assume, after rotation and translation if necessary, that it is given by
\begin{equation}\langlebel{linecrack}
\gamma_\varepsilon = \{ (x_1, 0) :  -\varepsilon/2 \le x_1 \le \varepsilon/2 \}.
\end{equation}
In this case, one can check that
\begin{align*}
\partial_2v_{11}&= \left(a+b\frac{(x_1-y_1)^2}{|{\bf x}-{\bf y}|^2}\right)
\frac{1}{|{\bf x}-{\bf y}|^2}=\frac{a+b}{(x_1-y_1)^2},\\
\partial_1v_{21}&= a\left( \frac{1}{|{\bf x}-{\bf y}|^2 } -2\frac{x_1-y_1}{|{\bf x}-{\bf y}|^2}\frac{x_1-y_1}{|x_1-y_1|^2}\right) = -\frac{a}{(x_1-y_1)^2},\\
\partial_1v_{12}&= -a\left( \frac{1}{|{\bf x}-{\bf y}|^2 } -2\frac{x_1-y_1}{|{\bf x}-{\bf y}|^2}\frac{x_1-y_1}{|x_1-y_1|^2}\right) = \frac{a}{(x_1-y_1)^2},\\
\partial_2v_{22}&= \frac{a}{(x_1-y_1)^2} \\
\partial_1v_{11}&= \partial_2v_{21} = \partial_2v_{12} = \partial_1v_{22} = 0 .
\end{align*}
Since ${\bf e}^\partialerp=(n_1,n_2)=(0,1)$, we have
\begin{eqnarray*}
W({\bf x}-{\bf y})_{11} &=& \mu \partial_2v_{11}+\mu \partial_1v_{21}\\
&=& \mu(a+b) \frac{1}{(x_1-y_1)^2} - \mu a\frac{1}{(x_1-y_1)^2}
=  \frac{\mu b}{(x_1-y_1)^2},\\
W({\bf x}-{\bf y})_{12} &=& \mu \partial_2v_{12}+\mu \partial_1v_{22}
= 0,\\
W({\bf x}-{\bf y})_{21} &=& \langlembda  \partial_1v_{11}+ (\langlembda+2  \mu) \partial_2v_{21}= 0,\\
W({\bf x}-{\bf y})_{22} &=& \langlembda  \partial_1v_{12}+ (\langlembda+2\mu) \partial_2v_{22}\\
 &=& \langlembda  a \frac{1}{(x_1-y_1)^2}+(\langlembda+2\mu)a \frac{1}{(x_1-y_1)^2}
 =  \frac{2(\langlembda+\mu) a}{(x_1-y_1)^2},
\end{eqnarray*}
that is,
\begin{equation}
W({\bf x}-{\bf y})=  \frac{1}{(x_1-y_1)^2}
 \begin{bmatrix}
 -\frac{\mu (\langlembda+\mu)}{\partiali(\langlembda+2\mu)}& 0 \\  0 &  -\frac{\mu(\langlembda+\mu)}{\partiali(\langlembda+2\mu)}
 \end{bmatrix} .
\end{equation}
Note that
\begin{equation}\langlebel{c-lam}
\frac{\mu(\langlembda+\mu)}{\langlembda+2\mu}=\frac{E}{4}
\end{equation}
where $E$ is the Young's modulus in two dimensions. So, we have
\begin{equation}
W({\bf x}-{\bf y})=  -\frac{E}{4\partiali}\frac{1}{(x_1-y_1)^2} {\bf I}.
\end{equation}

So far we have shown that if $\gamma_\varepsilon$ is given by \eqnref{linecrack}, then
\begin{equation}\langlebel{ddb2}
 \frac{\partial } {\partial \nu_{\bf x}}  \mathcal{D}_{\varepsilon}[\varphi_\varepsilon]({\bf x})= - \frac{E}{4\partiali}
 \int_{-\varepsilon/2}^{\varepsilon/2} \frac{\varphi_{\varepsilon}(y)}{(x-y)^2}
   \,   dy,\quad {\bf x}=(x,0), \ -\varepsilon/2<x<\varepsilon/2.
  \end{equation}
Here the integral is hyper-singular and should be understood as a finite part in the sense of Hadamard, which will be defined in the next section. So the integral equation \eqnref{inteqn1} becomes
 \begin{equation}
 \frac{1}{\partiali} \int_{-\varepsilon/2}^{\varepsilon/2} \frac{\varphi_{\varepsilon}(y)}{(x-y)^2} \,dy= - \frac{4}{E} {\bf f}(x), \quad -\varepsilon/2<x<\varepsilon/2,
 \end{equation}
where
 \begin{equation}
 {\bf f}(x)=  \partiald{{\bf u}_0}{\nu}(x,0) + \partiald{}{\nu}\mathcal{D}_\Omega[{\bf w}_\varepsilon](x,0) .
 \end{equation}
Define
\begin{equation}
{\bf f}_\varepsilon(x) := {\bf f} (\frac{\varepsilon}{2} x),
\end{equation}
and
\begin{equation}\langlebel{psiGve}
\partialsi_\varepsilon(x):= \frac{2}{\varepsilon} \varphi_\varepsilon (\frac{\varepsilon}{2} x), \quad -1< x <1.
\end{equation}
Then the scaled integral equation is
\begin{equation}\langlebel{scaleeqn}
 \frac{1}{\partiali} \int_{-1}^{1} \frac{\partialsi_{\varepsilon}(y)}{(x-y)^2}\,dy = - \frac{4}{E} {\bf f}_\varepsilon (x), \quad -1<x<1,
 \end{equation}
which we solve in the next section.

\section{Asymptotic expansion} \langlebel{sect4}

The integral in \eqnref{scaleeqn} is understood as a finite-part in the
 sense of Hadamard \cite{Mart92, Mart96}: for $\partialsi \in \mathcal{C}^{1,\alpha}(-1,1)  ~(0<\alpha\le 1)$
\begin{equation}\langlebel{HM}
\Xint\times_{-1}^{1} \frac{\partialsi(y)}{(x-y)^2}\,dy = \lim_{\delta
\rightarrow 0} {\bf i}g[\int_{-1}^{x-\delta} \frac{\partialsi(y)}{(x-y)^2}dy  +
\int_{x+\delta}^1 \frac{\partialsi(y)}{(x-y)^2}dy -
\frac{2\partialsi(x)}{\delta}{\bf i}g].
\end{equation}
Define
\begin{equation}\langlebel{eq:pv}
\mathcal{A}[\partialsi](x): = \frac{1}{\partiali} \Xint\times_{-1}^{1} \frac{\partialsi(y) }{(x-y)^2}\,dy, \quad |x|<1.
\end{equation}
It is known (\cite{Ioak, Mart89}) that \begin{equation}\langlebel{aandh}
\mathcal{A}[\partialsi](x) =- \frac{d}{dx} \mathcal{H}[\partialsi](x), \end{equation} where $\mathcal{H}$
is the (finite) Hilbert transform, {\it i.e.}, \begin{equation} \mathcal{H}[\partialsi](x)
= \mbox{p.v.} \frac{1}{\partiali}\int_{-1}^{1}  \frac{\partialsi(y)}{x-y}\,dy.
\end{equation}
More properties of finite-part integrals and principal-value integrals
can be found in \cite{Kup, Mart92, Mart96, Musk,Tric, Zabr}.

If $\partialsi(-1)=\partialsi(1)=0$, we have from \eqnref{aandh} that
\begin{equation}\langlebel{hilbt}
\mathcal{A}[\partialsi](x)= -\mathcal{H}[\partialsi'](x).
\end{equation}
Thus we can invert the operator $\mathcal{A}$ using the properties of $\mathcal{H}$.
The set $\mathcal{Y}$, given by
\begin{equation} \langlebel{xeps}
  \mathcal{Y} \,=\,
\displaystyle \bigg\{ \varphi: \int_{-1}^1 \sqrt{1 - x^2}\,
|\varphi(x)|^2 \, dx < +\infty
  \bigg\},
  \end{equation}
 is a Hilbert space with the norm
$$ || \varphi ||_{\mathcal{Y}} = \bigg( \int_{-1}^1 \sqrt{1 - x^2}\,
|\varphi(x)|^2 \, dx  \bigg)^{1/2}.$$
It is well known (see, for example, \cite[section 5.2]{AKL2009}) that $\mathcal{H}$ maps $\mathcal{Y}$ onto itself and its null space is the one dimensional space generated by $1/\sqrt{1-x^2}$. Therefore, if we define
\begin{equation}\langlebel{yeps}
\begin{array}{lll}
 \mathcal{X} &=&\displaystyle
\bigg\{  \partialsi \in \mathcal{C}^0 \left(\,\left[\, -1,
1 \,\right]\,\right): \partialsi^{\partialrime} \in \mathcal{
Y},~ \partialsi(-1)=\partialsi(1)=0 \bigg\},
\end{array}
\end{equation}
where  $\partialsi^{\partialrime}$ is the distributional derivative of $\partialsi$, then $\mathcal{A}: \mathcal{X}\rightarrow \mathcal{Y}$ is invertible. We note that $\mathcal{X}$ is a Banach space with the norm
\[\displaystyle
|| \partialsi ||_{\mathcal{X}} \,=\,
||\partialsi||_{L^\infty} + ||
\partialsi^\partialrime||_{\mathcal{Y}} .
\]
Using the Hilbert inversion formula (see, for example,
\cite[section 5.2]{AKL2009}), we can check
\begin{align}
\mathcal{A}^{-1}[ 1 ](x)& =  -\sqrt{1-x^2},\langlebel{invA1} \\
\mathcal{A}^{-1}[y ](x)& = -\frac{x}{2}\sqrt{1-x^2}.\langlebel{invAy}
\end{align}

The equation \eqref{scaleeqn} can be written as
\begin{equation}\langlebel{a_equation}
 \mathcal{A}[ \partialsi_{\varepsilon}](x) = -\frac{4}{E} {\bf f}_\varepsilon (x), \quad -1<x<1.
\end{equation}
The Taylor expansion yields
\begin{equation}
\partiald{{\bf u}_0}{\nu}(\frac{\varepsilon}{2} x,0) =
\frac{\partial {\bf u}_{0}} {\partial \nu}(0) +
\frac{\varepsilon x}{2} \frac{\partial^2  {\bf u}_{0}}  {\partial t \partial \nu}(0)
+ e_1(x),
\end{equation}
where ${\partial}/{\partial t}$ denotes the tangential derivative on $\gamma_{\varepsilon}$. The remainder term $e_1$ satisfies
 $$
 |e_1(x)| \le C \varepsilon^2 |x|^2,
 $$
and in particular,
 \begin{equation}
 \| e_1 \|_{\mathcal{Y}} \le C \varepsilon^2.
 \end{equation}
On the other hand, since
$$
\partiald{}{\nu}\mathcal{D}_\Omega[{\bf w}_\varepsilon](\frac{\varepsilon}{2} x,0) = \int_{\partial\Omega}
\frac{\partial^2 {\bf N}}{\partial \nu_{\bf x}\partial \nu_{\bf y}}((\frac{\varepsilon}{2} x,0),{\bf y})
{\bf w}_\varepsilon({\bf y}) \,d\sigma({\bf y}),
$$
and $\gamma_\varepsilon$ is away from $\partial\Omega$, one can see that \begin{equation} \left\|
\partiald{}{\nu}\mathcal{D}_\Omega[{\bf w}_\varepsilon](\frac{\varepsilon}{2} \cdot,0)
\right\|_{\mathcal{Y}} \le C \| {\bf w}_\varepsilon \|_{L^\infty(\partial\Omega)}.
\end{equation} Therefore, we have \begin{equation} {\bf f}_\varepsilon (x)= \frac{\partial {\bf u}_{0}} {\partial
\nu}(0) + \frac{\varepsilon x}{2} \frac{\partial^2  {\bf u}_{0}}  {\partial t \partial \nu}(0)
+ e(x), \end{equation} where $e$ satisfies \begin{equation} \| e \|_{\mathcal{Y}} \le C
\bigg(\varepsilon^2 + \| {\bf w}_\varepsilon \|_{L^\infty(\partial\Omega)} \bigg). \end{equation}

We now obtain from \eqnref{a_equation} that
\begin{equation}
\partialsi_\varepsilon (x)= -\frac{4}{E} \left[ \frac{\partial {\bf u}_{0}} {\partial \nu}(0) \mathcal{A}^{-1}[1](x) + \frac{\varepsilon}{2} \frac{\partial^2  {\bf u}_{0}}  {\partial t \partial \nu}(0) \mathcal{A}^{-1}[y](x) + \mathcal{A}^{-1}[e](x) \right].
\end{equation}
Note that $\mathcal{E}_1(x) = \mathcal{A}^{-1}[e](x)$ satisfies
$$
\| \mathcal{E}_1 \|_{\mathcal{X}} \le C \| e \|_{\mathcal{Y}} \le C
\bigg(\varepsilon^2 + \| {\bf w}_\varepsilon \|_{L^\infty(\partial\Omega)} \bigg),
$$
and in particular,
\begin{equation}
\| \mathcal{E}_1 \|_{L^\infty(-1,1)} \le C
\bigg(\varepsilon^2 + \| {\bf w}_\varepsilon \|_{L^\infty(\partial\Omega)} \bigg).
\end{equation}
It then follows from \eqnref{invA1} and \eqnref{invAy} that
\begin{equation}
\partialsi_\varepsilon (x)= \frac{4}{E} \left[ \frac{\partial {\bf u}_{0}} {\partial \nu}(0)
\sqrt{1-x^2} + \frac{\varepsilon}{4} \frac{\partial^2  {\bf u}_{0}}  {\partial t \partial
\nu}(0) x \sqrt{1-x^2} + \mathcal{E}_1(x) \right].
\end{equation}
Thus we have from \eqnref{psiGve} that
\begin{equation}\langlebel{GvfGve}
\varphi_\varepsilon (x)= \frac{2}{E} \left[ \frac{\partial {\bf u}_{0}} {\partial \nu}(0) \sqrt{\varepsilon^2 - 4
x^2} + \frac{1}{2} \frac{\partial^2  {\bf u}_{0}}  {\partial t \partial \nu}(0) x
\sqrt{\varepsilon^2 - 4 x^2} + \mathcal{E}(x) \right], \quad (x,0) \in
\gamma_\varepsilon,
\end{equation}
where $\mathcal{E}(x)= \frac{\varepsilon}{2} \mathcal{E}_1(\frac{2}{\varepsilon} x)$ satisfies
\begin{equation}\langlebel{Eest}
\| \mathcal{E} \|_{L^\infty(\gamma_\varepsilon)} \le C \varepsilon \bigg(\varepsilon^2 + \| {\bf w}_\varepsilon
\|_{L^\infty(\partial\Omega)} \bigg).
\end{equation}

Substituting \eqnref{GvfGve} into \eqnref{Neumann_repr} we obtain
\begin{align}
{\bf w}_\varepsilon({\bf x}) &= \frac{2}{E} \int_{\gamma_\varepsilon} \frac{\partial} {\partial \nu_{\bf y}} {\bf N}({\bf x}, (y,0)) \sqrt{\varepsilon^2 - 4 y^2} \,dy \frac{\partial {\bf u}_{0}} {\partial \nu}(0) \nonumber \\
& \quad + \frac{1}{E} \int_{\gamma_\varepsilon} \frac{\partial} {\partial \nu_{\bf y}} {\bf N}({\bf x}, (y,0)) y \sqrt{\varepsilon^2 - 4 y^2} \,dy \frac{\partial^2  {\bf u}_{0}}  {\partial t \partial \nu}(0) \nonumber \\
& \quad + \frac{2}{E} \int_{\gamma_\varepsilon} \frac{\partial} {\partial \nu_{\bf y}} {\bf N}({\bf x}, (y,0)) \mathcal{E}(y) \,dy := I+ II+III , \quad {\bf x} \in \partial \Omega . \langlebel{BuGveBx}
\end{align}
Since
$$
\frac{\partial}{\partial \nu_{\bf y}} {\bf N}({\bf x}, (y,0))= \frac{\partial}{\partial \nu_{\bf y}} {\bf N}({\bf x}, 0) + \frac{\partial^2}{\partial t_{\bf y} \partial \nu_{\bf y}} {\bf N}({\bf x}, 0) y + O(y^2),
$$
we have
\begin{align*}
\int_{\gamma_\varepsilon} \frac{\partial} {\partial \nu_{\bf y}} {\bf N}({\bf x}, (y,0)) \sqrt{\varepsilon^2 - 4 y^2} \,dy
&= \frac{\partial}{\partial \nu_{\bf y}} {\bf N}({\bf x}, 0) \int_{\gamma_\varepsilon} \sqrt{\varepsilon^2 - 4 y^2} \,dy \\
& \quad + \frac{\partial^2}{\partial t_{\bf y} \partial \nu_{\bf y}} {\bf N}({\bf x}, 0) \int_{\gamma_\varepsilon} y \sqrt{\varepsilon^2 - 4 y^2} \,dy + O(\varepsilon^4) \\
&= \frac{\partiali\varepsilon^2}{2} \frac{\partial}{\partial \nu_{\bf y}} {\bf N}({\bf x}, 0)  + O(\varepsilon^4),
\end{align*}
and hence \begin{equation} I= \frac{\partiali\varepsilon^2}{E} \frac{\partial}{\partial \nu_{\bf y}}
{\bf N}({\bf x}, 0) \frac{\partial {\bf u}_{0}} {\partial \nu}(0) + O(\varepsilon^4). \end{equation} Here
and throughout this paper, $O(\varepsilon^4)$ is in the sense of the
uniform norm on $\partial\Omega$. Similarly, one can show that \begin{equation} II=
O(\varepsilon^4). \end{equation} So we obtain that  \begin{equation}\langlebel{BwGve} {\bf w}_\varepsilon({\bf x})
= \frac{\partiali\varepsilon^2}{E} \frac{\partial}{\partial \nu_{\bf y}} {\bf N}({\bf x}, 0) \frac{\partial
{\bf u}_{0}} {\partial \nu}(0) + O(\varepsilon^4) + III. \end{equation} In particular, we
have
$$
\| {\bf w}_\varepsilon \|_{L^\infty(\partial\Omega)} \le C (\varepsilon^2+ |III|).
$$
But, because of \eqnref{Eest}, we arrive at
$$
|III| \le C \varepsilon \| \mathcal{E} \|_{L^\infty(\gamma_\varepsilon)} \le C \varepsilon^2 (\varepsilon^2 + \| {\bf w}_\varepsilon \|_{L^\infty(\partial\Omega)}),
$$
and hence
$$
\| {\bf w}_\varepsilon \|_{L^\infty(\partial\Omega)} \le C \varepsilon^2 (1 + \| {\bf w}_\varepsilon \|_{L^\infty(\partial\Omega)}).
$$
So, if $\varepsilon$ is small enough, then
\begin{equation}
\| {\bf w}_\varepsilon \|_{L^\infty(\partial\Omega)} \le C \varepsilon^2.
\end{equation}
It then follows from \eqnref{BwGve} that
\begin{equation}\langlebel{BwGve2}
{\bf w}_\varepsilon({\bf x}) = \frac{\partiali\varepsilon^2}{E} \frac{\partial}{\partial \nu_{\bf y}} {\bf N}({\bf x}, 0) \frac{\partial {\bf u}_{0}} {\partial \nu}(0) + O(\varepsilon^4) .
\end{equation}

We obtain the following theorem.
\begin{theorem} \langlebel{thm-neuman} 
Suppose that $\gamma_\varepsilon$ is a linear crack of size $\varepsilon$ and ${\bf z}$ is the center of $\gamma_\varepsilon$. Then the solution to \eqnref{prob1} has the following asymptotic expansion:
\begin{equation}\langlebel{exp-u}
  ({\bf u}_{\varepsilon}- {\bf u}_0)({\bf x})
  = \frac{\partiali\varepsilon^2}{E} \frac{\partial {\bf N}} {\partial \nu_{\bf y}} ({\bf x},{\bf y}){\bf i}g|_{{\bf y}={\bf z}}
\frac{\partial {\bf u}_0}{\partial \nu}({\bf z})
  + O(\varepsilon^{4})
     \end{equation}
uniformly on ${\bf x}\in \partial \Omega$. Here $E$ is the Young's modulus.
\end{theorem}

It is worth emphasizing that in \eqnref{exp-u} the error is $O(\varepsilon^4)$ and
the $\varepsilon^3$-term vanishes. One can see from the derivation of \eqnref{exp-u} that
the $\varepsilon^3$-term vanishes because $\gamma_\varepsilon$ is a line segment. If it is a curve,
then we expect that the $\varepsilon^3$-term does not vanish. We also emphasize that \eqnref{exp-u} is a point-wise asymptotic formula, and it can be used to design algorithms to reconstruct cracks from boundary measurements. We can also integrate this formula against the traction ${\bf g}$ to obtain the asymptotic formula for the perturbation of the elastic energy as we do in the next section.

Similarly, if we consider the Dirichlet problem
 \begin{equation}\langlebel{prob21}
\begin{cases}
\nabla \cdot \sigma({\bf u}_{\varepsilon}) = 0 \quad  & \mbox{in}\; \Omega \setminus \bar{\gamma}_{\varepsilon}, \\
{\bf u}_{\varepsilon} = {\bf f}   & \mbox{on} \; \partial \Omega,\\
\sigma({\bf u}_{\varepsilon})\,{\bf e}^\partialerp = 0   & \mbox{on} \;
\gamma_{\varepsilon},
\end{cases}
\end{equation}
and denote the Green function of Lam\'{e} system in $\Omega$ by
${\bf G}$, then we get the following asymptotic expansion of its
solution ${\bf u}_{\varepsilon}$.

\begin{theorem} \langlebel{thm-dirich} 
Suppose that $\gamma_\varepsilon$ is a linear crack of size $\varepsilon$, located at ${\bf z}$. Then the solution to \eqnref{prob21} has the following asymptotic expansion:
\begin{equation}\langlebel{exp-Du}
  \partiald{}{\nu}({\bf u}_{\varepsilon}- {\bf u}_0)({\bf x})
  =  \frac{\partiali\varepsilon^2}{E}
  \frac{\partial^2  {\bf G}}{\partial \nu_{\bf x} \partial \nu_{\bf y}}({\bf x},{\bf y}){\bf i}g|_{{\bf y}={\bf z}}
  \frac{\partial {\bf u}_0}{\partial \nu}({\bf z})
   + O(\varepsilon^{4})
     \end{equation}
uniformly on ${\bf x}\in \partial \Omega$.
\end{theorem}

\section{Topological derivative of the potential energy} \langlebel{sect5}

The elastic potential energy functional of the cracked body is given by \eqnref{energy}, while without the crack the energy functional is given by
\begin{equation}
J[{\bf u}_{0}]
= -\frac{1}{2}\int_{\Omega} \sigma({\bf u}_{0}) :  \nabla^s {\bf u}_{0}.
\end{equation}
By the divergence theorem we have
\begin{equation}
J_{\varepsilon}[{\bf u}_{\varepsilon}] - J[{\bf u}_{0}] = - \frac{1}{2} \int_{\partial\Omega}  ({\bf u}_{\varepsilon}-{\bf u}_0) \cdot {\bf g} \, d\sigma.
\end{equation}
Thus we obtain from \eqref{exp-u}
 \begin{align*}
 J_{\varepsilon}[{\bf u}_{\varepsilon}] - J[{\bf u}_{0}] = - \frac{\partiali\varepsilon^2}{2E}
    \frac{\partial {\bf u}_0}{\partial \nu}({\bf z})\frac{\partial }{\partial \nu_{\bf y}}
\int_{\partial \Omega} {\bf N} ({\bf x},{\bf y})
    {\bf g}({\bf x}) \,d\sigma({\bf x}){\bf i}g|_{{\bf y}={\bf z}}  + O(\varepsilon^{4}).
  \end{align*}
Since
$$
{\bf u}_0({\bf y})= \int_{\partial \Omega} {\bf N} ({\bf x},{\bf y})
    {\bf g}({\bf x}) \,d\sigma({\bf x}), \quad {\bf y} \in \Omega,
$$
we have
\begin{equation}\langlebel{energyform}
 J_{\varepsilon}[{\bf u}_{\varepsilon}] - J[{\bf u}_{0}] = - \frac{\partiali\varepsilon^2}{2E}
{\bf i}g|\frac{\partial {\bf u}_0}{\partial \nu}({\bf z}){\bf i}g|^2 + O(\varepsilon^{4}).
\end{equation}

We may write \eqnref{energyform} in terms of the stress intensity factors. The (normalized) stress intensity factors $K_I$ and $K_{II}$ are defined by
\begin{equation}
K_{I}({\bf u}_0,{\bf e}):= \sigma({\bf u}_0){\bf e}^{\partialerp}\cdot {\bf e}^{\partialerp}\quad \mbox{and}\quad
K_{II}({\bf u}_0,{\bf e}):= \sigma({\bf u}_0){\bf e}^\partialerp \cdot {\bf e}.
\end{equation}
So, we have
\begin{equation}
\sigma ({\bf u}_0){\bf e}^{\partialerp}= K_{I}{\bf e}^{\partialerp}+K_{II}{\bf e},
\end{equation}
and hence
\begin{equation}
{\bf i}g|\frac{\partial {\bf u}_0}{\partial \nu}({\bf z}){\bf i}g|^2 = |\sigma ({\bf u}_0){\bf e}^{\partialerp}|^2 = K_{I}^2 + K_{II}^2.
\end{equation}
We obtain the following result.
\begin{theorem}
We have
\begin{equation}\langlebel{energyform2}
 J_{\varepsilon}[{\bf u}_{\varepsilon}] - J[{\bf u}_{0}] = - \frac{\partiali\varepsilon^2}{2E}
(K_{I}^2 + K_{II}^2) + O(\varepsilon^{4})
\end{equation}
as $\varepsilon \to 0$.
\end{theorem}

The topological derivative $D_T J_{\varepsilon}({\bf z})$ of the potential
energy is defined by \cite{SZ99, SZ01} 
\begin{equation}\langlebel{topd} 
D_T J_{\varepsilon}({\bf z}) :=\lim_{\varepsilon \rightarrow
0}\left(\frac{1}{\rho'(\varepsilon)} \frac{d}{d\varepsilon}J_{\varepsilon}\right), 
\end{equation}
where $\rho(\varepsilon)= \partiali\varepsilon^2$. So, one can immediately see from
\eqnref{energyform2} that 
\begin{equation}\langlebel{plainstress}
D_T J_{\varepsilon}({\bf z}) = - \frac{1}{2E} (K_{I}^2 + K_{II}^2). 
\end{equation} 
This formula is in accordance with the one obtained by Novotny {\it et al} in \cite{NDTP} (see also \cite{GN}).
In fact, in those papers the plane strain and the plain stress problems are considered, and  \eqnref{plainstress} is the formula for the latter problem.

\section{Derivation of \eqnref{repDij}} \langlebel{appendix}

The Kelvin matrix \eqref{lamefund} can be rewritten as
\begin{equation}\langlebel{lamefund2}
   \Phi_{ij}({\bf x}-{\bf y})=\langlembda' \delta_{ij}\log|{\bf x}-{\bf y}|
   + \mu'(x_i-y_i)\frac{\partial \log|{\bf x}-{\bf y}| }{\partial y_j},\quad
   i,j=1,2,
\end{equation}
where
$$
\langlembda' = \frac{\langlembda+3\mu}{4\partiali \mu(\langlembda+2\mu)}, \quad
\mu'= \frac{\langlembda+\mu}{4\partiali\mu(\langlembda+2\mu)}.
$$
Using the operator $T(\partial)$ defined by \eqnref{op-T} one can see that
\begin{equation}\langlebel{derfund-rule}
 \frac{\partial \Phi}{\partial \nu_{\bf y}}({\bf x}-{\bf y})
= (T (\partial_{\bf y})\Phi({\bf x}-{\bf y}))^T,
\end{equation}
or
\begin{equation}\langlebel{defopr}
\displaystyle{\bf i}g(\frac{\partial \Phi}{\partial \nu_{\bf y}}({\bf x}-{\bf y}){\bf i}g)_{kj}:=
\sum_{l=1}^2 T_{jl}(\partial_{\bf y})\Phi_{lk}({\bf x}-{\bf y}), \quad k,j=1,2.
\end{equation}
We use the formulas
\begin{eqnarray*}
 \frac{\partial}{\partial y_i}\log|{\bf x}-{\bf y}| &=&-\frac{x_i-y_i}{|{\bf x}-{\bf y}|^2},\\
\frac{\partial^2}{\partial y_i^2}\log|{\bf x}-{\bf y}| &=&-2\frac{(x_i-y_i)^2}{|{\bf x}-{\bf y}|^4}+ \frac{1}{|{\bf x}-{\bf y}|^2},\\
\frac{\partial^2}{\partial y_i\partial y_j}\log|{\bf x}-{\bf y}| &=&-2\frac{(x_i-y_i)(x_j-y_j)}{|{\bf x}-{\bf y}|^4} \quad\mbox{if } i \neq j.
  \end{eqnarray*}
By \eqref{defopr}, we have
\begin{eqnarray*}
{\bf i}g(\frac{\partial \Phi}{\partial \nu_{\bf y}}{\bf i}g)_{11}
&=& \left((\langlembda +2\mu) n_1 \frac{\partial }{\partial {y_1}} + \mu n_2\frac{\partial }{\partial {y_2}}\right)\left(\langlembda' \log|{\bf x}-{\bf y}| + \mu'(x_1-y_1)\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_1}}\right)\\
&& + \left( \mu n_2 \frac{\partial }{\partial {y_1}} + \langlembda n_1 \frac{\partial }{\partial {y_2}}\right)\mu'(x_1-y_1)\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_2}}\\
&=& \langlembda'(\langlembda +2\mu) n_1 \frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_1}} + \langlembda'\mu n_2\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_2}}\\
&& +\,  \mu'(\langlembda +2\mu) n_1 \left(-\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_1}}+(x_1-y_1)\frac{\partial^2 \log|{\bf x}-{\bf y}|}{\partial {y_1}^2}\right)\\
&& - \mu\mu'n_2\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_2}} + 2\mu\mu'n_2(x_1-y_1)\frac{\partial^2 \log|{\bf x}-{\bf y}|}{\partial {y_1}\partial{y_2}} \\ &&+\langlembda\mu'(x_1-y_1)\frac{\partial^2 \log|{\bf x}-{\bf y}|}{\partial {y_2}^2}.
  \end{eqnarray*}
Since $\Delta \log|{\bf x}-{\bf y}|=0$ for ${\bf x} \neq {\bf y}$, we have
\begin{eqnarray*}
{\bf i}g(\frac{\partial \Phi}{\partial \nu_{\bf y}}{\bf i}g)_{11}
 &=& (\langlembda +2\mu)(\langlembda'-\mu') n_1 \frac{\partial \log|{\bf x}-{\bf y}| }{\partial y_1}
+\mu (\langlembda'-\mu') n_2\frac{\partial\log|{\bf x}-{\bf y}| }{\partial y_2} \\
&& +2\mu\mu' \left( n_1(x_1-y_1)\frac{\partial^2 \log|{\bf x}-{\bf y}|}{\partial {y_1}^2}+n_2(x_1-y_1)\frac{\partial^2 \log|{\bf x}-{\bf y}|}{\partial {y_1}\partial{y_2}}\right).
 \end{eqnarray*}
Since
$$
(\langlembda +2\mu)(\mu'-\langlembda')+2\mu\mu'=-\frac{\mu}{2\partiali (\langlembda +2\mu) }=\mu(\mu'-\langlembda'),
$$
we obtain
\begin{eqnarray*}
{\bf i}g(\frac{\partial \Phi}{\partial \nu_{\bf y}}{\bf i}g)_{11}
=\left[ \mu(\mu'-\langlembda') - 4\mu\mu'\frac{(x_1-y_1)^2}{|{\bf x}-{\bf y}|^2}\right]
 \sum_{l=1}^2 n_l\frac{x_l-y_l}{|{\bf x}-{\bf y}|^2}.
 \end{eqnarray*}

Similarly, we can compute
\begin{eqnarray*}
{\bf i}g(\frac{\partial \Phi}{\partial \nu_{\bf y}}{\bf i}g)_{12}
&=& \left(\langlembda n_2 \frac{\partial }{\partial {y_1}} + \mu n_1\frac{\partial }{\partial {y_2}}\right)\left(\langlembda' \log|{\bf x}-{\bf y}|  + \mu'(x_1-y_1)\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_1}}\right)\\
&& + \left(\mu n_1 \frac{\partial }{\partial {y_1}}
 + (\langlembda +2\mu) n_2 \frac{\partial }{\partial {y_2}}\right)\mu'(x_1-y_1)\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_2}}\\
&=&\langlembda(\langlembda' - \mu')n_2 \frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_1}}
+\mu( \langlembda' -\mu')n_1\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_2}}\\
&&+2\mu\mu'n_1(x_1-y_1) \frac{\partial ^2\log|{\bf x}-{\bf y}|}{\partial {y_1} \partial{y_2}}
+ 2\mu\mu' n_2 (x_1-y_1)\frac{\partial ^2\log|{\bf x}-{\bf y}|}{\partial {y_2}^2}.
  \end{eqnarray*}
Since $\langlembda(\mu'-\langlembda')+2\mu\mu'=\mu(\langlembda'-\mu')$, we have
\begin{eqnarray*}
{\bf i}g(\frac{\partial \Phi}{\partial \nu_{\bf y}}{\bf i}g)_{12}
&=&[\langlembda(\mu'-\langlembda')+2\mu\mu']n_2 \frac{x_1-y_1}{|{\bf x}-{\bf y}|^2}
+\mu(\mu'- \langlembda')n_1\frac{x_2-y_2}{|{\bf x}-{\bf y}|^2} \\
&&-4\mu\mu'n_1 \frac{(x_1-y_1)^2(x_2-y_2)}{|{\bf x}-{\bf y}|^4}
-4\mu\mu' n_2 \frac{(x_1-y_1)(x_2-y_2)^2}{|{\bf x}-{\bf y}|^4}\\
 &=& -4\mu\mu'\frac{(x_1-y_1)(x_2-y_2)}{|{\bf x}-{\bf y}|^2}
 \sum_{l=1}^2 n_l\frac{x_l-y_l}{|{\bf x}-{\bf y}|^2} \\
&&  \quad -\,\mu(\mu'-\langlembda') \left[\frac{n_2(x_1-y_1) - n_1(x_2-y_2)}{|{\bf x}-{\bf y}|^2}\right].
  \end{eqnarray*}

We also have
\begin{align*}
{\bf i}g(\frac{\partial \Phi}{\partial \nu_{\bf y}}{\bf i}g)_{22}
& = \left(\langlembda n_2 \frac{\partial }{\partial {y_1}} + \mu n_1\frac{\partial }{\partial {y_2}}\right)\mu'(x_2-y_2)\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_1}}\\
&+\left(\mu n_1\frac{\partial }{\partial {y_1}}+(\langlembda +2\mu)n_2\frac{\partial }{\partial {y_2}}\right) \left(\langlembda' \log|{\bf x}-{\bf y}|  + \mu'(x_2-y_2)\frac{\partial \log|{\bf x}-{\bf y}|}{\partial {y_2}}\right)\\
& = \left[ \mu(\mu'-\langlembda') - 4\mu\mu'\frac{(x_2-y_2)^2}{|{\bf x}-{\bf y}|^2}\right]
 \sum_{l=1}^2 n_l\frac{x_l-y_l}{|{\bf x}-{\bf y}|^2}.
 \end{align*}
 This proves \eqref{repDij}.

\end{document}