\documentclass[a4paper, 10pt]{article}

\usepackage[utf8]{inputenc}
\usepackage{a4wide}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{bm}

\newcommand{\te}{\mathrm{e}}
\newcommand{\ti}{\mathrm{i}}
\newcommand{\td}{\mathrm{d}}

\newcommand{\tn}{\bm{\mathnormal{n}}}
\newcommand{\tr}{\bm{\mathnormal{r}}}
\newcommand{\tx}{\bm{\mathnormal{x}}}
\newcommand{\tT}{^{\mathrm{T}}}
\newcommand{\ty}{\bm{\mathnormal{y}}}
\newcommand{\sa}{\bm{\mathnormal{a}}}
\newcommand{\sx}{\bm{\mathnormal{x}}}
\newcommand{\sy}{\bm{\mathnormal{y}}}
\newcommand{\sn}{\bm{\mathnormal{n}}}
\newcommand{\bN}{\mathbf{N}}
\newcommand{\bp}{\mathbf{p}}
\newcommand{\bq}{\mathbf{q}}
\newcommand{\bH}{\mathbf{H}}
\newcommand{\bG}{\mathbf{G}}
\newcommand{\bI}{\mathbf{I}}
\newcommand{\limeps}{\lim\limits_{\epsilon \rightarrow 0}}
\newcommand{\intsph}{\int\limits_{0}^{2\pi} \int\limits_{0}^{\pi/2}}
\newcommand{\dsph}{\mathrm{d}\theta \mathrm{d}\phi}
\newcommand{\drph}{\mathrm{d} r \mathrm{d}\theta}
\newcommand{\summ}{\sum\limits_{m=1}^3}
\newcommand{\intme}{\int\limits_{\theta_1^m}^{\theta_2^m}\int\limits_\epsilon^{R(\theta)}}
\newcommand{\intmn}{\int\limits_{\theta_1^m}^{\theta_2^m}\int\limits_0^{R(\theta)}}
\newcommand{\intms}{\int\limits_{\theta_1^m}^{\theta_2^m}}


\author{Péter Rucz}
\title{Burton-Miller}

\begin{document}

\maketitle

\section{The Green function and its derivatives}

The Green-function reads as
%
\begin{equation}
	G_k(\tx, \ty) = \frac{\te^{-\ti k r}}{4 \pi r}, \quad
	r = |\ty - \tx|
\end{equation}

\begin{align}
	\frac{\partial G_k(\tx, \ty)}{\partial \tn_y}
	&= -G_k(\tx, \ty) \frac{1 + \ti kr}{r}\frac{\partial r}{\partial \tn_y} \qquad \text{with} \qquad \frac{\partial r}{\partial \tn_y} = \frac{\tr \cdot \tn_y}{r} \\
	\frac{\partial G_k(\tx, \ty)}{\partial \tn_x}
	&= -G_k(\tx, \ty) \frac{1 + \ti kr}{r}\frac{\partial r}{\partial \tn_x} \qquad \text{with} \qquad \frac{\partial r}{\partial \tn_x} = -\frac{\tr \cdot \tn_x}{r}
\end{align}

\begin{align}
	\frac{\partial^2 G_k(\tx, \ty)}{\partial \tn_x \partial \tn_y}
	&= \frac{\te^{- \ti kr}}{4 \pi r^3}\left[\left(3 + 3 \ti kr - k^2 r^2\right) \frac{\partial r}{\partial \tn_x}\frac{\partial r}{\partial \tn_y} + (1 + \ti kr) \tn_x \cdot \tn_y \right] \\
	&= G_k(\tx,\ty)\left[\left(\frac{3}{r^2} + \frac{3\ti k}{r} - k^2\right) \frac{\partial r}{\partial \tn_x}\frac{\partial r}{\partial \tn_y}
	+ \frac{1}{r^2}(1 + \ti kr) \tn_x \cdot \tn_y \right]
\end{align}

The static part is given as:
%
\begin{equation}
	G_0(\tx, \ty) = \frac{1}{4\pi r}
\end{equation}

\begin{equation}
	\frac{\partial G_0(\tx, \ty)}{\partial \tn_x} = -\frac{G_0(\tx, \ty)}{r} \frac{\partial r}{\partial \tn_x}
\end{equation}

\begin{equation}
	\frac{\partial G_0(\tx, \ty)}{\partial \tn_x \partial \tn_y} = \frac{G_0(\tx, \ty)}{r^2} \left[3 \frac{\partial r}{\partial \tn_x}\frac{\partial r}{\partial \tn_y} + \tn_x \cdot \tn_y \right]
\end{equation}

Taking the difference of the second derivatives:
%
\begin{multline}
	\frac{\partial G_k(\tx, \ty)}{\partial \tn_x \partial \tn_y} - \frac{\partial G_0(\tx, \ty)}{\partial \tn_x \partial \tn_y} = \\
	\frac{\te^{- \ti kr}}{4 \pi r^3}\left[\left(3 \cdot (1 - \te^{\ti kr}) + 3 \ti kr - k^2 r^2\right) \frac{\partial r}{\partial \tn_x}\frac{\partial r}{\partial \tn_y} + (1 - \te^{\ti kr} + \ti kr) \tn_x \cdot \tn_y \right]
\end{multline}
%
Taking the limit $r \rightarrow 0$, the difference can be expressed as a Taylor-series
%
\begin{equation}
	\lim\limits_{\sy \rightarrow \sx} \left(\frac{\partial G_k(\tx, \ty)}{\partial \tn_x \partial \tn_y} - \frac{\partial G_0(\tx, \ty)}{\partial \tn_x \partial \tn_y}\right) = \frac{k^2}{8 \pi r} + \frac{\ti k^3}{12 \pi} + o(r)
\end{equation}

\section{The Burton-Miller formulation}

\begin{equation}
	C(\sx) p (\sx) = \int_{\Gamma}\left[\frac{\partial{G}_k (\sx, \sy)}{\partial \sn_y} p(\sy)
	- G_k(\sx,\sy)\frac{\partial p(\sy)}{\partial \sn_y}\right] \td \Gamma_y + p_{in}(\sx)
\end{equation}

The normal derivative ($ \partial / \partial \sn_x$):

\begin{equation}
	C(\sx) \frac{\partial p(\sx)}{\partial \sn_x}
	=
	\int_{\Gamma}\left[\frac{\partial^2{G}_k (\sx, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) - \frac{\partial G_k(\sx,\sy)}{\partial \sn_x}\frac{\partial p(\sy)}{\partial \sn_y}\right] \td \Gamma_y + \frac{\partial p_{in}(\sx)}{\partial \sn_x}
\end{equation}

Using the Burton-Miller formulation for surface points ($C(\sx) = 1/2$), the linear combination is calculated with the coupling constant $\alpha$

\begin{multline}
	C(\sx) p (\sx) - p_{in}(\sx)
	+
	\alpha C(\sx) \frac{\partial p(\sx)}{\partial \sn_x}
	-
	\alpha \frac{\partial p_{in}(\sx)}{\partial \sn_x}
	= \\
	\int_{\Gamma}\left[\frac{\partial{G}_k (\sx, \sy)}{\partial \sn_y} p(\sy)
	- G_k(\sx,\sy)\frac{\partial p(\sy)}{\partial \sn_y}\right] \td \Gamma_y \\
	+
	\alpha \int_{\Gamma}
	\left[
	\frac{\partial^2{G}_k (\sx, \sy)}{\partial \sn_y \partial \sn_x} p(\sy)
	- \frac{\partial G_k(\sx,\sy)}{\partial \sn_x}\frac{\partial p(\sy)}{\partial \sn_y}
	\right] \td \Gamma_y
\end{multline}

Using shape functions $N_i(\sx)$ for the pressure and its normal derivative:

\begin{align}
	p(\sx) & = \sum\limits_{i=1}^P N_i(\sx)p_i \\
	\frac{\partial p(\sx)}{\partial \sn_x} & = q(\sx) = \sum\limits_{i=1}^P N_i(\sx)q_i
\end{align}

Rearranging:

\begin{multline}
	\frac{1}{2} p (\sx)
	-
	p_{in}(\sx)
	+
	\frac{\alpha}{2} q(\sx)
	-
	\alpha q_{in}(\sx)
	= \\
	\left(
	\int_{\Gamma} \frac{\partial{G}_k (\sx, \sy)}{\partial \sn_y} p(\sy) \td \Gamma_y
	+
	\alpha \int_{\Gamma}
	\frac{\partial^2{G}_k (\sx, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) \td \Gamma_y
	\right)
	\\
	-
	\left(
	\int_{\Gamma} G_k(\sx,\sy) q(\sy) \td \Gamma_y
	+
	\alpha \int_{\Gamma}
	\frac{\partial G_k(\sx,\sy)}{\partial \sn_x}q(\sy) \td \Gamma_y
	\right)
\end{multline}

In matrix form:

\begin{align}
	\frac{1}{2} \mathbf{p} - \mathbf{p}_{in} + \frac{\alpha}{2} \mathbf{q} - \alpha \mathbf{q}_{in} = \bH \bp - \bG \bq
\end{align}

With

\begin{align}
	\bH_{ij}
	& = \int_{\Gamma}\frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_y} N_j(\sy) \td \Gamma(\sy)
	+ \alpha \int_\Gamma \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} N_j(\sy) \td \Gamma(\sy) \nonumber \\
	\bG_{ij}
	& = \int_{\Gamma} G_k(\sx_i,\sy) N_j(\sy) \td \Gamma(\sy)
	+ \alpha \int_\Gamma \frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_x} N_j(\sy) \td \Gamma(\sy)
\end{align}




\section{Treatment of hypersingular integrals}

In the sequel, the following four identities are used:
%
% \begin{align}
% 	\alpha & \int_\Gamma \frac{\partial^2{G}_k (\sx, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) \td \Gamma_y -
% 	\alpha \int_\Gamma \frac{\partial{G}_k (\sx, \sy)}{\partial \sn_x} \frac{\partial p(\sy)}{\partial \sn_y} \td \Gamma_y = \nonumber \\
% 	& = \int_{\Gamma}\left[-\frac{\partial{G}_k (\sx, \sy)}{\partial \sn_y} p(\sy) + G_k(\sx,\sy)\frac{\partial p(\sy)}{\partial \sn_y}\right] \td \Gamma_y - \frac{1}{2} p(\sx) - \frac{1}{2}\alpha \frac{ \partial p(\sx)}{\partial \sn_x} + p_{in}(\sx) + \alpha \frac{\partial p_{in}(\sx)}{\partial \sn_x}
% \end{align}
%
\begin{align}
	\int_\Gamma \frac{\partial G_0(\sx, \sy)}{\partial \sn_y} \td \Gamma_y & = 0 \\
	\int_\Gamma \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y & = 0 \\
	\int_\Gamma \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y}(\sy - \sx) \td \Gamma_y & = \int_\Gamma \frac{\partial G_0(\sx, \sy)}{\partial \sn_x} \sn_y \td \Gamma_y \\
	\int_\Gamma \sa \cdot (\sy - \sx) \frac{\partial^2 G_k(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y & = \int_\Gamma \sa \cdot \sn_y \frac{G_k(\sx, \sy)}{\partial \sn_x} \td \Gamma_y \nonumber \\
	& + \int_\Omega - k^2 \sa \cdot (\sy - \sx) \frac{\partial G_k(\sx, \sy)}{\partial \sn_x} \td \Omega_y - \frac{\sa \sn_x}{2}
\end{align}

The hypersingular integral in matrix $\bH$
%
\begin{multline}
	\int_\Gamma \frac{\partial^2 G_k(\sx, \sy)}{\partial \sn_x \partial \sn_y} \phi(\sy) \td \Gamma_y = 
	\int_\Gamma \left[\frac{\partial^2 G_k(\sx, \sy)}{\partial \sn_x \partial \sn_y} - \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y}\right] \phi(\sy) \td \Gamma_y \\
	+ \int_\Gamma \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \phi(\sy) \td \Gamma_y
\end{multline}

The difference of the second derivative of the two Green functions can be expressed by the Taylor-series:
%
\begin{equation}
	\frac{\partial^2 G_k(\sx, \sy)}{\partial \sn_x \partial \sn_y} - \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} = 
	\frac{k^2}{8 \pi r} + \frac{\ti k^3}{12 \pi} + o(r)
\end{equation}

This shows that the first integral on the right hand side is only weakly singular.
The second integral needs to be taken care of. By adding and subtracting the same terms:
%
\begin{multline}
	\int_\Gamma \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \phi(\sy) \td \Gamma_y =
	\int_\Gamma \left[\phi(\sy) - \phi(\sx) - \nabla \phi(\sx) (\sy - \sx) \right] \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y \\
	+ \phi(\sx) \int_\Gamma \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y +
	\int_\Gamma \nabla \phi(\sx) (\sy - \sx) \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y
\end{multline}

Making use of identity 2, the second term on the right hand side vanishes, which leads to
%
\begin{multline}
	\int_\Gamma \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \phi(\sy) \td \Gamma_y =
	\int_\Gamma \left[\phi(\sy) - \phi(\sx) - \nabla \phi(\sx) (\sy - \sx) \right] \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y \\
	+ \int_\Gamma \nabla \phi(\sx) (\sy - \sx) \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y
\end{multline}

The terms in brackets in the first integral of the right hand side are the difference of $\phi(\sy)$ and its Taylor series, which means that the difference in the brackets is $o(r^2)$ here and therefore the first integral is now weakly singular.

The second term on the RHS should be treated by making use of identity 4, with $\sa = \nabla\phi(\sx)$ and $k = 0$.

\begin{align}
	\int_\Gamma \nabla \phi(\sx) (\sy - \sx) \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y & = 
	\int_\Gamma \nabla \phi(\sx) \cdot \sn_y \frac{\partial G_0(\sx,\sy)}{\partial \sn_x} \td \Gamma_y \nonumber \\
	&+ \int_\Omega -k^2 \nabla \phi(\sx) \cdot (\sy - \sx) \frac{\partial G_0(\sx, \sy)}{\partial \sn_x} \td \Omega_y
	-\frac{1}{2} \nabla \phi(\sx) \cdot \sn_x
\end{align}

Since $k = 0$ in the static case, the volume integral on the RHS disappears and we get
%
\begin{equation}
	\int_\Gamma \nabla \phi(\sx) (\sy - \sx) \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y = 
	\int_\Gamma \nabla \phi(\sx) \cdot \sn_y \frac{\partial G_0(\sx,\sy)}{\partial \sn_x} \td \Gamma_y -\frac{1}{2} \nabla \phi(\sx) \cdot \sn_x
\end{equation}

The surface integral of the static green function is obtained
%
\begin{multline}
	\int_\Gamma \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \phi(\sy) \td \Gamma_y
= \int_\Gamma \left[\phi(\sy) - \phi(\sx) - \nabla \phi(\sx) (\sy - \sx) \right] \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y \\
	+ \int_\Gamma \nabla \phi(\sx) \cdot \sn_y \frac{\partial G_0(\sx,\sy)}{\partial \sn_x} \td \Gamma_y -\frac{1}{2} \nabla \phi(\sx) \cdot \sn_x
\end{multline}

And finally, the hypersingular integral can be written as
%
\begin{multline}
	\int_\Gamma \frac{\partial^2 G_k(\sx, \sy)}{\partial \sn_x \partial \sn_y} \phi(\sy) \td \Gamma_y = \int_\Gamma \left[\frac{\partial^2 G_k(\sx, \sy)}{\partial \sn_x \partial \sn_y} - \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y}\right] \phi(\sy) \td \Gamma_y \\
	+ \int_\Gamma \left[\phi(\sy) - \phi(\sx) - \nabla \phi(\sx) (\sy - \sx) \right] \frac{\partial^2 G_0(\sx, \sy)}{\partial \sn_x \partial \sn_y} \td \Gamma_y \\
	+ \int_\Gamma \nabla \phi(\sx) \cdot \sn_y \frac{\partial G_0(\sx,\sy)}{\partial \sn_x} \td \Gamma_y -\frac{1}{2} \nabla \phi(\sx) \cdot \sn_x
\end{multline}

\section{Constant triangular elements}

The boundary integrals in the Burton-Miller formulation:

\begin{multline}
	\frac{1}{2} p (\sx)
	-
	p_{in}(\sx)
	+
	\frac{\alpha}{2} q(\sx)
	-
	\alpha q_{in}(\sx)
	= \\
	\left(
	\int_{\Gamma} \frac{\partial{G}_k (\sx, \sy)}{\partial \sn_y} p(\sy) \td \Gamma_y
	+
	\alpha \int_{\Gamma}
	\frac{\partial^2{G}_k (\sx, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) \td \Gamma_y
	\right)
	\\
	-
	\left(
	\int_{\Gamma} G_k(\sx,\sy) q(\sy) \td \Gamma_y
	+
	\alpha \int_{\Gamma}
	\frac{\partial G_k(\sx,\sy)}{\partial \sn_x}q(\sy) \td \Gamma_y
	\right)
\end{multline}

The integrals 
%

\begin{align}
	H = \int_{\Gamma}\frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_y} p(\sy) \td \Gamma(\sy)
	+ \alpha \int_\Gamma \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) \td \Gamma(\sy) \nonumber \\
	G = \int_{\Gamma} G_k(\sx_i,\sy) \frac{\partial p(\sy)}{\partial \sn_y} \td \Gamma(\sy)
	+ \alpha \int_\Gamma \frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_x} \frac{\partial p(\sy)}{\partial \sn_y}  \td \Gamma(\sy)
\end{align}

become hypersingular, which is treated by extending the integration area by a hemishphere and taking the limit.
By taking a hemisphere of radius $\epsilon$ in the element center the integration domain for a~triangular element can be expanded as
%
\begin{align}
	\int_{\Gamma_\Delta} f(\sy)\td\Gamma(\sy) = \lim\limits_{\epsilon \rightarrow 0 }\left[\int_{\Gamma_\epsilon} f(\sy) \td \Gamma(\sy) + \int_{S_{\epsilon}} f(\sy) \td \Gamma (\sy) \right] .
\end{align}
%

% \begin{align}
% 	\bH_{ii}
% 	& = \int_{\Gamma}\frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_y} N_i(\sy) \td \Gamma(\sy)
% 	+ \alpha \int_\Gamma \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} N_i(\sy) \td \Gamma(\sy) \nonumber \\
% 	\bG_{ii}
% 	& = \int_{\Gamma} G_k(\sx_i,\sy) N_i(\sy) \td \Gamma(\sy)
% 	+ \alpha \int_\Gamma \frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_x} N_i(\sy) \td \Gamma(\sy)
% \end{align}
%
\begin{multline}
	H_\Delta
	=
	\underbrace{\limeps \int_{\Gamma_\epsilon}\frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_y} p(\sy) \td \Gamma(\sy)}_{H_1}
    +
	\underbrace{\limeps \int_{S_\epsilon}\frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_y} p(\sy) \td \Gamma(\sy) \nonumber}_{H_2} \\
	+
	\alpha \underbrace{\limeps \int_{\Gamma_\epsilon} \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) \td \Gamma(\sy)}_{H_3} 
    +
	\alpha \underbrace{\limeps \int_{\S_\epsilon} \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) \td \Gamma(\sy)}_{H_4}
\end{multline}
%
and
%
\begin{multline}
	G_\Delta
	=
	\underbrace{\limeps \int_{\Gamma_\epsilon} G_k(\sx_i,\sy) \frac{\partial p(\sy)}{\partial \sn_y}\td \Gamma(\sy)}_{G_1}
	+
	\underbrace{\limeps \int_{S_\epsilon} G_k(\sx_i, \sy) \frac{\partial p(\sy)}{\partial \sn_y} \td \Gamma(\sy) \nonumber}_{G_2} \\
	+
	\alpha \underbrace{\limeps \int_{\Gamma_\epsilon} \frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_x} \frac{\partial p(\sy)}{\partial \sn_y} \td \Gamma(\sy)}_{G_3}
	+
	\alpha \underbrace{\limeps \int_{S_\epsilon} \frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_x} \frac{\partial p(\sy)}{\partial \sn_y} \td \Gamma(\sy)}_{G_4}
\end{multline}

Taking the first integral of matrix $\bG$ on the triangle and making use of the identity:
%
\begin{align}
	\int \frac{\te ^{-\ti kr}}{4 \pi} \mathrm{d} r = \frac{\ti}{ k} \frac{\te ^{-\ti kr}}{4 \pi}
\end{align}

The first integral, $G_1$, making use of the identity on the constant triangle $\partial p(\sy) / \partial \sn_y = q(\sx)$ reads as
%
\begin{align}
	G_1 = \limeps \int_{\Gamma_\epsilon} G_k(\sx_i,\sy) \frac{\partial p(\sy)}{\partial \sn_y} \td \Gamma(\sy) 
	& = \limeps \int_{\Gamma_\epsilon} G_k(\sx_i,\sy) \td \Gamma(\sy) q(\sx) \nonumber \\
	& = \limeps \summ \intme \frac{\te ^{- \ti k r}}{4 \pi r} r \drph q(\sx) \nonumber \\
	& = \summ \intmn \frac{\te^{-\ti kr}}{4 \pi} \drph q(\sx) \nonumber \\
	& = \summ \intms \frac{\ti}{4\pi k}\left(\te^{-\ti k R(\theta)} -1\right) \mathrm{d} \theta q(\sx) \nonumber \\
    & = -\frac{\ti}{2k} q(\sx) + \frac{\ti}{4\pi k}\summ \intms \te ^{-\ti kR(\theta)} \mathrm{d} \theta q(\sx)
\end{align}

For the first term of matrix $\bH$ first integral over the triangle, it can be used that $\partial r / \partial \sn_y$ over the planar element, and therefore the term vanishes
%
\begin{align}
	H_1 = \limeps \int_{\Gamma_\epsilon}\frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_y} p(\sy) \td \Gamma(\sy) =
	\limeps \int_{\Gamma_\epsilon}-\frac{\te ^{-\ti k r}}{4 \pi r^2}(1 + \ti kr) \frac{\partial r}{\partial \sn_y} p(\sy) \mathrm{d} \Gamma_y = 0
\end{align}

Introducing the coordinate system $\tau(\sx)$, $s(\sx)$ and $n(\sx)$ over the Hemisphere, the second integral in $G_\Delta$ can be expressed as:

\begin{align}
	G_2 = \int_{S_\epsilon} G_k(\sx_i, \sy) \frac{\partial p}{\partial \sn_y} \td \Gamma(\sy) & = \intsph \frac{\te^{-\ti k \epsilon}}{4 \pi \epsilon} \frac{\partial p}{\partial \sn_y} \epsilon^2 \sin \theta \mathrm{d} \theta \mathrm{d} \phi \nonumber \\
	& = \frac{\epsilon \te ^{-\ti k \epsilon}}{4 \pi} \intsph \left(\frac{\partial p(\sy)}{\partial \sn_y} -\frac{\partial p(\sx)}{\partial \sn_y} \right) \sin \theta \mathrm{d} \theta \mathrm{d} \phi \nonumber \\
	& + \frac{\epsilon \te ^{-\ti k \epsilon}}{4 \pi} \intsph \sn_y \sin \theta \mathrm{d} \theta \mathrm{d} \phi \cdot \nabla p(\sx)
\end{align}

Taking the limit $\epsilon \rightarrow 0$

\begin{align}
	\limeps G_2 = 0
\end{align}

Term $H_2$ can be expressed by making use of the fact that $\partial r / \partial n = 1$ on the hemishpere surface:

\begin{align}
	\int_{S_\epsilon}\frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_y} p(\sy) \td \Gamma(\sy) & = 
	\intsph  \left[-\frac{\te^{-\ti k \epsilon}}{4 \pi \epsilon^2}(1 + \ti k \epsilon) p(\sy) \right] \epsilon^2 \sin \theta \mathrm{d} \theta \mathrm{d} \phi = \nonumber \\
	& = -\frac{\te ^{- \ti k\epsilon}}{4\pi}(1 + \ti k\epsilon) \intsph \left(p(\sy) - p(\sx)\right) \sin \theta \dsph \nonumber \\
	& + -\frac{\te ^{ -\ti k\epsilon}}{4\pi}(1 + \ti k\epsilon) \intsph \sin \theta \dsph p(\sx)
\end{align}

The first integral on the RHS vanishes as $\epsilon \rightarrow 0$ and the second integral is constant (equals $2\pi$), therefore

\begin{equation}
	\limeps H_2 = -\frac{1}{2}p(\sx)
\end{equation}

Now, looking at the components arising from the Burton-Miller formulation:

The third part in matrix $\bG$ vanishes since $\partial r / \partial \sn_x = 0$ 

\begin{align}
	G_3 = \limeps \int_{\Gamma_\epsilon} \frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_x} \frac{\partial p(\sy)}{\partial \sn_y}\td \Gamma(\sy) = 
	\int_{\Gamma_\epsilon} \frac{\te ^{-\ti kr}}{4 \pi r^2}(1 + \ti kr) \frac{\partial r}{\partial \sn_x} \frac{\partial p(\sy)}{\partial \sn_y} \mathrm{d} \Gamma_y = 0
\end{align}

The fourth part in matrix $\bG$:

\begin{align}
	G_4 = \int_{S_\epsilon} \frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_x} \frac{\partial p}{\partial \sn_y}\td \Gamma(\sy) & = 
	\intsph - \frac{\te ^{-\ti k\epsilon}}{4 \pi \epsilon^2} (1 + \ti k \epsilon) \frac{\partial r}{\partial \sn_x} \frac{\partial p}{\partial \sn_y} \epsilon^2 \sin \theta \dsph \nonumber \\
	& = \intsph- \frac{\te ^{-\ti k\epsilon}}{4 \pi \epsilon^2} (1 + \ti k \epsilon) \frac{\partial r}{\partial \sn_x}\left[ \frac{\partial p(\sy)}{\partial \sn_y} - \frac{\partial p(\sx)}{\partial \sn_y}\right] \epsilon^2 \sin \theta \dsph \nonumber \\
	& + \intsph - \frac{\te ^{-\ti k\epsilon}}{4 \pi \epsilon^2} (1 + \ti k \epsilon) \frac{\partial r}{\partial \sn_x} \frac{\partial p(\sx)}{\partial \sn_y}  \epsilon^2 \sin \theta \dsph 
\end{align}

The first integral on the RHS vanishes and by using the relation $\partial r / \partial \sn_x = - \cos \theta$ on the hemisphere, then we have:

\begin{align}
	G_4 & = \frac{\te ^{-\ti k\epsilon}}{4 \pi} (1 + \ti k \epsilon) \intsph \frac{\partial p(\sx)}{\partial \sn_y} \sin \theta \cos \theta \dsph \nonumber \\
	& =  \frac{\te ^{-\ti k\epsilon}}{4 \pi} (1 + \ti k \epsilon) \intsph \left[ \sin \theta \cos \theta \frac{\partial p(\sx)}{\partial \tau} + \sin \theta \sin \phi \frac{\partial p(\sx)}{\partial s} + \cos\theta \frac{\partial p(\sx)}{\partial n} \right] \sin\theta \cos\theta \dsph
\end{align}

The results of the three integral terms on the RHS are $0$, $0$ and $2\pi / 3$, respectively, which gives:

\begin{align}
	G_4 & = \frac{1}{4\pi}\frac{2\pi}{3}\frac{\partial p(\sx)}{\partial n} = \frac{1}{6} \frac{\partial p(\sx)}{\partial n} = \frac{1}{6} q(\sx)
\end{align}

The third term in matrix $\bH$ can be integrated by making use of $\sn(\sx) = \sn(\sy)$ and $\partial r / \partial \sn_y = 0$ on $\Gamma_\epsilon$ the double derivative of the Green function can be simplified as:

\begin{align}
	\frac{\partial^2 G_k(\tx, \ty)}{\partial \tn_x \partial \tn_y}
	&= \frac{\te^{- \ti kr}}{4 \pi r^3}\left[\left(3 + 3 \ti kr - k^2 r^2\right) \frac{\partial r}{\partial \tn_x}\frac{\partial r}{\partial \tn_y} + (1 + \ti kr) \tn_x \cdot \tn_y \right] = \nonumber\\
	& = \frac{\te^{- \ti kr}}{4 \pi r^3}\left(1 + \ti kr\right),  \qquad \text{if} \quad \sx, \sy \in \Gamma_\epsilon
\end{align}

Therefore

\begin{align}
	H_3 = \limeps \int_{\Gamma_\epsilon} \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} \td \Gamma(\sy) & =
	\limeps \summ \intme \frac{\te^{- \ti kr}}{4 \pi r^3}\left(1 + \ti kr\right) r \mathrm{d} r \mathrm{d} \theta \nonumber \\
\end{align}

Since the term appearing in the integral on the RHS is the negative derivative of the Green's function

\begin{align}
	\frac{\te^{- \ti kr}}{4 \pi r^3}\left(1 + \ti kr\right) r = -\frac{\partial}{\partial r}\frac{\te ^ {-\ti kr}}{4\pi r}
\end{align}

The primitive function can be substituted into the integral

\begin{align}
	H_3 = \limeps \summ \intme \frac{\te^{- \ti kr}}{4 \pi r^3}\left(1 + \ti kr\right) r \mathrm{d} r \mathrm{d} \theta & =
	\limeps \summ \intms \left[\frac{\te ^{-\ti k \epsilon}}{4 \pi \epsilon} - \frac{\te ^{-\ti k R(\theta)}}{4 \pi R(\theta)}\right] \mathrm{d} \theta \nonumber\\ 
	& = \frac{1}{2\epsilon} - \frac{\ti k}{2}- \summ \intms \frac{\te ^{-\ti k R(\theta)}}{4 \pi R(\theta)} \mathrm{d} \theta
\end{align}

(Note: the term $-\frac{\ti k}{2}$ is from the Taylor-series approximation.)

Finally, term $H_4$

\begin{align}
	H_4 = \int_{S_\epsilon} \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) \td \Gamma(\sy)
	& = \int_{S_\epsilon} \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} \left[p(\sy) - p(\sx) - \nabla p(\sx) (\sy - \sx)\right]\nonumber \\
	& + p(\sx) \int_{S_\epsilon} \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} \mathrm{d} \Gamma_y
	+ \nabla p(\sx) \cdot \int_{S_\epsilon} \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} (\sy - \sx) \Gamma_y
\end{align}

On the hemisphere surface $\partial r / \partial n(\sx) = -\cos(\theta)$, $\partial r / \partial n(\sy) = 1$ and $r = \epsilon$ and  $\sn_x \cdot \sn_y = \cos \theta$:

\begin{align}
	\frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} & = 
	\frac{\te^{- \ti kr}}{4 \pi r^3}\left[\left(3 + 3 \ti kr - k^2 r^2\right) \frac{\partial r}{\partial \tn_x}\frac{\partial r}{\partial \tn_y} + (1 + \ti kr) \tn_x \cdot \tn_y \right] \nonumber \\
	& = \frac{\te^{- \ti k\epsilon}}{4 \pi \epsilon^3}\left[\left(3 + 3 \ti k\epsilon - k^2 r^2\right) (-\cos \theta) + (1 + \ti k\epsilon) \cos \theta) \right] \nonumber \\
	& = - \frac{\te^{-\ti k \epsilon}}{4 \pi \epsilon^3} \left[2 + 2\ti k \epsilon -k^2 \epsilon^2\right] \cos \theta , \qquad \text{if} \quad \sx, \sy \in S_\epsilon
\end{align}

The first integral on the RHS for $H_4$ gives zero ($H_{41} = 0$) in the limit $\epsilon \rightarrow 0$, because of Hölder 1,$\alpha$ continuity of $p(\sy)$. 

The second term $H_{42}$ can be treated as:

\begin{align}
	H_{42} &  = \limeps \intsph \frac{\te ^{-\ti k \epsilon}}{4 \pi \epsilon^3}\left[-2 -2\ti k \epsilon + k^2 \epsilon^2\right] \cos \theta \epsilon^2 \sin \theta \dsph \nonumber \\
	& = \limeps \frac{\te ^{-\ti k \epsilon}}{4 \pi \epsilon}\left[-2 -2\ti k \epsilon + k^2 \epsilon^2\right] \intsph \sin \theta \cos \theta \dsph \nonumber \\
	& = \limeps \frac{\te ^{-\ti k \epsilon}}{4 \pi \epsilon}\left[-2 -2\ti k \epsilon + k^2 \epsilon^2\right] \pi = - \frac{1}{2 \epsilon} + \frac{\ti k}{2} - \frac{\ti k}{2} = -\frac{1}{2 \epsilon} 
\end{align}

Now, the following is applied:

\begin{align}
	(\sy - \sx) \cdot \nabla p(\sx) & = \left[(\sin \theta \cos \phi \bm{\tau} + (\sin\theta \sin\phi)\bm{s} + \cos \theta \bm{n}
	\right] \frac{\partial p(\sx)}{\partial n} \bm{n} = r \frac{\partial p(\sx)}{\partial n} \cos \theta \nonumber \\
	& = \frac{\partial p(\sx)}{\partial n} \cos^2 \theta \qquad \text{if} \quad \sx, \sy \in S_\epsilon
\end{align}

The third term $H_{43}$ can be evaluated as:

\begin{align}
	H_{43} & = \limeps\nabla p(\sx) \cdot \int_{S_\epsilon} \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} (\sy - \sx) \Gamma_y \nonumber \\
	&= \limeps \intsph \frac{\te^{-\ti k \epsilon}}{4 \pi \epsilon^3}\left[-2 -2\ti k \epsilon + k^2 \epsilon^2\right] \epsilon \frac{\partial p(\sx)}{\partial \sn_x} \cos^2 \theta \epsilon^2 \sin \theta \dsph \nonumber \\
	&= \limeps \frac{\te^{-\ti k \epsilon}}{4 \pi}\left[-2 -2\ti k \epsilon + k^2 \epsilon^2\right] \frac{\partial p(\sx)}{\partial n} \intsph \sin \theta \cos^2 \theta \dsph \nonumber \\
	& = \limeps \frac{\te^{-\ti k \epsilon}}{4 \pi}\left[-2 -2\ti k \epsilon + k^2 \epsilon^2\right] \frac{\partial p(\sx)}{\partial n} \frac{2 \pi}{3} = -\frac{1}{3} \frac{\partial{p(\sx)}}{\sn_x}
\end{align}

Summing up the terms:

\begin{align}
	H & = H_1 + H_2 + \alpha H_3 + \alpha \left(H_{41} + H_{42} + H_{43}\right) = \nonumber \\ 
	&= 0 - \frac{1}{2} p(\sx) + \alpha \left[\frac{1}{2\epsilon} - \frac{\ti k}{2}- \summ \intms \frac{\te ^{-\ti k R(\theta)}}{4 \pi R(\theta)} \mathrm{d} \theta\right] p(\sx) + \alpha \left[0 - \frac{1}{2\epsilon} p(\sx) - \frac{1}{3} \frac{\partial{p(\sx)}}{\partial n}\right]
\end{align}

\begin{align}
	G & = G_1 + G_2 + \alpha G_3 + \alpha G_4 \nonumber \\
	& = \frac{i}{2k}\left[\frac{1}{2\pi}\summ \intms \te ^{-\ti kR(\theta)} \mathrm{d} \theta -1 \right]\frac{\partial p(\sx)}{\partial n} + 0 + \alpha \cdot 0 + \alpha \frac{1}{6} \frac{\partial p(\sx)}{\partial n}
\end{align}

Substituting into the Burton-Miller formulation for the surface matrices:

\begin{multline}
	\frac{1}{2} p (\sx)
	-
	p_{in}(\sx)
	+
	\frac{\alpha}{2} q(\sx)
	-
	\alpha q_{in}(\sx)
	= \\
	\left(
	\int_{\Gamma} \frac{\partial{G}_k (\sx, \sy)}{\partial \sn_y} p(\sy) \td \Gamma_y
	+
	\alpha \int_{\Gamma}
	\frac{\partial^2{G}_k (\sx, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) \td \Gamma_y
	\right)
	\\
	-
	\left(
	\int_{\Gamma} G_k(\sx,\sy) q(\sy) \td \Gamma_y
	+
	\alpha \int_{\Gamma}
	\frac{\partial G_k(\sx,\sy)}{\partial \sn_x}q(\sy) \td \Gamma_y
	\right)
\end{multline}

\begin{multline}
	\frac{1}{2} p (\sx)
	-
	p_{in}(\sx)
	+
	\frac{\alpha}{2} q(\sx)
	-
	\alpha q_{in}(\sx)
	= \\
	\int_{\Gamma-\Gamma_\Delta} \frac{\partial{G}_k (\sx, \sy)}{\partial \sn_y} p(\sy) \td \Gamma_y
	+
	\left[-\frac{1}{2}p(\sx) - \alpha \summ \intms \frac{\te ^{-\ti k R(\theta)}}{4 \pi R(\theta)} \mathrm{d} \theta p(\sx)\right] 
	\\ + 
	\alpha \left[\int_{\Gamma-\Gamma_\Delta}
	\frac{\partial^2{G}_k (\sx, \sy)}{\partial \sn_y \partial \sn_x} p(\sy) \td \Gamma_y
	- \frac{\ti k}{2} p(\sx) - \frac{1}{3}q(\sx)\right]
	\\
	-
	\int_{\Gamma-\Gamma_\Delta} G_k(\sx,\sy) q(\sy) \td \Gamma_y
	-
	\frac{\ti}{2k}\left[\frac{1}{2\pi}\summ \intms \te ^{-\ti kR(\theta)} \mathrm{d} \theta -1 \right] q(\sx)
	\\
	-
	\alpha \int_{\Gamma-\Gamma_\Delta}
	\frac{\partial G_k(\sx,\sy)}{\partial \sn_x}q(\sy) \td \Gamma_y
	- \alpha \frac{1}{6}q(\sx)
\end{multline}

Introducing matrix notations, making use of constant shape functions:

\begin{align}
	\bH_{ij}
	& = \int_{\Gamma_j}\frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_y} \td \Gamma_j(\sy)
	+ \alpha \int_\Gamma \frac{\partial^2{G}_k (\sx_i, \sy)}{\partial \sn_y \partial \sn_x} \td \Gamma_j(\sy)
	& \text{if} \quad i \ne j \nonumber \\
	\bH_{ii}
	& = - \alpha \summ \intms \frac{\te ^{-\ti k R(\theta)}}{4 \pi R(\theta)} \mathrm{d} \theta 
	& \text{if} \quad i = j \nonumber \\
	\bG_{ij}
	& = \int_{\Gamma_j} G_k(\sx_i,\sy) \td \Gamma_j(\sy)
	+ \alpha \int_{\Gamma_j} \frac{\partial{G}_k (\sx_i, \sy)}{\partial \sn_x} \td \Gamma_j(\sy)
	& \text{if} \quad i \ne j \nonumber \\
	\bG_{ii}
	& = \frac{\ti}{4\pi k}\summ \intms \te ^{-\ti kR(\theta)} \mathrm{d} \theta 
	& \text{if} \quad i = j
\end{align}
%
The matrix formulation reads as:
%
\begin{align}
	\frac{1}{2} \mathbf{p} - \mathbf{p}_{in} + \frac{\alpha}{2} \mathbf{q} - \alpha \mathbf{q}_{in} = - \alpha \frac{\ti k}{2} \bp 
	+ \frac{\ti }{2k} \bq + \bH \bp - \bG \bq
\end{align}
%
Rearranging leads to:
%
\begin{align}
	\left(\frac{1}{2} + \frac{\alpha \ti k}{2}\right)\bp - \bp_{in}+ \left(\alpha - \frac{\ti }{2 k}\right) \bq - \alpha \bq_{in} =\bH \bp - \bG \bq
\end{align}
%
Solution without incident field is obtained either for given $\bp$ or $\bq$ on the surface as:
%
\begin{align}
	\left[\bH - \bI\left(\frac{1}{2} + \frac{\alpha \ti k}{2}\right) \right]^{-1} \bp = \left[\bG + \bI \left(\alpha - \frac{\ti}{2k}\right) \right] \bq .
\end{align}



\end{document}
