1 Introduction

The Dirichlet-to-Neumann map for the Laplacian is important in many areas of basic analysis (elliptic boundary value problems [10, 23, 27, 32], inverse problems [5, 20]), as well as physics (e.g., fluid mechanics [9], electromagnetic theory [6], electrical impedance tomography [19, 22, 35], electrical transmission [11]). The map relates the values of a function \(f:\partial \Omega \rightarrow {\mathbb {R}}\) (Dirichlet datum) on the boundary \(\partial \Omega \) (\(\Omega \subseteq {\mathbb {R}}^3\)) to the normal derivative (Neumann datum) of the unique harmonic function \(u:\Omega \rightarrow {\mathbb {R}}\) whose boundary values are given by f. The Neumann problem is to recover f (or u) from the normal derivative.

When \(\Omega \) is a sphere, solving the Dirichlet problem amounts to copying the boundary function’s Fourier coefficients in terms of spherical harmonics to obtain the coefficients of the interior harmonic function in terms of solid spherical harmonics. Thus it would make no sense to apply a general-purpose method (such as finite elements) for the Dirichlet problem on the sphere; the same holds for the Neumann problem (cf. Remark 3.7 below).

When \(\Omega \) is a torus, harmonic functions defined on \(\Omega \) are naturally expressed as series based on a doubly-indexed collection of toroidal harmonics (involving associated Legendre functions of the first and second kinds of half-integer degree) [17], which are orthogonal with respect to a certain weighted \(L^2\)-inner product. As for the sphere, the coefficients of this expression provide a solution of the Dirichlet problem for the Laplacian quite directly. However, unlike the case for a sphere, the Dirichlet-to-Neumann mapping for a torus turns out to be far more complicated, and the solution of the Neumann problem involves solving an underdetermined infinite upper-triangular system of linear equations. We express the well-known necessary and sufficient condition for the solvability of the Neumann problem (compatibility condition), as well as the normalization condition, in terms of the Fourier coefficients. The solution of the Neumann problem turns out to involve a special twist in that the free (i.e., undetermined) parameter which appears in the linear system cannot be found algebraically, as far as we know (it is not connected to the normalization condition). However, it can be expressed as a limit of easily calculated algebraic expressions. The solution procedure is illustrated through numerical examples.

A similar phenomenon was found by Hicks [18] in calculating the velocity potential for the fluid motion of a torus moving parallel to a fluid at rest at infinity, in which the unique free parameter leading to convergence of the series was determined in terms of another infinite series. Love [26] approached an inhomogeneous Neumann problem by means of Green’s functions for second-order difference equations. In both cases, the construction was based on an expansion of the potential functions in terms of toroidal harmonics.

Of course, one can solve the numerical problem by approximating the normal derivative numerically, but this approach provides little insight into the structure of the solutions (cf. Proposition 5.4 and Corollary 5.5 below.) As far as we know, our explicit formula for the normal derivative in toroidal coordinates has not appeared previously. We have seen applications (e.g., [6]) where the potential in a toroidal conductor is modeled in spherical coordinates, which are not ideally suited for such problems. A method for solving the Neumann problem in a more general class of regions termed “toroidal” appeared in [28], but the notion of “toroidal coordinates” which we use here is not applicable to that method.

The paper is organized as follows. Some properties of toroidal harmonic functions are summarized in Sect. 2. In Sect. 3, a computation of the toroidal Neumann derivative leads to a formula for the Dirichlet-to-Neumann mapping. The expansion coefficients of the normal derivative are linear expressions in the Fourier coefficients on the surface of the torus. While the mapping exists in the context of the appropriate Sobolev function spaces associated with the Laplace equation and the boundary data, as well as the trace operator, we only justify the derivation of our formulas for the normal derivative under stronger smoothness assumptions. Sections 4–6 present the algorithm and numerical examples to show the accuracy of the procedure. The computational complexity is shown to be O(NM) in a specific sense, far faster than a general-purpose method such as finite elements. Section 7 explains how to solve the Neumann problem in a toroidal shell using the results for the interior and exterior.

2 Toroidal coordinates and toroidal harmonics

We begin with the notation and several well-known facts to be used throughout the paper.

2.1 Normal derivatives

Toroidal coordinates \((\eta ,\theta ,\varphi )\) for a point \(x=(x_0,x_1,x_2)\) in three-dimensional Euclidean space are defined [17, 30] by the relations

$$\begin{aligned} x_0&= \frac{\sin \theta }{\cosh \eta -\cos \theta }, \ \ x_1 = \frac{\sinh \eta \, \cos \varphi }{\cosh \eta -\cos \theta }, \ \ x_2 = \frac{\sinh \eta \, \sin \varphi }{\cosh \eta -\cos \theta } \end{aligned}$$
(1)

in the range \(\eta \in (0,\infty )\), \(\theta \in [-\pi ,\pi ]\), \(\varphi \in (-\pi ,\pi ]\). For any fixed \(\eta _0>0\), these coordinates define the interior domain \(\Omega _{\eta _0}= \{x:\ \eta _0<\eta \le \infty \} \cup S^1\) where the limiting value \(\eta =\infty \) corresponds respectively to the singular subset \(S^1=\{x\in {\mathbb {R}}^3:\ x_0=0,\ x_1^2+x_2^2=1\}\). Any open solid torus in \({\mathbb {R}}^3\) can be rotated, shifted, and rescaled to a torus of the form \(\Omega _{\eta _0}\).

The unit normal vector \({\textbf{n}}=(x_\theta \times x_\Phi )/\vert x_\theta \times x_\varphi \vert \) on \({\partial \Omega _{\eta _0}}\) where \(x_\theta \)= \(\partial x/\partial \theta ~\text {and}\)  \(x_\phi \)= \(\partial x/\partial \varphi \), is equal to

$$\begin{aligned} {{\textbf{n}}} =&\, \frac{-1}{(\cosh \eta _{0}-\cos \theta )} \\&\big ( \sinh \eta _{0}\sin \theta ,\ \cos \varphi (\cosh \eta _{0}\cos \theta -1),\ \sin \varphi (\cosh \eta _{0}\cos \theta -1) \big ) .\nonumber \end{aligned}$$
(2)

Since toroidal coordinates form an orthogonal coordinate system, we have also \({{\textbf{n}}} = x_\eta /\vert x_\eta \vert \), and since \(\eta \) increases as points move towards the interior of \(\Omega _{\eta _0}\), we see that \({{\textbf{n}}}\) is in fact the inward pointing normal vector on \({\partial \Omega _{\eta _0}}\). The normal derivative of a function f defined in a neighborhood V of a point \(x\in {\partial \Omega _{\eta _0}}\) (or a half-neighborhood \(V\cap {\overline{\Omega }}_{\eta _0}\)) is by definition \({{\,\mathrm{{\text{ nor }}}\,}}f(x) = \frac{d}{dt}f(x+t\,{{\textbf{n}}}) \vert _{t=0^+}\) \(= ({{\,\mathrm{{\text{ grad }}}\,}}f(x))\!\!\cdot {{\textbf{n}}}\), and by orthogonality of the coordinate system the normal derivative is also equal to the following, which is often more convenient for calculations:

$$\begin{aligned} {{\,\mathrm{{\text{ nor }}}\,}}f (x)= \frac{1}{\vert x_\eta \vert }\,\frac{d}{d\eta }f(\eta ,\theta ,\varphi )\bigg \vert _{\eta =\eta _0^+} . \end{aligned}$$
(3)

2.2 Toroidal harmonics

The associated Legendre functions (Ferrer’s functions) of the first and second kinds for \(t>1\) are defined for integer values of \(n,m\ge 0\), respectively, as

$$\begin{aligned} P_n^m(t)&=(t^{2}-1)^{m/2}\,\frac{d^{m}P_{n}(t)}{dt^{m}}, \nonumber \\ Q_n^m(t)&= (t^{2}-1)^{m/2}\,\frac{d^{m}Q_{n}(t)}{dt^{m}}, \end{aligned}$$
(4)

where \(P_n(t),Q_n(t)\) denote the classical Legendre polynomials of degree n [1, 2, 4, 8, 12, 14, 17, 21, 34, 36]. Extended analytically in the complex plane, these functions are branched at \(t=\pm 1\). However, they are entire functions of n and m regarded as complex variables [17]. In this sense, (4) can be taken as a definition of the associated Legendre functions for half-integer values of n as we will need here.

We will abbreviate \(\Phi {}^{\hspace{.1ex}\nu }_{n}(\theta ) = \cos n\theta \) for \(\nu =1\), and \(\Phi {}^{\hspace{.1ex}\nu }_{n}(\theta ) = \sin n\theta \) for \(\nu =-1\). The (normalized) interior toroidal harmonic functions are

$$\begin{aligned} I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}(x)=I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}[\eta _{0}](x) = \sqrt{\cosh \eta -\cos \theta }\,\frac{Q{}^{\hspace{.1ex}m}_{n-1/2}(\cosh \eta )}{Q{}^{\hspace{.1ex}m}_{n-1/2}(\cosh \eta _{0})} \, \Phi {}^{\hspace{.1ex}\nu }_{n}(\theta ) \, \Phi {}^{\hspace{.1ex}\mu }_{m}(\varphi ) \end{aligned}$$
(5)

for integers n, m with

$$\begin{aligned} n\ge 0,\quad m\ge 0,\quad \nu \in \{-1,1\} ,\quad \mu \in \{-1,1\} . \end{aligned}$$
(6)

A derivation of the Laplacian equation in toroidal coordinates and the verification that \(I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\) is harmonic in \(\Omega _{\eta _0}\) can be found in [17, p. 434]; standard results on removable sets for harmonic functions [3] cover the behavior at the singular set.

One may define a weighted \(L^2\) inner product on real-valued functions by

$$\begin{aligned} \langle f,g \rangle _{\eta _0} = \iiint _{\Omega _{\eta _0}} f g\, w \,dV \end{aligned}$$
(7)

with the weight function

$$\begin{aligned} w(\eta ,\theta ,\varphi ) = \frac{(\cosh \eta -\cos \theta )^2}{\sinh \eta }. \end{aligned}$$
(8)

One deduces the following using the fact that the collection of products \(\{\Phi _{n}^{v}\Phi _{m}^{\mu }\}\) is a complete orthogonal set in \(L^{2}(\Pi , \Pi ^{2})\) (cf. [31, Sec. 25]).

Proposition 2.1

The interior toroidal harmonics \(\{I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\}\) (for \(n,m,\nu ,\mu \) as in (6)) form a complete orthogonal system in the harmonic functions in \(L^2(\Omega _{\eta _0},w)\). Their norms are

$$\begin{aligned} \Vert I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\Vert ^2_{\eta _0} = \varepsilon _n\varepsilon _m\int \limits _{n_0}^\infty \bigg (\frac{Q{}^{\hspace{.1ex}m}_{n-1/2}(\cosh \eta )}{Q{}^{\hspace{.1ex}m}_{n-1/2}(\cosh \eta _{0})}\bigg )^2\,d\eta , \end{aligned}$$

where \(\varepsilon _n= 1+\delta _{n,0}\) and \(\delta _{n,m}\) is the Kronecker delta function. Further, the restrictions to the boundary \(\{\left. I{}^{\hspace{.1ex}\nu ,\mu }_{n,m} \right| _{{\partial \Omega _{\eta _0}}}\}\) are complete in \(L^2({\partial \Omega _{\eta _0}})\) and \(L^2({\partial \Omega _{\eta _0}},w)\).

3 Expression of the Dirichlet-to-Neumann map in toroidal coordinates

3.1 Background on the Neumann problem

Given a suitable \(f:{\partial \Omega _{\eta _0}}\rightarrow {\mathbb {R}}\), the Dirichlet-to-Neumann mapping is described by finding the unique harmonic function \(u\in {{\,\mathrm{{\text{ Har }}}\,}}(\Omega _{\eta _0})\) with boundary values \(f=u\vert _{\partial \Omega _{\eta _0}}\), and then taking the normal derivative \(h={{\,\mathrm{{\text{ nor }}}\,}}u\). The mapping is thus \(\Lambda f = h\).

A common setting [25] is for u to be in the Sobolev space \(H^1(\Omega _{\eta _0})\) and f in the boundary space \(H^{1/2}({\partial \Omega _{\eta _0}})\). Roughly speaking, \(H^1(\Omega _{\eta _0})\) consists of \(L^2\) functions having \(L^2\) derivatives, and \(H^{1/2}({\partial \Omega _{\eta _0}})=H^1(\Omega _{\eta _0})/H^1_0(\Omega )\) is identified with the space of boundary values of elements of \(H^1(\Omega _{\eta _0})\), where \(H^1_0(\Omega )\) denotes the closure of the subspace of functions of compact support. The Neumann problem is to find a function f such that \(\Lambda f=h\), for a given function h defined on \(\partial \Omega _{\eta _0}\). The solution of the Neumann problem, in general, is guaranteed by the following result [10, 13, 29], valid for domains \(\Omega \) with sufficiently smooth boundary.

Proposition 3.1

Let \(h\in H^{-1/2}(\partial \Omega )\) satisfy the compatibility condition

$$\begin{aligned} \int \limits _{\partial \Omega } h\, dS = 0. \end{aligned}$$
(9)

Then there exists an \(f\in H^{1/2}(\partial \Omega )\) such that \(\Lambda f = h\). This solution is unique up to an additive constant. If \(h\in L^2({\partial \Omega _{\eta _0}})\), then \(f\in L^2({\partial \Omega _{\eta _0}})\). If h is continuous, then f is continuous.

The solution f can made unique by requiring the normalization condition

$$\begin{aligned} \int \limits _{\partial \Omega _{\eta _0}}f\, dS = c \end{aligned}$$
(10)

for a chosen constant c.

3.2 Dirichlet-to-Neumann mapping in toroidal coordinates

The coordinate expression of the Dirichlet-to-Neumann mapping will be based on the following calculation.

Lemma 3.2

The normal derivatives of the interior toroidal harmonics are given by the formula

$$\begin{aligned} {{\,\mathrm{{\text{ nor }}}\,}}I {}^{\hspace{.1ex}\nu ,\mu }_{n,m}&= \Bigg (\Phi _{n-1}^{\nu }(\theta ) \bigg ( (1+2n)\cosh \eta _{0} -\Big (2(n-m)+1\Big ) \frac{Q_{n+1/2}^{m}(\cosh \eta _{0})}{Q_{n-1/2}^{m}(\cosh \eta _{0})}\bigg ) \\&\quad +\Phi _{n}^{\nu }(\theta )\bigg ( -2\Big ( (2n\cosh ^{2}\eta _{0}+1) \\&\quad -\Big (2(n-m)+1\Big )\cosh \eta _{0} \, \frac{Q_{n+1/2}^{m}(\cosh \eta _{0})}{Q_{n-1/2}^{m}(\cosh \eta _{0})}\Big )\bigg )\\&\quad +\Phi _{n+1}^{\nu }(\theta )\bigg ((1+2n)\cosh \eta _{0} -\Big (2(n-m)+1\Big ) \frac{Q_{n+1/2}^{m}(\cosh \eta _{0})}{Q_{n-1/2}^{m}(\cosh \eta _{0})}\bigg )\Bigg ) \\&\quad \times \frac{(\cosh \eta _{0}-\cos \theta )^{1/2}}{4\sinh \eta _{0}} \, \Phi _{m}^{\mu }(\varphi ) . \end{aligned}$$

Proof

Apply (3)–(5) and obtain a formula for \({{\,\mathrm{{\text{ nor }}}\,}}I {}^{\hspace{.1ex}\nu ,\mu }_{n,m}\), Now eliminate all occurrences of the derivative \((Q^m_{n+1})^{\prime }\) by means of the recurrence formula [4, pp. 161–162]

$$\begin{aligned} \quad \sinh ^2\!\eta \, (Q^m_{n+1})^{\prime }(\cosh \eta ) = (n+m+1) Q^m_{n}(\cosh \eta ) - (n+1)\cosh \eta \, Q^m_{n+1}(\cosh \eta ). \qquad \qquad \end{aligned}$$

Finally, regroup the terms as multiples of the trigonometric expressions of the form \(\Phi _n^{\nu }\). \(\square \)

The boundary function f may be represented conveniently for our purposes via \( f(\theta ,\varphi )/ \sqrt{\cosh \eta _0-\cos \theta }= \sum _{n,m,\nu ,\mu } a{}^{\hspace{.1ex}\nu ,\mu }_{n,m} \Phi _n^\nu (\theta ) \Phi _m^\mu (\varphi ),\) i.e.,

$$\begin{aligned} f(\theta ,\varphi ) = \sum _{n,m,\nu ,\mu } a{}^{\hspace{.1ex}\nu ,\mu }_{n,m} I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}(\eta _0,\theta ,\varphi ) \end{aligned}$$
(11)

as per Proposition 2.1, with \(a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\in {\mathbb {R}}\). Unless otherwise specified, the indices of summation will always be as in (6) but excluding the cases of \((n,m,\nu ,\mu )\) being \((0,m,-1,\mu )\) or \((n,0,\nu ,-1)\), taking into account that \(\Phi {}^{\hspace{.1ex}-}_{0}=0\) identically. For convenience, we will often write superscripts as “\(+\)” in place of 1 and “−” in place of \(-1\).

With this notation, the solution of the Dirichlet problem on \(\Omega _{\eta _0}\) is trivial: \(f=u\vert _{\partial \Omega _{\eta _0}}\) where \(u\in {{\,\mathrm{{\text{ Har }}}\,}}(\Omega _{\eta _0})\) is given by

$$\begin{aligned} u=\sum _{n,m,\nu ,\mu } a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m} . \end{aligned}$$
(12)

Then by Lemma 3.2, we have \(h = \Lambda f\) is in turn given by

$$\begin{aligned} h = \sum _{n,m,\nu ,\mu } a{}^{\hspace{.1ex}\nu ,\mu }_{n,m} {{\,\mathrm{{\text{ nor }}}\,}}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m} , \end{aligned}$$
(13)

assuming that f is sufficiently well-behaved to justify the exchange of summation and differentiation. For example, since the trace operator \(\text{ tr }:H^1(\Omega _{\eta _0})\rightarrow H^{1/2}({\partial \Omega _{\eta _0}}) \) and \(\Lambda :H^{1/2}({\partial \Omega _{\eta _0}})\rightarrow H^{-1/2}({\partial \Omega _{\eta _0}})\) are continuous [33], it is easily seen that the sum \({{\,\mathrm{{\text{ nor }}}\,}}u = \sum a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}{{\,\mathrm{{\text{ nor }}}\,}}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\) converges in the dual space \(H^{-1/2}({\partial \Omega _{\eta _0}})\) of \(H^{1/2}({\partial \Omega _{\eta _0}})\). Alternatively, the convergence is valid under the assumption that the sum in (12) together with the sums of the partial derivatives of the terms converge uniformly on compact subsets of \(\Omega _{\eta _0}\).

We will show how \({{\,\mathrm{{\text{ nor }}}\,}}u\) can be expressed in terms of the coefficients of f and certain constants defined in terms of Legendre functions. We will use the abbreviations

$$\begin{aligned} t_0=\cosh \eta _0,\quad s_0=\sinh \eta _0,\quad q_{n,m} = Q_{n-1/2}^m(t_0). \end{aligned}$$
(14)

Definition 3.3

The toroidal Neumann constants, \(\rho _{n,m} = \rho _{n,m}(\eta _{0})\), \(\sigma _{n,m} = \sigma _{n,m}(\eta _{0})\), and \(\tau _{n,m}= \tau _{n,m}(\eta _{0})\) associated to \(\Omega _{\eta _0}\) are defined as follows:

$$\begin{aligned} \rho _{1,m}&= \frac{1}{2s_0}\left( t_0 + (2m-1)\frac{q_{1,m}}{q_{0,m}}\right) , \nonumber \\ \rho _{n,m}&= \frac{1}{4s_0}\Big ((2n-1) t_0+ (2(m-n)+1) \frac{q_{n,m}}{q_{n-1,m}} \Big ) \ \ (n\ge 2), \nonumber \\ \sigma _{n,m}&= \frac{-1}{2s_0}\Big ((2nt_{0}^{2}+1) + (2(m-n)-1)t_0 \frac{q_{n+1,m}}{q_{n,m}}\Big ) \ \ (n\ge 0) ,\nonumber \\ \tau _{n,m}&= \frac{1}{4s_0}\Big ((2n+3) t_0 + (2(m-n)-3) \frac{q_{n+2, m}}{q_{n+1,m}}\Big ) \ \ (n\ge 0), \end{aligned}$$
(15)

for all \(m\ge 0\).

We will use the following well known facts about the growth of Fourier coefficients and the convergence of Fourier series.

Proposition 3.4

([15, Corollary 3.3.10, Proposition 3.3.12]) If \(f(\theta ,\varphi )\) is of class \(C^r\), then

$$\begin{aligned} \vert a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\vert \le \frac{C}{(m+n+1)^r} \end{aligned}$$
(16)

for some constant \(C>0\). Conversely, if (16) holds and \(r\ge 2\), then f is of class \(C^{r-2}\).

The main result of this section is that the coefficients (15) are entries in an infinite upper-triangular matrix which defines the Dirichlet-to-Neumann mapping \(\Lambda \).

Theorem 3.5

For a fixed \(\eta _0\), let \(f:{\partial \Omega _{\eta _0}}\rightarrow {\mathbb {R}}\) given by (11) and suppose that the Fourier coefficients satisfy (16) for some constant C where \(r\ge 4\). Then the coefficients in the absolutely convergent series

$$\begin{aligned} h(\theta ,\varphi ) = \sqrt{\cosh \eta _{0}-\cos \theta } \sum _{n,m,\nu ,\mu } b{}^{\hspace{.1ex}\nu ,\mu }_{n,m} \Phi {}^{\hspace{.1ex}\nu }_{n}(\theta ) \Phi {}^{\hspace{.1ex}\mu }_{m}(\varphi ) \end{aligned}$$
(17)

for the Dirichlet-to-Neumann mapping \(h=\Lambda f\in C^{r-3}\) are given by

$$\begin{aligned} b{}^{\hspace{.1ex}+,\mu }_{0,m}&= \sigma _{0,m}\, a{}^{\hspace{.1ex}+,\mu }_{0,m} + \tau {}^{\hspace{.1ex}\,}_{0,m} a{}^{\hspace{.1ex}+,\mu }_{1,m}, \nonumber \\ b{}^{\hspace{.1ex}+,\mu }_{n,m}&= \rho _{n,m}\, a{}^{\hspace{.1ex}+,\mu }_{n-1,m} + \sigma _{n,m}\, a{}^{\hspace{.1ex}+,\mu }_{n,m} + \tau _{n,m}\, a{}^{\hspace{.1ex}+,\mu }_{n+1,m} \ \ (n\ge 1), \nonumber \\ b{}^{\hspace{.1ex}-,\mu }_{1,m}&= \sigma _{1,m} \, a{}^{\hspace{.1ex}-,\mu }_{1,m} + \tau _{1,m}\, a{}^{\hspace{.1ex}-,\mu }_{2,m}, \nonumber \\ b{}^{\hspace{.1ex}-,\mu }_{n,m}&= \rho _{n,m}\, a{}^{\hspace{.1ex}-,\mu }_{n-1,m} + \sigma _{n,m}\, a{}^{\hspace{.1ex}-,\mu }_{n,m} + \tau _{n,m}\, a{}^{\hspace{.1ex}-,\mu }_{n+1,m} \ \ (n\ge 2), \end{aligned}$$
(18)

for \(\mu =\pm 1\) and all \(m\ge 0\).

Proof

From (5) it follows that \(\big \vert I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}(\eta _0,\theta ,\varphi ) \big \vert \le \sqrt{t_0+1}\), so the expansion (11) converges absolutely. In [17, p. 305] it is shown that for fixed \(\eta _0\) and m,

$$\begin{aligned} Q_n^m(\cosh \eta _0) \ \sim \ (-1)^m \frac{\Gamma (n+m+1)}{\Gamma (n+1)} \big (\frac{\pi }{n}\big )^{1/2} \frac{e^{-(n+1/2)\eta _0}}{(2\sinh \eta _0)^{1/2}} \end{aligned}$$
(19)

as \(n\rightarrow \infty \), where \(\sim \) means that the ratio of the two expressions tends to 1. From this it is seen that

$$\begin{aligned} \lim _{n\rightarrow \infty } \frac{ q_{n,m} }{ q_{n-1,m} } = e^{-\eta _0} \end{aligned}$$
(20)

independently of the value of m. From (20) and Lemma 3.2 one verifies that

$$\begin{aligned} \big \vert {{\,\mathrm{{\text{ nor }}}\,}}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m} \big \vert \le C_1( n+m+1) \end{aligned}$$

for some constant \(C_1\) (which depends on \(\eta _0\)), and hence by (16),

$$\begin{aligned} \big \vert a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}{{\,\mathrm{{\text{ nor }}}\,}}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m} \big \vert \le \frac{CC_1}{(n+m+1)^{r-1}}. \end{aligned}$$

This is enough to guarantee that (13), a double series in \(\theta \) and \(\varphi \), also converges absolutely. (Our hypothesis on the value of r takes into account that \(\sum 1/(m\!+n\!+1)^2\!=\infty \).) This in turn permits us to substitute the formula for \({{\,\mathrm{{\text{ nor }}}\,}}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\) into (13) and then reindex \(\Phi _{n-1}^{\nu }(\theta )\) and \(\Phi _{n+1}^{\nu }(\theta )\) into \(\Phi _n^{\nu }(\theta )\) to obtain after a routine calculation that

$$\begin{aligned} h(\theta ,\varphi ) = \sqrt{t_{0}-\cos \theta }\sum _{m,n,\nu ,\mu } \Phi _{n}^{\nu }(\theta )\Phi _{m}^{\mu }(\varphi )\big (\rho _{n,m}a_{n-1,m}^{\nu ,\mu }+ \sigma _{n,m}a_{n,m}^{\nu ,\mu }+\tau _{n,m}a_{n+1,m}^{\nu ,\mu }\big ) , \end{aligned}$$
(21)

which is (17). By Definition 3.3 and with observation that \(t_0 - s_0 = e^{-\eta _0}\) it also follows from (20) that

$$\begin{aligned} \lim _{n\rightarrow \infty } \frac{\rho _{n,m}}{n} = \frac{1}{2}, \quad \lim _{n\rightarrow \infty } \frac{\sigma _{n,m}}{n} = -t_0, \quad \lim _{n\rightarrow \infty } \frac{\tau _{n,m}}{n} = \frac{1}{2}. \end{aligned}$$
(22)

Thus by (22), the Fourier coefficients \(b{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\) defining h are of order no greater than \(1/((m+n+1)^ {r-1})\), so by Proposition 3.4 we are done. \(\square \)

Remark 3.6

While the differentiability requirement in Theorem 3.5 may appear quite restrictive, we will see that numerically it is often not necessary. It is not difficult to verify that when h in Proposition 3.1 is real analytic, the solution f to the Neumann problem is also real analytic. Since we are mainly interested in the formal relationships among the coefficients we will not go deeper into relaxing the differentiability conditions.

Remark 3.7

When one applies the same principle underlying Theorem 3.5 to the Neumann problem on the sphere \(\vert x\vert <1\) in spherical coordinates \(y_0=\rho \cos \theta \), \(y_1=\rho \sin \theta \cos \phi \), \(y_2=\rho \sin \theta \sin \phi \), and uses the standard solid spherical harmonics \(Y_{n,m}^\pm =\rho ^n P_n^m(\cos \theta )\Phi _m^\pm (\varphi )\), \(0\le m\le n\) as the basis for the harmonic functions, it is natural to represent the Dirichlet and Neumann functions \(f(\theta ,\varphi )\) and \(h(\theta ,\varphi )\) with the well-known basis \(\{Y_{n,m}^\pm (1,\theta ,\varphi )\)} for \(L^2\) functions on the sphere [32]. This yields certain coefficients \(a_{n,m}^\pm \) and \(b_{n,m}^\pm \) for f and h respectively. Since the normal derivative in this situation is equal to the radial derivative, i.e. \({{\,\mathrm{{\text{ nor }}}\,}}u = (\partial u/\partial \rho )\vert _{\rho =1}\), one sees immediately that \(b_{n,m}^\pm =na_{n,m}^\pm \), so the spherical analogue of the system (18) is rather trivial. (In the text [24, p. 218], this approach is suggested in a remark after expressing the solution of the Neumann problem for a spheroid, but only for zonal harmonics: functions constant with respect to the angular coordinate, i.e., involving \(P_n\) but not general \(P_n^m\). Apart from this, we have not seen this type of solution presented in the literature.)

4 Algebraic solutions for the Fourier coefficients in the Neumann problem

We now focus on the problem of solving for the coefficients \(a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\) of f in terms of the coefficients \(b{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\) of h. The upper triangular system (18) can be solved almost trivially by back-substitution. It separates naturally into a separate system for each fixed set of parameters \((m,\nu ,\mu )\), corresponding to the separate modes of the problem. For \(\nu =+1\), the solution is

$$\begin{aligned}&a{}^{\hspace{.1ex}+,\mu }_{1,m} = \frac{1}{\tau _{0,m}} (b{}^{\hspace{.1ex}+,\mu }_{0,m} - \sigma _{0,m}\, a{}^{\hspace{.1ex}+,\mu }_{0,m} ) ,\\&a{}^{\hspace{.1ex}+,\mu }_{n+1,m} = \frac{1}{\tau _{n,m}}( b{}^{\hspace{.1ex}+,\mu }_{n,m} -\rho _{n,m}\, a{}^{\hspace{.1ex}+,\mu }_{n-1,m} - \sigma _{n,m}\, a{}^{\hspace{.1ex}+,\mu }_{n,m} ); \end{aligned}$$

while for \(\nu =-1\), the relations are

$$\begin{aligned}&a{}^{\hspace{.1ex}-,\mu }_{2,m} = \frac{1}{\tau _{1,m}} (b{}^{\hspace{.1ex}-,\mu }_{1,m} - \sigma _{1,m}\, a{}^{\hspace{.1ex}-,\mu }_{1,m} ) ,\\&a{}^{\hspace{.1ex}-,\mu }_{n+1,m} = \frac{1}{\tau _{n,m}}( b{}^{\hspace{.1ex}-,\mu }_{n,m} - \rho _{n,m}\, a{}^{\hspace{.1ex}-,\mu }_{n-1,m} - \sigma _{n,m}\, a{}^{\hspace{.1ex}-,\mu }_{n,m} ). \end{aligned}$$

Thus one takes the values \(a{}^{\hspace{.1ex}+,\mu }_{0,m}\) or \(a{}^{\hspace{.1ex}-,\mu }_{1,m}\), respectively, to be the only free parameters. Of course, the following fact is necessary for the justification of these statements.

Lemma 4.1

The values \(\tau _{n,m}\) are never zero.

Proof

In [17, p. 195] it is shown that

$$\begin{aligned} Q_n^m(t) = \frac{(-1)^m}{2^{n + 1}} \frac{(n + m)!}{n!} \, (t^2-1)^{m/2} \int \limits _{-1}^1 \frac{(1 - s^2)^n}{(t - s)^{n+m+1}} \,ds \end{aligned}$$

(in fact, this rather than (4) is used to give a definition of \(Q_n^m(t)\)). From this it follows that

$$\begin{aligned} (-1)^m Q_n^m(\cosh \eta ) > 0 \end{aligned}$$
(23)

for all n, m, and \(\eta \in {\mathbb {R}}^+\). From (15), we need to show that the value

$$\begin{aligned} 16 s_{0}\,q_{n+1,m}\,\tau _{n,m}={4}\big ((2n+3)t_{0}\,q_{n+1,m}+\big (2(n-m)-3\big )q_{n+2,m}\big ) \end{aligned}$$
(24)

does not vanish. Consider the recursion formulas from [4, pp. 161–162] and [17, p. 108]

$$\begin{aligned}&(n-m+1)Q_{n+1}^{m}(t) = (2n+1)\,t\, Q_{n}^{m}(t)-(n+m)Q_{n-1}^{m}(t), \\&(t^{2}-1)^{1/2} Q_{n+1}^{m}(t) = \frac{1}{2n+3}\big (Q_{n+2}^{m+1}(t)-Q_{n}^{m+1}(t)\big ); \end{aligned}$$

i.e.,

$$\begin{aligned} 2(n+1)\,t_{0}\,q_{n+1,m} =\left( n-m+\frac{3}{2}\right) q_{n+2,m}+\left( n+m+\frac{1}{2}\right) q_{n,m}, \end{aligned}$$
(25)
$$\begin{aligned} 2(n+1)\,s_{0}\,q_{n+1,m-1}=q_{n,m}-q_{n+2,m}. \end{aligned}$$
(26)

Applying (25) to (24), we find

$$\begin{aligned} 16(n+1)\,s_{0}\,q_{n+1,m}\,\tau _{n,m}&= (2n+3)\big (2(m+n)+1\big )q_{n,m} \\&\quad -(1+2n)\big (2n-2m+3\big )q_{n+2,m}. \end{aligned}$$

Now add and subtract \((2n+3)(2(m+n)+1)q_{n+2,m}\) and use (25), yielding

$$\begin{aligned} 16\,s_{0}\,q_{n+1,m}\,\tau _{n,m}=-2(2n+3)(2n+2m+1)q_{n+1,m-1}+8\,m\,q_{n+2,m}, \end{aligned}$$

which by (23) is never zero. \(\square \)

Definition 4.2

We will say that a collection of real numbers \(\{a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\}\) is an algebraic solution of the Neumann problem posed by \(\{b{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\}\) when all of the equations (18) are satisfied.

An algebraic solution will give rise to a solution of the Neumann problem if it defines a convergent series (11).

5 Convergent solutions for the toroidal Neumann problem

We observed above that the Dirichlet problem is effectively trivial when data is expressed in terms of the basis of toroidal harmonics. We will now see that for the Neumann problem, some perhaps surprising particularities occur in connection with the convergence properties of algebraic solutions.

The question is when the algebraic solution determined by choice of the initial coefficients \(a{}^{\hspace{.1ex}+,\mu }_{0,m}\) and \(a{}^{\hspace{.1ex}-,\mu }_{1,m}\) will provide a convergent series in (11). That is to say, one wants to know when the function

$$\begin{aligned} f{}^{\hspace{.1ex}\nu ,\mu }_{m} = \sum _{n=0}^\infty a{}^{\hspace{.1ex}\nu ,\mu }_{n,m} I{}^{\hspace{.1ex}\nu ,\mu }_{n,m} \end{aligned}$$
(27)

corresponding to the mode \((m,\nu ,\mu )\) will exist. The principal mode will be considered separately from the remaining modes.

5.1 Convergence of the algebraic solution for the principal mode

For the particular indices \((m,\nu ,\mu )=(0,+1,+1)\), the Neumann constants satisfy some special relations which we will need.

Lemma 5.1

We have

$$\begin{aligned}&\sigma _{0,0}\, q_{0,0} + 2\,\tau _{0,0}\, q_{1,0}=0, \\&\rho _{1,0}\, q_{0,0} + 2(\sigma _{1,0} \, q_{1,0}+ \tau _{1,0}\, q_{2,0})=0, \end{aligned}$$

while for \(n\ge 2\),

$$\begin{aligned} \rho _{n,0}\, q_{n-1,0}+ \sigma _{n,0} \, q_{n,0}+ \tau _{n,0} \, q_{n+1,0}= 0. \end{aligned}$$

Proof

Via the recursion formula (25) as well as the following one [4, pp. 161–162],

$$\begin{aligned} (n-m)t\, Q_{n}^{m}(t) = (t^2-1)^{1/2} \, Q_{n}^{m+1}(t)+(n+m)Q_{n-1}^{m}(t), \end{aligned}$$

i.e.,

$$\begin{aligned} \left( n-m-\frac{1}{2}\right) t_{0}q_{n,m}=s_{0}\,q_{n,m+1}+\left( n+m-\frac{1}{2}\right) \,q_{n-1,m} , \end{aligned}$$
(28)

direct computations employing (15) show that

$$\begin{aligned} \sigma _{0,0}\, q_{0,0} + 2\,\tau _{0,0}\, q_{1,0} =&\, \frac{-1}{2s_{0}}\big (q_{0,0}-t_{0} q_{1,0} \big )+\frac{3}{2s_{0}}\big (t_{0}q_{1,0}-q_{2,0}\big )\\ =&\, \frac{3}{8s_{0}}\big (q_{2,0}-q_{0,0}\big )-\frac{3}{8s_{0}}\big (q_{2,0}-q_{0,0}\big )=0. \end{aligned}$$

Similarly, referring again to Definition 3.3 we find where the last step follows from (25).

Finally,

$$\begin{aligned}&\rho _{n,0}\, q_{n-1,0}+ \sigma _{n,0} \, q_{n,0}+ \tau _{n,0} \, q_{n+1,0}\\&\quad = (2n-1)t_{0} q_{n-1,0} +(-2n+1-4nt_{0}^{2}-2)q_{n,0} \\&\qquad \ +(6n+5)t_{0}q_{n+1,0}-(2n+3)q_{n+2,0}\\&\quad =0. \end{aligned}$$

\(\square \)

When \(b{}^{\hspace{.1ex}+,+}_{n,0}=0\) for all n (i.e., when h vanishes in the principal mode), the corresponding equations (18) are linear homogeneous, and Lemma 5.1 implies that

$$\begin{aligned} \frac{a{}^{\hspace{.1ex}+,+}_{n,0}}{q_{n,0}} = 2 \frac{a{}^{\hspace{.1ex}+,+}_{0,0}}{q_{0,0}}\ \ (n\ge 1); \end{aligned}$$

i.e., \( a{}^{\hspace{.1ex}+,+}_{n,0}/q_{n,0} = \varepsilon _n a{}^{\hspace{.1ex}+,+}_{0,0} /q_{0,0} \) for all n. On comparing the formula of Proposition 5.2 with the exponent determined by \(\alpha =1/2\), one sees that the resulting solution

$$\begin{aligned} f_{0}^{+,+}(\theta ,\varphi ) = \frac{a{}^{\hspace{.1ex}+,+}_{0,0}}{q_{0,0}} \sum _{n=0}^\infty \varepsilon _n q_{n,0}\cos n\theta \,\sqrt{t_0-\cos \theta } \end{aligned}$$

is indeed equal to a constant function on \({\partial \Omega _{\eta _0}}\) (as must be the case according to Proposition 3.1) with value \((\pi /\sqrt{2})(a{}^{\hspace{.1ex}+,+}_{0,0}/q_{0,0})\). The solution of the Dirichlet problem in \(\Omega _{\eta _0}\) is the same constant. Thus, for every \(a{}^{\hspace{.1ex}+,+}_{0,0}\in {\mathbb {R}}\), the algebraic solution gives a convergent series \(\sum _{n=0}^\infty a{}^{\hspace{.1ex}+,+}_{n,0}I{}^{\hspace{.1ex}+,+}_{n,0}(\eta _0,\theta ,\varphi )\).

In the rest of this subsection we give some results which shed some further light on the special nature of the principal mode. First, we show that the compatibility condition and the normalization depend only on the coefficients of the principal mode. We will need the following series expansion.

Proposition 5.2

([7]) For all \(\alpha \in {\mathbb {C}}\),

$$\begin{aligned} (\cosh \eta -\cos \theta )^{-\alpha } = \frac{1}{\Gamma (\alpha )}\sqrt{\frac{2}{\pi }} \frac{e^{-i\pi (\alpha -1/2)}}{(\sinh \eta )^{\alpha -1/2}} \sum _{n-0}^\infty \epsilon _n \cos (n\theta ) \, Q_{n-1/2}^{\alpha -1/2}(\cosh \eta ) . \end{aligned}$$

Lemma 5.3

 

  1. (a)

    The compatibility condition (9) applied to h of the form (17) is equivalent to

    $$\begin{aligned} \sum _{n=0}^\infty \varepsilon _n^2\, Q_{n-1/2}^1(\cosh \eta _0)\, b{}^{\hspace{.1ex}+,+}_{n,0} = 0. \end{aligned}$$
    (29)
  2. (b)

    The normalization condition (10) applied to f of the form (11) is equivalent to

    $$\begin{aligned} \sum _{n=0}^\infty \varepsilon _n^2\, Q_{n-1/2}^1(\cosh \eta _0)\, a{}^{\hspace{.1ex}+,+}_{n,0} = -\frac{c}{4\pi \sqrt{2} }. \end{aligned}$$
    (30)

Proof

(a) Since \(\int _{0}^{2\pi }\Phi _n^-(\theta )\,d\theta =0\), while \(\int _{0}^{2\pi }\Phi _0^+ (\varphi )\, d\varphi =2\pi \), \(\int _{0}^{2\pi }\Phi _m^+ (\varphi )\, d \varphi =0\) for \(m\ge 1\),

$$\begin{aligned} \iint _{\partial \Omega _{\eta _{0}}}h (\theta ,\varphi ) \, dS&= \int \limits _{0}^{2\pi } \int \limits _{0}^{2\pi }h(\theta ,\varphi ) \, \frac{s_{0}}{(t_{0}-\cos \theta )^{2}} \, d\theta d\varphi \\&= s_{0}\sum _{m}\sum _{n}b_{n,m}^{\nu ,\mu }\int \limits _{0}^{2\pi }(t_{0}-\cos \theta )^{-3/2}\Phi _{n}^{\nu }(\theta ) \, d\theta \int \limits _{0}^{2\pi }\Phi _{m}^{\mu }(\varphi ) \, d\varphi \\&=-4\sqrt{2}\,\pi \sum _{n=0}^{\infty }\varepsilon _{n}^{2}b_{n,0}^{+,+}Q_{n-\frac{1}{2}}^{1}(\cosh \eta _{0}), \end{aligned}$$

with the last equality following from Proposition 5.2 with \(\alpha =3/2\). The proof of (b) follows the same lines as (a). \(\square \)

The following application is of independent interest.

Proposition 5.4

The area of \({\partial \Omega _{\eta _0}}\) is equal to

$$\begin{aligned} \alpha (\eta _0) = -8 \sum _{n=0}^\infty \varepsilon _n^3 q_{n,0}q_{n,1}. \end{aligned}$$

Proof

Take \(a{}^{\hspace{.1ex}+,+}_{0,0}=(\sqrt{2}/\pi )q_{0,0}\), which gives \(f{}^{\hspace{.1ex}+,+}_{0}=1\) identically. Then apply (30) to evaluate \(\alpha (\eta _0)=\int _{\partial \Omega _{\eta _0}}f{}^{\hspace{.1ex}+,+}_{0}\,dS\). \(\square \)

Thus we obtain from (30) a formula for normalizing any given solution which has been found for the toroidal Neumann problem:

Corollary 5.5

Let f be a particular solution of \(\Lambda f=h\) and set \(c_1=\int _{\partial \Omega _{\eta _0}}f\,dS\). Let \({{\hat{f}}}\) be obtained by replacing the coefficients \(a{}^{\hspace{.1ex}+,+}_{n,0}\) for f with

$$\begin{aligned} {{\hat{a}}}{}^{\hspace{.1ex}+,+}_{n,0} = a{}^{\hspace{.1ex}+,+}_{n,0} + \varepsilon _n \frac{\sqrt{2}}{\pi } \frac{q_{n,1}}{\alpha (\eta _0)}(c-c_1). \end{aligned}$$

Then \({{\hat{f}}}\) is the unique solution of the Neumann problem which satisfies the normalization condition (10).

5.2 Determination of parameter for convergence for non-principal modes

We assume now that \((m,\nu ,\mu )\not =(0,+1,+1)\). According to Proposition 3.1, there is exactly one algebraic solution of (18) for which (11) converges. The determination of the corresponding unique value \(a_\text {opt}\) of \(a{}^{\hspace{.1ex}+,\mu }_{0,m}\) or \(a{}^{\hspace{.1ex}-,\mu }_{1,m}\) is not an algebraic question; in fact, we believe that \(a_\text {opt}\) is not the solution of any algebraic or transcendental equation associated with the Neumann data. We will assume that \(\mu =+1\) since the case \(\mu =-1\) is analogous, the only difference being the initial indexing from \(n=1\) instead of \(n=0\).

To simplify the notation, the mode being fixed we will write \(a_n=a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\) and let \(A_n(a)\) denote the value of \(a_n\) in the solution of equations (18) determined by setting the arbitrary parameter \(a_0 = a{}^{\hspace{.1ex}+,\mu }_{0,m}\) (or \(a{}^{\hspace{.1ex}-,\mu }_{1,m}\)) equal to a. Thus \(A_0(a)=a\), and by a simple induction, one finds recursively defined linear expressions

$$\begin{aligned} A_n(a)&= C_n a + D_n \ (n\ge 0), \end{aligned}$$
(31)

where

$$\begin{aligned} \begin{aligned} C_0&= 1,&D_0&=0 , \\ C_1&= \frac{-\sigma _0}{\tau _0},&D_1&= \frac{b_0}{\tau _0} ,\\ C_{n+1}&= \frac{-1}{\tau _n}(\rho _ nC_{n-1} + \sigma _nC_n), \quad&D_{n+1}&= \frac{-1}{\tau _n}(\rho _nD_{n-1} + \sigma _nD_n - b_n) \ \ (n\ge 1). \end{aligned} \end{aligned}$$
(32)

By construction, every value of \(a\in {\mathbb {R}}\) gives an algebraic solution \(\{A_n(a)\}\) of the system (18). According to (11) and Theorem 3.5, we want to find the unique value \(a_\text {opt}\) provided by Theorem 5.6 for which the series

$$\begin{aligned} \sum _{n=0}^\infty A_n(a_\text {opt}) \, \Phi {}^{\hspace{.1ex}\nu }_{n}(\theta )\Phi {}^{\hspace{.1ex}\mu }_{m}(\varphi ) \end{aligned}$$
(33)

converges absolutely and thus gives a boundary function \(f_{m}^{\nu ,\mu }(\theta ,\varphi )\). In particular, it is necessary that \( A_n(a_\text {opt})\rightarrow 0\) as \(n\rightarrow \infty \). By (31), this says \(C_n a_\text {opt}+ D_n \rightarrow 0\).

Note that the \(C_n\) do not depend on the data \(\{b_n\}\). It is clear that two consecutive terms \(C_n\), \(C_{n+1}\) can never vanish. Under the assumption that the \(C_n\) do not tend to zero, i.e., \(C_n>\epsilon >0\) for infinitely many n, we have

$$\begin{aligned} -\frac{D_n}{C_n} \rightarrow a_\text {opt}\end{aligned}$$
(34)

as \(n\rightarrow \infty \) on that subsequence. We will look further into the question of the asymptotics of \(C_n\) in the following section.

We summarize our conclusions as follows.

Theorem 5.6

Let the coefficients \(\{b{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\}\) be such that the series (17) converges absolutely, defining \(h:{\partial \Omega _{\eta _0}}\rightarrow {\mathbb {R}}\). Assume that \(b{}^{\hspace{.1ex}+,+}_{n,0}\) satisfy (29), so h satisfies the compatibility condition (9). Suppose further that the continuous solution \(f:{\partial \Omega _{\eta _0}}\rightarrow {\mathbb {R}}\) of \(\Lambda f=h\) specified in Proposition 3.1 has a double Fourier series which converges absolutely. Then (i) for every value of \(a{}^{\hspace{.1ex}+,+}_{0,0}\in {\mathbb {R}}\), the resulting algebraic solution for the sequence \(\{a{}^{\hspace{.1ex}+,+}_{n,0}\}\) produces an absolutely convergent series \(\sum _{n} a{}^{\hspace{.1ex}+,+}_{n,0}I{}^{\hspace{.1ex}+,+}_{n,0}(\eta _0,\theta ,\phi )\) whose value is \(f{}^{\hspace{.1ex}+,+}_{0}\) plus a constant. Further, (ii) for \((m,\nu ,\mu )\) different from \((0,+1,+1)\), there exists a unique value of \(a{}^{\hspace{.1ex}\nu ,\mu }_{0,m}\) (when \(\nu =1\)) or \(a{}^{\hspace{.1ex}\nu ,\mu }_{1,m}\) (when \(\nu =-1\)) for which the resulting algebraic solution gives a convergent series \(\sum _{n} a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}(\eta _0,\theta ,\phi )\). The sum of this series is the mode \(f{}^{\hspace{.1ex}\nu ,\mu }_{m}\) given by Eq. (27).

The contrast of the behavior of the principal mode in this theorem is striking, since Eq. (18) show no apparent structural difference between the principal and other modes. The algorithm is essentially as follows:

Algorithm 5.7

(Solution to the Neumann problem) Given the coefficients \(\{b{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\}\), calculate the coefficients \(C_n,D_n\), \(0\le n\le N\), by (32) which give the approximation \(a_\text {opt}=a_\text {opt}(m,\nu ,\mu )=-D_N/C_N\) and then also \(\{a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\}\) by (31) and (34).

6 Numerical results

6.1 Calculation of \(q_{n,m}\)

From (19) it is clear that \(q_{n,m}\) decreases exponentially as \(n\rightarrow \infty \). In fact, for many values of \(\eta _0\) in the range we will be considering, \(|q_{100,0}|\) is less than \(10^{-300}\) and thus effectively becomes zero when represented in standard IEEE double-precision format. For many problems values of n as large as 100 would not be necessary, but if one wants to use a large number of terms in a series solution, two ways out of this difficulty present themselves.

The first way is via asymptotic expansions. For large values of n we can estimate \(q_{n,m}\) to high accuracy by means of (19) and Stirling’s formula for the \(\Gamma \) function. However, it turns out there is a gap between small values of n where machine precision is applicable and the large values of n for which Stirling’s estimate is sufficiently accurate. For example, for \(\eta _0=1.5\) and \(m=2\), we find in machine precision that \(q_{131,2}\) vanishes, while \(q_{130,2}\) suffers a relative error of 0.0006. The estimate based on (19) gives a relative error of 0.003. The formula (19) is in fact the first term of an asymptotic expansion for \(Q_n^m\) given in [17], and can be made much more accurate by multiplying it by the factor

$$\begin{aligned} 1 - \frac{1}{8n} -\frac{m^2}{n} + \frac{4m^2-1}{n}\frac{1}{1-e^{-2\eta _0}}. \end{aligned}$$

Note that by (15), we only need to know the ratios \(q_{n+1,m}/q_{n,m}\). After cancelling several factors, one derives an approximation which can be substituted into (15) to find, for example, that \(\tau _{136,2}\) vanishes in machine precision, while the approximations of \(\tau _{n,2}\) already for \(n\ge 20\) have relative error less than 0.00001. Similar statements hold for the coefficients \(\rho _{n,m}\) and \(\sigma _{n,m}\). The optimal cutoff value of n at which to switch from a direct calculation to the asymptotic formula would depend on the values of m and \(\eta _0\) as well as a consideration of the acceptable error in the coefficients. It is likely that other more convenient asymptotic expressions could be developed but this would lead us rather far afield.

Multiple-precision arithmetic was used in Mathematica to verify the correctness of the above statements. Multiple precision is also a convenient second alternative for calculating all of the Neumann constants when one is planning to solve many Neumann problems on the same torus: once the coefficients have been obtained and rounded to machine precision, the computation will then easily move forward without further recourse to higher precision arithmetic.

6.2 Numerical behavior of \(C_n\)

We observed that \(a_\text {opt}\) is given by (34) unless it happens that \(C_n\rightarrow 0\). (Recall that \(C_n\) does not depend on the Neumann data \(\{b{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\}\).) For small values of n, we have a priori little control over even the sign of the coefficients defined in (15). However, from (22), \(\rho _{n,m}/\tau _{n,m}\rightarrow 1\) and \(\sigma _{n,m}/\tau _{n,m}\rightarrow -2t_0\). Therefore if for a single large n we have

$$\begin{aligned} C_n \approx e^{\eta _0}C_{n-1}, \end{aligned}$$

then by (32) it would follow that

$$\begin{aligned} C_{n+1} \approx -C_{n-1}+2t_0C_n = -e^{-\eta _0}C_n+2t_0C_n = e^{\eta _0}C_n; \end{aligned}$$

i.e., the sequence \(\{C_n\}\) grows exponentially. Table 1 lists calculated values of \(C_n\) corresponding to \(m=1\) and a range of values of \(\eta _0\). Other values of m are shown in Table 2. Even though the initial values can decrease, in all cases that we have examined it appears that \(C_n\rightarrow \infty \) exponentially as \(n\rightarrow \infty \), and we conjecture that this always holds. Then Algorithm 5.7 is applicable.

Table 1 Sample values of \(C_n\) for \(m=1\), \((\nu ,\mu )=(+,+)\). 200-digit precision was used to avoid underflow in the calculations
Full size table
Table 2 Sample values of \(C_n\) for \(\eta =0.4\), \((\nu ,\mu )=(+,+)\)
Full size table

6.3 Examples of Dirichlet-to-Neumann and Neumann problem calculations

We illustrate the solution of the Neumann problem with numerical examples. For the principal mode \((m,\nu ,\mu )=(0,+1,+1)\) there is not much to be said. Every choice of \(a_0=a{}^{\hspace{.1ex}++}_{0,0}\) gives an algebraic solution which defines a series which solves the Neumann problem. Corollary 5.5 may be used for normalization. It may be worth noting that when \(b_n=0\) for all n (i.e., the coefficients for a vanishing normal derivative in the mode under consideration), since the sum is a constant, one can dispense with the series altogether for this case. But if one does wish to follow the procedure, the formulas (32) give \(D_n=0\) always, and by (31), we have \(A_n(a)=aC_n\). Choosing \(a_0=1\) without any loss of generality and fixing \(\eta _0\), one obtains the values \(C_n\) by (32) and then the initial coefficients \(a_n\) by (31). This calculation amounts simply to finding values of the associated Legendre functions of the second kind via the classical recursion formulas, and the only numerical error is that which accumulates due to roundoff.

The numerical examples were calculated in Mathematica on a household laptop computer. All calculations were in machine precision except a few of the coefficients as noted. The calculations took a fraction of a second apart from the graphics.

Example 1

To illustrate the calculation of the Dirichlet to Neumann mapping, consider \(a{}^{\hspace{.1ex}+,+}_{n,m}=a_n\) where \(a_n = (-1)^{(n-1)/2} n^{-2}\) for n odd, and \(a_n=0\) for n even, \(0\le n\le 50\). Figure 1 shows the function f and the image \(h=\Lambda f\) for a few values of m. While f is smooth, it does not satisfy the regularity condition of Theorem 3.5, and h has jump discontinuities. Except at the jumps, the values of h agree to within \(10^{-14}\) to \(10^{-15}\) with the normal derivative of the f series obtained by numerical derivation.

As another example, let instead \(a_n = (-1)^{(n-1)/2} n^{-1}\) for n odd. As Fig. 2 shows, now f has a jump discontinuity while the series for the resulting h does not converge, further illustrating the tendency of the Dirichlet to Neumann mapping to reduce the degree of differentiability.

Fig. 1
figure 1

Given function \(f(\theta ,\varphi )\) determined by \(a_n = (-1)^{(n-1)/2} n^{-2}\) for n odd (above) and calculated D-to-N mapping \(h(\theta ,\varphi )\) (below). \(\eta _0=1.5\). All calculations carried out with machine precision

Full size image
Fig. 2
figure 2

Solution determined by \(a_n = (-1)^{(n-1)/2} n^{-1}\) for n odd, \(m=2\), \(\eta _0=1.5\). The profile of h for \(\varphi =0\) (right) results from an attempt to graph a nonconvergent series

Full size image

Example 2

Let

$$\begin{aligned} u = \left( \frac{\sinh \eta }{\cosh \eta -\cos \theta } \right) ^m \cos m \varphi . \end{aligned}$$
(35)

It is readily checked that u is harmonic and

$$\begin{aligned} {{\,\mathrm{{\text{ nor }}}\,}}u = m \left( \frac{\sinh \eta _0}{ \cosh \eta _0-\cos \theta }\right) ^m \left( (\cosh \eta _0-\cos \theta )\coth \eta _0 +\sinh \eta _0\right) \cos m \varphi . \end{aligned}$$
(36)

(Clearly one also would obtain a harmonic function with \(\sin m\varphi \) in place of \(\cos m\varphi \) in (35).) By Proposition 5.2, the coefficients in the series for u are equal to

$$\begin{aligned} a{}^{\hspace{.1ex}\nu ,\mu }_{n,m} = (-1)^m \frac{\sqrt{2/\pi }}{\Gamma (m + 1/2)} \, \varepsilon _n q_{n,m}. \end{aligned}$$
(37)

We substitute these coefficients into (18) to obtain numerical values for the \(b{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\). Then we compare truncations of the series (21) with the true values of \(h={{\,\mathrm{{\text{ nor }}}\,}}u\) according to (36). Figure 3 displays the base-10 logarithm of the absolute error for different combinations of m and \(\eta _0\). As is expected, the error is reduced when the number of terms in the series increases. It is also seen that the error increases steadily when larger values of m and \(\eta _0\) are used.

Fig. 3
figure 3

Base-10 logarithm indicating number of significant figures of approximation of the Dirichlet-to-Neumann mapping given by Eq. (18) truncating the series (21) to \(0\le n\le N\) for varying values of N. Accuracy is lost as m or \(\eta _0\) increases. For comparison purposes 100-digit precision was used to obtain the Neumann constants. It was found that machine precision was sufficient for \(\eta _0>0.5\), which is applicable for most “reasonably-shaped” toroidal domains

Full size image

Example 3

We illustrate our algorithm for solving the Neumann problem using the same function u as in the previous example. The Fourier coefficients \(b{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\) are obtained by numerical integration. Then the auxiliary coefficients \(C_n\), \(D_n\) are obtained recursively by (32), and then \(a_\text {opt}\) is approximated by the last value of \(-C_n/D_n\) according to (34). One would expect that the values \(a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}=A_n(a_\text {opt})\) of (34) provide a convergent series, while for \(a\not =a_\text {opt}\), \(\{A_n(a)\}\) would not. This is confirmed by Fig. 4, which shows the values of \(A_n(a+\epsilon )\) for small values of \(\epsilon \). The error in a particular series solution h of the Neumann problem compared to (36) is shown in Fig. 5. Maximum errors for combinations of \(\eta _0\), m are shown in Table 3.

Fig. 4
figure 4

Rapid growth of the first 50 coefficients in nonconvergent algebraic solutions generated by to \(a_\text {opt}+\epsilon \), illustrated for \(\eta _0=0.4\) and \(m=2\), with \(a_\text {opt}\) approximated by \(-D_{50}/C_{50}\). (The graphic is truncated: for \(\epsilon =.1\), the coefficients reach approximately \(10^7\). Even at this scale, the coefficients for \(\epsilon =0\) are virtually indistinguishable from the horizontal axis

Full size image
Fig. 5
figure 5

Error in solution for Neumann problem for \(\eta _0=0.4\), \(m=2\) and 50 terms, distributed over the range \(0\le \theta \le 2\pi \), with \(\varphi =0\)

Full size image

6.4 Computational complexity

In order to calculate the Dirichlet-to-Neumann mapping or solve the Neumann problem in the context of the Fourier representation, one needs the coefficients \(\{a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\}\) or \(\{b{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\}\) respectively. Methods for calculating these coefficients given the values of f or h on \({\partial \Omega _{\eta _0}}\) are well known, depend on a choice of numerical integration procedure, can be found in many standard software packages, and will not be discussed here. In this sense, Algorithm 5.7 can be considered as one part of a larger computation.

The amount of calculation required for Algorithm 5.7 is simple to estimate. Take N terms in the truncated Fourier series and use the values \(C_N\!,D_N\) to determine \(a_\text {opt}\). Assume that the Neumann constants \(\rho _{n,m},\sigma _{n,m},\tau _{n,m}\) have been calculated previously. This is also reasonable when one wishes to solve many Neumann problems on a single surface. Then by (32), approximately 4N multiplications in machine-precision are needed to find \(a_\text {opt}\), and since we then have all of the \(C_n,D_n\) for \(n<N\), by (31) we need N more multiplications to find the Taylor coefficients for each mode. The calculation of the Neumann constants is also of order O(NM), but one must take into account that either multiple precision or a slightly more complicated asymptotic formula may be necessary as described in Sect. 4 for a certain collection of coefficients of low index.

Table 3 Significant figures in the numerical solution of the Neumann problem on the torus showing the increase in accuracy with the number of terms
Full size table

7 Exterior toroidal domain and toroidal shells

7.1 Exterior domain

The formula for the normal derivative of an exterior harmonic function \(u\in {{\,\mathrm{{\text{ Har }}}\,}}(\Omega _{\eta _0}^*) = \{x:\ 0\le \eta <\eta _0\},\) and the solution of the corresponding Neumann problem are quite analogous to that of the interior domain \(\Omega _{\eta _0}\). The exterior harmonics \(E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\in {{\,\mathrm{{\text{ Har }}}\,}}(\Omega _{\eta _0}^*)\), defined by

$$\begin{aligned} E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}(x)=E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}[\eta _{0}](x) = \sqrt{\cosh \eta -\cos \theta }\,\frac{P{}^{\hspace{.1ex}m}_{n-1/2}(\cosh \eta )}{P{}^{\hspace{.1ex}m}_{n-1/2}(\cosh \eta _{0})} \, \Phi {}^{\hspace{.1ex}\nu }_{n}(\theta ) \, \Phi {}^{\hspace{.1ex}\mu }_{m}(\varphi ), \end{aligned}$$
(38)

are obtained from the interior harmonics (5) by writing \(P_{n-1/2}^{m}(\cosh \eta )\) in place of \(Q_{n-1/2}^{m}(\cosh \eta )\) and are orthogonal with the same weight function (8) formula but applied in \(\Omega _{\eta _0}^*\).

These Legendre functions of the first and second kinds satisfy identical recurrence relationships [17]. For this reason, one finds that \({{\,\mathrm{{\text{ nor }}}\,}}E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\) is obtained from the formula of Lemma 3.2 by replacing similarly \(Q_{sn-1/2}^{m}(\cosh \eta )\) with \(P_{n-1/2}^{m}(\cosh \eta )\). Since the solution of the Dirichlet problem in \(\Omega _{\eta _0}^*\) with boundary condition f given by (11) is

$$\begin{aligned} u=\sum _{n,m,\nu ,\mu } a{}^{\hspace{.1ex}\nu ,\mu }_{n,m}E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}(\eta ,\theta ,\varphi ), \end{aligned}$$
(39)

one finds that the normal derivative of f will be given by equations (18) when \(q_{n,m}\) is replaced in (15) with

$$\begin{aligned} p_{n,m} = P_{n-1/2}^m(\cosh \eta _0). \end{aligned}$$
(40)

The method we have described is then applicable with no essential changes for solving the Dirichlet-to-Neumann problem in \(\Omega _{\eta _0}^*\). It is worth noting that parallel to (20) we have [17, p. 305] that

$$\begin{aligned} \lim _{n\rightarrow \infty }\frac{p_{n-1,m}}{p_{n,m}}=e^{\eta _{0}}. \end{aligned}$$
(41)

7.2 Toroidal shell

We now combine our results in the interior and exterior domains for solving the Neumann problem in a toroidal shell. Let \(\eta _{\textrm{int}}<\eta _{{\textrm{ext}}}\). Common to an interior and an exterior domain are the points of the toroidal shell

$$\begin{aligned} \Omega = \Omega _{\eta _{\textrm{int}},\eta _{\textrm{ext}}} = \Omega _{\eta _{\textrm{ext}}} \cap \Omega _{\eta _{\textrm{int}}}^* . \end{aligned}$$

A general harmonic function u in \(\Omega _{\eta _{\textrm{ext}},\eta _{\textrm{int}}}\) and continuous in the closure can be expressed via an integral of its boundary values over \(\partial \Omega _{\eta _{\textrm{ext}},\eta _{\textrm{int}}}\) using the Poisson kernel for the torus [16, Ch. 1]. This integral is the difference of the integrals over \(\partial \Omega _{\eta _{\textrm{ext}}}\) and \(\partial \Omega _{\eta _{\textrm{int}}}\), which give a decomposition \(u=u_0+u_1\) with \(u_0\in {{\,\mathrm{{\text{ Har }}}\,}}\Omega _{\eta _{\textrm{int}}}\) and \(u_1\in {{\,\mathrm{{\text{ Har }}}\,}}\Omega _{\eta _{\textrm{ext}}}^*\). Consequently, we may express u as the sum of two series

$$\begin{aligned} u = \sum _{n,m,\nu ,\mu } c{}^{\hspace{.1ex}{\textrm{int}}\,\nu ,\mu }_{n,m}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m} + \sum _{n,m,\nu ,\mu } c{}^{\hspace{.1ex}{\textrm{ext}}\,\nu ,\mu }_{n,m}E{}^{\hspace{.1ex}\nu ,\mu }_{n,m} , \end{aligned}$$
(42)

analogous to the Laurent series for holomorphic functions in an annular domain in the complex plane, converging uniformly in proper closed subdomains. (Note, however, that the inner and outer harmonics together do not form an orthogonal system in \(\Omega _{\eta _{\textrm{ext}},\eta _{\textrm{int}}}\).)

A boundary function \(f:\partial \Omega \rightarrow {\mathbb {R}}\) is given collectively by its values for \(\eta =\eta _{\textrm{int}}\) and \(\eta =\eta _{\textrm{ext}}\) collectively, let us say

$$\begin{aligned} f_{\textrm{int}}(\theta ,\varphi )= f(\eta _{\textrm{int}},\theta ,\varphi )&= \sum a{}^{\hspace{.1ex}{\textrm{int}}\,\nu ,\mu }_{n,m} I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}[\eta _{1}] , \nonumber \\ f_{\textrm{ext}}(\theta ,\varphi )= f(\eta _{\textrm{ext}},\theta ,\varphi )&= \sum a{}^{\hspace{.1ex}{\textrm{ext}}\,\nu ,\mu }_{n,m} E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}[\eta _{0}]. \end{aligned}$$
(43)

For u to be the solution of the Dirichlet problem for f, we combine (42) with (43) to find

$$\begin{aligned} q{}^{\hspace{.1ex}\textrm{int}}_{n,m} c{}^{\hspace{.1ex}{\textrm{int}}\,\nu ,\mu }_{n,m} + p{}^{\hspace{.1ex}\textrm{int}}_{n,m} c{}^{\hspace{.1ex}{\textrm{ext}}\,\nu ,\mu }_{n,m}&= a{}^{\hspace{.1ex}{\textrm{int}} \,\nu ,\mu }_{n,m}, \nonumber \\ q{}^{\hspace{.1ex}\textrm{ext}}_{n,m} c{}^{\hspace{.1ex}{\textrm{int}}\,\nu ,\mu }_{n,m} + p{}^{\hspace{.1ex}\textrm{ext}}_{n,m} c{}^{\hspace{.1ex}{\textrm{ext}}\,\nu ,\mu }_{n,m}&= a{}^{\hspace{.1ex}{\textrm{ext}} \,\nu ,\mu }_{n,m}, \end{aligned}$$
(44)

where

$$\begin{aligned} \begin{aligned} q{}^{\hspace{.1ex}\textrm{int}}_{n,m}&= Q_{n-1/2}^m(\cosh \eta _{\textrm{int}}) , \quad q{}^{\hspace{.1ex}\textrm{ext}}_{n,m}&\,=\,&Q_{n-1/2}^m(\cosh \eta _{\textrm{ext}}) , \\ p{}^{\hspace{.1ex}\textrm{int}}_{n,m}&= P_{n-1/2}^m(\cosh \eta _{\textrm{int}}) , \quad p{}^{\hspace{.1ex}\textrm{ext}}_{n,m}&\,=\,&P_{n-1/2}^m(\cosh \eta _{\textrm{ext}}) . \end{aligned} \end{aligned}$$

This might be written symbolically as

$$\begin{aligned} \begin{pmatrix} q^{\textrm{int}} &{} p^{\textrm{int}} \\ q^{\textrm{ext}} &{} p^{\textrm{ext}} \end{pmatrix} \begin{pmatrix} c^{\textrm{int}} \\ c^{\textrm{ext}} \end{pmatrix} = \begin{pmatrix} a^{\textrm{int}} \\ a^{\textrm{ext}} \end{pmatrix} . \end{aligned}$$

To solve this system, one needs to verify that it is nonsingular. Instead of a direct verification as in Lemma 4.1, we simply note that if for even one combination of \((n,m,\nu ,\mu )\) there were more than one solution, one could easily construct a Dirichlet problem in the shell \(\Omega \) with more than one solution.

We see that \({{\,\mathrm{{\text{ nor }}}\,}}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\vert _{\partial \Omega _{\textrm{int}}}\) is obtained from the formula of Lemma 3.2 with \(\eta _0\) replaced with \(\eta _{\textrm{int}}\), while \({{\,\mathrm{{\text{ nor }}}\,}}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\vert _{\partial \Omega _{\textrm{ext}}}\) is obtained by using \(\eta _{\textrm{ext}}\) instead. The boundary values \({{\,\mathrm{{\text{ nor }}}\,}}E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\vert _{\partial \Omega _{\textrm{int}}}\) and \({{\,\mathrm{{\text{ nor }}}\,}}E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\vert _{\partial \Omega _{\textrm{ext}}}\) are then obtained by replacing \(Q_{n-1/2}^m\) with \(P_{n-1/2}^m\). Once we have the harmonic function u as in (42), we have then

$$\begin{aligned} {{\,\mathrm{{\text{ nor }}}\,}}u\bigg \vert _{\partial \Omega _{\textrm{int}}} = \sum _{n,m,\nu ,\mu } c{}^{\hspace{.1ex}{\textrm{int}}\,\nu ,\mu }_{n,m} {{\,\mathrm{{\text{ nor }}}\,}}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m} \bigg \vert _{\partial \Omega _{\textrm{int}}}+ \sum _{n,m,\nu ,\mu } c{}^{\hspace{.1ex}\textrm{ext}\,\nu ,\mu }_{n,m} {{\,\mathrm{{\text{ nor }}}\,}}E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\bigg \vert _{\partial \Omega _{\textrm{int}}} , \\ {{\,\mathrm{{\text{ nor }}}\,}}u\bigg \vert _{\partial \Omega _{\textrm{ext}}} = \sum _{n,m,\nu ,\mu } c{}^{\hspace{.1ex}\textrm{int}\,\nu ,\mu }_{n,m} {{\,\mathrm{{\text{ nor }}}\,}}I{}^{\hspace{.1ex}\nu ,\mu }_{n,m} \bigg \vert _{\partial \Omega _{\textrm{ext}}}+ \sum _{n,m,\nu ,\mu } c{}^{\hspace{.1ex}{\textrm{ext}}\,\nu ,\mu }_{n,m} {{\,\mathrm{{\text{ nor }}}\,}}E{}^{\hspace{.1ex}\nu ,\mu }_{n,m}\bigg \vert _{\partial \Omega _{\textrm{ext}}} . \end{aligned}$$

When the convergence of the series is absolute, one may apply the same rearranging and reindexing as described in the proof of Theorem 3.5 to obtain the coefficients in the Dirichlet-to-Neumann mapping \(h=\Lambda f\),

$$\begin{aligned} h(\eta _{\textrm{int}},\theta ,\varphi )&= \sqrt{\cosh \eta _0-\cos \theta }\sum b{}^{\hspace{.1ex}{\textrm{int}}\,\nu ,\mu }_{n,m} \Phi {}^{\hspace{.1ex}\nu }_{n}(\theta ) \Phi {}^{\hspace{.1ex}\mu }_{m}(\varphi ) ,\nonumber \\ h(\eta _{\textrm{ext}},\theta ,\varphi )&= \sqrt{\cosh \eta _0-\cos \theta }\sum b{}^{\hspace{.1ex}{\textrm{ext}}\,\nu ,\mu }_{n,m} \Phi {}^{\hspace{.1ex}\nu }_{n}(\theta ) \Phi {}^{\hspace{.1ex}\mu }_{m}(\varphi ). \end{aligned}$$
(45)

As in the solution of the Neumann problem for the interior domain, the equations for a fixed value of \((m,\nu ,\mu )\) are independent of those for another value of these indices. They can be solved recursively. The only difference will be that one must solve a pair of equations at each step.

8 Conclusions

We have presented an approach to studying the Dirichlet-to-Neumann mapping and solving the Neumann problem for the Laplacian on a torus. It is shown how the Dirichlet-to-Neumann mapping may be expressed by means of certain infinite series based on toroidal harmonics. We describe the well-known necessary and sufficient condition for the solvability of the Neumann problem (compatibility condition), as well as the normalization condition in terms of the Fourier coefficients. These results show that in this context the Neumann problem involves an infinite upper-triangular system of linear equations. The solution of the problem involves a special twist in that the unique value of the free parameter in this underdetermined linear system which truly gives a solution, cannot be found algebraically. Therefore we express it as a limit of easily calculated algebraic expressions. Numerical results are displayed for the accuracy of the algorithm. While the theoretical justification requires a fairly strict smoothness assumption, the method appears applicable in great generality. The Neumann problem is also solved in a toroidal shell.