On applications of Mathematica Package"FAPT"in QCD

We consider computational problems in the framework of nonpower Analityc Perturbation Theory and Fractional Analytic Perturbation Theory that are the generalization of the standard QCD perturbation theory. The singularity-free, finite couplings ${\cal A}_{\nu}(Q^2),{\mathfrak A}_{\nu}(s)$ appear in these approaches as analytic images of the standard QCD coupling powers $\alpha_s^{\nu}(Q^2)$ in the Euclidean and Minkowski domains, respectively. We provide a package"FAPT"based on the system Mathematica for QCD calculations of the images ${\mathcal A}_{\nu}(Q^2)$, ${\mathfrak A}_{\nu}(s)$ up to N$^3$LO of renormalization group evolution. Application of these approaches to Bjorken sum rule analysis and $Q^2$-evolution of higher twist $\mu_4^{p-n}$ is considered.


Introduction
The QCD perturbation theory (PT) in the region of space-like momentum transfer Q 2 = −q 2 > 0 is based on expansions in a series in powers of the running coupling α s (µ 2 = Q 2 ) which in the one-loop approximation is given by α (1) s (Q 2 ) = (4π/b 0 )/L with b 0 being the first coefficient of the QCD beta function, L = ln(Q 2 /Λ 2 ), and Λ is the QCD scale. The one-loop solution α (1) s (Q 2 ) has a pole singularity at L = 0 called the Landau pole. The ℓ-loop solution α (ℓ) s (Q 2 ) of the renormalization group (RG) equation has an ℓ-root singularity of the type L −1/ℓ at L = 0, which produces the pole as well in the ℓ-order term d ℓ α ℓ s (Q 2 ). This prevents the application of perturbative QCD in the low-momentum space-like regime, Q 2 ∼ Λ 2 , with the effect that hadronic quantities, calculated at the partonic level in terms of a power-series expansion in α s (Q 2 ), are not everywhere well defined.
In 1997, Shirkov and Solovtsov discovered couplings A 1 (Q 2 ) free of unphysical singularities in the Euclidean region [1], and Milton and Solovtsov discovered couplings A 1 (s) in the Minkowski region [2]. Due to the absence of singularities of these couplings, it is suggested to use this systematic approach, called Analytic Perturbation Theory (APT), for all Q 2 and s. The APT yields a sensible description of hadronic quantities in QCD (see reviews [3][4][5]), though there are alternative approaches to the singularity of effective charge in QCD -in particular, with respect to the deep infrared region Q 2 < Λ 2 . One of the main advantages of the APT analysis is much faster convergence of the APT nonpower series as compared with the standard PT power series (see [6]). Recently, the analytic and numerical methods, necessary to perform calculations in two-and three-loop approximations, were developed [7][8][9]. The APT approach was applied to calculate properties of a number of hadronic processes, including the width of the inclusive τ lepton decay to hadrons [10][11][12][13][14], the scheme and renormalization-scale dependencies in the Bjorken [15,16] and Gross-Llewellyn Smith [17] sum rules, the width of Υ meson decay to hadrons [18], meson spectrum [19], etc.
The generalization of APT for the fractional powers of an effective charge was done in [20,21] and called the Fractional Analytic Perturbation Theory (FAPT). The important advantage of FAPT in this case is that the perturbative results start to be less dependent on the factorization scale. This reminds the results obtained with the APT and applied to the analysis of the pion form factor in the O(α 2 s ) approximation, where the results also almost cease to depend on the choice of the renormalization scheme and its scale (for a detailed review see [22] and references therein). The process of the Higgs boson decay into a bb pair of quarks was studied within a FAPT-type framework in the Minkowski region at the one-loop level in [23] and within the FAPT at the three-loop level in [21]. The results on the resummation of nonpower-series expansions of the Adler function of scalar D S and a vector D V correlators within the FAPT were presented in [24]. The interplay between higher orders of the perturbative QCD expansion and higher-twist contributions in the analysis of recent Jefferson Lab data on the lowest moment of the spindependent proton structure function, Γ p 1 (Q 2 ), was studied in [25] using both the standard PT and APT/FAPT. The FAPT technique was also applied to analyse the structure function F 2 (x) behavior at small values of x [26,27] and calculate binding energies and masses of quarkonia [28]. All these successful applications of APT/FAPT necessitate to have a reliable mathematical tool for extending the scope of these approaches. In this paper, we present the theoretical background which is necessary for the running of A ν [L] and A ν [L] in the framework of APT and its fractional generalization, FAPT, and which is collected in the easy-to-use Mathematica package "FAPT" [29]. This task has been partially realized for APT as the Maple package QCDMAPT in [30] and as the Fortran package QCDMAPT F in [31]. We have organized "FAPT" in the same manner as the well-known package "RunDec" [32]. A few examples of APT and FAPT applications are given.

Theoretical framework
Let us start with the standard definitions used in "FAPT" for standard PT calculations. The QCD running coupling, α s (µ 2 ) = α s [L] with L = ln[µ 2 /Λ 2 ], is defined through where n f is the number of active flavours. The β-function coefficients are given by (see [33]) ζ is Riemann's zeta function. We introduce the following notation: ( Then Eq. (1) in the l-loop approximation can be rewritten as: In the one- with the Landau pole singularity at L → 0. In the two-loop (ℓ = 2) approximation where W −1 [z] is the appropriate branch of the Lambert function.
The three-(c k (n f ) = b k (n f ) = 0 for all k ≥ 3) and higher-loop solutions a (ℓ) [L; n f ] can be expanded in powers of the two-loop one, a (2) [L; n f ], as has been suggested and investigated in [8,9,14]: The coefficients C (ℓ) n can be evaluated recursively. As has been shown in [9], this expansion has a finite radius of convergence, which appears to be sufficiently large for all values of n f of practical interest. Note here that this method of expressing the higher-ℓ-loop coupling in powers of the two-loop one is equivalent to the 't Hooft scheme, where one puts by hand all coefficients of the β-function, except b 0 and b 1 , equal to zero and effectively takes into account all higher coefficients b i by redefining perturbative coefficients d i (see for more detail [35]).
The basic objects in the Analytic approach are the analytic couplings in the Euclidian A It is convenient to use the following representation for spectral functions: In the one-loop approximation the corresponding functions have the simplest form whereas at the two-loop order they have a more complicated form with W 1 [z] being the appropriate branch of Lambert function.
In the three-(ℓ = 3) and four-loop (ℓ = 4) approximation we use Eq. (7) and then obtain Here we do not show explicitly the n f dependence of the corresponding quantities -it goes inside through R (2) The package "FAPT" performs the calculations of the basic required objects: α

APT and FAPT applications
As an example of the APT application, we present the Bjorken sum rule (BSR) analysis (see for more details [39]). The BSR claims that the difference between the proton and neutron structure functions integrated over all possible values of the Bjorken variable x in the limit of large momentum squared of the exchanged virtual photon at Q 2 → ∞ is equal to g A /6, where the nucleon axial charge g A = 1.2701 ± 0.0025 [33]. Commonly, one represents the Bjorken integral in Eq. (16) as a sum of perturbative and higher twist contributions The perturbative QCD correction ∆ Bj (Q 2 ) has a form of the power series in the QCD running coupling α s (Q 2 ). At the up-to-date four-loop level in the massless case in the modified minimal subtraction (MS) scheme, for three active flavors, n f = 3, it looks like [36] ∆ The perturbative representation (18) violates analytic properties due to the unphysical singularities of α s (Q 2 ). To resolve the issue, we apply APT. In particular, the four-loop APT expansion for the perturbative part ∆ PT Bj (Q 2 ) is given by the formal replacement Clearly, at low Q 2 a value of α s is quite large, questioning the convergence of perturbative QCD series (18). The qualitative resemblance of the coefficients pattern to the factorial growth did not escape our attention although more definite statements, if possible, would require much more efforts. This observation allows one to estimate the value of α s ∼ 1/3 providing a similar magnitude of three-and four-loop contributions to the BSR. To test that, we present in Figs  dominant contribution to the pQCD correction ∆ Bj (Q 2 ) comes from the four-loop term ∼ α 4 s . Moreover, its relative contribution increases with decreasing Q 2 . In the region Q 2 > 2 GeV 2 the situation changes -the major contribution comes from one-and two-loop orders there. Analogous curves for the APT series given by Eq. (19) are presented in Fig. 2. Figures 1 and 2 demonstrate the essential difference between the PT and APT cases, namely, the APT expansion obeys much better convergence than the PT one. In the APT case, the higher order contributions are stable at all Q 2 values, and the one-loop contribution gives about 70 %, two-loop -20 %, three-loop -not exceeds 5%, and four-loop -up to 1 %.
One can see that the four-loop PT correction becomes equal to the three-loop one at Q 2 = 2 GeV 2 and noticeably overestimates it (note that the slopes of these contributions are quite close in the relatively wide Q 2 region) for Q 2 ∼ 1 GeV 2 which may be considered as an extra argument supporting an asymptotic character of the PT series in this region. In the APT case, the contribution of the higher loop corrections is not so large as in the PT one. The four-loop order in APT can be important, in principle, if the theoretical accuracy to better than 1 % will be required. Now we briefly discuss how the APT applications affect the values of the higher-twist coefficients µ p−n 2i in Eq. (17) extracted from Jlab data. Previously, a detailed higher-twist analysis of the four-loop expansions in powers of α s was performed in [39]. In Figs. 3 and 4 we present the results of 1-and 3-parametric fits in various orders of the PT and APT. The corresponding fit results for higher twist terms µ p−n 2i , extracted in different orders of the PT and APT, are given in Table 1 (all numerical results are normalized to the corresponding powers of the nucleon mass M ). From these figures and Table 1 one can see that APT allows one  to move down up to Q 2 ∼ 0.1 GeV 2 in description of the experimental data [39]. At the same time, in the framework of the standard PT the lower border shifts up to higher Q 2 scales when increasing the order of the PT expansion. This is caused by extra unphysical singularities in the higher-loop strong coupling. It should be noted that the magnitude of µ p−n 4 /M 2 decreases with an order of the PT and becomes compatible to zero at the four-loop level. It is interesting to mention that a similar decreasing effect has been found in the analysis of the experimental data for the neutrino-nucleon DIS structure function xF 3 [37] and for the charged lepton-nucleon DIS structure function F 2 [38].
Consider the application of the FAPT approach by the example of the RG-evolution of the non-singlet higher-twist µ p−n 4 (Q 2 ) in Eq. (17). The evolution of the higher-twist terms µ p−n 6,8, ... is still unknown. The RG-evolution of µ p−n 4 (Q 2 ) in the standard PT reads In the framework of FAPT the corresponding expression reads as follows: We present in Table 2 the best fits for µ p−n 4 (Q 2 0 ) taking into account the corresponding RGevolution with Q 2 0 = 1 GeV 2 as a normalization point and without the RG-evolution. We do not Table 2. Results of higher twist extraction from the JLab data on BSR with inclusion and without inclusion of the RG-evolution of µ p−n 4 (Q 2 ) normalized at Q 2 0 = 1 GeV 2 . for the standard PT calculations and compare with FAPT since the only effect of that would be the enhancement of the Landau singularities by extra divergencies at Q 2 ∼ Λ 2 , whereas at higher Q 2 ∼ 1 GeV 2 the evolution is negligible with respect to other uncertainties. We see from Table 2 that the fit results become more stable with respect to Q min variations, which reduces the theoretical uncertainty of the BSR analysis.

Summary
To summarize, APT and FAPT are the closed theoretical schemes without unphysical singularities and additional phenomenological parameters which allow one to combine RGinvariance, Q 2 -analyticity, compatibility with linear integral transformations and essentially incorporate nonperturbative structures. The APT provides a natural way for the coupling constant and related quantities. These properties of the coupling constant are the universal loop-independent infrared limit and weak dependence on the number of loops. At the same time, FAPT provides an effective tool to apply the Analytic approach for RG improved perturbative amplitudes. This approaches are used in many applications. In particular, in this paper we consider the application of APT and FAPT to the RG-evolution of nonsinglet structure functions and Bjorken sum rule higher-twist analysis at the scale Q 2 ∼ Λ 2 considered.
The singularity-free, finite couplings A ν (Q 2 ), A ν (s) appear in APT/FAPT as analytic images of the standard QCD coupling powers α ν s (Q 2 ) in the Euclidean and Minkowski domains, respectively. In this paper, we presented the theoretical background, used in a package "FAPT" [29] based on the system Mathematica for QCD calculations in the framework of APT/FAPT, which are needed to compute these couplings up to N 3 LO of the RG running. We hope that this will expand the use of these approaches.