Orthogonality, Lommel integrals and cross product zeros of linear combinations of Bessel functions

The cylindrical Bessel differential equation and the spherical Bessel differential equation in the interval \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R \le r \le \gamma R$$\end{document}R≤r≤γR with Neumann boundary conditions are considered. The eigenfunctions are linear combinations of the Bessel function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Phi _{n,\nu }(r)=Y_{\nu }^{\prime }(\lambda _{n,\nu }) J_{\nu }(\lambda _{n,\nu } r/R)-J_{\nu }^{\prime }(\lambda _{n,\nu }) Y_{\nu }(\lambda _{n,\nu } r/R)$$\end{document}Φn,ν(r)=Yν′(λn,ν)Jν(λn,νr/R)-Jν′(λn,ν)Yν(λn,νr/R) or linear combinations of the spherical Bessel functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\psi _{m,\nu }(r)=y_{\nu }^{\prime }(\lambda _{m,\nu }) j_{\nu }(\lambda _{m,\nu } r/R)-j_{\nu }^{\prime }(\lambda _{m,\nu }) y_{\nu }(\lambda _{m,\nu } r/R)$$\end{document}ψm,ν(r)=yν′(λm,ν)jν(λm,νr/R)-jν′(λm,ν)yν(λm,νr/R). The orthogonality relations with analytical expressions for the normalization constant are given. Explicit expressions for the Lommel integrals in terms of Lommel functions are derived. The cross product zeros \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_{\nu }^{\prime }(\lambda _{n,\nu }) J_{\nu }^{\prime }(\gamma \lambda _{n,\nu })-J_{\nu }^{\prime }(\lambda _{n,\nu }) Y_{\nu }^{\prime }(\gamma \lambda _{n,\nu }) = 0$$\end{document}Yν′(λn,ν)Jν′(γλn,ν)-Jν′(λn,ν)Yν′(γλn,ν)=0 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y_{\nu }^{\prime }(\lambda _{m,\nu }) j_{\nu }^{\prime }(\gamma \lambda _{m,\nu })-j_{\nu }^{\prime }(\lambda _{m,\nu }) y_{\nu }^{\prime }(\gamma \lambda _{m,\nu }) = 0$$\end{document}yν′(λm,ν)jν′(γλm,ν)-jν′(λm,ν)yν′(γλm,ν)=0 are considered in the complex plane for real as well as complex values of the index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document}ν and approximations for the exceptional zero \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{1,\nu }$$\end{document}λ1,ν are obtained. A numerical scheme based on the discretization of the two-dimensional and three-dimensional Laplace operator with Neumann boundary conditions is presented. Explicit representations of the radial part of the Laplace operator in form of a tridiagonal matrix allow the simple computation of the cross product zeros.

Cylindrical Bessel functions and their properties are well described, as for example in the textbook of Watson (1966) where orthogonality relations and Lommel integrals are analyzed in detail. However, for a linear combination of spherical Bessel functions such orthogonality relations or Lommel integrals are not given explicitly, although they occur in many problems of diffusion or heat conduction in spherical geometries (for a review, see Carslaw and Jaeger 1959). This is mainly due to the fact that spherical Bessel functions with an integer index can be expressed in terms of trigonometric functions, and, thus, the remaining orthogonality relations and Lommel integrals can be given explicitly. Yet, the case for arbitrary complex-valued indices of cylindrical and spherical Bessel functions has not been analyzed and there are no generally valid expressions for the orthogonality relations and the Lommel integral. Moreover, the cross product zeros of Bessel functions that provide the eigenvalues in the case of an eigenfunction expansion are highly relevant for inverse eigenvalue problems and go back to the famous question "Can one hear the shape of a drum?" asked by Kac (1966). Cross product zeros were first analyzed by McMahon (1894) who derived an expression for large zeros. In addition, Cochran examined the asymptotic nature and analyticity of cross product Bessel functions (Cochran 1964(Cochran , 1966a. Recurrence relations for the Bessel function cross products are given by Goodwin (1949). But, until now, the cross products of cylindrical Bessel functions or spherical Bessel functions for an arbitrary complex index of the Bessel function have not yet been discussed in depth, although its application in physics becomes increasingly important, e.g. in optics or quantum mechanics, where non-hermitean potentials are involved.
In this work, the cylindrical and spherical Bessel differential equation, respectively, are considered in the radial interval R ≤ r ≤ γ R where γ ≥ 1. In "Preliminary facts", we investigate the orthogonality relation between eigenfunctions and the Lommel integral for a linear combination of cylindrical Bessel functions J and Y. In an analogous procedure, general expressions for the orthogonality relation and Lommel integral for a linear combination of spherical Bessel functions j and y are derived by using similarity relations between cylindrical and spherical Bessel functions. In "Cylindrical Bessel functions", the cylindrical Bessel differential equation R 2 [ν 2 /r 2 − � r ]� n,ν (r) = 2 n,ν � n,ν (r) is considered for Neumann boundary conditions on both ends of the radial interval. Explicit expressions for the orthogonality relation and Lommel integral of the according eigenfunctions � n,ν (r) = Y ′ ν ( n,ν )J ν ( n,ν r/R) − J ′ ν ( n,ν )Y ν ( n,ν r/R) are given where the eigenvalues n,ν are determined by the cross product zeros Y ′ ν ( n,ν )J ′ ν (γ n,ν ) − J ′ ν ( n,ν )Y ′ ν (γ n,ν ) = 0. To analyze cross product zeros for a complex valued index ν, the function is considered in the complex -plane and symmetry relations of the function f ν ( ) are provided. The first exceptional zero is approximated by a Taylor expansion of the function f ν ( ). A discretization scheme for the radial part of two-dimensional Laplace-operator with Neumann boundary conditions is proposed and given in terms of a tridiagonal matrix. Consequently, the cross-product zeros can be obtained by solving a simple matrix eigenvalue problem. In "Spherical Bessel functions", similar results are derived for the spherical Bessel differential equation R 2 [ν[ν + 1]/r 2 − � r ]ψ m,ν (r) = 2 m,ν ψ m,ν (r) and the corresponding eigenfunctions ψ m,ν (r) = y ′ ν ( m,ν )j ν ( m,ν r/R) − j ′ ν ( m,ν )y ν ( m,ν r/R). In "Numerical implementation", numerical algorithms are provided that correctly implement discretization schemes for the radial part of the two-dimensional and three-dimensional Laplace-operator and that allow a direct determination of cross product zeros for both cylindrical and spherical Bessel functions. Applications to the diffusion process around dipole fields are described in "Application to the diffusion process around dipole fields". A summary of the results and conclusions are given in "Summary and conclusions".

Cylindrical Bessel functions
For a linear combination of Bessel functions the following orthogonality relation holds [see Eq. (11.4.2) on page 485 in Abramowitz and Stegun (1972) or Eq. (11) on page 135 in Watson (1966) or Eqs. (19) and (20)  where the index ν denotes the order of the cylindrical Bessel function and n numerates the eigenvalues n,ν which obey the boundary conditions The normalization constant is The Lommel integral involving the linear combination of Bessel functions and a power function is given in the following form: which can be given in terms of the function � n,ν (z) and its derivative at the boundaries [see Eq. (5) on page 350 in paragraph 10.74 in Watson (1966) in combination with Eq.
(3) on page 83 in paragraph 3.9 in Watson (1966)]. It can be evaluated using the Lommel functions S κ,ν (z) = s κ,ν (z) (2) on page 285 and on page 286 in paragraph 9.3 in Watson (1966)]. The exceptional case S −1,0 (z) is treated in Watson (1966) and Glasser (2010). The most general case in which Lommel functions S κ,ν can be expressed in terms of power functions occurs when the sum of the indices κ + ν = 2p + 1 is an odd integer [see Eq. (8) in section 10.74 on page 351 in Watson (1966)]: where the Gegenbauer polynomials A n,ν are defined as [see Eq. (1) in section 9.2 on page 283 in Watson (1966)]: From this definition it is easy to obtain the relation and, further, to give the following explicit expressions: The term for A 0,ν (z) given in Eq. (25) is in agreement with the results of Watson [see Eq. (7) in section 9.2 on page 283 in Watson (1966)]. Other special cases of the Lommel functions can be found in section 3.10.2 on pages 111-113 in Magnus et al. (1966).

Spherical Bessel functions
To obtain similar relations for spherical Bessel functions, it is advantageous to consider the linear combination and to use the following general relation between cylindrical Bessel functions and spherical Bessel functions: Introducing the relations (29) and (30) into the orthogonality relation (2) of cylindrical Bessel functions, the following orthogonality relation for spherical Bessel functions can be obtained: where the same notation as in the case of cylindrical Bessel functions is used, i.e. the index µ denotes the order of the spherical Bessel functions and the index m enumerates the eigenvalues m,µ which obey the boundary condition The corresponding normalization constant is given by The Lommel integral can be evaluated by using Eqs. (29) and (30):

Cross-product zeros
The cylindrical Bessel differential equation with the radial part of two-dimensional Laplace operator � r = ∂ rr + r −1 ∂ r inside an annular ring R ≤ r ≤ γ R is considered. Reflecting boundary conditions at r = R and r = γ R are assumed: The respective eigenfunctions are given by which obviously obey the reflecting boundary condition at r = R.
Reflecting boundary conditions at the outer boundary at r = γ R lead to the conditions that have to be evaluated numerically. The function f ν ( ) possesses the following symmetry properties that translates to the properties of eigenvalues and eigenfunctions For real values of the Bessel function index ν the function f ν ( ) always takes real values and the eigenvalues n,ν are real. The respective eigenvalues n,ν are tabulated (Bauer 1964;Bridge and Angrist 1962;Truell 1943) or they can be found by using numerical routines (Sorolla et al. 2013) and empirical approximations (Laslett and Lewish 1962). The computer algebra system MATHEMATICA ® (Wolfram Research, Inc., Champaign, IL, USA, Wolfram 1999) provides the command for numerical computation of the eigenvalues. However, these standard numerical routines based on search algorithms for zeros can run into problems for inappropriate initial conditions of the algorithm. Consequently eigenvalues can be missed or found multiple. Furthermore, for complex values of the index of the Bessel functions, these numerical routines often fail. To circumvent these problems, an algorithm based on solving the corresponding original eigenvalue problem will be introduced in "Numerical implementation".
For complex values of the Bessel function index ν, however, it is helpful to consider the complex valued function f ν ( ) in the complex -plane as shown in Fig. 1. Similar considerations are performed in Figure 2 in .
For small values of the parameter , the function f ν ( ) given in Eq. (41) can be approximated by and the first zero of this approximated function is In the limit γ → 1, the first eigenvalue tends to In addition, for small values of the index ν, the first eigenvalue can be obtained from a Taylor expansion: In the same limit, Gottlieb derived [see Eq. (A.5) in ]:

Discretization
For complex values of the index ν, it is cumbersome to find the respective eigenvalues in the complex plane. To avoid a global search, a simple numerical procedure in terms of a discretization scheme is introduced. For this, the interval R ≤ r ≤ γ R is divided into p intervals with length h in the following way: As shown in Ziener et al. (2009Ziener et al. ( , 2015, the radial part of the two-dimensional Laplace operator can be written in the form Fractional indices are meant as r j 1 Additionally, the main diagonal of this discretized radial part of the Laplace operator is given by −2/h 2 except for the first and the last row. Furthermore, the sum of all elements in each row must vanish. Thus, the discretized form of the Bessel differential equation (36) can be expressed as an eigenvalue equation: where the eigenvector contains the values of the discretized eigenfunction � n,ν (r j ).

Orthogonality
Due to Abel's identity of the Wronski-determinant, the eigenfunctions take the value at the inner boundary at r = R. At the outer boundary at r = γ R, the eigenfunctions can be evaluated using the Wronski-determinant (63) and the eigenvalue equation (41): The orthogonality relation of the eigenfunctions (39) can be obtained from the general expressions (2) and (5) by replacing the eigenfunctions with the expressions from Eqs. (63) and (64):

Lommel integral
Applying eigenfunction � n,ν (r) from Eq. (39) to the general Lommel integral in Eq. (6), we arrive at where Eqs. (63) and (64) were used. The special case α = 1 and ν = 0 can be evaluated using the relation Eq. (12) with the respective Neumann polynomial given in Eq. (14). Since the corresponding Lommel function takes the constant value S 1,0 (z) = 1, this special case unfolds as i.e. the constant function 1 and the function � n,0 (r) are orthogonal in agreement with Thambynayagam (2011) (chapter 2.5, page 36). This result can be generalized by integrating the original cylindrical Bessel differential equation (36):

Cross-product zeros
The spherical Bessel differential equation with the radial part of the three-dimension Laplace operator � r = ∂ rr + 2r −1 ∂ r inside a spherical shell R ≤ r ≤ γ R is considered. Reflecting boundary conditions at r = R and r = γ R are assumed: γ R R dr r� n,0 (r) = 0 , The respective eigenfunctions are given by which already fulfill the reflecting boundary condition at r = R. The reflecting boundary conditions at the outer boundary at r = γ R lead to an equation for the determination of the eigenvalues: For real values of the Bessel function index ν the function g ν ( ) always takes real values and the eigenvalues m,ν are real. For complex values of the Bessel function index ν it is advantageous to consider the complex valued function g ν ( ) in the complex -plane as visualized in Fig. 3. The function g ν ( ) possesses similar symmetry properties as in the two-dimensional case: λ 1,ν λ 2,ν λ 3,ν λ 4,ν λ 1,ν λ 2,ν λ 3,ν λ 4,ν Fig. 3 Contour plot of the function g ν ( ) in the complex -plane. a As in Fig. 1, the real and the imaginary part of the function g ν ( ) vanish on the blue lines and red lines, respectively, and the eigenvalues can be found at the intersection points of blue and red lines. For a real index (for ν = 2 and γ = 5 as an example), all intersection points are located on the real axis. b For a complex index (exemplarily shown for ν = 2 + i and γ = 5), the geometry of the contour plots appears again more complicated as for real indices.

This translates to the symmetry properties to eigenvalues and eigenfunctions
For small values of the parameter , the function g ν ( ) given in Eq. (77) can be approximated by and, thus, the first eigenvalue can be found at In the limit γ → 1, the first eigenvalue tends to For small values of the index ν, the approximate value of first eigenvalue from Eq. (87) can be approximated further through a Taylor expansion: In the same limit for small values of the parameter γ, Grebenkov found [see Eq. (34) in Grebenkov (2007)]: and Gottlieb obtained [see Eq. (2.14) in  or Eq. (A.10) in ]: A comparison of the different approximations and the exact numerical solution is shown in Fig. 4. For m ≥ 2, we find with the approximation of McMahon (1894): that is in agreement with Eq. (36) in Grebenkov (2007).

Discretization
In analogy to the two-dimensional case, the radial interval R ≤ r ≤ γ R can be discretized in the same form as given in Eq. (59): As shown in Ziener et al. (2009), the radial part of the three-dimensional Laplace operator can be represented as   Thus, the discretized form of the spherical Bessel differential equation (72) can be written in the form of the eigenvalue equation: where the eigenvector contains the values of the discretized eigenfunction ψ m,ν (r j ).

Orthogonality
Due to Abel's identity of the Wronski-determinant, the eigenfunctions take the value at the inner boundary at r = R. The value of the eigenfunction at the outer boundary at r = γ R can be evaluated using the Wronski-determinant (97) and the eigenvalue equation (77): The orthogonality relation of the eigenfunctions (75) can be obtained from the general expressions (31) and (34) by replacing the respective eigenfunctions with the expressions in Eqs. (97) and (99) as in the case of cylindrical Bessel functions:

Lommel integral
The Lommel integral can be obtained from the general expression given in Eq. (35), respecting the Neumann boundary conditions, the Wronski-determinant from Eq. (97) and expression (98): For α + ν = 2p + 1, the Lommel functions can be expressed in terms of Gegenbauer polynomials as above [see Eq. (22)].
In analogy to Eq. (70), the special case for α = 2 and ν = 0 can be determined, since in this case the prefactors of the Lommel functions in Eq. (103) vanish: Therefore, the functions 1 and ψ m,0 (r) are orthogonal. In analogy to the cylindrical case in Eq. (71) this result can be generalized by integrating the original spherical Bessel differential equation (72) Numerical implementation

Cylindrical Bessel functions
The discretization scheme provided for the radial part of the two-dimension Laplace operator with Neumann boundary conditions is implemented in the following algorithm: In lines 1 and 2, the length of the discretization interval h and the discretized radius according to the scheme (59) is computed. In lines 3-9, the discretized form of radial part of the two-dimensional Laplace operator given in Eq. (60) is implemented and, in line 10, the diagonal matrix ν 2 diag(1/r 2 1 , . . . , 1/r 2 p ) represents the part ν 2 /r 2 of the original cylindrical Bessel differential equation. Finally, in line 11, the eigenvalue equation (61) is solved and the eigenvalues are given in a list sorted in ascending order. The eigenvalue n,ν can be obtained by the command

Spherical Bessel functions
In analogy to the implementation for cylindrical Bessel functions, the eigenvalues of the radial part of the three-dimensional discretization scheme according to Eq. (94) and (95) can be implemented as follows: Finally, the eigenvalue m,ν of the matrix eigenvalue problem (95) can be obtained by the command

Application to the diffusion process around dipole fields
To visualize the applicability of the results obtained in the previous sections, the frequency autocorrelation function of diffusing spins with the local resonance frequency is considered as an example. The function f (r) describes the shape of the local resonance frequency and δω its strength. The same analysis has already been performed in Ziener et al. (2008) for the diffusion between two concentric cylinders (two-dimensional case) and between two concentric spheres (three-dimensional case). The frequency autocorrelation function can be used to analyze magnetic resonance pulse sequences and to determine properties of red blood cells as shown in Ziener et al. (2010).
However, using the explicit expressions for the normalization and the Lommel integral for linear combinations of cylindrical or spherical Bessel functions, respectively, it is possible to obtain considerably simpler expressions as will be demonstrated below.
The frequency autocorrelation function for spins diffusing in the volume V is defined as (Ziener et al. 2006) (106) ω(r) = δωf (r) where p(r, r 0 , t) is the probability of a spin to diffuse during the time span t from the position r 0 to the position r and p(r 0 ) = 1/V is the probability of finding a spin at position r 0 . This transition probability is a solution of the diffusion equation where D is the diffusion coefficient inside the volume V in which the diffusion occurs (Bauer et al. 2005). The transition probability p(r, r 0 , t) can be expressed in terms of an eigenfunction expansion with eigenvalues κ l and eigenfunctions φ l that fulfil the orthogonality condition and are solutions of the eigenvalue equation with R being the radius of the inner cylinder or the inner sphere, respectively. The eigenfunction corresponds to the lowest eigenvalue

Cylinders
In the two-dimensional case, the diffusion occurs in the space between two concentric circles with radius R and γ R: The eigenfunctions of the eigenvalue equation (112) are This expression for the expansion coefficients is evidently a simpler expression than its equivalent given in Eq. (B1) in Ziener et al. (2008). The corresponding eigenvalues n,2 can be obtained from Eq. (41) for ν = 2:

Fig. 5
Diffusion propagator according to Eq. (119) around a cylinder with radius R = 1 µm, diffusion coefficient D = 1 µm 2 /ms and γ = 8. The initial position is at r 0 = 4 µm and φ 0 = π/4. For small times (exemplary for t = 1 ms) the transition probability is concentrated around the initial position. For large times (exemplary for t = 100 ms) the diffusion propagator takes the constant value lim t→∞ p(r, φ, r 0 , φ 0 , t) = and the following sums can be used to check if a sufficient numerical accuracy is obtained:

Spheres
In the three-dimensional case the diffusion occurs between two concentric spheres, and, thus, the diffusion volume corresponds to The orthogonal eigenfunctions of the eigenvalue equation (112)   (131) V = 4 3 πR 3 γ 3 − 1 .
Introducing this shape function and the three-dimensional propagator from Eq. (138) into the definition of the autocorrelation function (108), leads to with the expansion coefficients The remaining Lommel integral can be solved by using the expression (103) for α = −1 and ν = 2. The occurring Lommel functions S − 5 2 ,+ 7 2 (z) and S − 3 2 ,+ 5 2 (z) can explicitly be given by using the relation (22): where A 0, 7 2 (z) and A 0, 5 2 (z) can be obtained from Eq. (25). Thus, the integration yields: Introducing the normalization constant from Eq. (102) yields: For ν = 2 the spherical Bessel functions can be expressed in terms of the sine and cosine function:  j ′ 2 ( m,2 ) j ′ 2 (γ m,2 ) 2 γ 2 2 m,2 −6 γ 3 j ′ 2 ( m,2 ) j ′ 2 (γ m,2 ) 2 + 6 − 2 m,2 which will be necessary for further simplifications. The derivative of the spherical Bessel functions can be simplified by using the relation (151): Finally, the expansion coefficients can be written as or in the form This expression for the expansion coefficients is much more simpler than the equivalent expression given in Eq. (B4) in Ziener et al. (2008) or in Eq. (7) in Ziener et al. (2010).