# Closed-form evaluation of potential integrals in the Boundary Element Method
111This work was partially carried out under the Horizon 2020 project
AERIALIST, project identifier 723367

Michael Carley<EMAIL_ADDRESS>Department of Mechanical Engineering,
University of Bath, Bath, BA2 7AY, United Kingdom

###### Abstract

A method is presented for the analytical evaluation of the singular and near-
singular integrals arising in the Boundary Element Method solution of the
Helmholtz equation. An error analysis is presented for the numerical
evaluation of such integrals on a plane element, and used to develop a
criterion for the selection of quadrature rules. The analytical approach is
based on an optimized expansion of the Green’s function for the problem,
selected to limit the error to some required tolerance. Results are presented
showing accuracy to tolerances comparable to machine precision.

##  1 Introduction

A central part of the Boundary Element Method (BEM) is the evaluation of
potential integrals, to compute the contribution of an element to the
potential field, or to the entries of the solution matrix. It is thus a key
factor in the accuracy and efficiency of any implementation, and one which has
attracted great interest over many decades. In this paper we develop a method
for the evaluation of integrals which arise in the three-dimensional BEM for
the wave equation, in particular in acoustics, where the acoustic potential
$\phi$, external to a surface $A$, is given by the integral formulation:

$\displaystyle\phi(\mathbf{x})$
$\displaystyle=\int_{A}\frac{\partial\phi_{1}}{\partial
n}G(\mathbf{x},\mathbf{x}_{1})-\frac{\partial
G(\mathbf{x},\mathbf{x}_{1})}{\partial n}\phi_{1}\,\mathrm{d}A,$ (1)

where $\mathbf{x}$ indicates position, subscript $1$ variables of integration
on the surface $A$, and $n$ the outward pointing normal to the surface. The
Green’s function $G$ is:

$\displaystyle G(\mathbf{x};\,\mathbf{x}_{1})$
$\displaystyle=\frac{\mathrm{e}^{\mathrm{j}kR}}{4\pi R},$ (2) $\displaystyle
R$ $\displaystyle=|\mathbf{x}-\mathbf{x}_{1}|,$

where $k$ is acoustic wavenumber. Given the surface potential $\phi$ and
gradient $\partial\phi/\partial n$, the potential, and, after differentiation,
its gradient(s), can be evaluated at any point in the field. Also, given a
boundary condition for $\phi$ and/or $\partial\phi/\partial n$ on $A$, the
integral equation can be solved for $\phi(\mathbf{x})$ and/or
$\partial\phi/\partial n(\mathbf{x})$, $\mathbf{x}\in A$.

If the boundary integral equation is solved using a collocation method, the
surface $A$ is divided into elements, here taken to be plane triangles, and
suitable shape functions are used to interpolate the potential on these
elements. The integral equation is transformed to a linear system in the
element potentials, with the influence coefficients determined by the
potential generated by each element at each node of the surface mesh. This
leads to the requirement to evaluate integrals $I$ and $\partial I/\partial n$
where:

$\displaystyle I$
$\displaystyle=\iint_{A_{e}}f(\xi,\eta)G(\mathbf{x},\mathbf{x}_{1}(\xi,\eta))\,\mathrm{d}A_{e},$
(3)

with $A_{e}$ the surface of an element and $(\xi,\eta)$ a coordinate system
local to $A_{e}$. The requirement then is to evaluate integrals of
$\exp[\mathrm{j}kR]/R$ and its derivatives over a triangular element. This is
especially challenging when the field point $\mathbf{x}$ is on, or near, the
element, and the $1/R$ singularity must be accommodated in the integration
scheme.

There are numerous numerical schemes for the evaluation of the surface
integrals, which mainly vary in their approach to dealing with the
singularity. There are also a number of analytical schemes for the equivalent
integral in the Laplace equation [4, 8, 7, 3, 10, 12, for example], some of
which can be used to deal with the singular terms in the acoustic problem and
thus ease numerical integration, but there are few analytical methods for the
Helmholtz problem. Clearly, given the absence of an exact analytical solution
for the retarded potential from a plane element, any closed-form solution is
an approximation, but it should be possible to approximate the integral to any
required accuracy, in a form amenable to analytical manipulation, so that the
result can be used as if it were an analytical formula for the potential. This
is especially important for the case of a field point on or near the element,
where the ability to handle singularities analytically offers an advantage
over purely numerical schemes.

To the author’s knowledge there are two published methods for closed-form or
analytical evaluation of the Helmholtz potential from a planar element [13,
9]. These use two different approaches to the problem. In one [13], an
expression is derived in the Fourier domain resulting in an expression based
on a series of terms defined by integrals of Hankel functions. These integrals
can be evaluated analytically in terms of Struve functions, yielding a closed-
form solution for the potential from a planar element, but at the expense of
using special functions not routinely available in numerical libraries.

The second approach [9], which is similar in spirit to the method of this
paper, makes use of results derived for the Laplace problem [7] and
approximates $\exp[\mathrm{j}kR]$ as a polynomial over the element. This is
justified by noting that the element size is limited by the requirement to
avoid aliasing in the representation of the surface potential, so that a
relatively low-order approximation containing five or six terms of the Taylor
series for $\exp[\mathrm{j}kR]$ is adequate for evaluation of the integrals to
the tolerance specified.

The method of this paper uses a similar approach, in that it replaces the
exponential with an approximation of controlled error, and uses results from
an analysis of the Laplace problem [4] to compute the terms in the resulting
expansion. It differs in the form of Laplace solution used, and in the choice
of expansion for the exponential, to give a systematic control on quadrature
accuracy optimized to require a minimum number of terms. Additionally, a
criterion is provided for choosing when, and when not, to use the analytical
approach or a purely numerical method, based on an error analysis of
integration using a polar coordinate transformation. To the author’s
knowledge, this error analysis is novel and may have applications more
generally.

##  2 Integration of the $1/R$ potential

In order to motivate the development of the closed-form expression for the
acoustic potential, we begin by analyzing the numerical evaluation of the
Laplace potential, which corresponds to the leading-order, singular, part of
the Helmholtz potential, which gives rise to the difficulties in numerical
integration.

Figure 1: Integration of Laplace potential over a triangle

The model problem is shown in Figure 1 and consists of the evaluation of

$\displaystyle I$
$\displaystyle=\int_{A}\frac{r}{R}\,\mathrm{d}r\mathrm{d}\theta,$
$\displaystyle R^{2}$ $\displaystyle=r^{2}+z^{2},$

over the area of the triangle shown, which lies in the plane $z=0$, with the
usual transformation to polar coordinates $(r,\theta)$ for the integration.

The error in the evaluation of this integral, especially at small values of
$z$ arises from the singular, or near-singular, term $1/R$. Here we develop an
approximate error analysis for the evaluation of this term, which can be used
in determining the required order of integration for $r/R$ or when to switch
to some other quadrature approach, such as that in the next section. An error
analysis for integration using the polar coordinate transformation has been
published previously [11] but the analysis presented here appears to be novel
and is simple enough for use as an a priori estimator in determining
quadrature order in applications.

The analysis depends on an error estimate for the $1/R$ term in a numerical
polar integration, such as (8a). If such an integration is performed using
Gaussian quadrature, the integrand is being approximated by a polynomial over
the interval of integration and the accuracy of the approximation is
determined by the number of terms required to approximate the integrand
accurately. We perform the analysis by estimating the error in the polynomial
expansion of $1/R$ and use it to give an approximation of the order of
polynomial required to approximate $1/R$ to a given tolerance. Given this
polynomial order, a Gaussian quadrature of sufficiently high order can be
selected, or if the order required is too great, the analytical method of the
following section can be used.

Figure 2: Notation for error analysis of $1/R$ expansion

From Figure 2, we write

$\displaystyle r$ $\displaystyle=r_{\text{mid}}-tr_{\text{mid}},\quad
r_{\text{mid}}=r_{\text{max}}/2,\quad-1\leq t\leq 1,$ $\displaystyle R^{2}$
$\displaystyle=r^{2}+z^{2}=R_{\text{mid}}^{2}\left[1-2\cos^{2}\phi
t+(t\cos\phi)^{2}\right],$ $\displaystyle R_{\text{mid}}^{2}$
$\displaystyle=r_{\text{mid}}^{2}+z^{2},\quad\cos\phi=r_{\text{mid}}/R_{\text{mid}},$

and neglect the case of $\phi=0$ as in this case $r/R\equiv 1$ and the
polynomial representation of the integrand raises no difficulties.

Using the generating function for Legendre polynomials [5, 8.921],

$\displaystyle\frac{1}{R}$
$\displaystyle=\frac{1}{R_{\text{mid}}}\sum_{q=0}^{\infty}(t\cos\phi)^{q}P_{q}(\cos\phi).$
(4)

If the expansion is truncated at $q=Q$, the error at any value of $t$ is given
by the remainder

$\displaystyle\epsilon_{Q}$
$\displaystyle=\frac{1}{R}-\frac{1}{R_{\text{mid}}}\sum_{q=0}^{Q}(t\cos\phi)^{q}P_{q}(\cos\phi)=\frac{1}{R_{\text{mid}}}\sum_{q=Q+1}^{\infty}(t\cos\phi)^{q}P_{q}(\cos\phi),$

which can be rewritten using the large-order asymptotic form of the Legendre
polynomial [5, 8.918],

$\displaystyle P_{q}(\cos\phi)\sim\left(\frac{2}{\pi
q\sin\phi}\right)^{1/2}\cos\left[(q+1/2)\phi-\pi/4\right],$

so that

$\displaystyle\epsilon_{Q}$
$\displaystyle\approx\frac{1}{R_{\text{mid}}}\frac{1}{(\pi\sin\phi)^{1/2}}\sum_{q=Q+1}^{\infty}t^{q}\frac{\cos^{q}\phi}{q^{1/2}}\left[\cos(q+1/2)\phi+\sin(q+1/2)\phi\right].$

An upper bound for the sum can be found by replacing $q^{1/2}$ with
$(Q+1)^{1/2}$ and, upon rearrangement,

$\displaystyle\epsilon_{Q}$
$\displaystyle\approx\Im\frac{1+j}{R_{\text{mid}}}\frac{\mathrm{e}^{\mathrm{j}\phi/2}}{[\pi(Q+1)\sin\phi]^{1/2}}\sum_{q=Q+1}^{\infty}\left(t\cos\phi\mathrm{e}^{\mathrm{j}\phi}\right)^{q},$
$\displaystyle=\Im\frac{1+j}{R_{\text{mid}}}\frac{\mathrm{e}^{\mathrm{j}\phi/2}}{[\pi(Q+1)\sin\phi]^{1/2}}\frac{\left(t\cos\phi\mathrm{e}^{\phi}\right)^{Q+1}}{1-t\cos\phi\mathrm{e}^{\mathrm{j}\phi}}.$
(5)

As will be seen, this is an accurate estimate of the remainder in the
polynomial expansion of $1/R$ but it is oscillatory as a function of $q$, so
we adopt the more convenient measure of the magnitude rather than the
imaginary part,

$\displaystyle E_{Q}$
$\displaystyle=\frac{1}{R_{\text{mid}}}\left(\frac{2}{\pi\sin\phi}\right)^{1/2}\frac{|t|^{Q+1}}{(Q+1)^{1/2}}\frac{\cos^{Q+1}\phi}{\left[(1-t)^{2}\cos^{2}\phi+\sin^{2}\phi\right]^{1/2}}.$
(6)

We note that (5) could be integrated over $t$ to give an estimate of the total
error in the integral of $r/R$ but this gives an unwieldy expression with
little advantage in implementations. Instead we adopt as error criterion the
absolute value given by (6) with the value of $t$ given by the nearest point
on the element. In particular, when the projection of the field point lies on
the element, i.e. when the triangle in Figure 1 encloses the origin, $t=1$ and
the error estimate for $1/R$ is

$\displaystyle E_{Q}$
$\displaystyle=\frac{1}{R_{\text{mid}}}\left(\frac{2}{\pi\sin^{3}\phi}\right)^{1/2}\frac{\cos^{Q+1}\phi}{(Q+1)^{1/2}}.$
(7)

From the form of the error estimate, the reason for the difficulty in
evaluating near-singular integrals is clear: near the element plane where
$\phi\to 0$, approximation of the integrand by a polynomial, implicit in the
use of Gaussian quadratures, incurs a very large error, even for quite high
order quadratures with large $Q$.

Figure 3: Error in polynomial approximation of $1/R$, $Q=32$,
$r_{\text{mid}}=1/2$, $t=1$: solid line exact error; boxes estimate from (5);
dashed line absolute value from (6) Figure 4: Error in polynomial
approximation of $1/R$, $r_{\text{mid}}=1/2$, $z=0.1$, $t=1,7/8,3/4$: solid
line exact error; boxes estimate from (5); dashed line absolute value from (6)

To minimize the computational burden of using the criterion, it is applied in
the following manner. Given the transformation into coordinates based on the
element plane, the minimum distance $r_{\min}$ from the origin to the triangle
can be determined (see Figure 5) with $r_{\min}\equiv 0$ when the triangle
encloses the origin. The maximum distance to a vertex $r_{\max}$ is found
similarly. Then we set $r_{\text{mid}}=r_{\max}/2$,
$t=(r_{\text{mid}}-r_{\min})/\overline{r}$ and other quantities as above. The
criterion is then applied by computing $E_{Q}$ for $Q=1,2,\ldots$ until
$E_{Q}$ falls below some specified tolerance, and returning the resulting
value of $Q$, the order of polynomial required to compute $1/R$ to the
specified tolerance over the range of the integral. We note that the error
measure here is the maximum error in $1/R$ at any point in the range of
integration, which is quite a stringent, though conservative, measure, but it
will be found that $E_{Q}$ is a useful assessment of the accuracy of
quadrature.

Figure 3 shows the error estimates as a function of $z$ for a test case with a
32nd order polynomial, equivalent to a 16 point Gaussian quadrature. The error
estimate $\epsilon_{Q}$ is seen to be very reliable, and the magnitude $E_{Q}$
does indeed match the amplitude of the error quite closely. Figure 4 shows the
error as a function of $Q$ for fixed $z$ and again the error behavior is
accurately captured by the estimators. Despite the relative simplicity of the
error measures, they give reliable indicators of the accuracy of the
quadrature or of the order of quadrature required for a given tolerance. We
note finally that the quantities used in the error measure are typically
computed as part of the geometric transformations required in generating a
quadrature on an element, so that there is very little overhead in applying
the error estimate.

##  3 Analysis

The problem to be considered is evaluation of the Helmholtz single- and
double-layer potential integrals on a planar triangular element. Integration
is performed after transformation of coordinates such that the triangular
element is defined by vertices $(x_{i},y_{i},0)$ and the field point lies at
$(0,0,z)$. The triangle is then decomposed into up to three triangles each
having a vertex at $(0,0,0)$. The process is shown in Figure 5. The approach
is similar to that taken in a previous analysis for the Laplace potential [4],
though some changes are required to make it suitable for the Helmholtz
problem.

Figure 5: Integration over a general triangle (left) by subdivision into three
triangles centred at the origin (right). The triangle shown dashed in the
exploded view on the right has negative orientation and its contribution is
subtracted from that of the others. The distance $r_{\min}$ is used in
applying the criterion of Section 2. Figure 6: Reference triangle for
integration

Figure 6 shows the basic triangle which is used for the evaluation of the
contributions from the subtriangles of Figure 5. It has one vertex at the
origin, i.e. at the projection of the field point onto the element plane, and
is defined by the lengths of the two sides which meet at the origin, $r_{1}$
and $r_{2}$, and by the angle $\Theta$ between them.

In developing the analysis, we assume that the triangular element conforms to
some reasonable standards of quality, in particular that the edge length is no
greater than some specified fraction of a wavelength, typically between one
sixth and one eighth. This translates into a limit on $k\ell$, where $\ell$ is
a typical edge length. Taking into account the need to deal with triangles
which are larger than the element proper, such as triangle $023$ in Figure 5,
we assume that $k\ell<\pi/2$, which allows us to limit the size of the
expansions which will be employed in evaluating the potential integrals, while
retaining the required accuracy. If necessary, this limit can be increased, at
the expense of extra computational effort and a small increase in stored data.

###  3.1 Basic integrals

Integration is performed on the reference triangle of Figure 6, using the
polar coordinate system $(r,\theta)$. Geometric parameters are defined,

$\displaystyle\phi$
$\displaystyle=\tan^{-1}\frac{r_{1}-r_{2}\cos\Theta}{r_{2}\sin\Theta},\,\overline{r}(\theta)=\frac{s}{\cos(\theta-\phi)},$
$\displaystyle s$ $\displaystyle=r_{1}\cos\phi,\,S^{2}=s^{2}+z^{2},$

and auxiliary variables used in performing the integrations are

$\displaystyle\alpha^{2}=z^{2}/S^{2},\,\overline{R}$
$\displaystyle=\left(\overline{r}^{2}+z^{2}\right)^{1/2}=S\Delta/\cos\theta,\quad\Delta^{2}=1-\alpha^{2}\sin^{2}\theta.$

The integrals to be evaluated are the zeroth and first order derivatives with
respect to $z$ of

$\displaystyle I_{0}=\mathrm{e}^{\mathrm{j}k|z|}I_{0}^{\prime},\quad
I_{x}=\mathrm{e}^{\mathrm{j}k|z|}I_{x}^{\prime},\quad
I_{y}=\mathrm{e}^{\mathrm{j}k|z|}I_{y}^{\prime},$

where the basic integrals are

$\displaystyle I^{\prime}_{0}$
$\displaystyle=\int_{-\phi}^{\Theta-\phi}\int_{0}^{\overline{r}}\frac{\mathrm{e}^{\mathrm{j}k(R-|z|)}}{R}r\,\mathrm{d}r\,\mathrm{d}\theta,$
(8a) $\displaystyle I^{\prime}_{x}$
$\displaystyle=\int_{-\phi}^{\Theta-\phi}\int_{0}^{\overline{r}}\frac{\mathrm{e}^{\mathrm{j}k(R-|z|)}}{R}r^{2}\,\mathrm{d}r\,\cos\theta\,\mathrm{d}\theta,$
(8b) $\displaystyle I^{\prime}_{y}$
$\displaystyle=\int_{-\phi}^{\Theta-\phi}\int_{0}^{\overline{r}}\frac{\mathrm{e}^{\mathrm{j}k(R-|z|)}}{R}r^{2}\,\mathrm{d}r\,\sin\theta\,\mathrm{d}\theta,$
(8c) $\displaystyle\text{with}\,R$
$\displaystyle=\left(r^{2}+z^{2}\right)^{1/2},$

and correspond to the zero and first order source terms required for linear
shape functions on a plane element. The normal derivatives are given by
differentiation with respect to $z$, for example,

$\displaystyle\frac{\partial I_{0}}{\partial z}$
$\displaystyle=\pm\mathrm{j}k\mathrm{e}^{\mathrm{j}k|z|}I_{0}^{\prime}+\mathrm{e}^{\mathrm{j}k|z|}\frac{\partial
I_{0}^{\prime}}{\partial z}=-\frac{\partial I_{0}}{\partial n}$

where the element is oriented such that the normal lies in the positive $z$
direction, and the upper (lower) signs are taken for positive (negative) $z$.

The integrals are evaluated by expanding the complex exponential in a
polynomial approximation, and evaluating term-by-term using analytical
formulae defined by recursion relations, as in previous work [4]. The form of
the approximation for $\exp[\mathrm{j}kx]$ will be considered in Section 3.3,
but for now we write

$\displaystyle\mathrm{e}^{\mathrm{j}k(R-|z|)}$
$\displaystyle\approx\sum_{q=0}^{Q}e_{q}k^{q}(R-|z|)^{q},$ (9) $\displaystyle
e_{q}$ $\displaystyle=c_{q}+\mathrm{j}s_{q},$ $\displaystyle\text{where}\,\sin
x$ $\displaystyle\approx\sum_{q=0}^{Q}s_{q}x^{q},\quad\cos
x\approx\sum_{q=0}^{Q}c_{q}x^{q}.$

Expanding in powers of $k(R-|z|)$ has the advantages of ensuring that the
expansion remains valid for large values of $z$ as $(R-|z|)\to 0$ as
$z\to\infty$, and providing a natural reduction in the number of terms in (9)
for increasing $z$.

Substituting (9) into (8), yields

$\displaystyle I_{0}^{\prime}$
$\displaystyle\approx\sum_{q=0}^{Q}e_{q}K_{q,0},$ (10a) $\displaystyle
I_{x}^{\prime}$ $\displaystyle\approx\sum_{q=0}^{Q}e_{q}K_{q,x},$ (10b)
$\displaystyle I_{y}^{\prime}$
$\displaystyle\approx\sum_{q=0}^{Q}e_{q}K_{q,y},$ (10c)

where

$\displaystyle K_{q,0}$
$\displaystyle=k^{q}\int_{-\phi}^{\Theta-\phi}\int_{0}^{\overline{r}}\left(R-|z|\right)^{q}\frac{r}{R}\,\mathrm{d}r\,\mathrm{d}\theta,$
(11a) $\displaystyle K_{q,x}$
$\displaystyle=k^{q}\int_{-\phi}^{\Theta-\phi}\int_{0}^{\overline{r}}\left(R-|z|\right)^{q}\frac{r^{2}}{R}\,\mathrm{d}r\cos\theta\,\mathrm{d}\theta,$
(11b) $\displaystyle K_{q,y}$
$\displaystyle=k^{q}\int_{-\phi}^{\Theta-\phi}\int_{0}^{\overline{r}}\left(R-|z|\right)^{q}\frac{r^{2}}{R}\,\mathrm{d}r\sin\theta\,\mathrm{d}\theta.$
(11c)

The integrals of (11) can be evaluated analytically using a combination of
recursions and tabulated integrals. The inner integrals are given by,

$\displaystyle
k^{q}\int_{0}^{\overline{r}}\left(R-|z|\right)^{q}\frac{r}{R}\,\mathrm{d}r$
$\displaystyle=\frac{k^{q}(R-|z|)^{q+1}}{q+1}=\frac{S(kS)^{q}}{q+1}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1},$
(12) $\displaystyle
k^{q}\int_{0}^{\overline{r}}\left(R-|z|\right)^{q}\frac{r^{2}}{R}\,\mathrm{d}r$
$\displaystyle=\frac{k^{q}\overline{r}}{q+2}\left(\overline{R}-|z|\right)^{q+1}+\frac{2|z|}{q+2}J_{q},$
$\displaystyle=\frac{sS(kS)^{q}}{q+2}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{1}{\cos\theta}+\frac{2|z|}{q+2}J_{q},$
(13) $\displaystyle J_{q}$
$\displaystyle=k^{q}\int_{(2|z|)^{1/2}}^{\left(\overline{R}+|z|\right)^{1/2}}\left(t^{2}-2|z|\right)^{q+1/2}\,\mathrm{d}t.$

The integral $J_{q}$ can be evaluated using the recursion

$\displaystyle J_{q}$
$\displaystyle=\frac{k^{q}\overline{r}}{2(q+1)}\left(\overline{R}-|z|\right)^{q}-k|z|\frac{2q+1}{q+1}J_{q-1},$
$\displaystyle=\frac{s}{2}\frac{(kS)^{q}}{q+1}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q}\frac{1}{\cos\theta}-k|z|\frac{2q+1}{q+1}J_{q-1},$
(14) $\displaystyle J_{0}$
$\displaystyle=\frac{\overline{r}}{2}-\frac{|z|}{2}\log\frac{\overline{R}+\overline{r}}{|z|},$
$\displaystyle=\frac{s}{2}\frac{1}{\cos\theta}+\frac{|z|}{4}\log\frac{\Delta-\alpha^{\prime}}{\Delta+\alpha^{\prime}},$
(15)

so that all required terms are written in a form suitable for the application
of standard formulae for trigonometric integrals [5, 2.58],

$\displaystyle K_{q,0}$
$\displaystyle=\frac{S(kS)^{q}}{q+1}\int_{-\phi}^{\Theta-\phi}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\,\mathrm{d}\theta,$
(16a) $\displaystyle K_{q,x}$
$\displaystyle=\frac{sS(kS)^{q}}{q+2}\int_{-\phi}^{\Theta-\phi}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\,\mathrm{d}\theta+\frac{2|z|}{q+2}I_{q,c},$
(16b) $\displaystyle K_{q,y}$
$\displaystyle=\frac{sS(kS)^{q}}{q+2}\int_{-\phi}^{\Theta-\phi}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta+\frac{2|z|}{q+2}I_{q,s},$
(16c) $\displaystyle I_{q,c}$
$\displaystyle=\int_{-\phi}^{\Theta-\phi}J_{q}\cos\theta\,\mathrm{d}\theta,\,I_{q,s}=\int_{-\phi}^{\Theta-\phi}J_{q}\sin\theta\,\mathrm{d}\theta.$

The normal derivatives of the integrals can be evaluated by differentiating
terms, yielding,

$\displaystyle\frac{\partial K_{q,0}}{\partial z}$
$\displaystyle=\mp(kS)^{q}\int_{-\phi}^{\Theta-\phi}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{\cos\theta}{\Delta}\,\mathrm{d}\theta,$
(17) $\displaystyle\frac{\partial K_{q,x}}{\partial z}$ $\displaystyle=\mp
s(kS)^{q}\frac{q+1}{q+2}\int_{-\phi}^{\Theta-\phi}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{\cos\theta}{\Delta}\,\mathrm{d}\theta\pm\frac{2}{q+2}I_{q,c}+\frac{2|z|}{q+2}\frac{\partial
I_{q,c}}{\partial z},$ (18) $\displaystyle\frac{\partial K_{q,y}}{\partial z}$
$\displaystyle=\mp
s(kS)^{q}\frac{q+1}{q+2}\int_{-\phi}^{\Theta-\phi}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{\sin\theta}{\Delta}\,\mathrm{d}\theta\pm\frac{2}{q+2}I_{q,s}+\frac{2|z|}{q+2}\frac{\partial
I_{q,s}}{\partial z}.$ (19)

All integrals can then be evaluated using the results of Section A. This gives
a means of evaluating all required expressions for the integrals on the
triangular element, which can then be summed to give the integral over the
initial general triangle.

###  3.2 Hypersingular integral

The results of (3.1) may be used to solve boundary integral problems using a
standard Helmholtz equation. It is often desirable to employ a Burton and
Miller approach [2] to avoid the well-known problem of fictitious resonances
when the wavenumber $k$ in the exterior problem coincides with an eigenvalue
of the interior problem. In this approach, the Helmholtz equation is combined
with its normal derivative to yield a formulation which is numerically valid
for all real wavenumbers, at the expense of requiring the evaluation of
hypersingular integrals of the form

$\displaystyle\frac{\partial^{2}}{\partial
n^{2}}\iint_{A_{e}}f(\xi,\eta)G(\mathbf{x},\mathbf{x}_{1}(\xi,\eta))\,\mathrm{d}\xi\,\mathrm{d}\eta.$

In order to meet continuity requirements, the collocation points in a
hypersingular method must lie strictly within elements, though discontinuous
elements offer a way around this [6], and so we give a result for the zero-
order (constant) element only:

$\displaystyle\frac{\partial^{2}K_{q,0}}{\partial
z^{2}}=\frac{(kS)^{q}}{S}\biggr{[}\alpha$
$\displaystyle\int_{-\phi}^{\Theta-\phi}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{\cos^{3}\theta}{\Delta^{3}}\,\mathrm{d}\theta$
(20) $\displaystyle+(q+1)$
$\displaystyle\int_{-\phi}^{\Theta-\phi}\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{\cos^{2}\theta}{\Delta^{2}}\,\mathrm{d}\theta\biggr{]}.$

###  3.3 Approximation of exponentials

In order to efficiently evaluate the formulae of Section 3.1, we require a
means of selecting the polynomial approximation to the exponential, (9). The
most obvious choice is to truncate the Taylor series for $\mathrm{e}^{x}$ at
some point where the estimated remainder is smaller than a specified tolerance
$\epsilon$. For reasons of efficiency, however, we adopt an “economized”
series which replaces the truncated Taylor series with a polynomial
approximation with minimum deviation and a minimized error over the range
where the polynomial is used. Given that the integral terms are evaluated
using recursion relations, by reducing the number of terms, we also reduce the
chance of numerical error accumulating in moving from term to term.

The economization algorithm is that given by Acton [1, p291–296] and is used
to generate a set of polynomial approximations of $\sin x$ and $\cos x$ over a
range $0\leq x<\Delta x$, to a tolerance $\epsilon$ where $\epsilon$ is the
maximum difference between $\exp[\mathrm{j}x]$ and the polynomial
approximation over the range $0\leq x<\Delta x$. For the calculations of this
paper, $\Delta x=\pi/16,\pi/8,\pi/4,\pi/2$, and $\epsilon=10^{-n}$,
$n=3,6,9,12,15$. In the implementation, a polynomial approximation is chosen
which has the required maximum error less than $\epsilon$ and $k\ell<\Delta
x$. In the case of $\Delta x=\pi/2$, for example, this gives a reduction in
the number of terms required from fifteen for the truncated Taylor series to
eight for the economized polynomial when $\epsilon=10^{-9}$.

###  3.4 Summary of method

The quadrature method of the previous sections can be summarized as follows,
for a triangle $(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3})$ which has been
rotated into the plane $z=0$ and field point $\mathbf{x}=(0,0,z)$, Figure 5:

  1. 1.

determine the closest and furthest points on the triangle boundary and their
radial distances $r_{\min}$ (Figure 5) and $r_{\max}$;

  2. 2.

compute the required order of quadrature $Q$ for polynomial approximation of
$1/R$, Section 2;

  3. 3.

if $Q$ falls below the set limit:

    1. (a)

evaluate the integrals numerically and terminate;

otherwise

    1. (a)

decompose the triangle into up to three sub-triangles centered at the origin;

    2. (b)

for each sub-triangle compute the contribution using the formulae of Section
3.1 and accumulate, taking account of sub-triangle orientation.

##  4 Numerical testing

Figure 7: Sample triangle and field points

As a numerical test of the performance of the method, we use the same test
case as Pourahmadian and Mogilevskaya [9], Figure 7. Four points are selected
in the element plane, as indicated, and we evaluate $I_{0}$ for $k=1$ as a
function of $z$, vertical displacement from the element; results for $I_{x}$
and $I_{y}$ are similar. As a reference for error estimation, we follow
Pourahmadian and Mogilevskaya and use a polar transformation and a $50\times
50$ point Gaussian quadrature, which is accurate to eight significant figures
[9]. The reported error is

$\displaystyle\epsilon(z)$ $\displaystyle=|I_{0}^{(a)}(z)-I_{0}^{(c)}(z)|,$
(21)

with superscripts ‘a’ and ‘c’ denoting ‘analytical’ and ‘computed’ values
respectively. Error is evaluated by specifying the required tolerance
$\epsilon^{(a)}$ in the analytical method and computing the resulting
$\epsilon$. A second set of error calculations are presented by fixing
$\epsilon^{(a)}=10^{-12}$, varying the order of numerical quadrature in the
polar transformation, and computing the resulting error $\epsilon^{(c)}$.
Figure 8 gives $\epsilon$ as a function of $z$ for varying $\epsilon^{(a)}$,
and $\epsilon^{(c)}$ for varying order of Gaussian quadrature, plotted with
$Q$ computed for varying values of $E_{Q}$. Figure 8 shows error data for the
evaluation of $I_{0}$ and Figure 9 for the normal derivative $\partial
I_{0}/\partial z$.

---
Figure 8: Error for $I_{0}(k,z)$ on element of Figure 7 at points 1–4 (top to
bottom); left-hand column: $|I_{0}^{(a)}-I_{0}^{(c)}|$ for
$\epsilon^{(a)}=10^{-3,-6,-9,-12}$ and $50\times 50$ point polar quadrature;
right-hand column $|I_{0}^{(a)}-I_{0}^{(c)}|$ for $\epsilon^{(a)}=10^{-12}$
and $4\times 4$, $8\times 8$, $16\times 16$, $32\times 32$ polar quadrature;
gray curves: $Q$ for $E_{Q}=10^{-3,-6,-9,-12}$

The left-hand column of Figure 8 shows the error estimate for points 1, 2, 3,
and 4 in Figure 7 which correspond respectively to field points whose
projections lie on a vertex, in the interior, on an edge, and outside the
element. The reference integral $I_{0}^{(c)}$ for points 1–3 is computed using
the $50\times 50$ Gaussian quadrature after transformation to polar
coordinates, and that for point 4 using the 175 point symmetric quadrature of
Wandzura and Xiao [14]. Errors are computed with a requested tolerance
$\epsilon^{(a)}=10^{-3,-6,-9,-12}$ and the computed errors reflect both the
accuracy of the analytical method and the conformity to the requested
tolerance. In applications, there is reasonable confidence that the error will
be approximately equal to that requested, without excessive computation.

The right-hand column of Figure 8 presents data relevant to the error estimate
$E_{Q}$ and the accuracy of Gaussian quadrature in this problem. The darker
curves show an error estimate computed as the difference between the
analytical method with a requested tolerance of $10^{-12}$ and numerical
quadrature of varying order. As expected the low-order methods, e.g. $4\times
4$ points, give a larger error and the high-order approach, $32\times 32$
points, gives accuracy comparable to the analytical technique, except for
small values of $z$. In each case, the breakdown of the polynomial
approximation for $1/R$ is apparent in the increase of error as $z\to 0$, most
clearly for the $16\times 16$ quadrature where the error increases markedly
from $z\approx 0.3$. Of interest here is the value of $E_{Q}$ as a criterion
for selecting quadrature rules. The lighter curves show the value of $Q$ found
from (6) with varying values of $E_{Q}$. The curves do indeed predict quite
well the point at which the polynomial approximation to $1/R$ is no longer
accurate and the Gaussian quadrature begins to fail, confirming the
reliability of the measure as a criterion for the selection of quadrature
rules.

---
Figure 9: Error for $\partial I_{0}(k,z)/\partial n$ on element of Figure 7 at
points 1–4 (top to bottom), notation as in Figure 8

Figure 9 gives similar results but for the evaluation of the normal derivative
of the layer potential, also required in BEM calculations. The results are
similar to those in Figure 8 and the discussion of those data carries over to
here, but it is worth noting that though the error behavior of the Gaussian
quadratures is different from that in Figure 8 (compare the results for point
3, for example), the curves of $Q$ still function as a reliable criterion for
selecting a quadrature method.

##  5 Conclusions

An analytical method for the evaluation of potential integrals in boundary
element codes for the Helmholtz equation has been presented and tested. An
error estimator for purely numerical quadrature has been derived and used to
establish a criterion for quadrature method selection. The quadrature method
has been tested and found to be accurate and reliable; the error criterion is
a reliable technique for quadrature selection. We believe that the quadrature
method proposed is a suitable plug-in replacement in BEM codes for the wave
equation where an a priori error estimate for element integrals and an
economical integration are required.

##  Appendix A Basic integrals

The evaluation of the potential integrals requires a number of elementary
integrals which can be computed using results from standard tables combined
with recursions. This appendix contains the results required for the
evaluation of the trigonometric integrals of the main paper, written in terms
of the parameter $\alpha$, $0\leq\alpha<1$, and
$\Delta^{2}=1-\alpha^{2}\sin^{2}\theta$. The results are given as the
indefinite integral, with a separate result where necessary for the in-plane
case $\alpha=0$.

The first basic term is

$\displaystyle\int\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q}\left(\frac{\Delta}{\cos\theta}\right)^{-s}\,\mathrm{d}\theta$
$\displaystyle=\sum_{u=0}^{q}\binom{q}{u}\left(-\alpha\right)^{q}\int\left(\frac{\Delta}{\cos\theta}\right)^{q-s}\,\mathrm{d}\theta,$
(22)

where $s=0,1,2,3$. The terms in the summation are pseudo-elliptic integrals
which can be evaluated using elementary functions and recursion relations [5,
2.58].

Using the transformation $u=\tan\theta$ and noting that
$\Delta/\cos\theta=(1+{\alpha^{\prime}}^{2}\tan^{2}\theta)^{1/2}$

$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{n}\,\mathrm{d}\theta$
$\displaystyle=\alpha^{2}\int\left(\frac{\Delta}{\cos\theta}\right)^{n-2}\,\mathrm{d}\theta+{\alpha^{\prime}}^{2}\int\left(1+{\alpha^{\prime}}^{2}u^{2}\right)^{(n-2)/2}\,\mathrm{d}u,$
(23)

with $\alpha^{\prime}=(1-\alpha^{2})^{1/2}$.

The integral term can be evaluated using the recursion

$\displaystyle\int\left(1+{\alpha^{\prime}}^{2}u^{2}\right)^{(n+2)/2}\,\mathrm{d}u$
$\displaystyle=\frac{u}{n+3}\left(1+{\alpha^{\prime}}^{2}u^{2}\right)^{(n+2)/2}+\frac{n+2}{n+3}\int\left(1+{\alpha^{\prime}}^{2}u^{2}\right)^{n/2}\,\mathrm{d}u,$
(24)

seeding the recursion with

$\displaystyle\int\left(1+{\alpha^{\prime}}^{2}u^{2}\right)^{-2/2}\,\mathrm{d}u$
$\displaystyle=\frac{\tan^{-1}(\alpha^{\prime}u)}{\alpha^{\prime}},$ (25)
$\displaystyle\int\left(1+{\alpha^{\prime}}^{2}u^{2}\right)^{-1/2}\,\mathrm{d}u$
$\displaystyle=\frac{1}{\alpha^{\prime}}\log\left[\left(1+{\alpha^{\prime}}^{2}u^{2}\right)^{1/2}+\alpha^{\prime}u\right],$
(26)

and using

$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-3}\,\mathrm{d}\theta$
$\displaystyle=-\frac{{\alpha^{\prime}}^{2}}{\alpha^{2}}\frac{\sin\theta}{\Delta}+\frac{\sin^{-1}(\alpha\sin\theta)}{\alpha},$
(27)
$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-2}\,\mathrm{d}\theta$
$\displaystyle=\frac{\theta}{\alpha^{2}}-\frac{\alpha^{\prime}}{\alpha^{2}}\tan^{-1}(\alpha^{\prime}\tan\theta),$
(28)
$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-1}\,\mathrm{d}\theta$
$\displaystyle=\frac{\sin^{-1}(\alpha\sin\theta)}{\alpha}.$ (29)

For $\alpha=0$,

$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-3}\,\mathrm{d}\theta$
$\displaystyle=\sin\theta-\frac{\sin^{3}\theta}{3},$ (30)
$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-2}\,\mathrm{d}\theta$
$\displaystyle=\frac{\sin\theta\cos\theta}{2}+\frac{\theta}{2}$ (31)
$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-1}\,\mathrm{d}\theta$
$\displaystyle=\sin\theta.$ (32)

A second, similar, integral is

$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{n}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta$
$\displaystyle=\alpha^{2}\int\left(\frac{\Delta}{\cos\theta}\right)^{n-2}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta+\frac{1}{n}\left(\frac{\Delta}{\cos\theta}\right)^{n},$
(33)

which can be seeded with [5, 2.584]

$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-3}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta$
$\displaystyle=\frac{\cos\theta}{\alpha^{2}\Delta}-\frac{1}{\alpha^{3}}\log\left(\alpha\cos\theta+\Delta\right),$
$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-2}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta$
$\displaystyle=-\frac{1}{\alpha^{2}}\log\Delta,$
$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-1}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta$
$\displaystyle=-\frac{1}{\alpha}\log(\alpha\cos\theta+\Delta),$
$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{0}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta$
$\displaystyle=-\log\cos\theta.$

For $\alpha=0$,

$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-3}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta$
$\displaystyle=-\frac{\cos^{3}\theta}{3},$
$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-2}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta$
$\displaystyle=\frac{\sin^{2}\theta}{2},$
$\displaystyle\int\left(\frac{\Delta}{\cos\theta}\right)^{-1}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta$
$\displaystyle=-\cos\theta.$

In an implementation of the method of this paper, when the required geometric
parameters have been calculated for the reference triangle, and the
appropriate expansion for $\exp[\mathrm{j}kx]$ has been selected, the first
step is to compute the required elementary integrals (23) and (33) using the
initial values and the recursion relations. The computed terms can then be
used in the summations of (3.1) to evaluate the potential integrals.

For convenience, we define

$\displaystyle L_{c}$
$\displaystyle=\int\cos\theta\log\frac{\Delta-\alpha^{\prime}}{\Delta+\alpha^{\prime}}\,\mathrm{d}\theta=\sin\theta\log\frac{\Delta-\alpha^{\prime}}{\Delta+\alpha^{\prime}}+\log\frac{\Delta+\alpha^{\prime}\sin\theta}{\Delta-\alpha^{\prime}\sin\theta}-2\frac{\alpha^{\prime}}{\alpha}\sin^{-1}(\alpha\sin\theta),$
(34) $\displaystyle L_{s}$
$\displaystyle=\int\sin\theta\log\frac{\Delta-\alpha^{\prime}}{\Delta+\alpha^{\prime}}\,\mathrm{d}\theta=-\cos\theta\log\frac{\Delta-\alpha^{\prime}}{\Delta+\alpha^{\prime}}+2\frac{\alpha^{\prime}}{\alpha}\log(\alpha\cos\theta+\Delta),$
(35)

which are readily evaluated using integration by parts. Then,

$\displaystyle I_{q,c}$ $\displaystyle=\int
J_{q}\cos\theta\,\mathrm{d}\theta,$ (36) $\displaystyle I_{q+1,c}$
$\displaystyle=\frac{s(kS)^{q+1}}{2(q+2)}\int\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\,\mathrm{d}\theta-k|z|\frac{2q+3}{q+2}I_{q,c},$
(37) $\displaystyle I_{0,c}$
$\displaystyle=\frac{s}{2}\theta+\frac{|z|}{4}L_{c},$ (38)
$\displaystyle\frac{\partial I_{q+1,c}}{\partial z}$
$\displaystyle=\mp\frac{s}{2S}(kS)^{q+1}\frac{q+1}{q+2}\int\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{\cos\theta}{\Delta}\,\mathrm{d}\theta$
$\displaystyle\mp\frac{2q+3}{q+2}kI_{q,c}-k|z|\frac{2q+3}{q+2}\frac{\partial
I_{q,c}}{\partial z},$ (39) $\displaystyle\frac{\partial I_{0,c}}{\partial z}$
$\displaystyle=\pm
L_{c}\pm\frac{s}{2S}\int\frac{\cos\theta}{\Delta}\,\mathrm{d}\theta,$ (40)

and similarly

$\displaystyle I_{q+1,s}$
$\displaystyle=\frac{s(kS)^{q+1}}{2(q+2)}\int\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{\sin\theta}{\cos\theta}\,\mathrm{d}\theta-k|z|\frac{2q+3}{q+2}I_{q,s},$
(41) $\displaystyle I_{0,s}$
$\displaystyle=-\frac{s}{2}\log\cos\theta+\frac{|z|}{4}L_{s},$ (42)
$\displaystyle\frac{\partial I_{q+1,s}}{\partial z}$
$\displaystyle=\mp\frac{s}{2S}(kS)^{q+1}\frac{q+1}{q+2}\int\left(\frac{\Delta}{\cos\theta}-\alpha\right)^{q+1}\frac{\sin\theta}{\Delta}\,\mathrm{d}\theta$
$\displaystyle\mp\frac{2q+3}{q+2}kI_{q,s}-k|z|\frac{2q+3}{q+2}\frac{\partial
I_{q,s}}{\partial z},$ (43) $\displaystyle\frac{\partial I_{0,s}}{\partial z}$
$\displaystyle=\pm
L_{s}\pm\frac{s}{2S}\int\frac{\sin\theta}{\Delta}\,\mathrm{d}\theta$ (44)

## References

  * [1] F. S. Acton, Numerical methods that work, Mathematical Association of America, 1990.
  * [2] A. J. Burton and G. F. Miller, The application of integral equation methods to the numerical solution of some exterior boundary-value problems, Proceedings of the Royal Society of London. A., 323 (1971), pp. 201–210.
  * [3] A. Carini and A. Salvadori, Analytical integrations in 3D BEM: preliminaries, Computational Mechanics, 28 (2002), pp. 177–185, https://doi.org/10.1007/S00466-001-0278-7.
  * [4] M. Carley, Analytical formulae for potential integrals on triangles, ASME Journal of Applied Mechanics, 80 (2013), https://doi.org/10.1115/1.4007853.
  * [5] I. Gradshteyn and I. M. Ryzhik, Table of integrals, series, and products, Academic, London, 5th ed., 1980.
  * [6] S. Marburg and S. Schneider, Influence of element types on numeric error for acoustic boundary elements, Journal of Computational Acoustics, 11 (2003), pp. 363–386, https://doi.org/10.1142/S0218396X03001985.
  * [7] S. G. Mogilevskaya and D. V. Nikolskiy, The use of complex integral representations for analytical evaluation of three-dimensional BEM integrals—potential and elasticity problems, Quarterly Journal of Mechanics and Applied Mathematics, 67 (2014), pp. 505–523, https://doi.org/10.1093/qjmam/hbu015.
  * [8] J. N. Newman, Distributions of sources and normal dipoles over a quadrilateral panel, Journal of Engineering Mathematics, 20 (1986), pp. 113–126, https://doi.org/10.1007/BF00042771.
  * [9] F. Pourahmadian and S. G. Mogilevskaya, Complex variables-based approach for analytical evaluation of boundary integral representations of three-dimensional acoustic scattering, Engineering Analysis with Boundary Elements, 53 (2015), pp. 9–17, https://doi.org/10.1016/j.enganabound.2014.11.019.
  * [10] A. Salvadori, Analytical integrations in 3D BEM for elliptic problems: Evaluation and implementation, International Journal for Numerical Methods in Engineering, 84 (2010), pp. 505–542, https://doi.org/10.1002/nme.2906.
  * [11] C. Schwab and W. L. Wendland, On numerical cubatures of singular surface integrals in boundary element methods, Numerische Mathematik, 62 (1992), pp. 343–369.
  * [12] J.-C. Suh, The evaluation of the Biot–Savart integral, Journal of Engineering Mathematics, 37 (2000), pp. 375–395, https://doi.org/10.1023/A:1004666000020.
  * [13] A. Tadeu and J. António, 3D acoustic wave simulation using BEM formulations: Closed form integration of singular and hypersingular integrals, Engineering Analysis with Boundary Elements, 36 (2012), pp. 1389–1396, https://doi.org/10.1016/j.enganabound.2012.03.011.
  * [14] S. Wandzura and H. Xiao, Symmetric quadrature rules on a triangle, Computers and Mathematics with Applications, 45 (2003), pp. 1829–1840.


††thanks: A footnote to the article title

# Polarization-based cyclic weak value metrology for angular velocity
measurement

Zi-Rui Zhong  Yue Chen  Wei-Jun Tan  Xiang-Ming Hu  Qing-Lin Wu
<EMAIL_ADDRESS>_Department of Physics, Central China Normal University,
Wuhan 430079, China_ _Key Laboratory of Quark & Lepton Physics (MOE) and
Institute of Particle Physics, Central China Normal University, Wuhan 430079,
China_

###### Abstract

Weak value has been proved to amplify the detecting changes of the meters at
the cost of power due to post-selection. Previous power-recycling schemes
enable the failed post-selection photons to be reselected repeatedly, thus
surpassing the upper noise limit and improving the precision of
interferometric systems. Here we introduce three cyclic methods to improve the
sensitivity of polarization-based weak-value-based angular velocity
measurement: power-, signal- and dual-recycling schemes. By inserting one or
two partially transmitting mirrors inside the system, both the power and
precision of detected signals are greatly enhanced, and the dual-recycling
scheme has wider optimal region than that of power- or signal-recycling
schemes. Compared to non-polarization schemes, polarization-based schemes
enjoy lower optical loss and unique cyclic directions. These reduce the
crosstalk among different paths of light and, theoretically, eliminate the
walk-off effect, thus towering in both theoretical performance and
application.

††preprint: APS/123-QED

##  I Introduction

Since first introduced by Aharonov, Albert and Vaidman(AAV) in Ref.[1], weak
measurement has shown its numerous potential in various precise measurements.
Unlike the classical(or strong) measurements set forth by von Neumann[2], the
coupling between the probe and system is very weak, leading to the famous
measurement disturbance guaranteed by Heisenberg limit invalidated[3].
Therefore, it can be used to reconsider many interesting quantum phenomena
such as Hardy’s paradox[4, 5, 6, 7], three-box problem[8, 9, 10] and quantum
Cheshire cats[11, 12, 13, 14, 15]. By preparing appropriate pre- and post-
selected states, a weak measurement enables a so called ‘weak value’ to record
the information of the weak interaction process. Generally, the weak value is
defined as $A_{w}=\left\langle
f\right|\hat{A}\left|i\right\rangle/\left\langle f|i\right\rangle$, where
$\left|i\right\rangle$, $\left|f\right\rangle$ are the pre-, post-selected
states and $\hat{A}$ is the measured observable. Because $\left\langle
f|i\right\rangle$ exists in the denomination of the formula, $A_{w}$ can be
very large if $\left|i\right\rangle$ and $\left|f\right\rangle$ are nearly
orthogonal. Thus, it has the potential to detect many small physical effects
such as the spin Hall effect[16, 17, 18], Goos-Hänchen shift[19, 20], beam
deflection[21], velocity[22], phase shifts[23, 24, 25], temperature[26],
angular velocity[27, 28, 29] and resonance[30], to name a few.

However, the weak-value-amplification (WWA) effect comes at the sacrifice of
the post-selection probability, which is defined as $P=\left|\left\langle
f\middle|i\right\rangle\right|^{2}$. This results in only a small part of the
incident light being detected while discarding the rest most. To solve this
problem, the power-recycling technique[31, 32, 33, 34, 35] is introduced by
placing a partially transmitting mirror(PTM) at the bright port of the
interferometer. It reuses the failed post-selection photons and can permit all
the input light to be detected if providing ideal experimental conditions.
Besides, this enhancement enables the amplification of SNR itself by large
weak value factor, thus breaking the upper limits of classical measurements.
The similar conclusion is obtained for the signal-recycling weak
measurement[36], which works by placing the PTM at the dark port of the
interferometer. In the power- and signal-recycling system above, the failed
and successful post-selected photons are reused, respectively. Furthermore,
these two parts can be combined in one system for further improvement, denoted
by dual-recycling scheme[37, 38, 39, 40].

Previous dual-recycled interferometric WWA setup obtains large precision
improvement while sacrificing some of WWA effect of pointer due to the walk-
off effect. In addition, the walkway of cyclic photons inside the
interferometer is intricate, resulting in inevitable crosstalk, which
increases the system loss. Here we introduce above-mentioned cyclic schemes of
polarization-based WWA set up based on the angular velocity measurement of
[41] . Compared with the non-polarization cyclic schemes, we substitute the
polarization-beam-splitter (PBS) for beam-splitter(BS). This simplifies the
direction of light paths, which is only clockwise $\circlearrowright$ in this
article, and reduces the optical loss. In addition, the ’only-direction’
permits a filter to refresh all cyclic photons before their last weak
interaction, thus eliminating the walk-off effect.

Figure 1: (Color online) (a)Schematic of weak-value-based angular velocity
measurement. A laser wave is generated by an acoustic optical modulator(AOM)
and enters a polarization-dependent angular velocity measurement system
consisting of two Glan prisms(G1 and G2) and the PHWP. The small angular
velocity $\omega$ is induced by the PHWP and finally measured by the detector.
(b): The power-recycling scheme. The PBS, which distinguishes the polarization
between H and V, is combined with the PTM to reuse the failed post-selection
photons repeatedly. Two QWPs near the PTM provide the initial polarization and
rotate the polarization direction of cyclic light back to V. The filter in
front of the PHWP refreshes the beam profile on each pass. (c)The signal-
recycling scheme. The he PTM at the output is combined with the previous
system to form a signal-recycling cavity, thus improving the detected signal.
(d)The dual-recycling scheme. Combining the power- and signal-recycling PTMs
into one scheme to further enhance the precision of detecting. QWP: quarter-
wave plate. HWP: half-wave plate. PBS: polarization-beam-splitter. PHWP:
piezo-driven half-wave plate. PTM: partially transmitting mirror. H:
horizontal. V: vertical.

##  II Standard WWA setup

We first review the standard WWA setup for angular velocity measurement in
[41]. As shown in Fig. 1(a), a non-Fourier limit Gaussian pulse
$I_{0}\left(t\right)=\left(N^{2}/2\pi\tau^{2}\right)^{1/2}exp\left(-t^{2}/2\tau^{2}\right)$
, where $N$ is the number of photons and $\tau$ is the length of pulse, is
sent to a polarization-dependent system. The first Glan prism(G1) combined
with half-wave plate(HWP) provide the pre-selection, where the axis of G1 is
vertical and the angle between G1 and HWP is $\phi_{1}$. The second Glan
prism(G2), whose original orientation of axis is the horizon, provides to the
post-selected state. The weak interaction is expressed by a unitary operator
${\hat{U}}_{w}=exp\left(i2\omega t\hat{A}\right)$ , where $\omega$ is angular
velocity by rotating HWP and $\hat{A}$ is a Hermitian operator
$\hat{A}=\left|\left.L\right\rangle\right.\left\langle\left.L\right|\right.-\left|\left.R\right\rangle\right.\left\langle\left.R\right|\right.$
($\left|\left.L\right\rangle\right.$, $\ \left|\left.R\right\rangle\right.$
are the left, right circularly polarized states). Here we introduce a unitary
operator ${\hat{U}}_{\phi}=exp\left(i2\phi_{1}\hat{A}\right)$ to represent the
polarization rotation $\phi_{1}$ produced by the piezo-driven half-wave
plate(PHWP). In this way, the pre-selected state of system is
$\left|\left.\psi_{pre}\right\rangle\right.={\hat{U}}_{\phi}\left|\left.V\right\rangle\right.=\frac{i}{\sqrt{2}}\left[exp\left(-i2\phi_{1}\right)\left|\left.R\right\rangle-
exp\left(i2\phi_{1}\right)\left|\left.L\right\rangle\right.\right.\right]]$
and the post-selected state is
$\left.\left|\psi_{pos}\right.\right\rangle=cos\phi_{2}\left|\left.H\right\rangle\right.-sin\phi_{2}\left.\left|V\right.\right\rangle=\frac{i}{\sqrt{2}}\left[exp\left(-i2\phi_{2}\right)\left|\left.R\right\rangle+exp\left(i2\phi_{2}\right)\left|\left.L\right\rangle\right.\right.\right]$
. We assume $\left|\left.\phi_{0}\right\rangle\right.$ is the infinite state
of probe. Therefore, the intensity of detected light is given by

$\displaystyle I_{d}\left(t\right)$
$\displaystyle=\left|\left\langle\psi_{pos}\right|\hat{U}_{w}\left|\psi_{pre}\right\rangle\otimes\left|\left.\varphi_{0}\right\rangle\right.\right|^{2}$
(1) $\displaystyle\approx\frac{N}{\sqrt{2\pi\tau^{2}}}{sin}^{2}\phi
exp\left[-\frac{1}{2\tau^{2}}\left(t-\frac{4\omega\tau^{2}}{\phi}\right)^{2}\right],$

where we introduce the angle $\phi=2\phi_{1}-\phi_{2}$ and assume
$2\omega\tau\ll\phi\ll 1$. The corresponding weak value is given by
$A_{w}=\left\langle\psi_{pos}\right|\hat{A}\left|\psi_{pre}\right\rangle/\left\langle\psi_{pos}|\psi_{pre}\right\rangle$,
which is related to the time shift. Compared with the incident light, the time
shift induced by the PHWP is $\delta
t=\frac{4\omega\tau^{2}}{\phi}=4\omega\tau^{2}\left|A_{w}\right|$. Based on
Fisher information theory[42, 43], the fisher information of time shift
$\delta t$ is $F\left(\delta t\right)=\int{dtI_{d}\left|\frac{d}{d\delta
t}lnI_{d}\right|^{2}}\approx\frac{N\phi^{2}}{\tau^{2}}$. So the minimum
uncertainty of angular velocity determined by Cramér-Rao bound(CRB) satisfies

$\Delta\omega_{CRB}=\frac{\phi}{4\tau^{2}}\Delta\left(\delta
t\right)\approx\frac{1}{4\sqrt[]{N}\tau}.$ (2)

Therefore, the corresponding SNR is given by

$SNR=\frac{\omega}{\left(\Delta\omega\right)}=4\omega\tau\sqrt{N}=\frac{\sqrt{N}\phi}{\tau}\delta
t.$ (3)

##  III Power-recycling

The power-recycled weak-value setup is shown in Fig. 1(b). The initial state
$\left|\left.V\right\rangle\right.$ is provided by the G and two QWP, and we
use the combination ‘HWP1 PBS HWP2’ to replace the G2 for the post-selection.
The PTM, whose reflection and transmission coefficients are $r$ and $p$
$\left(r^{2}+p^{2}=1\right)$, is placed between two QWP, thus reflecting the
failed post-selected light while rotating the light polarization to
$\left|\left.V\right\rangle\right.$. In this angular velocity measurement
scheme, the initial light state can be expressed as
$\left|\left.\varphi_{0}\right\rangle\right.=\sqrt{I_{0}\left(t\right)}\left|\left.\alpha\right\rangle\right.=\left(N^{2}/2\pi\tau^{2}\right)^{1/4}exp\left(-t^{2}/4\tau^{2}\right)\left|\left.\alpha\right\rangle\right.$,
where $\left|\left.\alpha\right\rangle\right.$ is coherent light state. Here
we define two orthogonal states $\left|\left.\psi_{1}\right\rangle\right.$ and
$\left|\left.\psi_{2}\right\rangle\right.$ to represent the input and output
system states, where $\left|\left.\psi_{1}\right\rangle\
=\left|\left.V\right\rangle\right.\right.=\frac{i}{\sqrt{2}}\left(\left|\left.R\right\rangle-\left|\left.L\right\rangle\right.\right.\right)$
and $\left|\left.\psi_{2}\right\rangle\
=\left|\left.H\right\rangle\right.\right.=\frac{1}{\sqrt{2}}\left(\left|\left.R\right\rangle+\left|\left.L\right\rangle\right.\right.\right)$.
Post-selected by the input and output ends, the meter states become
$\left|\left.\varphi_{ref}\right\rangle\right.=\left\langle
V\right|\hat{U_{w}}\hat{U_{\phi}}\left|V\right\rangle\left|\left.\varphi_{0}\right\rangle\right.$
and $\left|\left.\varphi_{out}\right\rangle\right.=\left\langle
H\right|\hat{U_{w}}\hat{U_{\phi}}\left|V\right\rangle\left|\left.\varphi_{0}\right\rangle\right.$
, respectively. This produces two measurement operators
$M_{11}=\left\langle\psi_{1}\right|{\hat{U}}_{w}{\hat{U}_{\phi}}\left|\psi_{1}\right\rangle=\cos{\left(\phi+2\omega
t\right)}$ and
$M_{12}=\left\langle\psi_{2}\right|{\hat{U}}_{w}{\hat{U}_{\phi}}\left|\psi_{1}\right\rangle=\sin{\left(\phi+2\omega
t\right)}$. We introduce the non-unitary operator $\hat{L}=\sqrt{1-\gamma}$,
where $\gamma$ is the single-pass power loss, to express the loss of optical
imperfection in one return. Assuming the length of one traversal is $l_{cav}$,
the pulse transition time of per traversal is given by $t_{cav}=2l_{cav}/c$.
Generally, both the measurement operators and the meter state are related to
the number of traversals. For example, $M_{11}$ should be written as
${M_{11}}^{n}=\cos{\left[\phi+2\omega\left(t-nt_{cav}\right)\right]}$. [35]
proved that this change is small and only induces a constant delay, which can
be eliminated. Therefore, with the resonance cavity, the amplitude of the
detected signal is given by the sum of amplitude from all traversal numbers,

$\left|\left.\varphi_{d}\right\rangle\right._{pow}=pM_{12}\sum_{n=0}^{\infty}{\left(rLM_{11}\right)^{n}\left|\left.\varphi_{0}\right\rangle\right.}$
(4)

It is a summation of the convergence series so that there is a maximum value
of $n$, denoted by $n_{max}$. Therefore, the formula above can be simplified
as

$\displaystyle\left|\left.\varphi_{d}\right\rangle\right._{pow}$
$\displaystyle=pM_{12}\sum_{n=0}^{n_{max}}{\left(rLM_{11}\right)^{n}\left|\left.\varphi_{0}\right\rangle\right.}$
(5)
$\displaystyle\approx\left(\frac{N^{2}}{2\pi\tau^{2}}\right)^{\frac{1}{4}}exp\left(-\frac{t^{2}}{4\tau^{2}}\right)\frac{p\sin{\left(\phi+2\omega
t\right)}}{1-rL\cos{\left(\phi+2\omega t\right)}}$

Next, we do a Taylor expansion on the function
$f\left(x\right)=p\sin{\left(\phi+2\omega
t\right)}/\left[1-rL\cos{\left(\phi+2\omega t\right)}\right]$ and make an
approximation $f\left(x\right)\approx
f\left(0\right)+xf^{\prime}\left(0\right)\approx
exp\left(-f^{\prime}\left(0\right)x/f\left(0\right)\right)$. Then the
amplitude of detected state is

$\left|\left.\varphi_{d}\right\rangle\right._{pow}\approx
A\left(\frac{N^{2}}{2\pi\tau^{2}}\right)^{\frac{1}{4}}\sin^{2}{\left(\phi+2\omega
t\right)}exp\left[-\frac{\left(t-\delta t_{p}\right)^{2}}{4\tau^{2}}\right]$
(6)

where

$A=\frac{p}{1-r\sqrt{1-\gamma}\cos{2\phi}}$ (7)

and

$\delta
t_{p}=\frac{2\omega\tau^{2}\left(\cos{\phi}-r\sqrt{1-\gamma}\right)}{\sin{\phi}\left(1-r\sqrt{1-\gamma}\cos{2\phi}\right)}.$
(8)

Due to the walk-off effect, the time shift changes from $\delta t$ to $\delta
t_{p}$. But if we place a filter in front of the PHWP, each time reflected by
the PTM, the light pass through the filter and is projected into
$\left|\left.\varphi_{0}\right\rangle\right.$. This leaves the pre-filter
state as
$\left|\left.\varphi\prime\right\rangle\right._{pow}=M_{11}\left|\left.\varphi_{0}\right\rangle\right./\sqrt{\left|M_{11}\left|\left.\varphi_{0}\right\rangle\right.\right|^{2}}$.
So the probability of surviving the filter is

$\displaystyle p_{f}=$
$\displaystyle\left|\left\langle\varphi_{0}\middle|\varphi\prime\right\rangle\right|^{2}=\frac{\cos^{2}{\phi}}{\sinh{4\omega^{2}\tau^{2}}+\cos^{2}{\phi}e^{-4\omega^{2}\tau^{2}}}$
(9) $\displaystyle\approx
1-\phi^{2}\left(4\omega^{2}\tau^{2}\right)-\left(4\omega^{2}\tau^{2}\right)^{2}/2+\cdots,$

where we make the approximation in the weak value range,
$2\omega\tau\ll\phi\ll 1$. In this way, the time shift is refreshed every
cycle, eliminating the walk-off effect while adding a minimum ‘filter’ loss
$\gamma_{min}\approx 4\omega^{2}\tau^{2}\phi^{2}$ to the system[31].
Therefore, the power of the detected signal is given by

$I_{pow}\approx\left(\frac{N^{2}}{2\pi\tau^{2}}\right)^{\frac{1}{2}}A^{2}\sin^{2}{\left(\phi+2\omega
t\right)}exp\left[-\frac{\left(t-\delta t\right)^{2}}{2\tau^{2}}\right].$ (10)

Similarly, the corresponding SNR determined by Cramér-Rao bound(CRB) is

${\rm SNR}_{pow}\approx A\frac{\sqrt{N}\phi}{\tau}\delta t,$ (11)

which is A times to the SNR of standard weak measurement.

Figure 2: (Color online) Comparison of power-, signal- and dual-recycling
schemes. (a), (b) and (c) correspond to the case of A and B varying with $r$
under different value of $\gamma$: 0, 0.1 and 0.2 where $\phi=0.1$. In (d),
(e) and (f), we assume $r=0.9$ and plot A, B varying with $\phi$ under
$\gamma=0,\ 0.1,\ 0.2$, respectively. A and B: the improvement factor of power
(or signal) and dual recycling schemes. $r$: the reflection coefficient of
PTM. $\gamma$: optical loss.

##  IV Signal-recycling

The similar methods are obtained in the signal-recycled scheme. As shown in
Fig. 1(c), the optical axis of Glan prism is vertical, promising the input
state $\left|\left.V\right\rangle\right.$. The post-selection is provided by
the combination ‘HWP PBS QWPs’. The PTM is placed between the two QWPs to
reflect the output signal. This post-selection processes provide two
measurement operators:
$M_{12}=\left\langle\psi_{2}\right|{\hat{U}}_{w}{\hat{U}_{\phi}}\left|\psi_{1}\right\rangle=\sin{\left(\phi+2\omega
t\right)}$ and
$M_{22}=\left\langle\psi_{2}\right|{\hat{U}}_{w}{\hat{U}_{\phi}}\left|\psi_{2}\right\rangle=\cos{\left(\phi+2\omega
t\right)}$. In this signal-recycled cavity, the amplitude of the detected
signal is given by

$\displaystyle\left|\left.\varphi_{d}\right\rangle\right._{sig}$
$\displaystyle=pM_{12}\sum_{n=0}^{n_{max}}{\left(rLM_{22}\right)^{n}\left|\left.\varphi_{0}\right\rangle\right.}$
(12)
$\displaystyle\approx\left(\frac{N^{2}}{2\pi\tau^{2}}\right)^{\frac{1}{4}}exp\left(-\frac{t^{2}}{4\tau^{2}}\right)\frac{p\sin{\left(\phi+2\omega
t\right)}}{1-rL\cos{\left(\phi+2\omega t\right)}}$

which is the same to $\left|\left.\varphi_{d}\right\rangle\right._{pow}$. This
means the power- and signal-recycling are of the same status in the weak-
value-based signal improvement. Different from previous interferometric
signal-recycling scheme, the only clockwise path permits all cyclic photons to
be refreshed before the next post-selection. With the filter, the pre-filter
state is
$\left|\left.\varphi\prime\right\rangle\right._{sig}=M_{22}\left|\left.\varphi_{0}\right\rangle\right./\sqrt{\left|M_{22}\left|\left.\varphi_{0}\right\rangle\right.\right|^{2}}=\left|\left.\varphi\prime\right\rangle\right._{pow}$.
We can predict that the same conclusion is taken in the calculation of the
detected power and SNR, $I_{sig}=I_{pow}$ and ${SNR}_{sig}={SNR}_{pow}$

Figure 3: (Color online) Comparison of polarization-based and interferometric
dual recycling schemes. (a), (b) and (c) correspond to the case of $B$ varying
with $r$ and $\phi$ under different value of $\gamma$: 0, 0.1 and 0.2 in 3D,
respectively. (d), (e) and (f) respresent that $B_{non}$ varies with $r$ and
$\phi$ under different value of $\gamma$: 0, 0.1 and 0.2 in 3D. $B$ and
$B_{non}$: the improvement factor of polarization-based and non-polarization
dual recycling schemes. $r$: the reflection coefficient of PTM. $\phi$: post-
selected angle. $\gamma$: optical loss.

##  V Dual-recycling

Because the power- and signal-recycling mirrors reuse different parts of the
photon respectively, it seems that the two PTMs can be compatible in one
system. In Fig. 1(d), the Glan prism together with two QWPs provide the input
state $\left|\left.\psi_{1}\right\rangle\ \right.$, and the output state
$\left|\left.\psi_{2}\right\rangle\right.$ is provided by the combination ‘HWP
PBS QWPs’. In this cyclic process, all possible post-selections are available.
Similarly, the filter in front of the PHWP projects the meter state into
$\left|\left.\varphi_{0}\right\rangle\right.$, thus eliminating the walk-off
effect and maintaining the large point shift associated with the WVA. This
also results in a minimum optical loss $\gamma_{min}\approx
4\omega^{2}\tau^{2}\phi^{2}$, which can be ignored in the weak value range
$2\omega\tau\ll\phi\ll 1$. For simple calculation, we assume the parameters of
PTMs are the same and introduce the measurement matrix

$U=\begin{bmatrix}M_{11}&M_{12}\\\ M_{21}&M_{22}\\\
\end{bmatrix}=\begin{bmatrix}\cos{\left(\phi+2\omega
t\right)}&\sin{\left(\phi+2\omega t\right)}\\\ -\sin{\left(\phi+2\omega
t\right)}&\cos{\left(\phi+2\omega t\right)}\\\ \end{bmatrix}$ (13)

which is formed by four measurement operators and arranged in the order
corresponding to the subscripts. $\left(U^{n}\right)_{12}$ represents the
physical process that the incident light travels through the dual recycling
cavity $n$ times and finally reaches the detector. Therefore, the steady state
amplitude detected by the meter is given by the sum of amplitude of all
traversal numbers,

$\displaystyle\left\langle\alpha\right.\left|\left.\varphi_{d}\right\rangle\right.=p\sum_{n=0}^{n_{max}}{\left(\sqrt{1-\gamma}\right)^{n+1}\left(U^{n+1}\right)_{12}p\left\langle\alpha\right.\left|\left.\varphi_{0}\right\rangle\right.}$
(14) $\displaystyle\approx
p^{2}\left(\frac{N^{2}}{2\pi\tau^{2}}\right)^{\frac{1}{4}}exp\left(-\frac{t^{2}}{4\tau^{2}}\right)\hat{L}\left(\frac{U}{I-r\sqrt{1-\gamma}U}\right)_{12}$
$\displaystyle=-\frac{p^{2}\left(\frac{N^{2}}{2\pi\tau^{2}}\right)^{\frac{1}{4}}exp\left(-\frac{t^{2}}{4\tau^{2}}\right)\sqrt{1-\gamma}\sin{\left(\phi+2\omega
t\right)}}{1+\left(1-\gamma\right)r^{2}-2\sqrt{1-\gamma}rcos\left(\phi+2\omega
t\right)}$
$\displaystyle\approx-\frac{p^{2}\left(\frac{N^{2}}{2\pi\tau^{2}}\right)^{\frac{1}{4}}exp\left(-\frac{t^{2}}{4\tau^{2}}\right)\sqrt{1-\gamma}\sin{\left(\phi+2\omega
t\right)}}{1+\left(1-\gamma\right)r^{2}-2\sqrt{1-\gamma}r\cos{\phi}},$

Where the last approximation is taken with the minimum ‘filter’ loss
$\gamma_{min}\approx 4\omega^{2}\tau^{2}\phi^{2}$. In this way, the intensity
of the detected signal is given by

$I_{dua}\approx\left(\frac{N^{2}}{2\pi\tau^{2}}\right)^{\frac{1}{2}}B^{2}\sin^{2}{\left(\phi+2\omega
t\right)}exp\left[-\frac{\left(t-\delta t\right)^{2}}{2\tau^{2}}\right]$ (15)

where

$B=\frac{p^{2}}{1+\left(1-\gamma\right)r^{2}-2\sqrt{1-\gamma}r\cos{\phi}}.$
(16)

Thus, the corresponding SNR is

${\rm SNR}_{pow}\approx B\frac{\sqrt{N}\phi}{\tau}\delta t.$ (17)

##  VI Comparison

Here we define $A$ and $B$ are the improvement factors of power(or signal) and
dual recycling, respectively. It is clear that the power-, signal- and dual-
recycling schemes improve the SNR of standard WWA setup $A$, $A$ and $B$
times. $A^{2}$ and $B^{2}$ also corresponds to the improvement of the detected
power. Therefore, as shown in Fig. 2, we plot $A$ and $B$ varying with
$r$(Fig. 2(a), (b) and (c)) or $\phi$(Fig. 2(d), (e) and (f)) under different
loss $\gamma=0,\ 0.1,\ 0.2$, which correspond to ideal, low and regular loss.

As expected, both $A$ and $B$ can reach the maximum value $1/\phi=10$, as
shown in Fig. 2(a), (b) and (c). However, the peak of $A$ decreases faster
that of $B$, corresponding to a larger limitation to the improvement of
detecting. Therefore, the dual-recycling cavity has tolerance for a wider
range of $r$ and $\gamma$, which applies more to real experiments. In Fig.
2(d), (e) and (f), we set $r=0.9$, a common parameter of the PTM, and draw the
curves of $A$ and $B$ varying with $\phi$ where $\phi\in\left[0.01,\
0.2\right]$. We can see that the improvement factor of dual recycling is
larger in most weak value range, thus suppressing the power or signal
recycling.

In previous dual-recycled interference-based WWA system[40], the amplification
effect of pointer is reduced by the walk-off effect, leading to a limitation
of the precision gain. From the equations (25) and (26) in [40], without a
filter, the improvement factor $B$ changes to $B_{non}$

$B_{non}=\xi\frac{p^{2}}{1+\left(1-\gamma\right)r^{2}-2\sqrt{1-\gamma}r\cos{\phi}},$
(18)

where

$\xi=\phi\frac{\cos{\phi\left[1+r^{2}\left(1-\gamma\right)\right]}-2r\sqrt{1-\gamma}}{\sin{\phi}[1+r^{2}\left(1-\gamma\right)-2r\sqrt{1-\gamma}\cos{\phi}]}.$
(19)

Due to the proper using of the filter, a minimum filter loss replaces original
performance reduction in the polarization-based dual recycling scheme. For
clear comparison, in Fig. 3, we similarly set $\phi=0.1$ and plot $B$,
$B_{non}$ vary with $r$ and $\phi$ under different losses. The polarization-
based scheme has obvious improvement and the gaps between $B_{non}$ and $B$
decrease as the loss increases. This is established on the assumption that
both systems enjoy the same loss $\gamma$. Actually, replacing the BS with the
PBS can effectively reduce the optical loss. The probability of photons
surviving the PBS ($\geq 95\%$) is known to be larger than that of the BS
($\geq 90\%$). In addition, the PBS simplifies the propagation paths of cyclic
photons, which reduces the crosstalk among photons. All these reasons make
this polarization-based scheme advantageous in both theoretical performance
and experimental application.

##  VII Conclusion

In summary, we have proposed three polarization-based cyclic weak measurement
schemes based on the angular velocity weak measurement setup. By inserting one
or two PTMs into the system to form a resonant cavity, all the incident light
can be detected in principle. In our analysis, this polarization-based schemes
can suppress the pervious interferometric schemes by their lower theoretical
loss and improved cyclic paths. This permits a filter to refresh the cyclic
photons, thus eliminating the walk-off effect. In addition, we get the similar
conclusion about the improvement factor, all three schemes can greatly improve
the power and SNR of the detected signal. Especially, the dual-recycling
scheme enjoys the wider optimal region.

However, the beam will be a diffracting Gaussian beam with a waist instead of
the parallel beam treated above. Similar to the solutions of [32], several
lenses should be placed well-designed to form stable self-reproduction of
Gaussian mode. In addition, the cavity is unstable due to the influence of
optical platform jitter, temperature, pressure and so on. Therefore, the
Pound-Drever-Hall(PDH) system is also essential to feedback the instability,
locking the lengths of cavity.

This cyclic mode is also applied to other WVA experimental realizations
because the post-selection absolutely exists in all weak value setup. In
addition, the quantum resources can be used to increase the precision beyond
the standard quantum limit[43, 44], which is a predictable way to further
improve the performance of weak value metrology.

##  VIII acknowledgements

This work was supported by the National Natural Science Foundation of China
(Grants No. 11734015).

## References

  * Aharonov  _et al._ [1988] Y. Aharonov, D. Z. Albert, and L. Vaidman, How the result of a measurement of a component of the spin of a spin-1/2 particle can turn out to be 100, Phys. Rev. Lett. 60, 1351 (1988).
  * von Neumann [2018] J. von Neumann,  _Mathematical Foundations of Quantum Mechanics_, edited by N. A. Wheeler (Princeton University Press, Princeton, 2018).
  * Steinberg [2010] A. Steinberg, A light touch, Nature 463, 890–891 (2010).
  * Hardy [1992] L. Hardy, Quantum mechanics, local realistic theories, and lorentz-invariant realistic theories, Phys. Rev. Lett. 68, 2981 (1992).
  * Aharonov  _et al._ [2002] Y. Aharonov, A. Botero, S. Popescu, B. Reznik, and J. Tollaksen, Revisiting hardy’s paradox: counterfactual statements, real measurements, entanglement and weak values, Physics Letters A 301, 130 (2002).
  * Lundeen and Steinberg [2009] J. S. Lundeen and A. M. Steinberg, Experimental joint weak measurement on a photon pair as a probe of hardy’s paradox, Phys. Rev. Lett. 102, 020404 (2009).
  * Yokota  _et al._ [2009] K. Yokota, T. Yamamoto, M. Koashi, and N. Imoto, Direct observation of hardy’s paradox by joint weak measurement with an entangled photon pair, New Journal of Physics 11, 033011 (2009).
  * Aharonov and Vaidman [1991] Y. Aharonov and L. Vaidman, Complete description of a quantum system at a given time, Journal of Physics A: Mathematical and General 24, 2315 (1991).
  * Ravon and Vaidman [2007] T. Ravon and L. Vaidman, The three-box paradox revisited, Journal of Physics A: Mathematical and Theoretical 40, 2873 (2007).
  * Resch  _et al._ [2004] K. Resch, J. Lundeen, and A. Steinberg, Experimental realization of the quantum box problem, Physics Letters A 324, 125 (2004).
  * Aharonov  _et al._ [2013] Y. Aharonov, S. Popescu, D. Rohrlich, and P. Skrzypczyk, Quantum cheshire cats, New Journal of Physics 15, 113015 (2013).
  * Aharonov  _et al._ [2021] Y. Aharonov, E. Cohen, and S. Popescu, Nature Communications 12, 4770 (2021).
  * Denkmayr  _et al._ [2014] T. Denkmayr, H. Geppert, S. Sponar, _et al._ , Nature Communications 5, 4492 (2014).
  * Duprey  _et al._ [2018] Q. Duprey, S. Kanjilal, U. Sinha, D. Home, and A. Matzkin, The quantum cheshire cat effect: Theoretical basis and observational implications, Annals of Physics 391, 1 (2018).
  * Quach [2019] J. Q. Quach, Dual to the anomalous weak-value effect of photon-polarization separation, Phys. Rev. A 100, 052117 (2019).
  * Hosten and Kwiat [2008] O. Hosten and P. Kwiat, Science 319, 787 (2008).
  * Zhou  _et al._ [2012] X. Zhou, Z. Xiao, H. Luo, and S. Wen, Phys. Rev. A 85, 043809 (2012).
  * Bai  _et al._ [2020] X. Bai, Y. Liu, L. Tang, Q. Zang, J. Li, W. Lu, H. Shi, X. Sun, and Y. Lu, Opt. Express 28, 15284 (2020).
  * Xu  _et al._ [2013a] X.-Y. Xu, Y. Kedem, K. Sun, L. Vaidman, C.-F. Li, and G.-C. Guo, Phys. Rev. Lett. 111, 033604 (2013a).
  * Li  _et al._ [2018a] L. Li, Y. Li, Y.-L. Zhang, S. Yu, C.-Y. Lu, N.-L. Liu, J. Zhang, and J.-W. Pan, Phys. Rev. A 97, 033851 (2018a).
  * Dixon  _et al._ [2009] P. B. Dixon, D. J. Starling, A. N. Jordan, and J. C. Howell, Phys. Rev. Lett. 102, 173601 (2009).
  * Viza  _et al._ [2013] G. I. Viza, J. Martínez-Rincón, G. A. Howland, H. Frostig, I. Shomroni, B. Dayan, and J. C. Howell, Optics letters 38 16, 2949 (2013).
  * Xu  _et al._ [2013b] X.-Y. Xu, Y. Kedem, K. Sun, L. Vaidman, C.-F. Li, and G.-C. Guo, Phys. Rev. Lett. 111, 033604 (2013b).
  * Brunner and Simon [2010] N. Brunner and C. Simon, Phys. Rev. Lett. 105, 010405 (2010).
  * Qiu  _et al._ [2017] X. Qiu, L. Xie, X. Liu, L. Luo, Z. Li, Z. Zhang, and J. Du, Applied Physics Letters 110, 10.1063/1.4976312 (2017), 071105.
  * Li  _et al._ [2018b] H. Li, J.-Z. Huang, Y. Yu, Y. Li, C. Fang, and G. Zeng, Applied Physics Letters 112, 10.1063/1.5027117 (2018b), 231901.
  * Pfeifer and Fischer [2011] M. Pfeifer and P. Fischer, Weak value amplified optical activity measurements, Opt. Express 19, 16508 (2011).
  * Magaña Loaiza  _et al._ [2014] O. S. Magaña Loaiza, M. Mirhosseini, B. Rodenburg, and R. W. Boyd, Amplification of angular rotations using weak measurements, Phys. Rev. Lett. 112, 200401 (2014).
  * de Lima Bernardo  _et al._ [2014] B. de Lima Bernardo, S. Azevedo, and A. Rosas, Ultrasmall polarization rotation measurements via weak value amplification, Physics Letters A 378, 2029 (2014).
  * Qu  _et al._ [2020] W. Qu, S. Jin, J. Sun, L. Jiang, J. Wen, and Y. Xiao, Nature Communications 11 (2020).
  * Lyons  _et al._ [2015] K. Lyons, J. Dressel, A. N. Jordan, J. C. Howell, and P. G. Kwiat, Phys. Rev. Lett. 114, 170801 (2015).
  * Wang  _et al._ [2016] Y.-T. Wang _et al._ , Phys. Rev. Lett. 117, 230801 (2016).
  * Dressel  _et al._ [2013] J. Dressel, K. Lyons, A. N. Jordan, T. M. Graham, and P. G. Kwiat, Phys. Rev. A 88, 023821 (2013).
  * Krafczyk  _et al._ [2021] C. Krafczyk, A. N. Jordan, M. E. Goggin, and P. G. Kwiat, Phys. Rev. Lett. 126, 220801 (2021).
  * Fang  _et al._ [2020] S.-Z. Fang, L.-L. Zhu, R.-B. Jin, H.-T. Tan, G.-X. Li, and Q.-L. Wu, Optics Communications 460, 125117 (2020).
  * Meers and Strain [1991] B. J. Meers and K. A. Strain, Phys. Rev. D 43, 3117 (1991).
  * Strain and Meers [1991] K. A. Strain and B. J. Meers, Phys. Rev. Lett. 66, 1391 (1991).
  * Abbott  _et al._ [2016a] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016a).
  * Abbott  _et al._ [2016b] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 131103 (2016b).
  * Zhong  _et al._ [2023] Z.-R. Zhong, W.-J. Tan, Y. Chen, and Q.-L. Wu, Dual-recycled interference-based weak value metrology, Phys. Rev. A 108, 032608 (2023).
  * Fang  _et al._ [2021] S.-Z. Fang, H.-T. Tan, G.-X. Li, and Q.-L. Wu, Weak value amplification for angular velocity measurements, Appl. Opt. 60, 4335 (2021).
  * Care [1983] C. M. Care, Probabilistic and statistical aspects of quantum theory: North-holland series in statistics and probability vol 1, Physics Bulletin 34, 395 (1983).
  * Pang  _et al._ [2014] S. Pang, J. Dressel, and T. A. Brun, Entanglement-assisted weak value amplification, Phys. Rev. Lett. 113, 030401 (2014).
  * Chen  _et al._ [2019] J.-S. Chen, B.-H. Liu, M.-J. Hu, X.-M. Hu, C.-F. Li, G.-C. Guo, and Y.-S. Zhang, Realization of entanglement-assisted weak-value amplification in a photonic system, Phys. Rev. A 99, 032120 (2019).
  * Helstrom and Kennedy [1974] C. Helstrom and R. Kennedy, Noncommuting observables in quantum detection and estimation theory, IEEE Transactions on Information Theory 20, 16 (1974).