On Roots of the Macdonald Function

by

Kaylee Robert Tejeda

B.S., Pure Mathematics, University of New Mexico, 2003

THESIS

Submitted in Partial Fulfillment of the Requirements for the Degree of

Master of Science Mathematics

The University of New Mexico

Albuquerque, New Mexico

July, 2014 \begin{matrix}c \\ ◯ \end{matrix}2014, Kaylee Robert Tejeda

iii

Dedication

This is for all seekers of the truth, regardless of consequence.

“In the twentieth century it has been understood that all knowledge is dependent upon the

question asked; and the relationship of mathematics to nature is one of the profound indicators, I think, that truth can be known. Maybe not ‘the’ truth, but I always think of the positivist philosopher Wittgenstein who was once asked, in a classroom situation about a certain proposition, ‘Is it the truth?’, and he said ‘Well, it’s certainly true enough!’; and that is where we are with our modelling of the world and with our mathematics. It is the

truest truth we know. It is true enough.” – Terence McKenna

Acknowledgments

I would like to thank my advisor Stephen Lau for his generous guidance, in addition to all remaining UNM Department of Mathematics and Statistics faculty and staff, particularly Ana Parra Lombard, for their assistance and mentorship over the years. I thank my family, for without them I would not exist. I would also like to thank those who believed in me, for their moral support; as well as those who doubted me, for the inspiration to prove otherwise. Some of this work was supported by NSF grant No. PHY 0855678.

On Roots of the Macdonald Function

by

Kaylee Robert Tejeda

B.S., Pure Mathematics, University of New Mexico, 2003

M.S., Mathematics, University of New Mexico, 2014

Abstract

An overview is given for the Dirichlet-to-Neumann map for outgoing solutions to the “radial wave equation” in the context of nonreflecting radiation boundary conditions on a spherical domain. We then consider the Macdonald function K_{ℓ+1/2}(z) for ℓ ∈ ℤ_{≥0}, a solution to the half-integer order modified Bessel equation. This function can be expressed as K_{ℓ+1/2}(z) = \sqrt{\frac{π}{2z}}e^{−z}z^{−ℓ}p_{ℓ}(z), where p_{ℓ}(z) is a degree-ℓ monic polynomial with simple roots in the left-half plane. By exploiting radiation boundary conditions for the “radial wave equation”, we show that the root set of p_{ℓ}(z) also obeys ℓ additional polynomial constraints. These constraints are in fact Newton’s identities which relate a polynomial’s coefficients to the power sums of its roots. We follow this with numerical verification up to order ℓ = 20.

Contents

List of Figures ix List of Tables x Glossary xi 1 Introduction 1 2 Radiation boundary conditions and the radial wave equation 2

2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 The radial wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.3 Macdonald function as an outgoing solution . . . . . . . . . . . . . . . . . . 5

3 Nonreflecting Boundary Conditions 7 3.1 Proof of Theorem 3.0.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2 Alternative proof. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4 Newton’s Identities 14 Contents

4.1 Main result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 4.2 Proof of Theorem 4.1.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

5 Numerical root evaluation and identity check for several values of ℓ 19 5.1 Root finding and identity verification for ℓ = 20 to 27 digits of precision … 20

6 Conclusion 30

A Location and count of Macdonald roots 31 A.1 Root count for K_{ℓ+1/2}(z) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 A.2 Root count for K_{ν}(z) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

References 36

List of Figures

2.1 Scaled roots of K_{ℓ+1/2}(z). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.1 Radial axis with support of Ψ_{ℓ}(0, r) and point at r_{B} − ϵδ (the dashed line). 8 3.2 Retarded time axis and support of profile function f. . . . . . . . . . . . . 8 A.1 Keyhole contour with arcs Γ : |z| = R and γ : |z| = δ . . . . . . . . . . . . . 34

List of Tables

5.1 Roots of K_{ℓ+1/2}(z) for ℓ = 1, ..., 20 . . . . . . . . . . . . . . . . . . . . . . . 22

Glossary

ψ_{ℓ} General solution to the radial wave equation Ψ_{ℓ} General solution defined to remove \frac{1}{r} decay Ψ̂_{ℓ} Laplace transform of the general non-decaying solution above Ω_{ℓ} Radiation boundary kernel for the radial wave equation Y_{ℓm}(θ, ϕ) Spherical harmonics I_{ν}(z), K_{ν}(z) Modified Bessel functions p_{ℓ}(z) Bessel polynomial b_{ℓn} Roots of Bessel polynomial W(r) Wronskian G(r, ξ; s) Green’s function

Chapter 1

Introduction

The overall goal of this thesis is to explore a set of identities that arise in the context of outgoing waves and the half-integer order modified Bessel equation. We derive identities for the roots of the Macdonald function, a solution to the modified Bessel equation. These identities prove to be Newton’s identities which relate a polynomial’s coefficients to the power sums of its roots. Our derivation of these identities exploits the concept of nonreflecting boundary conditions for the wave equation in the presence of a spherical “computational” domain. These identities are then verified numerically using Mathematica up to a specified order. The actual code used, along with the tables of roots generated and numerical verification of the identity for each root set are presented in the final chapter. Arguments [2] by Watson concerning the location of the roots of Macdonald’s function are presented in the Appendix.

Both subscript style notation ∂_{r} and ratio notation \frac{∂}{∂r} are used interchangeably in this paper to represent the partial derivative with respect to the radial variable r (for example), with the choice of notation based on readability and typesetting. Sometimes we also use, for example, the notation Ψ_{tt} = ∂_{t}^{2}Ψ to denote derivatives.

Radiation boundary conditions and

the radial wave equation

2.1 Overview

In mathematical physics, wave phenomena are typically modelled by hyperbolic partial differential equations on an unbounded domain. Examples include acoustic, electromagnetic, elastic, and gravitational phenomena. In numerical situations, this unbounded domain is usually truncated to produce a computational domain, with “artificial” boundary conditions (BCs) placed at the domain’s boundary, an artificial boundary. One approach for producing BCs which achieve exact domain reduction (i.e. the solution on the computational domain agrees with the restriction thereon of the solution on the unbounded domain) involves the Dirichlet-to-Neumann (DtN) map [3].

Consider the case of the ordinary wave equation on ℝ^{3},

Δψ − \frac{1}{c^{2}} \frac{∂^{2}ψ}{∂t^{2}} = 0. (2.1)

Chapter 2. Radiation boundary conditions and the radial wave equation

The initial value problem for the wave equation can be solved using spherical coordinates, chosen for symmetry considerations. We will also assume unit speed c = 1 in the rest of this paper for simplicity.

Let D = {x||x| ≤ r_{B}} contain the support of the initial data. Then BCs are constructed on the boundary ∂D with the purpose of approximating the exact solution for the wholespace problem restricted to D. We want BCs which are transparent (i.e. nonreflecting). An otherwise approximate BC which does not behave in this manner will generate a spurious reflection, leading to numerical error in the computational domain D. Such BC’s are specified in the next chapter. Here we lay some groundwork for their formulation.

2.2 The radial wave equation

In spherical coordinates the Laplace operator has the form

Δψ = \frac{1}{r^{2}} \frac{∂}{∂r}r^{2} \frac{∂}{∂r}ψ + \frac{1}{r^{2}} Δ_{S^{2}}ψ,  Δ_{S^{2}} = \frac{1}{sin θ} \frac{∂}{∂θ} sin θ\frac{∂}{∂θ} + \frac{1}{sin^{2} θ} \frac{∂^{2}}{∂ϕ^{2}} ,  (2.2)

where Δ_{S^{2}} is the Laplace operator on the unit sphere.

The general solution to (2.1) can then be expanded in terms of the eigenfunctions Y_{ℓm}(θ, ϕ) of Δ_{S^{2}}:

ψ = ∑_{ℓ=0}^{∞} ∑_{m=−ℓ}^{ℓ} ψ_{ℓm}(t, r)Y_{ℓm}(θ, ϕ) = ∑_{ℓ=0}^{∞} ∑_{m=−ℓ}^{ℓ} \frac{Ψ_{ℓm}(t, r)}{r}Y_{ℓm}(θ, ϕ). (2.3)

Here the Y_{ℓm}(θ, ϕ) are the standard spherical harmonics which obey

Δ_{S^{2}}Y_{ℓm}(θ, ϕ) = −ℓ(ℓ + 1)Y_{ℓm}(θ, ϕ). (2.4)

Plugging (2.3) into the wave equation (2.1) we find, using orthogonality, that the multipole coefficients ψ_{ℓm} and Ψ_{ℓm} satisfy

−\frac{∂^{2}ψ_{ℓ}}{∂t^{2}} + \frac{∂^{2}ψ_{ℓ}}{∂r^{2}} + \frac{2}{r} \frac{∂ψ_{ℓ}}{∂r} − \frac{ℓ(ℓ + 1)}{r^{2}}ψ_{ℓ} = 0, (2.5)

and

−\frac{∂^{2}Ψ_{ℓ}}{∂t^{2}} + \frac{∂^{2}Ψ_{ℓ}}{∂r^{2}} − \frac{ℓ(ℓ + 1)}{r^{2}}Ψ_{ℓ} = 0, (2.6)

where now we suppress the m indices. Note that (2.6) is similar to the simple 1 + 1 wave equation −∂_{t}^{2}Ψ + ∂_{r}^{2}Ψ = 0 except for the potential term V = ℓ(ℓ + 1)/r^{2}.

The 1 + 1 wave equation has solutions g(t + r), f(t − r) for general profile functions g(v), f(u). In fact, for the radial wave equation (2.6),

Ψ_{ℓ}(t, r) = ∑_{k=0}^{ℓ} \frac{1}{r^{k}} c_{ℓk}f^{(ℓ−k)}(t − r),  c_{ℓk} = \frac{1}{2^{k}k!} \frac{(ℓ + k)!}{(ℓ − k)!},  (2.7)

is the analog of the outgoing solution f(t − r) to the 1 + 1 wave equation, which we now formulate as a lemma with proof.

Lemma 2.2.1. Given a sufficiently smooth profile function f(u), the“multipole” Ψ_{ℓ}(t, r) given in (2.7) is a solution to the radial wave equation (2.6).

Proof. To prove the lemma, note that f(t − r) is clearly a solution to (2.6) for ℓ = 0, and then proceed by induction on ℓ. First, appealing to the explicit formula (2.7) for c_{ℓk}, it can be established that

Ψ_{ℓ} = −(\frac{∂}{∂r} − \frac{ℓ}{r})Ψ_{ℓ−1},  Ψ_{ℓ−1}(t, r) = ∑_{k=0}^{ℓ−1} \frac{1}{r^{k}} c_{ℓ−1,k}f^{(ℓ−1−k)}(t − r).  (2.8)

Next, assume that Ψ_{ℓ−1} in (2.8) obeys (2.7) with ℓ replaced by ℓ − 1 in the potential:

−\frac{∂^{2}Ψ_{ℓ−1}}{∂t^{2}} + \frac{∂^{2}Ψ_{ℓ−1}}{∂r^{2}} − \frac{ℓ(ℓ − 1)}{r^{2}}Ψ_{ℓ−1} = 0. (2.9)

Finally, apply the raising operator D_{ℓ}^{+} = −∂_{r} + ℓ/r to the last equation, thereby finding

−(D_{ℓ}^{+}Ψ_{ℓ−1})_{tt} + (D_{ℓ}^{+}Ψ_{ℓ−1})_{rr} − \frac{ℓ(ℓ + 1)}{r^{2}}D_{ℓ}^{+}Ψ_{ℓ−1} = 0. (2.10)

These steps prove the lemma.

2.3 Macdonald function as an outgoing solution

The Macdonald function (see p. 374 in [1], also [2]) K_{ν}(z) is a solution to the modified Bessel equation

z^{2}w′′ + zw′ − (z^{2} + ν^{2})w = 0, (2.11)

and it plays a prominent role in the representation of outgoing solutions to the ordinary wave equation. As an example, Ψ_{ℓ,s}(t, r) = e^{st}r^{1/2}K_{ℓ+1/2}(sr) obeys (2.6) and this solution is the one in (2.7) with f(t − r) = e^{s(t−r)} for the underlying profile function. Of course this does not have compact support and is not used in the following.

For half-integer order, we have [1]

K_{ℓ+1/2}(z) = \sqrt{\frac{π}{2z}} e^{−z}W_{ℓ}(z),  W_{ℓ}(z) = ∑_{k=0}^{ℓ} \frac{c_{ℓk}}{z^{k}} ,  c_{ℓk} = \frac{1}{2^{k}k!} \frac{(ℓ + k)!}{(ℓ − k)!},  (2.12)

and the first few Macdonald functions are (see 10.2.15-17 in [1])

\begin{aligned}+K_{3/2}(z) & = \sqrt{\frac{π}{2z}} e^{−z} (1 + \frac{1}{z}) \\ ⋄K_{5/2}(z) & = \sqrt{\frac{π}{2z}} e^{−z} (1 + \frac{3}{z} + \frac{3}{z^{2}} ) \\ ∘K_{7/2}(z) & = \sqrt{\frac{π}{2z}} e^{−z} (1 + \frac{6}{z} + \frac{15}{z^{2}} + \frac{15}{z^{3}}) \\ ∗K_{9/2}(z) & = \sqrt{\frac{π}{2z}} e^{−z} (1 + \frac{10}{z} + \frac{45}{z^{2}} + \frac{105}{z^{3}} + \frac{105}{z^{4}}) . \end{aligned}  (2.13)

The K_{ℓ+1/2}(z) can also be expressed as

K_{ℓ+1/2}(z) = \sqrt{\frac{π}{2z}} \frac{e^{−z}}{z^{ℓ}} p_{ℓ}(z),  p_{ℓ}(z) = ∑_{k=0}^{ℓ} c_{ℓk}z^{ℓ−k},  (2.14)

where here we refer to the polynomial p_{ℓ}(z) as the Bessel polynomial.

Let the roots of p_{ℓ}(z) be denoted by {b_{ℓj} : j = 1, …, ℓ}. Watson [2] shows that each b_{ℓj} lies in the left-half plane and is simple. Watson’s analysis is presented in the Appendix.

Figure 2.1: Scaled roots of K_{ℓ+1/2}(z).

The set {b_{ℓj}/(ℓ + 1/2) : j = 1, …, ℓ} is the collection of roots scaled by the Bessel order ν = ℓ + 1/2. For ℓ = 1, 2, 3, 4 the scaled roots are depicted in Fig. 2.1 using the symbols +, ⋄, ∘, ∗ respectively. Olver [4] showed that for large ℓ the scaled roots tend to lie on the curve given by

z(λ) = −\sqrt{λ^{2} − λ tanh λ} ± i\sqrt{λ coth λ − λ^{2}}

for λ ∈ [0, λ_{0}], where λ_{0} ≃ 1.19967864025773 solves tanh λ_{0} = 1/λ_{0}. Note that z(0) = ±i, giving the two points at the end of the arc. To the eye the roots are on this curve for even small ℓ, which also illustrates that often asymptotics are good when the parameter is not large.

Nonreflecting Boundary Conditions

Our goal in this section is to specify an exact transparent BC for the wave equation in the presence of a spherical boundary. Again, we will work “multipole by multipole” and specify a condition obeyed by (2.7). Technically, the analysis will be easier if we work with a spherical annulus centered at r = 0, thereby avoiding issues at the origin r = 0 which are not relevant for the outer BCs we are considering here. Therefore, we introduce inner r_{0} and outer r_{B} radii (B for “boundary”).

In this chapter we prove the following. Theorem 3.0.1. Consider 1 ≤ r_{0} < r_{B} and the associated interval D = [−r_{B} + δ, −r_{0} − δ], where 0 < δ < (r_{B} − r_{0})/2. Suppose that f ∈ C_{0}^{∞}(D), and consider the “wave field”

Ψ_{ℓ}(t, r) = ∑_{k=0}^{ℓ} \frac{1}{r^{k}} c_{ℓk}f^{(ℓ−k)}(t − r). (3.1)

Then at r = r_{B} the wave field obeys

Ψ_{t}(t, r_{B}) + Ψ_{r}(t, r_{B}) = \frac{1}{r_{B}} ∫_{0}^{t} Ω_{ℓ}(t − t′, r_{B})Ψ_{ℓ}(t′, r_{B})dt′, (3.2)

where the kernel

Ω_{ℓ}(t, r) = ∑_{n=1}^{ℓ}(b_{ℓn}/r) exp(b_{ℓn}t/r) = ∂_{t} ∑_{n=1}^{ℓ} exp(b_{ℓn}t/r). (3.3)

Chapter 3. Nonreflecting Boundary Conditions

Figure 3.1: Radial axis with support of Ψ_{ℓ}(0, r) and point at r_{B} − ϵδ (the dashed line).

Figure 3.2: Retarded time axis and support of profile function f.

We care only about (3.2) at r_{B} but in fact it holds for r ≥ r_{B}.

3.1 Proof of Theorem 3.0.1

Our first proof of the theorem relies on the following lemma.

Lemma 3.1.1. For r ∈ (r_{B} − ϵδ, ∞) with ϵ ≤ 1 the Laplace transform Ψ̂_{ℓ}(s, r) of the expression Ψ_{ℓ}(t, r) given in (3.1) with f ∈ C_{0}^{∞}(D) has form [cf. Eq. (2.12)]

Ψ̂_{ℓ}(s, r) = a(s)s^{ℓ}e^{−sr}W_{ℓ}(sr),  a(s) = ∫_{−r_{B}+δ}^{−r_{0}−δ} e^{−su}f(u)du.  (3.4)

Notice that a(s) is independent of r, and that this holds for any complex number s.

Proof. To obtain (3.4), we proceed as follows. First consider the Laplace transform of f^{(ℓ−k)}(t − r) which appears in (3.1). Repeated integration by parts generates t = 0 boundary terms of the form f^{(ℓ−k−p)}(−r) for 1 ≤ p ≤ ℓ − k; all such terms vanish. Indeed, for r ∈ (r_{B} − ϵδ, ∞), we have −r < −r_{B} + δ and so −r ∉ D. Moreover, the profile function f and all of its derivatives vanish on ℝ\D. The repeated integration by parts also generates t = ∞ boundary terms which are similarly shown to vanish.

The previous argument establishes that

∫_{0}^{∞} e^{−st}f^{(ℓ−k)}(t − r)dt = s^{ℓ−k}e^{−sr} ∫_{−r}^{∞} e^{−su}f(u)du. (3.5)

The lower limit −r of integration can now be replaced with −r_{B} + δ, because f(u) = 0 for u ∈ [−r, −r_{B} + δ) which implies u ∉ D, and similarly the upper limit by −r_{0} − δ.

The proof of Theorem 3.0.1 follows from Lemma 3.1.1.

Proof. The explicit expression (3.4) for Ψ̂_{ℓ}(s, r) determines an exact frequency-domain boundary condition (the Dirichlet-to-Neumann map)

sΨ̂_{ℓ}(s, r) + ∂_{r}Ψ̂_{ℓ}(s, r) = \frac{1}{r}Ω̂_{ℓ}(s, r)Ψ̂_{ℓ}(s, r), (3.6)

where the frequency-domain radiation kernel Ω̂_{ℓ}(s, r) defines the Sommerfeld residual. Indeed, the operator on the left hand of (3.6) corresponds to the Sommerfeld operator ∂_{t} + ∂_{r} in the time-domain. If ℓ = 0, then sΨ̂_{ℓ}(s, r) + ∂_{r}Ψ̂_{ℓ}(s, r) = 0. In (3.6) Ω̂_{ℓ}(s, r) is essentially the logarithmic derivative of W_{ℓ}(sr),

Ω̂_{ℓ}(s, r) ≡ sr\frac{W′_{ℓ}(sr)}{W_{ℓ}(sr)} = ∑_{k=1}^{ℓ} \frac{b_{ℓk}/r}{s − b_{ℓk}/r}. (3.7)

The last equality arises from the following consideration. We recognize W_{ℓ}(z) as the expression given in (2.12). With this result and (3.6), we get (3.2) - (3.3) upon inverse Laplace transformation.

Remark: The proof of Lemma 3.1.1 shows that Theorem 3.0.1 actually holds for any Ψ_{ℓ}(t, r) of the form (3.1). That is the c_{ℓk} need not be the special expressions in (2.12). Of course for general c_{ℓk} the b_{ℓk} in Theorem 3.0.1 must be the corresponding roots of p_{ℓ}(z) = ∑_{k=1}^{ℓ} c_{ℓk}z^{ℓ−k}. Moreover, if the c_{ℓk} are not those given in (2.12), then Ψ_{ℓ}(t, r) will not generally solve the radial wave equation (2.6).

Let us now consider what the boundary kernels for a spherical boundary look like for the cases ℓ = 1 and ℓ = 2. For the case ℓ = 1, we refer to the general radial wave equation (2.6) to get (suppressing the ℓ = 1 subscript)

−s^{2}Ψ̂ + Ψ̂_{rr} − \frac{2}{r^{2}} Ψ̂ = 0 (3.8)

for the Laplace-transformed equation. In order to peel off the exponential dependence of the outgoing solution, we let Ψ̂ = e^{−sr}Φ̂ to obtain

Φ̂_{rr} − 2sΦ̂_{r} − \frac{2}{r^{2}} Φ̂ = 0. (3.9)

Letting z = sr gives us

Φ̂_{zz} − 2Φ̂_{z} − \frac{2}{z^{2}} Φ̂ = 0. (3.10)

Linearly independent solutions to this equation are

W_{1}(z) = 1 + \frac{1}{z}  Z_{1}(z) = e^{2z}(1 − \frac{1}{z}).  (3.11)

The product Ω_{1}(s, r) of z and the logarithmic derivative of W_{1} for z = sr then takes the form

Ω_{1}(s, r) = \frac{zW′_{1}(z)}{W_{1}(z)} = \frac{z(−\frac{1}{z^{2}} )}{1 + \frac{1}{z}} = −\frac{1}{z + 1}. (3.12)

Now, for the case ℓ = 2, we again refer to the general radial wave equation (2.6) to get (suppressing the ℓ = 2 subscript)

Φ̂_{rr} − 2sΦ̂_{r} − \frac{6}{r^{2}} Φ̂ = 0 (3.13)

for the Laplace-transformed equation. Letting z = sr gives us

Φ̂_{zz} − 2Φ̂_{z} − \frac{6}{z^{2}} Φ̂ = 0. (3.14)

The linearly independent solutions to this equation are then

W_{2}(z) = 1 + \frac{3}{z} + \frac{3}{z^{2}}  Z_{2}(z) = e^{2z}(1 − \frac{3}{z} + \frac{3}{z^{2}} ),  (3.15)

and the product Ω_{2}(s, r) of z and the logarithmic derivative of W_{2} for z = sr then takes the form

Ω_{2}(s, r) = \frac{zW′_{2}(z)}{W_{2}(z)} = \frac{z(−\frac{3}{z^{2}} − \frac{6}{z^{3}} )}{1 + \frac{3}{z} + \frac{3}{z^{2}}} = \frac{−3z − 6}{z^{2} + 3z + 3}. (3.16)

3.2 Alternative proof.

This subsection sketches an outline of a different formal proof of the convolution result (3.2) in Theorem 3.0.1 by solving (2.6) on (r_{0}, ∞) as an initial value problem assuming (i) Ψ(t, r_{0}) = 0 and (ii) data of compact support on [r_{0} +δ, r_{B} −δ], 0 < δ < (r_{B} −r_{0})/2. Unlike the previous proof, this one relies on the wave equation and so assumes ultimately that the c_{ℓk} are those given in (2.12). However, this proof allows for different initial data than that associated with (3.1). In the proof we continue to suppress the ℓ on Ψ.

Proof. The Laplace transform of (2.6) is

[\frac{d^{2}}{dr^{2}} − \frac{ℓ(ℓ + 1)}{r^{2}} − s^{2}]Ψ̂(s, r) = −Ψ̇(0, r) − sΨ(0, r). (3.17)

Since Ψ(t, r_{0}) = 0, we have Ψ̂(s, r_{0}) = 0. Furthermore, we require Ψ̂(s, r) ∼ e^{−sr} as r → ∞ (the latter “boundary condition at ∞” is the outgoing assumption).

We solve (3.17) subject to the specified boundary conditions via the method of Green’s functions. The following are solutions to the homogeneous equation [i.e. (3.17) with the right-hand side set to zero]:

\sqrt{2πsr}I_{ℓ+1/2}(sr),  \sqrt{\frac{2sr}{π}}K_{ℓ+1/2}(sr),  (3.18)

where I_{ν}(z) is the other modified Bessel function. Choosing linearly independent solutions u, v (for Re(s) > 0) such that u(r_{0}) = 0 and v(r) ∼ e^{−sr} we obtain

u(r) = I_{ℓ+1/2}(r_{0}s)\sqrt{\frac{πr}{2s}}K_{ℓ+1/2}(sr) − K_{ℓ+1/2}(r_{0}s)\sqrt{\frac{πr}{2s}}I_{ℓ+1/2}(sr),  (3.19)
v(r) = \sqrt{\frac{2sr}{π}}K_{ℓ+1/2}(sr),  (3.20)

where of course these solutions also depend on s. Clearly, u(r_{0}) = 0 and v(r) ∼ e^{−sr} as r → ∞. The Wronskian W(r) = u(r)v′(r) − v(r)u′(r) is then

\begin{aligned}W(r) & = −srK_{ℓ+1/2}(sr_{0}) [I_{ℓ+1/2}(sr)K′_{ℓ+1/2}(sr) − K_{ℓ+1/2}(sr)I′_{ℓ+1/2}(sr)] \\ & = −srK_{ℓ+1/2}(sr_{0}) (−\frac{1}{sr}) = K_{ℓ+1/2}(sr_{0}), \end{aligned}  (3.21)

where the prime denotes differentiation in the argument. The boundary value problem could be solved using variation of parameters, however we take a Green’s function approach.

The Green’s function G(r, ξ; s) obeys

[\frac{d^{2}}{dr^{2}} − \frac{ℓ(ℓ + 1)}{r^{2}} − s^{2}] G(r, ξ; s) = −δ(r − ξ), (3.22)

as well as the boundary conditions G(r_{0}, ξ; s) = 0 and G(r, ξ; s) ∼ e^{−sr} as r → ∞. The Green’s function is then [5]

G(r, ξ; s) = \{ \begin{matrix}−\frac{u(r)v(ξ)}{K_{ℓ+1/2}(sr_{0})} & for r_{0} < r ≤ ξ \\ −\frac{u(ξ)v(r)}{K_{ℓ+1/2}(sr_{0})} & for ξ ≤ r < ∞. \end{matrix} (3.23)

The solution Ψ̂(s, r) to the boundary problem (3.17) is then the spatial convolution

Ψ̂(s, r) = ∫_{r_{0}+δ}^{r_{B}−δ} G(r, ξ, s)[Ψ̇(0, ξ) + sΨ(0, ξ)]dξ. (3.24)

For r > r_{B} − δ, we have

\begin{aligned}Ψ̂(s, r) & = − ∫_{r_{0}+δ}^{r_{B}−δ} \frac{u(ξ)v(r)}{K_{ℓ+1/2}(sr_{0})}[Ψ̇(0, ξ) + sΨ(0, ξ)]dξ \\ & = −\frac{v(r)}{K_{ℓ+1/2}(sr_{0})} ∫_{r_{0}+δ}^{r_{B}−δ} u(ξ)[Ψ̇(0, ξ) + sΨ(0, ξ)]dξ \\ & = −\frac{v(r)}{K_{ℓ+1/2}(sr_{0})}A(s). \end{aligned}  (3.25)

Using (2.12) and (3.20) then gives

\begin{aligned}Ψ̂(s, r) & = −A(s)\sqrt{\frac{2sr}{π}} \frac{K_{ℓ+1/2}(sr)}{K_{ℓ+1/2}(sr_{0})} \\ & = −A(s)\sqrt{\frac{2sr}{π}} \frac{1}{K_{ℓ+1/2}(sr_{0})}\sqrt{\frac{π}{2sr}}e^{−sr}W_{ℓ}(sr) \\ & = b(s)e^{−sr}W_{ℓ}(sr). \end{aligned}  (3.26)

Chapter 3. Nonreflecting Boundary Conditions Here

b(s) = −\frac{1}{K_{ℓ+1/2}(sr_{0})} ∫_{r_{0}+δ}^{r_{B}−δ} u(ξ)[Ψ̇(0, ξ) + sΨ(0, ξ)]dξ. (3.27)

For general compact initial data b(s) cannot be expressed as s^{ℓ}a(s) as in (3.4). From this point the proof follows exactly the steps of the first proof after Lemma 3.1.1.

Newton’s Identities

4.1 Main result

Our goal is to prove the following result.

Theorem 4.1.1. The roots b_{ℓj} : j = 1, ..., ℓ of the ℓth degree polynomial

p_{ℓ}(z) = ∑_{k=0}^{ℓ} c_{ℓk}z^{ℓ−k} (4.1)

also obey the following set of ℓ algebraic equations

−kc_{ℓk} = ∑_{n=1}^{ℓ} b_{ℓn} ∑_{q=1}^{k} c_{ℓ,q−1}(b_{ℓn})^{k−q},  k = 1, ..., ℓ,  (4.2)

where here we assume ℓ ≥ 1.

Remark: All of the factors in this system (roots b_{ℓj} and coefficients c_{ℓj}) carry a lead “name” index ℓ. We now suppress this index, and reorganize the last expressions using c_{ℓ0} ≡ c_{0} = 1 and the re-indexing p = q − 1 to reach

−kc_{k} = ∑_{n=1}^{ℓ}(b_{n})^{k} + ∑_{p=1}^{k−1} c_{p} ∑_{n=1}^{ℓ}(b_{n})^{k−p},  k = 1, ..., ℓ.  (4.3)

Chapter 4. Newton’s Identities

These formulas involve the power sums s_{k} = ∑_{n=1}^{ℓ}(b_{n})^{k} of the roots. In terms of the s_{k} and a_{ℓ−k} ≡ c_{k}, we then have

−ka_{ℓ−k} = s_{k} + ∑_{p=1}^{k−1} a_{ℓ−p}s_{k−p},  k = 1, ..., ℓ,  (4.4)

which we recognize as a subset of Newton’s identities [6].

In what remains of this section, we prove Theorem 4.1.1 assuming Theorem 3.0.1 which was established in the previous chapter. The idea of the proof is fairly straightforward; using the expressions (3.1) and (3.2), we find the identity

− ∑_{k=1}^{ℓ} \frac{k}{r_{B}^{k+1}} c_{ℓk}f^{(ℓ−k)}(t−r_{B}) = \frac{1}{r_{B}} ∫_{0}^{t} Ω_{ℓ}(t−t′, r_{B})[ ∑_{k=0}^{ℓ} \frac{1}{r^{k}} c_{ℓk}f^{(ℓ−k)}(t′ −r_{B})]dt′, t > 0.

(4.5)

The proof amounts to a careful comparison of both sides in the last equation.

Lemma 4.1.2. Let r = r_{B} and f ∈ C_{0}^{∞}(D) as in Theorem 3.0.1, and define

I^{(p)}[b, r, f] = ∫_{0}^{t} e^{b(t−t′)/r}f^{(p)}(t′ − r)dt′. (4.6)

Then, we have

I^{(p)}[b, r, f] = (b/r)^{p}I^{(0)}[b, r, f] + ∑_{j=1}^{p}(b/r)^{p−j}f^{(j−1)}(t − r). (4.7)

Proof. We prove the lemma by induction, showing that

I^{(p)}[b, r, f] = (b/r)^{p}I^{(0)}[b, r, f] + ∑_{j=1}^{p}(b/r)^{p−j}[f^{(j−1)}(t − r) − e^{bt/r}f^{(j−1)}(−r)], (4.8)

where the second term within the square brackets vanishes since r = r_{B} ∉ D. Integration by parts establishes the last formula for p = 1. Similarly, integration by parts yields

I^{(p)}[b, r, f] = (b/r)I^{(p−1)}[b, r, f] + f^{(p−1)}(t − r) − e^{bt/r}f^{(p−1)}(−r). (4.9)

Assuming now that (4.8) holds with p replaced by p − 1, we insert the p − 1 result into the last equation, thereby recovering (4.8) and verifying the induction step.

We now use Theorem 3.0.1 to prove Theorem 4.1.1.

4.2 Proof of Theorem 4.1.1

Proof. First, taking r = r_{B} and t > 0 and then appealing to the formulas (3.1) and (3.3), we get

\frac{1}{r} ∫_{0}^{t} Ω_{ℓ}(t − t′, r)Ψ_{ℓ}(t′, r)dt′ = ∑_{n=1}^{ℓ}(b_{ℓn}/r^{2}) ∫_{0}^{t} exp(b_{ℓn}(t − t′)/r)Ψ_{ℓ}(t′, r)dt′  (4.10)

= ∑_{n=1}^{ℓ}(b_{ℓn}/r^{2}) ∑_{k=0}^{ℓ} r^{−k}c_{ℓk}I^{(ℓ−k)}[b_{ℓn}, r, f].

We will now focus on reducing this last expression to a desirable form. Lemma 4.1.2 gives

I^{(ℓ−k)}[b_{ℓn}, r, f] = (b_{ℓn}/r)^{ℓ−k}I^{(0)}[b_{ℓn}, r, f] + ∑_{j=1}^{ℓ−k}(b_{ℓn}/r)^{ℓ−k−j}f^{(j−1)}(t − r). (4.11)

Combination of the last two equations yields

\frac{1}{r} ∫_{0}^{t} Ω_{ℓ}(t − t′, r)Ψ_{ℓ}(t′, r)dt′ = r^{−(ℓ+2)} ∑_{n=1}^{ℓ} b_{ℓn}I^{(0)}[b_{ℓn}, r, f] ∑_{k=0}^{ℓ} c_{ℓk}(b_{ℓn})^{ℓ−k}  (4.12)

+ ∑_{n=1}^{ℓ}(b_{ℓn}/r^{2}) ∑_{k=0}^{ℓ} r^{−k}c_{ℓk} ∑_{j=1}^{ℓ−k}(b_{ℓn}/r)^{ℓ−k−j}f^{(j−1)}(t − r).

Since ∑_{k=0}^{ℓ} c_{ℓk}(b_{ℓn})^{ℓ−k} = p_{ℓ}(b_{ℓn}) is the Bessel polynomial evaluated at one of its roots, the first term on the right-hand side of the last expression vanishes. Whence, up to now

\frac{1}{r} ∫_{0}^{t} Ω_{ℓ}(t − t′, r)Ψ_{ℓ}(t′, r)dt′ = ∑_{n=1}^{ℓ}(b_{ℓn}/r^{2}) ∑_{k=0}^{ℓ} r^{−k}c_{ℓk} ∑_{j=1}^{ℓ−k}(b_{ℓn}/r)^{ℓ−k−j}f^{(j−1)}(t − r). (4.13)

Within the sum over k, the sum over j is empty when k = ℓ (no boundary term is generated by the term in (3.1) corresponding to the undifferentiated profile function f). Therefore, here we may replace ∑_{k=0}^{ℓ} by ∑_{k=0}^{ℓ−1}. This replacement, followed by the reindexing q = k + 1, yields

\frac{1}{r} ∫_{0}^{t} Ω_{ℓ}(t−t′, r)Ψ_{ℓ}(t′, r)dt′ = ∑_{n=1}^{ℓ}(b_{ℓn}/r^{2}) ∑_{q=1}^{ℓ} r^{−(q−1)}c_{ℓ,q−1} ∑_{j=1}^{ℓ−q+1}(b_{ℓn}/r)^{ℓ−q−j+1}f^{(j−1)}(t−r).

(4.14)

Now, the double sum ∑_{q=1}^{ℓ} ∑_{j=1}^{ℓ−q+1}(terms)_{jq} is equivalent to ∑_{j=1}^{ℓ} ∑_{q=1}^{ℓ−j+1}(terms)_{jq}. As a result, we may group the two inner sums and make this exchange, thereby reaching

\frac{1}{r} ∫_{0}^{t} Ω_{ℓ}(t−t′, r)Ψ_{ℓ}(t′, r)dt′ = ∑_{n=1}^{ℓ} b_{ℓn} ∑_{j=1}^{ℓ} r^{j−(ℓ+2)} ∑_{q=1}^{ℓ−j+1} c_{ℓ,q−1}(b_{ℓn})^{ℓ−q−j+1}f^{(j−1)}(t−r). (4.15)

Finally, re-indexing by k = ℓ − j + 1 leads to

\frac{1}{r} ∫_{0}^{t} Ω_{ℓ}(t − t′, r)Ψ_{ℓ}(t′, r)dt′ = ∑_{n=1}^{ℓ} b_{ℓn} ∑_{k=1}^{ℓ} r^{−(k+1)} ∑_{q=1}^{k} c_{ℓ,q−1}(b_{ℓn})^{k−q}f^{(ℓ−k)}(t − r), (4.16)

which is the aforementioned desirable form.

Equation (4.16) can now be used to finish the proof. Recalling that we have set r = r_{B}, with (4.16) the identity (4.5) becomes

− ∑_{k=1}^{ℓ} kr^{−(k+1)}c_{ℓk}f^{(ℓ−k)}(t − r) = ∑_{k=1}^{ℓ} r^{−(k+1)} ∑_{n=1}^{ℓ} b_{ℓn} ∑_{q=1}^{k} c_{ℓ,q−1}(b_{ℓn})^{k−q}f^{(ℓ−k)}(t − r), (4.17)

and we may express that last equation as

0 = ∑_{k=1}^{ℓ} r^{−(k+1)}E_{k}f^{(ℓ−k)}(u),  E_{k} = kc_{ℓk} + ∑_{n=1}^{ℓ} b_{ℓn} ∑_{q=1}^{k} c_{ℓ,q−1}(b_{ℓn})^{k−q}.  (4.18)

Here u = t − r is retarded time.

Introducing the operator Q ≡ (∂_{t} + ∂_{r})r^{2}, we then use induction to show that

Q^{p} ∑_{k=1}^{ℓ} r^{−(k+1)}E_{k}f^{(ℓ−k)}(u) = ∑_{k=p+1}^{ℓ}(−1)^{p} \frac{(k − 1)!}{(k − p − 1)!}r^{−(k−p+1)}E_{k}f^{(ℓ−k)}(u).  (4.19)

Chapter 4. Newton’s Identities For p = ℓ − 1, this identity implies E_{ℓ} = 0. Then assuming E_{p+2} = ⋯ = E_{ℓ} = 0, we find

Q^{p} ∑_{k=1}^{ℓ} r^{−(k+1)}E_{k}f^{(ℓ−k)}(u) = (−1)^{p} \frac{p!}{r^{2}} E_{p+1}f^{(ℓ−p−1)}(u),  (4.20)

yielding E_{p+1} = 0. Therefore, backwards iteration p = ℓ − 1, ℓ − 2, …, 0 establishes that E_{k} = 0 for k = 1, …, ℓ.

Numerical root evaluation and

identity check for several values of ℓ

As a means of verifying the result from the previous chapter, we now proceed to check these identities numerically. One might ask why there is need to check Newton’s identities, classical formulas known for centuries. The point is that by checking them for the particular case of the Bessel polynomial we may (i) understand their conditioning in this setting and (ii) further confirm the accuracy of our numerically generated Macdonald roots. Using Mathematica, we first generate the roots of the polynomial p_{ℓ}(z) in (2.14) to as much precision as needed for several values of ℓ. These root sets for each value of ℓ are then input into the set of ℓ identities (4.2), specifically by evaluating the expression

kc_{ℓk} + ∑_{n=1}^{ℓ} b_{ℓn} ∑_{q=1}^{k} c_{ℓ,q−1}(b_{ℓn})^{k−q},  k = 1, ..., ℓ  (5.1)

Chapter 5. Numerical root evaluation and identity check for several values of ℓ

to the desired precision. The accuracy with which this expression (5.1) evaluates to zero can be seen as a verification of the identity (4.5). For each ℓ value, this generates a set of ℓ errors for each value of k in (4.2). The maximum absolute error of the identity in this set of ℓ errors is then found as |Err_{max}| = max{|Err_{k}|}, where Err_{k} is the error in the k^{th} identity for a particular value of ℓ. The accuracy to which this error is known is then output. One observes that the accuracy of the errors reduces as ℓ grows. The concluding section discusses this observation.

5.1 Root finding and identity verification for ℓ = 20 to

27 digits of precision

The Mathematica code and associated output follow, showing that the identity is indeed verified for these ℓ values at this level of precision. This code can be altered to increase both with satisfying results. Since complex roots appear as conjugate pairs, only the upper conjugate is listed for the sake of brevity.

(* set precision of calculations and highest l value *) (* set precision of calculations and highest l value *) (* set precision of calculations and highest l value *)
numdigits:=28 numdigits:=28 numdigits:=28
lmax:=20 lmax:=20 lmax:=20
$MaxExtraPrecision:=∞ $MaxExtraPrecision:=∞ $MaxExtraPrecision:=∞
(*define c {l, k}*) (*define c {l, k}*) (*define c {l, k}*)
cee[l , k ]:=(1/(2^{∧}k ∗ k!)) ∗ (l + k)!/(l − k)! cee[l , k ]:=(1/(2^{∧}k ∗ k!)) ∗ (l + k)!/(l − k)! cee[l , k ]:=(1/(2^{∧}k ∗ k!)) ∗ (l + k)!/(l − k)!
(* define polynomial to find roots of *) (* define polynomial to find roots of *) (* define polynomial to find roots of *)
polyroots:=NSolve[Sum[cee[l, k] ∗ z^{∧}(l − k), {k, 0, l}] == 0, z, numdigits] polyroots:=NSolve[Sum[cee[l, k] ∗ z^{∧}(l − k), {k, 0, l}] == 0, z, numdigits] polyroots:=NSolve[Sum[cee[l, k] ∗ z^{∧}(l − k), {k, 0, l}] == 0, z, numdigits]
(* define identity from main theorem *) (* define identity from main theorem *) (* define identity from main theorem *)
identity:=k ∗ cee[l, k] + SetPrecision[ identity:=k ∗ cee[l, k] + SetPrecision[ identity:=k ∗ cee[l, k] + SetPrecision[
Sum[Part[{z}/.polyroots, n] ∗ Sum[cee[l, q − 1] ∗ Part[{z}/.polyroots, n]^{∧}(k − q), {q, 1, k}], Sum[Part[{z}/.polyroots, n] ∗ Sum[cee[l, q − 1] ∗ Part[{z}/.polyroots, n]^{∧}(k − q), {q, 1, k}], Sum[Part[{z}/.polyroots, n] ∗ Sum[cee[l, q − 1] ∗ Part[{z}/.polyroots, n]^{∧}(k − q), {q, 1, k}],
{n, 1, l}], numdigits] {n, 1, l}], numdigits] {n, 1, l}], numdigits]

(* loop through each value of l *) (* loop through each value of l *) (* loop through each value of l *)
For[l = 1, l < lmax + 1, l++, For[l = 1, l < lmax + 1, l++, For[l = 1, l < lmax + 1, l++,
Print[] Print[] Print[]
Print[“real and imaginary parts of roots for l=”, l] Print[“real and imaginary parts of roots for l=”, l] Print[“real and imaginary parts of roots for l=”, l]
Print[] Print[] Print[]

(* print out formatted tables of polynomial roots *) (* print out formatted tables of polynomial roots *) (* print out formatted tables of polynomial roots *)
For[j = 1, j < l + 1, j++ For[j = 1, j < l + 1, j++ For[j = 1, j < l + 1, j++
If[Im[Part[Part[{z}/.polyroots, j − 1], 1]] < 0, , If[Im[Part[Part[{z}/.polyroots, j − 1], 1]] < 0, , If[Im[Part[Part[{z}/.polyroots, j − 1], 1]] < 0, ,
Print[ScientificForm[Re[Part[Part[{z}/.polyroots, j − 1], 1]], numdigits − 1, Print[ScientificForm[Re[Part[Part[{z}/.polyroots, j − 1], 1]], numdigits − 1, Print[ScientificForm[Re[Part[Part[{z}/.polyroots, j − 1], 1]], numdigits − 1,
NumberFormat → NumberFormat → NumberFormat →
(Row[{#1, E, If[#3 == “”, “+000”, If[Part[Characters[#3], 1] == “-”, (Row[{#1, E, If[#3 == “”, “+000”, If[Part[Characters[#3], 1] == “-”, (Row[{#1, E, If[#3 == “”, “+000”, If[Part[Characters[#3], 1] == “-”,
StringJoin[“-00”, StringDrop[#3, 1]], StringJoin[“+00”, #3]]]}]&)], StringJoin[“-00”, StringDrop[#3, 1]], StringJoin[“+00”, #3]]]}]&)], StringJoin[“-00”, StringDrop[#3, 1]], StringJoin[“+00”, #3]]]}]&)],
“ ”, ScientificForm[Im[Part[Part[{z}/.polyroots, j − 1], 1]], numdigits − 1, “ ”, ScientificForm[Im[Part[Part[{z}/.polyroots, j − 1], 1]], numdigits − 1, “ ”, ScientificForm[Im[Part[Part[{z}/.polyroots, j − 1], 1]], numdigits − 1,
NumberFormat → NumberFormat → NumberFormat →
(Row[{#1, E, If[#3 == “”, “+000”, If[Part[Characters[#3], 1] == “-”, (Row[{#1, E, If[#3 == “”, “+000”, If[Part[Characters[#3], 1] == “-”, (Row[{#1, E, If[#3 == “”, “+000”, If[Part[Characters[#3], 1] == “-”,
StringJoin[“-00”, StringDrop[#3, 1]], StringJoin[“+00”, #3]]]}]&)]] StringJoin[“-00”, StringDrop[#3, 1]], StringJoin[“+00”, #3]]]}]&)]] StringJoin[“-00”, StringDrop[#3, 1]], StringJoin[“+00”, #3]]]}]&)]]
]] ]] ]]

(* print error norm of all roots for a given value of l *) (* print error norm of all roots for a given value of l *) (* print error norm of all roots for a given value of l *)
Print[] Print[] Print[]
Print[“identity error accuracy ∼ ”, Print[“identity error accuracy ∼ ”, Print[“identity error accuracy ∼ ”,
ScientificForm[Max[Abs[Part[Part[Delete[Reap[For[k = 1, k < l + 1, k++, Sow[identity]]; ScientificForm[Max[Abs[Part[Part[Delete[Reap[For[k = 1, k < l + 1, k++, Sow[identity]]; ScientificForm[Max[Abs[Part[Part[Delete[Reap[For[k = 1, k < l + 1, k++, Sow[identity]];
], 1], 1], 1]]], numdigits, ], 1], 1], 1]]], numdigits, ], 1], 1], 1]]], numdigits,
NumberFormat → NumberFormat → NumberFormat →
(Row[{“1.”, E, If[#3 == “”, “+000”, If[Part[Characters[#3], 1] == “-”, (Row[{“1.”, E, If[#3 == “”, “+000”, If[Part[Characters[#3], 1] == “-”, (Row[{“1.”, E, If[#3 == “”, “+000”, If[Part[Characters[#3], 1] == “-”,
StringJoin[“-0”, StringDrop[#3, 1]], StringJoin[“+0”, #3]]]}]&)]]] StringJoin[“-0”, StringDrop[#3, 1]], StringJoin[“+0”, #3]]]}]&)]]] StringJoin[“-0”, StringDrop[#3, 1]], StringJoin[“+0”, #3]]]}]&)]]]

Table 5.1: Roots of K_{ℓ+1/2}(z) for ℓ = 1, ..., 20

real and imaginary parts of roots for l=1
-1.00000000000000000000000000E+000  0E+000
identity error accuracy ~ 1.E-028
real and imaginary parts of roots for l=2
-1.50000000000000000000000000E+000  8.66025403784438646763723171E-001
identity error accuracy ~ 1.E-027
real and imaginary parts of roots for l=3
-2.32218535462608559291147071E+000  0E+000
-1.83890732268695720354426464E+000  1.75438095978372166095183060E+000
identity error accuracy ~ 1.E-026

real and imaginary parts of roots for l=4
-2.89621060282037216839439393E+000  8.67234128934503751818973215E-001
-2.10378939717962783160560607E+000  2.65741804185675271685832210E+000
identity error accuracy ~ 1.E-025
real and imaginary parts of roots for l=5
-3.64673859532964325973516964E+000  0E+000
-3.35195639915353314301641425E+000  1.74266141618319772272632347E+000
-2.32467430318164522711600093E+000  3.57102292033797640038610261E+000
identity error accuracy ~ 1.E-024
real and imaginary parts of roots for l=6
-4.24835939586336394493607975E+000  8.67509673231365606386441571E-001
-3.73570835632581466794136197E+000  2.62627231144712564049355159E+000
-2.51593224781082138712255828E+000  4.49267295365394253591824392E+000
identity error accuracy ~ 1.E-023
real and imaginary parts of roots for l=7
-4.97178685852793567786117785E+000  0E+000
-4.75829052815462894523746281E+000  1.73928606113053654289381190E+000
-4.07013916363813747170592771E+000  3.51717404770975316581518078E+000
-2.68567687894326574412602056E+000  5.42069413071674889584849422E+000
identity error accuracy ~ 1.E-022

real and imaginary parts of roots for l=8
-5.58788604326308519899909296E+000  8.67614445352786459816302823E-001
-5.20484079063688191825036523E+000  2.61617515264252742873877731E+000
-4.36828921720240240703092253E+000  4.41444250047153908355026146E+000
-2.83898394889763047571961928E+000  6.35391129860487682208469639E+000
identity error accuracy ~ 1.E-020

real and imaginary parts of roots for l=9
-6.29701918171496853775919603E+000  0E+000
-6.12936790427427278778116357E+000  1.73784838348086250370031346E+000
-5.60442181950778141621337163E+000  3.49815691788609357601465626E+000
-4.63843988718039029665812495E+000  5.31727167543565114230342731E+000
-2.97926079818007123046774183E+000  7.29146368834218209076005832E+000
identity error accuracy ~ 1.E-019

real and imaginary parts of roots for l=10
-6.92204490542724611542320812E+000  8.67665195451221438446134323E-001
-6.61529096547687025894515467E+000  2.61156792080008987349400899E+000
-5.96752832858778584036140538E+000  4.38494718894193206818405874E+000
-4.88621956685899957989604215E+000  6.22498548247156710752681539E+000
-3.10891623364909820537418967E+000  8.23269945907358745046141002E+000
identity error accuracy ~ 1.E-018

real and imaginary parts of roots for l=11
-7.62233984579642949207664512E+000  0E+000
-7.48422986073193879164098577E+000  1.73710282075340382415004181E+000
-7.05789238766995313732999335E+000  3.48901450355582971973580520E+000
-6.30133745487130838170915958E+000  5.27619174369676828439756131E+000
-5.11564828390827892036268160E+000  7.13702075889336672017633185E+000
-3.22972208992030602291885715E+000  9.17711156870857844433181315E+000
identity error accuracy ~ 1.E-016

real and imaginary parts of roots for l=12
-8.25342201141208096388242751E+000  8.67693572009768815999260468E-001
-7.99727059960143416395771680E+000  2.60906653694579815898358031E+000
-7.46557124035177041635434551E+000  4.37016959335456516328553217E+000
-6.61100424995635186370036445E+000  6.17153499303722969677556413E+000
-5.32970859087582917461993816E+000  8.05290686425703263930869928E+000
-3.34302330780253341748520757E+000  1.01242968072408198876735301E+001
identity error accuracy ~ 1.E-015

real and imaginary parts of roots for l=13
-8.94770967439179101800795175E+000  0E+000
-8.83025208414490416096754110E+000  1.73666640030763054334032015E+000
-8.47059177147718456305833984E+000  3.48386845066099304960125471E+000
-7.84438027706259622338290970E+000  5.25490340661196188422886280E+000
-6.90037282614665984409350267E+000  7.07064431215294895763843768E+000
-5.53068098334403653612614328E+000  8.97224777515578769951113656E+000
-3.44986722062872316336758754E+000  1.10739285522161970404484911E+001
identity error accuracy ~ 1.E-013

real and imaginary parts of roots for l=14
-9.58317139364696666545939559E+000  8.67711028864253165359599786E-001
-9.36314585160955222518717854E+000  2.60755332438166661179716945E+000
-8.91100055537504426431713428E+000  4.36160417830244736590463016E+000
-8.19884696998847462160089744E+000  6.14304107147079673955607766E+000
-7.17239596217181731547627714E+000  7.97321735418496849451244718E+000
-5.72035238382751889004599461E+000  9.89470759748915927072561405E+000
-3.55108688338062601791312240E+000  1.20257380322545244924387964E+001
identity error accuracy ~ 1.E-012

real and imaginary parts of roots for l=15
-1.02731096663224778078893516E+001  0E+000
-1.01709139964400681778642607E+001  1.73638891945045594842184913E+000
-9.85956722839627922331631815E+000  3.48067121143276632621161317E+000
-9.32359932060896942898832732E+000  5.24225889523761694360563458E+000
-8.53245905229834060184246438E+000  7.03439362551704559270114646E+000
-7.42939699294215379659889775E+000  8.87898262112151570044737878E+000
-5.90015171366464763005831176E+000  1.08199991377535732063444092E+001
-3.64735686248830223738674416E+000  1.29795010707604194609117519E+001
identity error accuracy ~ 1.E-011

real and imaginary parts of roots for l=16
-1.09118860786775030650439450E+001  8.67722527435720452392272192E-001
-1.07189858189780134799231120E+001  2.60656700725828947911621149E+000
-1.03251196023414622533772726E+001  4.35616338060960859540158749E+000
-9.71232633256350378195566778E+000  6.12576089102176697237096423E+000
-8.84796819650278498065891895E+000  7.92877285588937164588993504E+000
-7.67324079086715998439086145E+000  9.78769743836906830811183539E+000
-6.07124138290869981857329635E+000  1.17478749384808882640804452E+001
-3.73923179716087263607692588E+000  1.39350284758133825846001679E+001
identity error accuracy ~ 1.E-09

real and imaginary parts of roots for l=17
-1.15985294923395501618999203E+001  0E+000
-1.15080767771397596143175561E+001  1.73620153790806340947249102E+000
-1.12334368172695439283853762E+001  3.47854389076469685688536076E+000
-1.07641341775628432905926096E+001  5.23407490203687631572176595E+000
-1.00802944448577807095976715E+001  7.01200998269376838973610942E+000
-9.14758867760315452597282895E+000  8.82599830149333386118061796E+000
-7.90544959593734230266452814E+000  1.06991450754651671694245213E+001
-6.23458097836041372958967294E+000  1.26781202290665046041369318E+001
-3.82717378509938681792979640E+000  1.48921589246642892796739770E+001
identity error accuracy ~ 1.E-07

real and imaginary parts of roots for l=18
-1.22399021368250294751608712E+001  8.67730500530609374492680708E-001
-1.20681358449367730276078721E+001  2.60588788173345421215154585E+000
-1.17189487956552908233983590E+001  4.35247975429983523855180469E+000
-1.11800390165370403880516547E+001  6.11439409303699536109075213E+000
-1.04300129653021452464250164E+001  7.90089310331303521101324559E+000
-9.43313222080871295101261560E+000  9.72590031412845700845172164E+000
-8.12728394509562492375307106E+000  1.16131317511959929859269619E+001
-6.39097278368397486896633901E+000  1.36105473490914326733665756E+001
-3.91157229115540829562420106E+000  1.58507535969377336782227290E+001
identity error accuracy ~ 1.E-06

real and imaginary parts of roots for l=19
-1.29239630554237278712291506E+001  0E+000
-1.28428277968952223042629424E+001  1.73606905090273320111409963E+000
-1.25970628096647620973052159E+001  3.47705490010667713280513635E+000
-1.21792312603829384923664555E+001  5.22845054839089777949190898E+000
-1.15756010654031831963984511E+001  6.99707637470125717969882230E+000
-1.07635384400032784129991877E+001  8.79229302167302653136215767E+000
-9.70610240075824858168220378E+000  1.06283211002462878627136475E+001
-8.33980071913673608260638162E+000  1.25294838239446239254763941E+001
-6.54109506216141403075300463E+000  1.45449913030212108164540177E+001
-3.99275891788235286601158213E+000  1.68106920601116228524404469E+001
identity error accuracy ~ 1.E-04

real and imaginary parts of roots for l=20
-1.35674242831533132904409805E+001  8.67736254955798277563901221E-001
-1.34125971436066018984302504E+001  2.60540014717949754240528038E+000
-1.30988224745771632490213023E+001  4.34986491179146240148265878E+000
-1.26172813166098517250117887E+001  6.10647987005239646124623649E+000
-1.19530908024999872535298649E+001  7.88205843424744989157438123E+000
-1.10825803337311520272082869E+001  9.68609324182857851119637990E+000
-9.96776247886039112461620511E+000  1.15331147285162466386412696E+001
-8.54389572685003190873213213E+000  1.34480452734196997247112442E+001
-6.68552687829519020092395369E+000  1.54813061879236185543330417E+001
-4.07101856181631732208523537E+000  1.77718690688854561225973949E+001
identity error accuracy ~ 1.E-03

Conclusion

We conclude with some comments on the significance of the analysis done here as it relates to the larger context. The main result can be viewed from two perspectives regarding Newton’s identities (4.2) and their relation to the convolution result (3.2).

Recall that the first proof of Theorem 3.0.1 does not rely on the particular form of the coefficients c_{ℓk} in (2.12). See Lemma 3.1.1 and the following remark. Thus, Chapter 4 may be viewed as a derivation of general Newton’s identities (or at least a subset thereof, see [6]), with the derivation assuming the convolution result provided earlier in Theorem 3.0.1. Alternatively, by first assuming Newton’s identities, the analysis in Chapter 4 can also be viewed as a time-domain proof of the convolution result in Theorem 3.0.1 without the need for the Laplace transform.

Newton’s identities for the particular case of the roots and coefficients of the Macdonald function appear ill-conditioned. In other words, they seem difficult to verify accurately for large values of ℓ. This is most likely due to the coefficients c_{ℓk} which for a fixed large ℓ (ℓ = 20 for example) vary over a large range, being factorially defined as they are in (2.12). This combined with the relative proximity of the roots to the origin (see Fig 2.1) leads to a necessity of many digits of precision for accurate numerical verification.

Appendix A

Location and count of Macdonald

roots

For any positive ν ≥ 0 Watson [2] has shewn that the roots of K_{ν}(z) are simple and lie in the left-half z plane. This appendix presents some of his analysis, often tailoring it to our specific case ν = ℓ + 1/2, that is half-integer order. In what follows we often follow [2] quite closely. Since K_{ν}(z) obeys the differential equation (2.11), any root b ≠ 0 must be simple. Indeed, were b ≠ 0 a root of multiplicity, then we could conclude that K′′_{ν}(b) = K′_{ν}(b) = K_{ν}(b) = 0, and by repeated differentiation of (2.11) that K_{ν}^{(p)}(b) = 0 for all p. Since all nonzero points of the Bessel equation are regular points, this would imply K_{ν}(z) = 0 in a neighborhood of b. Additionally, the small z asymptotic expansion for K_{ν}(z) (see [1]) shows that z = 0 is a singular point of K_{ν}(z) and not a root.

We begin by noting that the generalization of Bessel’s integral (see 9.6.24 [1])

K_{ν}(z) = ∫_{0}^{∞} e^{−z cosh(t)} cosh(νt)dt,  |arg z| < \frac{π}{2},  (A.1)

makes it clear that K_{ν}(z) has no positive real roots. Watson has extended this result to show that K_{ν}(z) has no roots such that |arg z| ≤ \frac{π}{2}.

Appendix A. Location and count of Macdonald roots

In order to demonstrate this, we must first consider a result of Macdonald (and proven in [2]) which represents the product K_{ν}(Z)K_{ν}(z) as an integral involving a single function of type K_{ν}:

K_{ν}(Z)K_{ν}(z) = \frac{1}{2} ∫_{0}^{∞}exp[ − \frac{v}{2} − \frac{Z^{2} + z^{2}}{2v} ]K_{ν}(\frac{Zz}{v} )\frac{dv}{v} , (A.2)

which is valid for all values of ν ≥ 0 when |arg Z| ≤ π, |arg z| ≤ π, and |arg (Z + z)| ≤ \frac{π}{4}. Now, if z = re^{−iα} were such a root (r > 0, −\frac{π}{2} < α < \frac{π}{2}), then Z = re^{iα} would be another root by (A.1). However, the integral shows us that

K_{ν}(re^{iα})K_{ν}(re^{−iα}) = \frac{1}{2} ∫_{0}^{∞}exp[ − \frac{v}{2} − \frac{r^{2} cos(2α)}{v}]K_{ν}(\frac{r^{2}}{v} )\frac{dv}{v}

> 0, (A.3)

which is contrary to the hypothesis that these are roots of K_{ν}(z). If α = ±\frac{π}{2}, then we have (again, see [2])

|K_{ν}(re^{± \frac{π}{2} i})| = \frac{π}{2} \sqrt{J_{ν}^{2}(r) + Y_{ν}^{2}(r)}, (A.4)

and so K_{ν}(z) has no purely imaginary roots, because J_{ν}(r) and Y_{ν}(r) cannot both vanish simultaneously.

A.1 Root count for K_{ℓ+1/2}(z)

For ν = ℓ + 1/2 (half-integer order) we have formula (2.13). Since p_{ℓ}(0) = c_{ℓℓ} ≠ 0 and \sqrt{π/2}e^{−z}z^{−(ℓ+2)} is singular at z = 0, then the origin is clearly not a root. This is an explicit check of the general result stated above. Moreover, from (2.14) we conclude that K_{ℓ+1/2}(z) has ℓ roots, since p_{ℓ}(z) is degree ℓ. Notice that p_{ℓ}(z) may have a root on the negative real axis, which is technically not in the domain of analyticity of K_{ℓ+1/2}(z) due to the branch associated with z^{−(ℓ+2)} in the prefactor \sqrt{π/2}e^{−z}z^{−(ℓ+2)}. Nevertheless, we shall consider such a root of p_{ℓ}(z) also as a root of K_{ℓ+1/2}(z). With the above analysis, we conclude that K_{ℓ+1/2}(z) has ℓ simple roots, each lying in the left-half z plane.

A.2 Root count for K_{ν}(z)

Now we turn our attention to general ν > 0 and the roots for which Re(z) < 0 and with \frac{π}{2} < |arg z| < π. We show that the number of roots in these two quadrants is the even integer nearest to ν − \frac{1}{2}, unless ν − \frac{1}{2} itself is an integer, in which case (as seen above) it is also the number of roots.

If ν − \frac{1}{2} is not an integer, then there are no roots on the negative real axis, where arg z = ±π. This follows from

K_{ν}(re^{±πi}) = e^{∓νπi}K_{ν}(r) ∓ πiI_{ν}(r), (A.5)

and for both the real and imaginary parts to vanish, then we must have

cos(νπ)K_{ν}(r) = 0, sin(νπ)K_{ν}(r) + πI_{ν}(r) = 0. (A.6)

The Wronskian of the functions on the left of the equations is \frac{π}{r}cos(νπ). Therefore, K_{ν}(r) and I_{ν}(r) cannot vanish simultaneously unless cos(νπ) = 0. We now use the Argument principle to count the roots in the second and third quadrants.

Now we consider the change in phase of z^{ν}K_{ν}(z) as z describes a contour consisting of small and large circles terminated by the lines arg z = ±π, together with the parts of these lines terminated by the circular arcs [cf. Fig. A.1]. Let us refer to the circles as Γ and γ, the equations of which are |z| = R and |z| = δ, respectively. It is clear then that the number of roots of K_{ν}(z) in the pair of quadrants under consideration is equal to the number of roots of z^{ν}K_{ν}(z) inside the contour, and therefore is equal to \frac{1}{2π} times the change in phase of z^{ν}K_{ν}(z) as z traverses the contour. The change in phase is now

[arg(z^{ν}K_{ν}(z))]_{Γ} − [arg(z^{ν}K_{ν}(z))]_{γ} + [arg(z^{ν}K_{ν}(z))]_{R exp(πi)}^{δ exp(πi)} + [arg(z^{ν}K_{ν}(z))]_{δ exp(−πi)}^{R exp(−πi)}. (A.7)

As R → ∞ and δ → 0, the first two terms tend respectively towards 2π(ν − \frac{1}{2}) and zero. This is evident from the consideration that the following asymptotic expansions are valid

Figure A.1: Keyhole contour with arcs Γ : |z| = R and γ : |z| = δ

when |arg z| < π [2]. Respectively, when |z| is large or small on the contour,

z^{ν}K_{ν}(z) ∼ z^{ν− \frac{1}{2}} e^{−z}\sqrt{\frac{π}{2}} ,  z^{ν}K_{ν}(z) ∼ 2^{ν−1}Γ(ν).  (A.8)

The second expansion requires modification when ν = 0.

Introduce the phase (angle) function ph(u, v) which is defined as ph(u, v) = arg(w), where w = u + iv and arg is the argument function in (A.6). (In this appendix u, v clearly have a different meaning than in the main text.) The phase function ph(u, v) = arctan(v/u) if u > 0, but otherwise gives the correct values in the second and third quadrants. Moreover, −π < ph(u, v) < π and the expression is undefined for v = 0, u ≤ 0. Now, above and below the negative real z-axis z^{ν}K_{ν}(z) can be expressed as

(re^{+πi})^{ν}K(re^{+πi}) = r^{ν}[K_{ν}(r) + π sin(νπ)I_{ν}(r)] − ir^{ν}π cos(νπ)I_{ν}(r)  (A.9)
(re^{−πi})^{ν}K(re^{−πi}) = r^{ν}[K_{ν}(r) + π sin(νπ)I_{ν}(r)] + ir^{ν}π cos(νπ)I_{ν}(r).  (A.10)

Therefore, the combination of the last two terms in (A.6) is

lim_{{δ→0^{+} \atop R→∞}} 2[ph(K_{ν}(r) + π sin(νπ)I_{ν}(r), π cos(νπ)I_{ν}(r))]_{δ}^{R}. (A.11)

Now, K_{ν}(r) is a positive decreasing function of r that vanishes at infinity, while I_{ν}(r) is a positive increasing function of r that vanishes at r = 0. Therefore, K_{ν}(r) + π sin(νπ)I_{ν}(r) has one root if sin(νπ) is negative, and no root if sin(νπ) is non-negative. In the latter case the last expression is

lim_{{δ→0^{+} \atop R→∞}} 2 [arctan \frac{π cos(νπ)I_{ν}(r)}{K_{ν}(r) + π sin(νπ)I_{ν}(r)}]_{δ}^{R} . (A.12)

Taking the δ → 0^{+}, R → ∞ limits in (A.11), we then find the phase

2 ph( sin(νπ), cos(νπ)) (A.13)

for the combination of the last two terms in (A.6). The phase ph( sin(νπ), cos(νπ)) has the same sign as cos(νπ), and, moreover, is less than π in absolute value.

Therefore, by the Argument Principle, the total number of roots of K_{ν}(z) in the pair of quadrants for which Re(z) < 0 and |argz| < π is

# of roots of K_{ν}(z) = ν − \frac{1}{2} + \frac{1}{π} ph( sin(νπ), cos(νπ))

= ν − \frac{1}{2} + \frac{1}{π} ph( cos(\frac{π}{2} − νπ), sin(\frac{π}{2} − νπ))

= ν − \frac{1}{2} + \frac{1}{π} (\frac{π}{2} − νπ + 2πk)

= 2k,

(A.14)

where k is chosen so that 2k is closest to ν − \frac{1}{2}.

References

[1] Milton Abramowitz and Irene A. Stegun Handbook of Mathematical Functions 1964: National Bureau of Standards.

[2] G. N. Watson A Treatise on the Theory of Bessel Functions, 2nd. ed 1966: Cambridge University Press.

[3] Marcus J. Grote and Joseph B. Keller “On Nonreflecting Boundary Conditions” Journal of Computational Physics, 122, 231 - 243, 1995.

[4] F. W. J. Olver, The Asymptotic Expansion of Bessel Functions of Large Order, Phil. Trans. Roy. Soc. Lond. A 247 (1954) 328-368.

[5] P. K. Chattopadhyay, Mathematical Physics, 1990: New Age International (P) Limited, Publishers.

[6] D. Kalman, A Matrix Proof of Newton’s Identities, Mathematics Magazine, vol. 73, Nov. 4, October 2000.