%~Mouliné par MaN_auto v.0.32.0 2024-04-22 12:35:26
\documentclass[CRMATH, Unicode, XML, biblatex, published]{cedram}

\TopicFR{\'Equations aux dérivées partielles, Mécanique}
\TopicEN{Partial differential equations, Mechanics}

\addbibresource{CRMATH_LeLouer_20210627.bib}

\usepackage{mathrsfs}
\usepackage{upgreek}
\usepackage{mathtools}

%\usepackage{ulem}

\newcommand*{\dd}{\mathrm{d}}
\newcommand*{\dxi}{\dd \xi}
\newcommand{\ds}{\dd s}
\newcommand{\dx}{\dd x}

\DeclareMathOperator{\Div}{div}
\DeclareMathOperator{\Rot}{curl}
\DeclareMathOperator{\rot}{curl}
\DeclareMathOperator{\Id}{Id}
\DeclareMathOperator{\Imop}{Im}
\DeclareMathOperator{\Reop}{Re}
\DeclareMathOperator{\Cap}{cap}
\DeclareMathOperator{\Trace}{Trace}
\DeclareMathOperator{\Ker}{Ker}
%\DeclareMathOperator{}{}

\newcommand{\vnabla}{\operatorname{\boldsymbol{\nabla}}}
\newcommand{\loc}{\mathrm{loc}}
\newcommand{\inc}{\mathrm{inc}}
\renewcommand{\Im}{\Imop}
\renewcommand{\Re}{\Reop}

\newcommand{\R}{\mathbb{R}}
\newcommand{\N}{{\mathbb N}}
\newcommand{\LL}{\boldsymbol{L}}
\newcommand{\HH}{\boldsymbol{H}}
\newcommand{\TT}{\boldsymbol{T}}
\newcommand{\VV}{\boldsymbol{V}}
\newcommand{\WW}{\boldsymbol{W}}

\newcommand{\nn}{\boldsymbol{n}}
\newcommand{\xb}{\boldsymbol{x}}
\newcommand{\yb}{\boldsymbol{y}}
\newcommand{\zb}{\boldsymbol{z}}
\newcommand{\uu}{\boldsymbol{u}}
\newcommand{\ttt}{\boldsymbol{t}}
\newcommand{\vv}{\boldsymbol{v}}
\newcommand{\ki}{\boldsymbol{\xi}}
\newcommand{\zero}{\boldsymbol{0}}
\newcommand{\vpsi}{\boldsymbol{\psi}}
\newcommand{\vphi}{\boldsymbol{\varphi}}

\newcommand{\0}{\boldsymbol{0}}

\renewcommand*{\to}{\mathchoice{\longrightarrow}{\rightarrow}{\rightarrow}{\rightarrow}}

\newcommand{\transposee}[1]{{\vphantom{#1}}^{\sf T}{#1}}
%Better widebar

\makeatletter
\let\save@mathaccent\mathaccent
\newcommand*\if@single[3]{%
  \setbox0\hbox{${\mathaccent"0362{#1}}^H$}%
  \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}%
  \ifdim\ht0=\ht2 #3\else #2\fi
  }
%The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath:
\newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna}
%If there's a superscript following the bar, then no negative kern may follow the bar;
%an additional {} makes sure that the superscript is high enough in this case:
\newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}}
%Use a separate algorithm for single symbols:
\newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}}
\newcommand*\wide@bar@[3]{%
  \begingroup
  \def\mathaccent##1##2{%
%Enable nesting of accents:
    \let\mathaccent\save@mathaccent
%If there's more than a single symbol, use the first character instead (see below):
    \if#32 \let\macc@nucleus\first@char \fi
%Determine the italic correction:
    \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}%
    \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}%
    \dimen@\wd\tw@
    \advance\dimen@-\wd\z@
%Now \dimen@ is the italic correction of the symbol.
    \divide\dimen@ 3
    \@tempdima\wd\tw@
    \advance\@tempdima-\scriptspace
%Now \@tempdima is the width of the symbol.
    \divide\@tempdima 10
    \advance\dimen@-\@tempdima
%Now \dimen@ = (italic correction / 3) - (Breite / 10)
    \ifdim\dimen@>\z@ \dimen@0pt\fi
%The bar will be shortened in the case \dimen@<0 !
    \rel@kern{0.6}\kern-\dimen@
    \if#31
      \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}%
      \advance\dimen@0.4\dimexpr\macc@kerna
%Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows:
      \let\final@kern#2%
      \ifdim\dimen@<\z@ \let\final@kern1\fi
      \if\final@kern1 \kern-\dimen@\fi
    \else
      \overline{\rel@kern{-0.6}\kern\dimen@#1}%
    \fi
  }%
  \macc@depth\@ne
  \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar
  \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}%
  \macc@set@skewchar\relax
  \let\mathaccentV\macc@nested@a
%The following initialises \macc@kerna and calls \mathaccent:
  \if#31
    \macc@nested@a\relax111{#1}%
  \else
%If the argument consists of more than one symbol, and if the first token is
%a letter, use that letter for the computations:
    \def\gobble@till@marker##1\endmarker{}%
    \futurelet\first@char\gobble@till@marker#1\endmarker
    \ifcat\noexpand\first@char A\else
      \def\first@char{}%
    \fi
    \macc@nested@a\relax111{\first@char}%
  \fi
  \endgroup
}

\def\@setafterauthor{%
  \vglue3mm%
%  \hspace*{0pt}%
\begingroup\hsize=12cm\advance\hsize\abstractmarginL\raggedright
\noindent
%\hspace*{\abstractmarginL}\begin{minipage}[t]{10cm}
   \leftskip\abstractmarginL
  \normalfont\Small
  \@afterauthor\par
\endgroup
\vskip2pt plus 3pt minus 1pt
}
\makeatother

\let\oldbar\bar
\renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}}

\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\renewcommand{\labelenumi}{(\theenumi)}}

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}


\title{Boundary integral equation methods for Lipschitz domains in linear elasticity}

\author{\firstname{Fr{\'e}d{\'e}rique} \lastname{Le Lou{\"e}r}}
\address{Alliance Sorbonne Universit\'e, Universit\'e de Technologie de Compi\`egne, LMAC EA2222 Laboratoire de Math\'ematiques Appliqu\'ees de Compi\`egne - CS 60319 - 60203 Compi\`egne cedex, France}
\email{frederique.le-louer@utc.fr}

\keywords{Boundary integral equation, Linear elasticity, Lipschitz domains, G\r{a}rding's inequality, Eigenvalues}
\subjclass{35J25, 35P25, 35C15, 74B05, 65N12}

\begin{abstract}
A review of stable boundary integral equation methods for solving the Navier equation with either Dirichlet or Neumann boundary conditions in the exterior of a Lipschitz domain is presented. The conventional combined-field integral equation (CFIE) formulations, that are used to avoid spurious resonances, do not give rise to a coercive variational formulation for nonsmooth geometries anymore. To circumvent this issue, either the single layer or the double layer potential operator is composed with a compact or a Steklov--Poincar\'e type operator. The later can be constructed from the well-know elliptic boundary integral operators associated to the Laplace equation and G\r{a}rding's inequalities are satisfied. Some Neumann interior eigenvalue computations for the unit square and cube are presented for forthcoming numerical investigations.
\end{abstract}

\begin{document}
\maketitle

\section{Introduction and problem settings}

Let $\Omega$ denote a bounded domain in $\R^{d}$, $d=2$ or $3$ and let $\Omega^c$ denote the exterior domain $\R^d\backslash\bar{\Omega}$. We assume that the boundary $\Gamma$ of $\Omega$ is a Lipschitz continuous (in the sense of~\cite[Definition~2.4.5]{HenrotPierre} or~\cite[p.~89]{McLean}), simply connected and oriented closed curve or surface. Let $\nn$ denote the outward unit normal vector on the boundary $\Gamma$. We assume that the Lam\'e parameters $\mu$ and $\lambda$ and the density $\rho$ are positive constants. The propagation of time-harmonic elastic waves in the three-dimensional isotropic and homogeneous elastic medium (with $e^{-i\omega t}$ time-dependence) is described by the Navier equation
\begin{equation}\label{NE}
\Updelta^*\uu+ \rho\omega^2\uu=\mu\Delta\uu+(\lambda+\mu)\vnabla\Div\uu+ \rho\omega^2\uu=\zero,
\end{equation}
where $\omega>0$ is the frequency. In this work we discuss solution methods based on boundary integral equation formulations for the following exterior boundary value problems in elastodynamics: Given an incident field $\uu^{\inc}$ which is assumed to solve the Navier equation in the absence of scatterer, find the solution $\uu$ to the Navier equation~\eqref{NE} in $\Omega^c$ which satisfies either a Dirichlet boundary condition
\begin{equation}\label{dirichlet}
\uu+\uu^{\inc}=\zero\quad\text{on }\Gamma
\end{equation}
or a Neumann boundary condition
\begin{equation}\label{neumann}
\TT(\uu+\uu^{\inc})=\zero\quad\text{on }\Gamma,
\end{equation}
where the traction operator is defined by
\begin{equation}\label{traction3}
\TT\uu=2\mu\,\nn\cdot(\nabla \uu)_{|_{\Gamma}} +\lambda \,\nn (\Div\uu)_{|_{\Gamma}} +\mu\,\nn\times(\Rot\uu)_{|_{\Gamma}} \text{ when $d=3$}.
\end{equation}
When $d=2$, by setting $\uu=(u_{1},u_{2},0)$ and $\nn=(n_{1},n_{2},0)$, we get
\begin{equation}\label{traction3b}
\TT\uu=2\mu\,\nn\cdot(\nabla \uu)_{|_{\Gamma}} +\lambda \,\nn (\Div\uu)_{|_{\Gamma}} -\mu\,\nn^{\bot}(\rot_{||}\uu)_{|_{\Gamma}}\text{ when $d=2$},
\end{equation}
where $\nn^{\bot}=(-n_{2},n_{1})$ and $\rot_{||}\uu=\partial_{x_{1}}u_{2}-\partial_{x_{2}}u_{1}$. Incident plane waves satisfy the homogeneous Navier equation in $\R^d$. Other incident fields are assumed to be generated by the incident source $\ttt^{\inc}$ with compact support in $\Omega^c$, so that $\Updelta^*\uu^{\inc}+ \rho\omega^2\uu^{\inc}=\ttt^{\inc}$. In addition the scattered field $\uu$ has to satisfy the Dirichlet-to-Neumann type radiation condition~\cite{GachterGrote}
\begin{equation}\label{radiation}
\lim_{|\xb|\to\infty}|\xb|^{\frac{d-1}{2}}\left|\TT\uu-i\mathcal{K}_{\omega}\uu\right|=0,
\end{equation}
which has to be satisfied uniformly for all unitary directions $\hat{\xb}=\frac{\xb}{|\xb|}$ and where
\[
\mathcal{K}_{\omega}=\kappa_{p}(\lambda+2\mu)\Id_{\hat{\xb}} +\kappa_{s}\mu \big(\Id-\Id_{\hat{\xb}}\big),
\]
$\Id_{\hat{\xb}}=\hat{\xb}\transposee{\hat{\xb}}$, $\kappa_{p}=\omega\sqrt{\frac{\rho}{\lambda+2\mu}}$ and $\kappa_{s}=\omega\sqrt{\frac{\rho}{\mu}}$ are the $P$- and $S$-wave numbers associated to longitudinal and transverse wave propagation, respectively. This radiation condition~\eqref{radiation} is equivalent to the Kupradze radiation conditions and avoids the splitting of scattered waves as the sum of share ($\Div\uu=0$) and pressure $(\rot_{||}\uu=0$ for $d=2$ or $\Rot\uu=0$ for $d=3$) waves in the analysis. Uniqueness of the solution to problems~\eqref{NE}--\eqref{dirichlet}--\eqref{radiation} and~\eqref{NE}--\eqref{neumann}--\eqref{radiation} is guaranted by~\cite[Lemma~2.1]{IvanyshynLeLouer} for the two dimensional case. The analogous uniqueness results for the three-dimensional case will be detailed in~\cite{LeLouerRais}.

A classical tool for establishing existence results for the solution to such elliptic boundary value problems is the reduction to coercive combined-field boundary integral equations. However, if the boundary $\Gamma$ is not smooth, the conventional linear combination of the single and double layer potential operators (required to avoid spurious frequencies) do not satisfy G\r{a}rding's inequalities anymore. This issue has been the subject of several works in acoustic and electromagnetic scattering~\cite{BuffaHiptmair,EnglederSteinbach1, EnglederSteinbach2,SteinbachWindisch}. From these three-dimensional studies, two different techniques emerge. Either we \emph{regularize} one of the potential operators by a compact operator or else we \emph{modify} them by the introduction of a Steklov--Poincar\'e type operator. Both of the ideas extend to the elastic case in two and three dimensions.

The sequel is organized as follows. In Section~\ref{Sec2}, collecting theories presented in~\cite{Kupradze,McLean}, we recall some standard results about trace mappings and regularity properties of the boundary integral operators in suitable Sobolev spaces for Lipschitz boundaries. With the help of the radiation condition~\eqref{radiation}, we then extend from acoustics to elastodynamics some well-known positivity results satisfied by the weakly and hypersingular boundary integral operators. Instead of constructing the Steklov--Poincar\'e type operators from the purely elastic potential theory we suggest to reuse the elliptic operators associated to the Laplace equation case. The regularized or modified CFIEs are presented and analysed in Section~\ref{MCFIE}. We only describe the so-called indirect approach (with non physical unknown density) since the direct one leads to the $L^2$-adjoint form of the proposed modified CFIEs. Finally, we outline concluding remarks and discuss possible research lines in the final section of the paper. The content is completed by two appendices providing the principal symbols of the boundary integral operators and explicit Neumann eigenvalue computations.


\section{Potential theory for linear elastic waves scattering}\label{Sec2}

\subsection{Trace mapping properties for Lipschitz domains}

We summarize some well known results about traces of vector fields and variational formulations for time-harmonic elastodynamic waves in a Lipschitz domain.

\begin{nota*}
For $G\subset\R^d$, we denote by $H^s(G)$ the usual $L^2$-based Sobolev space of order $s\in\R$, and by $H^s_{\loc}(\bar G)$ the space of functions whose restrictions to any bounded subdomain $B$ of $G$ belong to $H^s(B)$, with the convention $H^0\equiv L^2$. Spaces of vector functions will be denoted by boldface letters; thus $\HH^s=(H^s)^d$. We set :
\begin{align*}
\HH^1(G,\Updelta^*)&:=\left\{\uu\in \HH^1(G):\;\Updelta^*\uu\in \LL^2(G)\right\},\\
\HH^1_{\loc}(G,\Updelta^*)&:=\left\{\uu\in \HH^1_{\loc}(\bar{G}):\;\Updelta^*\uu\in \LL_{\loc}^2(\bar{G})\right\}.
\end{align*}
The space $\HH^1(\Omega, \Updelta^{*})$ is a Hilbert space endowed with the natural graph norm.
\end{nota*}

\begin{defi}\label{DefTrace}
For a vector function $\uu\in (\mathscr{C}^{\infty}(\bar{\Omega}))^d$ we define the
traces :
\begin{align*}
\gamma_{0}^- \uu&=\uu_{|_{\Gamma}}\text{ (Dirichlet)}\\
\gamma_{1}^-\uu&=\TT\uu\text{ (Neumann)}\\
\gamma_{g}^-\uu&=\nn\cdot(\nabla \uu)_{|_{\Gamma}} - \nn (\Div\uu)_{|_{\Gamma}} +\nn\times(\Rot\uu)_{|_{\Gamma}}, \text{ when $d=3$,}\\
&= \nn\cdot(\nabla \uu)_{|_{\Gamma}} - \nn (\Div\uu)_{|_{\Gamma}} -\nn^{\bot}(\rot_{||}\uu)_{|_{\Gamma}},\text{ when $d=2$},\text{ (G\"unter derivative).}
\end{align*}
\end{defi}

We use standard Sobolev spaces $H^t(\Gamma)$, $t\in[-1,1]$, endowed with standard norms $||\cdot||_{H^t(\Gamma)}$ and with the convention $H^0(\Gamma)=~L^2(\Gamma)$. Spaces of vector densities are denoted by boldface letters, thus $\HH^t(\Gamma)=\big(H^t(\Gamma)\big)^d$. By density arguments we extend the Definition~\ref{DefTrace} to Sobolev spaces. The trace maps
\begin{align*}
\gamma_{0}^-:\HH^{s+\frac{1}{2}}(\Omega)&\to \HH^{s}(\Gamma),
\end{align*}
are continuous for all $s>0$, if the domain is smooth (see~\cite[Theorem~3.37]{McLean}). Whereas for a bounded Lipschitz domain in general, the validity is only given for $s\in(0,1)$ (see~\cite[Lemma~3.6]{Costabel} or~\cite[Theorem~3.38]{McLean}). The dual spaces of $\HH^t(\Gamma)$ with respect to the $\LL^2$ scalar product are denoted by $\HH^{-t}(\Gamma)$. The brackets $\langle\,\cdot\,,\cdot\,\rangle_{\frac{1}{2}}$ represents the duality product defined for $(\vphi,\vpsi)\in\HH^{-\frac{1}{2}}(\Gamma)\times\HH^{\frac{1}{2}}(\Gamma)$~as
\[
\langle\vphi,\vpsi\rangle_{\frac{1}{2}}=\int_{\Gamma}\vphi\cdot\overline{\vpsi} \,\ds=\overline{\langle\vpsi,\vphi\rangle_{\frac{1}{2}}}.
\]

By the Gauss divergence theorem~\cite[p.~115]{McLean}, we derive the mapping properties of $\gamma^-_{1}$ and $\gamma^-_{g}$. Let set
\[
\begin{aligned} \upsigma(\uu)&=\lambda(\Div\uu)\Id+2\mu\upxi(\uu),\quad\text{with }\upxi(\uu)=\tfrac{1}{2}([\nabla \uu]+\transposee{[\nabla\uu]}),\\
\upchi(\uu)&=\transposee{[\nabla \uu]}-(\Div\uu)\Id,
\end{aligned}
\]
where $\Id$ is the $d$-by-$d$ identity matrix, $\transposee{M}$ is the tranpose of matrix $M$ and $[\nabla\uu]$ is the matrix whose the $j$th column is the gradient of the $j$-th component of $\uu$.

For two $d$-by-$d$ matrices $A$ and $B$ whose columns are denoted by $(a_{1},\hdots,a_{d})$ and $(b_{1},\hdots,b_{d})$, respectively, we set $A:B= a_{1}\cdot b_{1}+\cdots+a_{d}\cdot b_{d}$. From the identities
\begin{equation}\label{idtt1}
\begin{aligned}
\Div\big(\upsigma(\uu)\vv\big)&=\Updelta^*\uu\cdot\vv+\upsigma(\uu):\upvarepsilon(\vv), & \transposee{[\nn\cdot\sigma(\uu)_{|\Gamma}]}&=\TT\uu,\\
\Div\big(\upchi(\uu)\vv\big)&=\upchi(\uu):[\nabla\vv], & \transposee{[\nn\cdot\upchi(\uu)_{|\Gamma}]}&=\gamma^-_{g}(\uu),
\end{aligned}
\end{equation}
we deduce the following two Lemmas.

\begin{lemm}\label{Green}
For vector functions $\uu$ and $\vv$ in $\HH^1(\Omega,\Updelta^*)$, we have the following first Green formula,
\begin{equation}\label{IPP1}
\int_{\Omega}\Updelta^*\uu\cdot\vv \,\dx+\int_{\Omega}\upsigma(\uu):\upxi(\vv)\,\dx=\int_{\Gamma}\TT\uu\cdot \vv\, \ds.
\end{equation}
The symmetry of the product $\upsigma(\uu):\upvarepsilon(\vv)=\lambda(\Div\uu)(\Div\vv)+2\mu\,\upxi(\uu):\upxi(\vv)=\upsigma(\vv):\upxi(\uu)$ yields the second Green formula,
\begin{equation}\label{IPP2}
\int_{\Omega}\left(\uu\cdot\Updelta^*\vv-\Updelta^*\uu\cdot\vv\right)\dx=\int_{\Gamma}\big(\uu\cdot \TT\vv-\TT\uu\cdot\vv\big)\ds.
\end{equation}
If $\uu$ and $\vv$ solve the Navier equation in $\Omega$, then each term in~\eqref{IPP2} vanishes.
\end{lemm}

\begin{lemm}\label{DefGunter}
For vector functions $\uu$ and $\vv$ in $\HH^1(\Omega,\Updelta^*)$, we have the following formula,
\begin{equation}\label{IPP3}
\int_{\Omega}\upchi(\uu):[\nabla\vv]\,\dx=\int_{\Gamma}\gamma^-_{g}(\uu)\cdot \vv\, \ds.
\end{equation}
The symmetry of the product $\upchi(\uu):[\nabla\vv]=\upchi(\vv):[\nabla\uu]$ yields the symmetry property of the G\"unter derivative,
\begin{equation}\label{IPP4}
\int_{\Gamma}\big(\uu\cdot \gamma^-_{g}(\vv)-\gamma^-_{g}(\uu)\cdot\vv\big)\ds=0.
\end{equation}
\end{lemm}

The first Green formul{\ae}~\eqref{IPP1}-\eqref{IPP3} yield the continuous trace mappings
\begin{align*}
\gamma^-_{1}:\HH^{1}(\Omega,\Updelta^*)&\to \HH^{-\frac{1}{2}}(\Gamma), \\
\gamma^-_{g}:\HH^{1}(\Omega,\Updelta^*)&\to \HH^{-\frac{1}{2}}(\Gamma).
\end{align*}
For $\uu\in \HH^1_{\loc}({\Omega^c},\Delta^*)$ we define $\gamma_{0}^+\uu$, $\gamma_{1}^+\uu $ and $\gamma_{g}^+\vv $ in the same way and the same mapping properties hold true.

Now we introduce some surface differential operators: The tangential gradient denoted by $\nabla_{\Gamma}$ and the surface divergence denoted by $\Div_{\Gamma}$ for non tangential vector fields. For their definitions we refer to~\cite[paragraph 5.4.3]{HenrotPierre} and we have $\Div_{\Gamma}\vpsi=(\Trace[\nabla_{\Gamma}\vpsi])$. Their extensions to the Sobolev space of fractional order $t=\frac{1}{2}$ are both derived from~\cite[Proposition~3.6]{BuffaCostabelSheen}. From the decomposition of the gradient and divergence operators given in~\cite[Definitions 4.5.5 and 4.5.6]{HenrotPierre}, we deduce the following identity first derived in~\cite[Chapter~$V$, Eq.~(1.6)]{Kupradze}
\[
\gamma_{g}^{\pm}(\uu)=[\nabla_{\Gamma}\uu_{|\Gamma}]\nn-(\Div_{\Gamma}\uu_{|\Gamma})\nn:=\mathcal{M}(\uu_{|\Gamma}),
\]
so that the tangential derivative $\mathcal{M}$ is symmetric and continuous from $\HH^{\frac{1}{2}}(\Gamma)$ to $\HH^{-\frac{1}{2}}(\Gamma)$. In particular, we will use the following rewritting of the traction trace
\begin{equation}\label{Tu2}
\begin{aligned}
\TT\uu&=2\mu\mathcal{M}\uu_{|\Gamma} + (\lambda+2\mu)\nn (\Div\uu)_{|_{\Gamma}} +\mu(\Rot\uu)_{|_{\Gamma}}\times\nn, &&\text{when $d=3$,}\\
&=2\mu\mathcal{M}\uu_{|\Gamma} + (\lambda+2\mu)\nn (\Div\uu)_{|_{\Gamma}} +\mu\nn^{\bot}(\rot_{||}\uu)_{|_{\Gamma}},&&\text{when $d=2$}.
\end{aligned}
\end{equation}

\subsection{Boundary integral reformulation of the elastodynamic system}\label{subsec2}

In this paragraph, we introduce the classical potential theory associated to elastodynamic wave scattering problems and we establish some positivity criteria for the two boundary integral operators with either a weakly or a hypersingular kernel.

Denoting by $G(\kappa,\zb)$ the fundamental solution to the Helmholtz equation
\[
(\Delta+\kappa^2\Id)G(\kappa,\zb)=-\delta_{0}(|\zb|),
\]
the fundamental solution to the Navier equation is then given by
\[
\Phi_{\omega}(\zb)=\frac{1}{\mu}\left(G(\kappa_{s},\zb)\Id_{\R^d}+\frac{1}{\kappa_{s}^2}\vnabla_{\zb}\transposee{\vnabla}_{\zb}\Big(G(\kappa_{s},\zb)-G(\kappa_{p},\zb)\Big)\right).
\]
It is a $d\times d$ symmetric matrix-valued function and behaves as $\Phi_{\omega}(\zb)=O\left(\frac{1}{|\zb|^{\frac{d-1}{2}}}\right)$ as $|\zb|\to+\infty$.

For a solution to the Navier equation~\eqref{NE} in $\Omega^c$ that satisfies the Kupraze radiation condition, one can derive from~\eqref{IPP2} the Somigliana integral representation formula for $\xb\in\Omega^c$:
\begin{equation}\label{intrepext}
\uu(\xb)=\int_{\Gamma}\left(\transposee{\left[\TT_{\!\yb}\Phi_{\omega}(\xb-\yb)\right]}\uu(\yb)-\Phi_{\omega}(\xb-\yb)\TT_{\!\yb}\uu(\yb)\right)\ds(\yb),
\end{equation}
where the lower subscript $_{\yb}$ indicates differentiation with respect to $\yb$ and $\TT_{\!\yb}\Phi_{\omega}(\xb-\yb)$ is the matrix obtained by applying the traction operator $\TT_{\!\yb}$ to each column of $\Phi_{\omega}(\xb-\yb)$. Interior solutions to the Navier equation satisfy for $\xb\in\Omega$
\begin{equation}\label{intrepextn}
\uu(\xb)=-\int_{\Gamma}\left(\transposee{\left[\TT_{\!\yb}\Phi_{\omega}(\xb-\yb)\right]}\uu(\yb)-\Phi_{\omega}(\xb-\yb)\TT_{\!\yb}\uu(\yb)\right)\ds(\yb).
\end{equation}
Conversely, boundary integral equation methods consider ansatz based on the following result.

\begin{lemm}\label{potentials}\ 
\begin{enumerate}\romanenumi
\item \label{4i} The single layer potential operator $\Psi^{\!S}_{\omega}$ defined by
\[
\Psi^{\!S}_{\omega}\vphi(\xb)=\int_{\Gamma}\Phi_{\omega}(\xb-\yb)\vphi(\yb)\ds(\yb)
\]
is
continuous from $\HH^{-\frac{1}{2}}(\Gamma)$ to $\HH^1(\Omega,\Updelta^*)\cup\HH^1_{\loc}(\Omega^c,\Updelta^*)$. For $\vphi\in \HH^{-\frac{1}{2}}(\Gamma)$ we have
\[
(\Updelta^* +\rho\omega^2\Id)\Psi^{\!S}_{\omega}\vphi = 0\quad \text{in\/ }\R^d\backslash\Gamma,
\]
and $\Psi^{\!S}_{\omega}\vphi$
satisfies the radiation condition~\eqref{radiation} and $\big(\Psi^{\!S}_{\omega}\vphi\big)(\xb)=O\left(\frac{1}{|\xb|^{\frac{d-1}{2}}}\right)$ as $|\xb|\to+\infty$.

\item \label{4ii} Similar results hold true for the double layer potential operator $\Psi^{\!D}_{\omega}$ defined for $\vpsi\in\HH^{\frac{1}{2}}(\Gamma)$~by
\[
\Psi^{\!D}_{\omega}\vpsi(\xb)=\int_{\Gamma}\transposee{\left[\TT_{\!\yb}\Phi_{\omega}(\xb-\yb)\right]}\vpsi(\yb)\ds(\yb).
\]
\end{enumerate}
\end{lemm}

The Calder\'on projectors for the time-harmonic Navier equation are
\[
P^{\pm}_{\omega}=
\begin{pmatrix}\pm\tfrac{1}{2}\Id+D_{\omega}&-S_{\omega}\\N_{\omega}&\pm\tfrac{1}{2}\Id-D'_{\omega}
\end{pmatrix},
\]
where the boundary integral operators are defined by
\begin{align}\label{Sdef}
S_{\omega}\vphi(\xb)&=\frac{1}{2}\{\gamma_{0}^++\gamma_{0}^-\}\Psi^{\!S}_{\omega}\vphi=\int_{\Gamma}\Phi_{\omega}(\xb,\yb)\vphi(\yb)\,\ds(\yb),\\\label{Ddef}
D_{\omega}\vpsi(\xb)&=\frac{1}{2}\{\gamma_{0}^++\gamma_{0}^-\}\Psi^{\!D}_{\omega}\vpsi=\int_{\Gamma}\transposee{[\TT_{\!\yb}\Phi_{\omega}(\xb,\yb)]}\vpsi(\yb)\,\ds(\yb),\\\label{Dpdef}
D'_{\omega}\vphi(\xb)&=\frac{1}{2}\{\gamma_{1}^++\gamma_{1}^-\}\Psi^{\!S}_{\omega}\vphi=\int_{\Gamma}\TT_{\!\xb}\left\{\Phi_{\omega}(\xb,\yb)\vphi(\yb)\right\}\ds(\yb),\\\label{Ndef}
N_{\omega}\vpsi(\xb)&=\frac{1}{2}\{\gamma_{1}^++\gamma_{1}^-\}\Psi^{\!D}_{\omega}\vpsi=\int_{\Gamma}\TT_{\!\xb}\left\{\transposee{[\TT_{\!\yb}\Phi_{\omega}(\xb,\yb)]}\vpsi(\yb)\right\}\ds(\yb).
\end{align}
The operator $S_{\omega}$ has a pseudo-homogeneous kernel of order $-1$ (see Appendix~\ref{Ap-A}) and is thus bounded from $\HH^{-\frac{1}{2}}(\Gamma)$ to $\HH^{\frac{1}{2}}(\Gamma)$. The operators $D_{\omega}$ and $D'_{\omega}$ have a pseudo-homogeneous kernel of order 0 and are bounded on $\HH^{\frac{1}{2}}(\Gamma)$ and $\HH^{-\frac{1}{2}}(\Gamma)$ respectively. The operator $N_{\omega}$ has a pseudo-homogeneous kernel of order $+1$ and is bounded from $\HH^{\frac{1}{2}}(\Gamma)$ to $\HH^{-\frac{1}{2}}(\Gamma)$. More precisely we have
\begin{equation}\label{calderon}
P^{\pm}
\begin{pmatrix}\vpsi\\\vphi
\end{pmatrix}=
\begin{pmatrix}\gamma_{0}^{\pm}\Psi_{\omega}^{\!D}\vpsi\\\gamma_{1}^{\pm}\Psi_{\omega}^{\!D}\vpsi
\end{pmatrix}-
\begin{pmatrix}\gamma_{0}^{\pm}\Psi_{\omega}^{\!S}\vphi\\\gamma_{1}^{\pm}\Psi_{\omega}^{\!S}\vphi
\end{pmatrix},
\end{equation}
so that the potential operators $\Psi_{\omega}^{\!S}$ and $\Psi_{\omega}^{\!D}$ satisfy the following jump relations across $\Gamma$
\begin{equation}\label{jump}
\begin{aligned}
\big(\gamma^+_{0}-\gamma_{0}^-\big)\Psi_{\omega}^{\!S}&=\zero, &\big(\gamma^+_{0}-\gamma_{0}^-\big)\Psi_{\omega}^{\!D}&=\Id,\\
\big(\gamma^+_{1}-\gamma_{1}^-\big)\Psi_{\omega}^{\!S}&=-\Id, &\big(\gamma^+_{1}-\gamma_{1}^-\big)\Psi_{\omega}^{\!D}&=\zero.
\end{aligned}
\end{equation}

For numerical implementation purpose, the integration by part formula~\eqref{IPP4} (initially established in~\cite[Chapter~$V$ Theorem~1.3]{Kupradze} for smooth boundaries and still valid for Lipschitz boundaries) permits to express the above four boundary integral operators in terms of other operators with weakly singular kernels and surface differential derivatives. Adapting the computations in~\cite[p.~47-49]{HsiaoWendland} to the dynamic case $\omega\not=0$, the double layer boundary integral operator and the traction derivative of the single layer potential operator can be rewritten for $\xb\in\Gamma$ as:
\begin{multline}\label{D-ipp}
(D_{\omega}\vpsi)(\xb)=2\mu(S_{\omega}\mathcal{M}\vpsi)(\xb)-{\int_{\Gamma}G(\kappa_{s},\xb-\yb)\mathcal{M}_{\yb}\vpsi(\yb)\,\ds(\yb)}+{\int_{\Gamma}\dfrac{\partial}{\partial\nn(\yb)}G(\kappa_{s},\xb-\yb)\vpsi(\yb)\ds(\yb)}\\
+\int_{\Gamma}\vnabla_{\yb} \Big(G(\kappa_{p},\xb-\yb)-G(\kappa_{s},\xb-\yb)\Big)\big(\nn(\yb)\cdot\vpsi(\yb)\big)\,\ds(\yb),
\end{multline}
and
\begin{multline}\label{Dp-ipp}
(D'_{\omega}\vphi)(\xb)=2\mu(\mathcal{M}S_{\omega}\vphi)(\xb)-{\,\mathcal{M}_{\xb}\int_{\Gamma}G(\kappa_{s},\xb-\yb)\vphi(\yb)\,\ds(\yb)} +\int_{\Gamma}\dfrac{\partial}{\partial\nn(\xb)}G(\kappa_{s},\xb-\yb)\vphi(\yb)\ds(\yb)\\
+\nn(\xb)\int_{\Gamma}\vphi(\yb)\cdot\vnabla_{\xb} \Big(G(\kappa_{p},\xb-\yb)-G(\kappa_{s},\xb-\yb)\Big)\,\ds(\yb).
\end{multline}
The treatment of the hypersingularity derived in~\cite{LeLouer1} for smooth boundaries extends to the Lipschitz case without any changes and we have:
\begin{equation}\label{N-ipp}
\begin{aligned}
N_{\omega}=2\mu D'_{\omega}\mathcal{M} +2\mu\mathcal{M}(D_{\omega}-2\mu S_{\omega}\mathcal{M})+C_{\omega},
\end{aligned}
\end{equation}
with
\begin{multline}\label{C3}
C_{\omega}\vpsi(\xb)=\rho\omega^2\,\nn(\xb)\times\int_{\Gamma}G(\kappa_{s},\xb-\yb)(\vpsi(\yb)\times\nn(\yb))\ds(\yb)-\mu\Rot_{\Gamma}\int_{\Gamma}G(\kappa_{s},\xb-\yb)\rot_{\Gamma}\vpsi(\yb) \ds(\yb)\\
+\rho\omega^2\nn(\xb)\int_{\Gamma}G(\kappa_{p},\xb-\yb)\big(\nn(\yb)\cdot\vpsi(\yb)\big)\ds(\yb), \qquad\text{ when $d=3$,}
\end{multline}
where we have $\Rot_{\Gamma}=-\,\nn\times\vnabla_{\Gamma}$ and $\rot_{\Gamma}=\Div_{\Gamma}(\,\cdot\times\nn)$\,, and
\begin{multline}\label{C2}
C_{\omega}\vpsi(\xb)=\rho\omega^2\,\nn^{\bot}(\xb)\int_{\Gamma}G(\kappa_{s},\xb-\yb)(\vpsi(\yb)\cdot\nn^{\bot}(\yb))\ds(\yb)\\
+\rho\omega^2\nn(\xb)\int_{\Gamma}G(\kappa_{p},\xb-\yb)\big(\nn(\yb)\cdot\vpsi(\yb)\big)\ds(\yb), \qquad\text{ when $d=2$.}
\end{multline}

In virtue of~\cite[Theorems~7.6 and 7.8]{McLean}, the boundary integral operators $S_{\omega}$ and $N_{\omega}$ satisfy G\r{a}rding's inequalities. By strong ellipticity of $(-\Updelta^*)$ (see~\cite[Chapter~10]{McLean}) we can formulate the following theorem

\begin{theo}\label{G\r{a}rdings}\ 
\begin{enumerate}\romanenumi
\item \label{5i} There exists a positive constant $c_{S}$ and a compact operator $L_{S}:\HH^{-\frac{1}{2}}(\Gamma)\to\HH^{\frac{1}{2}}(\Gamma)$ such that
\[
\Re\left(\langle S_{\omega}\vphi,\vphi\rangle_{\frac{1}{2}}+\langle L_{S}\vphi,\vphi\rangle_{\frac{1}{2}}\right)\ge c_{S}\|\vphi\|^2_{\HH^{-\frac{1}{2}}(\Gamma)},\qquad \text{ for all }\vphi\in\HH^{-\frac{1}{2}}(\Gamma).
\]
\item \label{5ii} There exists a positive constant $c_{N}$ and a compact operator $L_{N}:\HH^{\frac{1}{2}}(\Gamma)\to\HH^{-\frac{1}{2}}(\Gamma)$ such that
\[
\Re\left(\langle -N_{\omega}\vpsi,\vpsi\rangle_{\frac{1}{2}}+\langle L_{N}\vpsi,\vpsi\rangle_{\frac{1}{2}}\right)\ge c_{N}\|\vpsi\|^2_{\HH^{\frac{1}{2}}(\Gamma)},\qquad \text{ for all }\vpsi\in\HH^{\frac{1}{2}}(\Gamma).
\]
\end{enumerate}
\end{theo}

\begin{remark}
As for the Helmholtz and Maxwell equation cases~\cite{BuffaHiptmair,CostabelLeLouer,EnglederSteinbach1,EnglederSteinbach2,SteinbachWindisch}, the boundary integral operators arising from the elastodynamic case $\omega\not=0$ are compact perturbations of their counterpart arising from the static case $\omega=0$, even in the Lipschitz case. Therefore, when $d=3$, the first item for the operator $S_{\omega}=S_{0}+(S_{\omega}-S_{0})$ can also be justified by the ellipticity of $S_{0}$ (see~\cite[Theorem~10.7]{McLean} or~\cite{CostabelStephan2}). The two dimensional case seems more involved~\cite{Steinbach}.
\end{remark}

The operators $S_{\omega}$ and $N_{\omega}$ are then Fredholm operators of index zero. However, we have a loss of injectivity when $\rho\omega^2$ is either a Dirichlet or a Neumann eigenvalue of $(-\Updelta^*)$. Some Neumann eigenvalues for the unit square and cube are shown in Appendix~\ref{Ap-B}. Then combined field integral equation methods are usually considered to circumvent this issue. The following theorem established in~\cite[Theorem~2.11]{AgranovichAmosovLevitin} (see also~\cite[Theorems~3.17 and 3.22]{ColtonKress}) will be help full to get injectivity of the modified combined-field boundary integral equations presented in Section~\ref{MCFIE}.

\begin{theo}\label{Eigenspace}
For any $\omega>0$, we have
\[
\Ker(\tfrac{1}{2}\Id-D'_{\omega})=\left\{\gamma_{1}^-(\uu)\;;\; \uu\in\HH^{1}(\Omega),\; \Updelta^*\uu+\rho\omega^2\uu=0,\; \uu_{|\Gamma}=\zero\right\},
\]
\[
\Ker(\tfrac{1}{2}\Id+D_{\omega})=\left\{\gamma_{0}^-(\uu)\;;\; \uu\in\HH^{1}(\Omega),\; \Updelta^*\uu+\rho\omega^2\uu=0,\; \TT\uu=\zero\right\}.
\]
\end{theo}

In the next two Lemmas we extend~\cite[Lemma~3.1]{EnglederSteinbach1} also required for the injectivity analysis of the proposed boundary integral equation operators.

\begin{lemm}\label{Im-neg}
Let $\uu$ be a radiating solution to the Navier equation in $\Omega^c$. Then
\[
\lim_{R\to+\infty}\Imop\left(\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}\right)\ge0,
\]
in which the normal vector is taken outward the boundary ${|\xb|=R}$, namely it is $\hat{\xb}=\frac{\xb}{|\xb|}$.
\end{lemm}

\begin{proof}
Radiating solutions satisfy the same asymptotic behavior as the potential operators introduced in Lemma~\ref{potentials} so that their G\"unter deivative on the circle $|\xb|=R$ behaves as
\[
\mathcal{M}\uu(\xb)=O\left(\frac{1}{R^{\frac{d+1}{2}}}\right)\qquad\text{ as }R\to+\infty.
\]
It follows from~\eqref{Tu2} that
\[
\lim_{R\to+\infty}\int_{|\xb|=R}|\TT\uu-i\mathcal{K}_{\omega}\uu|^2=0\quad\Longleftrightarrow\;\; \lim_{R\to+\infty}\int_{|\xb|=R}\left|(\TT-2\mu\mathcal{M})\uu-i\mathcal{K}_{\omega}\uu\right|^2=0.
\]
Next, we restrict the computations to the two-dimensional case. The three-dimensional computations are identical~\cite{LeLouerRais}. Spliting the above integral according to the normal and tangential components of the traction derivative~\eqref{Tu2}, we get equivalently
\begin{equation}\label{Rellich}
\lim_{R\to+\infty}\int_{|\xb|=R}\lvert\Div\uu-i\kappa_{p}\hat{\xb}\cdot\uu|^2=0\quad \text{ and }\quad \lim_{R\to+\infty}\int_{|\xb|=R}\lvert\rot_{||}\uu-i\kappa_{s}\hat{\xb}^{\bot}\cdot\uu|^2=0.
\end{equation}
The symmetry property~\eqref{IPP4} of the G\"unter derivative leads to
\[
\int_{|\xb|=R}\mathcal{M}\uu\cdot\overline{\uu}=\int_{|\xb|=R}\uu\cdot\mathcal{M}\overline{\uu}=\overline{ \int_{|\xb|=R}\mathcal{M}\uu\cdot\overline{\uu}}\;,
\]
so that
\[
\begin{aligned}\Im\left(\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}\right)&=(\lambda+2\mu)\Im\left(\int_{|\xb|=R}(\Div\uu)(\hat{\xb}\cdot\overline{\uu})\right)+\mu\Im\left(\int_{|\xb|=R}(\rot_{||}\uu)(\hat{\xb}^{\bot}\cdot\overline{\uu})\right).
\end{aligned}
\]
Expending the square modulus in~\eqref{Rellich}, we finaly get
\begin{multline*}
\lim_{R\to+\infty}\Im\left(\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}\right)\\
=\lim_{R\to+\infty}\frac{\lambda+2\mu}{2\kappa_{p}}\int_{|\xb|=R}\lvert\Div\uu|^2+|\kappa_{p}(\hat{\xb}\cdot\uu)|^2+\frac{\mu}{2\kappa_{s}}\int_{|\xb|=R}\lvert\rot\uu|^2+|\kappa_{s}(\hat{\xb}^{\bot}\cdot\uu)|^2\ge0,
\end{multline*}
which ends the proof.
\end{proof}

\begin{lemm}\label{ImS-ImN}
For any $\vphi\in\HH^{-\frac{1}{2}}(\Gamma)$ and $\vpsi\in\HH^{\frac{1}{2}}(\Gamma)$, we have
\[
\Imop\langle S_{\omega}\vphi,\vphi\rangle_{\frac{1}{2}}\ge0\quad\text{ and }\quad \Imop\langle N_{\omega}\vpsi,\vpsi\rangle_{\frac{1}{2}}\ge0.
\]
\end{lemm}

\begin{proof}
%$\bullet$
Let $\vphi\in\HH^{-\frac{1}{2}}(\Gamma)$ and $\uu=\Psi^{\!S}_{\omega}\vphi$. We apply the first Green formula~\eqref{IPP1} to $\uu$ and $\vv=\overline{\uu}$ in $\Omega$ and in $B_{R}\backslash\Omega^c$ where $B_{R}$ is a sufficiently large ball of radius $R>0$ centered at origin and containing $\bar{\Omega}$. We have
\[
\int_{\Omega}\lambda\lvert\Div\uu|^2+2\mu|\upxi(\uu)|^2-\rho\omega^2\int_{\Omega}|\uu|^2= <\gamma_{1}^-(\uu),\gamma_{0}^-(\uu)\rangle_{\frac{1}{2}}
\]
and
\[
\int_{B_{R}\backslash\bar{\Omega}}\lambda\lvert\Div\uu|^2+2\mu|\upxi(\uu)|^2-\rho\omega^2\int_{B_{R}\backslash\bar{\Omega}}|\uu|^2=\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}-\langle\gamma_{1}^{+}(\uu),\gamma_{0}^{+}(\uu)\rangle_{\frac{1}{2}}.
\]
From~\eqref{calderon}--\eqref{jump} and adding the two equalities, we get
\begin{equation}\label{elliptic}
\begin{aligned} \int_{B_{R}}\lambda\lvert\Div\uu|^2+2\mu|\upxi(\uu)|^2-\rho\omega^2\int_{B_{R}}|\uu|^2&=\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}-\big\langle(\gamma_{1}^{+}-\gamma_{1}^{-})(\uu),\gamma_{0}^{+}(\uu)\big\rangle_{\frac{1}{2}}\\&=\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}-\big\langle(-\vphi),{S_{\omega}\vphi}\big\rangle_{\frac{1}{2}}.
\end{aligned}
\end{equation}
Taking the imaginary part and passing to the limit as $R\to+\infty$ we get from Lemma~\ref{Im-neg}
\[
\Imop\big\langle{S_{\omega}\vphi},\vphi\big\rangle_{\frac{1}{2}}=- \Imop\big\langle\vphi,{S_{\omega}\vphi}\big\rangle_{\frac{1}{2}}= \lim_{R\to+\infty}\Imop\left(\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}\right)\ge0.
\]

%$\bullet$
Let $\vpsi\in\HH^{\frac{1}{2}}(\Gamma)$. Substituing $\uu=\Psi^{\!D}_{\omega}\vpsi$ in the above equalities, we get in the same way
\[
\begin{aligned} \int_{B_{R}}\lambda\lvert\Div\uu|^2+2\mu|\upxi(\uu)|^2-\rho\omega^2\int_{B_{R}}|\uu|^2&=\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}-\big\langle \gamma_{1}^{+}(\uu),(\gamma_{0}^{+}-\gamma_{0}^{-})(\uu)\big\rangle_{\frac{1}{2}}\\&=\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}-\big\langle{N_{\omega}\vpsi},\vpsi\big\rangle_{\frac{1}{2}}.
\end{aligned}
\]
Taking the imaginary part and passing to the limit as $R\to+\infty$ we get from Lemma~\ref{Im-neg}
\[
\Imop\big\langle{N_{\omega}\vpsi},\vpsi\big\rangle_{\frac{1}{2}}= \lim_{R\to+\infty}\Imop\left(\int_{|\xb|=R}\TT\uu\cdot\overline{\uu}\right)\ge0.\qedhere
\]
\end{proof}


\section{The modified field boundary integral equations}\label{MCFIE}

For smooth boundaries, uniquely solvable boundary integral equations equivalent to the scattering problems~\eqref{NE}--\eqref{dirichlet}--\eqref{radiation} and~\eqref{NE}--\eqref{neumann}--\eqref{radiation} rely on the following layer ansatz for any radiating solution to the Navier equation:
\[
\uu(\xb)=a(\Psi_{\omega}^{\!D}\vphi)(\xb)-b(\Psi_{\omega}^{\!S}\vphi)(\xb),\qquad \text{ for all }\xb\in\Omega^c,
\]
where $a$, $b$ are non zero complex constants that have to be chosen such that the associated interior boundary value problem:
\begin{equation}\label{IBVP}
\begin{aligned}
\Updelta^{\!*}\uu+\rho\omega^2\uu&=0&&\text{in }\Omega\\
a\gamma_{1}^-\uu-b\gamma_{0}^-\uu&=0&&\text{on }\Gamma,
\end{aligned}
\end{equation}
do not admit any spurious frequencies $\rho\omega^2$. Choosing $\vv=\overline{\uu}$ in the first Green's formula~\eqref{Green} leads to the choice $\frac{b}{a}=i\eta$ where $\eta$ is a non zero real number for positive Lam\'e parameters and $\eta>0$ for complex valued Lamé parameters (viscoelastic model with $e^{-i\omega t}$ time dependence). As for the acoustic case, one can even replace the scalar number $\eta$ by the matrix $\mathcal{K}_{\omega}$ given in the radiation condition~\eqref{radiation} for an improved eigenvalue clustering for an application of the GMRES solver~\cite{DarbasLeLouer}. Even with the lack of compactness of the operators $D_{\omega}$ and its adjoint $D'_{\omega}$, on can prove that the resulting boundary integral equations discussed in~\cite{LeLouer1} are Fredholm of index zero.

In the case of Lipschitz boundaries, we have to \emph{regularize} or \emph{modify} either the single or the double layer potential operator. Extending earlier work in the acoustic case~\cite{BuffaHiptmair,EnglederSteinbach1,EnglederSteinbach2}, we present various approaches for deriving well-posed boundary integral equations in linear elasticity.

\subsection{Positivity results for the Laplace equation case}

In this paragraph we introduce the symmetric positive definite boundary integral operators arising from the Laplace equation studies that may define a norm on $\HH^{\pm\frac{1}{2}}(\Gamma)$. More precisely we vectorize several results presented in~\cite{EnglederSteinbach1,EnglederSteinbach2,NedelecPlanchard,SteinbachOf} and~\cite[Chapter~8]{McLean}.

A fundamental solution for the Laplace equation is
\[
G(\kappa=0,\zb)=
\begin{cases}
\frac{1}{2\pi}\log\frac{r}{|\zb|}&\text{ when }d=2\text{ for any }r>0,\\
\frac{1}{4\pi|\zb|}&\text{ when }d=3.
\end{cases}
\]
Let us introduce first the single layer operator
\begin{equation}\label{V0}
\begin{aligned}
\VV_{\!0}:\HH^{-\frac{1}{2}}(\Gamma)&\longrightarrow \HH^{\frac{1}{2}}(\Gamma)\\
\vphi&\longmapsto {\left\{\xb\mapsto\int_{\Gamma}G(0,\xb-\yb)\vphi(\yb)\,\ds(\yb)\right\}},
\end{aligned}
\end{equation}

\begin{theo}\ 
\begin{enumerate}\romanenumi
\item \label{10i} When $d=2$, under the asumption $r>\Cap_{\Gamma}$ the operator $\VV_{\!0}$ is positive and bounded below. Then it admits a bounded inverse and there exists a positive constant $c^2_{V}$ such that
\[
\langle \VV_{\!0}^{-1}\vpsi,\vpsi\rangle_{\frac{1}{2}}\ge c^2_{V}\|\vpsi\|^2_{\HH^{\frac{1}{2}}(\Gamma)} \qquad \text{ for all }\vpsi\in\HH^{\frac{1}{2}}(\Gamma).
\]
The logarithmic capacity $\Cap_{\Gamma}$ is defined in~\cite[p.~262]{McLean}.

\item \label{10ii} When $d=3$, the operator $\VV_{\!0}$ is positive and bounded below. Then it admits a bounded inverse and there exists a positive constant $c^3_{V}$ such that
\[
\langle \VV_{\!0}^{-1}\vpsi,\vpsi\rangle_{\frac{1}{2}}\ge c^3_{V}\|\vpsi\|^2_{\HH^{\frac{1}{2}}(\Gamma)} \qquad \text{ for all }\vpsi\in\HH^{\frac{1}{2}}(\Gamma).
\]
\end{enumerate}
\end{theo}

Now we turn to the hypersingular operator for the three dimensional case~\cite{SteinbachOf}. Let consider
\[
\begin{aligned}
\WW_{\!0}:\HH^{\frac{1}{2}}(\Gamma)&\longrightarrow \HH^{-\frac{1}{2}}(\Gamma)\\
\vpsi&\longmapsto {\left\{\xb\mapsto\int_{\Gamma}\frac{\partial^2G(0,\xb-\yb)}{\partial\nn(\xb)\partial\nn(\yb)}\vpsi(\yb)\,\ds(\yb)\right\}}.
\end{aligned}
\]
We perturbe $-\WW_{0}$ as follows
\begin{equation}\label{Well}
\widetilde{\WW_{0}}\vpsi=\big(\langle\psi_{i},1\rangle_{\frac{1}{2}}\big)_{1\le i\le 3}-\WW_{0}\vpsi\qquad\text{ where }\vpsi=\transposee{(\psi_{1},\psi_{2},\psi_{3})}.
\end{equation}

\begin{theo}
The operator $\widetilde{\WW_{0}}$ defines an isomorphism from $\HH^{\frac{1}{2}}(\Gamma)$ to $\HH^{-\frac{1}{2}}(\Gamma)$. Moreover $\widetilde{\WW_{0}}$ is positive and bounded below. Then it admits a bounded inverse and there exists a positive constant $c_{W}$ such that
\[
\langle \widetilde{\WW_{0}}^{-1}\vphi,\vphi\rangle_{\frac{1}{2}}\ge c_{W}\|\vphi\|^2_{\HH^{-\frac{1}{2}}(\Gamma)}, \qquad \text{ for all }\vphi\in\HH^{-\frac{1}{2}}(\Gamma).
\]
\end{theo}

\subsection{The Dirichlet boundary condition case}

For an equivalent reduction of the volume problem~\eqref{NE}--\eqref{dirichlet}--\eqref{radiation} as a well-posed integral equation on the Lipschitz boundary $\Gamma$ one can use the layer ansatz
\begin{equation}\label{ansatz1}
\uu(\xb)=i\eta\big(\Psi_{\omega}^{\!D}(\boldsymbol{\mathcal{R}}\vphi)\big)(\xb)+(\Psi_{\omega}^{\!S}\vphi)(\xb),\qquad \text{ for all }\xb\in\Omega^c,
\end{equation}
where the operator $\mathcal{R}$ might be chosen as follows
\begin{alignat}{2}
\mathcal{R}&=(\Id-\Delta_{\Gamma})^{-1}, && \text{when $d=2$ or $3$},
\\
\mathcal{R}&= \widetilde{\WW}_0^{-1}(\tfrac{1}{2}\Id+D'_{-\omega}), \quad && \text{when $d=3$}.
\end{alignat}
The first choice is proposed in~\cite{BuffaHiptmair} and here we consider the two- and three-dimensional modified Laplace--Beltrami operator variationally defined as
\[
\int_{\Gamma} \big((\Id-\Delta_{\Gamma})\vphi\big)\cdot\overline{\tilde{\vphi}}=\int_{\Gamma}\vphi\cdot\overline{\tilde{\vphi}}+\int_{\Gamma}[\nabla_{\Gamma}\vphi]:[\nabla_{\Gamma}\overline{\tilde{\vphi}}],\qquad\text{ for all }\tilde{\vphi}\in\HH^{1}(\Gamma).
\]
This operator is a symmetric definite positive operator bounded from $\HH^{1}(\Gamma)$ to $\HH^{-1}(\Gamma)$ and its inverse has a regularizing property since it is bounded from $\HH^{-\frac{1}{2}}(\Gamma)$ to $\HH^{1}(\Gamma)$ and is thus compact from $\HH^{-\frac{1}{2}}(\Gamma)$ to $\HH^{\frac{1}{2}}(\Gamma)$.
The second choice is a Steklov--Poincar\'e type operator bounded from $\HH^{-\frac{1}{2}}(\Gamma)$ to $\HH^{\frac{1}{2}}(\Gamma)$. In all cases, we deduce the following theorem

\begin{theo}
For all frequencies $\omega>0$, the solution to the Dirichlet scattering problem can be equivalently reduced to the solution to the boundary integral equation
\begin{equation}\label{Iop1}
S_\omega\vphi+i\eta\left(\tfrac{1}{2}\Id +D_\omega\right)\mathcal{R}\vphi=-\uu^{\inc}, \quad\text{for any\/ }\eta>0.
\end{equation}
More precisely, the operator $\left[S_\omega +i\eta(\tfrac{1}{2}\Id +D_\omega)\mathcal{R}\right]:\HH^{-\frac{1}{2}}(\Gamma)\to\HH^{\frac{1}{2}}(\Gamma)$ satisfies a G\r{a}rding's inequality and is injective.
\end{theo}

\begin{proof}
\begin{proof}[(i)]
Let us start with $\mathcal{R}=(\Id-\Delta_{\Gamma})^{-1}$.

%$\bullet$
The compactness of $\mathcal{R}$ yieds the compactness of $(\tfrac{1}{2}\Id +D_\omega)\mathcal{R}:\HH^{-\frac{1}{2}}(\Gamma)\to\HH^{\frac{1}{2}}(\Gamma)$. In virtue of the first item \eqref{5i} in Theorem~\ref{G\r{a}rdings}, the boundary integral equation operator $\left[S_\omega +i\eta(\tfrac{1}{2}\Id +D_\omega)\mathcal{R}\right]$ satisfies a G\r{a}rding's inequality.

%$\bullet$
Now let us assume that $\uu^{\inc}=0$. The uniqueness of the solution to the exterior scattering problem~\eqref{NE}--\eqref{dirichlet}--\eqref{radiation} ensure that the integral representation~\eqref{ansatz1} vanishes in $\Omega^c$ so that the jump relations~\eqref{jump} gives
\[
\gamma_0^-\uu(x)=-i\eta\mathcal{R}\vphi\quad\text{ and } \gamma_1^-\uu(x)=\vphi.
\]
Then the ansatz~\eqref{ansatz1} solves the interior problem
\[
\begin{aligned}
\Updelta^{\!*}\uu+\rho\omega^2\uu&=0 &&\text{in }\Omega\\
\gamma_{1}^-\uu+[i\eta\mathcal{R}]^{-1}\gamma_{0}^-\uu&=0 &&\text{on }\Gamma.
\end{aligned}
\]
The positivity of modified the Laplace--Beltrami operator ensures that the above impedance-like interior problem do not admit any real eigenvalues $\rho\omega^2$. Then $\uu=0$ in $\Omega$ and $\vphi=0$ on $\Gamma$. We deduce the injectivity of $\left[S_\omega +i\eta(\tfrac{1}{2}\Id +D_\omega)\mathcal{R}\right]$.
\let\qed\relax
\end{proof}

\begin{proof}[(ii)] Now let us consider the second choice $\mathcal{R}=\widetilde{\WW}_0^{-1}(\Id+D'_\omega)$.

%$\bullet$
We have
\[
\begin{aligned}
\Re\left\langle \left[S_\omega +i\eta(\tfrac{1}{2}\Id +D_\omega)\mathcal{R}\right]\vphi,\vphi\right\rangle_{\frac{1}{2}}&=\Re\left(\langle S_\omega \vphi,\vphi\rangle_{\frac{1}{2}} +i\eta\left\langle\widetilde{\WW}_0^{-1}(\tfrac{1}{2}\Id +D'_{-\omega})\vphi,(\tfrac{1}{2}\Id +D'_{-\omega})\vphi\right\rangle_{\frac{1}{2}}\right)\\
&=\Re\langle S_\omega\vphi,\vphi\rangle_{\frac{1}{2}}.
\end{aligned}.
\]
In virtue of the first item \eqref{5i} in Theorem~\ref{G\r{a}rdings}, the boundary integral equation operator $\left[S_\omega +i\eta(\tfrac{1}{2}\Id +D_\omega)\mathcal{R}\right]$ satisfies again a G\r{a}rding's inequality then it is Fredholm of index zero.

%$\bullet$
Now let us assume that $\uu^{\inc}=0$. Taking the imaginary part of~\eqref{Iop1} we obtain
\[
\Imop\left(\langle S_\omega\vphi,\vphi\rangle_{\frac{1}{2}}\right)+ \eta\left\langle \widetilde{\WW}_0^{-1}(\tfrac{1}{2}\Id +D'_{-\omega})\vphi,(\tfrac{1}{2}\Id +D'_{-\omega})\vphi\right\rangle_{\frac{1}{2}}=0.
\]
From Lemma~\ref{ImS-ImN} and the positivity of $\widetilde{\WW}_0^{-1}$, both of the two above terms vanish and we have
\[
(\tfrac{1}{2}\Id +D'_{-\omega})\vphi=\zero.
\]
Coming back to the BIE~\eqref{Iop1} we also get $S_\omega\vphi=\zero$. We deduce that the potential $\Psi_{\omega}^{\!S}\vphi$ is a Dirichlet eigenfunction of the Navier equation. In virtue of Theorem~\ref{Eigenspace}, since $-\omega$ leads to the same Navier equation then we also have
\[
(\tfrac{1}{2}\Id -D'_{-\omega})\vphi=\zero,
\]
then $\vphi=\zero$. We deduce again the injectivity of $\left[S_\omega +i\eta(\tfrac{1}{2} \Id+D_\omega)\mathcal{R}\right]$.
\end{proof}
\let\qed\relax
\end{proof}

\subsection{The Neumann boundary condition case}

For an equivalent reduction of the volume problem~\eqref{NE}-\eqref{neumann}-\eqref{radiation} as a well-posed integral equation on the Lipschitz boundary $\Gamma$ one can use the layer ansatz
\begin{equation}\label{ansatz1b}
\uu(\xb)=\big(\Psi_{\omega}^{\!D}\vpsi\big)(\xb)-i\eta(\Psi_{\omega}^{\!S}(\boldsymbol{\mathcal{R}}\vpsi))(\xb),\qquad \text{ for all }\xb\in\Omega^c,
\end{equation}
where
\begin{equation}
\mathcal{R}= {\VV}_0^{-1}(\tfrac{1}{2}\Id-D_{-\omega}), \qquad \text{when $d=2$ or $3$},
\end{equation}
which is a Steklov--Poincar\'e type operator too and is bounded from $\HH^{\frac{1}{2}}(\Gamma)$ to $\HH^{-\frac{1}{2}}(\Gamma)$.

\begin{theo}
For all frequencies $\omega>0$, the solution to the Neumann scattering problem can be equivalently reduced to the solution to the boundary integral equation
\begin{equation}\label{Iop2}
N_\omega\vpsi+i\eta\left(\tfrac{1}{2}\Id -D'_\omega\right)\mathcal{R}\vpsi=-\TT\uu^{\inc}, \quad\text{for any }\eta>0.
\end{equation}
More precisely, the operator $\left[N_\omega +i\eta(\tfrac{1}{2}\Id -D'_\omega)\mathcal{R}\right]:\HH^{-\frac{1}{2}}(\Gamma)\to\HH^{\frac{1}{2}}(\Gamma)$ satisfies a G\r{a}rding's inequality and is injective.
\end{theo}

\begin{proof}
%$\bullet$
We have
\[
\begin{aligned}
\Re\left\langle \left[N_\omega +i\eta(\tfrac{1}{2}\Id -D'_\omega)\mathcal{R}\right]\vpsi,\vpsi\right\rangle_{\frac{1}{2}}&=\Re\left(\langle N_\omega \vpsi,\vpsi\rangle_{\frac{1}{2}} +i\eta\left\langle\VV_0^{-1}(\tfrac{1}{2}\Id -D_{-\omega})\vpsi,(\tfrac{1}{2}\Id -D_{-\omega})\vpsi\right\rangle_{\frac{1}{2}}\right)\\
&=\Re\left\langle N_\omega\vpsi,\vpsi\right\rangle_{\frac{1}{2}}.
\end{aligned}.
\]
In virtue of the second item \eqref{5ii} in Theorem~\ref{G\r{a}rdings}, the boundary integral equation operator $\left[N_\omega +i\eta(\tfrac{1}{2}\Id -D'_\omega)\mathcal{R}\right]$ satisfies a G\r{a}rding's inequality then it is Fredholm of index zero.

%$\bullet$
Now let us assume that $\uu^{\inc}=0$. Taking the imaginary part of~\eqref{Iop2} we obtain
\[
\Imop\left(\langle N_\omega\vpsi,\vpsi\rangle_{\frac{1}{2}}\right)+ \eta\langle \VV_0^{-1}(\tfrac{1}{2}\Id -D_{-\omega})\vpsi,(\tfrac{1}{2}\Id -D_{-\omega})\vpsi\rangle_{\frac{1}{2}}=0.
\]
From Lemma~\ref{ImS-ImN} and the positivity of $\VV_0^{-1}$, both of the two above terms vanish and we have
\[
(\tfrac{1}{2}\Id -D_{-\omega})\vpsi=\zero.
\]
Coming back to the BIE~\eqref{Iop2} we also get $N_\omega\vpsi=\zero$. We deduce that the potential $\Psi_{\omega}^{\!D}\vpsi$ is a Neumann eigenfunction of the Navier equation. In virtue of Theorem~\ref{Eigenspace}, we also have
\[
(\tfrac{1}{2}\Id +D_{-\omega})\vpsi=\zero,
\]
then $\vpsi=\zero$. We deduce again the injectivity of $\left[N_\omega +i\eta(\tfrac{1}{2}\Id -D'_\omega)\mathcal{R}\right]$.
\end{proof}


\section{Conclusion}

The aim of establishing G\r{a}rding's inequalities is two fold. On the one hand, one obtains existence proofs for the solution to the modified combined field boundary integral equations, and consequently for the equivalent exterior elliptic boundary value problems. On the other hand it guarantees stability of the aproximate Galerkin scheme for any choice of finite dimensional subspaces of the Sobolev spaces $\HH^{\pm\frac{1}{2}}(\Gamma)$. A numerical comparison of the proposed formulations~\eqref{Iop1} and~\eqref{Iop2} with the classical ones discussed in~\cite{LeLouer1} is the subject of a forthcoming work.

For readers interested in the analysis of the elastodynamic transmission problem we refer to~\cite{LeLouer2} for smooth boundaries and~\cite{CostabelStephan1} for the case of Lipschitz boundaries. In both cases, the volume problem is reduced to a couple of boundary integral equations for the unknown Dirichlet and Neumann traces. However, it is possible to extend the approach of combined field integral equations for a single unknown to transmission problems that also satisfies G\r{a}rding's inequalities even for Lipschitz boundaries~\cite{CostabelLeLouer}. In this work, a key point relied on the fact that the double layer boundary integral operators associated to two different wavenumbers admit the same principal part. The principal symbols obtained in Appendix~\ref{Ap-A} reveal that this property is not valid anymore and the extension of the work to elasticity is not immediate.

\appendix
\setcounter{equation}{0}
\renewcommand\theequation{A\arabic{equation}}
\section{Pseudo-differential calculus for the Navier system}\label{Ap-A}

The principal symbols of the single and the double layer boundary integral operators for the Lam\'e system (i.e. $\omega=0$) where obtained in~\cite{AgranovichAmosovLevitin}. In this section, we extend the computations to the dynamic case in $\R^3$ (i.e. $\omega\not=0$).

At a fixed regular point $\xb=(x_{1},x_{2},x_{3})\in\Gamma$, there exists a neighborhood $\mathcal{O}\subset\R^3$ of $\xb$ such that the portion of surface $\Gamma\cap\mathcal{O}$ can be defined by the explicit equation $y_{3}=\phi(y_{1},y_{2})$. We translate the origin of the orthonormal system of $\R^3$ to the point $\xb$ by setting the new coordinates
\[
(\tilde{y}_{1},\tilde{y}_{2},\tilde{y}_{3})=\psi^{-1}\big(y_{1}-x_{1},y_{2}-x_{2},y_{3}-\phi(y_{1},y_{2})\big),
\]
where $\psi^{-1}$ rotates the axis in such a way that $\tilde{y}_{3}>0$ indicates the direction of the exterior normal to $\Gamma$ and $\tilde{y_{3}}=0$ indicates the tangent plane to $\Gamma$ at the point $\xb$. In this new coordinate system, the tangential direction is $\tilde{\yb}_{\|} := \transposee{(\tilde{y}_{1}, \tilde{y}_{2})}$ and the outward normal vector to $\Gamma$ is $\transposee{(0,0,1)}$. The associated Fourier variable of $\tilde{\yb}_{||}$ is $\ki_{||} := \transposee{(\xi_{1}, \xi_{2})}$ and we set $\|\ki\| := \sqrt{\ki\cdot\ki}=(|\ki_{||}|^2+|\xi_{3}|^2)^{\frac{1}{2}}$.

\begin{prop}
For all $\omega>0$, the principal symbol of the single layer operator is given by
\[
\sigma_{S}(\ki_{||})=\frac{1}{2\rho\omega^2}\left(\frac{1}{\big|{|\ki_{_{||}}|^2-\kappa_{p}^2}\big|^{\frac{1}{2}}}
\begin{pmatrix}\ki_{_{||}}\transposee{\ki_{_{||}}}&0\\0&{\kappa_{p}^2-|\ki_{_{||}}|^2}
\end{pmatrix}+\frac{1}{\big|{|\ki_{_{||}}|^2-\kappa_{s}^2}\big|^{\frac{1}{2}}}
\begin{pmatrix}\kappa_{s}^2\Id_{_{||}}-\ki_{_{||}}\transposee{\ki_{_{||}}}&0\\0&{|\ki_{_{||}}|^2}
\end{pmatrix}\right),
\]
where $\Id_{||}$ is the $(2\times2)$ identity matrix.
\end{prop}

\begin{proof}
We first compute the symbol of the volume potential with kernel $\Phi_{\omega}(\xb-\,\cdot\,)$ where $\Phi_{\omega}$ is the fundamental solution to the Navier system : $\Updelta^*\Phi_{\omega}+\rho\omega^2\Phi_{\omega}=-\delta_{0}\Id$.
\[
\Delta=\vnabla\Div-\Rot\Rot\quad \Longrightarrow\quad { \Updelta^*+\rho\omega^2\Id=(\lambda+2\mu)\vnabla\Div-\mu\Rot\Rot+\rho\omega^2\Id}.
\]
The Navier operator has the symbol $\ell(\ki)=-(\lambda+2\mu)\ki\transposee{\ki}-\mu(|\ki|^2-\ki\transposee{\ki})+\rho\omega^2\Id$ \;, that can be rewritten:
\[
{ \ell(\ki)=(\lambda+2\mu)\left(\kappa_{p}^2-|\ki|^2\right)\frac{1}{|\ki|^2}\ki\transposee{\ki}+\mu\left(\kappa_{s}^2-|\ki|^2\right)\left(\Id-\frac{1}{|\ki|^2}\ki\transposee{\ki}\right)}.
\]
We deduce that the volume potential with kernel ${\Phi_{\omega}(\xb-\,\cdot\,)}$ has the symbol
\[
{ -\ell^{-1}(\ki)=(\lambda+2\mu)^{-1}\frac{1}{|\ki|^2-\kappa_{p}^2}\frac{1}{|\ki|^2}\ki\transposee{\ki}+\mu^{-1}\frac{1}{|\ki|^2-\kappa_{s}^2}\left(\Id-\frac{1}{|\ki|^2}\ki\transposee{\ki}\right)}.
\]
The boundary integral operator with kernel ${\Phi_{\omega}(\xb-\,\cdot\,)}$ has the principal symbol
\begin{equation}\label{IFT}
{ \sigma_{S}(\ki_{||})={\frac{1}{2\pi}\int_{-\infty}^{\infty}-\ell^{-1}(\ki)\dxi_{3}}}\qquad\text{where $\ki=(\xi_{1},\xi_{2},\xi_{3})=(\ki_{_{||}},\xi_{3})$}.
\end{equation}
We have
\[
\ki\transposee{\ki}=
\begin{pmatrix}\ki_{_{||}}\transposee{\ki_{_{||}}}&\ki_{_{||}}\xi_{3}\\\xi_{3}\transposee{\ki_{_{||}}}&|\xi_{3}|^2
\end{pmatrix}.
\]
The non diagonal terms are odd so that their partial inverse Fourier transform~\eqref{IFT} vanishes. For the diagonal term associated to $\ki_{||}\transposee{\ki_{||}}$, we simplify
\begin{multline*}
\frac{1}{|\ki|^2-\kappa^2}\frac{1}{|\ki|^2}=\frac{1}{\kappa^2}\left(\frac{1}{|\ki|^2-\kappa^2}-\frac{1}{|\ki|^2}\right)\\
\Longrightarrow\quad { \bigl[-\ell^{-1}(\ki)\Big]_{||,||}=\frac{1}{\rho\omega^2}\left(\frac{1}{|\ki|^2-\kappa_{p}^2}\ki_{||}\transposee{\ki_{||}}+\frac{1}{|\ki|^2-\kappa_{s}^2}\left(\kappa_{s}^2\Id_{||}-\ki_{||}\transposee{\ki_{||}}\right)\right)}
\end{multline*}
and we use the well-known result from the Helmholtz equation case:
\[
\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{1}{|\ki|^2-\kappa^2}\dxi_{3}=\frac{1}{2\big||\ki_{||}|^2-\kappa^2\big|^{\frac{1}{2}}}.
\]
For the last diagonal term associated to $|\xi_{3}|^2$, we simplify
\[
\frac{1}{|\ki|^2-\kappa^2}\frac{|\xi_{3}|^2}{|\ki|^2}=\frac{-1}{\kappa^2}\left(\frac{|\ki_{||}|^2-\kappa^2}{|\ki|^2-\kappa^2}-\frac{|\ki_{||}|^2}{|\ki|^2}\right)\quad\Longrightarrow\quad {\bigl[-\ell^{-1}(\ki)\Big]_{3,3}=\frac{1}{\rho\omega^2}\left(\frac{\kappa_{p}^2-|\ki_{||}|^2}{|\ki|^2-\kappa_{p}^2}+\frac{|\ki_{||}|^2}{|\ki|^2-\kappa_{s}^2}\right)}.\qedhere
\]
\end{proof}

The principal symbol of the tangential G\"unter derivative $\mathcal{M}$ is
\[
\sigma_{\!\!\mathcal M}(\ki_{||})=
\begin{pmatrix}0&i\ki_{||}\\-i\transposee{\ki_{||}}&0
\end{pmatrix}.
\]
One can derive the principal symbols for the three other operators $D_{\omega}$, $D'_{\omega}$ and $N_{\omega}$ using product rules in the integral representations~\eqref{D-ipp}--\eqref{Dp-ipp}--\eqref{N-ipp} and the well-known principal symbols of the acoustic and electromagnetic boundary integral operators.

For a fixed frequency $\omega$ and fixed Lam\'e parameters $\lambda$ and $\mu$, by using the following asymptotic expansion when $|\ki_{||}|\to+\infty$,
\[
\frac{1}{(|\ki_{||}|^2-\kappa^2)^{\frac{1}{2}}}=\frac{1}{|\ki_{||}|}\left(1+\frac{\kappa^2}{2|\ki_{||}|^2}+o\big(\frac{1}{|\ki_{||}|^2}\big) \right),
\]
we retrieve the principal symbols of the four boundary integral operators associated to the static case $\omega=0$ that are in agreement with the results presented in~\cite[Proposition~4.1]{DarbasLeLouer}:
\[
\sigma_{\!S_{0}}(\ki_{||})=
\begin{dcases}
\frac{1}{2\mu|\ki_{||}|}
\begin{pmatrix}\Id_{||}-\tfrac{\ki_{||}\transposee{\ki_{||}}}{|\ki_{||}|^2}&&0\\0&&0
\end{pmatrix}+\frac{\lambda+3\mu}{4\mu(\lambda+2\mu)|\ki_{||}|}
\begin{pmatrix}\tfrac{\ki_{||}\transposee{\ki_{||}}}{|\ki_{||}|^2}&0\\0&1
\end{pmatrix} & \text{when $d=3$},\\
\frac{\lambda+3\mu}{4\mu(\lambda+2\mu)|\ki_{||}|}\Id &\text{when $d=2$,}
\end{dcases}
\]
\[
\begin{aligned} \sigma_{\!D_{0}}(\ki_{||})=\left(2\mu\sigma_{\!S_{0}}(\ki_{||})-\frac{1}{2|\ki_{||}|}\right)\sigma_{\!\!\mathcal M}(\ki_{||})&=\frac{\mu}{2(\lambda+2\mu)|\ki_{||}|}\sigma_{\!\!\mathcal M}(\ki_{||}),
\end{aligned}
\]
and
\[
\sigma_{\!N_{0}}(\ki_{||})=
\begin{dcases}
 -\frac{\mu|\ki_{||}|}{2}
\begin{pmatrix}\Id_{||}-\tfrac{\ki_{||}\transposee{\ki_{||}}}{|\ki_{||}|^2}&&0\\0&&0
\end{pmatrix}-\frac{\mu(\lambda+\mu)|\ki_{||}|}{\lambda+2\mu}
\begin{pmatrix}\tfrac{\ki_{||}\transposee{\ki_{||}}}{|\ki_{||}|^2}&0\\0&1
\end{pmatrix} &\text{when $d=3$},\\
=-\frac{\mu(\lambda+\mu)|\ki_{||}|}{\lambda+2\mu}\Id &\text{when $d=2$.}
\end{dcases}
\]

\renewcommand\theequation{B\arabic{equation}}
\section{Some Neumann interior eigenvalues for elastic square and cube}\label{Ap-B}

For forthcoming numerical investigations, we provide in this appendix some eigenvalues and the associated eigenfunctions for the Neumann interior elastodynamic problem inspired by the work presented in~\cite{GrebenkovNguyen}.

\begin{prop}
Let $\Omega=(0,1)^d$. For any $n\in\N^*$, some Neumann eigenvalues are given by $\rho\omega^2=2\mu n^2\pi^2$ associated to the eigenvector field
\[
\begin{aligned}
\uu_{s}(x_{1},x_{2})&=\Rot_{||}\big(\cos(n\pi x_{1})\cos(n\pi x_{2})\big), &\text{ for }d=2.
\\
\uu_{s}(x_{1},x_{2},x_{3})&=
\begin{pmatrix}\Rot_{||}\big(\cos(n\pi x_{1})\cos(n\pi x_{2})\big)\\0
\end{pmatrix}, &\text{ for }d=3,
\end{aligned}
\]
where $\Rot_{||}=\transposee{(\partial_{x_{2}},-\partial_{x_{1}})}$.
\end{prop}

\begin{proof}
\begin{proof}[(i)]
We start by the two-dimensional case.

%$\bullet$
Let $\uu_{s}=\transposee{(u_{s}^1,u_{s}^2)}$ be a share wave (i.e. $\Div\uu_{s}=0$) and $\uu_{p}=\transposee{(u_{p}^1,u_{p}^2)}$ be a pressure wave (i.e. $\rot_{||}\uu_{p}=0$). Then the boundary condition $\TT(\uu_{s}+\uu_{p})=0$ takes the form
\[
\left\{
\begin{aligned}
(\lambda+2\mu)\partial_{x_{1}}u_{p}^1+\lambda\partial_{x_{2}}u_{p}^2+2\mu\partial_{x_{1}}u_{s}^1&=0\\ 2\mu\partial_{x_{1}}u_{p}^2+\mu(\partial_{x_{1}}u_{s}^2+\partial_{x_{2}}u_{s}^1)&=0
\end{aligned}\right.
\quad \text{on }\{0;1\}\times[0,1],
\]
and
\[
\left\{
\begin{aligned}
(\lambda+2\mu)\partial_{x_{2}}u_{p}^2+\lambda\partial_{x_{1}}u_{p}^1+2\mu\partial_{x_{2}}u_{s}^2&=0\\ -2\mu\partial_{x_{2}}u_{p}^1-\mu(\partial_{x_{1}}u_{s}^2+\partial_{x_{2}}u_{s}^1)&=0
\end{aligned}\right.
\quad \text{on }[0,1]\times \{0;1\}.
\]
One can check that $\uu_{p}=\zero$ and
\[
\uu_{s}(x_{1},x_{2})=\Rot_{||}\big(\cos(n\pi x_{1})\cos(n\pi x_{2})\big)=
\begin{pmatrix}-\cos(n\pi x_{1})\sin(n\pi x_{2})\\\sin(n\pi x_{1})\cos(n\pi x_{2})
\end{pmatrix}
\]
solves the two systems and $-\Updelta^*\uu_{s}=\mu\Rot_{||}\big(-\Delta(\cos(n\pi x_{1})\cos(n\pi x_{2}))\big)=\mu (2n^2\pi^2)\uu_{s}$.

%$\bullet$
The converse, that is $\uu_{s}=0$ and
\[
\uu_{p}(x_{1},x_{2})=\nabla\big(\phi_{1}(x_{1})\phi_{2}(x_{2})\big)=
\begin{pmatrix}\phi'_{1}(x_{1})\phi_{2}(x_{2})\\\phi_{1}(x_{1})\phi'_{2}(x_{2})
\end{pmatrix},
\]
where $\phi_{i}(x_{i})=\cos(n_{i}\pi x_{i})$ or $\phi_{i}(x_{i})=\sin(n_{i}\pi x_{i})$ for $i=1,2$ does not solve the homogeneous Neumann boundary condition.
\let\qed\relax
\end{proof}

\begin{proof}[(ii)]
Now we turn to the three-dimensional case. Let consider
\[
\uu_{s}(x_{1},x_{2},x_{3})=
\begin{pmatrix}\Rot_{||}\big(\cos(n\pi x_{1})\cos(n\pi x_{2})\big)\\0
\end{pmatrix}\phi(x_{3})
\]
be a share wave (i.e. $\Div\uu_{s}=0$). We have
\[
\Rot\uu_{s}(x_{1},x_{2},x_{3})=
\begin{pmatrix}\phi'(x_{3})\nabla\big(\cos(n\pi x_{1})\cos(n\pi x_{2})\big)\\-\Delta\big(\cos(n\pi x_{1})\cos(n\pi x_{2})\big)\phi(x_{3})
\end{pmatrix}.
\]
Then the boundary condition $\TT\uu_{s}=2\mu\frac{\partial \uu_{s}}{\partial\nn}+\mu\nn\times\Rot\uu_{s}=0$ leads to the additional equations on the third component
\[
\sin(n\pi x_{2})\phi'(x_{3})=0\quad \text{when $x_{1}=0$ or $1$ and }(x_{2},x_{3})\in[0,1]^2,
\]
and
\[
\sin(n\pi x_{1})\phi'(x_{3})=0\quad \text{when $x_{2}=0$ or $1$ and }(x_{1},x_{3})\in[0,1]^2,
\]
and on the first two components
\[
\left\{
\begin{aligned}
\cos(n\pi x_{1})\sin(n\pi x_{2})\phi'(x_{3})&=0\\
\sin(n\pi x_{1})\cos(n\pi x_{2})\phi'(x_{3})&=0
\end{aligned}\right.
\quad \text{when $x_{3}=0$ or $1$ and }(x_{1},x_{2})\in[0,1]^2.
\]
The choice $\phi(x_{3})=\cos(n_{3}\pi x_{3})$ cancels the last system. However, the first two equations require that $n_{3}=0$.
\end{proof}
\let\qed\relax
\end{proof}

\vspace{-1.5em}
\section*{Declaration of interests}
\vspace{-0.5em}
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.
\vspace{-1em}

\printbibliography
\end{document}