S. K. H. Auluck

## On the Integral of the Product of Three Bessel Functions over an Infinite Domain NB   CDF   PDF

### Fourier-Space Representation of Nonlinear Dynamics of Continuous Media in Cylindrical Geometry

Fourier-space representation of the partial differential equations describing nonlinear dynamics of continuous media in cylindrical geometry can be achieved using Chandrasekhar-Kendall (C-K) functions defined over infinite domain as an orthogonal basis for solenoidal vector fields and their generating function and its gradient as orthogonal bases for scalar and irrotational vector fields, respectively. All differential and integral operations involved in translating the partial differential equations into transform space are then carried out on the basis functions, leaving a set of time evolution equations, which describe the rate of change of the spectral coefficient of an evolving mode in terms of an aggregate effect of pairs of interacting modes computed as an integral over a product of spectral coefficients of two physical quantities along with a kernel, which involves the following integral:

involving the product of three Bessel functions of the first kind of integer order. This article looks at this integral’s properties using a semi-empirical approach supported by numerical experiments. It is shown that this integral has well-characterized singular behavior. Significant reduction in computational complexity is possible using the proposed empirical approximation to this integral.

### Introduction

Fourier-space representation of the nonlinear dynamics of continuous media is useful in studies on turbulent or self-organizing behavior that looks at the transport of energy across a wide range of scale lengths. Many problems in physics possess a native cylindrical geometry, and for such cases Chandrasekhar-Kendall (C-K) functions [1] (curl eigenfunctions) provide a complete orthogonal basis [2] for solenoidal fields, under appropriate boundary conditions. C-K functions defined over a finite cylinder have been used as an orthogonal basis for spatial-spectral expansion methods in plasma physics, usually in the incompressible regime [3, 4, 5], which neglects the role of irrotational fields in plasma dynamics. However, there are important physical effects [6] that involve a nonlinear interaction between solenoidal and irrotational fields, which can be studied only if both irrotational and solenoidal fields are taken into account on an equal footing.

The limitation to incompressible flows (or solenoidal fields in general) in spatial-spectral methods referred to above arises from the fact that the radial component of C-K functions and the radial derivative of their generating function (scalar potential solution of the Helmholtz equation) cannot both be zero on the same finite cylinder, so that the gradient of the generating function does not provide an orthogonal basis for irrotational flow. In contrast, when defined over an infinite domain, the generating function and its gradient do serve as orthogonal bases for scalar and irrotational vector fields, respectively. This is equivalent to Fourier-space representation in cylindrical geometry.

Although the mathematics of Fourier-space representation deals with spatial dimensions of the system and Fourier-space coordinates ranging from zero to infinity, no physical system has infinite size, and physical models of continuous media often stipulate a lower bound on scale lengths as a condition for their applicability. An infinite domain thus has a utilitarian meaning within the context of a physical theory: the upper and lower limits of integral transforms used in a physical theory are physically meaningful large and small numbers; mathematical zero and infinity are used to obtain the limit of a sequence of calculations with increasingly smaller and larger limits of integration in the interest of obtaining a tractable model and experimentally verifiable predictions [7].

The C-K functions, eigenfunctions of the curl operator in cylindrical geometry (defined by unit vectors , , and with coordinates ) are defined as , . The set of eigenfunctions, labeled by the four parameters , , , , is

The parameter labels circular polarizations and is the azimuthal mode number; the radial mode number and axial mode number are related to the magnitude of the eigenvalue by the relation . The orthogonality relations over an infinite cylinder are:

For the infinite domain, the generating function of C-K functions

satisfies the relation

and its gradient satisfies the following relations

The orthogonality between the gradient of the generating function and the C-K functions is an exact consequence of their structure and does not depend on their asymptotic behavior. Partial differential equations of the nonlinear dynamics of continuous media can be translated into cylindrical Fourier space by performing all differential and integral operations with respect to spatial variables on the basis functions, leaving a set of time evolution equations. As an example, the equation of continuity,

can be written in transform space as

with the following shorthand notation: . The spectral coefficient of density is , that of solenoidal velocity is , and that of irrotational velocity is . The functions of mode numbers used in equation (7) are defined as

Time-evolution equations like (7) describe the rate of change of the spectral coefficient of an evolving mode on the left-hand side (LHS) in terms of an aggregate effect of pairs of interacting modes on the right-hand side (RHS) computed as an integral over a product of spectral coefficients of two physical quantities along with a kernel. The physics contained in the vector and differential operators of the partial differential equations is transferred to the functions of mode numbers similar to , which involve the following integral over an infinite domain consisting of products of three Bessel functions of the first kind of integer order,

This integral, referred to in the next section as the triple-Bessel (or 3B) integral, belongs to a class of integrals discussed by Watson [9, 10], whose analytical theory is not known and is the subject of this article. From (7), its properties as a distribution are seen to be relevant, but not its numerical value as a conventional Riemann integral, which may not exist for some values of its arguments. The comment made above concerning the utilitarian meaning of infinity in physical theory applies in this context, where properties of a sequence of integrals similar to (11) with increasingly larger but finite upper limit can be investigated usefully. The upper limit of integration is therefore mentioned as an argument in (11). Whenever the argument is not mentioned, it is implied to be infinity.

Transformation of the equations of continuity and momentum conservation in the two-fluid model of plasma into cylindrical Fourier space involves 115 instances of the 3B integral in 19 expressions similar to those described in equations (8)-(10) [8]. Investigations of plasma phenomena involving nonlinear interaction between solenoidal and irrotational fields, such as self-organization or turbulence in the two-fluid plasma model, would benefit considerably if an adequate approximation of the 3B integral becomes available. This is the motivation of the present work.

### Properties of the Triple-Bessel Integral

The 3B integral is trivially invariant under permutation of pairs of arguments , , . If the arguments , , are scaled as , , , then the integral scales as

Since azimuthal mode numbers can be positive or negative integers or zero,

reflecting the property of Bessel functions of integer orders: .

Using the known relations for Bessel functions

the following properties of the triple-Bessel integral can be derived as

The 3B integral is thus the Hankel transform of of order , and its inverse Hankel transform is also seen to exist.

Relations (12) and (13) yield

For , . If in addition ,

For , , the following result is known [10, p. 412] in terms of Legendre functions :

The following simpler case is also known, corresponding to the case , equation 3 in Watson’s treatise [10, p. 411]:

where

Relations (15)-(18) point to the existence of singular behavior of the 3B integral at specific values of its arguments. While (17) and (18) are formulas for the 3B integrals in the sense of Riemann integrals, (15) and (16) refer to them in the sense of distributions.

### Approximation to the Triple-Bessel Integral

This integral occurs as part of a kernel in an integral over the mode numbers , over the interval , which itself is a term in a summation over azimuthal mode numbers , (integers) from to . Although the exact evaluation of this integral in the sense of a Riemann integral is probably not possible [9], the existence of the singular behavior noted above suggests that its behavior in the neighborhood of singularities provides the dominant contribution. In view of known difficulties in the analytical approach [10] and the availability of the sophisticated numerical integrator in Mathematica 8.0.4, a semi-empirical approach [11] is chosen. The upper limit of integration is taken as a large but finite number , so that the required integral can be treated as the limit of a sequence for increasing .

For large values of argument, the Bessel function can be approximated as [10]:

The scaling property of the 3B integral discussed above suggests that the radial mode numbers can always be scaled so that approximation (19) becomes a good approximation in most of the infinite domain of integration (except perhaps in a small region around ). This yields

Each integral on the right-hand side of (21) is of the form

For ,

For ,

The functions and are the Fresnel cosine and sine integrals. Applying (24) to (21),

This can be simplified as follows: When is an odd integer, the cosine term in the ratio

is zero, so that the ratio is . Similarly, if this number is an even integer or zero, the sine term is zero, so that the ratio must be . Hence this ratio can be represented as . Applying this logic to (25),

Equation (25) has the advantage that both its sides can be independently computed, and their numerical equivalence and computational time can be compared, particularly at the singular points. Each term in (25) has a factor of the form

The functions , which are mirror images of each other, are plotted in Figure 1 and some numerical values are shown in Table 1.

Figure 1. Plot of functions .

Table 1. Numerical values of functions .

For a large but finite number , only one term out of four in (25) dominates when , and then the integral has the value

The full width at half maximum (FWHM) of the peaks of as a function of any one of the radial mode number arguments is seen to be , irrespective of other arguments. The semi-empirical relations (27) and (28), based on a large but finite upper limit of integration, are numerically validated in the next section, showing that (20) is a reasonably good approximation near the upper and lower singularities at and .

When and are equal, the last two terms in (25) become comparable, and (28) breaks down. However, the question of how close and can be before (28) breaks down can form the basis of a numerical experiment. Set , where is a sufficiently large positive number. Then (28) takes the form

The ratio of the numerically evaluated LHS and RHS for randomly selected parameters can be plotted against by varying over many decades. This is considered in the next section. The value of above which the formula breaks down is seen to be of the order of , irrespective of the upper limit of integration, which is varied over five decades. The numerical experiments show that (28) remains valid for increasingly close values of and as the limit of integration is increased.

The limit requires some care. It is seen from (27) and (28) that the maximum of at shifts to and tends to infinity. At the same time, the FWHM of the peak tends to zero. The weighted area of the peak at for a smooth weight function is defined in terms of an arbitrarily small positive number and a function where as

For smooth weight functions, in the limit as , and the integral on the RHS of (29) can be written as

Observing that the integral does not converge, choose , in order to ensure convergence in the infinite limit of . The integral on the RHS of (30) is then evaluated with infinite limits

The weighted area under the peak is therefore given by the limit of (31)

Similar logic can be applied to (28) when is treated as a fixed nonzero number

The region away from singularities does not need a finite upper limit of integration. For this case, a better approximation for the 3B integral can be constructed by applying (19) to only two out of three Bessel functions:

This can be computed by recognizing that (for positive )

The Heaviside -function in (36) and (37) excludes the singular points. Some general features of this approximate expression require attention. For all pairs of primed azimuthal mode numbers whose sum is even (or odd), the expression is unchanged. For , the third term in (35) varies as , and the fourth term is zero. For , the expression oscillates rapidly as and at the same time falls off in magnitude as . For , the last terms of (36) and (37) provide a sharp fall. Significant magnitude of the expression is thus expected only in the region satisfying the triangle inequality: as in the analytical results (17) and (18).

A computational survey of the behavior of the above approximation to the 3B integral can be graphically presented in terms of scaled parameters , , making use of the scaling relation mentioned above:

This may be conveniently displayed in terms of the scaled :

This parameter is negative below the lower singularity and greater than 1 above the upper singularity.

This is considered in the next section.

### Numerical Exploration of the Triple-Bessel Integral

Numerical exploration of the triple-Bessel integral, treated as a function of for given values of , , uses Mathematica 8.0.4, which implements sophisticated algorithms for numerical integration of highly oscillating functions with user-defined precision. Since the integral is expected to diverge at , the upper limit of numerical integration is kept at a high but finite number, which is also a parameter varied in the study.

Although the NIntegrate function in Mathematica is well tested, it is nevertheless prudent to benchmark it against the known result (18) in order to confirm that its options and parameters are correctly chosen for the intended application. The result of benchmark calculations is displayed in Figure 2, which summarizes the results of thousands of calculations with randomly generated arguments satisfying the triangle inequality and located near the lower singularity, upper singularity, or between the singularities. It is seen that there is quite good, but not perfect, agreement between the numerically calculated result and the formula for the following options: , , , and , which were chosen after a few trials.

Outside the triangle inequality region, NIntegrate gives very small values, consistent with the zero value from the analytical formula.

The upper limit of integration used in these calculations is symbolic Infinity, which is interpreted by NIntegrate as a large number that delivers the requested performance. The upper limit thus may not be the same in all calculations of the triple-Bessel integral. Using arbitrary fixed large numbers for the upper limit, however, gives absurd results.

Figure 2. Result of 3600 calculations of with randomly chosen and in various decade ranges and randomly chosen to lie in the region between and .

Figures 3 and 4 show the detail of the variation (magnified times) of the triple-Bessel integral, treated as a function of for , , , , in the neighborhood of and , with integration limit . The relevant Mathematica code is shown in Appendix 2. Figure 3 took over six hours to generate on a Core-2 Duo, 2 GHz, 4 GB ram machine. The resemblance of the shape of the plot to the functions is evident. Note the glitches at values of away from the singularity, showing that the numerical integrator does not work well for some values of the parameters. Since , the sign of the singularity in Figure 4 (time of computation: nearly nine hours) is negative for the given values of azimuthal mode numbers, as predicted by (29).

Figure 3. magnification detail in the neighborhood of : , , , , . Time required: about 22000 seconds.

Figure 4. magnification detail in the neighborhood of : , , , , . Time required: about 31000 seconds.

A randomized parameter space survey of (27) was carried out using the following protocol: 200 values of and were chosen randomly in each of four decades from to , was chosen randomly in the range to , and integers , , were chosen randomly in the range 0 to 100. For each set of parameters, the result of numerical integration over the limit , in the range of 9-12, divided by the value calculated using (27), was plotted against the latter. The result is summarized in Figure 5. A similar exercise for equation (28) is summarized in Figure 6. The code used is described in Appendix 3. It is seen that a higher integration limit gives less deviation from the formula.

Figure 5. Randomized parameter space survey of formula (27).

Figure 6. Randomized parameter space survey of formula (28).

As discussed in the previous section, (28) is expected to break down when and are sufficiently close, since the assumption concerning domination of only one term of (25) fails. This phenomenon is explored in the following test. (See Appendix 4).

The closeness of and is quantified in terms of a multiple of the parameter , which scans the width of the peak by setting . For a set of randomly selected parameters , , , , and , the ratio of values obtained from numerical integration and (28) is plotted against , which varies over 11 decades from to . This plot is generated for five values of upper integration limit , , , , and for the same set of randomly selected parameters, This is repeated for many randomly selected sets of , , in the range 0 to 100, 50 random values of in each of four decades from to , and in the range to The result is shown in Figure 7. The formula is seen to hold good when is above , irrespective of the limit of integration.

Figure 7. Randomized parameter space survey of formula (28) with .

The same dataset is displayed as a function of in Figure 8 for different values of the upper limit of integration.

Figure 8. Randomized parameter space survey of formula (28) with . The agreement between numerical integration and formula becomes better at closer values of as the limit of integration is increased.

Formula (28) is seen to be valid for increasingly smaller values of difference in and as the limit of integration is increased: less than a few parts per million when the integration limit is . For infinite limit of integration, the formula should therefore hold for arbitrarily close and .

The significance of numerical validation of (27) and (28) demonstrated above is that the dependence of the polarity of the peak on the three azimuthal mode numbers is correctly reproduced without a single error in several thousand calculations (where and differ more than a few parts per million), in addition to a reasonable agreement with peak amplitude and shape, for a large number of randomly selected parameters.

A similar exercise for the better approximation for regions away from singularities, formula (35), gives a very different picture: the formula overestimates the numerically calculated integral by several orders of magnitude. In order to understand this, a randomized survey of the approximation (35) and the numerical integration with upper limits , , was carried out and is summarized below. The program used in this survey is given in Appendix 5. The parameters used in the survey were: , , , . NIntegrate did not report any errors. Consistency of the trend in the survey suggests that the 3B integral tends to zero for the infinite upper limit of integration when : a (tentative) conclusion, which is predicated on the assumption that NIntegrate works correctly in these kinds of calculations. That can, in principle, be further examined by repeating this test for a much larger number of calculations and many different sets of options.

Using symbolic Infinity for the upper limit in these calculations is found to give absurd results. There is thus a distinct possibility that NIntegrate sometimes behaves in an anomalous manner while calculating the triple-Bessel integral.

Figure 9. Randomized survey of as functions of scaled , which are for various upper limits compared with (35).

The reasonable degree of agreement between the results of three theoretical formulas ((18), (27), and (28)) and numerical integration (with symbolic Infinity as the upper limit for (18) and a large number as the upper limit in (27) and (28)), with randomly selected parameters, tends to establish the credibility of numerical integration. At the same time, the fact that absurd results are produced by using a large number as the upper limit in the test for (18) and symbolic Infinity as the upper limit in the above test suggests that there is scope for improvement in the implementation of NIntegrate.

### Proposed Approximate Formula for the Triple-Bessel Integral

The above discussion and numerical experimentation lead to the following observations, subject to the caveat concerning the correct behavior of NIntegrate:

1. tends to infinity as .
2. The FWHM of the peak of at tends to zero as .
3. The weighted integrals (32) and (33) over this peak converge to finite values for smooth weight functions.
4. Both outside and within the triangle inequality region , is negligibly small.

These observations suggest the following approximate semi-empirical representation of the 3B integral in the sense of distribution with respect to for smooth weight functions:

Here, for , and for . Note the inclusion of the case : this is proposed as the generalization of the results displayed in Figure 8 to the infinite integration limit. Putting , , in (38) recovers the usual orthogonality formula (12) for Bessel functions, supporting the discussion concerning equation (33) and also the inclusion of the case.

Its only claim to merit is that it is a less restrictive approximation than the total neglect of irrotational and scalar fields in spatial-spectral treatments of nonlinear dynamics of continuous media. Adoption of (38) as a kind of tentative semi-empirical model for the triple-Bessel integral, at least until a better replacement becomes available, would enable construction of a theoretical framework for studying effects arising out of the nonlinear interaction of solenoidal and irrotational electron modes in the two-fluid plasma model in cylindrical geometry [6], which is currently an intractable problem. This framework can be improved if and when further research, hopefully motivated by this discussion, leads to a better representation of (38).

### Discussion

The insight gained above about the properties and numerical behavior of the 3B integral has significant implications concerning the usefulness of the Fourier space representation of nonlinear dynamics of continuous media in cylindrical geometry. The question, how far is infinity?, in the context of a physical theory, is seen to have an answer in terms of the error with which a physical quantity, such as radial mode number, might be defined using an experimental procedure. Two radial mode numbers may then be considered equivalent if their difference lies within the error band. The reciprocal of this error band would then represent the order of magnitude of the limit of integration in coordinate space, which reproduces, within experimental errors, the result that nonlinear interaction of these modes generates a near-zero radial-mode-number mode.

The origin of the singularities in this integral, treated as a function of three radial and three azimuthal mode numbers, can be seen from (20) to lie in resonant interaction between three waves of radial mode numbers , , , having phases governed by the azimuthal mode numbers , , , when any radial mode number equals the sum or difference of the other two radial mode numbers.

The triple-Bessel integral occurs in expressions that describe how the rate of change of the spectrum of a given quantity (the evolving mode), labeled by mode numbers without primes, depends on the interaction between the spectra of two other quantities, labeled by primed mode numbers (the interacting modes) and is accompanied by the factor . The semi-empirical model (38) then takes on a much simpler form, expressed below as a distribution with respect to using the invariance with respect to a permutation of pairs of arguments , , :

According to the model formula (40), the contribution of the 3B integral to the evolution of the spectrum is governed by two singularities. One evolves the spectrum toward higher radial mode numbers than the radial mode numbers of interacting quantities, generating smaller radial scale lengths and steeper gradients as in shock formation or turbulence. The other evolves the spectrum toward radial scale lengths larger than the radial scale lengths of the interacting modes, indicating the formation of a larger structure, as in self-organization. The tendencies toward turbulence/shock formation or self-organization are influenced by azimuthal mode numbers.

The computational advantage of the semi-empirical model of the 3B integral presented above can be judged by comparison with the approach of Chen, Shan, and Montgomery (CSM) [4, 5]. The first advantage over the method of CSM is that the spatial-spectral representation of irrotational vector fields as well as scalar fields becomes possible in addition to solenoidal fields: while the CSM method is limited to solenoidal fields only, over the infinite domain, an orthogonal basis becomes simultaneously available for scalar, irrotational, and solenoidal fields, facilitating the investigation of effects due to their nonlinear interaction [6]. The second advantage is that rather than computing a table of numerical values of coupling coefficients for an infinity (or a suitably large number) of azimuthal mode numbers, the dependence on azimuthal mode numbers is reduced to a simple computation of sign of a term. The semi-empirical formula (40) takes on the role of the table of coupling coefficients in the Galerkin approach of CSM [4, 5].

The spatial-spectral representation of the equation of continuity then becomes

The third advantage is the significant reduction in the computational complexity of equation (7) noticed in its reduced form (40): reduction in the number of integrations, simpler dependence on azimuthal mode numbers, and simpler form of the coupling factors. This opens up the possibility of doing computations similar to those of CSM without the restriction to incompressible flow and with a much larger number of modes.

### Summary and Conclusion

This article concerns the integral of the product of three Bessel functions of the first kind of integer orders over the semi-infinite domain , which is encountered when transforming a system of nonlinear partial differential equations containing solenoidal vector fields, scalar fields, and irrotational vector fields, using Chandrasekhar-Kendall functions, their generating function, and its gradient as orthogonal bases over the semi-infinite domain. Known properties of this integral, which show singular behavior, are summarized. Using the approximation of Bessel functions for large arguments, it is shown that the triple-Bessel integral represents resonant three-wave interaction that leads to singularities when one argument equals the sum or difference of the other two arguments. Numerical investigations using Mathematica have been used to validate a scaling formula for the polarity, absolute magnitude, and shape of the singular peak with a large finite upper limit of integration. A semi-empirical approximate formula of the integral has been suggested in terms of the Dirac delta function for regions near the singularities for smooth weight functions. It is shown that the computational complexity of a cylindrical Fourier-space representation of the equation of continuity is reduced considerably using this formulation. This opens up the possibility of computational efforts similar to those of CSM without restriction to incompressible dynamics for a much larger number of modes.

### Appendix 1

This generates Figure 2.

### Appendix 2

This generates Figures 3 and 4.

### Appendix 3

This generates Figure 5.

This generates Figure 6.

### Appendix 4

This generates Figures 7 and 8.

### Appendix 5

This generates Figure 9.

### Appendix 6

The following outlines an argument in support of the assumption that the 3B integral is zero outside the triangle inequality region. Using the integral representation of the Bessel function

the 3B integral can be written as

The integral would be zero if has no zeros. This would correspond to the region outside the triangle inequality.

### Acknowledgments

This article owes its origin to the following insightful question from an anonymous referee of Physics of Plasmas: “How far is infinity?” The help given by William Rummler of Wolfram Technical Support is gratefully acknowledged.

### References

 [1] S. Chandrasekhar and P. C. Kendall, “On Force-Free Magnetic Fields,” Astrophysical Journal, 126, 1957 pp. 457-460. doi:10.1086/146413. [2] Z. Yoshida, “Eigenfunction Expansions Associated with the Curl Derivatives in Cylindrical Geometries: Completeness of Chandrasekhar-Kendall Eigenfunctions,” Journal of Mathematical Physics, 33(4), 1992 pp. 1252-1256. doi:10.1063/1.529703. [3] D. Montgomery, L. Turner, and G. Vahala, “Three Dimensional Magnetohydrodynamic Turbulence in Cylindrical Geometry,” Physics of Fluids, 21(5), 1978 pp. 757-764. doi:10.1063/1.862295. [4] H. Chen, X. Shan, and D. Montgomery, “Galerkin Approximations for Dissipative Magnetohydrodynamics,” Physical Review A, 42(10), 1990 pp. 6158-6165. doi:10.1103/PhysRevA.42.6158. [5] X. Shan, D. Montgomery, and H. Chen, “Nonlinear Magnetohydrodynamics by Galerkin-Method Computation,” Physical Review A, 44(10), 1991 pp. 6800-6818. doi:10.1103/PhysRevA.44.6800. [6] S. K. H. Auluck, “Role of Electron-Inertia-Linked Current Source Terms in the Physics of Cylindrically Symmetric Imploding Snowplow Shocks,” Physics of Plasma, 9(11), 2002 pp. 4488-4494. doi:10.1063/1.1508775. [7] P. M. Morse and H. Feshbach, Methods of Theoretical Physics, Part I, New York: McGraw-Hill, 1953 p. 3. [8] S. K. H. Auluck, “Coherent Effects in the Stochastic Electrodynamics of Two-Fluid Plasma.” arxiv.org/abs/1208.1573. [9] G. N. Watson, “An Infinite Integral Involving Bessel Functions,” Journal of the London Mathematical Society, s1-9(1), 1934 pp. 16-22.jlms.oxfordjournals.org/content/s1-9/1/16.extract. [10] G. N. Watson, A Treatise on the Theory of Bessel Functions, 2nd ed., London: Cambridge University Press, 1944. [11] D. H. Bailey, J. M. Borwein, N. J. Calkin, R. Girgensohn, D. R. Luke, and V. H. Moll, Experimental Mathematics in Action, Wellesley, MA: A. K. Peters, 2007. S. K. H. Auluck, “On the Integral of the Product of Three Bessel Functions over an Infinite Domain,” The Mathematica Journal, 2012. dx.doi.org/doi:10.3888/tmj.14-15.