EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 15, No. 3, 2022, 856-863 ISSN 1307-5543 – ejpam.com Published by New York Business Global Double Integral involving the Product of the Bessel Function of the First Kind and modified Bessel Function of the Second Kind: Derivation and Evaluation Robert Reynolds1,∗, Allan Stauffer1 1 Department of Mathematics and Statistics, Faculty of Science, York University, Toronto, Ontario, Canada, M3J1P3 Abstract. A double integral whose kernel involves the Bessel functions Kv(xβ) and Jv(yα) is derived. This integral is expressed in terms of the Hurwitz-Lerch zeta function and evaluated for various values of the parameters involved. Some examples are evaluated and expressed in terms of fundamental constants. All the results in this work are new. 2020 Mathematics Subject Classifications: 30E20, 33-01, 33-03, 33-04, 33-33B Key Words and Phrases: Bessel functions, double integral, Cauchy integral 1. Introduction Integrals involving Bessel functions have been studied in the works by Glasser [3], where the study of wave propagation along a coaxial cable was investigated, Temme [7], where the mathematical discussion of the exchange processes, of heat or of matter (as in ion exchange or adsorption), that arise when a fluid flows through the pores or voids along a column containing matter in the solid state, was studied. Throughout these works the authors derive definite integrals involving the Bessel function or the product of Bessel functions for specific orders. in our present paper we will be expanding on the previous for- mulae by deriving a double integral of the product of Bessel functions over a general order. In this paper we derive the double definite integral given by (1) ∫ ∞ 0 ∫ ∞ 0 xm−1y1−mKv(xβ)Jv(yα) log k ( ax y ) dxdy where the parameters k, a, α, β, v,m are general complex numbers and Re(α, β, v,m) > 0, Re(v) < Re(m) < 3/2. This definite integral will be used to derive special cases in ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v15i3.4239 Email addresses: milver@my.yorku.ca (R. Reynolds), stauffer@yorku.ca (A. Stauffer) https://www.ejpam.com 856 © 2022 EJPAM All rights reserved. R. Reynolds, A. Stauffer / Eur. J. Pure Appl. Math, 15 (3) (2022), 856-863 857 terms of special functions and fundamental constants. The derivations follow the method used by us in [6]. This method involves using a form of the generalized Cauchy’s integral formula given by yk Γ(k + 1) = 1 2πi ∫ C ewy wk+1 dw. (2) where C is in general an open contour in the complex plane where the bilinear concomitant has the same value at the end points of the contour. We then multiply both sides by a function of x and y, then take a definite double integral of both sides. This yields a definite integral in terms of a contour integral. Then we multiply both sides of Equation (2) by another function of x and y and take the infinite sums of both sides such that the contour integral of both equations are the same. 2. Definite Integral of the Contour Integral We use the method in [6]. The variable of integration in the contour integral is t = w + m. The cut and contour are in the first quadrant of the complex t-plane. The cut approaches the origin from the interior of the first quadrant and the contour goes round the origin with zero radius and is on opposite sides of the cut. Using a generalization of Cauchy’s integral formula we form the double integral by replacing y by log ( ax y ) and multiplying by xm−1y1−mKv(xβ)Jv(yα) then taking the definite integral with respect to x ∈ [0,∞) and y ∈ [0,∞) to obtain (3) 1 Γ(k + 1) ∫ ∞ 0 ∫ ∞ 0 xm−1y1−mKv(xβ)Jv(yα) log k ( ax y ) dxdy = 1 2πi ∫ ∞ 0 ∫ ∞ 0 ∫ C aww−k−1xm+w−1y−m−w+1Kv(xβ)Jv(yα)dwdxdy = 1 2πi ∫ C ∫ ∞ 0 ∫ ∞ 0 aww−k−1xm+w−1y−m−w+1Kv(xβ)Jv(yα)dxdydw = 1 2πi ∫ C 1 2 πaww−k−1αm+w−2β−m−w csc ( 1 2 π(m− v + w) ) dw from equations (3.10.1.2) and (3.14.3) in [1] where Re(α) > 0, |Re(v)|< Re(w +m− v) < 3/2 and using the reflection formula (8.334.3) in [4] for the Gamma function. We are able to switch the order of integration over x and y using Fubini’s theorem since the integrand is of bounded measure over the space C× [0,∞)× [0,∞) 3. The Hurwitz-Lerch zeta Function and Infinite Sum of the Contour Integral In this section we use Equation (2) to derive the contour integral representations for the Hurwitz-Lerch zeta function. R. Reynolds, A. Stauffer / Eur. J. Pure Appl. Math, 15 (3) (2022), 856-863 858 3.1. The Hurwitz-Lerch zeta Function The Hurwitz-Lerch zeta function (25.14) in [2] has a series representation given by Φ(z, s, v) = ∞∑ n=0 (v + n)−szn (4) where |z|< 1, v ̸= 0,−1, .. and is continued analytically by its integral representation given by Φ(z, s, v) = 1 Γ(s) ∫ ∞ 0 ts−1e−vt 1− ze−t dt = 1 Γ(s) ∫ ∞ 0 ts−1e−(v−1)t et − z dt (5) where Re(v) > 0, and either |z|≤ 1, z ̸= 1, Re(s) > 0, or z = 1, Re(s) > 1. 3.2. Infinite sum of the Contour Integral Using equation (2) and replacing y by log(a) + log(α) − log(β) + 1 2 iπ(2y + 1) then multiplying both sides by −iπαm−2β−me 1 2 iπ(2y+1)(m−v) taking the infinite sum over y ∈ [0,∞) and simplifying in terms of the Hurwitz-Lerch zeta function we obtain (6) − 1 Γ(k + 1) iπk+1αm−2β−me 1 2 iπ(k+m−v) Φ ( eiπ(m−v),−k, −2i log(a)− 2i log(α) + 2i log(β) + π 2π ) = − 1 2πi ∞∑ y=0 ∫ C iπw−k−1αm−2β−m exp ( w(log(a) + log(α)− log(β)) + 1 2 iπ(2y + 1)(m− v + w) ) dw = − 1 2πi ∫ C ∞∑ y=0 iπw−k−1αm−2β−m exp ( w(log(a) + log(α)− log(β)) + 1 2 iπ(2y + 1)(m− v + w) ) dw = 1 2πi ∫ C 1 2 πaww−k−1αm+w−2β−m−w csc ( 1 2 π(m− v + w) ) dw from equation (1.232.2) in [4] where Im ( 1 2π(m− v + w) ) > 0 in order for the sum to converge. R. Reynolds, A. Stauffer / Eur. J. Pure Appl. Math, 15 (3) (2022), 856-863 859 4. Definite Integral in terms of the Hurwitz-Lerch zeta Function Theorem 1. For all k, a ∈ C, Re(α, β, v,m) > 0, Re(v) < Re(m) < 3/2, (7 ) ∫ ∞ 0 ∫ ∞ 0 xm−1y1−mKv(xβ)Jv(yα) log k ( ax y ) dxdy = −iπk+1αm−2β−me 1 2 iπ(k+m−v) Φ ( eiπ(m−v),−k, −2i log(a)− 2i log(α) + 2i log(β) + π 2π ) Proof. The right-hand sides of relations (3) and (6) are identical; hence, the left-hand sides of the same are identical too. Simplifying with the Gamma function yields the desired conclusion. Example 1. The degenerate case.∫ ∞ 0 ∫ ∞ 0 xm−1y1−mKv(xβ)Jv(yα)dxdy = 1 2 παm−2β−m csc ( 1 2 π(m− v) ) (8) Proof. Use equation (7) and set k = 0 and simplify using entry (2) in Table below (64:12:7) in [5]. Example 2. The Hurwitz zeta function ζ(s, v) (9 ) ∫ ∞ 0 ∫ ∞ 0 √ xeβ(−x) sin(αy) logk ( ax y ) √ y √ βx √ αy dxdy = − 1 √ αβ3/2 ie 1 2 iπ(k+1)πk+1 ( 2kζ ( −k, −2i log(a)− 2i log(α) + 2i log(β) + π 4π ) − 2kζ ( −k, 1 2 ( −2i log(a)− 2i log(α) + 2i log(β) + π 2π + 1 ))) Proof. Use equation (7) and set m = 3/2, v = 1/2 and simplify using entry (4) in Table below (64:12:7) in [5]. Example 3. (10 ) ∫ ∞ 0 ∫ ∞ 0 e−x sin(y) log ( x y ) y ( log2 ( x y ) + π2 )dxdy = 0 and (11) ∫ ∞ 0 ∫ ∞ 0 e−x sin(y) y ( log2 ( x y ) + π2 )dxdy = 2 π − 1 2 R. Reynolds, A. Stauffer / Eur. J. Pure Appl. Math, 15 (3) (2022), 856-863 860 Proof. Use equation (9) apply l’Hopital’s rule as k → −1 and set a = −1, α = β = 1 rationalize the denominator and simplify using entry (2) in Table below (64:7) in [5]. Example 4. ∫ ∞ 0 ∫ ∞ 0 √ xK 1 4 (x)J 1 4 (y) log ( x y ) √ y ( log2 ( x y ) + π2 ) dxdy = 1 4 (√ 2π − 8 sin (π 8 ) − 2 √ 2 tanh−1 ( sin (π 8 ))) (12) and (13 ) ∫ ∞ 0 ∫ ∞ 0 √ xK 1 4 (x)J 1 4 (y) √ y ( log2 ( x y ) + π2 )dxdy = 8 cos ( π 8 ) − √ 2 ( π + 2 tanh−1 ( sin ( π 8 ))) 4π Proof. Use equation (7) and set k = −1,m = 3/2, v = 1/4, α = β = 1, a = −1 and simplify using entry (3) in Table below (64:12:7) in [5]. Example 5. (14 ) ∫ ∞ 0 ∫ ∞ 0 xm−1y1−mKv(xβ)Jv(yα) log ( −βx αy ) dxdy = 2iαm−2β−me− 1 2 iπ(2m−2v+1) ( e 1 2 iπ(m−v) − tanh−1 ( e 1 2 iπ(m−v) )) Proof. Use equation (7) and set k = −1,m = 3/2, v = 1/4, a = β/α and simplify using entry (3) in Table below (64:12:7) in [5]. Example 6. The Polylogarithm function Lik(z),∫ ∞ 0 ∫ ∞ 0 xm−1y1−mKv(x)Jv(y) log k ( ix y ) dxdy = −iπk+1e 1 2 iπ(k−m+v)Li−k ( eiπ(m−v) ) (15) Proof. Use equation (7) and set a = i, β = α = 1 and simplify using equation (64:12:2) in [5]. Example 7. The Polylogarithm function Li2(z), (16 ) ∫ ∞ 0 ∫ ∞ 0 xm−1y1−mKv(x)Jv(y) log2 ( ix y ) dxdy = − ie− 1 2 iπ(m−v+2)Li2 ( eiπ(m−v) ) π Proof. Use equation (7) and set k = −2, a = i, β = α = 1 and simplify using equation (64:12:2) in [5]. R. Reynolds, A. Stauffer / Eur. J. Pure Appl. Math, 15 (3) (2022), 856-863 861 Example 8. Catalan’s constant G, (17 ) ∫ ∞ 0 ∫ ∞ 0 e−x√x sin(y) ( π2 − 4 log2 ( x y )) y3/2 ( 4 log2 ( x y ) + π2 )2 dxdy = 48G+ π2 192 √ 2π and (18 ) ∫ ∞ 0 ∫ ∞ 0 e−x√x sin(y) log ( x y ) y3/2 ( 4 log2 ( x y ) + π2 )2dxdy = −π2 − 48G 768 √ 2π2 Proof. Use equation (16) and set m = 2, v = 1/2 and simplify. Example 9.∫ ∞ 0 ∫ ∞ 0 y−m−p+1Kv(x)Jv(y) (y mxp − xmyp) x log ( x y ) dxdy = 2 ( tanh−1 ( e 1 2 iπ(m−v) ) − tanh−1 ( e 1 2 iπ(p−v) )) (19) Proof. Use equation (7) and form a second equation by replacing m → p and taking their difference. Next set k = −1, a = 1, α = β = 1 and simplify using entry (3) in Table below (64:12:7) in [5]. Example 10.∫ ∞ 0 ∫ ∞ 0 ( x− √ x √ y ) K 1 3 (x)J 1 3 (y) y log ( x y ) dxdy = 2 tanh−1 ( 1− √ 3 2 + 1√ 2 ) (20) Proof. Use equation (19) set v = 1/3,m = 3/2, p = 2 and simplify. Example 11. (21 ) ∫ ∞ 0 ∫ ∞ 0 x2/5 ( 10 √ y − 10 √ x ) K0(x)J0(y) √ y log ( x y ) dxdy = − tanh−1 ( 1 29 √ 2 ( 254− 31 √ 5− 2 √ 4505 + 1109 √ 5 )) Proof. Use equation (19) set v = 0,m = 3/2, p = 7/5 and simplify. R. Reynolds, A. Stauffer / Eur. J. Pure Appl. Math, 15 (3) (2022), 856-863 862 Example 12.∫ ∞ 0 ∫ ∞ 0 (√ x− x2/5 10 √ y ) K 1 5 (x)J 1 5 (y) √ y log ( x y ) dxdy = 1 2 ( log ( 5− 2 √ 5 ) + 4 tanh−1 ( tan ( 3π 40 ))) (22) Proof. Use equation (19) set v = 1/5,m = 3/2, p = 7/5 and simplify. Example 13. (23 ) ∫ ∞ 0 ∫ ∞ 0 5 √ x ( y3/10 − x3/10 ) K 2 9 (x)J 2 9 (y) √ y log ( x y ) dxdy = − tanh−1 ( sin ( π 90 )) − 2 tanh−1 ( tan ( 5π 72 )) Proof. Use equation (19) set v = 2/9,m = 3/2, p = 6/5 and simplify. Example 14. (24 ) ∫ ∞ 0 ∫ ∞ 0 ( x− √ x √ y ) K 1 3 ( 4 √ 3x ) J 1 3 ( 2 √ 2y ) y log ( x y ) dxdy = 1 48 ( (−1)7/12 4 √ 6Φ ( − 6 √ −1, 1, π + i log(6) 2π ) − (−1)5/6Φ ( −(−1)2/3, 1, π + i log(6) 2π )) Proof. Use equation (7) and form a second equation by replacing m → p and taking their difference. Next set k = −1, a = 1, α = √ 2, β = √ 3,m = 3/2, p = 2, v = 1/7 and simplify. Example 15. (25 ) ∫ ∞ 0 ∫ ∞ 0 √ xK 1 3 (3x)J 1 3 (2y) √ y log ( −2x y ) dxdy = − (−1)7/12Φ ( − 6 √ −1, 1, 3π−4i log(2)+2i log(3) 2π ) 3 √ 6 Proof. Use equation (7) set k = −1, a = −2,m = 3/2, v = 1/3, , α = 2, β = 3 and simplify. REFERENCES 863 Example 16. (26 ) ∫ ∞ 0 ∫ ∞ 0 √ xK 1√ 5 ( 5x√ 11 ) J 1√ 5 ( y√ 7 ) √ y √ log ( −3x y ) dxdy = 1 5 4 √ 7113/4e − iπ 2 √ 5 √ π 5 Φ ( −ie − iπ√ 5 , 1 2 , 3 2 + i log ( 175 99 ) 2π ) Proof. Use equation (7) set k = −1/2, a = −3,m = 3/2, v = 1/ √ 5, , α = 1/ √ 7, β = 5/ √ 11 and simplify. 5. Discussion In this paper, we have presented a novel method for deriving a new double integral involving the product of Bessel functions along with some interesting definite integrals using contour integration. The results presented were numerically verified for both real and imaginary and complex values of the parameters in the integrals using Mathematica by Wolfram. References [1] Yu. A. Brychkov, O. I. Marichev, and N. V. Savischenko. Handbook of Mellin tranforms. CRC Press., 2019. [2] Nist digital library of mathematical functions. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds. [3] M. L. Glasser. Integral representations for the exceptional univariate lommel functions. J. Phys. A, 43, 2010. [4] I. S. Gradshteyn and I. M. Ryzhik. Table of integrals, series, and products. Else- vier/Academic Press, Amsterdam, seventh edition, 2007. [5] Keith B. Oldham, Jan Myland, and Jerome Spanier. An Atlas of Functions: with Equator, the Atlas Function Calculator. Springer Science & Business Media, 07 2010. [6] Robert Reynolds and Allan Stauffer. A method for evaluating definite integrals in terms of special functions with examples. International Mathematical Forum, 15:235– 244, 2020. [7] N.M. Temme. A double integral containing the modified bessel function: Asymptotics and computation. MATHEMATICS OF COMPIJTATION, 47:683–691, 1986.