Asymptotics of the contour of the stationary phase and efficient evaluation of the Mellin-Barnes integral for the F_3 structure function

A new approximation is proposed for the contour of the stationary phase of the Mellin--Barnes integrals in the case of its finite asymptotic behavior as ${\rm Re} z\to -\infty$. The efficiency of application of the proposed contour and the quadratic approximation to the contour of the stationary phase is compared by the example of the inverse Mellin transform for the structure function $F_3$. It is shown that, although for a small number of terms $N$ in quadrature formulas used to calculate integrals along these contours, the quadratic contour is more efficient, but for $N>20$ the asymptotic stationary phase integration contour gives better accuracy. The case of the $Q^2$-dependence of the $F_3$ structure function is also considered.


I. INTRODUCTION
The Mellin-Barnes (MB) integrals [1] are widely used at high-energy physics. Recently, significant progress has been made in the numerical computation of these integrals. The review of digital packages is presented in the introduction to the paper [2]. The list of physical tasks solved using the MB integrals (two-loop massive Bhabha scattering in QED, in three-loop massless form factors and static potentials, in massive two-loop QCD form factors, in B-physics studies, in hadronic top-quark physics, and for angular integrations in phase-space integrals) can be supplemented by the problem of determining the structure functions and parton distributions in x-space on the basis of their Mellin moments.
The inverse Mellin transform method is widely used in calculations related to deep inelastic scattering (DIS) [3][4][5] for describing the scaling violation in polarized and unpolarized structure functions and fragmentation functions. The general expression for the inverse Mellin transform is written as a contour integral in the complex zplane in the form where the contour C usually runs parallel to the imaginary axis, to the right of the rightmost pole. In case of the DIS, a functionf on the right-hand side of expression (1) is the moments of the structure function at some fixed value of momentum transfer squared Q 2 0 and is usually expressed in terms of the ratio of Γ-functions, like the expression (27) for the Mellin moments of the F 3structure function, which will be investigated in Sec. IV.
Then, the integral (1) can be considered as a typical onedimensional MB integral.
The best efficiency in a numerical integration Eq. (1) can be achieved on the contour of the stationary phase, where the oscillations of the integrand are minimal. However, the solution of the differential equation for the stationary phase contour and its subsequent application to calculate the MB integral requires big computing expenses (see, e.g., Ref. [6]). Instead, it is proposed to build such approximations of the stationary phase contour that would allow the effective application of the quadrature integration formulas [2]. It can be said that the first attempt to construct an effective approximation for the contour using the expansion of the integrand at the saddle point was made by Kosower [7] as applied to the calculation of parton distributions in the x-space. We call this contour C K . Recently, it was proposed to construct the Padé approximation of the stationary phase contour, taking into account the presence of a saddle point in the integrand and its asymptotic behavior for large z [2]. The effectiveness of this approach is shown in the summary Table [2]. However, according to this table, for the integrand F 1 (z, s) = (−s) −z Γ 3 (−z)Γ(1 + z)/Γ(−2z) at s = −1/20 the relative accuracy 10 −8 is achieved by taking into account 16 polynomials in the quadrature formula, whereas for the contour C K , we achieve the same accuracy with 12 terms. In this concrete case, the use of the integration contour C K is more effective than the use of the Padé approximation of the stationary phase contour constructed in Ref. [2]. When does the contour C K , which does not take into account the asymptotic behavior, allows us to obtain a reasonable result? Will the advantages of this contour be preserved with increasing the required accuracy? The answer to these questions will prompt our consideration in this work.
Here, we propose a new approximation for the contour of the stationary phase. The asymptotic behavior of the constructed contour coincides with the contour of the stationary phase as Rez → −∞. We consider a particular case of a finite asymptotic behavior of the contour of the stationary phase in the limit Rez → −∞. The MB integral arising for the F 3 structure function corresponds to this case. First, the construction of the contour will be presented on a simpler example of the MB integral whose value is known exactly. The corresponding integral arising in Feynman diagrams was examined, for example, in Ref. [2], where it was denoted by I 5 (s). Further, we consider in detail the application of the new asymptotic contour in the calculation of the F 3 structure function on the basis of their Mellin moments, and give a comparison with the results of applying the integration contour C K .
It should be stressed that the choice of a contour of integration is dictated by a physical task. For example, in the calculation of massive diagrams it is usually enough to calculate the MB integral once but with a high relative accuracy of ∼ 10 −10 − 10 −16 . In the case of finding the shape of the structure functions, the accuracy of ∼ 10 −4 − 10 −5 is sufficient for the fit of experimental data. However, in the process of fitting, the integral has to be calculated many times (more than 10 6 times), so the computational speed that directly depends on the number of terms in the quadrature formula N is important. It is known that the application of the contour proposed by Kosower is effective for a small number of N terms in the quadrature formula when the nodes of the quadrature formula are located near the saddle point [7,8]. However, this contour considerably moves away from the contour of the stationary phase at large |z|. One can expect that the advantages of the contour coinciding with the asymptotic behavior of the exact contour of the stationary phase as Rez → −∞ should be manifested for sufficiently large values of N . We will try to find out at what value of N ⋆ occurs the "change of the mode", i.e. the use of the contour coinciding with the asymptotic behavior of the contour of a stationary phase becomes more efficient.
The paper is organized as follows. In Sec. II, we start with review of the basic expressions relating to the choice of the integration contour in the MB integrals, according to the method proposed in Ref. [7]. In Sec. III, we explain how a contour that coincides with the asymptotic behavior of the stationary phase contour is constructed. Here, we confine ourselves to integrals with finite asymptotic behavior of the stationary phase contour as Rez → −∞. We present the example of the application of constructed contour in the numerical calculation of well-known integral I(s) whose value is known exactly. The case of the implementation of the MB integrals for evaluation of the F 3 (x B ) structure function, where x B is the Bjorken variable, is given in Sec. IV. We also investigate the applicability of the asymptotic stationary phase integration contour if the Q 2 -dependence of the F 3 structure function is taken into account. Numerical estimates of the relative accuracy of the reconstruction of the functions considered above using different contours are given in Sec. V. Summarizing comments are presented in Sec. VI.
In Appendix A, using the previously considered integral I(s), we perform an additional study for the contour with finite asymptotic as Rez → +∞.

II. CONSTRUCTION OF THE CK INTEGRATION CONTOUR
Let us begin with the review of the main relations for creation of the effective contour proposed by Kosower [7] based on the saddle-point method of the integrand in Eq. (1). We use this contour, which recall that is denoted by C K , for comparison of the efficiency of the numerical evaluation of the integral (1) for different contours.
Following the saddle-point method one can rewrite expression (1) as where F (z) ≡ x −zf (z), C ′ means the contour running from the saddle point c 0 , where F ′ (c 0 ) = 0 (see Ref. [7] for additional details). The complex variable z is parameterized as follows: with the conditions x(t 0 = 0) = c 0 and y(t 0 ) = 0. Then, the expansion of F (z(t)) in a series around the saddle point looks like The requirement that the imaginary part of the function is equal to zero ImF (z(t)) = 0 (5) in any order of the F -expansion (4) determines the exact stationary-phase contour. According to Eq. (5), to order O(t 3 ) the contour is written as where the coefficient . Hence, to this order we may write the integral (2) in the form: with the use of the quadratic integration contour C K given by Eq. (6).
Next, putting t = c 2 √ u , one can introduce a new variable u, where c 2 = 2F (c 0 )/F ′′ (c 0 ), and the integration over u can be presented as where the function H 1 (u) up to order O(u) has the form and the contour integration C K in terms of u read as The function H 1 (u) up to order O(u 2 ) looks like and the integration contour is defined by the expression with the coefficients . We designate this contour by C K . It is obvious that if c 6 = 0, then expression (11) turns into Eq. (9).
Finally, one can see that the prefactor in front of the brackets in Eq. (8) corresponds exactly to the weight function for the generalized Laguerre polynomials L (−1/2) n (u). Therefore, one can expect that the application of the generalized Gauss-Laguerre quadrature formula (see, e. g., Ref. can give a fast numerical evaluation of the integral (8) in the lower orders of approximation. This was the key achievement of the method proposed by Kosower [7]. Indeed, already the first approximation (N = 1) gives relative accuracy of the parton distributions reconstruction about a few percent (see, e. g., Table in Ref. [7]). However, the number of points N required for evaluating of the integral (8) using Eq. (13) with desirable accuracy depends on the integrand function and remains the subject of a numerical study.

III. EXAMPLE FOR CONSTRUCTION OF ASYMPTOTIC STATIONARY PHASE CONTOUR
In order to get a more detailed idea about the construction of the asymptotic stationary phase integration contour, we begin with a simple illustrative example of the MB integral, which can arise in Feynman diagrams. The asymptotics of the stationary phase contour in this case tends to a finite limit, just as in the case which will be considered below for the F 3 structure function.
Let us consider the integral [14] with the fundamental strip δ ∈ (−1, 0) in the region 0 < −s < 4. This integral can be evaluated analytically with result: We present the integrand as where ω = − ln(−s), and using known asymptotic ex- we get the integrand asymptotic behavior Whence from the zero phase condition Im ln [Φ(z)] = 0 , we obtain that arg z = 2[(ω + ln 4) y + 3 2 π sign(y)]/5 and the asymptotic behavior of the zero-phase contour is determined by the equation From this equation it follows that in the asymptotic region as Rez → −∞ the zero-phase contour, which is a contour of stationary phase C st , tends to lines parallel to the real axis Thus, the asymptotic stationary phase contour C as in the complex z-plane takes the form z as (y) = x as (y) + iy .
Using this contour, we can represent the original integral (14) as where the function H 2 (y) is given by and, finally, we have where y j = |y as | 2 (x j +1), x j are the roots of the Legendre polynomials P n (x) with normalization P n (1) = 1, and  Comparing the shape of the contours, we also consider the exact contour of the stationary phase, which can be found from the differential equation with condition x(0) = c 0 . (The designation ' * ' means the complex conjugation.) In the case under consideration, we used x and y, since the solution x = x(y) is uniquely determined. In the general case, the equation for exact stationary-phase contour is written in the parametric form using Eq. (3), see Ref. [2] for more details. Figure 1 shows the shape of the contours for the integral I(s) given by Eq. (14) for a fixed value of −s = 1/20. The solid curve represents the asymptotic stationary phase contour C as determined by Eq. (24), the dotted curve corresponds to the exact contour of the stationary phase C st calculated by Eq. (25), and the dashed curve is the contour C K corresponding to Eq. (6). The dashdot-dotted horizontal lines denote the asymptotic limit of the stationary phase contour determined by Eq. (19). The main contribution to the integral I(s) is given by a region near the saddle point c 0 , where the curves are very close to each other. It can be found that the contours C st and C as are close also in the asymptotic region Rez → −∞, but the contours C K and C K quickly move away from the contour C st (see Ref. [8] for details). Note that for the considered value of s, the coefficient combination c 6 c 4 2 in Eq. (12) is a small negative number and the contour corresponding to Eq. (12) practically coincides with the contour C K . Therefore, we do not plot the contour C K in Fig. 1. The influence of the coefficient c 6 for other region of −s > 4 is shown in the Appendix A.

IV. ASYMPTOTIC STATIONARY PHASE INTEGRATION CONTOUR FOR THE F3 STRUCTURE FUNCTION
The reconstruction of the DIS structure function F 3 based on its Mellin moments can be performed using the inverse Mellin transformation (1), which reduces to the MB integral. For an optimal calculation of this integral we construct the efficient contours in exactly the same way as in the previous section for the integral I(s). It should be noted that the accuracy of calculations of the structure function values can be limited to 4-5 decimal places, which corresponds to the accuracy of the experimental data. At the same time, the method of calculating of the F 3 structure function using the asymptotic stationary phase contour makes it easy to obtain the F 3 values with high accuracy.

A. Contour of integration
Let us consider the following parametrization of the F 3 structure function [15] with values of the parameters found in Ref. [10]. We can write the Mellin moments of the structure function via the Γ-functions: = A Γ(β + 1)Γ(α + z) Γ(α + β + 1 + z) + γ Γ(β + 1)Γ(α + 1 + z) Γ(α + β + 2 + z) and present the structure function in the form Introducing the notation where ω B ≡ − ln(x B ), and using the relation (16) one can get the asymptotic behavior Φ(z) as z → ∞ Calculating the argument of the Φ-function and equating its imaginary part to zero, we arrive at the equation From this equation it follows that as z tends to infinity, the argument z tends to ±π and we get the asymptotic behavior x DIS as (y) = y ctg Hence, in the asymptotic region the contour C st (Rez → −∞) tends to the finite limit As a result, the asymptotics of the stationary phase contour are parallel to the real axis. Note that the r.h.s. of Eq. (33) depends on the Bjorken variable x B and the parameter β which relates to the shape of the structure function at large x B -values. There is no dependence on the parameters α and γ contained in Eq. (26). With the growth of x B and β, the width of the corridor between the asymptotic limits increases. In accordance with expression (24), we obtain that the asymptotic stationary phase contour is determined by the expression z DIS as = y ctg To calculate the F 3 -structure function numerically, we can use the expression (21), which is now read as and also Eqs. (22), (23), and (29), replacing y as ⇒ y DIS as and z as ⇒ z DIS as . Figure 2 shows the shape of a set of contours C as (solid curve), C st (dotted), C K (dashed) and C K (dash-dotted) for a fixed value of x B = 0.1. In the calculation, we used the parameter values in Eq. (26) that were found in Ref. [10] at fixed Q 2 = Q 2 0 = 3 GeV 2 . One can see that, qualitatively, the behavior of the curves in Fig. 2 is similar to the behavior in Fig. 1. There is some difference between the shapes of the contours C st and C as in the vicinity of the saddle point. However, in the limit Rez → −∞ the contours C st and C as coincide, while the contours C K and C K move away from the contour C st in this limit.

B. Q 2 evolution and the contour of integration
Let us discuss the change of the behavior of the contour C as with a change of the momentum transfer squared Q 2 .
The perturbative Q 2 -evolution of the Mellin moments is well known (see, e.g., Refs. [4,11]), and in the nonsinglet case in the leading order (LO) is given by the formula Here γ is the Euler constant, ψ is the usual logarithmic derivative of the Γ-function, n f is the number of active flavors, and α LO s is the LO running coupling. Finding the asymptotic behavior of the integrand Φ(z, Q 2 ) ≡ x −z B M 3 (z, Q 2 ) as |z| → ∞, analogously to how it was made before for a fixed of value Q 2 , and also taking into account an asymptotic behavior of the ψ function, we obtain Φ(z, Q 2 ) ∼ exp ω B z − (β + 1) + ∆(Q 2 ) ln z . Next, from the zero phase condition follows that ω B y − (β + 1) + ∆(Q 2 ) arg z = 0 . Therefore, the asymptotic behavior of the contour as Rez → −∞ is determined by the formula In the asymptotic region the contour C as tends to the finite limiting value Thus, the account of the Q 2 -evolution of the Mellin moments of the F 3 structure function comes down to replacement in the expressions without evolution β ⇒ β + ∆(Q 2 ) and using Eq. (35) x DIS as (y) ⇒ x DIS as (y, Q 2 ) and y DIS as ⇒ y DIS as (Q 2 ). It is important to note that the expression defining the contour C as in higher orders of the perturbation theory QCD will coincides with the LO expression with the replacement only of α LO s (Q 2 ) in Eq. (37) by an expression for the running coupling in the corresponding order of the perturbation theory. Figure 3 shows the influence of the Q 2 evolution on a values of the saddle point c 0 and the 'asymptotic' point c as = x as (y = 0) in the case of the F 3 structure function. In this figure c 0 and c as are shown as functions of the Bjorken variable x B at Q 2 = Q 2 0 = 3 GeV 2 (solid and dashed lines for c 0 and c as , respectively) and Q 2 = 100 GeV 2 (dotted and dash-dotted lines). One can see that there is no significant difference in the behavior of the point c 0 with changing of Q 2 as the solid and dotted curves are close to each other. The same behavior is observed also for the point c as . It is interesting to note a sharp increase in the values of c 0 and c as if x B > 0.5, as well as the convergence of all curves at large values of x B .

V. NUMERICAL ESTIMATIONS OF ACCURACY
In this section, by numerically calculating the MB integrals given by Eqs. (14) and (28), we investigate the question at what values of polynomials N in the corresponding quadrature formula the advantages of the asymptotic contour C as are manifested, in comparison with the quadratic approximation C K to the contour of the stationary phase.
The relative accuracy of a reconstruction is defined as usual ε(N ) = |(f N − f exact )/f exact | , where f N is the sum given by Eq. (13), when the contour C K is used, or Eq. (23) for the contour C as ; f exact is the exact value of the corresponding integral.
The relative accuracy ε(N ) of calculating the I(s) integral depending on number of terms N in the sums (13) and (23) is presented in Table I. The result is given for −s = 1/20 and −s = 2.0. The calculations for the F 3 structure function are presented in Table II for the values x B = 0.01, 0.1, and 0.5. In doing so, we use the values of the parameters in Eq. (26) that were found in Ref. [10] at fixed Q 2 0 = 3 GeV 2 . Tables I and II show similar results for the relative accuracy ε(N ). The contour C K works more efficiently for a small the number of terms (N < 20) in the quadrature formulae (13) and (23), as for a large N > 20, the result of the asymptotic contour C as is more accurate by 2-3 orders of magnitude. In addition, as can be seen from Table II, the high accuracy achieved at N * = 20 exceeds the experimental accuracy by several orders of magnitude. However, as can be seen from Table I, the "regime change" for the integral I(−s = 2.0) occurs at value of ε(N * ) ≃ 4 × 10 −6 . Such accuracy has practical importance in the QFT calculations. Figure 4 shows the effect of the contour Q 2 evolution according to Eq. (39) for the F 3 structure function. In this figure the relative accuracy ε(N ) is plotted as a function of number terms N in the Gauss-Legendre quadrature formula (23) at values of the Bjorken variable x B = 0.01 (triangles) and x B = 0.8 (squares) for the contours C as (x B , Q 2 0 ) (open triangles and squares) and C as (x B , Q 2 ) (full triangles and squares) for Q 2 = Q 2 0 = 3 GeV 2 and Q 2 = 100 GeV 2 . For the running coupling α LO s (Q 2 ) we used a value Λ QCD = 363 MeV and n f = 4  (see Ref. [10]). One can see that the contour C as (x B , Q 2 ) yields a more exact result; however, this advantage is compensated when using the contour C as (Q 2 0 ) if the number of terms N is increased by 2-4 units.

VI. CONCLUSIONS
We presented a new approximation for the construction of a contour close to the contour of the stationary phase. A special case of the finite asymptotics of the stationary-phase contour in the limit Rez → −∞ was considered as such an asymptotics arises in the calculation of the MB integral that represents the F 3 structure function in terms of its Mellin moments. The proposed approximation reproduces the behavior of the contour zero phase for large values of |z| and has the form like Eq. (34) in the case of the DIS structure functions.
It was compared the efficiency of application of the asymptotic stationary phase contour C as and the contour C K , determined by the saddle-point method, as described in Sec. II, in calculating of the MB integrals (14) and (28). As expected, the contour C K turned out to be more effective for a small number of N terms in the Gauss-Laguerre quadrature formula, when the nodes of the quadrature formulas are located near the saddle point [7,8]. The advantage of the contour C as is manifested at large values of N . The "regime change" occurs at N * ≈ 20. For the region 20 < N 35, the contour C as gives a relative error of 2-3 orders of magnitude better than the contour C K . With increasing N , this gap increases too.
When the Q 2 -evolution of the structure function F 3 is taken into account, the efficiency of the contours C as (Q 2 0 ) and C as (Q 2 ) was compared. Although the contour C as (Q 2 ) gives a more accurate result, but this advantage is compensated by using the contour C as (Q 2 0 ) if we increase the number of terms in the quadrature formula only by 2-3 units. The contour C as (Q 2 0 ) can be considered as a universal one, that is, applicable for other values of Q 2 .
It should be emphasized that since the necessary accuracy is achieved when using a small number of N in the quadrature formulas (13) or (21), the computer time for calculating the F 3 structure function by using Eq. (28) with the efficient contours is significantly less than if one uses the linear contours that are usually parallel to the imaginary axis, to the right of the rightmost pole in the integrand, or a straight line at an angle (see, e.g., Refs. [4,5,12,13]).
Here we have considered the F 3 structure function, which is the simplest among the DIS structure functions since it does not contain the contribution of gluons and sea quarks. The parameterization of the shape of this structure function (26) is typical and widely used in the QCD analysis of the structure functions. Our consideration can be applied to the polarized nonsinglet combination ∆q 3 , ∆q 8 and the nonsinglet combination fragmentation function D π + uv . The choice of the efficient contour for the singlet case can be performed along the same line, but requires more complicated formulas. This is the task for forthcoming investigation.
Our result can be useful in choosing an integration method in both the one-dimensional case and the case of few-dimensional MB integrals if it is required to achieve relative accuracy unattainable in the integration by the Monte Carlo method. This case is interesting for us, because the exact contour of zero phase C st , see Eq. (25), go away from the imaginary axis to the right. The asymptotic behavior of the contour C st is defined by Eq. (18), as before, however, from this equation it follows that in the asymptotic region as Rez → +∞ the contour tends to another limit value y as = 3 2 π |ω + ln 4| sign(y) .
Therefore, the numerical value of the integral I(s) can be found using the Gauss-Legendre quadrature formula (23) in which y as is given by Eq. (A2) and the contour C as by Eqs. (18) and (24).
The contours C K and C K are defined by the same way as before, and the numerical value of the integral I(s) is determined using the Gauss-Laguerre quadrature formula (13). Figure A shows the shape of contours for fixed −s = 5. The solid curve corresponds to the contour C as , the dotted curve is the exact contour C st , the dashed curve is the contour C K , and the dash-dotted is the contour C K . From this figure it is clear that in calculating the integral I(s) using the contour C K , difficulties will arise, as this contour has an incorrect direction. Only the first few terms in the sum (13) may give a reasonable result, since the contour C K is close to the contour C st near the saddle point c 0 . Thus, for the case −s > 4 the contour C K it is impossible to consider as efficient. At the same time, the contour C K works rather good, but applying it to achieve relative accuracy, for example, 10 −12 , can be problematic. The proposed contour C as , whose a construction is quite simple, reproduces the behavior of the exact contour C st well, and the use of this contour can provide the required high relative accuracy.