Fast spherical Bessel transform via fast Fourier transform and recurrence formula
Introduction
Numerical evaluation of integrals containing the spherical Bessel function is of importance in many fields of computational science and engineering since the spherical Bessel function is often used as the eigenfunction for spherical coordinate systems. The integrals are known as the spherical Bessel transform (SBT) which is classified into a more general family of the Hankel or Fourier–Bessel transforms. The SBT is involved in many physical problems such as the scattering in atomic or nuclear systems [1], [2], the simulation of the cosmic microwave background [3], and the interaction of electrons in molecules and crystals [4], [5].
Several computational methods have been developed for the Hankel transforms [6], [7], [8]. However, not many of them can be applied for SBT. There have been proposed three different approaches for SBT for general order: (a) recasting the transform integral as a convolution integral by changing the coordinate variables [9], [10], [11], (b) expansion in terms of a series of Fourier cosine and sine transforms by the trigonometric expansion of the spherical Bessel function [12], and (c) the discrete Bessel transform method which describes SBT as an orthogonal transform [13]. The approach (a) is quite fast since it utilizes the fast Fourier transform (FFT) algorithm. The computation time actually scales as , where N is the number of quadrature points. It, however, has a major drawback caused by the logarithmic grid that almost all the grid points are located in close proximity of the origin. The computation time of the approach (b) which also uses FFT scales as , where ℓ is the order of transform. Therefore, it becomes slower for the higher order transform. In addition to that, since it requires the integrand to be multiplied by inverse powers of the radial coordinate, the high order transforms may become unstable. The computation time of the approach (c) scales as . The quadrature points are located at each zero of the spherical Bessel function. The optimized selection of the quadrature points enables us to use a small number of N while keeping the accuracy of the computation. However, when consecutive transforms with different orders are required, it may become a minor trouble that the optimized quadrature points differ depending on the order of transform.
In this paper, we propose a new method for the numerical SBT which uses a linear coordinate grid. The transform is decomposed into the Fourier transforms and the numerical integrations which can be evaluated recursively. The computation time for the present method scales as with overhead for the numerical integration which scales as N. The linear coordinate grid prevents us from troubles caused by the non-uniform or order-dependent grid points. If the considered problem requires to transform a function with various orders, the present method has further the advantage that the results of the most time consuming calculations (i.e. the Fourier transforms and the integrations) for a transform with a certain order ℓ can also be used for transforms with any order less than ℓ.
Section snippets
Formulation
We are interested in the ℓ-th order SBT of a function which is defined as follows: where is the spherical Bessel function of the first kind. The integral representation of the spherical Bessel function is given by where is the Legendre polynomials Here, is the double factorial with an exceptional definition
Implementation
In order to avoid problems arising from at the origin, we define the r- and k-grid points at half-interval shifted positions as follows: At each point of , the integral is divided into its segments where the segments are defined as By defining the sum of the segments as , the following simple recurrence formula is obtained:
Computation results
The present method has been applied in transforming the atomic orbital wave functions which are used in quantum chemistry and condensed matter physics. Three types of atomic orbitals have been examined: the Gaussian-type orbital (GTO), the Slater-type orbital (STO), and the numerically defined pseudo-atomic orbital (PAO) functions. The first two functions are good examples to investigate in detail the numerical error accompanied by the present method since the exact form of the transformed
Summary
In conclusion, we have proposed and demonstrated a new method for the numerical evaluation of SBT. Sufficiently accurate results are obtained in transforming analytic atomic orbital functions. Even the non-analytic PAO functions can be transformed by the present method with a practically acceptable accuracy. Application of the present method to evaluate the electron–electron repulsion integrals is currently in progress by the authors. A similar framework could also be developed for the Hankel
Acknowledgements
This work was partly supported by CREST-JST and the Next Generation Super Computing Project, Nanoscience Program, MEXT, Japan.
References (15)
- et al.
Calculations for quasi-elastic scattering , and at 460 MeV
Nucl. Phys. A
(1972) An algorithms for the Fourier Bessel transform
Comput. Phys. Comm.
(1981)Numerical evaluation of the Hankel transform
Comput. Phys. Comm.
(1999)- et al.
Numerical evaluation of the Hankel transform by using linear Legendre multi-wavelets
Comput. Phys. Comm.
(2008) Numerical Fourier and Bessel transforms in logarithmic variables
J. Comput. Phys.
(1978)NumSBT: A subroutine for calculating spherical Bessel transforms numericaly
Comput. Phys. Comm.
(2009)- et al.
Numerical evaluation of spherical Bessel transforms via fast Fourier transforms
J. Comput. Phys.
(1992)
Cited by (18)
Simulation of atomically resolved elemental maps with a multislice algorithm for relativistic electrons
2019, Advances in Imaging and Electron PhysicsCitation Excerpt :As can be seen in Fig. 8 and 9, the accuracy of our implementation of the Hankel transform depends on the order of the Hankel transform l, the transformed function and its sampling. The obtained accuracy is comparable to the results reported by Toyoda and Ozaki (2010) and should be sufficient for our purposes as long as the transformed functions are sufficiently sampled and converge towards zero within the sampling interval. Usually, investigations with an electron microscope are carried out at room temperature.
Numerical calculation of the Spherical Bessel Transform from Gaussian quadrature in the complex-plane
2018, Computational and Theoretical ChemistrySpherical Bessel transform via exponential sum approximation of spherical Bessel function
2018, Journal of Computational PhysicsOrthogonal fast spherical Bessel transform on uniform grid
2017, Computer Physics CommunicationsCitation Excerpt :The discrete spherical Bessel transform (DSBT) arises in a number of applications, such as, e.g., the analysis of the cosmic microwave background [1], the numerical solution of the differential equations [2–4], and the numerical evaluation of multi-center integrals [5,6]. Many different SBT algorithms have been proposed so far [7–10]. But none of them possess all of the advantages of their trigonometric progenitor, namely the fast Fourier transform (FFT).
Automated segmentation and shape characterization of volumetric data
2014, NeuroImageCitation Excerpt :The angular spherical harmonic decomposition was implemented using a divide and conquer approach and a fast Legendre transform (Driscoll and Healy, 1994; Healy et al., 1996, 2004; Mohlenkamp, 1999; Rokhlin and Tygert, 2006). For the radial spherical Bessel transform evaluation several different fast implementations are also available (Bisseling and Kosloff, 1985; Koval and Talman, 2010; Pettitt et al., 1993; Sharafeddin et al., 1992; Talman, 1978, 2009; Toyoda and Ozaki, 2010). In our implementation we followed Toyoda and Ozaki (2010) but replaced their recurrence formula with asymptotic expansion of the integral (see Appendix A).
LIBERI: Library for numerical evaluation of electron-repulsion integrals
2010, Computer Physics CommunicationsCitation Excerpt :In the method, an SBT is evaluated through two consecutive Fourier transforms, by employing a logarithmic radial mesh. The other method for numerical SBT has recently been developed by the authors [14] which shows performance comparable to or even better in some specific cases than the Siegman–Talman method. This method uses a linearly spaced radial mesh.