A physicist's guide to explicit summation formulas involving zeros of Bessel functions and related spectral sums

TL;DR

This paper reviews methods for explicitly computing spectral sums involving Bessel function zeros, crucial for various diffusion-related applications, emphasizing practical strategies and their limitations.

Contribution

It provides a comprehensive pedagogic overview of analytical tools and summation formulas for spectral sums involving Bessel function zeros, with practical insights.

Findings

Summarizes three main strategies for spectral sum computation.

Highlights advantages and limitations of each method.

Provides a collection of summation formulas from literature.

Abstract

We summarize the mathematical basis and practical hints for the explicit analytical computation of spectral sums that involve the eigenvalues of the Laplace operator in simple domains. Such spectral sums appear as spectral expansions of heat kernels, survival probabilities, first-passage time densities, and reaction rates in many diffusion-oriented applications. As the eigenvalues are determined by zeros of an appropriate linear combination of a Bessel function and its derivative, there are powerful analytical tools for computing such spectral sums. We discuss three main strategies: representations of meromorphic functions as sums of partial fractions, Fourier-Bessel and Dini series, and direct evaluation of the Laplace-transformed heat kernels. The major emphasis is put on a pedagogic introduction, the practical aspects of these strategies, their advantages and limitations. The review summarizes many summation formulas for spectral sums that are dispersed in the literature.

Figures

Tables

Table 1. Table 1: Laplacian eigenfunctions and eigenvalues in a disk and a ball of radius R 𝑅 R . We recall that J n ​ ( z ) subscript 𝐽 𝑛 𝑧 J_{n}(z) and j n ​ ( z ) subscript 𝑗 𝑛 𝑧 j_{n}(z) are the Bessel and spherical Bessel functions of the first kind, P n l ​ ( z ) superscript subscript 𝑃 𝑛 𝑙 𝑧 P_{n}^{l}(z) are the associated Legendre polynomials, and h = R / ℓ ℎ 𝑅 ℓ h=R/\ell . Here α n ​ k subscript 𝛼 𝑛 𝑘 \alpha_{nk} are all strictly positive zeros of the function g n ​ ( z ) subscript 𝑔 𝑛 𝑧 g_{n}(z) , enumerated by k = 1 , 2 , 3 , … 𝑘 1 2 3 … k=1,2,3,\ldots , except for the Neumann boundary condition ( h = 0 ℎ 0 h=0 ), for which α 01 = 0 subscript 𝛼 01 0 \alpha_{01}=0 and β 01 = 1 subscript 𝛽 01 1 \beta_{01}=1 for the disk, as well as α 01 = 0 subscript 𝛼 01 0 \alpha_{01}=0 and β 01 = 3 / 2 subscript 𝛽 01 3 2 \beta_{01}=\sqrt{3/2} for the ball. Note that the eigenfunctions for the disk are often written in the real-valued form, in which e i ​ n ​ ϕ superscript 𝑒 𝑖 𝑛 italic-ϕ e^{in\phi} is replaced by cos ⁡ ( n ​ ϕ ) 𝑛 italic-ϕ \cos(n\phi) and sin ⁡ ( n ​ ϕ ) 𝑛 italic-ϕ \sin(n\phi) with n = 0 , 1 , 2 , … 𝑛 0 1 2 … n=0,1,2,\ldots .

	Disk	Ball
Domain $\Omega$
Coordinates	Polar $(r,\varphi )$	Spherical $(r,\theta ,\varphi )$
Laplacian $\Delta$	${\partial }_r^2+\frac{1}{r}{\partial }_r$	${\partial }_r^2+\frac{2}{r}{\partial }_r$
	$+\frac{1}{r^2}{\partial }_{\varphi }^2$	$+\frac{1}{r^2}\left(\frac{1}{\sin \theta }{\partial }_{\theta }\sin \theta {\partial }_{\theta }+\frac{1}{{\sin }^2\theta }{\partial }_{\varphi }^2\right)$
eigenfunctions $U$	$C_{nk}{J}_n\left({\alpha }_{nk}r/R\right)e^{in\varphi }$	$C_{nk}{j}_n\left({\alpha }_{nk}r/R\right)P_n^l\left(\cos \theta \right)e^{il\varphi }$
indices	$n\in \mathbb{Z}$	$n=0,1,2,\dots$,   $l=-n,\dots ,n$
function $g_n\left(z\right)$	$z{J}_n^{\prime }\left(z\right)+h{J}_n\left(z\right)$	$z{j}_n^{\prime }\left(z\right)+h{j}_n\left(z\right)$
eigenvalues	${\lambda }_{nk}={\alpha }_{nk}^2/R^2$	${\lambda }_{nk}={\alpha }_{nk}^2/R^2$
$C_{nk}$ for $\ell >0$	$\frac{{\beta }_{nk}}{\sqrt{\pi }{J}_n\left({\alpha }_{nk}\right)R}$	$\frac{\sqrt{2n+1}{\beta }_{nk}}{\sqrt{2\pi }{j}_n\left({\alpha }_{nk}\right)R^{3/2}}$
$C_{nk}$ for $\ell =0$	$\frac{1}{\sqrt{\pi }{J}_n^{\prime }\left({\alpha }_{nk}\right)R}$	$\frac{\sqrt{2n+1}}{\sqrt{2\pi }{j}_n^{\prime }\left({\alpha }_{nk}\right)R^{3/2}}$
coefficient ${\beta }_{nk}$	$\frac{{\alpha }_{nk}}{\sqrt{{\alpha }_{nk}^2-n^2+h^2}}$	$\frac{{\alpha }_{nk}}{\sqrt{{\alpha }_{nk}^2-n\left(n+1\right)+h\left(h-1\right)}}$

Table 2. Table 3: Summary of summation formulas for a disk and a ball of radius b 𝑏 b . Here 0 ≤ x ≤ x 0 ≤ 1 0 𝑥 subscript 𝑥 0 1 0\leq x\leq x_{0}\leq 1 , and α n ​ k = b ​ q n ​ k subscript 𝛼 𝑛 𝑘 𝑏 subscript 𝑞 𝑛 𝑘 \alpha_{nk}=bq_{nk} are the positive zeros of g n ​ ( z ) = z ​ J n ′ ​ ( z ) + h ​ J n ​ ( z ) subscript 𝑔 𝑛 𝑧 𝑧 subscript superscript 𝐽 ′ 𝑛 𝑧 ℎ subscript 𝐽 𝑛 𝑧 g_{n}(z)=zJ^{\prime}_{n}(z)+hJ_{n}(z) in two dimensions, and of g n ​ ( z ) = z ​ j n ′ ​ ( z ) + h ​ j n ​ ( z ) subscript 𝑔 𝑛 𝑧 𝑧 subscript superscript 𝑗 ′ 𝑛 𝑧 ℎ subscript 𝑗 𝑛 𝑧 g_{n}(z)=zj^{\prime}_{n}(z)+hj_{n}(z) in three dimensions, with h = b / ℓ b ℎ 𝑏 subscript ℓ 𝑏 h=b/\ell_{b} . Expressions for Neumann boundary condition are simply deduced by setting h = 0 ℎ 0 h=0 . Several misprints were found in the published version of this Table; they are corrected here and highligted by red color.

	Robin ($h_1=h$)	Neumann ($h_1=0$)	Dirichlet ($h_1=\infty$)
Robin $\left(h_0=h\right)$	$\tan {\alpha }_k=\frac{2h{\alpha }_k}{{\alpha }_k^2-h^2}$	${\alpha }_k\tan {\alpha }_k=h$	$\tan {\alpha }_k=-{\alpha }_k/h$
	$\sum _{k=1}^{\infty }\frac{u_k\left(x\right)u_k\left(x_0\right)}{\left(z^2+{\alpha }_k^2\right)\left({\alpha }_k^2+2h+h^2\right)}$	$\sum _{k=1}^{\infty }\frac{u_k\left(x\right)u_k\left(x_0\right)}{\left(z^2+{\alpha }_k^2\right)\left({\alpha }_k^2+h+h^2\right)}$	$\sum _{k=1}^{\infty }\frac{u_k\left(x\right)u_k\left(x_0\right)}{\left(z^2+{\alpha }_k^2\right)\left({\alpha }_k^2+h+h^2\right)}$
	$=\frac{v(1-x_0;z,h)v(x;z,h)}{2z\left[\left(z^2+h^2\right)\sinh z+2hz\cosh z\right]}$	$=\frac{v(x;z,h)\cosh \left(z\left(1-x_0\right)\right)}{2z\left[z\sinh z+h\cosh z\right]}$	$=\frac{\sinh \left(z\left(1-x_0\right)\right)v(x;z,h)}{2z\left[h\sinh z+z\cosh z\right]}$
	$Z=\frac{z\left(z^2+h^2\right)\cosh z+\left(z^2\left(1+2h\right)-h^2\right)\sinh z}{2z^2\left(\left(z^2+h^2\right)\sinh z+2zh\cosh z\right)}$	$Z=\frac{z\cosh z+\left(h+1\right)\sinh z}{2z\left(z\sinh z+h\cosh z\right)}$	$Z=\frac{hz\cosh z+\left(z^2-h\right)\sinh z}{2z^2\left(h\sinh z+z\cosh z\right)}$
Neumann $\left(h_0=0\right)$	${\alpha }_k\tan {\alpha }_k=h$	${\alpha }_k=\pi \left(k-1\right)$	${\alpha }_k=\pi \left(k-1/2\right)$
	$\sum _{k=1}^{\infty }\frac{\cos \left({\alpha }_{k}x\right)\cos \left({\alpha }_k{x}_0\right)\left({\alpha }_k^2+h^2\right)}{\left(z^2+{\alpha }_k^2\right)\left({\alpha }_k^2+h+h^2\right)}$	$\frac{1}{2z^2}+\sum _{k=2}^{\infty }\frac{\cos \left({\alpha }_{k}x\right)\cos \left({\alpha }_k{x}_0\right)}{z^2+{\alpha }_k^2}$	$\sum _{k=1}^{\infty }\frac{\cos \left({\alpha }_{k}x\right)\cos \left({\alpha }_k{x}_0\right)}{z^2+{\alpha }_k^2}$
	$=\frac{v(1-x_0;z;h)\cosh \left(zx\right)}{2z\left[z\sinh z+h\cosh z\right]}$	$=\frac{\cosh \left(z\left(1-x_0\right)\right)\cosh \left(zx\right)}{2z\sinh z}$	$=\frac{\sinh \left(z\left(1-x_0\right)\right)\cosh \left(zx\right)}{2z\cosh z}$
	$Z=\frac{z\cosh z+\left(h+1\right)\sinh z}{2z\left(z\sinh z+h\cosh z\right)}$	$Z=\frac{z\cosh z+\sinh z}{2z^2\sinh z}$	$Z=\frac{\sinh z}{2z\cosh z}$
Dirichlet $\left(h_0=\infty \right)$	$\tan {\alpha }_k=-{\alpha }_k/h$	${\alpha }_k=\pi \left(k-1/2\right)$	${\alpha }_k=\pi k$
	$\sum _{k=1}^{\infty }\frac{\sin \left({\alpha }_{k}x\right)\sin \left({\alpha }_k{x}_0\right)\left({\alpha }_k^2+h^2\right)}{\left(z^2+{\alpha }_k^2\right)\left({\alpha }_k^2+h+h^2\right)}$	$\sum _{k=1}^{\infty }\frac{\sin \left({\alpha }_{k}x\right)\sin \left({\alpha }_k{x}_0\right)}{z^2+{\alpha }_k^2}$	$\sum _{k=1}^{\infty }\frac{\sin \left({\alpha }_{k}x\right)\sin \left({\alpha }_k{x}_0\right)}{z^2+{\alpha }_k^2}$
	$=\frac{v(1-x_0;z,h)\sinh \left(zx\right)}{2z\left[h\sinh z+z\cosh z\right]}$	$=\frac{\cosh \left(z\left(1-x_0\right)\right)\sinh \left(zx\right)}{2z\cosh z}$	$=\frac{\sinh \left(z\left(1-x_0\right)\right)\sinh \left(zx\right)}{2z\sinh z}$
	$Z=\frac{hz\cosh z+\left(z^2-h\right)\sinh z}{2z^2\left(h\sinh z+z\cosh z\right)}$	$Z=\frac{\sinh z}{2z\cosh z}$	$Z=\frac{\cosh z}{2z\sinh z}-\frac{1}{2z^2}$

Table 3. Table 3: Summary of summation formulas for a disk and a ball of radius b 𝑏 b . Here 0 ≤ x ≤ x 0 ≤ 1 0 𝑥 subscript 𝑥 0 1 0\leq x\leq x_{0}\leq 1 , and α n ​ k = b ​ q n ​ k subscript 𝛼 𝑛 𝑘 𝑏 subscript 𝑞 𝑛 𝑘 \alpha_{nk}=bq_{nk} are the positive zeros of g n ​ ( z ) = z ​ J n ′ ​ ( z ) + h ​ J n ​ ( z ) subscript 𝑔 𝑛 𝑧 𝑧 subscript superscript 𝐽 ′ 𝑛 𝑧 ℎ subscript 𝐽 𝑛 𝑧 g_{n}(z)=zJ^{\prime}_{n}(z)+hJ_{n}(z) in two dimensions, and of g n ​ ( z ) = z ​ j n ′ ​ ( z ) + h ​ j n ​ ( z ) subscript 𝑔 𝑛 𝑧 𝑧 subscript superscript 𝑗 ′ 𝑛 𝑧 ℎ subscript 𝑗 𝑛 𝑧 g_{n}(z)=zj^{\prime}_{n}(z)+hj_{n}(z) in three dimensions, with h = b / ℓ b ℎ 𝑏 subscript ℓ 𝑏 h=b/\ell_{b} . Expressions for Neumann boundary condition are simply deduced by setting h = 0 ℎ 0 h=0 . Several misprints were found in the published version of this Table; they are corrected here and highligted by red color.

	Disk		Ball
Robin ()	$\sum _{k=1}^{\infty }\frac{2{\alpha }_{nk}^2{J}_n\left({\alpha }_{nk}x\right)J_n\left({\alpha }_{nk}{x}_0\right)}{\left({\alpha }_{nk}^2-n^2+h^2\right)\left(z^2+{\alpha }_{nk}^2\right)J_n^2\left({\alpha }_{nk}\right)}=$	(D1)	$\sum _{k=1}^{\infty }\frac{2{\alpha }_{nk}^2{j}_n\left({\alpha }_{nk}x\right)j_n\left({\alpha }_{nk}{x}_0\right)}{\left({\alpha }_{nk}^2-n\left(n+1\right)+h^2-h\right)\left(z^2+{\alpha }_{nk}^2\right)j_n^2\left({\alpha }_{nk}\right)}=$	(S1)
	$\left(K_n\left(z{x}_0\right)-I_n\left(z{x}_0\right)\frac{z{K}_n^{\prime }\left(z\right)+h{K}_n\left(z\right)}{z{I}_n^{\prime }\left(z\right)+h{I}_n\left(z\right)}\right)I_n\left(zx\right)$		$z\left(k_n\left(z{x}_0\right)-i_n\left(z{x}_0\right)\frac{z{k}_n^{\prime }\left(z\right)+h{k}_n\left(z\right)}{z{i}_n^{\prime }\left(z\right)+h{i}_n\left(z\right)}\right)i_n\left(zx\right)$
	$\sum _{k=1}^{\infty }\frac{1}{z^2+{\alpha }_{nk}^2}=\frac{\left(z^2+n^2\right)I_n\left(z\right)+zh{I}_n^{\prime }\left(z\right)}{2z^2\left(z{I}_n^{\prime }\left(z\right)+h{I}_n\left(z\right)\right)}-\frac{n}{2z^2}$	(D2)	$\sum _{k=1}^{\infty }\frac{1}{z^2+{\alpha }_{nk}^2}=\frac{\left(z^2+n\left(n+1\right)\right)i_n\left(z\right)+z\left(h-1\right)i_n^{\prime }\left(z\right)}{2z^2\left(z{i}_n^{\prime }\left(z\right)+h{i}_n\left(z\right)\right)}-\frac{n}{2z^2}$	(S2)
	$\sum _{k=1}^{\infty }\frac{2{\alpha }_{nk}^2{J}_n\left({\alpha }_{nk}x\right)}{\left({\alpha }_{nk}^2-n^2+h^2\right)\left(z^2+{\alpha }_{nk}^2\right)J_n\left({\alpha }_{nk}\right)}=\frac{I_n\left(zx\right)}{z{I}_n^{\prime }\left(z\right)+h{I}_n\left(z\right)}$	(D3)	$\sum _{k=1}^{\infty }\frac{2{\alpha }_{nk}^2{j}_n\left({\alpha }_{nk}x\right)}{\left({\alpha }_{nk}^2-n\left(n+1\right)+h^2-h\right)\left(z^2+{\alpha }_{nk}^2\right)j_n\left({\alpha }_{nk}\right)}=\frac{i_n\left(zx\right)}{z{i}_n^{\prime }\left(z\right)+h{i}_n\left(z\right)}$	(S3)
	$\sum _{k=1}^{\infty }\frac{2{\alpha }_{nk}^2{J}_n\left({\alpha }_{nk}x\right)J_n\left({\alpha }_{nk}{x}_0\right)}{\left({\alpha }_{nk}^2-n^2+h^2\right)\left(z^2-{\alpha }_{nk}^2\right)J_n^2\left({\alpha }_{nk}\right)}=$	(D4)	$\sum _{k=1}^{\infty }\frac{2{\alpha }_{nk}^2{j}_n\left({\alpha }_{nk}x\right)j_n\left({\alpha }_{nk}{x}_0\right)}{\left({\alpha }_{nk}^2-n\left(n+1\right)+h^2-h\right)\left(z^2-{\alpha }_{nk}^2\right)j_n^2\left({\alpha }_{nk}\right)}=$	(S4)
	$\frac{\pi }{2}\left(Y_n\left(z{x}_0\right)-J_n\left(z{x}_0\right)\frac{z{Y}_n^{\prime }\left(z\right)+h{Y}_n\left(z\right)}{z{J}_n^{\prime }\left(z\right)+h{J}_n\left(z\right)}\right)J_n\left(zx\right)$		$z\left(y_n\left(z{x}_0\right)-j_n\left(z{x}_0\right)\frac{z{y}_n^{\prime }\left(z\right)+h{y}_n\left(z\right)}{z{j}_n^{\prime }\left(z\right)+h{j}_n\left(z\right)}\right)j_n\left(zx\right)$
	$\sum _{k=1}^{\infty }\frac{1}{z^2-{\alpha }_{nk}^2}=\frac{\left(n^2-z^2\right)J_n\left(z\right)+zh{J}_n^{\prime }\left(z\right)}{2z^2\left(z{J}_n^{\prime }\left(z\right)+h{J}_n\left(z\right)\right)}-\frac{n}{2z^2}$	(D5)	$\sum _{k=1}^{\infty }\frac{1}{z^2-{\alpha }_{nk}^2}=\frac{\left(n\left(n+1\right)-z^2\right)j_n\left(z\right)+z\left(h-1\right)j_n^{\prime }\left(z\right)}{2z^2\left(z{j}_n^{\prime }\left(z\right)+h{j}_n\left(z\right)\right)}-\frac{n}{2z^2}$	(S5)
	$\sum _{k=1}^{\infty }\frac{2{\alpha }_{nk}^2{J}_n\left({\alpha }_{nk}x\right)}{\left({\alpha }_{nk}^2-n^2+h^2\right)\left(z^2-{\alpha }_{nk}^2\right)J_n\left({\alpha }_{nk}\right)}=\frac{-J_n\left(zx\right)}{z{J}_n^{\prime }\left(z\right)+h{J}_n\left(z\right)}$	(D6)	$\sum _{k=1}^{\infty }\frac{2{\alpha }_{nk}^2{j}_n\left({\alpha }_{nk}x\right)}{\left({\alpha }_{nk}^2-n\left(n+1\right)+h^2-h\right)\left(z^2-{\alpha }_{nk}^2\right)j_n\left({\alpha }_{nk}\right)}=\frac{-j_n\left(zx\right)}{z{j}_n^{\prime }\left(z\right)+h{j}_n\left(z\right)}$	(S6)
Dirichlet ($h=\infty$)	$\sum _{k=1}^{\infty }\frac{2J_n\left({\alpha }_{nk}x\right)J_n\left({\alpha }_{nk}{x}_0\right)}{\left(z^2+{\alpha }_{nk}^2\right){\left[J_n^{\prime }\left({\alpha }_{nk}\right)\right]}^2}=$	(D7)	$\sum _{k=1}^{\infty }\frac{2j_n\left({\alpha }_{nk}x\right)j_n\left({\alpha }_{nk}{x}_0\right)}{\left(z^2+{\alpha }_{nk}^2\right){\left[j_n^{\prime }\left({\alpha }_{nk}\right)\right]}^2}=$	(S7)
	$\left(K_n\left(z{x}_0\right)-I_n\left(z{x}_0\right)\frac{K_n\left(z\right)}{I_n\left(z\right)}\right)I_n\left(zx\right)$		$z\left(k_n\left(z{x}_0\right)-i_n\left(z{x}_0\right)\frac{k_n\left(z\right)}{i_n\left(z\right)}\right)i_n\left(zx\right)$
	$\sum _{k=1}^{\infty }\frac{1}{z^2+{\alpha }_{nk}^2}=\frac{I_n^{\prime }\left(z\right)}{2z{I}_n\left(z\right)}-\frac{n}{2z^2}$	(D8)	$\sum _{k=1}^{\infty }\frac{1}{z^2+{\alpha }_{nk}^2}=\frac{i_n^{\prime }\left(z\right)}{2z{i}_n\left(z\right)}-\frac{n}{2z^2}$	(S8)
	$\sum _{k=1}^{\infty }\frac{{\alpha }_{nk}{J}_n\left({\alpha }_{nk}x\right)}{\left(z^2+{\alpha }_{nk}^2\right)J_n^{\prime }\left({\alpha }_{nk}\right)}=-\frac{I_n\left(zx\right)}{2I_n\left(z\right)}$	(D9)	$\sum _{k=1}^{\infty }\frac{{\alpha }_{nk}{j}_n\left({\alpha }_{nk}x\right)}{\left(z^2+{\alpha }_{nk}^2\right)j_n^{\prime }\left({\alpha }_{nk}\right)}=-\frac{i_n\left(zx\right)}{2i_n\left(z\right)}$	(S9)
	$\sum _{k=1}^{\infty }\frac{2J_n\left({\alpha }_{nk}x\right)J_n\left({\alpha }_{nk}{x}_0\right)}{\left(z^2-{\alpha }_{nk}^2\right){\left[J_n^{\prime }\left({\alpha }_{nk}\right)\right]}^2}=$	(D10)	$\sum _{k=1}^{\infty }\frac{2j_n\left({\alpha }_{nk}x\right)j_n\left({\alpha }_{nk}{x}_0\right)}{\left(z^2-{\alpha }_{nk}^2\right){\left[j_n^{\prime }\left({\alpha }_{nk}\right)\right]}^2}=$	(S10)
	$\frac{\pi }{2}\left(Y_n\left(z{x}_0\right)-J_n\left(z{x}_0\right)\frac{Y_n\left(z\right)}{J_n\left(z\right)}\right)J_n\left(zx\right)$		$z\left(y_n\left(z{x}_0\right)-j_n\left(z{x}_0\right)\frac{y_n\left(z\right)}{j_n\left(z\right)}\right)j_n\left(zx\right)$
	$\sum _{k=1}^{\infty }\frac{1}{z^2-{\alpha }_{nk}^2}=\frac{J_n^{\prime }\left(z\right)}{2z{J}_n\left(z\right)}-\frac{n}{2z^2}$	(D11)	$\sum _{k=1}^{\infty }\frac{1}{z^2-{\alpha }_{nk}^2}=\frac{j_n^{\prime }\left(z\right)}{2z{j}_n\left(z\right)}-\frac{n}{2z^2}$	(S11)
	$\sum _{k=1}^{\infty }\frac{{\alpha }_{nk}{J}_n\left({\alpha }_{nk}x\right)}{\left(z^2-{\alpha }_{nk}^2\right)J_n^{\prime }\left({\alpha }_{nk}\right)}=+\frac{J_n\left(zx\right)}{2J_n\left(z\right)}$	(D12)	$\sum _{k=1}^{\infty }\frac{{\alpha }_{nk}{j}_n\left({\alpha }_{nk}x\right)}{\left(z^2-{\alpha }_{nk}^2\right)j_n^{\prime }\left({\alpha }_{nk}\right)}=+\frac{j_n\left(zx\right)}{2j_n\left(z\right)}$	(S12)

Equations

$$\Delta U_m+{\lambda }_m{U}_m$$

$$\ell \frac{\partial }{\partial n}{U}_m+U_m$$

$$P(\mathbf{x},t|{\mathbf{x}}_0)=\sum _{m=1}^{\infty }{U}_m(\mathbf{x})U_m^{\ast }({\mathbf{x}}_0)e^{-Dt{\lambda }_m}$$

$$\frac{\partial }{\partial t}P(\mathbf{x},t|{\mathbf{x}}_0)=D\Delta P(\mathbf{x},t|{\mathbf{x}}_0)\text{in}\Omega ,$$

$$\tilde{P}(\mathbf{x},p|{\mathbf{x}}_0)=\underset{0}{\overset{\infty }{\int }}dt{e}^{-pt}P(\mathbf{x},t|{\mathbf{x}}_0)=\sum _{m=1}^{\infty }\frac{U_m(\mathbf{x})U_m^{\ast }({\mathbf{x}}_0)}{p+D{\lambda }_m},$$

$$\bigl{(}p-D\Delta\bigr{)}\tilde{P}(\bm{x},p|\bm{x}_{0})=\delta(\bm{x}-\bm{x}_{0})\quad\textrm{in}~{}\Omega,$$

$$Z(t)=\underset{\Omega }{\int }d{\mathbf{x}}_{0}P({\mathbf{x}}_0,t|{\mathbf{x}}_0)=\sum _{m=1}^{\infty }{e}^{-D{\lambda }_{m}t},$$

$$\tilde{Z}(p)=\underset{0}{\overset{\infty }{\int }}dt{e}^{-pt}Z(t)=\sum _{m=1}^{\infty }\frac{1}{p+D{\lambda }_m},$$

$$S(t|{\mathbf{x}}_0)={\mathbb{P}}_{{\mathbf{x}}_0}\{\tau >t\}=\underset{\Omega }{\int }d\mathbf{x}P(\mathbf{x},t|{\mathbf{x}}_0)=\sum _{m=1}^{\infty }{e}^{-Dt{\lambda }_m}{U}_m^{\ast }({\mathbf{x}}_0)\underset{\Omega }{\int }d\mathbf{x}{U}_m(\mathbf{x}),$$

$$\overline{S(t)}=\frac{1}{|\Omega |}\underset{\Omega }{\int }d{\mathbf{x}}_{0}S(t|{\mathbf{x}}_0)=\frac{1}{|\Omega |}\sum _{m=1}^{\infty }{e}^{-Dt{\lambda }_m}{\left|\underset{\Omega }{\int }d\mathbf{x}{U}_m(\mathbf{x})\right|}^2,$$

$$\overline{H(t)}=\frac{1}{|\Omega |}\underset{\Omega }{\int }d{\mathbf{x}}_{0}H(t|{\mathbf{x}}_0)=\frac{D}{|\Omega |}\sum _{m=1}^{\infty }{e}^{-Dt{\lambda }_m}{\lambda }_m{\left|\underset{\Omega }{\int }d\mathbf{x}{U}_m(\mathbf{x})\right|}^2.$$

$$\tilde{S}(p|{\mathbf{x}}_0)=\sum _{m=1}^{\infty }\frac{1}{p+D{\lambda }_m}{U}_m^{\ast }({\mathbf{x}}_0)\underset{\Omega }{\int }d\mathbf{x}{U}_m(\mathbf{x}).$$

$${\mathbb{E}}_{{\mathbf{x}}_0}\{\tau \}=\tilde{S}(0|{\mathbf{x}}_0)=\frac{1}{D}\sum _{m=1}^{\infty }{\lambda }_m^{-1}{U}_m^{\ast }({\mathbf{x}}_0)\underset{\Omega }{\int }d\mathbf{x}{U}_m(\mathbf{x}),$$

$$\overline{{\mathbb{E}}\{\tau\}}=\frac{1}{D\,|\Omega|}\sum\limits_{m=1}^{\infty}\lambda_{m}^{-1}\,\biggl{|}\int\limits_{\Omega}d\bm{x}\,U_{m}(\bm{x})\biggr{|}^{2}.$$

$$D(t)$$

$$C(t_1,\dots ,t_j|{\mathbf{x}}_0)={\mathbb{E}}_{{\mathbf{x}}_0}\{X_{t_1}{X}_{t_2}\dots X_{t_j}\}=\underset{\Omega }{\int }d{\mathbf{x}}_1\dots \underset{\Omega }{\int }d{\mathbf{x}}_j{\mathbf{x}}_j$$

$$\times P({\mathbf{x}}_j,t_j-t_{j-1}|{\mathbf{x}}_{j-1}){\mathbf{x}}_{j-1}P({\mathbf{x}}_{j-1},t_{j-1}-t_{j-2}|{\mathbf{x}}_{j-2})\dots {\mathbf{x}}_{1}P({\mathbf{x}}_1,t_1|{\mathbf{x}}_0)$$

$$C(t_1,\dots ,t_j|{\mathbf{x}}_0)=\sum _{m_1=1}^{\infty }\dots \sum _{m_j=1}^{\infty }{e}^{-D{\lambda }_{m_1}{t}_1-D{\lambda }_{m_2}(t_2-t_1)-\dots -D{\lambda }_{m_j}(t_j-t_{j-1})}$$

$$\times U_{m_1}^{\ast }({\mathbf{x}}_0)B_{m_1,m_2}{B}_{m_2,m_3}\dots B_{m_{j-1},m_j}\underset{\Omega }{\int }d{\mathbf{x}}_j{\mathbf{x}}_j{U}_{m_j}({\mathbf{x}}_j),$$

$$B_{m,m^{\prime }}=\underset{\Omega }{\int }d\mathbf{x}{U}_m(\mathbf{x})\mathbf{x}{U}_{m^{\prime }}^{\ast }(\mathbf{x}).$$

$$\tilde{C}(p_1,\dots ,p_j|{\mathbf{x}}_0)=\sum _{m_1=1}^{\infty }\dots \sum _{m_j=1}^{\infty }\frac{1}{p_1+D({\lambda }_{m_1}-{\lambda }_{m_2})}\dots \frac{1}{p_j+D{\lambda }_{m_j}}$$

$$\times U_{m_1}^{\ast }({\mathbf{x}}_0)B_{m_1,m_2}{B}_{m_2,m_3}\dots B_{m_{j-1},m_j}\underset{\Omega }{\int }d{\mathbf{x}}_j{\mathbf{x}}_j{U}_{m_j}({\mathbf{x}}_j).$$

$$\begin{array}[]{r l l}U_{k}(r)&=A_{k}\,\cos(\sqrt{\lambda_{k}}r)+B_{k}\,\sin(\sqrt{\lambda_{k}}r)&(d=1),\\ U_{nk}(r,\phi)&=\bigl{(}A_{nk}\,J_{n}(\sqrt{\lambda_{nk}}r)+B_{nk}\,Y_{n}(\sqrt{\lambda_{nk}}r)\bigr{)}\,e^{in\phi}&(d=2),\\ U_{nkl}(r,\theta,\phi)&=\bigl{(}A_{nk}\,j_{n}(\sqrt{\lambda_{nk}}r)+B_{nk}\,y_{n}(\sqrt{\lambda_{nk}}r)\bigr{)}\,P_{n}^{l}(\cos\theta)\,e^{il\phi}&(d=3),\\ \end{array}$$

$$\underset{\Omega }{\int }d\mathbf{x}|U_m(\mathbf{x}){|}^2=1,$$

$$\bigl{(}\ell_{a}\partial_{r}U_{m}-U_{m}\bigr{)}_{r=a}=\bigl{(}\ell_{b}\partial_{r}U_{m}+U_{m}\bigr{)}_{r=b}=0,$$

$${\mathcal{L}}_1={\partial }_r^2,{\mathcal{L}}_2={\partial }_r^2+\frac{1}{r}{\partial }_r-\frac{n^2}{r^2},{\mathcal{L}}_3={\partial }_r^2+\frac{2}{r}{\partial }_r-\frac{n(n+1)}{r^2}.$$

$$\frac{(z-1{)}^4}{z^2(z+1)}=z-5-\frac{5}{z}+\frac{1}{z^2}+\frac{16}{z+1}=z-5+\sum _{k=1}^2\sum _{j=1}^{n_k}\frac{c_{k,j}}{(z-z_k{)}^j},$$

$$F(z)=\underbrace{F_{0}(z)}_{\textrm{regular part}}+\sum\limits_{k}\Biggl{[}\underbrace{\sum\limits_{j=1}^{n_{k}}\frac{c_{k,j}}{(z-z_{k})^{j}}}_{\textrm{principal parts}}\hskip 5.69054pt+\underbrace{P_{k}(z)}_{\textrm{regularization}}\Biggr{]},$$

$$F(z)=\sum _k\frac{e_k}{z-z_k},$$

$$F(z)=\frac{f(z)}{g(z)}=\frac{f(z_k)}{g^{\prime }(z_k)(z-z_k)}+\mathrm{regular}\mathrm{part}(z\to z_k),$$

$$\frac{f(z)}{g(z)}=F(z)=F_0(z)+\sum _k\left[\frac{f(z_k)}{g^{\prime }(z_k)}\frac{1}{z-z_k}+P_k(z)\right],$$

$$F(z)=\frac{f(z)}{g(z)}=\frac{b_0+b_1(z-z_k)+b_2(z-z_k{)}^2+\dots }{(z-z_k{)}^{n_k}(a_0+a_1(z-z_k)+a_2(z-z_k{)}^2+\dots )}$$

Full text

{history}

Keywords: Laplacian eigenvalues, Bessel functions, spectral sums, summation formulas, heat kernel, diffusion

1 Introduction

The Laplace operator, $\Delta=\nabla^{2}$, plays the central role in mathematical physics and found countless applications in quantum mechanics, acoustics, chemical physics, biology, heat transfer, hydrodynamics, nuclear magnetic resonance, and stochastic processes [1, 2, 3, 4, 5, 6, 7, 8, 9]. In a bounded Euclidean domain $\Omega\subset{\mathbb{R}}^{d}$, the Laplace operator has a discrete spectrum, i.e., a countable number of eigenvalues $\lambda_{m}$ and eigenfunctions $U_{m}$ satisfying

[TABLE]

subject to Dirichlet ($\ell=0$), Neumann ($\ell=\infty$) or Robin ($0<\ell<\infty$) boundary condition on a smooth boundary $\partial\Omega$, with $\partial/\partial n$ being the normal derivative oriented outward the domain $\Omega$ (the Neumann case can formally be understood as the limit $\ell\to\infty$; in particular, the boundary condition (2) reads $\partial U_{m}/\partial n=0$). The eigenvalues are nonnegative and grow to infinity, $0\leq\lambda_{1}\leq\lambda_{2}\leq\ldots\nearrow+\infty$, whereas the eigenfunctions form an orthonormal complete basis in the space $L_{2}(\Omega)$ of square integrable functions. In particular, the eigenvalues and eigenfunctions determine the Green’s functions of homogeneous diffusion and wave equations and thus allow for solving most common boundary value problems [10, 11, 12]. For instance, the spectral decomposition

[TABLE]

determines the Green’s function of the diffusion (or heat) equation,

[TABLE]

subject to the same boundary condition as Eq. (2) and the initial condition $P(\bm{x},t=0|\bm{x}_{0})=\delta(\bm{x}-\bm{x}_{0})$, where $\delta$ is the Dirac distribution, $D$ is the diffusion coefficient, $\bm{x}_{0}$ is a fixed point in $\Omega$, and asterisk denotes the complex conjugate. This Green’s function, also known as heat kernel or propagator, describes various properties of diffusion-controlled processes in reactive media [13, 14, 15, 16]. In turn, the Laplace transform of the heat kernel,

[TABLE]

is the Green’s function of the modified Helmholtz equation,

[TABLE]

subject to the same boundary condition as Eq. (2), and tilde denotes the Laplace transform.

The spectral decomposition (3) and formulas derived from it determine most diffusion-based characteristics. An emblematic example is the heat trace,

[TABLE]

which can be interpreted as the probability density of the return to the starting point $\bm{x}_{0}$, averaged over the starting point. This quantity is also known as the return-to-the-origin probability [17, 18, 19]. Its Laplace transform,

[TABLE]

is one of the simplest spectral sums.

Another important example is the survival probability of a particle diffusing in a domain with partially absorbing boundaries:

[TABLE]

from which the probability density of the (random) reaction time $\tau$ (i.e., the first-passage time to a reaction event) follows as $H(t|\bm{x}_{0})=-\partial S(t|\bm{x}_{0})/\partial t$. If the starting point is uniformly distributed (or, equivalently, the initial concentration of particles is uniform), the volume average of the survival probability yields the NMR signal attenuation due to surface relaxation [20]:

[TABLE]

whereas the volume average of the probability density $H(t|\bm{x}_{0})$ gives the reaction rate on the boundary $\partial\Omega$ of the domain [21]:

[TABLE]

The Laplace transform of these spectral expansions yields several spectral sums that are often used to investigate the short-time asymptotic behavior of the above quantities. For instance, one gets from (9)

[TABLE]

In particular, the mean reaction time reads

[TABLE]

while its volume average is

[TABLE]

Yet another common quantity of interest in the effective diffusion coefficient in a confining domain:

[TABLE]

More generally, the $j$-point correlation function of the particle position $X_{t}$ is

[TABLE]

(with $0<t_{1}<\ldots<t_{j}$), which can be expressed as the $j$-fold spectral sum:

[TABLE]

where

[TABLE]

To evaluate this multiple sum, one needs to perform $j$ Laplace transforms, for each time variable $t_{i}$:

[TABLE]

Similar expressions appear in the perturbative expansion of the diffusion NMR signal attenuation in a linear magnetic field gradient [7, 22, 23]. As the eigenvalues and eigenfunctions are in general not known, only the asymptotic behavior ofspectral sums is usually accessible, e.g., from the Weyl’s asymptotic law for eigenvalues or from the analysis of the heat kernel [24, 25, 26].

However, for a very limited number of simple domains (see examples in Fig. 1), their symmetries allow for the separation of variables, and the eigenvalues and eigenfunctions are known explicitly [1, 2]. An interval, a disk, a circular annulus, a ball and a spherical shell are the most studied examples [27] that play the role of important toy models for understanding diffusion-controlled processes, e.g., in chemical physics and biology. The Laplacian eigenfunctions in these basic domains, defined generally as $\Omega=\{\bm{x}\in{\mathbb{R}}^{d}~{}:~{}a<|\bm{x}|<b\}$ with given radii $0\leq a<b$, are

[TABLE]

where $P_{n}^{l}(z)$ are the associated Legendre polynomials, $J_{n}(z)$ and $Y_{n}(z)$ are the Bessel functions of the first and second kind, and $j_{n}(z)$ and $y_{n}(z)$ are the spherical Bessel functions of the first and second kind. In two and three dimensions, polar and spherical coordinates are respectively used, while double ($n$ and $k$) and triple ($n$, $k$, and $l$) indices are employed to enumerate the eigenfunctions and eigenvalues. As this enumeration is formally equivalent to that by the single index $m$, we will use both schemes and choose one depending on the context. For each eigenfunction $U_{k}$ (or $U_{nk}$, or $U_{nkl}$), the three unknown parameters, the eigenvalue $\lambda_{k}$ (or $\lambda_{nk}$) and the coefficients $A_{k}$ and $B_{k}$ (or $A_{nk}$ and $B_{nk}$), are determined via one $L_{2}$ normalization condition,

[TABLE]

and two boundary conditions in Eq. (2), written explicitly as

[TABLE]

with two given constants $\ell_{a}$ and $\ell_{b}$ (note the sign change in the first relation due to opposite orientation of the radial and normal derivatives). The explicit form of the coefficients $A_{nk}$ and $B_{nk}$, as well as an explicit equation determining the eigenvalues $\lambda_{k}$ (or $\lambda_{nk}$) are well known and will be summarized below. Note that the radial functions in Eq. (18) satisfy the radial eigenvalue problem ${\mathcal{L}}_{d}u+\lambda u=0$ for the differential operators

[TABLE]

In this review, we discuss three strategies to compute explicitly spectral sums involving the eigenvalues $\lambda_{k}$ (or $\lambda_{nk}$): (i) representation of meromorphic functions as sums of partial fractions (Sec. 2), (ii) Fourier-Bessel and Dini series (Sec. 3), and (iii) direct evaluation of Laplace-transformed heat kernels (Sec. 4). The first strategy relies on the corollary of the classical Mittag-Leffler theorem from complex analysis that affirms the existence of a meromorphic function with prescribed poles and allows one to express any meromorphic function as a sum of partial fractions. In this way, one aims at identifying an appropriate meromorphic function to represent a given spectral sum. The second strategy is based on the theory of Fourier-Bessel and Dini series that generalize Fourier series to Bessel functions. Here, one searches for an appropriate function to be decomposed as the Dini series whose form would coincide with the studied spectral sum. The third strategy employs the spectral decomposition (5) of the Laplace-transformed heat kernel (or propagator). As the left-hand side can be evaluated independently by solving the modified Helmholtz equation (6) in simple domains, one obtains a powerful identity for evaluating spectral sums. Remarkably, even though the eigenvalues usually not known explicitly and require a numerical solution of a transcendental equation, spectral sums like Eq. (5) involving such unknown eigenvalues can be often computed exactly and explicitly. Moreover, the particular form (18) of the eigenfunctions often enables analytical computation of the integrals in (9 – 17) (or similar) and thus allows one to derive explicit formulas for such spectral sums.

The ability of computing spectral sums is particularly relevant for understanding the short-time limit of diffusion-controlled quantities. In fact, as each eigenmode in the spectral expansion (3) (or similar) is weighted by the exponential factor $e^{-Dt\lambda_{m}}$, more and more terms are needed to determine the behavior of the heat kernel and other quantities in this limit. As a consequence, the analysis is usually performed in Laplace domain by evaluating the spectral sum (5), determining its asymptotic behavior as $p\to\infty$, and applying Tauberian theorems to get back to the time domain.

At the same time, spectral sums should not be exclusively associated with the short-time limit. For instance, the integral of the heat trace in Eq. (7) yields the sum of inverse eigenvalues, which also corresponds to the limit $p\to 0$ of Eq. (8). In other words, such spectral sums correspond to the long-time limit. For instance, Neuman used the summation formula (60) involving the zeros of the derivative of Bessel and spherical Bessel functions to derive the long-time asymptotic behavior of the diffusion NMR signal [28]. This analysis was further extended in [7] (see also more detailed discussions in [8, 29]).

The explicit summation formulas also help simplifying some analytical solutions. For instance, the separation of variables for the diffusion equation in a right circular cylinder, $\Omega=\{\bm{x}\in{\mathbb{R}}^{2}~{}:~{}|\bm{x}|<R\}\times(a,b)\subset{\mathbb{R}}^{3}$, would naturally lead to double sums over the eigenvalues for an interval $(a,b)$ and the eigenvalues for a disk of radius $R$ (this is also true for a hollow cylinder when the disk is replaced by an annulus). In some cases, one of these sums can be calculated explicitly, thus simplifying expressions and speeding up numerical computations (see an example in [30]). The same argument is applicable to rectangles and rectangular cuboids: $(a_{1},b_{1})\times(a_{2},b_{2})\times\ldots\times(a_{d},b_{d})\subset{\mathbb{R}}^{d}$. For instance, summation formulas were intensively employed to study first-passage times for two particles with reversible target-binding kinetics in a one-dimensional setting, which can be mapped onto diffusion in a rectangle [31].

Rather than dwelling on mathematical proofs and technicalities, our goal is to provide an intuitive practical guide over powerful analytical tools for evaluating spectral sums that could be accessible to non-experts. While many common summation formulas often appearing in applications are summarized, the review does not pretend to be exhaustive.

2 Expansion of meromorphic functions

The first strategy for evaluating spectral sums is based on the corollary of the Mittag-Leffler theorem from complex analysis [32]. This classical theorem describes the properties of meromorphic functions, i.e., functions which are holomorphic (or analytic) at all points except for a discrete set of isolated points (their poles). Any meromorphic function $F(z)$ can be represented as a ratio of two holomorphic functions, $F(z)=f(z)/g(z)$, so that the poles of $F(z)$ are just the zeros of $g(z)$. While general theory refers to an open subset of the complex plane ${\mathbb{C}}$, we restrict discussion to meromorphic functions in ${\mathbb{C}}$, in which case holomorphic functions are called entire functions (i.e., they can be represented by a power series uniformly converging on compact sets). Polynomials, exponential function, sine and cosine, Bessel functions $J_{n}(z)$ of the first kind of an integer order $n$, spherical Bessel functions $j_{n}(z)$ of the first kind of an integer order $n$, their finite sums, products and compositions, derivatives and integrals are examples of entire functions. In contrast, square root, logarithm and functions with essential singularities (e.g., $e^{1/z}$) are not entire functions. The Mittag-Leffler theorem claims the existence of a meromorphic function with prescribed poles. As a consequence, any meromorphic function can be decomposed as a sum of partial fractions.

In this light, a meromorphic function is a natural extension of rational functions (i.e., ratio of two polynomials). Any rational function can be decomposed into a sum of partial fractions, e.g.,

[TABLE]

with the regular part $z-5$, the poles $z_{1}=0$ and $z_{2}=-1$ of order $n_{1}=2$ and $n_{2}=1$, and the weights $c_{1,1}=-5$, $c_{1,2}=1$ and $c_{2,1}=16$. The Mittag-Leffler theorem extends such decompositions to meromorphic functions [32] (p. 297): any meromorphic function $F(z)$ with a prescribed sequence of poles $\{z_{k}\}$ (with orders $n_{k}$) can be represented as

[TABLE]

where the regular part $F_{0}(z)$ is an entire function, each principal part is the sum of the partial fractions near the pole $z_{k}$ with weights $c_{k,j}$, and $P_{k}(z)$ are polynomials. The theorem states that a meromorphic function is essentially determined by its principal parts, up to an entire function. Since a meromorphic function may have infinite number of poles (i.e., the sum over $k$ can be infinite), adding polynomials $P_{k}(z)$ to the principal parts in the right-hand side of Eq. (23) may be needed to ensure the convergence of the series. For instance, the decomposition of the meromorphic function in Eq. (22) includes the entire function $F_{0}(z)=z-5$ and two principal parts.

The strategy for computing spectral sums consists in finding (or guessing) an appropriate meromorphic function. Let us first consider spectral sums that can be reduced to

[TABLE]

with a prescribed set of simple poles $\{z_{k}\}$ (i.e., $n_{k}=1$) and coefficients $\{e_{k}\}$. In the first step, one searches for an entire function $g(z)$ whose complete set of zeros is the set $\{z_{k}\}$. In practice, this is a simple step because the eigenvalues are often expressed in terms of zeros of some explicit entire functions. Next, one searches for another entire function $f(z)$ such that $f(z_{k})\neq 0$ for all $k$, and $F(z)=f(z)/g(z)$. The zeros $\{z_{k}\}$ are therefore the poles of the meromorphic function $F(z)=f(z)/g(z)$. Around $z\approx z_{k}$, one has $g(z)=g^{\prime}(z_{k})(z-z_{k})+O((z-z_{k})^{2})$, from which

[TABLE]

where prime denotes the derivative with respect to the argument of the function. The first term in (25) is the principal part of the function $F(z)$. The function $f(z)$ should satisfy $f(z_{k})/g^{\prime}(z_{k})=e_{k}$ for each $k$. Applying the corollary of the Mittag-Leffler theorem, one has

[TABLE]

where an entire function $F_{0}(z)$ (accounting for remaining regular parts) and polynomials $P_{k}(z)$ (ensuring the convergence) have to be found. In many cases, the “indefinite” ingredients $F_{0}(z)$ and $P_{k}(z)$ can be guessed or computed. A practical strategy may consist in computing numerically the truncated series in the right-hand side and checking this relation (paying attention to the convergence which may be critical for numerical analysis).

In the general case when a pole $z_{k}$ has the order $n_{k}$, the series expansion around $z\approx z_{k}$ reads

[TABLE]

where

[TABLE]

and the indices $i_{0}$, … $i_{m}$, $j$ are nonnegative. For instance, one has

[TABLE]

The coefficients $a_{m}$ and $b_{m}$ are related to the derivatives of the functions $f(z)$ and $g(z)$ at points $z_{k}$:

[TABLE]

In the particular case $e_{k}=n_{k}$, one formally gets an expansion of the logarithmic derivative of $g(z)$:

[TABLE]

In a typical case of an infinite sequence of simple zeros (i.e., $n_{k}=1$)111 This mathematical result can be significantly extended as described by [32]. Here we do not discuss convergence of the series in the right-hand side of Eq. (30), especially since the most attention will be focused on the function $\eta_{n}(z)$ defined by the convergent series in Eq. (43)., this relation can be made rigorous by requiring

[TABLE]

The relation (30) can be formally obtained by representing $g(z)$ as a product over its zeros,

[TABLE]

(with a constant $g_{0}$), taking derivative and dividing it by $g(z)$.

2.1 Basic example: interval with absorbing endpoints

To illustrate this strategy, let us compute the Laplace-transformed heat trace $\tilde{Z}(p)$ in the simple case of the Laplace operator on the interval of length $R$ with Dirichlet boundary conditions at both endpoints ($\ell=0$), for which $\lambda_{m}=\pi^{2}m^{2}/R^{2}$, with $m=1,2,3,\ldots$:

[TABLE]

where we set $-z^{2}=pR^{2}/D$. One needs therefore to compute the last sum with simple poles at $z_{k}=\pi k$ that are all the zeros of the entire function $g(z)=\sin z$. Applying Eq. (30), one gets

[TABLE]

which is still formal because the series in the right-hand side is divergent. To regularize this series, one can subtract $1/(\pi k)$ (for $k\neq 0$) whose contribution would be canceled anyway due to the similar term for $-k$:

[TABLE]

where the term with $k=0$ was moved to the left-hand side. The subtracted term $1/(\pi k)$ is precisely the polynomial $P_{k}(z)$ (here, just a constant) in the general representation (23). As the series in Eq. (2.1) is convergent, this regularization is actually not needed, and the substitution of Eq. (33) into Eq. (2.1) yields two equivalent classical summation formulas

[TABLE]

from which

[TABLE]

A similar derivation with the function $f(z)=1$ in (26) yields

[TABLE]

The explicit form (35) of the Laplace-transformed heat trace $\tilde{Z}(p)$ facilitates the analysis of the short-time and long-time asymptotic behaviors of the return-to-the-origin probability discussed in Sec. 1. The knowledge of $\tilde{Z}(p)$ also allows one to compute some other spectral series, in particular,

[TABLE]

For instance, one easily retrieves the even integer values of the Riemann zeta-function,

[TABLE]

as

[TABLE]

for instance

[TABLE]

As shown below, many other spectral series can be reduced to (35). While the same technique is applicable to a more general case of an interval with partially absorbing endpoints ($\ell>0$), we will discuss this case in Sec. 4.2.

2.2 Disk and ball

In the case of a disk or a ball of radius $R$, the separation of variables in polar or spherical coordinates leads again to an explicit form of the eigenfunctions of the Laplace operator. The eigenvalues $\lambda_{m}$ can be written as $\alpha_{nk}^{2}/R^{2}$ with double indices $n=0,1,2,\ldots$ and $k=1,2,\ldots$, where $\alpha_{nk}$ are strictly positive simple zeros of an explicit entire function:

[TABLE]

where $h=R/\ell$ is a nonnegative parameter (see Table 1). In this case, the computation should be performed separately for each $n$. Note that the zero eigenvalue appears only for Neumann boundary condition and will be treated separately.

To illustrate this technique, we focus on computing the Laplace-transformed heat kernel trace $\tilde{Z}(p)$, while more sophisticated summation formulas will be reported in Sec. 4 (see also Table 3). For convenience, we rewrite $\tilde{Z}(p)$ as

[TABLE]

where

[TABLE]

and $z=R\sqrt{p/D}$. According to the form (18) of Laplacian eigenfunctions, the function $\eta_{n}(z)$ can be associated with the $n$-th Fourier (or spherical harmonic) mode. Note that $-\alpha_{nk}$ are also zeros of the function $g_{n}(z)$ due to its symmetry. Moreover, the function $g_{n}(z)$ vanishes at $z=0$ for any $n>0$ but also for $n=0$ in the case of Neumann boundary condition. Since the set $\{0,\pm\alpha_{nk}\}$ contains all zeros of $g_{n}(z)$, one can apply the corollary of the Mittag-Leffler theorem to get

[TABLE]

where the second term accounts for the zero at [math], which is excluded from $\eta_{n}(z)$. For Neumann boundary condition, $z=0$ is also a simple zero of $g_{0}(z)$, and as it was accounted twice, it should be subtracted:

[TABLE]

Note that the condition (31) is satisfied thanks to the asymptotic behavior $\alpha_{nk}\propto k$ (for fixed $n$). As shown in [29], $\eta_{n}(z)$ from (44) with $n>0$ can also be written for a ball in any dimension and boundary condition as

[TABLE]

where the function $\psi_{n,d}(z)$ is defined by the recurrent relation

[TABLE]

and

[TABLE]

Note that the above representation with $d=1$ is also valid for a symmetric interval $(-R,R)$, see details in [29].

As an example, let us consider a disk with Dirichlet boundary condition ($\ell=0$), for which Eq. (44) implies the summation formula:

[TABLE]

This is a particular case of the summation formula (see [33], p. 498)

[TABLE]

where $\pm\alpha_{\nu k}$ are zeros of $z^{-\nu}J_{\nu}(z)$. Evaluating the derivatives of this identity at $z=0$, Rayleigh found other summation formulas [34], e.g.,

[TABLE]

(see also [33, 35, 36]). At $\nu=1/2$, $z^{-\nu}J_{\nu}(z)=\sqrt{2/\pi}\,\frac{\sin z}{z}$ so that Eqs. (50, 51) are reduced to Eqs. (2.1, 40) for the interval. Similarly, one gets a representation of the inverse of the Bessel function as

[TABLE]

where $\mathcal{P}_{n}$ is the polynomial of degree $n$ that accounts for the pole at [math]:

[TABLE]

The explicit form of this polynomial can be obtained by expanding $z^{n}/J_{n}(z)$ in a Taylor series around [math] and keeping the terms up to the $n$-th order.

Once $\eta_{n}(z)$ is known, one can establish a number of useful relations for more complicated sums involving the Laplacian eigenvalues. Such spectral sums naturally appear when the quantity of interest depends on several time variables, $t_{1},\ldots,t_{j}$, as, e.g., the $j$-point correlation function in (16). In this case, the Laplace transform is applied to each time variable, and the Laplace-transformed quantity in (17) depends on multiple conjugate parameters $p_{1},\ldots,p_{j}$. First, one can compute the sums

[TABLE]

iteratively, applying the identity

[TABLE]

with $\eta_{n}^{(1)}(z)\equiv\eta_{n}(z)$. In particular,

[TABLE]

Multiple differentiation of Eq. (54) with respect to the variables $s_{i}=z^{2}_{i}$ further extends the set of useful relations for any positive integers $\sigma_{1}$, …, $\sigma_{j}$:

[TABLE]

For instance, for $j=1$ and $\sigma=0,1,2,\ldots$, one gets

[TABLE]

from which the sum of inverse powers of eigenvalues (as in Eq. (51)) can be evaluated by setting $z\to 0$ (or $s\to 0$). For instance, setting $\sigma=0$, one retrieves (see, e.g., [38] and references therein):

[TABLE]

Although cumbersome, these expressions allow one to rigorously compute many spectral sums and related quantities involving the Laplacian eigenvalues (see [29, 37] for examples and applications).

As an example, we discuss the Calogero’s formula for the zeros $\alpha_{\nu k}$ of the function $g_{n}(z)=zJ^{\prime}_{\nu}(z)+hJ_{\nu}(z)$ with any real $\nu>-1$ [39, 40]

[TABLE]

This formula can be easily derived for integer $\nu=n$ by evaluating the limit:

[TABLE]

Using the explicit form (44) of the function $\eta_{n}(z)$, one can compute the limit by expanding $g_{n}(z)$ and $g^{\prime}_{n}(z)$ in a Taylor series around $z=\alpha_{nj}$. After simplifications, one gets

[TABLE]

Using the Bessel equation,

[TABLE]

one can express in two dimensions

[TABLE]

Since $\alpha_{nj}$ is a zero of $g_{n}(z)$, one finally gets

[TABLE]

from which the Calogero’s formula (62) follows for an integer $\nu=n$. Similarly, one can deduce the summation formulas for higher-order inverse differences between eigenvalues, see [41, 42, 38, 43]. Interestingly, Calogero and co-workers considered (62) as an infinite system of nonlinear equations that determine the zeros $\alpha_{\nu j}$, together with the normalization relation (61) extended to non-integer $\nu$.

In three dimensions, one uses the modified Bessel equation to express

[TABLE]

from which

[TABLE]

and thus

[TABLE]

This is an extension of the Calogero’s formula to the zeros $\alpha_{nk}$ of the function $g_{n}(z)=zj^{\prime}_{n}(z)+hj_{n}(z)$ involving spherical Bessel function and its derivative.

As the review is focused on the eigenvalues of the Laplace operator in rotation-invariant domains (here, a disk or a ball), we restricted the discussion to Bessel function $J_{n}(z)$ and spherical Bessel function $j_{n}(z)$ of integer index $n$. However, as illustrated above, many summation formulas of this section are valid if $n$ is replaced by a real index $\nu>-1$. In this more general setting, one needs to consider $z^{-\nu}J_{\nu}(z)$ or $z^{-\nu}j_{\nu}(z)$, which are again entire functions. Note that the zeros of Bessel and spherical Bessel functions with non-integer indices are naturally related to the Laplacian eigenvalues for circular and spherical sectors, respectively [1]. In spite of their practical interest, we do not discuss the summation formulas for these domains. Finally, for $d$-dimensional balls with $d>3$, one needs to consider the zeros of ultraspherical Bessel functions $z^{1-d/2}J_{d/2-1+n}(z)$.

3 Fourier-Bessel and Dini series

The second strategy for evaluating spectral sums relies on Fourier-Bessel and Dini series which extend Fourier series to Bessel functions. The Fourier-Bessel series involve the zeros of a Bessel function $J_{\nu}(z)$ or of its derivative $J^{\prime}_{\nu}(z)$, whereas Dini series involve the zeros of a linear combination $zJ^{\prime}_{\nu}(z)+hJ_{\nu}(z)$ with a constant $h$. As these zeros naturally appear from Dirichlet, Neumann and Robin boundary conditions on Laplacian eigenfunctions (see Table 1), we refer to them as Dirichlet, Neumann and Robin cases.

A continuous function $f(x)$ on $(0,1)$, integrable on $(0,1)$ with the weight $x^{1/2}$ and of limited total variation, admits [33, 35]:

• •

the Fourier-Bessel expansion (for $\nu\geq-1/2$)

[TABLE]

where $\alpha_{\nu k}$ are the positive zeros of $J_{\nu}(z)$ (Dirichlet case);

• •

the Fourier-Bessel expansion (for $\nu>0$)

[TABLE]

where $\alpha_{\nu k}$ are the positive zeros of $J^{\prime}_{\nu}(z)$ (Neumann case);

• •

the Dini expansion (for $\nu\geq-1/2$ and $\nu+h>0$)

[TABLE]

where $\alpha_{\nu k}$ are the positive zeros of $zJ^{\prime}_{\nu}(z)+hJ_{\nu}(z)$ (Robin case).

Note that the continuity assumption and some other constraints can be weakened (see details in [33]). If the Fourier series can be seen as spectral expansions on over Laplacian eigenfunctions on an interval, Fourier-Bessel and Dini series are essentially the spectral expansions over radial parts of the Laplacian eigenfunctions on the disk (see also Sec. 4).

As pointed out by Watson [33], the coefficients of these expansions, given by integrals with the function $f(y)$, are rarely known explicitly. Among such examples, one has

[TABLE]

for $\alpha_{\nu k}$ being zeros of $J_{\nu}(z)$, and

[TABLE]

for $\alpha_{\nu k}$ being zeros of $zJ^{\prime}_{\nu}(z)+hJ_{\nu}(z)$ (under the additional condition $\nu+h>0$).

Relying on these expansions, Sneddon derived several kinds of summation formulas for some positive integer $m$ [35]

• •

Dirichlet case: $\alpha_{\nu k}$ are positive zeros of $J_{\nu}(z)$

[TABLE]

and

[TABLE]

(note that (77) reproduces Rayleigh’s formulas (51)).

• •

Neumann case: $\alpha_{\nu k}$ are positive zeros of $J^{\prime}_{\nu}(z)$

[TABLE]

and

[TABLE]

• •

Robin case: $\alpha_{\nu k}$ are positive zeros of $zJ^{\prime}_{\nu}(z)+hJ_{\nu}(z)$

[TABLE]

Note that the Sneddon’s technique can be used to derive such formulas for larger values of $m$.

Sneddon also gave some examples of other series, e.g.,

[TABLE]

for the Dirichlet case, and

[TABLE]

for the Robin case.

In summary, Fourier-Bessel and Dini series present a powerful tool for deriving spectral sums over the zeros of Bessel functions. However, as such series involve only Bessel functions of the first kind, $J_{\nu}(z)$, they are not directly applicable to circular annuli or spherical shells, for which the Laplacian eigenfunctions also include Bessel functions of the second kind, $Y_{\nu}(z)$. In the next section, we discuss a different method based on the direct computation of the Laplace-transformed heat kernel. Here, the radial coordinate will be restricted on an interval $(a,b)$, and the Sturm-Liouville theory will ensure the completeness of a basis formed by linear combinations of $J_{\nu}(z)$ and $Y_{\nu}(z)$.

4 Laplace-transformed heat kernel

The two strategies discussed in the previous sections are well adapted for evaluating spectral sums involving the eigenvalues of the Laplace operator on an interval, a disk, and a ball. In particular, the associated eigenfunctions are expressed through entire functions: sine, cosine, Bessel and spherical Bessel functions of the first kind. However, these strategies are not suitable for circular annuli and spherical shells because Bessel and spherical Bessel functions of the second kind, $Y_{n}(z)$ and $y_{n}(z)$, are not entire functions. Similarly, Fourier-Bessel and Dini series do not operate with $Y_{n}(z)$ and $y_{n}(z)$. For this reason, we discuss here in more detail the third strategy for evaluating spectral sums that is also applicable to the Laplace operator in circular annuli and spherical shells. In these rotation-invariant domains, the separation of variables reduces the computation of the Laplacian eigenvalues to the analysis of radial second-order differential operators ${\mathcal{L}}_{d}$ in Eq. (21).

4.1 General case

Let us first recall the basic elements of the Sturm-Liouville theory can be found in most textbooks on differential equations and mathematical physics, e.g., [44, 45]. Let ${\mathcal{L}}$ be a self-adjoint second-order differential operator on an interval $(a,b)$, with prescribed boundary conditions at endpoints $a$ and $b$. As we will focus on the particular well-studied case of the operators ${\mathcal{L}}_{d}$ in Eq. (21), we skip mathematical details and eventual constraints known from the spectral theory of ordinary differential operators [46]. In particular, we postulate that the spectrum of the operator ${\mathcal{L}}$ is discrete that is clearly the case in our setting.

The third strategy relies on computing the Green’s function $G_{q}(r,r_{0})$, i.e., the resolvent of the operator ${\mathcal{L}}$ (acting on $r$) that obeys the equation

[TABLE]

where $\omega(r_{0})$ is a given positive weighting function (see below). On one hand, the standard construction of $G_{q}(r,r_{0})$ involves two linearly independent solutions $v_{q}^{a}$ and $v_{q}^{b}$ of the homogeneous equation $(q^{2}-{\mathcal{L}})v_{q}^{a,b}=0$ such that $v_{q}^{a}$ satisfies the boundary condition at $r=a$ and $v_{q}^{b}$ satisfies the boundary condition at $r=b$. From these two solutions, one can construct the Green’s function as

[TABLE]

where the coefficients $A$ and $B$ are determined by matching these parts at the point $r=r_{0}$ to ensure the continuity of $G_{q}(r,r_{0})$ and the unit jump of its derivative. One gets thus

[TABLE]

On the other hand, the function $v_{i\sqrt{\lambda}}^{a}(r)$ with a nonnegative $\lambda$ is a natural candidate to be an eigenfunction of the operator ${\mathcal{L}}$ (a similar construction is applicable to $v_{i\sqrt{\lambda}}^{b}(r)$). Indeed, this function obeys the eigenvalue equation, ${\mathcal{L}}u+\lambda u=0$, with the imposed boundary condition at $r=a$. The remaining boundary condition at $r=b$ determines the (infinite) set of eigenvalues of this operator that we denote as $\lambda_{k}$, with the index $k=1,2,\ldots$. The corresponding eigenfunctions are then

[TABLE]

where the prefactor $i$ was introduced for convenience. Denoting

[TABLE]

the $L_{2}$-normalization coefficient of $u_{k}$ with the weighting function $\omega(r)$, one can express the Green’s function via the spectral decomposition over the eigenfunctions of the operator ${\mathcal{L}}$

[TABLE]

Equating relations (86, 89), one gets an explicit summation formula over the eigenvalues $\lambda_{k}$:

[TABLE]

Multiplying both sides of this formula by an integrable function $f(r,r_{0})$ (with the weighting function $\sqrt{\omega(r)}$) and integrating over $r$ and $r_{0}$, one can deduce many other summation formulas.

In the following, we focus on Robin boundary condition (20), which reads for the Green’s function as

[TABLE]

If ${\mathcal{I}}(r)$ and ${\mathcal{K}}(r)$ denote two independent solutions of the differential equation ${\mathcal{L}}u=u$, then the functions $v_{q}^{a,b}$ can be expressed as their linear combinations:

[TABLE]

with coefficients

[TABLE]

In particular, the denominator in the right-hand side of Eq. (90) becomes

[TABLE]

where

[TABLE]

and

[TABLE]

is the Wronskian of two solutions.

The Robin boundary condition for an eigenfunction $-iv_{i\sqrt{\lambda}}^{a}(r)$ at $r=b$ yields an equation on the eigenvalues

[TABLE]

Denoting the solutions of this equation as $q_{k}$ (with $k=1,2,\ldots$), one gets the eigenvalues and eigenfunctions

[TABLE]

The general relation (90) can thus be written as

[TABLE]

(the other case $r\geq r_{0}$ is obtained in a similar way). As the operator ${\mathcal{L}}$ with Robin boundary conditions is self-adjoint, its eigenfunctions $u_{k}(r)$ form a complete basis in the space of square-integrable functions on $(a,b)$ with the weighting function $\omega(r_{0})$. The completeness relation can be formally written as

[TABLE]

4.2 One-dimensional case

Due to the translational invariance in one dimension, diffusion on the interval $(a,b)$ is fully equivalent to that on the shifted interval $(0,b-a)$. Without loss of generality, we set then $a=0$ and consider diffusion on $(0,b)$, which is described by the second-order derivative, ${\mathcal{L}}_{1}=\partial_{r}^{2}$. Two independent solutions of the equation $u^{\prime\prime}-u=0$ read

[TABLE]

so that

[TABLE]

and

[TABLE]

from which

[TABLE]

The eigenvalues are determined as $\lambda_{k}=q_{k}^{2}$, with $q_{k}$ being the positive solutions of the equation

[TABLE]

while the eigenfunctions from Eq. (87) have a form

[TABLE]

Integrating $u_{k}^{2}$ with the weighting function $\omega(r_{0})=1$ and using trigonometric identities, one finds

[TABLE]

where

[TABLE]

are dimensionless quantities. It is also convenient to rewrite Eq. (110) as an equation on $\alpha_{k}$:

[TABLE]

As the function in the left-hand side is continuous and monotonously decreasing on an interval $(\pi(k-1),\pi k)$ from $+\infty$ to $-\infty$ (that can be checked by computing its derivative), there is a single zero $\alpha_{k}$ on each such interval, implying $\pi(k-1)\leq\alpha_{k}\leq\pi k$ for $k=1,2,3,\ldots$.

Using Eq. (97) and ${\mathcal{W}}(z)=1$, the summation formula (102) becomes

[TABLE]

for $0\leq r\leq r_{0}\leq b$. Note that the left-hand side is precisely the Laplace-transformed heat kernel:

[TABLE]

so that Eq. (115) yields its analytical form, with $q=\sqrt{p/D}$. From the identity (115), one can derive various summation formulas. For instance, setting $r_{0}=r$ and integrating over $r$ from [math] to $b$, one finds

[TABLE]

which is proportional to the Laplace-transformed heat trace $\tilde{Z}(p)$.

The particular cases of the general identity (115) corresponding to nine combinations of boundary conditions at two endpoints are summarized in Table 4.2. For instance, in the Neumann-Neumann case ($\ell_{0}=\ell_{b}=\infty$), for which $\alpha_{k}=\pi(k-1)$ ($k=1,2,\ldots$), one has

[TABLE]

where we fixed $b=1$. Setting $x=x_{0}=0$, one retrieves Eq. (34b). In the particular cases of Dirichlet-Dirichlet and Neumann-Neumann conditions, many summation formulas can be deduced from standard Fourier series, e.g.,

[TABLE]