# Consistency between the attached-eddy model and the inner-outer interaction
model: a study of streamwise wall-shear stress fluctuations in a turbulent
channel flow

Cheng Cheng1  Lin Fu1,2,3,4<EMAIL_ADDRESS>1Department of Mechanical and
Aerospace Engineering, The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong 2Department of Mathematics, The Hong Kong
University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong 3
Center for Ocean Research in Hong Kong and Macau (CORE), The Hong Kong
University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
4Shenzhen Research Institute, The Hong Kong University of Science and
Technology, Shenzhen, China

###### Abstract

The inner-outer interaction model (Marusic, Mathis & Hutchins, Science, vol.
329, 2010, 193-196) and the attached-eddy model (Townsend, Cambridge
University Press, 1976) are two fundamental models describing the multi-scale
turbulence interactions and the organization of energy-containing motions in
the logarithmic region of high-Reynolds number wall-bounded turbulence,
respectively. In this paper, by coupling the additive description with the
attached-eddy model, the generation process of streamwise wall-shear
fluctuations, resulting from wall-attached eddies, is portrayed. Then, by
resorting to the inner-outer interaction model, the streamwise wall-shear
stress fluctuations generated by attached eddies in a turbulent channel flow
are isolated. Direct comparison between the statistics from these two models
demonstrates that they are consistent to and complement each other. Meanwhile,
we further show that the superpositions of attached eddies follow an additive
process strictly by verifying the validity of the strong and extended self
similarity. Moreover, we propose a Gaussian model to characterize the
instantaneous distribution of streamwise wall-shear stress, resulting from the
attached-eddy superpositions. These findings are important for developing an
advanced reduced-order wall model.

###### keywords:

MSC Codes  (Optional) Please enter your MSC Codes here

##  1 Introduction

Wall-shear stress fluctuation is a crucial physical quantity in wall-bounded
turbulence, as it is of importance for noise radiation, structural vibration,
drag generation, and wall heat transfer, among others (Diaz-Daniel et al.,
2017; Cheng et al., 2020). In the past two decades, ample evidence has shown
that the root mean squared value of streamwise wall-shear stress fluctuations
($\tau_{x,rms}^{{}^{\prime}}$) is sensitive to the flow Reynolds number (Abe
et al., 2004; Schlatter & Örlü, 2010; Yang & Lozano-Durán, 2017; Guerrero et
al., 2020). It indicates that large-scale energy-containing eddies populating
the logarithmic and outer regions in high-Reynolds-number wall turbulence have
non-negligible influences on the near-wall turbulence dynamics, and thus the
wall friction (de Giovanetti et al., 2016; Li et al., 2019).

Till now, several models have been proposed on the organization of motions in
logarithmic and outer regions and their interactions with the near-wall
dynamics. Marusic et al. (2010) have established that superposition and
modulation are the two basic mechanisms that large-scale motions (LSM) and
very-large-scale motions (VLSM) exert influences on the near-wall turbulence.
The former refers to the footprints of LSMs and VLSMs on the near-wall
turbulence, while the latter indicates the intensity amplification or
attenuation of near-wall small-scale turbulence by the outer motions. Mathis
et al. (2013) extended the model to interpret the generation of wall-shear
stress fluctuations in high-Reynolds number flows. They emphasized that
superposition and modulation are still two essential factors. This inner-outer
interaction model (IOIM) has also been successfully developed to predict the
near-wall velocity fluctuations with data inputs from the log layer (Marusic
et al., 2010; Baars et al., 2016; Wang et al., 2021).

On the other hand, the most elegant conceptual model describing the motions in
logarithmic region is the attached-eddy model (AEM) (Townsend, 1976; Perry &
Chong, 1982). It conjectures that the logarithmic region is occupied by an
array of self-similar energy-containing motions (or eddies) with their roots
attached to the near-wall region. Extensive validations support the existence
of attached eddies in high-Reynolds number turbulence, such as the logarithmic
decaying of streamwise velocity fluctuation intensities (Meneveau & Marusic,
2013), as originally predicted by Townsend (1976). The reader is referred to a
recent review work by Marusic & Monty (2019) for more details. Given the
existence of wall-attached energy-containing motions in the logarithmic
region, it would be quite natural to hypothesize that the near-wall part of
these motions would affect the generation of the wall-shear fluctuations to
some extent, maybe, via the superposition and modulation mechanisms. However,
some fundamental questions may be raised, e.g., whether the IOIM and AEM are
consistent with each other? There’s a possibility that the superposition
component of $\tau_{x}^{\prime}$ decomposed by the IOIM in physical space can
not fully follow the predictions made by the AEM quantitatively. If yes,
whether these two models can shed light on the mechanism of wall-shear
fluctuation generation and be indicative for modeling approaches?

Previous study (Yang & Lozano-Durán, 2017) verified that the generation of
wall-shear stress fluctuations can be interpreted as the outcomes of the
momentum cascade across momentum-carried eddies of different scales, and
modeled by an additive process. Here, we first aim to couple the additive
description with the AEM to portray the generation process of streamwise wall-
shear fluctuations, resulting from wall-attached eddies. Two scaling laws
describing their intensities and the linkages with the characteristic scales
of attached eddies can be derived (the characteristic scales of attached
eddies are their wall-normal heights according to AEM (Townsend, 1976)). Then,
we intend to isolate the streamwise wall-shear stress fluctuations generated
by attached eddies in a turbulent channel flow at $Re_{\tau}=2003$
($Re_{\tau}=hu_{\tau}/\nu$, $h$ denotes the channel half-height, $u_{\tau}$
the wall friction velocity and $\nu$ the kinematic viscosity) by resorting to
the IOIM (Marusic et al., 2010; Baars et al., 2016). Here, IOIM is employed as
a tool to estimate the streamwise wall-shear fluctuations generated by
attached eddies. The statistics from IOIM can be processed to verify the
scaling laws deduced by AEM, so as to demonstrate their consistency. Moreover,
a simple algebraic model describing the instantaneous distributions of the
streamwise wall-shear stress fluctuations generated by attached eddies will be
proposed.

##  2 Streamwise wall-shear stress fluctuations generated by attached eddies

According to Mandelbrot (1974) and Yang & Lozano-Durán (2017), the generation
of streamwise wall-stress fluctuations can be modeled as an additive process
within multifractal formalism, which takes the form of

$\tau_{x}^{{}^{\prime}+}=\sum_{i=1}^{n}a_{i},$ (1)

where $a_{i}$ are random addends, representing an increment in
$\tau_{x}^{{}^{\prime}+}$ due to eddies with wall-normal height $h/2^{i}$, and
superscript $+$ denotes the normalization with wall units. Here, we intend to
isolate the contributions from the eddies populating logarithmic region
($\tau_{x,o}^{{}^{\prime}+}$) and link to their wall-normal positions $y$.
$\tau_{x,o}^{{}^{\prime}+}$ can be expressed as

$\tau_{x,o}^{{}^{\prime}+}=\sum_{i=n_{s}}^{n_{o}}a_{i},$ (2)

where $n_{s}$ and $n_{o}$ represent the additives that correspond to the
eddies with the wall-normal height at $y_{s}$ and $y_{o}$, respectively. Here,
$y_{s}$ is the lower bound of logarithmic region, and generally believed to be
$80\leq y_{s}^{+}\leq 100$ (Jiménez, 2018; Baars & Marusic, 2020); $y_{o}$ is
the outer reference height. It can be found that $1\textless n_{s}\textless
n_{o}\textless n$. The addends $a_{i}$ are assumed to be identically and
independently distributed (i.i.d) and equal to $a$. The number of the addends
should be proportional to

$n_{o}-n_{s}+1\sim\int_{y_{s}}^{y_{o}}p(y)dy\sim\int_{y_{s}}^{y_{o}}\frac{1}{y}dy\sim\ln(\frac{y_{o}}{y_{s}}),$
(3)

where $p(y)$ is the eddy population density, which is proportional to $1/y$
according to AEM (Townsend, 1976; Perry & Chong, 1982). A momentum generation
function $\left\langle exp(q\tau_{x,o}^{{}^{\prime}+})\right\rangle$, where
$<>$ represents the averaging in the temporal and spatially homogeneous
directions, is defined to scrutinize the scaling behavior of
$\tau_{x,o}^{{}^{\prime}+}$ (Yang et al., 2016). $\left\langle
exp(q\tau_{x,o}^{{}^{\prime}+})\right\rangle$ can be evaluated as

$\left\langle exp(q\tau_{x,o}^{{}^{\prime}+})\right\rangle=\left\langle
exp(qa)\right\rangle^{n_{o}-n_{s}+1}\sim(\frac{y_{o}}{y_{s}})^{s(q)},$ (4)

where $q$ is a real number, $s(q)=C_{1}\ln\left\langle exp(qa)\right\rangle$
is called anomalous exponent, $C_{1}$ is a constant. Eq. (4) is called strong
self similarity (SSS). If $a$ is a Gaussian variable, the anomalous exponent
can be recast as

$s(q)=C_{2}q^{2},$ (5)

where $C_{2}$ is another constant. On the other hand, an extended self-
similarity (ESS) is defined to describe the relationship between $\left\langle
exp(q\tau_{x,o}^{{}^{\prime}+})\right\rangle$ and $\left\langle
exp(q_{0}\tau_{x,o}^{{}^{\prime}+})\right\rangle$ (fixed $q_{0}$) (Benzi et
al., 1993), i.e.,

$\displaystyle\left\langle
exp(q\tau_{x,o}^{{}^{\prime}+})\right\rangle=\left\langle
exp(q_{0}\tau_{x,o}^{{}^{\prime}+})\right\rangle^{\xi(q,q_{0})},$ (6)

where $\xi(q,q_{0})$ is a function of $q$ (fixed $q_{0}$). Note that ESS does
not strictly rely on i.i.d of the addends, but the additive process Eq. (2).

##  3 DNS database and scale decomposition method

The direct numerical simulation (DNS) database used in the present study is an
incompressible turbulent channel flow at $Re_{\tau}=2003$, which has been
extensively validated by previous studies (Hoyas & Jiménez, 2006; Jiménez &
Hoyas, 2008). The decomposition of $\tau^{\prime}_{x}$ is based on the IOIM
first proposed by Marusic et al. (2010). Baars et al. (2016) modified the
computational process by introducing spectral stochastic estimation to avoid
artificial scale decomposition. In this work, the modified version of IOIM is
adopted to investigate the multi-scale characteristics of $\tau^{\prime}_{x}$.
It can be expressed as

$u_{p}^{+}\left(y^{+}\right)=\underbrace{u^{*}\left(y^{+}\right)\left\\{1+\Gamma_{uu}u_{L}^{+}\left(y^{+}\right)\right\\}}_{u_{s}^{+}}+u_{L}^{+}\left(y^{+}\right),$
(7)

where $u_{p}^{+}$ denotes the predicted near-wall streamwise velocity
fluctuation, $u^{*}$ denotes the universal velocity signal without large-scale
impact, $u_{L}^{+}$ is the superposition component, $\Gamma_{uu}$ is the
amplitude-modulation coefficient, and $u_{s}^{+}$ denotes the amplitude
modulation of the universal signal $u^{*}$. $u_{L}^{+}$ is obtained by
spectral stochastic estimation of the streamwise velocity fluctuation at the
logarithmic region $y_{o}^{+}$, namely,

$u_{L}^{+}\left(x^{+},y^{+},z^{+}\right)=F_{x}^{-1}\left\\{H_{L}\left(\lambda_{x}^{+},y^{+}\right)F_{x}\left[u_{o}^{+}\left(x^{+},y_{o}^{+},z^{+}\right)\right]\right\\}$
(8)

where $u_{o}^{+}$ is the streamwise velocity fluctuation at $y_{o}^{+}$ in the
logarithmic region, and, $F_{x}$ and $F_{x}^{-1}$ denote FFT and inverse FFT
in the streamwise direction, respectively. $H_{L}$ is the transfer kernel,
which evaluates the correlation between $u^{+}(y^{+})$ and
$u_{o}^{+}(y_{o}^{+})$ at a given length scale $\lambda_{x}^{+}$, and can be
calculated as

$H_{L}\left(\lambda_{x}^{+},y^{+}\right)=\frac{\left\langle\hat{u}\left(\lambda_{x}^{+},y^{+},z^{+}\right)\overline{\hat{u}_{o}}\left(\lambda_{x}^{+},y_{o}^{+},z^{+}\right)\right\rangle}{\left\langle\hat{u}_{o}\left(\lambda_{x}^{+},y_{o}^{+},z^{+}\right)\overline{\hat{u}_{o}}\left(\lambda_{x}^{+},y_{o}^{+},z^{+}\right)\right\rangle},$
(9)

where $\hat{u}$ is the Fourier coefficient of $u$, and $\overline{\hat{u}}$ is
the complex conjugate of $\hat{u}$.

Figure 1: A schematic of the attached-eddy model (Hwang, 2015). Each circle
represents an individual attached eddy. $y_{s}^{+}$ and $y_{e}^{+}$ are the
lower and upper bound of the logarithmic region, respectively. $y_{0}^{+}$ is
the outer reference height, and varies from $y_{s}^{+}$ to $y_{e}^{+}$.

In this work, we mainly pay attention to the quantity, $\tau^{\prime}_{x}$
generated by the attached eddies. Thus, the predicted position $y^{+}$ is
fixed at $y^{+}=0.3$, and the outer reference height $y_{o}^{+}$ varies from
$100$ (namely $y_{s}^{+}$) to $0.2h^{+}$ (denoted as $y_{e}^{+}$), i.e., the
upper boundary of logarithmic region (Jiménez, 2018). We have checked that as
long as the predicted position is around $y^{+}\leq O(1)$, the results
presented below are insensitive to the choice of specific $y^{+}$. Once
$u_{L}^{+}$ is obtained, the superposition component of
$\tau_{x}^{{}^{\prime}+}$ can be calculated by definition (i.e.,
$\frac{\partial u_{L}^{{}^{\prime}+}}{\partial y^{+}}$ at the wall) and
denoted as $\tau_{x,L}^{{}^{\prime}+}(y_{o}^{+})$. According to the
hierarchical attached eddies in high-Reynolds number wall turbulence (see Fig.
1), $\tau_{x,L}^{{}^{\prime}+}(y_{o}^{+})$ represents the superposition
contributed from the wall-coherent motions with their height larger than
$y_{o}^{+}$. Thus, the difference value
$\tau_{x,L}^{{}^{\prime}+}(y_{s}^{+})-\tau_{x,L}^{{}^{\prime}+}(y_{o}^{+})$
can be interpreted as the superposition contribution generated by the wall-
coherent eddies with their wall-normal heights within $y_{s}^{+}$ and
$y_{o}^{+}$, i.e., $\tau_{x,o}^{{}^{\prime}+}$ in Eq. (2). Considering that
$y_{s}^{+}$ is the lower bound of the logarithmic region, the increase of
$y_{o}^{+}$ corresponds to the enlargement of the addends in the additive
description (see Eq. (2)). In this way, the connection between AEM and IOIM
are established, and the AEM predictions (see Eqs. (4)-(6)) can be verified
directly.

##  4 Results and discussion

###  4.1 Scaling laws of $\tau_{x,o}^{{}^{\prime}+}$

Figure 2: ($a$) $G$ as functions of $y_{o}/y_{s}$ for $q=\pm 5$ and $q=\pm 3$;
($b$) premultiplied $G$ as functions of $y_{o}/y_{s}$ for $q=\pm 5$ and $q=\pm
3$.

Here, we further define a moment generation function based on the IOIM, i.e.,

$G(q,y_{o}^{+})=\left\langle
exp(q(\tau_{x,L}^{{}^{\prime}+}(y_{s}^{+})-\tau_{x,L}^{{}^{\prime}+}(y_{o}^{+})))\right\rangle.$
(10)

Fig. 2($a$) shows the variations of $G$ as a function of $y_{o}/y_{s}$ for
$q=\pm 5$ and $q=\pm 3$. Power-law behaviours can be found in the interval
between $1.7\leq y_{o}/y_{s}\leq 2.9$ for positive $q$ and $1.7\leq
y_{o}/y_{s}\leq 4$ for negative $q$, justifying the validity of SSS, i.e., Eq.
(4). Fig. 2 is in aid of accessing the scalings by displaying the variations
of premultiplied $G$. This observation highlights that the superpositions of
wall-attached log-region motions on wall surface follow the additive process,
characterized by Eq. (2). It is also worth mentioning that the power-law
behaviour can be observed for larger wall-normal intervals for negative $q$.
As $G(q,y_{o}^{+})$ quantifies
$\tau_{x,L}^{{}^{\prime}+}(y_{s}^{+})-\tau_{x,L}^{{}^{\prime}+}(y_{o}^{+})$,
which features the same sign as $q$, this observation is consistent with the
work of Cheng et al. (2020), which showed that the footprints of the inactive
part of attached eddies populating the logarithmic region are actively
connected with large-scale negative $\tau^{\prime}_{x}$. Other $q$ values
yield similar results and are not shown here for brevity.

Figure 3: ($a$) Anomalous exponent $s(q)$ as a function of $q$. The black line
is a quadratic fit; ($b$) second- to sixth- order moments of
$\tau_{x}^{{}^{\prime}+}$ as functions of $Re_{\tau}$. The dashed lines are
the log-normal predictions from Eq. (11)-(13).

The anomalous exponent $s(q)$ can be obtained by fitting the range $2\leq
y_{o}/y_{s}\leq 2.9$, where both positive and negative $q$ display good power-
law scalings. Fig. 3($a$) displays the variation of the anomalous exponent
$s(q)$ as a function of $q$. The solid line denotes the quadratic fit within
$-0.5\leq q\leq 0.5$. It can be seen that the variation of $s(q)$ is very
close to the model prediction, i.e., the quadratic function as Eq. (5) with
$C_{2}=0.00629$. Only minor discrepancies between DNS data and model
predictions can be observed. As such, it is reasonable to hypothesize that the
streamwise wall-shear stress fluctuation $\tau^{\prime}_{x}$ generated by
attached eddies of a given size follows the Gaussian distribution. Moreover,
we can also estimate the statistical moments of $\tau_{x,o}^{+}$ by taking the
derivative of $G(q,y_{e}^{+})$ with respect to $q$ around $q=0$ (Yang et al.,
2016), i.e.,

$\left\langle\tau_{x,o}^{{}^{\prime}2+}\right\rangle=\left.\frac{\partial^{2}G(q;y_{o}^{+})}{\partial
q^{2}}\right|_{q=0}\sim 2C_{2}\ln(y_{o}/y_{s})\sim 2C_{2}\ln Re_{\tau},$ (11)
$\left\langle\tau_{x,o}^{{}^{\prime}4+}\right\rangle^{1/2}=(\left.\frac{\partial^{4}G(q;y_{o}^{+})}{\partial
q^{4}}\right|_{q=0})^{1/2}\sim 2\sqrt{3}C_{2}\ln(y_{o}/y_{s})\sim
2\sqrt{3}C_{2}\ln Re_{\tau},$ (12)
$\left\langle\tau_{x,o}^{{}^{\prime}6+}\right\rangle^{1/3}=(\left.\frac{\partial^{6}G(q;y_{o}^{+})}{\partial
q^{6}}\right|_{q=0})^{1/3}\sim 2\sqrt[3]{15}C_{2}\ln(y_{o}/y_{s})\sim
2\sqrt[3]{15}C_{2}\ln Re_{\tau}.$ (13)

Fig. 3($b$) shows the variations of second- ($p=1$) to sixth- ($p=3$) order
moments of $\tau^{\prime}_{x}$ calculated from DNS of channel flows (Iwamoto
et al., 2002; Del Álamo & Jiménez, 2003; Abe et al., 2004; Del Álamo et al.,
2004; Hu et al., 2006; Lozano-Durán & Jiménez, 2014; Lee & Moser, 2015; Cheng
et al., 2019; Kaneda & Yamamoto, 2021) and compares them with the model
prediction, i.e., Eq. (11)-(13). For the second- and fourth- order variances,
the model predictions are roughly consistent with the DNS results. The
comparisons also indicate a Reynolds-number dependence of
$\left\langle\tau_{x}^{{}^{\prime}2+}\right\rangle$, which has been reported
by vast studies (Schlatter & Örlü, 2010; Mathis et al., 2013; Guerrero et al.,
2020), and may be ascribed to the superposition effects of the wall-attached
log-region motions. Wang et al. (2020) speculated that the amplitude
modulation effect plays a more prominent role in affecting the statistic
characteristics of $\tau_{x,rms}^{{}^{\prime}+}$ than the superposition
effect, which contradicts the present findings. In fact, amplitude modulation
has been demonstrated to exert a negligible effect on the even-order moments
(Mathis et al., 2011; Blackman et al., 2019). Therefore, the deduction of Wang
et al. (2020) needs to be revisited. For sixth-order moments, the model
prediction displays substantial discrepancies with the DNS data. It is
expected since high-order moments are dominated by the rare events resulting
from the intermittent small-scale motions (Frisch & Donnelly, 1996), which can
not be captured by IOIM (see Fig. 5($a$)).

Figure 4: ($a$) $G(q)$ as functions of $G(-2)$ for $q=-1,-3,-5,-7$; ($b$)
$G(q)$ as functions of $G(2)$ for $q=1,3,5,7$. Both vertical and horizontal
axes in ($a$) and ($b$) are plotted in logarithmic form.

ESS (i.e., Eq. (6)) is another scaling predicted by the multifractal
formalism. Different from SSS, ESS does not rely on the i.i.d of the addends,
but the additive process (see Eq. (2)). Fig. 4($a$) and 4($b$) shows the ESS
scalings for $q_{0}=-2$ and $q_{0}=2$, respectively. ESS holds for the entire
logarithmic region. The observation suggests that the streamwise wall-shear
fluctuations generated by logarithmic motions obey the additive process,
though the streamwise wall shear fluctuations generated by attached eddies
with wall-normal heights at approximately $0.2h^{+}$ are not identically and
independently distributed due to the scale interactions (see Fig. 2), which
are not described by the attached-eddy model.

###  4.2 Instantaneous distribution of $\tau^{\prime}_{x}$

Furthermore, the instantaneous $\tau_{x}^{{}^{\prime}+}$ can be decomposed as

$\tau_{x}^{{}^{\prime}+}=\tau_{x,s}^{{}^{\prime}+}+\underbrace{\tau_{x,L}^{{}^{\prime}+}(y_{s}^{+})-\tau_{x,L}^{{}^{\prime}+}(y_{e}^{+})}_{\tau_{x,log}^{{}^{\prime}+}}+\underbrace{\tau_{x,L}^{{}^{\prime}+}(y_{e}^{+})}_{\tau_{x,out}^{{}^{\prime}+}},$
(14)

where $\tau_{x,s}^{{}^{\prime}+}$ denotes the amplitude modulation of the
universal signal $\tau_{x}^{{}^{\prime}*+}$, $\tau_{x,log}^{{}^{\prime}+}$ and
$\tau_{x,out}^{{}^{\prime}+}$ are the superposition components contributed
from the log region and the outer wall-coherent motions, respectively. The
methodology of removing modulation effects can be found in Mathis et al.
(2011) and Baars et al. (2016), whose details are out of the range of present
study. Fig. 5($a$) shows the probability density functions (p.d.f.s) of
$\tau_{x}^{{}^{\prime}*+}$, $\tau_{x,s}^{{}^{\prime}+}$,
$\tau_{x,log}^{{}^{\prime}+}$ and $\tau_{x,out}^{{}^{\prime}+}$, and compares
with the p.d.f for the full channel data. The p.d.f.s of
$\tau_{x,s}^{{}^{\prime}+}$ and $\tau_{x}^{{}^{\prime}*+}$ nearly coincide
with that of $\tau_{x}^{{}^{\prime}+}$ with asymmetric and positively skewed
shape, which demonstrates that removing the superposition and modulation
effects barely affects the instantaneous distributions. The asymmetries
between the positive and negative wall-shear fluctuations are the essential
characters of the near-wall small-scale turbulence, which may be associated
with the celebrated near-wall sustaining process (Schoppa & Hussain, 2002). In
contrast, the p.d.f.s of $\tau_{x,log}^{{}^{\prime}+}$ and
$\tau_{x,out}^{{}^{\prime}+}$ are more symmetric with rare events invisible,
suggesting that the superposition components of logarithmic and outer motions
are less intermittent than the small-scale universal signals. This also
explains the reason why the log-normal model describes the additive process
well (see Fig. 3($a$)), although the log-normal model is inapplicable for rare
events (Landau & Lifshitz, 1987). Moreover, the skewness and flatness of
$\tau_{x,log}^{{}^{\prime}+}$ are 0.05 and 2.91, which are very close to those
of a Gaussian distribution. It strongly supports the conclusion drawn above
that the streamwise wall-shear stress fluctuations generated by attached
eddies populating logarithmic region can be absolutely treated as Gaussian
variables with

$p(\xi)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\xi^{2}}{2\sigma^{2}}\right),$
(15)

where $p(\xi)$ denotes the p.d.f, and $\xi$ is the independent variable. It is
worth noting that the variation of variance can be well predicted by the log-
normal model, namely

$\sigma^{2}=\left.\frac{\partial^{2}G(q;y_{o}^{+})}{\partial
q^{2}}\right|_{q=0}=2C_{2}\ln(Re_{\tau})+C_{3},$ (16)

where $C_{2}\approx 0.00629$, and $C_{3}\approx-0.07959$ is a constant and
determined by the DNS data at $Re_{\tau}=2003$. Fig. 5($b$) shows the p.d.f.s
of $\tau_{x,log}^{{}^{\prime}+}$ and the model prediction by Eq. (15), results
of other two Reynolds numbers (Del Álamo et al., 2004; Lozano-Durán & Jiménez,
2014) are also included for comparison. It can be seen that the Gaussian model
proposed here works reasonably well and can cover a wide range of Reynolds
numbers. The model remains to be validated by higher-Reynolds number DNS data.

Figure 5: ($a$) P.d.f.s of $\tau_{x}^{{}^{\prime}*+}$,
$\tau_{x,s}^{{}^{\prime}+}$, $\tau_{x,log}^{{}^{\prime}+}$,
$\tau_{x,out}^{{}^{\prime}+}$, and $\tau_{x}^{{}^{\prime}+}$; ($b$) P.d.f.s of
$\tau_{x,log}^{{}^{\prime}+}$ in channel flows with $Re_{\tau}=934$, $2003$,
and $4179$. Dashed lines denote the Gaussian model predictions with Eqs.
(15)-(16).

##  5 Concluding remarks

In summary, the present study reveals that IOIM and AEM are consistent to each
other quantitatively. The statistical characteristics of the superpositions of
log-region eddies follow the predictions of AEM, namely, the SSS and ESS
scalings. Based on these observations, we conclude that the streamwise wall-
shear stress fluctuations generated by attached eddies populating the
logarithmic region can be treated as Gaussian variables. A Gaussian model is
then proposed to describe their instantaneous distributions and verified by
DNS data spanned broad-band Reynolds numbers. Considering the fact that the
intensity of wall-shear stress fluctuations is typically underpredicted by the
state-of-the-art wall-modelled large-eddy simulation (WMLES) approaches (Park
& Moin, 2016), the Gaussian model proposed in the present study may be
constructive for the development of the LES methodology, and the distribution
characteristics of $\tau_{x}^{{}^{\prime}*+}$ are helpful for developing more
accurate near-wall models of WMLES approaches.

It is noted that some previous works adopted IOIM to investigate the spectral
characteristics of the wall-coherent components of the signals in near-wall
region, such as the work of Marusic et al. (2017), but whether they are
consistent with the AEM predictions in physical space quantitatively has not
been verified in detail. The consistency of the two models demonstrated here
fills the gap and complements their works. Moreover, the findings in present
study indicate that we can isolate the footprints of attached eddies within a
selected wall-normal range by employing IOIM, i.e., by adjusting $y_{1}^{+}$
and $y_{2}^{+}$ in
$\tau_{x,L}^{{}^{\prime}+}(y_{1}^{+})-\tau_{x,L}^{{}^{\prime}+}(y_{2}^{+})$.
Here $y_{1}^{+}$ and $y_{2}^{+}$ are two selected wall-normal heights in the
logarithmic region, and $y_{1}^{+}<y_{2}^{+}$. In this regard, the present
study may provide a new perspective for analyzing some flow physics in wall-
bounded turbulence, such as the inner peak of the intensity of $u^{\prime}$,
and the streamwise inclined angles of attached eddies. All these are under
investigation currently and will be reported in separate forthcoming papers.

## Acknowledgments

We are grateful to the authors cited in Fig. 3($b$) for making their
invaluable data available. L.F. acknowledges the fund from CORE as a joint
research center for ocean research between QNLM and HKUST.

## Declaration of interests

The authors report no conflict of interest.

## References

  * Abe et al. (2004) Abe, H., Kawamura, H. & Choi, H. 2004 Very large-scale structures and their effects on the wall shear-stress fluctuations in a turbulent channel flow up to $\mathit{Re}_{{\it\tau}}=640$. J. Fluids Eng. 126 (5), 835–843.
  * Baars & Marusic (2020) Baars, W.J. & Marusic, I. 2020 Data-driven decomposition of the streamwise turbulence kinetic energy in boundary layers. part 1. energy spectra. J. Fluid Mech. 882.
  * Baars et al. (2016) Baars, W. J, Hutchins, N. & Marusic, I. 2016 Spectral stochastic estimation of high-reynolds-number wall-bounded turbulence for a refined inner-outer interaction model. Physical Review Fluids 1 (5), 054406.
  * Benzi et al. (1993) Benzi, R., Ciliberto, S., Tripiccione, R., Baudet, C., Massaioli, F. & Succi, S. 1993 Extended self-similarity in turbulent flows. Phys. Rev. E 48, R29–R32.
  * Blackman et al. (2019) Blackman, K., Perret, L. & Mathis, R. 2019 Assessment of inner–outer interactions in the urban boundary layer using a predictive model. J. Fluid Mech. 875, 44–70.
  * Cheng et al. (2019) Cheng, C., Li, W., Lozano-Durán, A. & Liu, H. 2019 Identity of attached eddies in turbulent channel flows with bidimensional empirical mode decomposition. J. Fluid Mech 870, 1037–1071.
  * Cheng et al. (2020) Cheng, C., Li, W., Lozano-Durán, A. & Liu, H. 2020 On the structure of streamwise wall-shear stress fluctuations in turbulent channel flows. J. Fluid Mech. 903, A29.
  * Del Álamo & Jiménez (2003) Del Álamo, J. C. & Jiménez, J. 2003 Spectra of the very large anisotropic scales in turbulent channels. Phys. Fluids 15 (6), L41–L44.
  * Del Álamo et al. (2004) Del Álamo, J. C., Jiménez, J., Zandonade, P. & Moser, R.D. 2004 Scaling of the energy spectra of turbulent channels. J. Fluid Mech. 500, 135–144.
  * Diaz-Daniel et al. (2017) Diaz-Daniel, C., Laizet, S. & Vassilicos, J.C. 2017 Wall shear stress fluctuations: Mixed scaling and their effects on velocity fluctuations in a turbulent boundary layer. Phys. Fluids 29 (5), 055102.
  * Frisch & Donnelly (1996) Frisch, U. & Donnelly, R.J. 1996 Turbulence: the legacy of AN Kolmogorov. AIP.
  * de Giovanetti et al. (2016) de Giovanetti, M., Hwang, Y. & Choi, H. 2016 Skin-friction generation by attached eddies in turbulent channel flow. J. Fluid Mech. 808, 511–538.
  * Guerrero et al. (2020) Guerrero, B., Lambert, M.F. & Chin, R.C. 2020 Extreme wall shear stress events in turbulent pipe flows: spatial characteristics of coherent motions. J. Fluid Mech. 904, A18.
  * Hoyas & Jiménez (2006) Hoyas, S. & Jiménez, J. 2006 Scaling of the velocity fluctuations in turbulent channels up to $\mathit{Re}_{{\it\tau}}=2003$. Phys. Fluids 18 (1), 011702.
  * Hu et al. (2006) Hu, Z.W., Morfey, C.L. & Sandham, N.D. 2006 Wall pressure and shear stress spectra from direct simulations of channel flow. AIAA Journal 44 (7), 1541–1549.
  * Hwang (2015) Hwang, Y. 2015 Statistical structure of self-sustaining attached eddies in turbulent channel flow. J. Fluid Mech. 767, 254–289.
  * Iwamoto et al. (2002) Iwamoto, K., Suzuki, Y. & Kasagi, N. 2002 Reynolds number effect on wall turbulence: toward effective feedback control. Int. J. Heat Fluid Flow 23 (5), 678–689.
  * Jiménez (2018) Jiménez, J. 2018 Coherent structures in wall-bounded turbulence. J. Fluid Mech. 842, P1.
  * Jiménez & Hoyas (2008) Jiménez, J. & Hoyas, S. 2008 Turbulent fluctuations above the buffer layer of wall-bounded flows. J. Fluid Mech. 611, 215–236.
  * Kaneda & Yamamoto (2021) Kaneda, Y. & Yamamoto, Y. 2021 Velocity gradient statistics in turbulent shear flow: an extension of kolmogorov’s local equilibrium theory. J. Fluid Mech. 929, A13.
  * Landau & Lifshitz (1987) Landau, L. & Lifshitz, E 1987 Fluid mechanics. In Course of Theoretical Physics, Pergamon.
  * Lee & Moser (2015) Lee, M. & Moser, R.D. 2015 Direct numerical simulation of turbulent channel flow up to $\mathit{Re}_{{\it\tau}}\approx 5200$. J. Fluid Mech. 774, 395–415.
  * Li et al. (2019) Li, W., Fan, Y., Modesti, D. & Cheng, C. 2019 Decomposition of the mean skin-friction drag in compressible turbulent channel flows. J. Fluid Mech. 875, 101–123.
  * Lozano-Durán & Jiménez (2014) Lozano-Durán, A. & Jiménez, J. 2014 Time-resolved evolution of coherent structures in turbulent channels: characterization of eddies and cascades. J. Fluid Mech. 759, 432–471.
  * Mandelbrot (1974) Mandelbrot, B.B 1974 Intermittent turbulence in self-similar cascades: divergence of high moments and dimension of the carrier. J. Fluid Mech. 62 (2), 331–358.
  * Marusic et al. (2017) Marusic, I., Baars, W.J. & Hutchins, N. 2017 Scaling of the streamwise turbulence intensity in the context of inner-outer interactions in wall turbulence. Phys. Rev. Fluids 2 (10), 100502.
  * Marusic et al. (2010) Marusic, I., Mathis, R. & Hutchins, N. 2010 Predictive model for wall-bounded turbulent flow. Science 329 (5988), 193–196.
  * Marusic & Monty (2019) Marusic, I. & Monty, J.P. 2019 Attached eddy model of wall turbulence. Annu. Rev. Fluid Mech. 51, 49–74.
  * Mathis et al. (2011) Mathis, R., Hutchins, N. & Marusic, I. 2011 A predictive inner–outer model for streamwise turbulence statistics in wall-bounded flows. J. Fluid Mech. 681, 537–566.
  * Mathis et al. (2013) Mathis, R., Marusic, I., Chernyshenko, S.I. & Hutchins, N. 2013 Estimating wall-shear-stress fluctuations given an outer region input. J. Fluid Mech. 715, 163–180.
  * Meneveau & Marusic (2013) Meneveau, C. & Marusic, I. 2013 Generalized logarithmic law for high-order moments in turbulent boundary layers. J. Fluid Mech. 719.
  * Park & Moin (2016) Park, G.I. & Moin, P. 2016 Space-time characteristics of wall-pressure and wall shear-stress fluctuations in wall-modeled large eddy simulation. Phys. Rev. Fluids 1 (2), 024404\.
  * Perry & Chong (1982) Perry, A. E. & Chong, M. S 1982 On the mechanism of wall turbulence. J. Fluid Mech. 119 (119), 173–217.
  * Schlatter & Örlü (2010) Schlatter, P. & Örlü, R. 2010 Assessment of direct numerical simulation data of turbulent boundary layers. J. Fluid Mech. 659, 116–126.
  * Schoppa & Hussain (2002) Schoppa, W & Hussain, Fazle 2002 Coherent structure generation in near-wall turbulence. J. Fluid Mech. 453, 57–108.
  * Townsend (1976) Townsend, A. A. 1976 The structure of turbulent shear flow, 2nd edn. Cambridge University Press.
  * Wang et al. (2020) Wang, J., Pan, C. & Wang, J. 2020 Characteristics of fluctuating wall-shear stress in a turbulent boundary layer at low-to-moderate Reynolds number. Phys. Rev. Fluids 5 (7), 074605.
  * Wang et al. (2021) Wang, L., Hu, R. & Zheng, X. 2021 A scaling improved inner–outer decomposition of near-wall turbulent motions. Physics of Fluids 33 (4), 045120.
  * Yang & Lozano-Durán (2017) Yang, X.I.A. & Lozano-Durán, A. 2017 A multifractal model for the momentum transfer process in wall-bounded flows. J. Fluid Mech. 824, R2.
  * Yang et al. (2016) Yang, X.I.A., Marusic, I. & Meneveau, C. 2016 Moment generating functions and scaling laws in the inertial layer of turbulent wall-bounded flows. J Fluid Mech 791, R2.


# Expansions for a fundamental solution of Laplace’s equation on
$\mathbb{R}^{3}$ in 5-cyclidic harmonics

Howard S. Cohl1 and Hans Volkmer2 1Applied and Computational Mathematics
Division, National Institute of Standards and Technology, Gaithersburg, MD
20899-8910, USA  2Department of Mathematical Sciences, University of
Wisconsin–Milwaukee, P. O. Box 413, Milwaukee, WI 53201, USA

###### Abstract.

We derive eigenfunction expansions for a fundamental solution of Laplace’s
equation in three-dimensional Euclidean space in 5-cyclidic coordinates. There
are three such expansions in terms of internal and external 5-cyclidic
harmonics of first, second and third kind. The internal and external
5-cyclidic harmonics are expressed by solutions of a Fuchsian differential
equation with five regular singular points.

##  1\. Introduction

Expansions for a fundamental solution of Laplace’s equation on
$\mathbb{R}^{3}$ in terms of solutions found by the method of separation of
variables in a suitable curvilinear coordinate system are known for a long
time. For example, when we choose spherical coordinates, we obtain the well-
known expansion [21]

(1.1)
$\frac{1}{\|\mathbf{r}-\mathbf{r}^{\prime}\|}=\sum_{\ell=0}^{\infty}\frac{r^{\ell}}{(r^{\prime})^{\ell+1}}\sum_{m=-\ell}^{\ell}\frac{(\ell-m)!}{(\ell+m)!}{\mathsf{P}}_{\ell}^{m}(\cos\theta){\mathsf{P}}_{\ell}^{m}(\cos\theta^{\prime})e^{im(\phi-\phi^{\prime})},$

where $\|\mathbf{r}\|<\|\mathbf{r}^{\prime}\|$ ($\|\mathbf{r}\|$ denotes the
Euclidian norm of $\mathbf{r}\in\mathbb{R}^{3}$), and $r,\theta,\phi$,
$r^{\prime},\theta^{\prime},\phi^{\prime}$ are the spherical coordinates of
$\mathbf{r}$ and $\mathbf{r}^{\prime}$, respectively. The expansion (1.1)
contains the Ferrers function of the first kind (associated Legendre function
of the first kind on-the-cut) ${\mathsf{P}}_{\ell}^{m}$ [22, (14.3.1)]. We may
write expansion (1.1) in the more concise form

(1.2)
$\frac{1}{\|\mathbf{r}-\mathbf{r}^{\prime}\|}=\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}G_{\ell}^{m}(\mathbf{r})\overline{H_{\ell}^{m}(\mathbf{r}^{\prime})},$

where $G_{\ell}^{m}:\mathbb{R}^{3}\to\mathbb{C}$ is the internal spherical
harmonic

(1.3)
$G_{\ell}^{m}(\mathbf{r}):=\left(\frac{(\ell-m)!}{(\ell+m)!}\right)^{1/2}r^{\ell}{\mathsf{P}}_{\ell}^{m}(\cos\theta)e^{im\phi},$

and $H_{\ell}^{m}:\mathbb{R}^{3}\to\mathbb{C}$ is the external spherical
harmonic

(1.4)
$H_{\ell}^{m}(\mathbf{r}^{\prime}):=\left(\frac{(\ell-m)!}{(\ell+m)!}\right)^{1/2}(r^{\prime})^{-\ell-1}{\mathsf{P}}_{\ell}^{m}(\cos\theta^{\prime})e^{im\phi^{\prime}}.$

In this paper we derive expansions analogous to (1.2) for the 5-cyclidic
coordinate system [20, (6.24)] in place of spherical coordinates. The
coordinate surfaces of 5-cyclidic coordinates are triply-orthogonal confocal
cyclides. There are three kinds of internal and external 5-cyclidic harmonics,
one for each family of coordinate surfaces, and three corresponding
expansions. The authors already introduced internal 5-cyclidic harmonics in
[13]. As far as we know, the definition of external 5-cyclidic harmonics and
the expansions analogous to (1.2) are given in this paper for the first time.
We also derive some needed additional properties of internal 5-cyclidic
harmonics. In the definitions of internal and external spherical harmonics
(1.3), (1.4) there appear only the associated Legendre functions apart from
elementary functions. In the case of 5-cyclidic coordinates the definition of
internal and external harmonics requires solutions of a Fuchsian differential
equation with 5 regular singularities. The particular solutions of interest
are eigenfunctions of two-parameter Sturm-Liouville eigenvalue problems; see
[13].

In Maxime Bôcher’s 1891 dissertation, Ueber die Reihenentwickelungen der
Potentialtheorie [3], it was shown that the 3-variable Laplace equation can be
solved using separation of variables in seventeen conformally distinct quadric
and cyclidic coordinate systems. These coordinates have coordinate surfaces
which are zero sets for polynomials in $x,y,z$ with degree at most two and
four respectively. The Helmholtz equation on $\mathbb{R}^{3}$ admits simply
separable solutions in the same eleven quadric coordinate systems that the
Laplace equation admits separable solutions [14]. The Laplace equation also
admits ${R}$-separable solutions in an additional six conformally distinct
coordinate systems [20, Table 17, page 210]. Unlike the Laplace equation, the
Helmholtz equation does not admit solutions via ${R}$-separation of variables.
The appearance of ${R}$-separation is intrinsic to the existence of conformal
symmetries for a linear partial differential equation (see Boyer, Kalnins &
Miller (1976) [4]), i.e. dilatations, special conformal transformations,
inversions and reflections. The theory of separation of variables from a Lie
group theoretic viewpoint has been treated in Miller (1977) [20]. In Miller’s
book, separation of variables for the Laplace equation on $\mathbb{R}^{3}$ was
treated and the general asymmetric ${R}$-separable 5-cyclidic coordinate
system was introduced (see [20, Table 17, System 12]). In regard to this
coordinate system, and the corresponding separable harmonic solutions, Miller
indicates that “Very little is known about the solutions.”

To the authors’ knowledge, eigenfunction expansions for the fundamental
solution (the $1/r$ potential) have been obtained for the following coordinate
systems. See [11, 16, 18, 21] for expansions in spherical,
circular/parabolic/elliptic cylinder, oblate/prolate spheroidal, parabolic,
bi-spherical and toroidal coordinates. The expansion in confocal ellipsoidal
coordinates is treated in [2, 15]. This paper is a stepping-stone for
derivations of eigenfunction expansions for a fundamental solution of
Laplace’s equation in coordinate systems where these expansions are not known
such as paraboloidal, flat-ring cyclide, flat-disk cyclidic, bi-cyclide, cap-
cyclide and 3-cyclide [20, Table 17, System 13]) coordinates.

The eigenfunction expansions are often connected with integral identities
(such as the integral of Lipschitz [23, Section 13.2] and the Lipschitz-Hankel
integral [23, Section 13.21] which appear in cylindrical coordinates),
addition theorems (such as Neumann’s and Graf’s generalization of Neumann’s
addition theorem [23, Section 11.1, Section 11.3] which appear in cylindrical
coordinates and the addition theorem for spherical harmonics [24] which
appears in spherical coordinates), generating functions for orthogonal
polynomials (such as the generating function for Legendre polynomials [22,
(18.12.4)] which appears in spherical coordinates), and special function
expansion identities (such as Heine’s reciprocal square root identity [7,
(3.11)] which appears in circular cylindrical coordinates). In this setting,
one may perform eigenfunction expansions for a fundamental solution of
Laplace’s equation in alternative separable coordinate systems to obtain new
special function summation and integration identities which often have
interesting geometrical interpretations (see for instance [5, 9, 10]).
Eigenfunction expansions for fundamental solutions of elliptic partial
differential equations have been extended to more general separable linear
partial differential equations [6] and to partial differential equations on
Riemannian manifolds of constant curvature [8].

The outline of this paper is as follows. The 5-cyclidic coordinate system
$s_{1},s_{2},s_{3}$ is discussed in Section 2. In Section 3, we consider
internal and external 5-cyclidic harmonics of the second kind which are
related to the coordinate surfaces $s_{2}={\rm const}$. We start with
functions of the second kind because they are slightly easier to treat than
the harmonics of the first and third kind related to the coordinate surfaces
$s_{1}={\rm const}$, $s_{3}={\rm const}$, respectively. In Section 4, as one
of our main results, we obtain the expansion of the fundamental solution of
Laplace’s equation in terms of internal and external 5-cyclidic harmonics of
the second kind. The proof is based on $(a)$ an integral representation of the
external harmonics in terms of internal harmonics given in Section 4, and
$(b)$ the completeness property of internal harmonics obtained in [13]. In
Sections 5,6 we treat 5-cyclidic harmonics of the first kind. In Sections 7,8
we treat 5-cyclidic harmonics of the third kind.

##  2\. 5-cyclidic coordinates

We work on $\mathbb{R}^{3}$ with Cartesian coordinates $x,y,z$, and we use the
notations $\mathbf{r}=(x,y,z)$ and $\|\mathbf{r}\|=(x^{2}+y^{2}+z^{2})^{1/2}$.
Fix $a_{0}<a_{1}<a_{2}<a_{3}$. The 5-cyclidic coordinates of a point
$\mathbf{r}\in\mathbb{R}^{3}$ are the solutions $s=s_{1},s_{2},s_{3}$ of the
equation

(2.1)
$\frac{(\|\mathbf{r}\|^{2}-1)^{2}}{s-a_{0}}+\frac{4x^{2}}{s-a_{1}}+\frac{4y^{2}}{s-a_{2}}+\frac{4z^{2}}{s-a_{3}}=0$

(strictly speaking, this equation is multiplied by the common denominator of
the left-hand side), where

$a_{0}\leq s_{1}\leq a_{1}\leq s_{2}\leq a_{2}\leq s_{3}\leq a_{3};$

see [13, Section 4]. On the set

(2.2) $R:=\\{\mathbf{r}:x,y,z>0,\,\|\mathbf{r}\|<1\\},$

the map $(x,y,z)\in
R\mapsto(s_{1},s_{2},s_{3})\in(a_{0},a_{1})\times(a_{1},a_{2})\times(a_{2},a_{3})$
is bijective. The inverse map is given by

(2.3) $x=\frac{x_{1}}{1+x_{0}},\quad y=\frac{x_{2}}{1+x_{0}},\quad
z=\frac{x_{3}}{1+x_{0}},\\\ $

where

(2.4) $x_{j}^{2}=\frac{\prod_{i=1}^{3}(s_{i}-a_{j})}{\prod_{j\neq
i=0}^{3}(a_{i}-a_{j})},\quad x_{j}>0.$

We note that each $s_{i}$ is a continuous function on $\mathbb{R}^{3}$. Of
particular interest are the sets

$\displaystyle A_{1}$ $\displaystyle:=$
$\displaystyle\\{\mathbf{r}:s_{1}=s_{2}\\}=\\{(0,y,z):g_{1}(y,z)=0\\},$
$\displaystyle A_{2}$ $\displaystyle:=$
$\displaystyle\\{\mathbf{r}:s_{2}=s_{3}\\}=\\{(x,0,z):g_{2}(x,z)=0\\},$

where

$\displaystyle g_{1}(y,z)$ $\displaystyle:=$
$\displaystyle\frac{(y^{2}+z^{2}-1)^{2}}{a_{1}-a_{0}}+\frac{4y^{2}}{a_{1}-a_{2}}+\frac{4z^{2}}{a_{1}-a_{3}},$
$\displaystyle g_{2}(x,z)$ $\displaystyle:=$
$\displaystyle\frac{(x^{2}+z^{2}-1)^{2}}{a_{2}-a_{0}}+\frac{4x^{2}}{a_{2}-a_{1}}+\frac{4z^{2}}{a_{2}-a_{3}}.$

Each set $A_{1},A_{2}$ consists of two closed curves; see Figures 1, 2. The
function $s_{1}$ is (real-)analytic on $\mathbb{R}^{3}\setminus A_{1}$,
$s_{2}$ is analytic on $\mathbb{R}^{3}\setminus(A_{1}\cup A_{2})$, and $s_{3}$
is analytic on $\mathbb{R}^{3}\setminus A_{2}$. We will also encounter the
sets

$\displaystyle K_{1}$ $\displaystyle:=$
$\displaystyle\\{\mathbf{r}:\|\mathbf{r}\|<1,s_{1}=a_{1}\\}=\\{(0,y,z):y^{2}+z^{2}<1,g_{1}(y,z)\geq
0\\},$ $\displaystyle L_{1}$ $\displaystyle:=$
$\displaystyle\\{\mathbf{r}:s_{2}=a_{1}\\}=\\{(0,y,z):g_{1}(y,z)\leq 0\\},$
$\displaystyle M_{1}$ $\displaystyle:=$
$\displaystyle\\{\mathbf{r}:\|\mathbf{r}\|>1,s_{1}=a_{1}\\}=\\{(0,y,z):y^{2}+z^{2}>1,g_{1}(y,z)\geq
0\\},$ $\displaystyle K_{2}$ $\displaystyle:=$
$\displaystyle\\{\mathbf{r}:z>0,s_{3}=a_{2}\\}=\\{(x,0,z):z>0,g_{2}(x,z)\leq
0\\},$ $\displaystyle L_{2}$ $\displaystyle:=$
$\displaystyle\\{\mathbf{r}:s_{2}=a_{2}\\}=\\{(x,0,z):g_{2}(x,z)\geq 0\\},$
$\displaystyle M_{2}$ $\displaystyle:=$
$\displaystyle\\{\mathbf{r}:z<0,s_{3}=a_{2}\\}=\\{(x,0,z):z<0,g_{2}(x,z)\leq
0\\}.$

The sets $A_{1},K_{1},L_{1},M_{1}$ are subsets of the plane $x=0$, and
$A_{2},K_{2},L_{2},M_{2}$ are subsets of the plane $y=0$; see Figures 1, 2.

Figure 1. Curves $A_{1}$ and regions $K_{1},L_{1},M_{1}$ for $a_{j}=j$. Figure
2. Curves $A_{2}$ and regions $K_{2},L_{2},M_{2}$ for $a_{j}=j$.

We denote the inversion at the unit sphere on $\mathbb{R}^{3}$ by

(2.5) $\sigma_{0}(\mathbf{r}):=\|\mathbf{r}\|^{-2}\mathbf{r},$

and the reflections at the coordinate planes by

(2.6)
$\sigma_{1}(x,y,z):=(-x,y,z),\,\sigma_{2}(x,y,z):=(x,-y,z),\,\sigma_{3}(x,y,z):=(x,y,-z).$

We note that the functions $s_{1},s_{2},s_{3}$ are invariant under
$\sigma_{j}$, $j=0,1,2,3$.

We define auxiliary functions $\chi_{j}:\mathbb{R}^{3}\to\mathbb{R}$,
$j=0,1,2,3$, by

$\displaystyle\chi_{0}(\mathbf{r})$ $\displaystyle:=$ $\displaystyle{\rm
sgn}(1-\|\mathbf{r}\|)(s_{1}-a_{0})^{1/2},$
$\displaystyle\chi_{1}(\mathbf{r})$ $\displaystyle:=$ $\displaystyle{\rm
sgn}(x)((s_{2}-a_{1})(a_{1}-s_{1}))^{1/2},$
$\displaystyle\chi_{2}(\mathbf{r})$ $\displaystyle:=$ $\displaystyle{\rm
sgn}(y)((s_{3}-a_{2})(a_{2}-s_{2}))^{1/2},$
$\displaystyle\chi_{3}(\mathbf{r})$ $\displaystyle:=$ $\displaystyle{\rm
sgn}(z)(a_{3}-s_{3})^{1/2}.$

######  Lemma 2.1.

The functions $\chi_{j}$, $j=0,1,2,3$, are continuous on $\mathbb{R}^{3}$.
$\chi_{0},\chi_{2}$ are analytic on $\mathbb{R}^{3}\setminus A_{1}$, and
$\chi_{1},\chi_{3}$ are analytic on $\mathbb{R}^{3}\setminus A_{2}$. Moreover,

(2.7) $\chi_{j}\circ\sigma_{i}=\begin{cases}\chi_{j}&\text{if $i\neq j$},\\\
-\chi_{j}&\text{if $i=j$.}\end{cases}$

###### Proof.

Consider first $\chi_{3}$. The function $s_{3}$ is continuous, and
$s_{3}=a_{3}$ if and only if $z=0$. Therefore, $\chi_{3}$ is continuous. In
order to prove that $\chi_{3}$ is analytic on $\mathbb{R}^{3}\setminus A_{2}$,
it is enough to show that $\chi_{3}$ is analytic at every point of the plane
$z=0$. Let $\mathbf{r}_{0}=(x_{0},y_{0},0)$. There is $\epsilon\in(0,1)$ such
that $s_{3}\neq a_{2}$ for $\mathbf{r}\in
B_{\epsilon}(\mathbf{r}_{0})=\\{\mathbf{r}:\|\mathbf{r}-\mathbf{r}_{0}\|<\epsilon\\}$.
Then (2.1) with $s=s_{3}$ implies

$a_{3}-s_{3}=\frac{4z^{2}}{f(\mathbf{r})}\quad\text{for $\mathbf{r}\in
B_{\epsilon}(\mathbf{r}_{0})$},$

where

$f(\mathbf{r}):=\frac{(\|\mathbf{r}\|^{2}-1)^{2}}{s_{3}-a_{0}}+\frac{4x^{2}}{s_{3}-a_{1}}+\frac{4y^{2}}{s_{3}-a_{2}}$

is positive and analytic on $B_{\epsilon}(\mathbf{r}_{0})$. Therefore, we
obtain

$\chi_{3}(\mathbf{r})=\frac{2z}{(f(\mathbf{r}))^{1/2}}\quad\text{for
$\mathbf{r}\in B_{\epsilon}(\mathbf{r}_{0})$},$

and this shows that $\chi_{3}$ is analytic at $\mathbf{r}_{0}$. $\chi_{0}$ is
treated similarly.

Consider next $\chi_{2}$. The functions $s_{2},s_{3}$ are continuous, and
$(a_{2}-s_{2})(s_{3}-a_{2})=0$ if and only if $y=0$. Thus $\chi_{2}$ is
continuous. In order to prove that $\chi_{2}$ is analytic on
$\mathbb{R}^{3}\setminus A_{1}$, it is enough to show that $\chi_{2}$ is
analytic at all points of the plane $y=0$ which do not lie in $A_{1}$. Suppose
first $\mathbf{r}_{0}=(x_{0},0,z_{0})\in(K_{2}\cup M_{2})\setminus A_{2}$.
There is $\epsilon>0$ such that $s_{3}\neq a_{3}$ and $s_{2}\neq a_{2}$ for
$\mathbf{r}\in B_{\epsilon}(\mathbf{r}_{0})$. Then, by (2.1) with $s=s_{3}$,
we obtain

$s_{3}-a_{2}=\frac{4y^{2}}{g(\mathbf{r})},$

where

$g(\mathbf{r}):=-\frac{(\|\mathbf{r}\|^{2}-1)^{2}}{s_{3}-a_{0}}-\frac{4x^{2}}{s_{3}-a_{1}}-\frac{4z^{2}}{s_{3}-a_{3}}$

is analytic on $B_{\epsilon}(\mathbf{r}_{0})$. Since
$g(\mathbf{r}_{0})=-g_{2}(x_{0},z_{0})>0$, $g$ is also positive on
$B_{\epsilon}(\mathbf{r}_{0})$ for sufficiently small $\epsilon>0$. Then

$\chi_{2}(\mathbf{r})=(a_{2}-s_{2})^{1/2}\frac{2y}{(g(\mathbf{r}))^{1/2}}\quad\text{for
$\mathbf{r}\in B_{\epsilon}(\mathbf{r}_{0})$.}$

This shows that $\chi_{2}$ is analytic at $\mathbf{r}_{0}$ provided that
$\mathbf{r}_{0}\notin A_{1}$. In a similar way, by using (2.1) with $s=s_{2}$,
we show that $\chi_{2}$ is analytic at all points $\mathbf{r}_{0}\in
L_{2}\setminus A_{2}$. Finally, by subtracting equations (2.1) with
$s=s_{2},s_{3}$ from each other, we show that $\chi_{2}$ is analytic at all
points $\mathbf{r}_{0}\in A_{2}$. $\chi_{1}$ is treated similarly.

The symmetries (2.7) follow from the definition of $\chi_{j}$. ∎

Solving the Laplace equation

(2.8) $\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial
y^{2}}+\frac{\partial^{2}u}{\partial z^{2}}=0$

by the method of separation of variables, we find solutions

(2.9)
$u(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}w_{1}(s_{1})w_{2}(s_{2})w_{3}(s_{3}),\quad
s_{i}\in(a_{i-1},a_{i}).$

Each function $w=w_{1},w_{2},w_{3}$ satisfies the Fuchsian equation

(2.10)
$\prod_{j=0}^{3}(s-a_{j})\left[w^{\prime\prime}+\frac{1}{2}\sum_{j=0}^{3}\frac{1}{s-a_{j}}w^{\prime}\right]+\left(\frac{3}{16}s^{2}+\lambda_{1}s+\lambda_{2}\right)w=0,$

where $\lambda_{1},\lambda_{2}$ are separation constants; see [13]. This
equation has five regular singularities at $a_{0},a_{1},a_{2},a_{3},\infty$.
The exponents at each finite singularity are $0$ or $\frac{1}{2}$.

The function $u(\mathbf{r})$ defined in (2.9) is harmonic for all choices of
solutions $w_{i}$ to (2.10). However, it is harmonic only in the open set
obtained from $\mathbb{R}^{3}$ by removing the coordinate planes $x=0,y=0,z=0$
and the unit sphere $\|\mathbf{r}\|=1$. In order to obtain globally defined
harmonic functions we have to select the Frobenius solutions $w$ at the finite
singularities, that is, solutions that are either analytic at $a_{j}$ or of
the form $(s-a_{j})^{1/2}g(s)$ with $g(s)$ analytic at $s=a_{j}$. It is
impossible to choose the parameters $\lambda_{1},\lambda_{2}$ in such a way
that each solution $w_{i}$, $i=1,2,3$, is a nontrivial Frobenius solution
belonging to either one of the exponents $0$ or $\frac{1}{2}$ at both end
points $a_{i-1}$, $a_{i}$. If this were possible (2.9) would define a function
which is harmonic in the whole space $\mathbb{R}^{3}$ (as we see later) and
converges to $0$ as $\|\mathbf{r}\|\to\infty$. But such a function would have
to be identically zero. However, as shown in [13], we can determine special
values of $\lambda_{1},\lambda_{2}$ (eigenvalues) such that two solutions
(either (1) $w_{2}$, $w_{3}$, or (2) $w_{1}$, $w_{3}$, or (3) $w_{1}$,
$w_{2}$) are nontrivial Frobenius solution at both end points simultaneously.
These cases lead to 5-cyclidic harmonics of the first, second and third kind.
If the remaining function $w_{i}$ in case $(i)$ is chosen appropriately, we
obtain internal or external 5-cyclidic harmonics.

##  3\. 5-cyclidic harmonics of the second kind

In [13, Section VII] we introduced special solutions
$w_{i}(s_{i})=E^{(2)}_{i,\mathbf{n},\mathbf{p}}(s_{i})$ to equation (2.10) for
eigenvalues $\lambda_{j}=\lambda^{(2)}_{j,\mathbf{n},\mathbf{p}}$, $j=1,2$,
for every $\mathbf{n}\in\mathbb{N}_{0}^{2}$,
$\mathbf{p}=(p_{0},p_{1},p_{2},p_{3})\in\\{0,1\\}^{4}$. If
$\mathbf{n}=(n_{1},n_{3})$ then $n_{i}$ denotes the number of zeros of
$E^{(2)}_{i,\mathbf{n},\mathbf{p}}$ in $(a_{i-1},a_{i})$ for $i=1,3$. The
subscript $p_{j}$ describes the behavior of the solutions at the endpoint
$a_{j}$: We have

$E^{(2)}_{i,\mathbf{n},\mathbf{p}}(s_{i})=(s_{i}-a_{i-1})^{p_{i-1}/2}(a_{i}-s_{i})^{p_{i}/2}{\mathcal{E}}^{(2)}_{i,\mathbf{n},\mathbf{p}}(s_{i}),\quad
s_{i}\in(a_{i-1},a_{i}),$

where ${\mathcal{E}}^{(2)}_{1,\mathbf{n},\mathbf{p}}$ is analytic on
$[a_{0},a_{1}]$, ${\mathcal{E}}^{(2)}_{2,\mathbf{n},\mathbf{p}}$ is analytic
on $[a_{1},a_{2})$ (but not at $a_{2}$), and
${\mathcal{E}}^{(2)}_{3,\mathbf{n},\mathbf{p}}$ is analytic on
$[a_{2},a_{3}]$.

According to (2.9) the function

(3.1)
$G^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}E^{(2)}_{1,\mathbf{n},\mathbf{p}}(s_{1})E^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2})E^{(2)}_{3,\mathbf{n},\mathbf{p}}(s_{3}),\quad\mathbf{r}\in
R,$

is harmonic on $R$. In order to analytically extend
$G^{(2)}_{\mathbf{n},\mathbf{p}}$ we use the functions $\chi_{j}$ introduced
in Section 2. We set

(3.2)
$G^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}\prod_{j=0}^{3}(\chi_{j}(\mathbf{r}))^{p_{j}}\prod_{i=1}^{3}{\mathcal{E}}^{(2)}_{i,\mathbf{n},\mathbf{p}}(s_{i})\quad\text{if
$s_{2}\neq a_{2}$}$

which is consistent with (3.1). The condition $s_{2}\neq a_{2}$ is equivalent
to $\mathbf{r}\in\mathbb{R}^{3}\setminus L_{2}$. We call
$G^{(2)}_{\mathbf{n},\mathbf{p}}$ an internal 5-cyclidic harmonic of the
second kind.

######  Theorem 3.1.

Let $\mathbf{n}\in\mathbb{N}_{0}^{2}$ and $\mathbf{p}\in\\{0,1\\}^{4}$. Then
$G^{(2)}_{\mathbf{n},\mathbf{p}}$ is harmonic on $\mathbb{R}^{3}\setminus
L_{2}$. Moreover,

(3.3)
$G^{(2)}_{\mathbf{n},\mathbf{p}}(\sigma_{j}(\mathbf{r}))=(-1)^{p_{j}}G^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})\quad\text{for
$j=1,2,3$},$

and

(3.4)
$G^{(2)}_{\mathbf{n},\mathbf{p}}(\sigma_{0}(\mathbf{r}))=(-1)^{p_{0}}\|\mathbf{r}\|G^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}).$

###### Proof.

By (3.2) and Lemma 2.1, $G^{(2)}_{\mathbf{n},\mathbf{p}}$ is a composition of
continuous functions, and thus it is continuous on $\mathbb{R}^{3}\setminus
L_{2}$. As a composition of analytic functions,
$G^{(2)}_{\mathbf{n},\mathbf{p}}$ is analytic and thus harmonic on
$\mathbb{R}^{3}\setminus(A_{1}\cup L_{2})$. The set $A_{1}$ is a removable
line singularity of $G^{(2)}_{\mathbf{n},\mathbf{p}}$. This can be seen in two
different ways. 1) We may appeal to the general theory of harmonic functions.
$A_{1}$ is a polar set, and we may apply [1, Cor. 5.2.3]. 2) We can show
directly that $G^{(2)}_{\mathbf{n},\mathbf{p}}$ is analytic at each point of
$A_{1}$ by the method used in the proof of [13, Lemma 6.1]. For example, take
the simplest case $\mathbf{p}=(0,0,0,0)$. Then (3.1) holds for all
$\mathbf{r}\in\mathbb{R}^{3}\setminus L_{2}$, and the product
$E^{(2)}_{1,\mathbf{n},\mathbf{p}}(s_{1})E^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2})$
is analytic at each point of $A_{1}$. This is because
$E^{(2)}_{1,\mathbf{n},\mathbf{p}}(s)$ and
$E^{(2)}_{2,\mathbf{n},\mathbf{p}}(s)$ are analytic extensions of each other,
and $s_{1},$ $s_{2}$ enter symmetrically. Note that $s_{1}s_{2}$ and
$s_{1}+s_{2}$ are analytic at each point of $A_{1}$ although $s_{1}$, $s_{2}$
are not analytic there.

The symmetry properties of $G^{(2)}_{\mathbf{n},\mathbf{p}}$ also follow from
(3.2) and Lemma 2.1. ∎

If $U(\mathbf{r})$ is a harmonic function then its Kelvin transformation

$V(\mathbf{r})=\|\mathbf{r}\|^{-1}U(\sigma_{0}(\mathbf{r}))$

is also harmonic [17, page 232]. Equation (3.4) states that
$G^{(2)}_{\mathbf{n},\mathbf{p}}$ is invariant or changes sign under the
Kelvin transformation if $p_{0}=0$ or $p_{0}=1$, respectively. We see that
$L_{2}$ is a “surface singularity” of $G^{(2)}_{\mathbf{n},\mathbf{p}}$ which
is not removable (it is not a polar set). In fact,
$G^{(2)}_{\mathbf{n},\mathbf{p}}$ cannot be harmonic on $\mathbb{R}^{3}$
because it would be identically zero otherwise.

Let $F^{(2)}_{2,\mathbf{n},\mathbf{p}}$ be the Frobenius solution to the
Fuchsian equation (2.10) (with
$\lambda_{j}=\lambda^{(2)}_{j,\mathbf{n},\mathbf{p}}$) on $(a_{1},a_{2})$
belonging to the exponent $\frac{p_{2}}{2}$ at $s_{2}=a_{2}$, uniquely
determined by the Wronskian condition

(3.5)
$\omega(s)\left(E^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2})\frac{d}{ds_{2}}F^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2})-F^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2})\frac{d}{ds_{2}}E^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2})\right)=1,$

where

(3.6) $\omega(s):=\left|(s-a_{0})(s-a_{1})(s-a_{2})(s-a_{3})\right|^{1/2}.$

This definition is possible because we know that
$E^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2})$ is not a Frobenius solution
belonging to the exponent $\frac{p_{2}}{2}$ at $s_{2}=a_{2}$. Now we define
external 5-cyclidic harmonics of the second kind by

(3.7)
$H^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}E^{(2)}_{1,\mathbf{n},\mathbf{p}}(s_{1})F^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2})E^{(2)}_{3,\mathbf{n},\mathbf{p}}(s_{3}),\quad\mathbf{r}\in
R.$

In order to analytically extend $H^{(2)}_{\mathbf{n},\mathbf{p}}$ we write

$F^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2})=(s_{2}-a_{1})^{p_{1}/2}(a_{2}-s_{2})^{p_{2}/2}{\mathcal{F}}^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2}),\quad
s_{2}\in(a_{1},a_{2}),$

where ${\mathcal{F}}^{(2)}_{2,\mathbf{n},\mathbf{p}}$ is analytic on
$(a_{1},a_{2}]$ (but not at $a_{1}$). Then we define

(3.8)
$H^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}\prod_{j=0}^{3}(\chi_{j}(\mathbf{r}))^{p_{j}}{\mathcal{E}}^{(2)}_{1,\mathbf{n},\mathbf{p}}(s_{1}){\mathcal{F}}^{(2)}_{2,\mathbf{n},\mathbf{p}}(s_{2}){\mathcal{E}}^{(2)}_{3,\mathbf{n},\mathbf{p}}(s_{3})\quad\text{if
$s_{2}\neq a_{1}$}.$

The condition $s_{2}\neq a_{1}$ is equivalent to
$\mathbf{r}\in\mathbb{R}^{3}\setminus L_{1}$.

######  Theorem 3.2.

Let $\mathbf{n}\in\mathbb{N}_{0}^{2}$ and $\mathbf{p}\in\\{0,1\\}^{4}$. Then
$H^{(2)}_{\mathbf{n},\mathbf{p}}$ is harmonic on $\mathbb{R}^{3}\setminus
L_{1}$. The functions $H^{(2)}_{\mathbf{n},\mathbf{p}}$ share the symmetries
(3.3), (3.4) with $G^{(2)}_{\mathbf{n},\mathbf{p}}$. Moreover,

(3.9)
$H^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})=O(\|\mathbf{r}\|^{-1})\quad\text{as
$\|\mathbf{r}\|\to\infty$},$

and

(3.10) $\|\nabla
H^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})\|=O(\|\mathbf{r}\|^{-2})\quad\text{as
$\|\mathbf{r}\|\to\infty$}.$

###### Proof.

The proof of analyticity and symmetry of $H^{(2)}_{\mathbf{n},\mathbf{p}}$ is
similar to that given for $G^{(2)}_{\mathbf{n},\mathbf{p}}$ in Theorem 3.1,
and is omitted. Estimates (3.9) and (3.10) follow easily from the observation
that the Kelvin transformation of $H^{(2)}_{\mathbf{n},\mathbf{p}}$ is $\pm
H^{(2)}_{\mathbf{n},\mathbf{p}}$ which is analytic at $\mathbf{0}\notin
L_{1}$. ∎

##  4\. Expansion of the reciprocal distance in 5-cyclidic harmonics of second
kind

For given $d_{2}\in(a_{1},a_{2})$ we consider the “5-cyclidic ring”

(4.1) $D_{2}:=\\{\mathbf{r}\in\mathbb{R}^{3}:s_{2}<d_{2}\\},$

or, equivalently,

(4.2)
$D_{2}=\\{\mathbf{r}:\frac{(\|\mathbf{r}\|^{2}-1)^{2}}{d_{2}-a_{0}}+\frac{4x^{2}}{d_{2}-a_{1}}+\frac{4y^{2}}{d_{2}-a_{2}}+\frac{4z^{2}}{d_{2}-a_{3}}<0\\}.$

Note that each internal 5-cyclidic harmonic $G^{(2)}_{\mathbf{n},\mathbf{p}}$
is harmonic in $D_{2}$ (and on its boundary), and each external 5-cyclidic
harmonic is harmonic on $\mathbb{R}^{3}\setminus D_{2}$ (and on its boundary).

We represent external harmonics in terms of internal harmonics by a surface
integral over the boundary $\partial D_{2}$ of the ring $D_{2}$ as follows.

######  Theorem 4.1.

Let $d_{2}\in(a_{1},a_{2})$, $\mathbf{n}\in\mathbb{N}_{0}^{2}$,
$\mathbf{p}\in\\{0,1\\}^{4}$. Then

(4.3)
$H^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}^{\prime})=\frac{1}{4\pi\omega(d_{2})\\{E^{(2)}_{2,\mathbf{n},\mathbf{p}}(d_{2})\\}^{2}}\int_{\partial
D_{2}}\frac{G^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})}{h_{2}(\mathbf{r})\|\mathbf{r}-\mathbf{r}^{\prime}\|}\,dS(\mathbf{r})$

for all $\mathbf{r}^{\prime}\in\mathbb{R}^{3}\setminus\bar{D}_{2}$. The scale
factor $h_{2}$ is given by

(4.4)
$16\\{h_{2}(\mathbf{r})\\}^{2}=\frac{(\|\mathbf{r}\|^{2}-1)^{2}}{(d_{2}-a_{0})^{2}}+\frac{4x^{2}}{(d_{2}-a_{1})^{2}}+\frac{4y^{2}}{(d_{2}-a_{2})^{2}}+\frac{4z^{2}}{(d_{2}-a_{3})^{2}}.$

###### Proof.

Let $D$ be an open bounded subset of $\mathbb{R}^{3}$ with smooth boundary.
For $u,v\in C^{2}(\bar{D})$, Green’s formula states that

(4.5) $\int_{D}(u\Delta v-v\Delta u)\,d\mathbf{r}=\int_{\partial
D}\left(u\frac{\partial v}{\partial\nu}-v\frac{\partial
u}{\partial\nu}\right)\,dS,$

where $\frac{\partial u}{\partial\nu}$ is the outward normal derivative of $u$
on the boundary $\partial D$ of $D$.

We apply (4.5) to the domain $D=D_{2}$, and functions
$u=G=G^{(2)}_{\mathbf{n},\mathbf{p}}$,
$v(\mathbf{r})=\frac{1}{4\pi\|\mathbf{r}-\mathbf{r}^{\prime}\|}$. Since $u,v$
are harmonic on an open set containing $\bar{D}_{2}$ we obtain

(4.6) $0=\int_{\partial D_{2}}\left(G\frac{\partial
v}{\partial\nu}-v\frac{\partial G}{\partial\nu}\right)dS.$

We now use (4.5) a second time. We choose $R>0$ so large that the ball
$B_{R}(\mathbf{0})$ contains $\mathbf{r}^{\prime}$ and $\bar{D}_{2}$. Then we
take $D=B_{R}(\mathbf{0})-\bar{D}_{2}-B_{\epsilon}(\mathbf{r}^{\prime})$ with
small radius $\epsilon>0$. Take $u=H=H^{(2)}_{\mathbf{n},\mathbf{p}}$ and $v$
as before. Note that $u,v$ are harmonic on an open set containing $\bar{D}$.
By a standard argument [19, Theorem 1, page 109], taking the limit
$\epsilon\to 0$, we obtain

(4.7) $H(\mathbf{r}^{\prime})=\int_{\partial
B_{R}(\mathbf{0})}\left(H\frac{\partial v}{\partial\nu}-v\frac{\partial
H}{\partial\nu}\right)dS-\int_{\partial D_{2}}\left(H\frac{\partial
v}{\partial\nu}-v\frac{\partial H}{\partial\nu}\right)dS,$

where, in the second integral, $\frac{\partial}{\partial\nu}$ denotes the same
derivative as in (4.6). The first integral in (4.7) tends to $0$ as
$R\to\infty$ by (3.9), (3.10). Therefore,

(4.8) $H(\mathbf{r}^{\prime})=-\int_{\partial D_{2}}\left(H\frac{\partial
v}{\partial\nu}-v\frac{\partial H}{\partial\nu}\right)dS.$

We now multiply (4.6) by $F_{2}(d_{2})$,
$F_{2}:=F^{(2)}_{2,\mathbf{n},\mathbf{p}}$, then multiply (4.8) by
$E_{2}(d_{2})$, $E_{i}:=E^{(2)}_{i,\mathbf{n},\mathbf{p}}$, and add these
equations. By (3.1) and (3.7) we have

$F_{2}(d_{2})G(\mathbf{r})=E_{2}(d_{2})H(\mathbf{r}),\quad\mathbf{r}\in\partial
D_{2},$

first for $\mathbf{r}\in\partial D_{2}\cap R$ but then for all
$\mathbf{r}\in\partial D_{2}$ by shared symmetries (3.3), (3.4) of $G,H$.
Therefore, we find

(4.9) $E_{2}(d_{2})H(\mathbf{r}^{\prime})=\int_{\partial
D_{2}}v\left(E_{2}(d_{2})\frac{\partial
H}{\partial\nu}-F_{2}(d_{2})\frac{\partial G}{\partial\nu}\right)dS.$

The normal derivative and the derivative with respect to $s_{2}$ are related
by

$\frac{\partial}{\partial\nu}=\frac{1}{h_{2}}\frac{\partial}{\partial s_{2}},$

where $h_{2}$ is the scale factor of the 5-cyclidic coordinate $s_{2}$ given
by (4.4); see [13, (22)]. Let $\mathbf{r}\in\partial D_{2}\cap R$ with
5-cyclidic coordinates $s_{1},s_{2}=d_{2},s_{3}$. Then

$\displaystyle\left(E_{2}(d_{2})\frac{\partial
H}{\partial\nu}-F_{2}(d_{2})\frac{\partial G}{\partial\nu}\right)(\mathbf{r})$
$\displaystyle\hskip
36.98866pt=E_{2}(d_{2})\frac{\partial(\|\mathbf{r}\|^{2}+1)^{-1/2}}{\partial\nu}E_{1}(s_{1})F_{2}(d_{2})E_{3}(s_{3})$
$\displaystyle\hskip
65.44142pt+E_{2}(d_{2})(\|\mathbf{r}\|^{2}+1)^{-1/2}h_{2}^{-1}E_{1}(s_{1})F_{2}^{\prime}(d_{2})E_{3}(s_{3})$
$\displaystyle\hskip 65.44142pt-
F_{2}(d_{2})\frac{\partial(\|\mathbf{r}\|^{2}+1)^{-1/2}}{\partial\nu}E_{1}(s_{1})E_{2}(d_{2})E_{3}(s_{3})$
$\displaystyle\hskip 65.44142pt-
F_{2}(d_{2})(\|\mathbf{r}\|^{2}+1)^{-1/2}h_{2}^{-1}E_{1}(s_{1})E_{2}^{\prime}(d_{2})E_{3}(s_{3})$
$\displaystyle\hskip
36.98866pt=h_{2}^{-1}(\|\mathbf{r}\|^{2}+1)^{-1/2}E_{1}(s_{1})\left\\{E_{2}(d_{2})F_{2}^{\prime}(d_{2})-E_{2}^{\prime}(d_{2})F_{2}(d_{2})\right\\}E_{3}(s_{3}).$

We now use (3.5) and obtain

(4.10) $\left(E_{2}(d_{2})\frac{\partial
H}{\partial\nu}-F_{2}(d_{2})\frac{\partial
G}{\partial\nu}\right)(\mathbf{r})=\frac{G(\mathbf{r})}{h_{2}(\mathbf{r})\omega(d_{2})E_{2}(d_{2})},$

which holds for all $\mathbf{r}\in\partial D_{2}$ because $G$ and $H$ share
the symmetries (3.3), (3.4). When we substitute (4.10) in (4.9) we arrive at
(4.3) ∎

We obtain the expansion of the reciprocal distance in 5-cyclidic harmonics.

######  Theorem 4.2.

Let $\mathbf{r},\mathbf{r}^{\prime}\in\mathbb{R}^{3}$ with 5-cyclidic
coordinates $s_{2},s_{2}^{\prime}$, respectively. If $s_{2}<s_{2}^{\prime}$
then

(4.11)
$\frac{1}{\|\mathbf{r}-\mathbf{r}^{\prime}\|}=\pi\sum_{\mathbf{n}\in\mathbb{N}_{0}^{2}}\sum_{\mathbf{p}\in\\{0,1\\}^{4}}G^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})H^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}^{\prime}).$

###### Proof.

We pick $d_{2}$ such that $s_{2}<d_{2}<s_{2}^{\prime}$, and consider the
domain $D_{2}$ defined in (4.1). The function
$f(\mathbf{q}):=\|\mathbf{q}-\mathbf{r}^{\prime}\|^{-1}$ is harmonic on an
open set containing $\bar{D}_{2}$. Therefore, by [13, (95),(97)], we have

(4.12)
$\frac{1}{\|\mathbf{r}-\mathbf{r}^{\prime}\|}=\sum_{\mathbf{n}\in\mathbb{N}_{0}^{2}}\sum_{\mathbf{p}\in\\{0,1\\}^{4}}d^{(2)}_{\mathbf{n},\mathbf{p}}G^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),$

where

$d^{(2)}_{\mathbf{n},\mathbf{p}}:=\frac{1}{4\omega(d_{2})\\{E^{(2)}_{2,\mathbf{n},\mathbf{p}}(d_{2})\\}^{2}}\int_{\partial
D_{2}}\frac{G^{(2)}_{\mathbf{n},\mathbf{p}}(\mathbf{q})}{h_{2}(\mathbf{q})\|\mathbf{q}-\mathbf{r}^{\prime}\|}\,dS(\mathbf{q}).$

Using Theorem 4.1, we obtain (4.11). ∎

##  5\. 5-cyclidic harmonics of the first kind

In [13, Section V] we introduced special solutions
$w_{i}(s_{i})=E^{(1)}_{i,\mathbf{n},\mathbf{p}}(s_{i})$ to equation (2.10) for
eigenvalues $\lambda_{j}=\lambda^{(1)}_{j,\mathbf{n},\mathbf{p}}$, $j=1,2$,
for every $\mathbf{n}\in\mathbb{N}_{0}^{2}$,
$\mathbf{p}=(p_{1},p_{2},p_{3})\in\\{0,1\\}^{3}$. These functions have the
form

$\displaystyle E^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})$ $\displaystyle=$
$\displaystyle(a_{1}-s_{1})^{p_{1}/2}{\mathcal{E}}^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1}),\quad
s_{1}\in(a_{0},a_{1}),$ $\displaystyle
E^{(1)}_{i,\mathbf{n},\mathbf{p}}(s_{i})$ $\displaystyle=$
$\displaystyle(s_{i}-a_{i-1})^{p_{i-1}/2}(a_{i}-s_{i})^{p_{i}/2}{\mathcal{E}}^{(1)}_{i,\mathbf{n},\mathbf{p}}(s_{i}),\quad
s_{i}\in(a_{i-1},a_{i}),i=2,3,$

where ${\mathcal{E}}^{(1)}_{1,\mathbf{n},\mathbf{p}}$ is analytic on
$(a_{0},a_{1}]$ (but not at $a_{0}$) while
${\mathcal{E}}^{(1)}_{i,\mathbf{n},\mathbf{p}}$ is analytic on
$[a_{i-1},a_{i}]$ for $i=2,3$. As in [13, Section VI] we define the internal
5-cyclidic harmonic of the first kind by

(5.1)
$G^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}E^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})E^{(1)}_{2,\mathbf{n},\mathbf{p}}(s_{2})E^{(1)}_{3,\mathbf{n},\mathbf{p}}(s_{3}),\quad\mathbf{r}\in
R.$

According to (2.9), $G^{(1)}_{\mathbf{n},\mathbf{p}}$ is a harmonic function
in the region $R$. In order to analytically extend
$G^{(1)}_{\mathbf{n},\mathbf{p}}$ to a larger domain of definition, some
preparations are necessary.

Let $P^{(1)}_{1,\mathbf{n},\mathbf{p}}$ be the solution to (2.10) (with
$\lambda_{j}=\lambda^{(1)}_{j,\mathbf{n},\mathbf{p}}$) on $(a_{0},a_{1})$
belonging to the exponent $0$ at $s=a_{0}$ and uniquely determined by the
condition $P^{(1)}_{1,\mathbf{n},\mathbf{p}}(a_{0})=1$. We write

$P^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})=(a_{1}-s_{1})^{p_{1}/2}{\mathcal{P}}^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1}),\quad
s_{1}\in(a_{0},a_{1}),$

where ${\mathcal{P}}^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})$ is analytic on
$[a_{0},a_{1})$. Then using the functions $\chi_{j}$ from Section 2 we define

(5.2)
$I^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}\prod_{j=1}^{3}(\chi_{j}(\mathbf{r}))^{p_{j}}{\mathcal{P}}^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})\prod_{i=2}^{3}{\mathcal{E}}^{(1)}_{i,\mathbf{n},\mathbf{p}}(s_{i})\quad\text{if
$s_{1}\neq a_{1}$}.$

The condition $s_{1}\neq a_{1}$ is equivalent to
$\mathbf{r}\in\mathbb{R}^{3}\setminus(K_{1}\cup M_{1})$; see Figure 1.

Similarly, let $Q^{(1)}_{1,\mathbf{n},\mathbf{p}}$ be the solution to (2.10)
(with $\lambda_{j}=\lambda^{(1)}_{j,\mathbf{n},\mathbf{p}}$) on
$(a_{0},a_{1})$ belonging to the exponent $\frac{1}{2}$ at $s=a_{0}$ and
uniquely determined by the condition $\lim_{s_{1}\to
a_{0}^{+}}\omega(s_{1})\frac{d}{ds_{1}}Q^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})=1$.
We write

$Q^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})=(s_{1}-a_{0})^{1/2}(a_{1}-s_{1})^{p_{1}/2}{\mathcal{Q}}^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1}),\quad
s_{1}\in(a_{0},a_{1}),$

where ${\mathcal{Q}}^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})$ is analytic on
$[a_{0},a_{1})$. Then we define

(5.3)
$J^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}\chi_{0}(\mathbf{r})\prod_{j=1}^{3}(\chi_{j}(\mathbf{r}))^{p_{j}}{\mathcal{Q}}^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})\prod_{i=2}^{3}{\mathcal{E}}^{(1)}_{i,\mathbf{n},\mathbf{p}}(s_{i})\quad\text{if
$s_{1}\neq a_{1}$}.$

######  Lemma 5.1.

The functions $I^{(1)}_{\mathbf{n},\mathbf{p}}$ and
$J^{(1)}_{\mathbf{n},\mathbf{p}}$ are harmonic on
$\mathbb{R}^{3}\setminus(K_{1}\cup M_{1})$. They have the symmetries

(5.4) $\displaystyle I^{(1)}_{\mathbf{n},\mathbf{p}}(\sigma_{0}(\mathbf{r}))$
$\displaystyle=$
$\displaystyle\|\mathbf{r}\|I^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),$
(5.5) $\displaystyle I^{(1)}_{\mathbf{n},\mathbf{p}}(\sigma_{j}(\mathbf{r}))$
$\displaystyle=$
$\displaystyle(-1)^{p_{j}}I^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),\quad
j=1,2,3,$ (5.6) $\displaystyle
J^{(1)}_{\mathbf{n},\mathbf{p}}(\sigma_{0}(\mathbf{r}))$ $\displaystyle=$
$\displaystyle-\|\mathbf{r}\|J^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),$
(5.7) $\displaystyle J^{(1)}_{\mathbf{n},\mathbf{p}}(\sigma_{j}(\mathbf{r}))$
$\displaystyle=$
$\displaystyle(-1)^{p_{j}}J^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),\quad
j=1,2,3.$

###### Proof.

By definition (5.2), $I^{(1)}_{\mathbf{n},\mathbf{p}}$ is a composition of
continuous functions provided $s_{1}\neq a_{1}$, that is,
$I^{(1)}_{\mathbf{n},\mathbf{p}}$ is continuous on
$\mathbb{R}^{3}\setminus(K_{1}\cup M_{1})$. $I^{(1)}_{\mathbf{n},\mathbf{p}}$
is also a composition of analytic functions provided $s_{1}\neq a_{1}$ and
$s_{2}\neq s_{3}$, that is, $I^{(1)}_{\mathbf{n},\mathbf{p}}$ is analytic on
$\mathbb{R}^{3}\setminus(K_{1}\cup M_{1}\cup A_{2})$. Thus it is also harmonic
on $\mathbb{R}^{3}\setminus(K_{1}\cup M_{1}\cup A_{2})$. By the same argument
as in the proof of Theorem 3.1, $A_{2}$ is a removable singularity of
$I^{(1)}_{\mathbf{n},\mathbf{p}}$. Thus $I^{(1)}_{\mathbf{n},\mathbf{p}}$ is
harmonic on $\mathbb{R}^{3}\setminus(K_{1}\cup M_{1})$. The proof that
$J^{(1)}_{\mathbf{n},\mathbf{p}}$ is harmonic on
$\mathbb{R}^{3}\setminus(K_{1}\cup M_{1})$ is analogous. The symmetry
properties follow from (5.2), (5.3) and Lemma 2.1. ∎

Since $P^{(1)}_{1,\mathbf{n},\mathbf{p}},Q^{(1)}_{1,\mathbf{n},\mathbf{p}}$
form a fundamental system of solutions to (2.10) (with
$\lambda_{j}=\lambda^{(1)}_{j,\mathbf{n},\mathbf{p}}$) on $(a_{0},a_{1})$,
there are (nonzero) scalars $\alpha^{(1)}_{\mathbf{n},\mathbf{p}}$,
$\beta^{(1)}_{\mathbf{n},\mathbf{p}}$ such that

$E^{(1)}_{1,\mathbf{n},\mathbf{p}}=\alpha^{(1)}_{\mathbf{n},\mathbf{p}}P^{(1)}_{1,\mathbf{n},\mathbf{p}}+\beta^{(1)}_{\mathbf{n},\mathbf{p}}Q^{(1)}_{1,\mathbf{n},\mathbf{p}}.$

This leads us to the global definition of internal 5-cyclidic harmonics of the
first kind

(5.8)
$G^{(1)}_{\mathbf{n},\mathbf{p}}:=\alpha^{(1)}_{\mathbf{n},\mathbf{p}}I^{(1)}_{\mathbf{n},\mathbf{p}}+\beta^{(1)}_{\mathbf{n},\mathbf{p}}J^{(1)}_{\mathbf{n},\mathbf{p}}$

which is consistent with (5.1). We also note that, if $\|\mathbf{r}\|<1$ and
$\mathbf{r}\not\in K_{1}$, then (5.2), (5.3), (5.8) imply that

(5.9)
$G^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})=(\|\mathbf{r}\|^{2}+1)^{-1/2}\prod_{j=1}^{3}(\chi_{j}(\mathbf{r}))^{p_{j}}\prod_{i=1}^{3}{\mathcal{E}}^{(1)}_{i,\mathbf{n},\mathbf{p}}(s_{i}).$

######  Theorem 5.2.

Let $\mathbf{n}\in\mathbb{N}_{0}^{2}$ and
$\mathbf{p}=(p_{1},p_{2},p_{3})\in\\{0,1\\}^{3}$. Then
$G^{(1)}_{\mathbf{n},\mathbf{p}}$ extends continuously to a harmonic function
on $\mathbb{R}^{3}\setminus M_{1}$. Moreover,

(5.10)
$G^{(1)}_{\mathbf{n},\mathbf{p}}(\sigma_{j}(\mathbf{r}))=(-1)^{p_{j}}G^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})\quad\text{for
$j=1,2,3$}.$

###### Proof.

By Lemma 5.1, $G^{(1)}_{\mathbf{n},\mathbf{p}}$ is harmonic on
$\mathbb{R}^{3}\setminus(K_{1}\cup M_{1})$. If $\|\mathbf{r}\|<1$ we have
$s_{1}\neq a_{0}$. Therefore, the right-hand side of (5.9) is continuous on
the ball $B_{1}(\mathbf{0})$ and harmonic on
$B_{1}(\mathbf{0})\setminus(A_{1}\cup A_{2})$. Thus it is harmonic on
$B_{1}(\mathbf{0})$ which proves the first part of the statement of the
theorem. The symmetries follow from (5.5), (5.7). ∎

It will be useful to introduce another solution to (2.10) by

(5.11)
$F^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1}):=\gamma^{(1)}_{\mathbf{n},\mathbf{p}}\left(\alpha^{(1)}_{\mathbf{n},\mathbf{p}}P^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})-\beta^{(1)}_{\mathbf{n},\mathbf{p}}Q^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})\right),\quad
s_{1}\in(a_{0},a_{1}).$

We determine $\gamma^{(1)}_{\mathbf{n},\mathbf{p}}$ from the Wronskian

(5.12)
$\omega(s_{1})\left(E^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})\frac{d}{ds_{1}}F^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})-F^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})\frac{d}{ds_{1}}E^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})\right)=1$

which is equivalent to

$\gamma^{(1)}_{\mathbf{n},\mathbf{p}}=\frac{-1}{2\alpha^{(1)}_{\mathbf{n},\mathbf{p}}\beta^{(1)}_{\mathbf{n},\mathbf{p}}}.$

We define external 5-cyclidic harmonics of the first kind by

(5.13)
$H^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=\gamma^{(1)}_{\mathbf{n},\mathbf{p}}\|\mathbf{r}\|^{-1}G^{(1)}_{\mathbf{n},\mathbf{p}}(\sigma_{0}(\mathbf{r}))\quad\text{for
$\mathbf{r}\in\mathbb{R}^{3}\setminus K_{1}$}.$

The reason to include the factor $\gamma^{(1)}_{\mathbf{n},\mathbf{p}}$ is
that we aim for a simple form of the expansion formula (6.4). In particular,
we have

(5.14)
$H^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})=(\|\mathbf{r}\|^{2}+1)^{-1/2}F^{(1)}_{1,\mathbf{n},\mathbf{p}}(s_{1})E^{(1)}_{2,\mathbf{n},\mathbf{p}}(s_{2})E^{(1)}_{3,\mathbf{n},\mathbf{p}}(s_{3})\quad\text{for
$\mathbf{r}\in R$}.$

We notice an important difference between 5-cyclidic harmonics of the first
and second kind (considered in Section 3). The external 5-cyclidic harmonics
of the first kind are simply the Kelvin transformations of the internal
5-cyclidic harmonics of the first kind up to a constant factor. There is no
such simple relationship between internal and external 5-cyclidic harmonics of
the second kind.

######  Theorem 5.3.

Let $\mathbf{n}\in\mathbb{N}_{0}^{2}$ and
$\mathbf{p}=(p_{1},p_{2},p_{3})\in\\{0,1\\}^{3}$. Then
$H^{(1)}_{\mathbf{n},\mathbf{p}}$ is harmonic on $\mathbb{R}^{3}\setminus
K_{1}$. The functions $H^{(1)}_{\mathbf{n},\mathbf{p}}$ share the symmetries
(5.10) with $G^{(1)}_{\mathbf{n},\mathbf{p}}$. Moreover,

(5.15)
$H^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})=O(\|\mathbf{r}\|^{-1})\quad\text{as
$\|\mathbf{r}\|\to\infty$},$

and

(5.16) $\|\nabla
H^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})\|=O(\|\mathbf{r}\|^{-2})\quad\text{as
$\|\mathbf{r}\|\to\infty$}.$

###### Proof.

The proof of analyticity and symmetry follows directly from (5.13) and Theorem
5.2. Estimates (5.15) and (5.16) follow from the fact that the Kelvin
transformation of $H^{(1)}_{\mathbf{n},\mathbf{p}}$ is analytic at the origin.
∎

##  6\. Expansion of the reciprocal distance in 5-cyclidic harmonics of first
kind

For fixed $s\in(a_{0},a_{1})$ the coordinate surface (2.1) consists of two
closed surfaces of genus $0$. One lies inside the unit ball
$B_{1}(\mathbf{0})$ and the other one is obtained from it by inversion
$\sigma_{0}$. We consider the region $D_{1}$ interior to the coordinate
surface $s=d_{1}$ which lies in $B_{1}(\mathbf{0})$:

(6.1) $D_{1}:=\\{\mathbf{r}\in\mathbb{R}^{3}:\|\mathbf{r}\|<1,s_{1}>d_{1}\\}.$

######  Theorem 6.1.

Let $d_{1}\in(a_{0},a_{1})$, $\mathbf{n}\in\mathbb{N}_{0}^{2}$,
$\mathbf{p}\in\\{0,1\\}^{3}$. Then

(6.2)
$H^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}^{\prime})=\frac{1}{4\pi\omega(d_{1})\\{E^{(1)}_{1,\mathbf{n},\mathbf{p}}(d_{1})\\}^{2}}\int_{\partial
D_{1}}\frac{G^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})}{h_{1}(\mathbf{r})\|\mathbf{r}-\mathbf{r}^{\prime}\|}\,dS(\mathbf{r})$

for all $\mathbf{r}^{\prime}\in\mathbb{R}^{3}\setminus\bar{D}_{1}$. The scale
factor $h_{1}$ is given by

(6.3)
$16\\{h_{1}(\mathbf{r})\\}^{2}=\frac{(\|\mathbf{r}\|^{2}-1)^{2}}{(d_{1}-a_{0})^{2}}+\frac{4x^{2}}{(d_{1}-a_{1})^{2}}+\frac{4y^{2}}{(d_{1}-a_{2})^{2}}+\frac{4z^{2}}{(d_{1}-a_{3})^{2}}.$

###### Proof.

The proof is similar to the proof of Theorem 4.1. We use (5.1), (5.14) and the
Wronskian (5.12). ∎

We obtain the expansion of the reciprocal distance in 5-cyclidic harmonics of
first kind.

######  Theorem 6.2.

Let $\mathbf{r},\mathbf{r}^{\prime}\in\mathbb{R}^{3}$ with 5-cyclidic
coordinates $s_{1},s_{1}^{\prime}$, respectively. If either (a)
$\|\mathbf{r}\|,\|\mathbf{r}^{\prime}\|\leq 1$, $s_{1}>s_{1}^{\prime}$, or (b)
$\|\mathbf{r}\|<1<\|\mathbf{r}^{\prime}\|$, or (c)
$\|\mathbf{r}\|,\|\mathbf{r}^{\prime}\|\geq 1$, $s_{1}<s_{1}^{\prime}$, then

(6.4)
$\frac{1}{\|\mathbf{r}-\mathbf{r}^{\prime}\|}=2\pi\sum_{\mathbf{n}\in\mathbb{N}_{0}^{2}}\sum_{\mathbf{p}\in\\{0,1\\}^{3}}G^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})H^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}^{\prime}).$

###### Proof.

Suppose (a) or (b) holds. Pick $d_{1}$ such that $s_{1}^{\prime}<d_{1}<s_{1}$
if (a) holds, or such that $a_{0}<d_{1}<s_{1}$ if (b) holds. Then consider the
domain $D_{1}$ defined in (6.1). The function
$f(\mathbf{q}):=\|\mathbf{q}-\mathbf{r}^{\prime}\|^{-1}$ is harmonic on an
open set containing $\bar{D}_{1}$. Therefore, by [13, (71),(73)], we have

(6.5)
$\frac{1}{\|\mathbf{r}-\mathbf{r}^{\prime}\|}=\sum_{\mathbf{n}\in\mathbb{N}_{0}^{2}}\sum_{\mathbf{p}\in\\{0,1\\}^{3}}d^{(1)}_{\mathbf{n},\mathbf{p}}G^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),$

where

$d^{(1)}_{\mathbf{n},\mathbf{p}}=\frac{1}{2\omega(d_{1})\\{E^{(1)}_{1,\mathbf{n},\mathbf{p}}(d_{1})\\}^{2}}\int_{\partial
D_{1}}\frac{G^{(1)}_{\mathbf{n},\mathbf{p}}(\mathbf{q})}{h_{1}(\mathbf{q})\|\mathbf{q}-\mathbf{r}^{\prime}\|}\,dS(\mathbf{q}).$

Using Theorem 6.1, we obtain (6.4).

Now suppose (c) holds. Then the points $\sigma_{0}(\mathbf{r}^{\prime})$,
$\sigma_{0}(\mathbf{r})$ in place of $\mathbf{r},\mathbf{r}^{\prime}$ satisfy
(a), so, by what we already proved,

$\frac{1}{\|\sigma_{0}(\mathbf{r})-\sigma_{0}(\mathbf{r}^{\prime})\|}=2\pi\sum_{\mathbf{n}\in\mathbb{N}_{0}^{2}}\sum_{\mathbf{p}\in\\{0,1\\}^{3}}G^{(1)}_{\mathbf{n},\mathbf{p}}(\sigma_{0}(\mathbf{r}^{\prime}))H^{(1)}_{\mathbf{n},\mathbf{p}}(\sigma_{0}(\mathbf{r})).$

This gives (6.4) by using (5.13) and observing that

$\|\mathbf{r}-\mathbf{r}^{\prime}\|=\|\mathbf{r}\|\|\mathbf{r}^{\prime}\|\|\sigma_{0}(\mathbf{r})-\sigma_{0}(\mathbf{r}^{\prime})\|.$

∎

##  7\. 5-cyclidic harmonics of the third kind

The 5-cyclidic harmonics of the third kind are treated analogously to the
harmonics of the first kind. Therefore, we will omit all proofs in the
following two sections.

In [13, Section IX] we introduced special solutions
$w_{i}(s_{i})=E^{(3)}_{i,\mathbf{n},\mathbf{p}}(s_{i})$ to equation (2.10) for
eigenvalues $\lambda_{j}=\lambda^{(3)}_{j,\mathbf{n},\mathbf{p}}$, $j=1,2$,
for every $\mathbf{n}\in\mathbb{N}_{0}^{2}$,
$\mathbf{p}=(p_{0},p_{1},p_{2})\in\\{0,1\\}^{3}$. These functions have the
form

$\displaystyle E^{(3)}_{i,\mathbf{n},\mathbf{p}}(s_{i})$ $\displaystyle=$
$\displaystyle(s_{i}-a_{i-1})^{p_{i-1}/2}(a_{i}-s_{i})^{p_{i}/2}{\mathcal{E}}^{(3)}_{i,\mathbf{n},\mathbf{p}}(s_{i}),\quad
s_{i}\in(a_{i-1},a_{i}),\ i=1,2,$ $\displaystyle
E^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})$ $\displaystyle=$
$\displaystyle(s_{3}-a_{2})^{p_{2}/2}{\mathcal{E}}^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3}),\quad
s_{3}\in(a_{2},a_{3}),$

where ${\mathcal{E}}^{(3)}_{i,\mathbf{n},\mathbf{p}}$ is analytic on
$[a_{i-1},a_{i}]$ for $i=1,2$ while
${\mathcal{E}}^{(3)}_{3,\mathbf{n},\mathbf{p}}$ is analytic on
$[a_{2},a_{3})$. As in [13, Section X] we define the internal 5-cyclidic
harmonic of the third kind by

(7.1)
$G^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}E^{(3)}_{1,\mathbf{n},\mathbf{p}}(s_{1})E^{(3)}_{2,\mathbf{n},\mathbf{p}}(s_{2})E^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3}),\quad\mathbf{r}\in
R.$

Let $P^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})$ be the solution to (2.10) (with
$\lambda_{j}=\lambda^{(3)}_{j,\mathbf{n},\mathbf{p}}$) on $(a_{2},a_{3})$
belonging to the exponent $0$ at $s=a_{3}$ and uniquely determined by the
condition $P^{(3)}_{3,\mathbf{n},\mathbf{p}}(a_{3})=1$. We write

$P^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})=(s_{3}-a_{2})^{p_{2}/2}{\mathcal{P}}^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3}),\quad
s_{3}\in(a_{2},a_{3}),$

where ${\mathcal{P}}^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})$ is analytic on
$(a_{2},a_{3}]$. Then we define

(7.2)
$I^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}\prod_{j=0}^{2}(\chi_{j}(\mathbf{r}))^{p_{j}}\prod_{i=1}^{2}{\mathcal{E}}^{(3)}_{i,\mathbf{n},\mathbf{p}}(s_{i}){\mathcal{P}}^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})\quad\text{if
$s_{3}\neq a_{2}$}.$

The condition $s_{3}\neq a_{2}$ is equivalent to
$\mathbf{r}\in\mathbb{R}^{3}\setminus(K_{2}\cup M_{2})$; see Figure 2.

Similarly, let $Q^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})$ be the solution to
(2.10) (with $\lambda_{j}=\lambda^{(3)}_{j,\mathbf{n},\mathbf{p}}$) on
$(a_{2},a_{3})$ belonging to the exponent $\frac{1}{2}$ at $s=a_{3}$ and
uniquely determined by the condition $\lim_{s_{3}\to
a_{3}^{-}}\omega(s_{3})\frac{d}{ds_{3}}Q^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})=1$.
We write

$Q^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})=(a_{3}-s_{3})^{1/2}(s_{3}-a_{2})^{p_{2}/2}{\mathcal{Q}}^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3}),\quad
s_{3}\in(a_{2},a_{3}),$

where ${\mathcal{Q}}^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})$ is analytic on
$(a_{2},a_{3}]$. Then we define

(7.3)
$J^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=(\|\mathbf{r}\|^{2}+1)^{-1/2}\chi_{3}(\mathbf{r})\prod_{j=0}^{2}(\chi_{j}(\mathbf{r}))^{p_{j}}\prod_{i=1}^{2}{\mathcal{E}}^{(3)}_{i,\mathbf{n},\mathbf{p}}(s_{i}){\mathcal{Q}}^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})\quad\text{if
$s_{3}\neq a_{2}$}.$

######  Lemma 7.1.

The functions $I^{(3)}_{\mathbf{n},\mathbf{p}}$ and
$J^{(3)}_{\mathbf{n},\mathbf{p}}$ are harmonic on
$\mathbb{R}^{3}\setminus(K_{2}\cup M_{2})$. They have the symmetries

(7.4) $\displaystyle I^{(3)}_{\mathbf{n},\mathbf{p}}(\sigma_{0}(\mathbf{r}))$
$\displaystyle=$
$\displaystyle(-1)^{p_{0}}\|\mathbf{r}\|I^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),$
(7.5) $\displaystyle I^{(3)}_{\mathbf{n},\mathbf{p}}(\sigma_{j}(\mathbf{r}))$
$\displaystyle=$
$\displaystyle(-1)^{p_{j}}I^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),\quad
j=1,2,$ (7.6) $\displaystyle
I^{(3)}_{\mathbf{n},\mathbf{p}}(\sigma_{3}(\mathbf{r}))$ $\displaystyle=$
$\displaystyle I^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),$ (7.7)
$\displaystyle J^{(3)}_{\mathbf{n},\mathbf{p}}(\sigma_{0}(\mathbf{r}))$
$\displaystyle=$
$\displaystyle(-1)^{p_{0}}\|\mathbf{r}\|J^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),$
(7.8) $\displaystyle J^{(3)}_{\mathbf{n},\mathbf{p}}(\sigma_{j}(\mathbf{r}))$
$\displaystyle=$
$\displaystyle(-1)^{p_{j}}J^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),\quad
j=1,2,$ (7.9) $\displaystyle
J^{(3)}_{\mathbf{n},\mathbf{p}}(\sigma_{3}(\mathbf{r}))$ $\displaystyle=-$
$\displaystyle J^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}).$

Since $P^{(3)}_{3,\mathbf{n},\mathbf{p}},Q^{(3)}_{3,\mathbf{n},\mathbf{p}}$
form a fundamental system of solutions to (2.10) (with
$\lambda_{j}=\lambda^{(3)}_{j,\mathbf{n},\mathbf{p}}$) on $(a_{2},a_{3})$,
there are (nonzero) scalars $\alpha^{(3)}_{\mathbf{n},\mathbf{p}}$,
$\beta^{(3)}_{\mathbf{n},\mathbf{p}}$ such that

$E^{(3)}_{3,\mathbf{n},\mathbf{p}}=\alpha^{(3)}_{\mathbf{n},\mathbf{p}}P^{(3)}_{3,\mathbf{n},\mathbf{p}}+\beta^{(3)}_{\mathbf{n},\mathbf{p}}Q^{(3)}_{3,\mathbf{n},\mathbf{p}}.$

This leads to the global definition of internal 5-cyclidic harmonics of the
third kind

(7.10)
$G^{(3)}_{\mathbf{n},\mathbf{p}}:=\alpha^{(3)}_{\mathbf{n},\mathbf{p}}I^{(3)}_{\mathbf{n},\mathbf{p}}+\beta^{(3)}_{\mathbf{n},\mathbf{p}}J^{(3)}_{\mathbf{n},\mathbf{p}}.$

If $z>0$, we can write $G^{(3)}_{\mathbf{n},\mathbf{p}}$ as follows

(7.11)
$G^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})=(\|\mathbf{r}\|^{2}+1)^{-1/2}\prod_{j=0}^{2}(\chi_{j}(\mathbf{r}))^{p_{j}}\prod_{i=1}^{3}{\mathcal{E}}^{(3)}_{i,\mathbf{n},\mathbf{p}}(s_{i}).$

######  Theorem 7.2.

Let $\mathbf{n}\in\mathbb{N}_{0}^{2}$ and
$\mathbf{p}=(p_{0},p_{1},p_{2})\in\\{0,1\\}^{3}$. Then
$G^{(3)}_{\mathbf{n},\mathbf{p}}$ extends continuously to a harmonic function
on $\mathbb{R}^{3}\setminus M_{2}$. Moreover

(7.12) $\displaystyle G^{(3)}_{\mathbf{n},\mathbf{p}}(\sigma_{0}(\mathbf{r}))$
$\displaystyle=$
$\displaystyle(-1)^{p_{0}}\|\mathbf{r}\|G^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),$
(7.13) $\displaystyle G^{(3)}_{\mathbf{n},\mathbf{p}}(\sigma_{j}(\mathbf{r}))$
$\displaystyle=$
$\displaystyle(-1)^{p_{j}}G^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}),\quad
j=1,2.$

We introduce another solution of (2.10) by

(7.14)
$F^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})=\gamma^{(3)}_{\mathbf{n},\mathbf{p}}\left(\alpha^{(3)}_{\mathbf{n},\mathbf{p}}P^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})-\beta^{(3)}_{\mathbf{n},\mathbf{p}}Q^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})\right),\quad
s_{3}\in(a_{2},a_{3}).$

We determine $\gamma^{(3)}_{\mathbf{n},\mathbf{p}}$ from the Wronskian

(7.15)
$\omega(s_{3})\left(E^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})\frac{d}{ds_{3}}F^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})-F^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})\frac{d}{ds_{3}}E^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})\right)=1,$

which is equivalent to

$\gamma^{(3)}_{\mathbf{n},\mathbf{p}}=\frac{-1}{2\alpha^{(3)}_{\mathbf{n},\mathbf{p}}\beta^{(3)}_{\mathbf{n},\mathbf{p}}}.$

We define external 5-cyclidic harmonics of the third kind by

(7.16)
$H^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}):=\gamma^{(3)}_{\mathbf{n},\mathbf{p}}G^{(3)}_{\mathbf{n},\mathbf{p}}(\sigma_{3}(\mathbf{r}))\quad\text{for
$\mathbf{r}\in\mathbb{R}^{3}\setminus K_{2}$}.$

In particular, we have

(7.17)
$H^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})=(\|\mathbf{r}\|^{2}+1)^{-1/2}E^{(3)}_{1,\mathbf{n},\mathbf{p}}(s_{1})E^{(3)}_{2,\mathbf{n},\mathbf{p}}(s_{2})F^{(3)}_{3,\mathbf{n},\mathbf{p}}(s_{3})\quad\text{for
$\mathbf{r}\in R$.}$

######  Theorem 7.3.

Let $\mathbf{n}\in\mathbb{N}_{0}^{2}$ and
$\mathbf{p}=(p_{0},p_{1},p_{2})\in\\{0,1\\}^{3}$. Then
$H^{(3)}_{\mathbf{n},\mathbf{p}}$ is harmonic on $\mathbb{R}^{3}\setminus
K_{2}$. The functions $H^{(3)}_{\mathbf{n},\mathbf{p}}$ share the symmetries
(7.12), (7.13) with $G^{(3)}_{\mathbf{n},\mathbf{p}}$. Moreover,

(7.18)
$H^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})=O(\|\mathbf{r}\|^{-1})\quad\text{as
$\|\mathbf{r}\|\to\infty$},$

and

(7.19) $\|\nabla
H^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})\|=O(\|\mathbf{r}\|^{-2})\quad\text{as
$\|\mathbf{r}\|\to\infty$}.$

##  8\. Expansion of the reciprocal distance in 5-cyclidic harmonics of third
kind

For fixed $s\in(a_{2},a_{3})$ the coordinate surface (2.1) consists of two
closed surfaces of genus $0$. One lies in the half-space $z>0$ and the other
one is obtained from it by reflection at the plane $z=0$. We consider the
region interior to the coordinate surface $s=d_{3}$ which lies in the half-
space $\\{\mathbf{r}:z>0\\}$:

(8.1) $D_{3}:=\\{\mathbf{r}\in\mathbb{R}^{3}:z>0,s_{3}<d_{3}\\}.$

######  Theorem 8.1.

Let $d_{3}\in(a_{2},a_{3})$, $\mathbf{n}\in\mathbb{N}_{0}^{2}$,
$\mathbf{p}\in\\{0,1\\}^{3}$. Then

(8.2)
$H^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}^{\prime})=\frac{1}{4\pi\omega(d_{3})\\{E^{(3)}_{3,\mathbf{n},\mathbf{p}}(d_{3})\\}^{2}}\int_{\partial
D_{3}}\frac{G^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})}{h_{3}(\mathbf{r})\|\mathbf{r}-\mathbf{r}^{\prime}\|}\,dS(\mathbf{r})$

for all $\mathbf{r}^{\prime}\in\mathbb{R}^{3}\setminus\bar{D}_{3}$. The scale
factor $h_{3}$ is given by

(8.3)
$16\\{h_{3}(\mathbf{r})\\}^{2}=\frac{(\|\mathbf{r}\|^{2}-1)^{2}}{(d_{3}-a_{0})^{2}}+\frac{4x^{2}}{(d_{3}-a_{1})^{2}}+\frac{4y^{2}}{(d_{3}-a_{2})^{2}}+\frac{4z^{2}}{(d_{3}-a_{3})^{2}}.$

We obtain the expansion of the reciprocal distance in 5-cyclidic harmonics of
the third kind.

######  Theorem 8.2.

Let
$\mathbf{r}=(x,y,z),\mathbf{r}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})\in\mathbb{R}^{3}$
with 5-cyclidic coordinates $s_{3},s_{3}^{\prime}$, respectively. If either
(a) $z,z^{\prime}\geq 0$, $s_{3}<s_{3}^{\prime}$, or (b) $z^{\prime}<0<z$, or
(c) $z,z^{\prime}\leq 0$, $s_{3}^{\prime}<s_{3}$, then

(8.4)
$\frac{1}{\|\mathbf{r}-\mathbf{r}^{\prime}\|}=2\pi\sum_{\mathbf{n}\in\mathbb{N}_{0}^{2}}\sum_{\mathbf{p}\in\\{0,1\\}^{3}}G^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r})H^{(3)}_{\mathbf{n},\mathbf{p}}(\mathbf{r}^{\prime}).$

## References

  * [1] D. H. Armitage and S. J. Gardiner.  Classical potential theory.  Springer Monographs in Mathematics. Springer-Verlag London Ltd., London, 2001.
  * [2] J. Blimke, J. Myklebust, H. Volkmer, and S. Merrill.  Four-shell ellipsoidal model employing multipole expansion in ellipsoidal coordinates.  Medical & Biological Engineering & Computing, 46(9):859–869, Sep 2008.
  * [3] M. Bôcher.  Ueber die Reihenentwickelungen der Potentialtheorie.  B. G. Teubner, Leipzig, 1894.
  * [4] C. P. Boyer, E. G. Kalnins, and W. Miller, Jr.  Symmetry and separation of variables for the Helmholtz and Laplace equations.  Nagoya Mathematical Journal, 60:35–80, 1976.
  * [5] H. S. Cohl.  Erratum: “Developments in determining the gravitational potential using toroidal functions”.  Astronomische Nachrichten, 333(8):784–785, 2012.
  * [6] H. S. Cohl.  Fourier, Gegenbauer and Jacobi expansions for a power-law fundamental solution of the polyharmonic equation and polyspherical addition theorems.  Symmetry, Integrability and Geometry: Methods and Applications, 9(042):26, 2013.
  * [7] H. S. Cohl and D. E. Dominici.  Generalized Heine’s identity for complex Fourier series of binomials.  Proceedings of the Royal Society A, 467:333–345, 2011.
  * [8] H. S. Cohl and E. G. Kalnins.  Fourier and Gegenbauer expansions for a fundamental solution of the Laplacian in the hyperboloid model of hyperbolic geometry.  Journal of Physics A: Mathematical and Theoretical, 45(14):145206, 2012.
  * [9] H. S. Cohl, A. R. P. Rau, J. E. Tohline, D. A. Browne, J. E. Cazes, and E. I. Barnes.  Useful alternative to the multipole expansion of $1/r$ potentials.  Physical Review A: Atomic and Molecular Physics and Dynamics, 64(5):052509, Oct 2001.
  * [10] H. S. Cohl, J. E. Tohline, A. R. P. Rau, and H. M. Srivastava.  Developments in determining the gravitational potential using toroidal functions.  Astronomische Nachrichten, 321(5/6):363–372, 2000.
  * [11] H. S. Cohl and H. Volkmer.  Eigenfunction expansions for a fundamental solution of Laplace’s equation on $\mathbb{R}^{3}$ in parabolic and elliptic cylinder coordinates.  Journal of Physics A: Mathematical and Theoretical, 45(35):355204, 2012.
  * [12] H. S. Cohl and H. Volkmer.  Publisher’s Note: “Separation of variables in an asymmetric cyclidic coordinate system”.  Journal of Mathematical Physics, 54(7):079904, 2013.
  * [13] H. S. Cohl and H. Volkmer.  Separation of variables in an asymmetric cyclidic coordinate system.  Journal of Mathematical Physics, 54(6):063513, 2013.
  * [14] L. P. Eisenhart.  Separable systems of Stackel.  Annals of Mathematics. Second Series, 35(2):284–305, 1934.
  * [15] E. Heine.  Handbuch der Kugelfunctionen, Theorie und Anwendungen (volume 2).  Druck und Verlag von G. Reimer, Berlin, 1881.
  * [16] E. W. Hobson.  The theory of spherical and ellipsoidal harmonics.  Chelsea Publishing Company, New York, 1955.
  * [17] O. D. Kellogg.  Foundations of potential theory.  Reprint from the first edition of 1929. Die Grundlehren der Mathematischen Wissenschaften, Band 31. Springer-Verlag, Berlin, 1967.
  * [18] T. M. MacRobert.  Spherical Harmonics. An Elementary Treatise on Harmonic Functions with Applications.  Methuen & Co. Ltd., London, second edition, 1947.
  * [19] R. C. McOwen.  Partial Differential Equations: Methods and Applications.  Prentice Hall, Upper Saddle River, New Jersey, 1996.
  * [20] W. Miller, Jr.  Symmetry and separation of variables.  Addison-Wesley Publishing Co., Reading, Mass.-London-Amsterdam, 1977.  With a foreword by Richard Askey, Encyclopedia of Mathematics and its Applications, Vol. 4.
  * [21] P. M. Morse and H. Feshbach.  Methods of theoretical physics. 2 volumes.  McGraw-Hill Book Co., Inc., New York, 1953.
  * [22] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors.  NIST handbook of mathematical functions.  Cambridge University Press, Cambridge, 2010.
  * [23] G. N. Watson.  A treatise on the theory of Bessel functions.  Cambridge Mathematical Library. Cambridge University Press, Cambridge, second edition, 1944.
  * [24] Z. Y. Wen and J. Avery.  Some properties of hyperspherical harmonics.  Journal of Mathematical Physics, 26(3):396–403, 1985.