Asymptotic expansion of the modified exponential integral involving the Mittag-Leffler function

We consider the asymptotic expansion of the generalised exponential integral involving the Mittag-Leffler function introduced recently by Mainardi and Masina [Fract. Calc. Appl. Anal. 21 (2018) 1156–1169]. We extend the definition of this function using the two-parameter Mittag-Leffler function. The expansions of the similarly extended sine and cosine integrals are also discussed. Numerical examples are presented to illustrate the accuracy of each type of expansion obtained. MSC: 30E15, 30E20, 33E20, 34E05


Introduction
The complementary exponential integral Ein(z) is defined by and is an entire function. Its connection with the classical exponential integral E 1 (z) = ∞ z t −1 e −t dt, valid in the cut plane | arg z| < π, is [4, p. 150] Ein(z) = log z + γ + E 1 (z), (1.2) where γ = 0.5772156 . . . is the Euler-Mascheroni constant. In a recent paper, Mainardi and Masina [2] proposed an extension of Ein(z) by replacing the exponential function in (1.1) by the one-parameter Mittag-Leffler function E α (z) = ∞ n=0 z n Γ(αn + 1) (z ∈ C, α > 0), which generalises the exponential function e z . They introduced the function for any α > 0 in the cut plane | arg z| < π (−) n z αn+1 (αn + 1)Γ(αn + α + 1) , (1.3) which when α = 1 reduces to the function Ein(z). A physical application of this function for 0 ≤ α ≤ 1 arises in the study of the creep features of a linear viscoelastic model; see [3] for details. An analogous extension of the generalised sine and cosine integrals was also considered in [2]. Plots of all these functions for α ∈ [0, 1] were given. In this note we consider a slightly more general version of (1.3) based on the two-parameter Mittag-Leffler function given by where β will be taken to be real. Then the extended complementary exponential integral we shall consider is upon replacement of n − 1 by n in the last summation. When β = 1 this reduces to (1.3) so that Ein α,1 (z) = Ein α (z). The asymptotic expansion of this function will be obtained for large complex z with the parameters α, β held fixed. We achieve this by consideration of the asymptotics of a related function using the theory developed for integral functions of hypergeometric type as discussed, for example, in [7, §2.3]. An interesting feature of the expansion of Ein α,β (x) for x → +∞ when α ∈ (0, 1] is the appearance of a logarithmic term whenever α = 1, 1 2 , 1 3 , . . . . Similar expansions are obtained for the extended sine and cosine integrals in Section 4. The paper concludes with the presentation of some numerical results that demonstrate the accuracy of the different expansions obtained.

The asymptotic expansion of a related function for |z| → ∞
To determine the asymptotic expansion of Ein α,β (z) for large complex z with the parameters α and β held fixed, we shall find it convenient to consider the related function defined by The parameter γ > 0, but will be chosen to have two specific values in Sections 3 and 4; namely, γ = 1 and γ = 1 + α. It will be shown that the asymptotic expansion of F (χ) consists of an algebraic and an exponential expansion valid in different sectors of the complex χ-plane. The function F (χ) in (2.1) is a case of the Wright function corresponding to p = q = 2. In (2.2) the parameters α r and β r are real and positive and a r and b r are arbitrary complex numbers. We also assume that the α r and a r are subject to the restriction α r n + a r = 0, −1, −2, . . . (n = 0, 1, 2, . . . ; 1 ≤ r ≤ p) so that no gamma function in the numerator in (2.2) is singular. The parameters associated 1 with p Ψ q (χ) are given by The algebraic expansion of F (χ) is obtained from the Mellin-Barnes integral representation [7, p. 56] where, with −γ/α < c < 0, the integration path lies to the left of the poles of Γ(−s) at s = 0, 1, 2, . . . but to the right of the poles at s = −γ/α and s = −k −1, k = 0, 1, 2, . . . . The upper or lower sign is taken according as arg χ > 0 or arg χ < 0, respectively. It is seen that when α = γ/m, m = 1, 2, . . . the pole at s = −m is double and its residue must be evaluated accordingly. Displacement of the integration path to the left when 0 < α < 2 and evaluation of the residues then produces the algebraic expansion H(χe ∓πi ), where (2.4) and ψ denotes the logarithmic derivative of the gamma function.
The exponential expansion associated with p Ψ q (χ) is given by [6, p. 299], [7, p. 57] where the coefficients A j are those appearing in the inverse factorial expansion The coefficients c j are independent of s and depend only on the parameters p, q, α r , β r , a r and b r . For the function F (χ), we have 1 Empty sums and products are to be interpreted as zero and unity, respectively.
We are in the fortunate position that the normalised coefficients c j in this case can be determined explicitly as c j = (α + β − γ) j . This follows from the well-known (convergent) expansion given in [1], [7, p. 41] to which, in the case of F (χ), the ratio of gamma functions appearing on the left-hand side of (2.6) reduces. Then, with X = χ 1/α we have from (2.5) the exponential expansion associated with F (χ) given by pp. 57-58], we then obtain the asymptotic expansion for |χ| → ∞ when 0 < α < 2 and, when α = 2, (2.10) The upper and lower signs are chosen according as arg χ > 0 or arg χ < 0, respectively. It may be noted that the expansions E(χe ∓2πi ) in (2.10) only become significant in the neighbourhood of arg χ = ±π. When α > 2, the expansion of F (χ) is exponentially large for all values of arg χ (see [7, p. 58]) and accordingly we omit this case as it is unlikely to be of physical interest.

Remark 2.1
The exponential expansion E(χ) in (2.9) continues to hold beyond the sector | arg χ| < 1 2 πα, where it becomes exponentially small in the sectors πα ≤ | arg χ| < 1 2 πα when 0 < α ≤ 1. The rays arg χ = ±πα are Stokes lines, where E(χ) is maximally subdominant relative to the algebraic expansion H(χe ∓πi ). On these rays, E(χ) undergoes a Stokes phenomenon, where the exponentially small expansion "switches off" in a smooth manner as | arg χ| increases [4, §2.11(iv)], with its value to leading order given by 1 2 E(χ); see [5] for a more detailed discussion of this point in the context of the confluent hypergeometric functions. We do not consider exponentially small contributions to F (χ) here, except to briefly mention the situation pertaining to the case α = 1.

The asymptotic expansion of Ein
The asymptotic expansion of Ein α,β (z) can now be constructed from that of F (χ) with the parameter γ = 1. It is sufficient, for real α, β, to consider 0 ≤ arg z ≤ π, since the expansion when arg z < 0 is given by the conjugate value. With χ = −z α = e −πi z α , the exponentially large sector | arg χ| < 1 2 πα becomes | − π + α arg z| < 1 2 πα; that is On the boundaries of this sector the exponential expansion is of an oscillatory character. When 0 < α < 2 3 , we note that the exponentially large sector (3.1) lies outside the sector of interest 0 ≤ arg z ≤ π.

3.1
The cases α = 1 and α = 2 The special case α = 1 deserves further consideration. From (3.2) and (3.7) we have the expansion If β = 1, the asymptotic sum in (3.9) vanishes and for large x. But we have the exact evaluation (compare (1.2)) by [4, (6.12.1)]. The additional sum appearing in (3.11) is exponentially small as x → +∞ and is consequently not accounted for in the result (3.10). From Remark 2.1, it is seen that there are Stokes lines at arg z = ±π(1 − α), which coalesce on the positive real axis when α = 1. In the sense of increasing arg z in the neighbourhood of the positive real axis, the exponential expansion E 1,β (z) is in the process of switching on across arg z = π(1 − α) and 2 E 1,β (z) is in the process of switching off across arg z = −π(1 − α). When α = 1, this produces the exponential contribution for large x. Thus, the more accurate version of (3.9) reads as x → +∞. When β = 1, this correctly reduces to (3.11). When β = 2, we have [8] (−x) n (n + 1) 2 (n + 2)n! = log x − ψ(2) + 1 x This can be seen to agree with (3.12) after a little rearrangement. When α = 2, it follows from (1.4) that so that it is sufficient to consider only the sector 0 ≤ arg z ≤ 1 2 π. From (2.10) we deduce that as |z| → ∞ in 0 ≤ arg z ≤ 1 2 π. Near arg z = 1 2 π, the dominant term is zE 2,β (z) = O(z −β−1 e z ). We note that the term E 2,β (ze πi ) is only significant in the neighbourhood of the positive real axis.
When z = x is real, we therefore have the expansion (3.14)
Theorem 5. When α = 1 and β is real the following expansions hold: and as x → +∞.

Numerical results
In this section we present numerical results confirming the accuracy of the various expansions obtained in this paper. In all cases we have employed optimal truncation (that is truncation at, or near, the least term in modulus) of the exponential and (when appropriate) the algebraic expansions. It is worth pointing out that when β = M α, where M = 1, 2, . . . , the infinite sum over k in the algebraic expansion in (3.2) will reduce to a finite sum consisting of M terms on account of the vanishing of the term 1/Γ(β − αk) for k ≥ M .   Table 1 shows the values 3 of the absolute relative error in the computation of Ein α,β (z) from the asymptotic expansions in Theorem 1 for complex z for values of α in the range 0 < α ≤ 2. It will noticed that there is a sudden reduction in the error when α = 1 and θ = π/4. In this case, the value of θ 0 = 1 2 π and a more accurate treatment would include the exponentially small contribution zE α,β (z). When this term is included we find the absolute relative error equal to 6.935 × 10 −11 . Table 2 presents the error resulting from the expansions for Ein α,β (x) in Theorem 2 for several values of x and different α.    Finally, in Table 3 we present the error associated with the expansions of Sin α,β (x) and Cin α,β (x) given in Theorems 3-6. For Sin α,β (x), the logarithmic expansion arises for α = 1 2 and α = 1 4 ; for Cin α,β (x) this type of expansion arises for α = 1 3 .