The QCD Static Energy using Optimal Renormalization and Asymptotic Pad\'e-approximant Methods

The perturbative QCD static potential and ultrasoft contributions, which together give the static energy, have been calculated to three- and four-loop order respectively, by several authors. Using the renormalization group, and Pad\'e approximants, we estimate the four-loop corrections to the static energy. We also employ the optimal renormalization method and resum the logarithms of the perturbative series in order to reduce sensitivity to the renormalization scale in momentum space. This is the first application of the method to results at these orders. The convergence behaviour of the perturbative series is also improved in position space using the Restricted Fourier Transform scheme. Using optimal renormalization, we have extracted the value of $\Lambda^{\overline{\textrm{MS}}}_{\textrm{QCD}}$ at different scales for two active flavours by matching to the static energy from lattice QCD simulations.


I. INTRODUCTION
The static potential energy of quantum chromodynamics (QCD) is the non-Abelian analogue of the well-known Coulomb potential energy of quantum electrodynamics. The short distance part of this quantity is calculated in non-relativistic QCD (NRQCD) [1,2] framework and involves the evaluation of Feynman diagrams. It has been studied extensively in recent years and analytical results are known to three-loop. At four-loop order contributions involving the ultrasoft gluons that start to contribute from three-loop order in perturbation series are known. Such contributions only appear for short distances when system is weakly coupled. They are calculated using weakly coupled potential non-relativistic QCD (pNRQCD) [3,4] formalism. Hence, we discuss only weakly coupled regime of pNRQCD in this article. The static potential depends on the renormalization scale µ, the ultrasoft factorization scale µ us , and magnitude of the three momentum transfer between the heavy sources p (= |p|). In this article MS renormalization scheme is used.
The first attempt to perform the full three-loop (numerical) calculations was by two independent groups [5,6]. The results were found to be in agreement. Analytical calculations were presented almost six years later in ref [7]. At three-loop order, ultrasoft gluons also appear which are capable of changing singlet to octet state of the system and vice-versa. Such contributions were first pointed in ref [8], and were calculated in ref [9] in pN-RQCD. The renormalization group (RG) improvement of the ultrasoft terms at three-loop order was first discussed in ref [10]. The next order calculations for the ultrasoft terms can be found in ref [11] and their resummation is discussed in ref [12]. * anant@iisc.ac.in † diganta99@gmail.com ‡ alam.khan1909@gmail.com QCD is known to be non-perturbative at long distance, and this fact is manifest in lattice QCD (LQCD) simulations. The static energy, that includes the static potential and the ultrasoft contributions, can also be computed in LQCD simulations, and hence it becomes a topic of great interest to extract various parameters of the theory. Some recent LQCD simulation results for the static energy can be found in refs [13][14][15][16][17][18], and they support the Cornell potential type behaviour for the heavy quarkantiquark system.
In this article, we have addressed the following issues: • At four-loop order, the constant term contributing to the perturbative series is as yet unknown. It requires calculation of four-loop Feynman diagram. In the absence of such a calculations we use RGE and Padé approximants [33] to get the estimate for the four-loop coefficients.
• The RG improvement of the static energy by resummation of all RG-accessible running logarithms following the method advocated in refs [34][35][36][37][38][39][40][41], and for the first time applied in this paper to the ultrasoft logarithms present in the static energy. All the leading and next-to-leading logarithms at each order in perturbation theory that can be accessed through the RG equation (RGE) are called RGaccessible logarithms. We call this renormalization group improved perturbative series as RG-summed or optimal renormalized series. It should be noted that our resummation scheme discuss about RGimprovement with respect to renormalization scale whereas ref [10,12] deals with the RG-improvement with respect to ultrasoft scale.
• The convergence of the static energy is improved to the four-loop order in position space using the arXiv:2007.10775v2 [hep-ph] 16 Sep 2023 Restricted Fourier Transform (RFT) proposed in ref [42].
• We use the RG-summed and unsummed forms of the static energy in momentum space to extract Λ MS QCD to four-loop from the LQCD inputs from ref [15].
• In comparison to the unsummed series, the RGsummed static energy in momentum space gives better fit to the static energy obtained from LQCD in ref [15].
• The four-loop RG-summed static energy is used as a trial case to show that these improvements (scale sensitivity and better fit) also persist at the higher order.
The scheme of this paper is as follows: In section II, we discuss the perturbative treatment to the static energy and various pieces associated with this quantity in the weak coupling limit. In section III, we calculate some contributions to the static energy at the four-loop using the RGE. In section IV, we use Padé approximants to estimate all the four-loop coefficients even the one that is not accessible with the RGE. Some of the estimates are found in agreement with RGE solutions. In section V, we perform all order RG-summation of certain running logarithms and show that this RG-improvement will bring down the sensitivity to the renormalization scale. In section VI, we discuss the improvement of the static energy to the four-loop order in position space by removing the pathological uncontrolled contributions using RFT. The inputs from the previous sections are applied in section VII to fit the LQCD inputs for the static energy to extract the Λ MS QCD for two active flavours in momentum space. A discussion is presented in section VIII and we summarize our results in section IX. Appendix A contains the QCD beta function coefficients. Appendix B contains the formula used in for running of the strong coupling constant with momentum and also in terms of Λ MS QCD . Appendix C contains the known contribution to the static energy. Appendix D contains the useful formula for calculating the restricted and unrestricted versions of position space static potential and the static energy. Appendix E contains the final result of uncontrolled contribution to the static energy in position space to the four-loop order.

II. THE PERTURBATIVE QCD-STATIC ENERGY
A heavy quark and anti-quark system, with heavy quark mass m Q and relative velocity v < 1, is nonrelativistic in nature and various scales [4,43] present in the system are hard scale ∼ O(m Q ), soft scale ∼ O(m Q v), and ultrasoft scale ∼ O(m Q v 2 ). If these scales are well separated then we can integrate them one by one to study only relevant degrees of freedom. Integrating out the hard scale from QCD gives the non-relativistic QCD (NRQCD) in which the soft and the ultrasoft degrees of freedom are dynamical. Contributions to the static energy at different orders were calculated using this framework. The one-loop perturbative calculations for massless quarks were first performed in the late 1970s and can be found in refs [8,[44][45][46] and in the massive case in ref [47]. Two-loop massless calculations appeared in refs [48][49][50][51] while the massive case were considered in refs [29,[52][53][54].
The pNRQCD is obtained from NRQCD by integrating out the soft scale and this formalism is best suitable for studying the threshold systems. The heavy quark and anti-quark systems in this formalism are described by color singlet fields S and color octet fields O. The gauge fields are multipole expanded about the interquark separation and the gauge invariant Lagrangian for pNRQCD [4,19], at leading order in 1/m Q is given by where L light is Lagrangian for light quarks, , O] and E are chromoelectric field strength. V s (r, µ us ) and V o (r, µ us ) are singlet and octet potential at leading order in r while V A and V B appear at sub-leading order in r. It should be noted that these potential will also depend on the renormalization scale µ for any finite order calculation. The potentials appearing at higher order in 1/m Q are hidden in ellipses which disappear in static limit m Q → ∞. The gauge fields are already multipole expanded about inter-quark separation (R) in eq(1) and therefore, F µνa ≡ F µνa (R, t).
In the weak coupling regime, the system has the hierarchy of scales: and in this limit, the static energy can be written as a perturbative series in the strong coupling constant α s . The singlet static energy to known orders can be written as: where V s (p, µ, µ us ) is the perturbative static potential and δ us (p, µ, µ us ) are contribution from ultrasoft gluons. The static potential encodes the interaction of the quarks and gluons degrees of freedom in the singlet state and is known to three loop. It takes the following form: In the above, L ≡ log µ 2 p 2 is running logarithm, and C F = 4/3 is color factor of the SU (3) representation. The expansion parameter in the above equation is defined as x ≡ (α s (µ)/π), and any argument of x indicates the scale at which it is evaluated. The coefficients T i,0,0 are the i th −loop perturbative contributions. The oneloop coefficient T 1,0,0 can be found in refs [46,47], the two-loop coefficient T 2,0,0 in refs [49][50][51] and the threeloop coefficient T 3,0,0 in refs [5][6][7]. All these coefficients are collected in appendix C. The coefficients of infrared logarithms, T us i,j,k , in the static potential and can be found in refs [9,11].
In effective field theory language, the static potential is a matching coefficient which depends upon the factorization scale µ us . The presence of the logarithmic terms in the static potential at three-loop were first pointed out in the ref [8] and they act as a source of infrared divergences. The ultrasoft part δ us is now known to next to the leading order (NLO) [9,11] but it contributes to the static energy from three-loop order(briefly discussed in section VI). They carry the information of the dynamical ultrasoft gluon degrees of freedom with the ultraviolet cut-off µ us . This scale acts as a source of ultra-violet divergences for these gluons. Both divergences, however, cancel with each other for the static energy which results in non-analytic dependence (∼ α n s log m (α s )) in terms of the expansion parameter and the total energy in eq(3) takes the form: where T i,0,0 i>2 = T pert i,0,0 + δT us i,0,0 and constant terms (δT us i,0,0 ) are the contributions from the ultrasoft gluons and can be found in ref [11].
The ultrasoft contributions to the static energy are now known to the four-loop order, but the perturbative correction to static potential are yet to be calculated at this order. Some of the contributions to the static energy at this order can be accessed using RGE and they are discussed in the next section III.

III. RG SOLUTIONS OF THE FOUR-LOOP CONTRIBUTIONS
We can rewrite eq(5) in terms of the coupling at the renormalization scale µ using eq(B1) as: Despite the explicit dependence, all-orders series should be independent of the renormalization scale µ. Mathematically, this implies that any perturbative series W (x, L) must satisfy the RGE: where The QCD beta function β is defined in terms of x as: where the coefficients β i up to five-loop order [55][56][57][58][59][60][61][62][63] are given in appendix A. The RGE can be used to solve iteratively for the RG-accessible terms in terms of the QCD β i coefficients and the lower order RG-inaccessible coefficients T i,0,k . The RGE for the W (x, L) along with known results to three-loop and QCD beta functions allow us to extract the RG-accessible four-loop coefficients T 4,1,0 , T 4,2,0 , T 4,3,0 , T 4,1,1 and T 4,4,0 . To obtain them, we note that Coefficients of x i L j log k (x) of the above equation gives the RG-accessible coefficients The coefficients T i,0,k can not be obtained using RGE and are known as the RG-inaccessible terms. The calculation of such terms involves the evaluation of all the Feynman diagrams relevant for that order. Such calculation are yet to be performed so we use asymptotic Padé approximant (APAP) to estimate unknown coefficient, T 4,0,0 , which is discussed in the next section.

IV. PADÉ ESTIMATE OF THE FOUR-LOOP CONTRIBUTION
Padé approximants are rational functions that can be used to estimate the higher order terms of a series from its lower order coefficients. Both the original series and the Padé approximants have the same Taylor expansion to a given order, the next term is taken as its prediction. They are used in refs [37,[64][65][66][67][68][69][70] in the past to improve the higher order perturbative results for QCD. It is worth to mention that they were first used for the static potential in ref [33] and we are extending the results using the same method for static energy to the four-loop order. The procedure used is explained in this section.
Note that the ultrasoft corrections to the static energy are already known to the four-loop order, so we do not need to predict them using Padé approximants. They can be added to the predicted values at the end of calculations. We rewrite the perturbative series in eq(6), without the ultrasoft corrections as W (x, L), in the following form: In general, if the series coefficients {R 1 , R 2 , R 3 , · · · , R N } are known then the Padé approximant for this series is denoted by W [N −M |M ] , and given as: can be used to estimate R N +1 . For example, if only the NLO term of eq(12), R 1 , is known then the next coefficient, R 2 , is estimated from W [0|1] by Taylor expanding it for small x as i.e., R 2 1 is the Padé approximant prediction for R 2 . For the static energy, the coefficients R 1 , R 2 , R 3 are known and can be read off by comparing the x i -th terms in eq(6) with eq (12). To predict the unknown four-loop coefficient R 4 in eq(12), we use the approximant W [2|1] and its series expansion, for small x, is given by: Now, first we solve for A 1 , A 2 and B 1 in terms of known R 1 , R 2 and R 3 of eq(12) and the solutions are: In the next step, the coefficient of x 4 in eq (16) is taken as the prediction for R 4 which, in terms of the lower order R i 's, can be written as: In large L limit, we get the the following form for R pred 4 compatible with eq(13) as: These predictions can be further improved by knowing the asymptotic behaviour of the series, and this extra piece of information is taken as input to get a more precise approximation. This improvement is termed as APAP and can be found in ref [68].
The error associated with such approximation in asymptotic limit [33,68] is given by: where D = N + M (1 + a p ) + b p and A, a p , b p are fitting parameters. We get APAP results in terms of R pred 4 by: Repeating the procedure for known lower order R ′ i s, we can fix the constants A, a p and b p for a fixed M . It is worth to mention that among different choice of Padé approximant for APAP, for a given order, M = 1 and a p = 0, b p = 0 gives best result compatible with RG for the static energy and hence this particular choice is used in this article.
Following the procedure explained above, we get the four-loop Padé prediction in the large L limit as: We can see that for T 4,4,0 and T 4,4,3 the predictions from Padé and the renormalization group are in perfect agreement. For the other RG-accessible coefficient T 4,1,0 and T 4,2,0 , the predictions are different. However, numerical difference for T 4,2,0 is not more than 2.2% for active quark flavours n f ≤ 6. However, T 4,1,0 has larger deviations ≥ 2% for n f > 2 from RGE prediction. For this reason, we will restrict our discussion in next sections only to two active flavours for the four-loop. The Padé prediction for the unknown constant term at the four-loop order for the static potential can be obtained by setting T i,j,k For RG-inaccessible coefficient In fact, similar pattern has been observed for known lower orders: Now we have an estimate for T 4,0,0 but the truncated perturbation series also suffers from scale dependence. The scale sensitivity of the perturbative series can be minimized using RG-summation of running logarithms and the procedure is discussed in the next section.

V. RG IMPROVEMENT IN THE MOMENTUM SPACE
The issue with the perturbative series in the QCD is to account for the RG running of all the parameters. The optimal renormalization method advocated in refs [38,39] accounts for the RG running by summation of all the RG-accessible logarithms. The RG-accessible logarithms at each order in the perturbation theory are defined as the leading and the next-to-leading logarithms that can be accessed through the processes-dependent the RGE. Resummation becomes interesting from the threeloop order due to presence of the ultrasoft terms and this issue is discussed for static energy for the first time in this article.
The perturbative series in question is where the series coefficients are T i,j,k . To obtain RGsummed perturbation series which we call W where intermediate quantities S i,k (xL) are the resummed series obtained by summing terms: We substitute eq (28) in eq(7) which leads to a recursion relation between the series coefficients. We multiply the recursion relation with (xL) k−1 with appropriate k and sum it from n = k to infinity, which following eq(30), give differential equations for S i,k (xL). The solution to these differential equation results into the closed form expression for S i,k (xL). The RG-summed solutions S i,0 (xL) are calculated to the two-loop order in ref [38] and are given below: where w = (1 − β 0 u) andB i = β i /β 0 and L w ≡ log(w).
The RG-summation for e + e − process to three-loop order is also discussed in ref [39] which is a special case of the results of this article in the limit where ultrasoft coefficients are taken zero. The static energy requires a new series representation, compatible with the RGE, to incorporate the ultrasoft logarithms. These logarithms form a separate recurrence relation among the coefficients. The RG-summation of ultrasoft terms at the three-loop order is obtained by collecting the coefficients of x n L n−4 log(x) in eq(10) which results in the following recurrence relation among the coefficients: (n − 3)T n,n−3,1 − nβ 0 T n−1,n−3,1 = 0 .
Collecting x n L n−4 terms, we get the following recurrence relation: Notice that presence of the β 0 T n−1,n−4,1 terms in the above equation is new and differs from e + e − case in ref [39].
Collecting x n L n−5 terms we get the following recurrence relation: (n − 4)T n,n−4,0 − (n − 4)β 4 T n−5,n−5,0 − (n − 3) × β 3 T n−4,n−5,0 − (n − 2)β 2 T n−3,n−5,0 − (n − 1) × β 1 T n−2,n−5,0 − β 1 T n−2,n−5,1 − nβ 0 T n−1,n−5,0 − β 0 T n−1,n−5,1 = 0 , multiplying u k−1 , where k = n − 3 for three-loop and k = n−4 for the four-loop, to the recurrence relations and then summing n from k to ∞ gives separate differential equations for the above recurrence relations. All these recursion relations can be written in a general-differential equation for the static energy as: Solution to different RG-summed series to the order we are interested in are given by: Similarly, other solutions are: The importance of resummation of all accessible logarithms can be seen in the figure(1). The scale dependence of the RG-summed series eq(29) around momen- while the unsummed series in eq(6) has significant µ dependence. Scale sensitivity of unsummed series decreases order by order but the advantage of RG-summed series provides results less sensitive to scale with the same available information. This theoretical improvement provide us an opportunity to extract various parameters from available experimental data with less scale sensitivity.

VI. RESTRICTED FOURIER TRANSFORM
Let us begin by noting that the ultrasoft part of the static energy is calculated in position space, whereas the perturbative part is carried out in momentum space. To make any phenomenological study, it is necessary to bring all the contributions of the static energy into the same space. Bringing the momentum space results in position space via a Fourier transform may seem natural, but such a transformation induces pathological contributions to the static energy. These undesirable contributions originate from the non-perturbative, small momentum modes, which must be removed explicitly. After removing such contributions, the convergence behaviour of the static energy improves drastically. This scheme was first addressed in ref [42] and was termed as the Restricted Fourier transform (RFT). This scheme has been already discussed in detail for the static energy to the three-loop. Here we discuss it very briefly and provide corresponding results for the four-loop order.
Note that the ultrasoft terms are already known to the four-loop order. Therefore, what we are providing is the complete calculation of the static energy in the position and momentum space to the four-loop order. The final result for the uncontrolled contribution to the position space static energy in RFT scheme is given in appendix E but an overview of the calculation is given here.
The position space version of the static potential, V (r, µ, µ us ), from momentum space potential, saỹ V (p, µ, µ us ) to distinguish between the two basis, in RFT scheme is defined by: where µ f is a perturbative scale chosen such that µ > µ f ≫ Λ MS QCD and uncontrolled terms, δV (r, µ, µ us ; µ f ), in the potential is (43) The static energy, given by eq(3), in RFT scheme is: The ultrasoft gluonic contribution to the static energy, at order r 2 in multipole expansion, is given by: where T F = 1/2, N c is number of colors, ϕ(t, 0) ad j ab is Wilson line in the adjoint representation connecting two points at temporal separation t, E a/b is choromoelectric field strength, and V A is matching coefficient which appear at order r 2 in multipole expansion. This quantity has been calculated to NLO in refs [9,11] and final result, sub-leading in r, is given by: where δ us 3−loop = 2 log and Coefficients C 1 , C 2 and D are given in ref [11] and we have presented them in appendix C. Next to the leading order expression for V o (r) − V s (r) is given by: and we can see that if we substitute eq (49) to eq(46) then ultrasoft contributions to the three-and four-loop order can be obtained from NLO results. Due to the presence of the extra scale in the problem µ us , it is desirable to expand all the quantities in terms of α s (µ). This expansion induces mixed logarithms at the four-loop of form 1 r log(µ 2 us r 2 ) log(µ 2 r 2 ) in the ultrasoft contributions. To simplify the calculation we can write: and log r 2 µ 2 us = log µ 2 r 2 − log so that the ultrasoft scale will appear only in running logarithms. The uncontrolled term for (V o − V s ) (r) in RFT scheme, in next to the leading order is given by: where L µ f = log( µ 2 µ 2 f ) and H ′ i s functions defined by eq(D7) and are proportional to generalized hypergeometric function. The only matching coefficient left is V A and it is taken as same as calculated in ref [42]. These contributions will enter in eq (46) and the final expression is given in appendix E. Padé estimates and the ultrasoft contributions give us numerical estimate for T 4,0,0 = T Padé 4,0,0 + δT us 4,0,0 . Now, the restricted version for the static potential, the ultrasoft contributions and the static energy in position space can be constructed using eq(D1) and eq(D2).
Using Λ MS QCD = 315 MeV, we can see in figure(2) that the four-loop corrections to the static energy in RFT scheme makes a very small contribution but, without removing pathological contributions, behaviour of the same quantity in unrestricted scheme has very bad for r > 0.05 fm at the four-loop.
In the next section, we have performed the extraction of Λ MS QCD in momentum space using LQCD inputs from ref [15]. In this section, we fit the perturbative static energy to the lattice data and extract the value of the Λ MS QCD from RG-summed and unsummed case. Due to absence of all order results, the extracted quantity Λ MS QCD depends on the choice of renormalization scale. To reduce the scale dependence, we exploit the optimal renormalized perturbative static energy to extract Λ MS QCD . The parametrization of lattice data to the Cornell potential is given in position space in ref [15] for two-flavour QCD where α = 0.326 ± 0.005 and σ = 7.52 ± 0.55fm −2 are constants for n f = 2 flavour. The parameter V 0 is a potential offset needed to match to the perturbative static energy. However, V 0 is not needed in case of momentum space analysis [15]. It is important to note that these parameters are correlated and have cor(α, σ) = −0.17 which has to be taken account in sampling from normal distribution. The Fourier transform of lattice static energy, in eq(52), to momentum space is given by: To extract Λ MS QCD , we minimize the square-deviation [15]: It should be noted that the strategy for sampling is analogous to the one discussed in ref [15]. The matching region is chosen 1500 MeV ≤ p ≤ 3000 MeV and momentum values are randomly sampled from a uniform distribution with p min ⊂ [1500, 2250] MeV and p max ⊂ [2250, 3000] MeV such that p max − p min ≥ 375 MeV. It is important to note that five-loop running of α s is used to extract Λ MS QCD . Following the procedure explained above, the scale dependence of the extracted (for 500 samples for {p, α, σ}) Λ MS QCD can be seen from figure (3) and the corresponding results are presented in Table II.
The RG-summed and unsummed series give very similar values of Λ MS QCD in case of renormalization scale is chosen in the middle of the matching region. However, the RG-summed static energy provides better fit to the lattice static energy than the unsummed one which can be seen in Table III below   The RG-summed static energy gives better fit and this improvement persist even at next higher order. For example, at the four-loop order, we have found that the RG-summed static energy produces significantly better fit to lattice static energy compared to unsummed perturbative static energy for different cases when we choose T 4,0,0 according to eq (55).

VIII. DISCUSSIONS
The static energy between a heavy quark and antiquark is an important quantity in QCD as it can be calculated in both perturbation theory and in LQCD simulations. Existence of a matching region, where both perturbation theory and lattice agree, provides an opportunity to extract the parameters of the theory. If we assume the light quarks are massless and the heavier one decouples then the only parameter left in the theory is strong coupling constant. Precise calculation of the static energy thus becomes very important goal for precision physics involving the strong interactions.
The perturbative static energy involves the computation of all the Feynman diagrams at that order and the computer codes were used in refs [5,6] for the three-loop numerical calculation. At this order, virtual emission of the ultrasoft gluons, with energy and momentum smaller than 1/r, capable for changing singlet state to octet state and vice versa also appear. Such emissions induce infrared divergence in the static potential.
The gluon fields are multipole expanded about interquark separation and hence their energy and momentum are restricted to ultrasoft scale. This scale acts as source of ultraviolet divergence for gluonic contributions. However the static energy, which is sum of perturbative potential and contribution of ultrasoft gluons, is divergence free quantity as two divergences get cancelled with each other. This cancellation induces non-analytic terms(∼ α n s (1/r) log m (α s (1/r))) [8] in the static energy. The ultrasoft terms [9] and their resummation is known to NLO [10,12] which contributes at the the three-and four-loop order but the perturbative singlet potential is known only to the three-loop [5][6][7].

A. Padé estimate
In this article, we have extended the results for static energy to the four-loop order using the renormalization group and Padé approximant. Among the seven coefficients at the four-loop, the estimates for T 4,3,0 and T 4,4,0 , using APAP, are in exact agreement with solutions of renormalization group. Deviation for coefficient T 4,2,0 is within 2.2% for n f ≤ 6 but T 4,1,0 deviates more than 2% for n f > 2. Since, the deviation is less than 2% for RGaccessible coefficient for light flavour, we have used APAP estimate for T 4,0,0 to extract the Λ MS QCD to the four-loop order with different choices given in eq (55). The large-n f limits for RG-accessible terms from both methods are in perfect agreement. There are also some contribution to T 4,0,0 from ultrasoft terms which demands this quantity must be Fourier transformed in order to get complete the RG-inaccessible terms in momentum space.

B. Position space improvement
The static energy from LQCD simulations are mostly parametrized in position space and hence the perturbative static energy is Fourier transformed to the position space. This quantity in position space suffers from pathological contributions stemming from the nonperturbative regions and has to be removed explicitly. This is achieved using the Restricted Fourier transform advocated in ref [42] which improves the convergence behaviour for r ∼ 0.12 fm. The ultrasoft and the static potential has explicit dependence on another scale µ us which is absent in total energy. This scale-dependence should also be cancelled for the static energy in the RFT scheme. Final expression for uncontrolled contributions to the static energy is provided in appendix E. The fourloop contribution to static energy in RFT scheme has very little effect and can be seen in figure (2). The static energy in RFT scheme provided in the section VI for the four-loop order and can be used in future studies. The RG-summed static energy in momentum space is used in this article to extract Λ MS QCD by fitting the static energy from perturbative to the static energy from LQCD from ref [15] for two active flavour. Its value for the threeloop from the RG-summed static energy is found to be 308.4±10.6(2.25 GeV) MeV, 313.7±10.9(4.17 GeV)) MeV and 317.0 ± 11.1(6.5 GeV) MeV where quantity in the parenthesis is the renormalization scale. For the unsummed static energy, this parameter is found to be 307.3 ± 10.9 MeV. The RG-summed version of static energy has been observed to provide not only the better fit to the the lattice energy but also giving less standard deviation if the renormalization scale is chosen in the middle of matching region. Similar trend also persist for next order but we have used less sample size since the exact calculations are not available. Our finding of Λ MS QCD from RG-summed and unsummed series agrees within error bars to the findings of ref [15].

IX. SUMMARY
To summarize, the QCD static potential is known to three-loop order, and the ultrasoft terms which first appear at the three-loop order are known to four-loop order. In section II, we describe the perturbative and the ultrasoft part of the static static energy. The main results of this paper are the following • In section III, using the RGE we determine the RGaccessible coefficients at four-loop order which is shown in eq(11).
• The constant term of the four-loop coefficient can not be determined using RGE. In section IV, we use the Padé approximant method to obtain this term and is given in equation eq (24).
• In section V, we apply for the first time the method of optimal renormalization to QCD static energy beyond two-loop order to sum up the RG-accessible running logarithms to all order in perturbation theory. The RG-summed series is defined in equation eq (29) and the subsequent quantities are given in eq (31)(32), eq (39)(40)(41)(42). The RG-summed series ensures the expected reduction in sensitivity to the renormalization scale as shown in figure(1).
• In section VI, we use the Restricted Fourier Transform scheme to improve the convergence behaviour of the static energy in the position space to fourloop order.
• Using the RG-summed series eq(29) (see definition in eq (6)) and the lattice QCD parametrization eq(53), we fit the the QCD scale Λ MS QCD to the lattice data. Our fit results at different loop orders can be found Table I-II and scale dependence at these orders in figure 3.
• The uncertainties associated with our extraction of the Λ MS QCD is discussed in section VIII. In summary, we have used a variety of techniques, theoretical and numerical, and have rendered the picture of the QCD static energy as a very useful tool to obtain a clear handle on Λ MS QCD which is one of the fundamental parameters of the QCD, there by confirming the results in a large number of other studies. We have also studied the consistency of the picture by invoking Padé approximants as well as renormalization group summation in order to achieve these ends. We also provide improvement in position space using RFT-scheme. Our findings in this article provide better control over variation of renormalization scale for finite order results available for the static energy. It also discusses scale dependence of the extracted Λ MS QCD at different orders of perturbation theory for the first time in this article as an application of the method.

X. ACKNOWLEDGMENTS
BA thanks D. Wyler and the University of Zurich, Switzerland for hospitality when this work started. BA was partly supported by the Mysore Sales International Ltd. Chair of the Division of Physical and Mathematical Sciences, Indian Institute of Science. DD would like to thank the DST, Government of India for the INSPIRE Faculty Award (grant no IFA16-PH170). AK is support by a fellowship from the Ministry of Human Resources Development, Government of India. We thank R. Sarkar for collaboration at an early stage of this investigation. In addition, we are grateful to N. Brambilla, P. Hegde, F. Karbstein, P. Lamba, H. Takaura and A. Vairo for invaluable discussions. We also thank S. Dey and S. Banik for help with numerical simulations. We are particularly grateful to N. Brambilla for reading the manuscript and making numerous suggestions that have improved it greatly.

Appendix A: The QCD β-functions
The QCD beta function is given by where L = log(µ 2 /p 2 ). The couplant used to extract the Λ MS QCD at scale µ is given by: where b i ≡ βi β0 , ℓ ≡ log(log(µ 2 /(Λ MS QCD ) 2 )) and y ≡