6_816_paris.dvi EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 3, No. 6, 2010, 1006-1031 ISSN 1307-5543 – www.ejpam.com SPECIAL ISSUE ON COMPLEX ANALYSIS: THEORY AND APPLICATIONS DEDICATED TO PROFESSOR HARI M. SRIVASTAVA, ON THE OCCASION OF HIS 70TH BIRTHDAY Asymptotic Expansion of n-dimensional Faxén-type Integrals Richard B. Paris University of Abertay Dundee, Dundee DD1 1HG, UK Abstract. The asymptotic expansion of n-dimensional extensions of Faxén’s integral In(z) are derived for large complex values of the variable z. The theory relies on the asymptotics of the generalised hy- pergeometric, or Wright, function. The coefficients in the exponential expansion are obtained by means of an algorithm applicable for arbitrary n. Numerical examples are given to illustrate the accuracy of the expansions. 2000 Mathematics Subject Classifications: 30E15, 33C70, 41A60, 41A63 Key Words and Phrases: Faxén’s integral, Asymptotic expansion, Wright function, generalised hyper- geometric functions 1. Introduction We consider the n-dimensional integral In(z) = λn ∫ ∞ 0 . . . ∫ ∞ 0 x ν1−1 1 . . . xνn−1 n e− f (x1 ,...,xn;z) d x1 . . . d xn, (1) where f (x1, . . . , xn; z) = n ∑ j=1 x µ j j − zx m1 1 . . . xmn n , λn = n ∏ j=1 µ j, (2) and the factor λn has been added for later convenience. We suppose that the exponents (not necessarily integers) satisfy µ j > m j > 0, Re (ν j)> 0 (1≤ j ≤ n) and that z denotes a complex Email address: r.paris�abertay.a .uk http://www.ejpam.com 1006 c© 2010 EJPAM All rights reserved. R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1007 variable. For convergence of In(z) we require that µ j and m j be further restricted so that the parameter κ, defined by κ= 1− n ∑ j=1 m j µ j , (3) should satisfy 0< κ < 1. The geometrical interpretation of this condition results from consid- eration of the Newton diagram associated with the phase function f . In the two-dimensional case n = 2, the Newton diagram is given by the boundary of the convex hull formed by the point (m1, m2) and the points (µ1, 0) and (0,µ2) situated on the coordinate axes, with the line joining these last two points being termed the back face. Extension to n ≥ 3 dimensions is straightforward with the back face being a hyperplane in n dimensions passing through the points µ j (1 ≤ j ≤ n) on the coordinate axes. The condition 0 < κ < 1 then corresponds to the internal point (m1, . . . , mn) being situated in front of the back face of the Newton diagram. In the case n= 1, we have (dropping the subscript 1 on the parameters) I1(z) = µ ∫ ∞ 0 xν−1e−xµ+zxm d x = ∫ ∞ 0 τ(ν/µ)−1e−τ+zτm/µ dτ, where κ= 1−(m/µ)< 1. This integral can be expressed as Fi(m/µ,ν/µ; z), where Fi denotes Faxén’s integral defined by [9, p. 332] Fi(a, b; z) = ∫ ∞ 0 τb−1e−τ+zτa dτ (0≤ Re(a)< 1, Re(b)> 0). Consequently, (1) can be considered as an n-dimensional extension of Faxén’s integral. A dif- ferent extension of Faxén’s integral as a one-dimensional integral with more than one internal point in the phase function f has been considered in [7]. Another integral that is related to (1), but with different domains of integration, is Jn(z) = λn ∫ ∞ −∞ . . . ∫ ∞ −∞ x ν1−1 1 . . . xνn−1 n e− f (x1,...,xn;z) d x1 . . . d xn, (4) in which the µ j are now all restricted to be positive even integers. When the exponents ν j are nonintegers, the integral Jn(z) is specified by taking the integration paths along the negative x j-axes to be along the upper side of the branch cuts on these axes. Variants of the integral Jn(z) can also be considered in which p < n of the integrals in (4) are evaluated over the interval (−∞,∞), with the remainder over the interval [0,∞). Special cases of the integral In(z) when n = 1 and ν = 1 were first studied asymptoti- cally in [4, 5], and more generally in [1], using the method of steepest descents. In [13] the asymptotic expansion of In(z) for large complex z was obtained by application of the asymptotic theory of the generalised hypergeometric, or Wright, function pΨq(z) defined in (5). The expansion was found to consist of an exponential expansion, which is dominant in the sector |arg z| < 1 2 πκ, together with an algebraic expansion dominant in the rest of the z- plane. An application of the integral In(z) in the particular case n= 2 has been given in [8] in R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1008 the discussion of two-dimensional Laplace integrals with more general phase and amplitude functions. More recently, Breen and Wood [3] have discussed an application of In(z) as a representa- tion of the solutions of certain high-order linear differential equations. One of these equations has the form y(n)(z)− p ∑ r=0 αrzr y(r)(z) = 0 (n> p > 0), where the αr are arbitrary coefficients. This equation has a basis of solutions given by y(z; s) = ∫ ∞ 0 . . . ∫ ∞ 0 x ν1−1 1 . . . x νp−1 p exp{szx1 . . . xp − (x n 1 + . . .+ xn p)/n} d x1 . . . d xp, where sn = νp and the exponents νr are related to the coefficients αr in a manner that we do not specify here. This integral is clearly related to Ip(sz) in (1) with the parameters µ j = n and m j = 1 (1≤ j ≤ p). The integral representation of solutions of the above differential equation when there are two lower-order derivatives (p = 2) was first given by Spitzer [18], with the general case of p < n lower-order derivatives being considered in [16]. These results are described in [14, pp. 130–133]. In [3] the asymptotics of the solutions y(z; s) were obtained using the theory developed in [13]. In this paper, we review the asymptotic expansion of the integral In(z) in (1) using the asymptotic theory of the Wright function. A recent account of the asymptotic theory of the lat- ter function has been presented in [11] and a discussion of the properties of 0Ψ1(z) (the gen- eralised Bessel function), together with its application to the solution of fractional diffusion- wave equations, can be found in [6]. It is shown how the expansion of In(z)may be employed to determine the asymptotic structure of the integral Jn(z) and its variants when some of the integrals in (4) are taken over [0,∞). 2. The Expansion of the Wright Function pΨq(z) for |z| →∞ The asymptotic expansion of the integrals In(z) and Jn(z) will be obtained by utilising the asymptotic theory of the Wright (or generalised hypergeometric) function which we present in this section. The Wright function pΨq(z) is defined by pΨq(z) ≡ pΨq � (α1, a1), . . . , (αp, ap) (β1, b1), . . . , (βq, bq) ; z � = ∞ ∑ k=0 g(k) zk k! , g(k) := ∏p r=1 Γ(αr k+ ar) ∏q r=1 Γ(βr k+ br) , (5) where p and q are nonnegative integers, the parameters αr and βr are real and positive and ar and br are arbitrary complex numbers. In addition, it is assumed that the αr and ar are subject to the restriction αr k+ ar 6= 0,−1,−2, . . . (k = 0,1,2, . . . ; 1≤ r ≤ p) (6) so that no gamma function in the numerator of (5) is singular. In the special case αr = βr = 1, the function pΨq(z) reduces to a multiple of the generalised hypergeometric function pFq((ap); (bq); z); see, for example, [17, p. 40]. R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1009 We summarise the asymptotic expansion of the Wright function pΨq(z) for |z| →∞ given in Wright [20, 21] and Braaksma [2]; for a summary, see also [12, §2.3] and [11]. We first introduce the parameters associated with g(k) given by κ= 1+ q ∑ r=1 βr − p ∑ r=1 αr , h= p ∏ r=1 ααr r q ∏ r=1 β−βr r , ϑ = p ∑ r=1 ar − q ∑ r=1 br + 1 2 (q− p), ϑ′ = 1− ϑ, (7) where, as usual, an empty product has unit value. If it is supposed that αr and βr are such that κ > 0, then pΨq(z) is uniformly and absolutely convergent for all finite z. It is clear that pΨq(z) is an entire function of z in this case. If κ = 0, the sum in (1) has a finite radius of convergence equal to h−1, whereas for κ < 0 the sum is divergent for all nonzero values of z. The parameter κ will be found to play a critical role in the asymptotic theory of pΨq(z) by determining the sectors in the z-plane in which its behaviour is either exponentially large, algebraic or exponentially small in character as |z| →∞. The exponential expansion Ep,q(z) is given by the formal asymptotic sum Ep,q(z) = ZϑeZ ∞ ∑ j=0 A j Z − j, Z = κ(hz)1/κ, (8) where the coefficients A j are those appearing in the inverse factorial expansion of g(s)/s! in the form g(s) Γ(s+ 1) = κ(hκκ)s �M−1 ∑ j=0 A j Γ(κs+ ϑ′ + j) + O(1) Γ(κs+ ϑ′ +M) � (9) for |s| → ∞ uniformly in |arg s| ≤ π− ε, ε > 0 and arbitrary positive integer M . The leading coefficient A0 is specified by A0 = (2π) 1 2 (p−q)κ− 1 2 −ϑ p ∏ r=1 α ar− 1 2 r q ∏ r=1 β 1 2 −br r . (10) The coefficients A j are independent of s and depend only on the parameters p, q, αr , βr , ar and br . An algorithm for their evaluation in specific cases when αr > 0, βr > 0 is described in Appendix A. The algebraic expansion Hp,q(z) follows from the Mellin-Barnes integral representation [12, §2.3] pΨq(z) = 1 2πi ∫ ∞i −∞i Γ(s)g(−s)(ze∓πi)−sds, |arg(−z)|< 1 2 π(2− κ), (11) where the upper or lower sign is chosen according as arg z > 0 or arg z < 0, respectively. The path of integration in (11) is indented near s = 0 to separate∗ the poles of Γ(s) situated at ∗This is always possible when the condition (6) is satisfied. R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1010 s = 0,−1,−2, . . . from those of g(−s) at sk,r = (ar + k)/αr , k = 0,1,2, . . . (1≤ r ≤ p). (12) In general there will be p such sequences of simple poles though, depending on the values of αr and ar , some of these poles could be multiple poles or even ordinary points if any of the Γ(βrs + br) are singular there. Displacement of the integration contour in (11) to the right over the poles of g(−s) followed by evaluation of the residues then generates the algebraic expansion of pΨq(z) valid as |z| →∞ in the sector in (11). If it is assumed that the parameters are such that the poles in (12) are all simple, we obtain the algebraic expansion given by Hp,q(ze∓πi), where Hp,q(z) = p ∑ j=1 α−1 j z−a j/α j Sp,q(z; j) (13) and Sp,q(z; j) denotes the formal asymptotic sum Sp,q(z; j) = ∞ ∑ k=0 (−)k k! Γ � k+ a j α j � ∏′ p r=1 Γ(ar −αrsk, j) ∏q r=1 Γ(br − βrsk, j) z−k/α j , (14) with the prime indicating the omission of the term corresponding to r = j in the product. This expression consists of p expansions with the leading behaviour z−a j/α j (1≤ j ≤ p). When the parameters αr and ar are such that some of the poles are of higher order, the expansion (13) is invalid and the residues must then be evaluated according to the multiplicity of the poles concerned; this will lead to terms involving log z in the algebraic expansion. We present the asymptotic expansion of pΨq(z) for large |z| only in the case when 0 < κ ≤ 2, since the value of this parameter associated with the integrals (1) and (4) must satisfy 0< κ < 1. A fuller list of expansion theorems is given in [2, 20, 21]; see also [11]. We have the following theorems, where throughout we let ε denote an arbitrarily small positive quantity. Theorem 1. If 0< κ < 2, then pΨq(z) ∼ ¨ Ep,q(z) +Hp,q(ze∓πi) in |arg z| ≤ 1 2 πκ Hp,q(ze∓πi) in |arg(−z)| ≤ 1 2 π(2− κ)− ε (15) as |z| → ∞. The upper or lower sign in Hp,q(ze∓πi) is chosen according as z lies in the upper or lower half-plane, respectively. It is seen that the z-plane is divided into two sectors, with a common vertex at z = 0, by the rays (the anti-Stokes lines) arg z = ±1 2 πκ. In the sector |arg z| < 1 2 πκ, the asymptotic character of pΨq(z) is exponentially large, whereas in the complementary sector |arg(−z)| < 1 2 π(2− κ), pΨq(z) is algebraic in character. The choice of signs in Hp,q(ze∓πi) results from the fact that the positive real axis arg z = 0 is a Stokes line, where the algebraic R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1011 expansion is maximally subdominant. Since pΨq(z) is an entire function of z, we may write pΨq(z) = pΨq(ze−2πi). Then, when π≤ arg z < 2π the algebraic expansion is (with the lower sign) Hp,q(ze−2πi eπi) = Hp,q(ze−πi), and so has the same form as when 0< arg z ≤ π. Hence the algebraic expansion associated with pΨq(z) can be written alternatively as Hp,q(ze−πi) in ε≤ arg z ≤ 2π− ε. (16) The above theorem does not take into account the presence of an exponentially small contribution beyond the sector |arg z| ≤ 1 2 πκ. This is covered by the more precise result in the following theorem [2, p. 331], [20, 21]. Theorem 2. If 2 3 ≤ κ≤ 2, then pΨq(z)∼ Ep,q(z) + Ep,q(ze∓2πi) +Hp,q(ze∓πi) (|arg z| ≤ π) (17) as |z| →∞. When 0< κ < 2 3 , we have pΨq(z)∼ ¨ Ep,q(z) +Hp,q(ze∓πi) in |arg z| ≤ 3 2 πκ− ε Hp,q(ze∓πi) in 3 2 πκ+ ε≤ |arg z| ≤ π (18) as |z| →∞. The upper or lower signs are chosen according as arg z > 0 or arg z < 0, respectively. Since, when 2 3 ≤ κ ≤ 2, Ep,q(z) is exponentially small in 1 2 πκ < |arg z| ≤ π then in the sense of Poincaré the expansion Ep,q(z) can be neglected. Similarly, E(ze−2πi) is exponentially small compared to Ep,q(z) in 0 ≤ arg z < π and consequently there is no inconsistency be- tween (17) and the second expansion in (15). However, in the neighbourhood of arg z = π, Ep,q(z) and Ep,q(ze∓2πi) are of comparable magnitude and, for real parameters, they combine to generate a real result on arg z = π. A similar remark applies to the expansion E(ze2πi) in −π < arg z ≤ 0. When κ < 2 3 , Ep,q(z) is exponentially small in the sectors 1 2 πκ < |arg z| < 3 2 πκ and the behaviour of pΨq(z) in the complementary sector 3 2 πκ < |arg z| ≤ π is then algebraic. An even more precise result can be given by recognising that the rays arg z = ±πκ are also Stokes lines, where Ep,q(z) is maximally subdominant with respect to Hp,q(ze∓πi). This will result in the expansion Ep,q(z) switching off (as |arg z| increases) across the Stokes lines arg z = ±πκ. Thus, when 0< κ < 1, (17) and (18) can be replaced by pΨq(z) ∼ ¨ Ep,q(z) +Hp,q(ze∓πi) in |arg z| ≤ πκ− ε Hp,q(ze∓πi) in πκ+ ε ≤ |arg z| ≤ π (19) as |z| → ∞. In Appendix B we present a numerical example for 2Ψ0(z) which demonstrates the truth of this assertion; a fuller discussion is given in [10]. Although the expansions in (15) and (18) are valid asymptotic descriptions, more accurate evaluation will result from using (19) which takes into account the Stokes phenomenon.† In the application to the integrals In(z) and Jn(z) we shall employ the expansion of pΨq(z) in the form given in (19). †The expansion in the neighbourhood of the Stokes lines arg z = 0 and arg z = ±πκ would necessitate a detailed treatment of the Stokes phenomenon that we do not consider here. R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1012 3. The Expansion of In(z) for |z| →∞ The series representation of the integral In(z) can be obtained by use of the Maclaurin expansion of the factor exp{zx m1 1 . . . x mn n } in (1) followed by termwise integration to yield In(z) = λn ∞ ∑ k=0 zk k! n ∏ r=1 ∫ ∞ 0 xνr+mr k−1 r exp{−xµr r } d xr = ∞ ∑ k=0 n ∏ r=1 Γ � νr +mr k µr � zk k! . (20) Comparison with (5) shows that the above series is a particular case of the Wright (or gener- alised hypergeometric) function given by In(z) = nΨ0 � (α1, a1), . . . , (αn, an) ; z � ≡ nΨ0(z), (21) where the parameters αr = mr µr , ar = νr µr (1≤ r ≤ n) (22) and the dash denotes the omission of a parameter sequence. The asymptotic expansion of In(z) for |z| → ∞ then follows from (19). With κ defined in (3) (which follows from the definition in (7)) we therefore have In(z)∼ ¨ En,0(z) +Hn,0(ze∓πi) in |arg z| ≤ πκ− ε Hn,0(ze∓πi) in πκ+ ε≤ |arg z| ≤ π, (23) where the exponential expansion En,0(z) and the algebraic expansion Hn,0(ze∓πi) are obtained from (8), (13) and (14) with p = n, q = 0, and the upper or lower signs are chosen according as arg z > 0 or arg z < 0, respectively. The leading coefficient A0 in En,0(z) is, from (10), given by A0 = (2π) n/2κ− 1 2 −ϑ n ∏ r=1 � mr µr �(νr/µr )− 1 2 . (24) The large |z| behaviour of In(z) is consequently exponentially large in the sector |arg z| < 1 2 πκ. Outside of this sector the behaviour is dominated by an algebraic expansion, with a subdominant exponentially small contribution being present in the sectors 1 2 πκ < |arg z| < πκ. We now give an example of the expansion of In(z) when n= 3; other numerical examples can be found in [13]. Consider the three-dimensional integral I3(z) = 36 ∫ ∞ 0 ∫ ∞ 0 ∫ ∞ 0 (x3/x2) 1 2 exp{−(x3 1 + x3 2 + x4 3 − z(x1 x2) 1 2 x3)} d x1d x2d x3, R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1013 which is associated with the parameters µ1 = µ2 = 3, µ3 = 4, m1 = m2 = 1 2 , m3 = 1 and ν1 = 1, ν2 = 1 2 , ν3 = 3 2 . From (21) we therefore have I3(z) = 3Ψ0 � (1 6 , 1 3 ), (1 6 , 1 6 ), (1 4 , 3 8 ) ; z � ≡ 3Ψ0(z). (25) From (3), (7) and (24) we obtain the parameters κ = 5 12 , h= 2−5/63−1/3, ϑ = −5 8 , A0 = 4π3/235/851/8. Then, from (8), E3,0(z) = Z−5/8eZ ∞ ∑ j=0 A j Z − j, Z = 5 48 (1 3 z3)4/5, (26) where, by use of the algorithm described in Appendix A, the first few normalised coefficients c j ≡ A j/A0 are found to be c1 = 67 144 , c2 = 23785 41472 , c3 = 106119923 89579520 , c4 = 181613304677 51597803520 , c5 = 102937183723339 7430083706880 , . . . . The poles in (12) are situated at sk,1 = 2+ 6k, sk,2 = 1+ 6k, sk,3 = 3 2 + 4k (k = 0,1,2, . . .) and so are all simple poles. Hence, from (13) and (14), we find the algebraic expansion given by H3,0(z) = 3 ∑ j=1 µ j m j z−ν j/m j S3,0(z; j), (27) where S3,0(z; j) = ∞ ∑ k=0 (−)k k! Γ � µ jk+ ν j m j � 3 ∏ r=1 ′ Γ � νr −mrsk, j µr � z−µ j k/m j (28) with the prime denoting the omission of the gamma function factor corresponding to r = j. Then, from (23) we obtain the expansion I3(z) ∼ ¨ E3,0(z) +H3,0(ze∓πi) in |arg z| ≤ 5 12 π− ε H3,0(ze∓πi) in 5 12 π+ ε≤ |arg z| ≤ π (29) as |z| → ∞. It follows that I3(z) is exponentially large in the sector |arg z| < 5 24 π with the dominant expansion being algebraic in the rest of the z-plane. In the sectors 5 12 π < |arg z| < 5 24 π the exponential expansion E3,0(z) is subdominant and switches off (as |arg z| increases) across the Stokes lines arg z = ± 5 12 π. We show in Table 1 the values of the absolute relative error in the computation of I3(z) as a function of θ = arg z when |z| = 15 using the optimally truncated asymptotic expansions (that is, truncated at or near the least term in modulus) in (29). R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1014 Table 1: Values of the absolute relative error in the computation of I3(z) in (25) when |z| = 15 as a function of θ = arg z using an optimal truncation of the expansions in (29). θ/π |Rel. Error| θ/π |Rel. Error| 0 5.816× 10−13 0.625 9.259× 10−14 0.125 9.262× 10−14 0.750 2.744× 10−13 0.250 1.612× 10−13 0.875 9.841× 10−14 0.375 3.534× 10−14 1.000 1.727× 10−13 0.500 1.177× 10−13 4. The integral Jn(z) We can apply a similar treatment to the integral Jn(z) = λn ∫ ∞ −∞ . . . ∫ ∞ −∞ x ν1−1 1 . . . xνn−1 n e− f (x1,...,xn;z) d x1 . . . d xn = λn ∞ ∑ k=0 zk k! n ∏ r=1 ∫ ∞ −∞ xνr+mr k−1 r exp{−xµr r } d xr , where the phase function f and the factor λn are defined in (2) and it is supposed that the parameters µ j (1 ≤ j ≤ n) appearing in f are positive even integers. In the evaluation of the above integrals when xr < 0, we shall write xr = |xr |e πi . We now introduce the notation e(x) for brevity in this section and define the quantities Br(k) by e(x) := eπi x , Br(k) := 1− e(νr +mr k). (30) Then we find µr ∫ ∞ −∞ xνr+mr k−1 r exp{−xµr r } d xr = µr Br(k) ∫ ∞ 0 xνr+mr k−1 r exp{−xµr r } d xr = Br(k)Γ � νr +mr k µr � and hence that Jn(z) = ∞ ∑ k=0 zk k! n ∏ r=1 � Br(k)Γ � νr +mr k µr �� . (31) 4.1. The Representation of Jn(z) in Terms of nΨ0(z) Functions Expansion of the product of exponential factors ∏n r=1 Br(k) in (31) can be achieved by making use of the standard expansion (1− z1)(1− z2) . . . (1− zn) = 1− n ∑ i=1 zi + n−1 ∑ i=1 n ∑ j=i+1 ziz j R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1015 − n−2 ∑ i=1 n−1 ∑ j=i+1 n ∑ k= j+1 ziz jzk + . . .+ (−)nz1z2 . . . zn. Substitution of this expansion into (31) with zr = e(νr +mr k) then shows that Jn(z) may be expressed as a linear combination of nΨ0(z) with rotated argument in the form Jn(z) = nΨ0(z)− n ∑ i=1 e(νi) nΨ0(ze(mi)) + n−1 ∑ i=1 n ∑ j=i+1 e(νi + ν j) nΨ0(ze(mi +m j)) − n−2 ∑ i=1 n−1 ∑ j=i+1 n ∑ k= j+1 e(νi + ν j + νk) nΨ0(ze(mi +m j +mk))+ . . . + (−)ne(N) nΨ0(ze(M)), (32) where we have defined N := ν1 + . . .+ νn, M := m1 + . . .+mn. (33) The parameters appearing in each nΨ0 function are those given in (22). The asymptotic expansion (19), or equivalently (23), can then be employed to deal with each nΨ0(zeπiω) with argument rotated by ω, it being remembered that nΨ0(z) is an integral function of z with arg z evaluated modulo 2π. Rather than attempt to present a complicated general result, we indicate how to proceed with the asymptotic expansion of Jn(z) as |z| →∞ in specific cases in Section 5. 4.2. The Algebraic Contribution to the Expansion of Jn(z) We consider the contribution to the asymptotic expansion of the integral Jn(z) as |z| → ∞ that results from the algebraic expansions associated with each nΨ0 function of rotated argument in (32). It will be shown that this combination of algebraic expansions cancels in the sector ε≤ arg z ≤ (2−M)π− ε, (34) where M is defined in (33), and that consequently the expansion of Jn(z) in this sector is purely exponential in character. We shall suppose in this section that all the poles sk, j in (12) are simple; a case when there are multiple poles present is discussed in Appendix C. Then, the algebraic expansion of nΨ0(z) associated with the parameters in (22) is, from (13), (14) and (16), given by Hn,0(ze−πi) = n ∑ j=1 µ j m j h j(z) valid as |z| →∞ in the sector ε ≤ arg z ≤ 2π− ε, where h j(z) := (ze−πi)−ν j/m j Sn,0(ze−πi; j) (1≤ j ≤ n) (35) R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1016 and Sn,0(ze−πi ; j) = ∞ ∑ k=0 (−)k k! Γ � µ jk+ ν j m j � n ∏ r=1 ′ Γ � νr −mrsk, j µr � (ze−πi)−µ j k/m j . From the representation of Jn(z) in (32) as a finite sum of nΨ0 functions of rotated argument, the contribution from the different algebraic expansions when arg z lies in the common sector (34) can then be written in the form n ∑ j=1 µ j m j n ∑ k=0 (−)kTk, j(z), (36) where T0, j(z) = h j(z), T1, j(z) = n ∑ j1=1 e(ν j1 )h j(ze(m j1 )), T2, j(z) = n−1 ∑ j1=1 n ∑ j2= j1+1 e(ν j1 + ν j2 )h j(ze(m j1 +m j2 )), T3, j(z) = n−2 ∑ j1=1 n−1 ∑ j2= j1+1 n ∑ j3= j2+1 e(ν j1 + ν j2 + ν j3 )h j(ze(m j1 +m j2 +m j3 )), . . . and in general, for 1≤ k ≤ n− 1, Tk, j(z) = n−k+1 ∑ j1=1 n−k+2 ∑ j2= j1+1 . . . n ∑ jk= jk−1+1 e(ν j1 + . . .+ ν jk )h j(ze(m j1 + . . .+m jk )) with Tn, j(z) = e(N)h j(ze(M)). (37) It will be sufficient to consider just one value of j and accordingly we choose j = 1. Since the µ j (1≤ j ≤ n) are even integers it follows from (35) that e(ν j)h j(ze(m j)) = h j(z). (38) We now separate off the term corresponding to j1 = 1 in the sums Tk, j(z) and make repeated use of (38) to find T1,1(z) = h1(z) + n ∑ j1=2 e(ν j1 )h1(ze(m j1 )), T2,1(z) = n ∑ j2=2 e(ν j2 )h1(ze(m j2 ))+ n−1 ∑ j1=2 n ∑ j2= j1+1 e(ν j1 + ν j2 )h1(ze(m j1 +m j2 )), T3,1(z) = n−1 ∑ j2=2 n ∑ j3= j2+1 e(ν j2 + ν j3 )h1(ze(m j2 +m j3 )) R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1017 + n−2 ∑ j1=2 n−1 ∑ j2= j1+1 n ∑ j3= j2+1 e(ν j1 + ν j2 + ν j3 )h1(ze(m j1 +m j2 +m j3 )) and so on. An obvious relabelling of the summation indices then shows that in the inner sum in (36), taken over 0≤ k ≤ n− 1, all the terms cancel except the last to yield n−1 ∑ k=0 (−)kTk,1(z) = (−) n−1 2 ∑ j1=2 3 ∑ j2= j1+1 . . . n ∑ jn= jn−1+1 e(ν j1 + . . .+ ν jn )h1(ze(m j1 + . . .+m jn )) = (−)n−1e(ν2 + . . .+ νn)h1(ze(m2 + . . .+mn)). = (−)n−1e(N)h1(ze(M)) = (−)n−1Tn,1(z) upon application of (38). It therefore follows that ∑n k=1(−) kTk,1(z) ≡ 0. An analogous proce- dure applies to other values of j ≤ n and so we obtain in the simple-pole case n ∑ k=0 (−)kTk, j(z)≡ 0 (1≤ j ≤ n) (39) in the sector (34). The asymptotic expansion of Jn(z) in this sector is therefore purely expo- nential in character. 4.3. The Integral Kn,p(z) The procedure described above can be applied to the variant of the integral Jn(z) obtained by taking p < n integrals evaluated over (−∞,∞) with the rest being evaluated over [0,∞). Thus, if we define Kn,p(z) = λn ∫ ∞ 0 . . . ∫ ∞ 0 �∫ ∞ −∞ . . . ∫ ∞ −∞ x ν1−1 1 . . . xνn−1 n e− f (x1,...,xn;z) d x1 . . . d xp � d xp+1 . . . d xn, (40) where it is now supposed that µr (1 ≤ r ≤ p) are even integers and µr > 0 (p + 1 ≤ r ≤ n), then we obtain following the procedure described in Section 4 the series expansion Kn,p(z) = ∞ ∑ k=0 zk k! n ∏ r=1 Γ � νr +mr k µr � p ∏ r=1 Br(k). (41) In the evaluation of the integrals when xr < 0 (1≤ r ≤ p) we have again taken xr = |xr |e πi . The product ∏p r=1 Br(k) may be expanded as a sum of exponentials so that Kn,p(z) can be written as a finite sum of nΨ0(z) functions of rotated argument and parameters given in (22) in an analogous manner to that in (32). The asymptotic expansion of nΨ0(z) in (19) can then be employed to obtain the expansion of Kn,p(z) as |z| → ∞. We give an example of the asymptotic structure of Kn,p(z) in Section 5. R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1018 5. Numerical Examples In this section we give some numerical examples to illustrate the application of the ex- pansion (19) to the construction of the asymptotic structure of the integrals Jn(z) and Kn,p(z) defined in (4) and (40). 5.1. Example 1 Let us consider the two-dimensional integral J2(z) = µ1µ2 ∫ ∞ −∞ ∫ ∞ −∞ x ν1−1 1 x ν2−1 2 exp{−x µ1 1 − x µ2 2 + zx m1 1 x m2 2 } d x1d x2, (42) where ν1, ν2 > 0 and µ1, µ2 are positive even integers. From (32) and (21), this integral can be expressed in terms of the Wright function 2Ψ0(z) of rotated argument, where 2Ψ0(z)≡ 2Ψ0 � ( m1 µ1 , ν1 µ1 ), ( m2 µ2 , ν2 µ2 ) ; z � . For simplicity in presentation we only consider the case when m1 = m2 = m. Then we find from (32) J2(z) = 2Ψ0(z)− B1 2Ψ0(zeπim) + B2 2Ψ0(ze2πim), (43) where B1 = (e πiν1 + eπiν2), B2 = eπi(ν1+ν2). (44) Provided ν1 and ν2 are real, it is readily shown that J2(z) possesses a basic symmetry about the half-rays arg z =ω and arg z =ω−π, where ω = π(1−m). For, upon recalling that the argument of pΨq(z) can be written modulo 2π, we find‡ from (43) with the above definition of ω that J2(zeiω) = eπi(ν1+ν2) ¦ 2Ψ0(zeiω+2πim)− (e−πiν1 + e−πiν2) 2Ψ0(zeiω+πim) +e−πi(ν1+ν2) 2Ψ0(zeiω) © = eπi(ν1+ν2) ¦ 2Ψ0(ze−iω)− (e−πiν1 + e−πiν2) 2Ψ0(ze−πi) +e−πi(ν1+ν2) 2Ψ0(zeiω−2πi) © = eπi(ν1+ν2)J2(zeiω), (45) where the bar denotes the complex conjugate. Hence it is sufficient in this case to restrict our attention to an appropriate half-plane. We display in Fig. 1 the large-|z| sectorial behaviour of J2(z) in the case µ1 = µ2 = 4 for different values of m, where we suppose that B1 6= 0. In Fig. 1(a), m = 1 4 (κ = 7 8 ) so that the symmetry line is arg z = 3 4 π, −1 4 π and, from (34), the algebraic expansions cancel in the ‡When m1 6= m2, it can be shown that the symmetry relation (45) still holds with ω = π− 1 2 π(m1 +m2). R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1019 sector (0, 3 2 π). There are three overlapping exponentially large sectors |arg(zeπir/4)| < 7 16 π (r = 0,1,2), with the expansion in the sector ( 7 16 π, 17 16 π) being exponentially small. It can be seen that there are sectors in which J2(z) consists of either one, two or three exponentially large expansions. In Fig. 1(b), m = 1 2 (κ = 3 4 ) so that the symmetry line is the imaginary z-axis and the algebraic expansions cancel in the upper half-plane. There are again three ex- ponentially large sectors |arg(zeπir/2)| < 3 8 π (r = 0,1,2), which overlap to produce sectors containing either one or two exponentially large expansions, with the exponentially small sec- tor now being (3 8 π, 5 8 π). In Fig. 1(c), m = 3 4 (κ = 5 8 ) so that the symmetry line is arg z = 1 4 π, −3 4 π and the algebraic expansions cancel in the sector (0, 1 2 π). The exponentially large sectors |arg(ze3πir/4)| < 5 16 π (r = 0,1,2) only overlap in the sector ( 3 16 π, 5 16 π), with the behaviour in the sectors (13 16 π, 15 16 π) and (− 7 16 π,− 5 16 π) consisting of algebraic and exponentially small expansions. Finally in Fig. 1(d), m = 1 (κ = 1 2 ) so that the symmetry line is the real z-axis and, by (34), there is no longer a sector in which the algebraic expansions cancel. In this case, there are two non-overlapping exponentially large sectors given by |arg(zeπir )| < 1 4 π (r = 0,1), with the behaviour of J2(z) in the sectors (1 4 π, 3 4 π) and (−3 4 π,−1 4 π) consisting of algebraic and exponentially small expansions. If the parameters ν1 and ν2 are such that B1 = 0, then the number of exponentially large sectors decreases by one. We now consider the asymptotic expansion of J2(z) in (42) when µ1 = µ2 = 4 in some specific cases in more detail. We first take m= 1 2 (κ = 3 4 ) so that from (43) J2(z) = 2Ψ0(z)− B1 2Ψ0(ze 1 2 πi) + B2 2Ψ0(ze∓πi), (46) where we choose the upper or lower sign according as arg z > 0 or arg z < 0, respectively. Recalling from (19) that the Stokes lines for the exponential expansion Ep,q(z) are given by arg z = ±πκ, we see that the Stokes lines associated with E2,0(ze 1 2 πir) (r = 0,1,2) are the rays arg z = ±3 4 π, ±1 4 π; see Fig. 1(b). Taking into account these Stokes lines and the fact that the algebraic expansions all cancel in the upper half-plane by (34), we find that the exponential expansion of J2(z) in (46) is then given by E2,0(z) + B2 E2,0(ze−πi) in (1 4 π, 1 2 π] E2,0(z)− B1 E2,0(ze 1 2 πi) in (−1 4 π, 1 4 π) E2,0(z)− B1 E2,0(ze 1 2 πi) + B2 E2,0(zeπi) in [−1 2 π,−1 4 π), as |z| → ∞ in the right-half plane. The expansion E2,0(z) is obtained from (8) with Z = 3 8 z4/3 and ϑ = 1 4 (ν1 + ν2)− 1. The coefficient A0 is specified by (24) with the coefficients A j ( j ≥ 1) being determined in specific cases by the algorithm described in Appendix A. The expansion of J2(z) is exponentially small in the sector (3 8 π, 5 8 π). Although there is an algebraic expansion present in the lower half-plane, we do not consider its contribution here as it is subdominant throughout this domain. The exponential expansion in the left-hand half-plane can be obtained via (45). For our second case, we take m = 3 4 (κ = 5 8 ) and ν1 = 1 2 , ν2 = 3 2 , so that from (44) we R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1020 7π/16 3π/16 −π/16 −7π/16−11π/16 −15π/16 (a) 3π/8 π/4 −π/8 −π/4 −3π/8−5π/8−3π/4 −7π/8 5π/8 3π/4 (b) 5π/16 3π/16 −π/8 −5π/16−5π/8 9π/8 5π/8 13π/16 15π/16 −7π/16 (c) π/4 −π/4−3π/4 3π/4 (d) Figure 1: The sectorial behaviour of J2(z) for µ1 = µ2 = 4 and m1 = m2 = m when it is sup- posed that B1 6= 0: (a) m = 1 4 , (b) m = 1 2 , (c) m = 3 4 and (d) m = 1. The sectors marked with a circular arc with arrows denote exponentially large sectors. The hatched regions denote exponentially small behaviour and the shaded regions denote mixed algebraic and exponen- tially small behaviour. The dashed lines are Stokes lines and the dash-dot line is the axis of basic symmetry. have B1 = 0, B2 = 1 and J2(z) = 2Ψ0(z) + 2Ψ0(ze 3 2 πi) = 2Ψ0(z) + 2Ψ0(ze− 1 2 πi). (47) The sector in which the algebraic expansions associated with J2(z) cancel is (0, 1 2 π). Refer- ring to Fig. 2(a), we see that the expansion of J2(z) as |z| → ∞ is exponentially large in the sector (− 5 16 π, 13 16 π), where E2,0(z) is given by (8) with Z = 5 64 (63/5z8/5) and ϑ = −1 2 . In the sectors (13 16 π, 9 8 π) and (−5 8 π,− 5 16 π) the expansion of J2(z) is mixed algebraic and exponen- tially small, whereas due to the presence of the Stokes lines associated with the exponential expansions on arg z = −5 8 π and arg z = 9 8 π the expansion in the sector |arg(ze−3πi/4)| < 1 8 π is purely algebraic. R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1021 From (12) and (13), the algebraic expansion H2,0(z) is controlled by the poles situated at sk,1 = 2 3 + 16 3 k, sk,2 = 2+ 16 3 k (k = 0,1,2, . . .) which are all simple, and consequently H2,0(z) = 16 3 2 ∑ j=1 ∞ ∑ k=0 (−)k k! Γ(sk, j)Γ(γ j − 3 16 sk, j)z −sk, j with γ1 = 3 8 and γ2 = 1 8 . Taking into account the Stokes lines on arg z = 0 and arg z = 1 2 π for the algebraic expansions associated with 2Ψ0(z) and 2Ψ0(ze−πi/2), respectively, we see from (15) and (16) that the algebraic expansion H(z) of J2(z) is given by H(z) = ( H2,0(ze−πi) +H2,0(ze 1 2 πi)≡ 0 in (0, 1 2 π) H2,0(ze−πi) +H2,0(ze− 3 2 πi) in (1 2 π, 2π) It is easily verified with the above form of H2,0(z) that H(z) ≡ 0 in the sector (0, 1 2 π), in accordance with (34). Then the asymptotic expansion of J2(z) in (47) has the form J2(z)∼        E2,0(z) + E2,0(ze− 1 2 πi) in (0, 1 2 π) E2,0(z) + E2,0(ze− 1 2 πi) +H(z) in (1 2 π, 5 8 π) E2,0(ze− 1 2 πi) +H(z) in (5 8 π, 9 8 π) H(z) in (9 8 π, 5 4 π] as |z| →∞, with that in the remainder of the plane being determined by the symmetry relation (45). 5π/16 3π/16 −π/8 −5π/16−5π/8 9π/8 5π/8 13π/16 (a) 5π/12 5π/24 −5π/24 −5π/12−7π/12 −19π/24 7π/12 19π/24 (b) Figure 2: The sectorial behaviour of (a) J2(z) for µ1 = µ2 = 4 and m1 = m2 = 3 4 when it is supposed that B1 = 0 and (b) of K3,1(z) for µ1 = µ2 = 3, µ3 = 4 and m1 = m2 = 1 2 , m3 = 1. The sectors marked with a circular arc with arrows denote exponentially large sectors. The shaded regions denote mixed algebraic and exponentially small behaviour. The dashed lines are Stokes lines and the dash-dot line is the axis of basic symmetry. R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1022 Finally, we consider m= 1 (κ = 1 2 ) and ν1 = ν2 = 1 4 , so that B1 = 2eπi/4 and B2 = i and J2(z) = (1+ i) 2Ψ0(z)− 2e 1 4 πi 2Ψ0(ze−πi), (48) where we have replaced the argument of the second Wright function by ze−πi . From Fig. 1(d), the Stokes lines coincide with the imaginary axis and the symmetry axis is the real axis. The exponential expansion E2,0(z) in (8) has Z = 1 8 z2 and ϑ = −7 8 . The poles in (12) are all double situated at sk = 4k+ 1 4 (k = 0,1,2, . . .) and, from (60), we obtain the algebraic expansion in this case given by H2,0(z) = −16 ∞ ∑ k=0 Γ(4k+ 1 4 ) (k!)2 {ψ(4k+ 1 4 )− 1 2 ψ(k+ 1)− log z}z−4k− 1 4 . Since, by (34), there is no sector in which the algebraic expansions cancel when m = 1, we obtain the expansion of J2(z) in (48) given by J2(z) ∼    (1+ i){E2,0(z) +H2,0(ze−πi)}− 2e 1 4 πiH2,0(z) in [0, 1 2 π) (1+ i)H2,0(ze−πi)− 2e 1 4 πi{E(ze−πi) +H2,0(z)} in (1 2 π,π] as |z| →∞. The expansion in the lower half-plane can be obtained by (45). The rays arg z = 0, π are Stokes lines for the algebraic expansions H2,0(ze−πi) and H2,0(z), respectively and the rays arg z = ±1 2 π are Stokes lines for the exponential expansions. In Table 1 we present the absolute relative errors in the asymptotic expansion of J2(z) in (46), (47) and (48) for a given value of |z| and varying θ = arg z. In each case the exponential and algebraic expansions have been optimally truncated, with the exact value of J2(z) being computed both by evaluation of the Wright functions and also high-precision numerical quadrature of the integral in (42). We remark that in the three cases considered, an accurate determination of the subdominant expansions on the Stokes lines would require a detailed treatment of the Stokes phenomenon. 5.2. Example 2 We consider an example of the integral Kn,p(z) defined in (40) with n= 3, p = 1, namely K3,1(z) = 36 ∫ ∞ −∞ �∫ ∞ 0 ∫ ∞ 0 (x3/x2) 1 2 exp{−(x3 1 + x3 2 + x4 3 − z(x1 x2) 1 2 x3} d x1d x2 � d x3, which is associated with the parameters µ1 = µ2 = 3, µ3 = 4, m1 = m2 = 1 2 , m3 = 1 and ν1 = 1, ν2 = 1 2 , ν3 = 3 2 . From (41), we therefore find K3,1(z) = ∞ ∑ k=0 zk k! 3 ∏ r=1 Γ � νr +mr k µr � (1− e(ν3 +m3k)) R. Paris / Eur. J. Pure Appl. Math, 3 (2010), 1006-1031 1023 Table 2: Values of the absolute relative error in the computation of J2(z) when µ1 = µ2 = 4 and m1 = m2 = m as a function of θ = arg z. m= 1 2 , ν1 = 1 2 , ν2 = 1 m= 3 4 , ν1 = 1 2 , ν2 = 3 2 m = 1, ν1 = ν2 = 1 4 |z| = 20 |z| = 20 |z| = 15 θ/π |Rel. Error| θ/π |Rel. Error| θ/π |Rel. Error| 0 2.172× 10−10 0.250 8.081× 10−14 0 4.869× 10−11 0.125 3.326× 10−10 0.375 2.528× 10−13 0.125 9.220× 10−12 0.250 1.228× 10−9 0.500 1.104× 10−12 0.250 4.427× 10−14 0.375 3.939× 10−8 0.625 2.528× 10−13 0.375 1.025× 10−15 0.500 1.756× 10−4 0.750 1.414× 10−13 0.500 5.061× 10−15 −0.125 3.326× 10−10 0.875 7.550× 10−14 0.625 1.046× 10−15 −0.250 3.684× 10−10 1.000 7.014× 10−14 0.750 5.306× 10−14 −0.375 3.326× 10−10 1.125 2.584× 10−13 0.875 9.220× 10−12 −0.500 2.172× 10−10 1.250 1.058× 10−13 1.000 4.869× 10−11 = 3Ψ0(z) + i 3Ψ0(ze∓πi), where the 3Ψ0(z) function is that defined in (25) with the same parameter values. It is easily shown that K3,1(ze±πi) = iK3,1(z), where the bar denotes the complex conju- gate, so that there is a basic symmetry about the imaginary axis. The sectorial behaviour of K3,1(z) is shown in Fig. 2(b) with κ = 5 12 . There are two exponentially large sectors, four sectors with mixed algebraic and exponentially small behaviour and two sectors straddling the imaginary axis of angular width 1 6 π (bounded by the Stokes lines) in which the large-|z| behaviour is algebraic. From the asymptotic expansion of 3Ψ0(z) given in (29) we then obtain K3,1(z) ∼      E3,0(z) +H3,0(ze−πi) + iH3,0(z) in (0, 5 12 π) H3,0(ze−πi) + iH3,0(z) in ( 5 12 π, 1 2 π] E3,0(z) +H3,0(zeπi) + iH3,0(z) in (− 5 12 π, 0) H3,0(zeπi) + iH3,0(z) in [−1 2 π,− 5 12 π) as |z| → ∞, where the expansions E3,0(z) and H3,0(z) are given in (26), (27) and (28). The expansion of K3,1(z) in the left-hand half-plane is described by the above symmetry relation. 6. Concluding Remarks We have shown how the n-dimensional analogues of Faxén’s integral in (1), (4) and (40) can be expressed in terms of either a single Wright function nΨ0(z), or a linear combination of such functions with rotated arguments. Knowledge of the asymptotic expansion of nΨ0(z) for |z| → ∞ then enables the asymptotic structure of these integrals to be determined. Not surprisingly, this asymptotic structure becomes more complicated the larger the value of n. The asymptotic behaviour of nΨ0(z) for large |z| consists of an exponential expansion and an algebraic expansion. The formal sum En,0(z) is a compact representation of the exponential REFERENCES 1024 expansion in the n-dimensional case. The evaluation of the coefficients in this expansion can be easily carried out in specific cases for low values of n, although the computational effort in- volved in their calculation rapidly increases with the dimension of the integrals. The algebraic expansion Hn,0(z) consists, in general, of n different expansions each with its own asymptotic scale. The occurrence of terms in log z (when n≥ 2) depends to a considerable degree on the symmetry in the associated Newton diagram of the phase function f (x1, . . . , xn; z) defined in (2). References [1] N G Bakhoom, Asymptotic expansions of the function Fk(x) = ∫∞ 0 exp(xu−uk)du, Proc. Lond. Math. Soc. 35:83–100, 1935. [2] B L J Braaksma, Asymptotic expansions and analytic continuations for a class of Barnes- integrals, Compos. Math. 15:239–341, 1963 [3] S Breen and A D Wood, On the asymptotic behaviour of Laplace-type multiple integral solutions of linear differential equations, J. Comp. Appl. Math. 171:103–112, 2004. [4] L Brillouin, Sur une méthode de calcul approcheé de certaines intégrales, dite méthode de col, Ann. Sci. Ecole Norm. Sup. 33:17–69, 1916. [5] W R Burwell, Asymptotic expansions of generalised hypergeometric functions, Proc. Lond. Math. Soc. 53:599–611, 1924. [6] R Gorenflo, Y Luchko and F Mainardi, Analytical properties and applications of the Wright function, Frac. Calc. Appl. Anal. 2:383–414, 1998. [7] D Kaminski and R B Paris, Asymptotics via iterated Mellin-Barnes integrals: Application to the generalised Faxén integral, Methods Appl. Anal. 4:311–325, 1997. [8] G V Liakhovetski and R B Paris, Asymptotic expansions of Laplace-type integrals III, J. Comp. Appl. Math. 132:409–429, 2001. [9] F W J Olver, Asymptotics and Special Functions, Academic Press, New York, 1974. Reprinted in A K Peters, Massachussets, 1997. [10] R B Paris, Smoothing of the Stokes phenomenon for high-order differential equations, Proc. Roy. Soc. London A 436:165–186, 1992. [11] R B Paris, Exponentially small expansions in the asymptotics of the Wright function, J. Comp. Appl. Math. 234:488–504, 2010. [12] R B Paris and D Kaminski, Asymptotics and Mellin-Barnes Integrals, Cambridge University Press, Cambridge, 2001. [13] R B Paris and G V Liakhovetski, Asymptotics of the multidimensional Faxén integral, Frac. Calc. Appl. Anal. 3:63–73, 2000. REFERENCES 1025 [14] R B Paris and A D Wood, Asymptotics of High Order Differential Equations, Pitman Re- search Notes in Mathematics, 129, Longman Scientific and Technical, Harlow, 1986. [15] T D Riney, On the coefficients in asymptotic factorial expansions, Proc. Amer. Math. Soc. 7:245–249, 1956. [16] R Saxton, An integral representation solution for a class of higher order linear ordinary differential equations, MSc Dissertation, Cranfield Institute of Technology, 1978. [17] L J Slater, Generalised Hypergeometric Functions, Cambridge University Press, Cam- bridge, 1966. [18] S Spitzer, Integration der linearen Differentialgleichung y(n) = Ax2 y ′′ + Bx y ′ + C y, in welcher n eine ganze positive Zahl und A, B, C constante Zahlen bezeichnen, mittelst bestimmter Integrale, Math. Ann. 3:453–455, 1871. [19] N M Temme, Special Functions: An Introduction to the Classical Functions of Mathematical Physics, Wiley, New York, 1996. [20] E M Wright, The asymptotic expansion of the generalized hypergeometric function, Proc. Lond. Math. Soc. (Ser. 2) 10:286–293, 1935. [21] E M Wright, The asymptotic expansion of the generalized hypergeometric function, Proc. Lond. Math. Soc. (Ser. 2) 46:389–408, 1940. [22] E M Wright, A recursion formula for the coefficients in an asymptotic expansion, Proc. Glasgow Math. Assoc. 4:38–41, 1958. Appendix A. An algorithm for the computation of the coefficients c j = A j/A0 We describe an algorithm for the computation of the normalised coefficients A j/A0 ap- pearing in the exponential expansion Ep,q(z) in (8). Methods of computing these coefficients by recursion in the case when αr = βr = 1 have been given by Riney [15] and Wright [22]. Here we describe an algebraic method valid for arbitrary αr > 0 and βr > 0; see also [12, pp. 46–49]. We rewrite the inverse factorial expansion (9) in the form g(s)Γ(κs+ ϑ′) Γ(s+ 1) = κA0(hκ κ)s �M−1 ∑ j=0 c j (κs+ ϑ′) j + O(1) (κs+ ϑ′)M � , (49) for |s| → ∞ uniformly in |arg s| ≤ π− ε, where g(s) is the ratio of gamma functions defined in (5), (a) j = Γ(a+ j)/Γ(a) and c j = A j/A0. Introduction of the scaled gamma function Γ∗(z) defined by Γ∗(z) := Γ(z)(2π)− 1 2 ezz 1 2 −z REFERENCES 1026 leads to the representation Γ(αs+ a) = Γ∗(αs+ a)(2π) 1 2 e−αs(αs)αs+a− 1 2 e(αs; a), where e(αs; a) = exp � (αs+ a− 1 2 ) log � 1+ a αs � − a � . Some straightforward algebra then shows that the left-hand side of (49) becomes g(s)Γ(κs+ ϑ′) Γ(s+ 1) = κA0(hκ κ)sR(s)Υ(s), (50) where Υ(s) = ∏p r=1 Γ ∗(αrs+ ar) ∏q r=1 Γ ∗(βrs+ br) Γ∗(κs+ ϑ′) Γ∗(s+ 1) and R(s) = ∏p r=1 e(αrs; ar) ∏q r=1 e(βrs; br) e(κs;ϑ′) e(s; 1) . Substitution of (50) into (49) finally produces R(s)Υ(s) = M−1 ∑ j=0 c j (κs+ ϑ′) j + O(1) (κs+ ϑ′)M (51) as |s| →∞ in |arg s| ≤ π− ε. Now let χ = (κs)−1 and expand R(s) and Υ(s) for χ → 0 making use of the well-known expansion [19, p. 71], [12, p. 32] Γ∗(z)∼ ∞ ∑ k=0 (−)kγkz−k (|z| →∞; |arg z| ≤ π− ε), where γk are the Stirling coefficients. The first few coefficients are given by γ0 = 1, γ1 = − 1 12 , γ2 = 1 288 , γ3 = 139 51840 , . . . . Some routine algebra then yields Γ∗(αs+ a) = 1− γ1κχ α +O(χ2), e(αs; a) = 1+ κχ 2α a(a− 1) +O(χ2), whence R(s) = 1+ κχ 2 ( p ∑ r=1 ar(ar − 1) αr − q ∑ r=1 br(br − 1) βr − ϑ κ (1− ϑ) ) +O(χ2), Υ(s) = 1+ κχ 12 ( p ∑ r=1 1 αr − q ∑ r=1 1 βr + 1 κ − 1 ) +O(χ2). REFERENCES 1027 Upon equating coefficients of χ in (51) we obtain c1 = 1 2 κ(A + 1 6 B), (52) where A = p ∑ r=1 ar(ar − 1) αr − q ∑ r=1 br(br − 1) βr − ϑ κ (1− ϑ), B = p ∑ r=1 1 αr − q ∑ r=1 1 βr + 1 κ − 1. The higher coefficients are then obtained by continuation of this expansion process applied to R(s) and Υ(s) in (51) with the help of Mathematica. In specific cases (i.e., with numerical values for the various parameters) it is possible to generate the coefficients in this manner quite easily. In our computations we have used up to a maximum of 40 coefficients. Appendix B. The Asymptotic Expansion of 2Ψ0(z) We demonstrate the validity of the assertion in (19) concerning the asymptotic expansion of the Wright function pΨq(z) as |z| → ∞ by considering a particular case. Let us take p = 2, q = 0 with the parameter values α1 = α2 = 1 4 and ν1 = ν2 = 1 8 ; that is, we consider the function 2Ψ0(z) ≡ 2Ψ0 � (1 4 , 1 8 ), (1 4 , 1 8 ) ; z � = ∞ ∑ k=0 zk k! Γ2(1 4 k+ 1 8 ). (53) From (7), this function is associated with the parameters κ = 1 2 , h = 1 2 and ϑ = −3 4 . The exponential expansion is, from (8) and (10), then given by E2,0(z) = Z−3/4eZ ∞ ∑ j=0 A j Z − j, Z = 1 8 z2, (54) where A0 = 29/4π; we have employed coefficients with j ≤ 28 in our computations (see Appendix A). The first ten coefficients c j ≡ A j/A0 for 2Ψ0(z) in (53) are listed in Table 3. Table 3: The coefficients c j for 1≤ j ≤ 10 associated with the function in (53). j c j j c j 1 13 16 2 729 512 3 31575 8192 4 7432635 524288 5 554191155 8388608 6 100179200205 268435456 7 10645956497295 4294967296 8 10406881110208275 549755813888 9 1437596137005803775 8796093022208 10 443063017349580803175 281474976710656 REFERENCES 1028 From (11), the Mellin-Barnes integral representation for 2Ψ0(z) in (53) is given by 2Ψ0(z) = 1 2πi ∫ ∞i −∞i Γ(s)Γ2(1 8 − 1 4 s)(ze∓πi)−sds (|arg(−z)| < 3 4 π), where the upper or lower sign is chosen according as arg z > 0 or arg z < 0, respectively and the integration path separates the poles of Γ(s) from the sequence of double poles of Γ2(1 8 − 1 4 s) situated at s = 4k+ 1 2 , k = 0,1,2, . . . . Displacement of the integration path over the sequence of double poles and evaluation of the residues leads to the algebraic expansion H2,0(ze∓πi) = −16 ∞ ∑ k=0 Γ(4k+ 1 2 ) (k!)2 (ze∓πi)−4k− 1 2 � ψ(4k+ 1 2 )− 1 2 ψ(k+ 1)− log(ze∓πi) � , (55) where ψ denotes the logarithmic derivative of the gamma function; compare (60). Theorems 1 and 2 show that 2Ψ0(z) is exponentially large given by (54) in the sector |arg z| < 1 4 π, with the dominant expansion in the rest of the z-plane being the algebraic expansion in (55). The exponential expansion E2,0(z), which is subdominant in the sectors 1 4 π < |arg z| < 3 4 π, is maximally subdominant on the rays (the Stokes lines) arg z = ±1 2 π. Accordingly, as |arg z| increases, the expansion E2,0(z) should undergo a Stokes phenomenon and switch off smoothly across the rays arg z = ±1 2 π, to leave the algebraic expansion H2,0(ze∓πi) in the sectors 1 2 π < |arg z| ≤ π, as stated in (19); that is 2Ψ0(z) ∼ ¨ E2,0(z) +H2,0(ze∓πi) in |arg z| ≤ 1 2 π− ε H2,0(ze∓πi) in 1 2 π+ ε≤ |arg z| ≤ π (56) as |z| →∞. To demonstrate this, we set z = |z|eiθ and define the Stokes multiplier S(θ) (at fixed |z|) for the Stokes line arg z = 1 2 π by 2Ψ0(z) = H opt 2,0 (ze∓πi) + A0Z−3/4eZ S(θ), where the superscript ‘opt’ denotes that the algebraic expansion is truncated at its optimal truncation point and Z is defined in (54). In the first half of Table 4 we show the values§ of Re (S) for varying θ in the neighbourhood of θ = 1 2 π when |z| = 15, where the value of 2Ψ0(z) has been computed by high-precision summation of (53). The second half of Table 4 displays the absolute error in the computation of 2Ψ0(z) using the asymptotic expansion in (56) in the sector 1 2 π≤ θ ≤ 3 4 π with the same value of |z|. The values in the column labelled (a) were obtained using the first expansion in (56), that is with the exponential expansion retained in the sector 1 2 π < θ ≤ 3 4 π, whereas those in the column labelled (b) were obtained using the second expansion in (56). Both asymptotic series were truncated at their respective optimal truncation points. The first half of Table 4 confirms that the exponential expansion E2,0(z) switches off (as arg z increases) across the Stokes line arg z = 1 2 π to leave the algebraic §The Stokes multiplier S(θ ) has a small imaginary part that we do not show. REFERENCES 1029 Table 4: The variation of the real part of the Stokes multiplier S(θ) and the absolute error in the computation of 2Ψ0(z) for different θ in the sector 1 2 π ≤ θ ≤ 3 4 π when |z| = 15: (a) with E2,0(z) and (b) without E2,0(z). θ/π Re(S) θ/π Re(S) θ/π |Error| (a) |Error| (b) 0.40 0.97776 0.51 0.32448 0.50 4.010× 10−13 3.662× 10−13 0.45 0.94397 0.52 0.20592 0.55 2.819× 10−12 1.324× 10−13 0.47 0.83063 0.53 0.11579 0.60 1.571× 10−10 7.241× 10−14 0.48 0.72932 0.54 0.05599 0.65 7.955× 10−8 4.808× 10−14 0.49 0.60226 0.55 0.02191 0.70 2.036× 10−4 3.569× 10−14 0.50 0.46177 0.60 0.00008 0.75 1.222× 10−0 2.854× 10−14 expansion H2,0(ze−πi) in the remainder of the upper half-plane; a similar behaviour applies across the Stokes line arg z = −1 2 π. The values of the absolute error in the second half of the table clearly indicate that a uniform accuracy over the sector 1 2 π < θ ≤ 3 4 π is achievable by discarding the exponential expansion in this sector in accordance with the second expansion in (56). We remark that a detailed analysis of the Stokes multiplier has been carried out in [10] for the more general function pΨ0(z) with the parameters αr = 1/n (1 ≤ r ≤ p), for positive integer n, and general¶ ar (1 ≤ r ≤ p). It was shown that for large |z| the leading behaviour of the Stokes multiplier S(θ) across the Stokes lines θ = ±πκ is given by S(θ) ≃ 1 2 ± 1 2 erf [(θ ∓πκ)(2κ/n)−1/2(|z|/n)1/(2κ)] (|z| →∞) respectively, where erf denotes the error function and κ = 1− (p/n). Specialisation to the values p = 2, n = 4, to correspond to (53), shows the smooth transition of S(θ) across the Stokes lines θ = ±1 2 π, thereby confirming the above viewpoint. Appendix C. The Algebraic Expansion of Jn(z) in the Case of Double Poles When some, or all, of the poles in (12) are multiple, the analysis of the algebraic con- tributions to the integral Jn(z) presented in Section 4.2 no longer applies. The treatment of the multiple-pole case in general would be very tedious. Accordingly, we demonstrate in the case n = 2 when double poles are present that the cancelation of the algebraic expansions associated with Jn(z) continues to hold in the sector (34), where m = m1 +m2. The algebraic expansion for the Wright function 2Ψ0(z) associated with the parameters αr , ar (r = 1,2) as |z| → ∞ is given by (16) with p = 2, q = 0. The form of the expansion H2,0(ze−πi) in (13) and (14) has to be modified to take into account the presence of the double poles. From the Mellin-Barnes representation in (11), we have 2Ψ0(z) = 1 2πi ∫ ∞i −∞i Γ(s)Γ(a1 −α1s)Γ(a2 −α2s)(ze±πi)−sds (57) ¶It is assumed the parameters are such that only simple poles arise in the corresponding integral (11). REFERENCES 1030 valid in |arg(−z)| < 1 2 π(1 + α1 + α2), where the upper or lower sign is chosen according as arg z > 0 or arg z < 0, respectively. The two sequences of poles that contribute to the algebraic expansion are, from (12), sm,r = (ar +m)/αr (r = 1,2; m = 0,1,2, . . .), where we suppose for some nonnegative integers k, ℓ that sk,1 = a1 + k α1 = a2 + ℓ α2 (58) for double poles to arise. The algebraic expansion H2,0(ze−πi), obtained by displacement of the integration path in (57) to the right over the above poles, then becomes H2,0(ze−πi) = 2 ∑ j=1 α−1 j (ze−πi)−a j/α j S′2,0(ze−πi ; j) + G(z) (59) where the prime denotes the deletion of the terms in the asymptotic sum (14) corresponding to the double poles. The contribution G(z) resulting from the double poles may be shown to be G(z) = 1 α1α2 ∑ k,ℓ (−)k+ℓ+1 k!ℓ! Γ(sk,1)(ze−πi)−sk,1{ψ(sk,1)−α1ψ(k+1)−α2ψ(ℓ+1)− log(ze−πi)}, (60) where summation is over the integers k, ℓ satisfying (58) and ψ(z) denotes the logarithmic derivative of the gamma function. The integral J2(z) is associated with the function 2Ψ0(z) with the parameters αr = mr/µr , ar = νr/µr (r = 1,2), where µr are positive even integers. The contribution to the algebraic expansion of J2(z) from the first series on the right-hand side of (59) (resulting from the simple poles) vanishes in the sector (34) by virtue of the discussion in Section 4.2. To deal with the contribution from the double poles, we note that G(z) may be written in the form G(z) = ∑ k,ℓ (ze−πi)−λk{ck,ℓ+ dk,ℓ log (ze−πi)}, where ck,ℓ and dk,ℓ are coefficients independent of z and λk := ν1 +µ1k m1 = ν2 +µ2ℓ m2 . (61) Recalling that e(x)≡ exp(πi x), we then see that e(ν1)G(ze(m1)) = G(z) +πim1 ∑ k,ℓ dk,ℓ(ze−πi)−λk . From (32) with n= 2, the contribution to the algebraic expansion of J2(z) resulting from the logarithmic series G(z) in the common sector (34) is then G(z) − e(ν1)G(ze(m1))− e(ν2)G(ze(m2)) + e(ν1 + ν2)G(ze(m1 +m2)) REFERENCES 1031 = −πim1 ∑ k,ℓ dk,ℓ(ze−πi)−λk − e(ν2)G(ze(m2)) +e(ν2) � G(ze(m2)) +πim1 ∑ k,ℓ dk,ℓ(ze(m2 − 1))−λk � = πim1 ∑ k,ℓ dk,ℓ(ze−πi)−λk {e(ν2 −m2λk)− 1} ≡ 0 since e(ν2 − m2λk) = 1 by (61). Thus, the algebraic expansion associated with J2(z) when double poles are present similarly vanishes in the sector (34).