/compile/output.dvi EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 9, No. 1, 2016, 3-18 ISSN 1307-5543 – www.ejpam.com The Asymptotic Expansion of a Generalisation of the Euler-Jacobi Series Richard B. Paris Division of Computing and Mathematics, University of Abertay Dundee, Dundee DD1 1HG, UK Abstract. We consider the asymptotic expansion of the sum Sp(a; w) = ∞ ∑ n=1 e−anp nw as a→ 0 in |arg a|< 1 2π for arbitrary finite p > 0 and w> 0. Our attention is concentrated mainly on the case when p and w are both even integers, where the expansion consists of a finite algebraic expan- sion together with a sequence of increasingly subdominant exponential expansions. This exponentially small component produces a transformation for Sp(a; w) analogous to the well-known Poisson-Jacobi transformation for the sum with p = 2 and w= 0. Numerical results are given to illustrate the accuracy of the expansion obtained. 2010 Mathematics Subject Classifications: 30E15, 33B10, 34E05, 41A30 Key Words and Phrases: Euler-Jacobi series, Poisson-Jacobi transformation, asymptotic expansion, inverse factorial expansion 1. Introduction We consider the asymptotic expansion of the sum Sp(a; w) = ∞ ∑ n=1 e−anp nw (1) as the parameter a→ 0 in |arg a|< 1 2π, where p > 0 and, for convenience, w will be supposed throughout to be real and positive. When w= 0, this sum is known as the Euler-Jacobi series and when a = 0 then Sp(0; w) reduces to the Riemann zeta function ζ(w) (providedℜ(w)> 1). Consequently, the series in (1) can also be viewed as a smoothed Dirichlet series for ζ(w). The asymptotics of Sp(a; w) as a→ 0+ for p a rational fraction and w< 0 was considered by Ramanujan and is discussed in [1, Chapter 15]. It was shown that the expansion in this Email address: r.paris@abertay.ac.uk http://www.ejpam.com 3 c© 2016 EJPAM All rights reserved. R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 4 case consisted of an asymptotic sum involving the Riemann zeta function. A hypergeometric function approach for the Euler-Jacobi series when p is a rational fraction and w= 0 has been discussed at length in the monograph [3], where great emphasis was placed on obtaining ex- ponentially small expansions. The case w= 0 and arbitrary p > 0, thus generalising the work in [3], has been investigated in [5, §8.1] also for a→ 0+ using both a Mellin-Barnes integral approach and also a saddle point analysis of a Laplace-type integral representation. A similar Mellin-Barnes integral approach for this latter case has also been independently considered. The saddle point approach is instructive for understanding the appearance of exponentially small terms in the expansion as the parameter p increases. It was established in [5, §8.1] that an additional exponentially small contribution appears when p = 2+ 4k, k = 0,1,2, . . . . Indeed, the appearance of the first exponentially small expansion as p passes through the ‘clas- sical’ value p = 2 was demonstrated to be associated with a Stokes phenomenon [5, §8.1.7]. In the case p = 2, w= 0 the behaviour of the sum S2(a; 0) as a→ 0 can be obtained from the classical Poisson-Jacobi transformation given by S2(a; 0) = ∞ ∑ n=1 e−an2 = 1 2 s π a − 1 2 + s π a ∞ ∑ n=1 e−π 2n2/a (2) valid for all values of a in |arg a| < 1 2π. This well-known transformation relates a sum of Gaussian exponentials involving the parameter a to a similar sum with parameter π2/a. In the small-a limit, the convergence of the sum on the left-hand side becomes slow, whereas the sum on the right-hand side converges rapidly. Various proofs of (2) exist in the literature; see, for example, [5, p. 120], [7, p. 60] and [8, p. 124]. The dominant asymptotic expansion of Sp(a; w) as a → 0 in |arg a| < 1 2π for general p > 0 and w > 0 is relatively straightforward and is found to consist, in general, of a single term proportional to a(w−1)/p, together with a series in ascending powers of a with coefficients involving the Riemann zeta function (the algebraic expansion). When 0 < p ≤ 1, the expan- sion is convergent and the result is exact; when p > 1 the expansion is asymptotic as a→ 0. The most interesting case arises when p and w are both even integers. The above-mentioned algebraic expansion then terminates after a finite number of terms, and it becomes essential for accurate estimation to also include a subdominant sequence of exponentially small expan- sions. This exponentially small component produces a transformation for Sp(a; w) analogous to the Poisson-Jacobi transformation in (2), but valid as a→ 0 in |arg a|< 1 2π. This similarly involves a finite sequence of series similar to (2) with a in the exponential replaced by an inverse power of a, but with each term decorated by an asymptotic series in ascending powers of a1/(p−1). The approach we employ in this paper is based on a Mellin-Barnes integral representation for Sp(a; w) and is similar to that described in [5, §8.1.4]. An algorithm for the determination of the coefficients in the exponentially small asymptotic series is described. The case p = 2 when w is an even integer has been recently discussed in [6], where the coefficients in the decorating asymptotic series can be given in closed form. An application of the series when p = 2, with w = 2 and w = 4, has arisen in the geological problem of thermochronometry in spherical geometry [9]. R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 5 2. An Expansion for Sp(a; w) as a→ 0 When w, p 6= 2, 4, . . . We examine the expansion of the series Sp(a; w) defined in (1) as a → 0 in the sector |arg a|< 1 2π, where p > 0 and for simplicity in presentation we shall assume throughout real values of w > 0. The case p = 1, w = 0 may be excluded from our consideration since the series in this case is summable as a geometric progression. Our starting point is the well-known Cahen-Mellin integral (see, for example, [5, §3.3.1]) zαe−z = 1 2πi ∫ c+∞i c−∞i Γ(α− s)zsds (z 6= 0, |arg z|< 1 2 π), (3) where c < ℜ(α) so that the integration path passes to the left of all the poles of Γ(α − s) situated at s = k+α (k = 0,1,2, . . .). Then, it follows that Sp(a; w) = ∞ ∑ n=1 e−anp nw = ∞ ∑ n=1 n−w 2πi ∫ −c+∞i −c−∞i Γ(−s)(anp)sds = 1 2πi ∫ −c+∞i −c−∞i Γ(−s)ζ(w− ps)asds (|arg a|< 1 2 π), (4) upon reversal of the order of summation and integration, which is justified when c >max{0, (w−1)/p}, and evaluation of the inner sum in terms of the Riemann zeta function. The integrand in (4) possesses simple poles at s = k (k = 0,1,2, . . .) and s = s0 ≡ (w−1)/p, except if s0 = M (that is, w = pM + 1), where M is a non-negative integer, when the pole at s = s0 is double. The residue at the double pole is obtained by making use of the fact that ζ(s)≃ 1/(s− 1) + γ in the neighbourhood of s = 1, where γ is Euler’s constant, to find (−a)M M ! § γ− 1 p log a+ 1 p ψ(M + 1) ª (M = 0,1,2, . . .), where ψ(x) is the logarithmic derivative of the gamma function. The case when w and p are even positive integers requires a separate treatment which is discussed in Section 3. 2.1. The Case 0< p < 1 We first consider the case 0 < p < 1. The integration path in (4) can be made to coincide with the imaginary s-axis together with a suitable indentation to lie to the left of the poles at s = 0 and s = s0 (when 0 < w < 1). Then use of the functional relation for ζ(s) given by [8, p. 269] ζ(s) = 2sπs−1ζ(1− s)Γ(1− s) sin 1 2 πs (5) shows that the integrand can be written as (2π)wζ(1−w+ ps) Γ(1−w+ ps) Γ(1+ s) sin 1 2π(ps−w) sinπs as (2π)ps . R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 6 With s = Reiθ , where R→∞ and is chosen so that the arc passes between the poles on the positive real axis, the logarithm of the dominant real part of the integrand is controlled by (p− 1)R cosθ log R+O(R). When |θ | < 1 2π, this last expression tends to −∞ when 0 < p < 1. Consequently, the inte- gration path can be bent back to enclose the poles of the integrand to yield the convergent result Sp(a; w) = Jp(a; w) + ∞ ∑ k=0 ′ (−)k k! ζ(w− kp)ak (0< p < 1), (6) where Jp(a; w) =      1 p Γ � 1−w p � a(w−1)/p (w 6= pM + 1) (−a)M M ! {γ− 1 p log a+ 1 p ψ(M + 1)} (w= pM + 1) and the prime on the sum over k denotes the omission of the term corresponding to k = M when w= pM + 1. When p = 1, the sum S1(a; w) is given by (6) but now the sum over k on the right-hand side converges when |a| < 2π; see [5, §4.2.2] for details. If we let p = 1, w = 0 then use of the facts that ζ(1− 2k) = −B2k/(2k) and ζ(−2k) = 0, where B2k are even-order Bernoulli numbers, shows that S1(a; 0) = 1 a − 1 2 + 1 a ∞ ∑ k=1 B2k (2k)! a2k (|a|< 2π), which correctly reduces to the trivial summation 1/(ea−1) by application of [4, Eq. (24.2.1)]. 2.2. The Case p > 1 When p > 1, the integration path in (4) cannot be bent back over the poles and we proceed in a similar manner to that described for the case w = 0 in [5, §8.1.4]. Consider the integral taken round the rectangular contour with vertices at −c ± iT , c′ ± iT , where c′ > 0. The contribution from the upper and lower sides s = σ ± iT , −c ≤ σ ≤ c′, vanishes as T →∞ provided |arg a|< 1 2π, since from the behaviour Γ(σ± i t) = O(tσ− 1 2 e− 1 2πt), ζ(σ± i t) = O(tµ(σ) logA t) (t →∞), where for σ and t real µ(σ) = 0 (σ > 1), 1 2 − 1 2 σ (0≤ σ ≤ 1), 1 2 −σ (σ < 0), A= 1 (0≤ σ ≤ 1), A= 0 otherwise, R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 7 the modulus of the integrand is controlled by O(Tσ+µ(σ)− 1 2 log Te−∆T ), with∆= 1 2π−|arg a|. Displacement of the integration path to the right over a finite set of poles then yields (provided w and p are not even integers) Sp(a; w) = Jp(a; w) + N−1 ∑ k=0 ′ (−)k k! ζ(w− kp)ak + RN , (7) where N is a positive integer such that N > s0 + 3 2 and the prime on the sum over k again denotes the omission of the term corresponding to k = M when w= pM + 1. The remainder RN is given by RN = 1 2πi ∫ c+∞i c−∞i Γ(−s)ζ(w− ps)asds = (2π)w 2πi ∫ c+∞i c−∞i ζ(1−w+ ps) Γ(1−w+ ps) Γ(1+ s) sin 1 2π(ps−w) sinπs as (2π)ps ds where c = N − 1 2 and the second expression follows from (5). Upon use of the result |ζ(x + i y)|< ζ(x) when x > 1, we obtain the bound |RN |< (2π)w−1ζ(pN−p(s0− 1 2 )) � a (2π)p �N− 1 2 ∫ ∞ −∞ e−φ t F(t) d t, (8) where φ = arg a and F(t) = � � � � Γ(1−w+ ps) Γ(1+ s) � � � � cosh 1 2πpt coshπt (s = N − 1 2 + i t). On the integration path, F(t) is regular and satisfies F(t) = O(e− 1 2π|t|) as t →±∞. Hence the integral in (8) is convergent and independent of |a| provided |φ|< 1 2π. It then follows that RN = O(aN− 1 2 ) (a→ 0 in |arg a|< 1 2 π). The expansion (7) is the dominant algebraic expansion associated with Sp(a; w) valid as a→ 0 in |arg a| < 1 2π, provided w (> 0) and p are not even integers when the sum in (7) is finite. The same analysis can be applied to the case with non-positive w to yield the Berndt- Ramanujan result [1, Theorem 3.1, p. 306] Sp(a;−w) = 1 p Γ � 1+w p � a−(1+w)/p + N−1 ∑ k=0 (−)k k! ζ(−w− kp) ak +O(aN− 1 2 ) (w≥ 0) as a→ 0 in |arg a|< 1 2π. The reflection formula (5) can be employed to convert the argument of the zeta function to a positive form. R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 8 3. The Expansion of Sp(a; w) When w and p are Even Integers Throughout this section we let w and p be even positive integers, with w = 2m where m = 1,2, . . . . In this case s0 = (2m− 1)/p, which cannot equal an integer and so no double pole can arise. More importantly, there is now only a finite set of poles of the integrand in (4) at s = s0 and s = 0,1,2, . . . , K , where K = ⌊w/p⌋, since the poles of Γ(−s) at s = K + k (k = 1,2, . . .) are cancelled by the trivial zeros of the zeta function ζ(s) at s = −2,−4, . . .. This has the consequence that the integrand is holomorphic in ℜ(s) > max{s0, K}, so that further displacement of the contour can produce no additional algebraic terms in the expansion of Sp(a; w). Thus, we find from (6) upon displacement of the integration path to the right over the poles of the integrand Sp(a; w) = 1 p Γ � 1−w p � a(w−1)/p + K ∑ k=0 (−)k k! ζ(w− pk) ak + (−)m(2π)w IL , (9) where IL = (−)m(2π)−w 2πi ∫ L Γ(−s)ζ(w− ps)asds (10) and L denotes a path parallel to the imaginary axis withℜ(s)> (w/p)+δ, with δ denoting an arbitrary positive quantity. This is easily seen to satisfy the requirement ℜ(s) > max{s0, K0} necessary for the validity of (9). We now employ the functional relation for ζ(s) in (5) to convert the argument of the zeta function in (10) into one with real part greater than unity. The integral in (10) can then be written in the form IL = 1 2πi ∫ L ζ(1−w+ ps) Γ(1−w+ ps) Γ(1+ s) sin 1 2πps sinπs χ−sds, where χ = (2π)pa−1. In [5, §8.1.4], the zeta function appearing in the above integrand was written as an infi- nite series. Here we follow a suggestion made by J. Boersma and retain this function in the integrand; see also [3, §8]. Making use of the expansion (see, for example, [5, p. 368]) sin 1 2πps sinπs = 2 N−1 ∑ r=0 cosπ( 1 2 p−2r−1)s+ ¨ 0 (p/2 even) 1 (p/2 odd) , N = [ 1 4 p], where square brackets denote the nearest integer part∗, we obtain IL = N−1 ∑ r=0 {J+r + J−r }+ ¨ 0 (p/2 even) J (p/2 odd) . (11) ∗The nearest integer part corresponds to [x] = N when x is in the interval (N − 1 2 , N + 1 2 ]. Note that when p = 2, we have N = 0 and the above expansion contains no information. R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 9 Here we have defined the integrals J±r and J by J±r = 1 2πi ∫ L ζ(1−w+ ps) Γ(1−w+ ps) Γ(1+ s) (χe∓πi(p/2−2r−1))−sds (12) and J = 1 2πi ∫ L ζ(1−w+ ps) Γ(1−w+ ps) Γ(1+ s) χ−sds. (13) 3.1. Asymptotic Evaluation of J± r and J The integrals J±r and J have no poles in the half-plane ℜ(s) > (w/p) + δ, so that we can displace the path L as far to the right as we please. On such a displaced path |s| is everywhere large. Let M denote an arbitrary positive integer. The ratio of gamma functions appearing in (12) and (13) may then be expanded by making use of the result (for p > 1) given in [5, p. 53] Γ(1−w+ ps) Γ(1+ s) = A 2π (hκκ)−s §M−1 ∑ j=0 (−) jc jΓ(κs+ ϑ− j) +ρM (s)Γ(κs+ ϑ−M) ª , (14) where c0 = 1, ρM (s) = O(1) as |s| →∞ in |arg s|< π and κ= p− 1, h= p−p, ϑ = 1 2 − w, A= (2π) 1 2 κ 1 2−ϑpϑ. (15) The coefficients c j ≡ c j(w, p) (0 ≤ j ≤ 4) are listed in [5, pp. 46–48] where an algorithm for their determination is described; see Section 4 for details. Substitution of the expansion (14) into the integrals J±r in (12) then produces J±r = A 2π M−1 ∑ j=0 (−) jc j 2πi ∫ L ζ(1−w+ ps)Γ(κs+ ϑ− j) (X e∓πiψr )−κsds+R±M ,r = A 2πκ M−1 ∑ j=0 (X e∓πiψr )ϑ− j (−) jc j 2πi ∫ L′ Γ(u)ζ(qu+λ j)(X e∓πiψr )−udu+R±M ,r . Here we have made the change of variable u → κs + ϑ − j, with L′ denoting the modified integration path, and have defined X := κ(hχ)1/κ, ψr := 1 2 p−2r−1 κ , q := p κ , λ j := 1+ 1 κ (w+ p( j − 1 2 )) (16) together with the remainders R±M ,r = A 4πi ∫ L ρM (s)ζ(1−w+ ps)Γ(κs+ ϑ−M) (X e∓πiψr )−κsds. (17) We note that λ j > 0 for j ≥ 0 when w > 0 and p ≥ 2. The above integrals appearing in J±r may now be evaluated by means of (4), when we replace s by −s and allow the integration R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 10 path (c−∞i, c+∞i) to coincide with the path L′, to yield the sum Sq(X e∓πiψr ;λ j) as defined in (1). This evaluation is valid provided that the variable X e∓πiψr satisfies the convergence condition |arg(X e∓πiψr )|< 1 2π; that is � � � arg a κ ∓πψr � � �< 1 2 π (0≤ r ≤ N − 1). It is routine to verify that these conditions are met when |arg a|< 1 2π. Thus we find J±r = A 2πκ M−1 ∑ j=0 (−) jc j(X e∓πiψr )ϑ− j Sq(X e∓πiψr ;λ j) +R±M ,r . (18) Bounds for the remainders of the type R±M ,r have been considered in [5, p. 71, Lemma 2.7]; see also [2, §10.1]. The integration path in (17) is such that ℜ(1− w+ ps) > 1, so that we may employ the bound |ζ(x + i y)| ≤ ζ(x) for x > 1. A slight modification of Lemma 2.7 in [5, p. 71] then shows that R±M ,r = O � X ϑ−M e−X e∓πiψr � (19) as a→ 0 in the sector |arg a|< 1 2π. An analogous procedure applied to J in (13) shows that J = A 2πκ M−1 ∑ j=0 (−) jc jX ϑ− jSq(X ;λ j) +O(X ϑ−M e−X ) (20) as a→ 0 in |arg a|< 1 2π. 3.2. The Expansion of Sp(a; w) The expansion of IL as a→ 0 in |arg a| < 1 2π then follows from (9), (11), (18) and (20). We obtain the following theorem. Theorem 1. Let m and M be positive integers. Then, when w= 2m and p is also an even positive integer, with K = ⌊w/p⌋ and N = [1 4 p] (with square brackets denoting the nearest integer part), we have the expansion valid as a→ 0 in |arg a|< 1 2π Sp(a; w) = 1 p Γ � 1−w p � a(w−1)/p + K ∑ k=0 (−)k k! ζ(w− pk) ak + (−)m(2π)w IL , (21) with IL = N−1 ∑ r=0 Er(a; w, p) +δpp∗ ÊN (a; w, p), (22) R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 11 where δpp∗ is the Kronecker symbol with p∗ = 4N + 2. The sums Er(a; w, p) are given by† Er(a; w, p) = A 2πκ ∑ ± M−1 ∑ j=0 (−) jc j(X e∓πiψr )ϑ− j Sq(X e∓πiψr ;λ j) + RM ,r (23) for 0 ≤ r ≤ N − 1, where X = κ(h(2π)p/a)1/κ, ψr = ( 1 2 p−2r−1)/κ (0 ≤ r ≤ N − 1), λ j = 1+ (w+ p( j − 1 2))/κ, q = p/κ and the parameters κ, h, ϑ and A are defined in (15). The leading coefficient c0 = 1 and c j ≡ c j(w, p) ( j ≥ 1) are discussed in Section 4. The sum ÊN (a; w, p) is also given by (21) when we put ψN ≡ 0 and omit the summation ∑ ±. The remainders RM ,r satisfy the bound RM ,r = O(max{X ϑ−M e−X e±πiψr }) (0≤ r ≤ N). It is seen from (21), (22) and (23) that the sum Sp(a; w) has been expressed in terms of the sums Sq(X e∓πiψr ;λ j), which involve the reciprocal power of the asymptotic variable a scaling like a−1/κ. Thus, as a → 0 the argument X → ∞. It is obvious from the definition in (1) (when q > 0) that Sq(z;λ j)∼ e−z (z→∞ in |arg z|< 1 2 π), so that the Er(a; w, p) represent a series of exponentially small expansions of increasing sub- dominance in the small-a limit. In addition, the number of exponentially small expansions increases by one each time p increases by 4. By means of a saddle-point analysis in the case w = 0, this was demonstrated to correspond to a Stokes phenomenon when p was allowed to vary continuously through the values p = 2,6,10, . . .; see [5, §§8.1.2, 8.1.7]. Finally, we remark that the exponents p and q are conjugate exponents, since 1 p + 1 q = 1. When a is a real parameter, the expansion in Theorem 1 can be expressed in a different form by using (1) to represent the Sq(X e∓πiψr ;λ j) as infinite sums. Then from (23) we obtain the following theorem: Theorem 2. Let w and p be even positive integers, N = [1 4 p] and M be a positive integer. Then, the exponentially small expansions in (23) valid as a→ 0+ can be written in the form Er(a; w, p) = A πκ ∞ ∑ n=1 nw−1X ϑn e−Xn cosπψr Υn,r (0≤ r ≤ N − 1), (24) ÊN (a; w, p) = A 2πκ ∞ ∑ n=1 nw−1X ϑn e−XnΥn,N . (25) The Υn,r (0≤ r ≤ N) have the asymptotic expansions Υn,r = M−1 ∑ j=0 (−) jc jX − j n cos[Xn sinπψr +π( j − ϑ)ψr] +O(X−M n ), (26) †The symbol ∑ ± signifies that the series with ± signs are to be added. R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 12 where Xn = X np/κ = κ(h(2πn)p/a)1/κ, ψr = ( 1 2 p−2r −1)/κ, ψN ≡ 0 and the other quantities are as defined in Theorem 1. The result in (21) and (22) is the analogue of the Poisson-Jacobi transformation in (2) corresponding to w = 0, p = 2. In this latter case, N = 0 and, from (22), IL = E0(a; 0, 2). The ratio of gamma functions in (14) is replaced by the single gamma function Γ(s + 1 2) by the duplication formula for the gamma function, with the result that c0 = 1, c j = 0 ( j ≥ 1) and consequently Υn,N = 1 for all n ≥ 1. Then, since K = 0 and ζ(0) = −1 2 , (21), (22) and (25) reduce to (2). The resulting expansion is valid for all values of the parameter a (not just a → 0) satisfying |arg a| < 1 2π. When w = 0, p = 2m, the expansions (24) and (25) reduce to those given in [6]; see Section 5. We observe that the n-dependence in the sums Er(a; w, p) from the factor nw−1X ϑn is given by nw−1 npϑ/κ = n−(2w+p−2)/(2κ). Since p ≥ 2 and w> 0, this is seen to correspond to a negative power of n. 4. The Coefficients c j We describe an algorithm for the computation of the coefficients c j ≡ c j(w, p) that appear in the exponentially small expansions Er(a; w, p) in (23) and (26). The expression for the ratio of two gamma functions in (14), with α≡ 1−w for convenience, takes the form Γ(α+ ps) Γ(1+ s)Γ(κs+ ϑ) = A 2π (hκκ)−s §M−1 ∑ j=0 c j (1− κs− ϑ) j + ρM (s) (1− κs− ϑ)M ª , where the parameters κ, h, ϑ and A are defined in (15) and (α) j = Γ(α + j)/Γ(α) is the Pochhammer symbol. If we introduce the scaled gamma function Γ∗(z) = Γ(z)/( p 2π zz− 1 2 e−z), then we have Γ(βs+ γ) = Γ∗(βs+ γ)(2π) 1 2 e−βs(βs)βs+γ− 1 2 e(βs;γ), where e(βs;γ) := exp � (βs+ γ− 1 2 ) log(1+ γ βs )− γ � . The above ratio of gamma functions may therefore be rewritten as R(s)G(s) = M−1 ∑ j=0 c j (1− κs− ϑ) j + ρM (s) (1− κs− ϑ)M (27) as |s| →∞ in |arg s|< π, where R(s) = e(ps;α) e(s; 1)e(κs;ϑ) , G(s) = Γ∗(α+ ps) Γ∗(1+ s)Γ∗(κs+ ϑ) . We now let ξ := (κs)−1 and expand R(s) and G(s) for ξ→ 0 making use of the well-known expansion [5, p. 71] Γ∗(z)∼ ∞ ∑ k=0 (−)kγkz−k (|z| →∞; |arg z|< π), R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 13 where γk are the Stirling coefficients, with γ0 = 1, γ1 = − 1 12 , γ2 = 1 288 , γ3 = 139 51840 , γ4 = − 571 2488320 , . . . . After some straightforward algebra we find that R(s) = 1+ ξ 2 § α(α− 1)κ p − ϑ(ϑ− 1) ª +O(ξ2), G(s) = 1+ ξ 12 � 1− p− 1 p � +O(ξ2), so that upon equating coefficients of ξ in (27)) we can obtain c1. The higher coefficients can be obtained by matching coefficients recursively with the aid of Mathematica to find [5, p. 47] c0 =1, c1 = 1 24p (2− 5p+ 2p2 − 12w+ 12pw+ 12w2), c2 = 1 1152p2 (4+ 28p− 87p2 + 28p3 + 4p4 + 48w− 216pw+ 24p2w+ 144p3w −96w2 − 120pw2 + 480p2w2 − 96w3 + 480pw3 + 144w4), . . . . (28) The rapidly increasing complexity of the coefficients with j ≥ 3 prevents their presentation. However, this procedure is found to work well in specific cases when the various parameters have numerical values, where up to a maximum of 100 coefficients have been so calculated. In Table 1 we present the values‡ of the coefficients c j for 1 ≤ j ≤ 8 in the specific examples considered in Section 5. Table 1: The Coefficients c j (1≤ j ≤ 8) for Different p and w. j p = 4, w= 2 p = 4, w= 4 p = 6, w= 2 p = 6, w= 4 1 1.395833 (0) 3.645833 (0) 1.472222 (0) 3.305556 (0) 2 3.495009 (0) 1.648980 (1) 3.861497 (0) 1.469946 (1) 3 1.230179 (1) 9.075366 (1) 1.380091 (1) 8.081628 (1) 4 5.555372 (1) 5.899040 (2) 6.207979 (1) 5.260968 (2) 5 3.060544 (2) 4.424055 (3) 3.387328 (2) 3.949570 (3) 6 1.990604 (3) 3.760330 (4) 2.188492 (3) 3.358058 (4) 7 1.493190 (4) 3.572267 (5) 1.639364 (4) 3.189927 (5) 8 1.269216 (5) 3.750863 (6) 1.396172 (5) 3.348999 (6) When p = 2, use of the duplication formula shows that the ratio of gamma functions in (14) becomes Γ(1 2− 1 2 w+s)Γ(1− 1 2 w+s) Γ(1+ s) = M−1 ∑ j=0 (−) jc jΓ(s+ 1 2 −w− j) +ρM (s)Γ(s+ 1 2 −w−M). ‡In the tables we write the values as x(y) instead of x × 10y . R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 14 The coefficients c j in this case can be expressed in closed form as [5, p. 53] c j = (1 2 w) j( 1 2 + 1 2 w) j j! = 2−2 j(w)2 j j! (p = 2). (29) Finally, we mention that when w= 0 (corresponding to the Euler-Jacobi series) the coefficients c j are listed for p ≥ 2 and j ≤ 8 in [5, p. 374]. 5. Numerical Results and Concluding Remarks We present some examples of the expansion of Sp(a; w) given in Theorem 1 when p and w= 2m are even integers. For convenience in presentation, we extract the factor e−z from the sum Sq(z;λ j) by writing Sq(z;λ j) = e−z Ŝ2(z;λ j), Ŝq(z;λ j) := ∞ ∑ n=1 e−z(nq−1) nλ j . It follows that, when q > 0, Ŝq(z;λ j) = O(1) as z→∞ in |arg z|< 1 2π. Example 1. In the case p = 2, we have κ = 1, q = 2, N = 0, K = m, p∗ = 2, ψ0 ≡ 0 and X = π2/a. The quantity δpp∗ = 1 so that the exponentially small component of S2(a; 2m) consists of the single term Ê0(a; 2m, 2). From (21), (22), (23) and (29) we therefore find S2(a; 2m)−1 2 Γ � 1− 2m 2 � am− 1 2 − m ∑ k=0 (−)k k! ζ(2m− 2k)ak =(−)m � a π �2m− 1 2 e−π 2/a §M−1 ∑ j=0 (−) j(2m)2 j j! � a 4π2 � j Ŝ2(π 2/a; 2m+ 2 j) +O(aM ) ª (30) as a→ 0 in the sector |arg a| < 1 2π. The expansion in this case has been given in an equivalent form in [6]. Example 2. When p = 4, we have κ= 3, q = 4 3 , N = 1, K = ⌊12 m⌋, p∗ = 6 and X = 3( 1 2 π)4/3 a−1/3, λ j = 1 3 (2m+ 4 j + 1). The quantity δpp∗ = 0 so that there is the single exponentially small term E0(a; 2m, 4) with ψ0 = 1 3 . Then we find the expansion as a→ 0 in |arg a|< 1 2π given by S4(a; 2m)−1 4 Γ � 1− 2m 4 � a(2m−1)/4 − K ∑ k=0 (−)k k! ζ(2m− 4k)ak R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 15 =(−)m � 2a π �(4m−1)/6∑ ± e−X e ∓ 1 3 πi∓ 1 3πiϑ §M−1 ∑ j=0 (−) jc j (X e∓ 1 3πi) j Ŝ 4 3 (X e∓ 1 3πi;λ j) +O(aM/3) ª (31) where the coefficients c j ≡ c j(2m, 4) can be obtained from (28) and [5, p. 47] as c0 =1, c1 = 1 48 (7+ 36m+ 24m2), c2 = 1 4608 (385+ 4392m+ 7104m2 + 3648m3 + 576m4), c3 = 1 663552 (39655+ 1191132m+ 2970936m2 + 2666880m3 + 1080000m4 + 200448m5 + 13824m6), . . . . These coefficients are listed in Table 1 for 1≤ j ≤ 8 when m= 1 and m= 2. Example 3. When p = 6, we have κ= 5, q = 6 5 , N = 1, K = ⌊13 m⌋ and X = 5( 1 3 π)6/5 a−1/5, λ j = 1 5 (2m+ 6 j + 2). In this case p∗ = 6, so that δpp∗ = 1 and there are now two exponentially small expansions E0(a; 2m, 6), with ψ0 = 2 5 , and Ê1(a; 2m, 6). Then, as a → 0 in |arg a| < 1 2π, we have the expansion S6(a; 2m)−1 6 Γ � 1− 2m 6 � a(2m−1)/6 − K ∑ k=0 (−)k k! ζ(2m− 6k)ak =(−)m � 3a π �(4m−1)/10∑ ± e−X e ∓ 2 5πi∓ 2 5πiϑ §M−1 ∑ j=0 (−) jc j (X e∓ 2 5πi) j Ŝ 6 5 (X e∓ 2 5πi;λ j) +O(aM/5) ª + (−)m � 3a π �(4m−1)/10 e−X §M−1 ∑ j=0 (−) jc j X j Ŝ 6 5 (X ;λ j) +O(aM/5) ª . (32) The first few coefficients c j ≡ c j(2m, 6) are c0 =1, c1 = 1 36 (11+ 30m+ 12m2), c2 = 1 2592 (517+ 3840m+ 4116m2 + 1392m3 + 144m4), c3 = 1 1399680 (−22253+ 426550m+ 8181720m2 + 5237640m3 + 1468800m4 + 185760m5 + 8640m6). These coefficients are listed in Table 1 for 1≤ j ≤ 8 when m= 1 and m= 2. R. Paris / Eur. J. Pure Appl. Math, 9 (2016), 3-18 16 Table 2: Values of the Absolute Error in the Computation of Sp(a; w) Defined in (33) Using the Expansions (30) and (31). The Value of the Index j0 Corresponds to Optimal Truncation of the Subdominant Expansion E0(a; w, p). p = 2, w= 2 p = 2, w= 4 a |Sp,w| |Sp,w − E0| j0 |Sp,w| |Sp,w − E0| j0 1.00 8.146(−06) 6.637(−09) 8 6.252(−07) 3.642(−08) 6 0.75 2.031(−07) 8.089(−12) 11 9.296(−09) 4.659(−11) 9 0.50 1.584(−10) 1.260(−17) 18 3.437(−12) 7.635(−17) 16 0.20 5.774(−24) 1.542(−43) 47 2.189(−26) 9.830(−43) 45 0.10 7.667(−46) 1.486(−86) 97 7.506(−49) 9.631(−86) 95 p = 4, w= 2 p = 4, w= 4 a |Sp,w| |Sp,w − E0| j0 |Sp,w| |Sp,w − E0| j0 0.200 3.473(−03) 1.329(−06) 7 3.919(−04) 8.742(−06) 6 0.100 4.863(−04) 2.749(−08) 11 4.805(−05) 4.879(−07) 8 0.050 2.737(−05) 2.156(−10) 14 8.456(−06) 1.420(−09) 11 0.010 4.221(−09) 4.621(−17) 23 7.982(−09) 3.041(−16) 21 0.001 1.064(−14) 1.033(−36) 53 1.876(−16) 6.799(−36) 51 We show the results of numerical calculations to demonstrate the achievable accuracy of the expansion in Theorem 1. We define the difference between Sp(a; w) and the finite algebraic expansion by Sp,w ≡ Sp,w(a) := Sp(a; w)− 1 p Γ � 1− w p � a(w−1)/p − K ∑ k=0 (−)k k! ζ(w− pk) ak. (33) In Table 2 we present the absolute error in the computation of Sp(a; w) for different values of the parameter a in the two cases p = 2 and p = 4, with w= 2 and w= 4 using the expansions given in (30) and (31). The first column in each entry displays the absolute value of Sp,w; that is, the accuracy achievable with just the algebraic expansion and no subdominant exponential terms. The second column shows the absolute error when the single optimally truncated expo- nential expansion E0(a; w, p) (denoted by E0 in the table) is included. The optimal truncation index j0, corresponding to truncation of the exponential expansion E0(a; w, p) at, or near, the least term in magnitude, is indicated in the final column. The situation when there is only a single subdominant exponentially small expansion present is straightforward: this sum is truncated at some suitable point thereby introducing a truncation error. If truncation is optimal, then the resulting error is exponentially more reces- sive than the parent exponential expansion. However, in the case of two, or more, exponential expansions of different degrees of subdominance (corresponding to p ≥ 6) the situation is not so obvious. It is not clear, without further investigation, how the error from the truncated leading exponential series compares with the contribution from the next series. We illustrate this by considering the case p = 6 and w = 2 given in (32). In Table 3 we present the absolute error in the computation of S6(a; 2) as a function of the parameter a. We show, in order, the value of |S6,2| and the absolute error in S6,2− E0(a; 2, 6) when the leading REFERENCES 17 subdominant exponential expansion E0(a; 2, 6) is optimally truncated at index j0. The fourth column gives the absolute error when the first few terms of the second exponential expansion Ê1(a; 2, 6) are included (for brevity in the table these exponential expansions are labelled E0 and E1, and their sum is denoted by E0,1). The final two columns show the values of the least term (including prefactors) in E0(a; 2, 6) at optimal truncation and the values of the leading term ( j = 0) of the sub-subdominant expansion Ê1(a; 2, 6). A cursory inspection of Table 3 shows that for a ≃ 0.1 the leading term of Ê1(a; 2, 6) is less than the minimum term of E0(a; 2, 6) and consequently that inclusion of Ê1(a; 2, 6) cannot improve the accuracy. For a<∼0.01, the reverse is true: the leading terms of Ê1(a; 2, 6) are greater than the minimum term of E0(a; 2, 6) and their inclusion therefore increases the overall accuracy. However, it is clear that in both cases the final accuracy achievable is limited by the optimal truncation of the leading subdominant expansion E0(a; 2, 6). Further improvement in the accuracy would require a hyperasymptotic treatment in order to deal with the divergent tails of E0(a; 2, 6) and Ê1(a; 2, 6). A possible hyperasymptotic scheme for the Euler-Jacobi series with p = 3 and w= 0 has been discussed in [3, §8]. Table 3: Values of the Absolute Error in the Computation of S6(a; 2) Defined by (33) Using the Expansion (32). The Value of the Index j0 Corresponds to Optimal Truncation of the Expansion E0(a; 2, 6). a |S6,2| |S6,2 − E0| j0 |S6,2 − E0,1| Min |E0| E1( j = 0) 1×10−1 2.935(−02) 3.780(−05) 6 −− 9.422(−05) 5.095(−05) 5×10−2 1.617(−03) 3.037(−05) 8 1.200(−05) 1.729(−05) 1.191(−05) 1×10−2 9.512(−04) 1.193(−07) 12 5.339(−08) 1.228(−07) 1.904(−07) 5×10−3 1.292(−03) 1.099(−08) 13 8.713(−09) 9.090(−09) 2.148(−08) 1×10−3 1.604(−04) 3.452(−11) 19 3.483(−12) 3.757(−12) 4.053(−11) 1×10−4 9.894(−07) 8.801(−17) 31 2.230(−19) 3.024(−19) 9.201(−17) 1×10−5 6.209(−10) 1.522(−25) 51 1.963(−30) 1.964(−30) 1.564(−25) Finally we remark that the asymptotics of the alternating version of (1) can be deduced from the result in Theorem 1 by making use of the identity ∞ ∑ n=1 (−)n e−anp nw = 21−wSp(2 pa; w)− Sp(a; w). ACKNOWLEDGEMENTS The author wishes to acknowledge helpful correspondence with J. Boersma who pointed out the conjugate nature of the exponents in Section 3.2. References [1] B.C. Berndt. Ramanujan’s Notebooks, Part II. Springer-Verlag, New York, 1985. [2] B.L.J. Braaksma. Asymptotic expansions and analytic continuations for a class of Barnes integrals, Compositio Mathematica 15, 239–341, 1963. REFERENCES 18 [3] V. Kowalenko, N.E. Frankel, M.L. Glasser, and T. Taucher. Generalised Euler-Jacobi Inver- sion Formula and Asymptotics Beyond All Orders, London Math. Soc. Lecture Notes Series 214, Cambridge University Press, Cambridge, 1995. [4] F.W.J. Olver, D.W. Lozier, R.F. Boisvert, and C.W. Clark (eds.). NIST Handbook of Mathe- matical Functions, Cambridge University Press, Cambridge, 2010. [5] R.B. Paris and D. Kaminski. Asymptotics and Mellin-Barnes Integrals, Cambridge Univer- sity Press, Cambridge, 2001. [6] R.B. Paris. A Poisson-Jacobi-type transformation for the sum ∑∞ n=1 n−2m exp(−an2) for pos- itive integer m, arXiv:1501.00685. http://arxiv.org/abs/1501.00685 2015. [7] E.C. Titchmarsh. Introduction to the Theory of Fourier Integrals, Oxford University Press, Oxford, 1975. [8] E.T. Whittaker and G.N. Watson. Modern Analysis, Cambridge University Press, Cam- bridge, 1952. [9] R.A. Wolf, K.A. Farley, and D.M. Kass. Modeling of the temperature sensitivity of the apatite (U-Th)/He thermochronometer, Chemical Geology, 148, 105–114, 1998.