The theory predictions of high-energy hadronic collider observables rely on the factorization theorem in Quantum Chromodynamics (QCD), which states that the hadronic cross section is factorized into universal parton distribution functions (PDFs) and partonic cross sections. The partonic cross sections are perturbatively calculable. The PDFs are non-perturbative quantities, but their scale evolution (the well-known DGLAP evolution [1][2][3]) is governed by splitting functions, which can be evaluated perturbatively in QCD.
Several benchmark partonic cross sections in QCD have been evaluated to next-to-next-to-next-to-leading order (N 3 LO), see for example [4][5][6][7]. To achieve the same accuracy for hadronic cross sections, it is necessary to know the N 3 LO PDFs, which require the knowledge of four-loop splitting functions. The splitting functions at three-loop accuracy in QCD were computed almost 20 years ago [8,9], and allowed the complete determination of NNLO PDFs. At four-loop order, results are available only for some specific color structures, including the leading-power 𝑁 𝑓 contributions to all channels [10][11][12], the 𝑁 2 𝑓 contribution to non-singlet splitting functions [12], leading color contribution to non-singlet splitting functions [13], and recently the 𝑁 2 𝑓 contributions to pure-singlet [14] and quark-to-gluon splitting functions [15]. Beyond these leading color, leading and sub-leading 𝑁 𝑓 contributions, a finite number of Mellin moments were computed for all splitting functions in [13,[16][17][18][19]. Those results were already used to obtain approximate N 3 LO PDFs [20,21].
The anomalous dimensions 𝛾(𝑛) with Mellin moments 𝑛 are related to splitting functions 𝑃 (𝑥) via the following Mellin transformation,
To go beyond the currently available all-𝑛 results, in this paper, we consider one of the simplest remaining contributions: the 𝑁 𝑓 𝐶 3 𝐹 contribution to the non-singlet splitting functions, that also appear in Quantum Electrodynamics (QED). Following closely references [22,14], we performed our computations in the framework of the operator product expansion (OPE), and extracted the splitting functions from the single pole of off-shell operator matrix elements (OMEs). The off-shell OMEs are defined as the off-shell matrix elements with an operator insertion, for the case of two partons in the external states it is T. Gehrmann, A. von Manteuffel, V. Sotnikov et al.
where 𝑂 𝑖 is a twist-two operator. In the current context, it is the following quark non-singlet operator,
(3)
In the above equation, 𝜆 𝑘 ∕2 denotes diagonal generators of the flavor group SU(𝑁 𝑓 ), and Δ is a light-like reference vector with Δ 2 = 0. The symbol 𝜓 represents the quark field, and 𝐷 𝜇 𝑖𝑗 = 𝜕 𝜇 𝛿 𝑖𝑗 -𝑖𝑔 𝑠 𝑇 𝑎 𝑖𝑗 𝐴 𝑎 𝜇 is the covariant derivative in the fundamental representation of a general gauge group. Compared with the conceptually complicated renormalization [23][24][25][26]22] in the singlet sector, the non-singlet sector is much easier and allows for a straight-forward multiplicative renormalization:
where superscripts B and R are used to represent the bare and renormalized operators, respectively. The renormalized operator satisfies the following renormalization group equation,
which defines the anomalous dimension 𝛾 ns of the non-singlet, twist-two operator. From equation ( 5) and the fact of the bare operator 𝑂 B ns does not depend on the scale 𝜇, it is easy to see that
It is not difficult to solve the above equation order by order in 𝑎 𝑠 = 𝛼 𝑠 ∕(4𝜋) with the help of the 𝑑-dimensional QCD 𝛽 function
where 𝜖 = (4 -𝑑)∕2. To four-loop order, the explicit result is given by
ns -2𝛽 0 𝛾 (1) ns + 3𝛾 (0) ns 𝛾 (1) ns ]
0 𝛾 (1) ns -14𝛽 0 𝛾 (0) ns 𝛾 (1) ns + 12𝛽 0 𝛽 1 𝛾 (0) ns + 6𝛾 (1) ns (𝛾 (0) (2) ns -6𝛽 1 𝛾 (1) ns -6𝛽 2 𝛾 (0) ns + 8𝛾 (0) ns 𝛾 (2) ns + 3(𝛾 (1) ns ) 2 ] + 6𝛾 (3) ns 𝜖 )
where 𝛾 (𝑙) ns is defined as
Therefore, the four-loop, non-singlet anomalous dimension 𝛾 (3) ns can be determined from the single pole in 𝜖 of the renormalization constant 𝑍 ns . By separating the even and odd moments, 𝛾 ns can be decomposed as
where 𝛾 s ns represents the flavor singlet but valence non-singlet contribution.
To extract the 𝑍 ns order by order in 𝑎 𝑠 , we insert equation ( 4) between two off-shell external quark states,
Here we also need to consider the renormalization of the wave function, the strong-coupling constant, and the gauge parameter, all for 𝜉 = 1 in Feynman gauge. In addition to the explicit expressions for various contributions to 𝑍 𝑞 , 𝑍 𝑔 and 𝑍 𝑎 𝑠 that were collected in the appendix of the reference [14], we need only one more contribution to these renormalization constants, the 𝑁 𝑓 𝐶 3 𝐹 part of 𝑍 (4) 𝑞 , which we document in appendix A.
In the following, we compute the 𝑁 𝑓 𝐶 3 𝐹 contribution to the four-loop corrections to the off-shell OME ⟨𝑞|𝑂 B ns |𝑞⟩ B .
We generated the relevant Feynman diagrams with an insertion of the operator 𝑂 ns by QGRAF [27]; some sample diagrams can be found in Fig. 1. The required Feynman rules for the operator 𝑂 ns exhibit some peculiar patterns, i.e., terms like (Δ ⋅ 𝑝) 𝑛-1 appear and thus prevent the application of standard integration-by-parts (IBP) [28][29][30] algorithms in moment 𝑛-space. A method first proposed in [31,32] was used to overcome this difficulty, by summing these peculiar terms into linear propagators using a tracing parameter 𝑡. As an example,
After the desired manipulations have been performed, one can reexpand in 𝑡 to obtain the result for some moment 𝑛. This method allows the applications of standard IBP algorithms and has been widely used to study the matching coefficient of heavy flavor quark contributions in deepinelastic scattering [32][33][34][35] and splitting function calculations [36,37,22] from off-shell OMEs.
We translated the Feynman rules in 𝑛-space to parameter 𝑡-space and then worked in 𝑡-space throughout. For this calculation, we used Mathematica to substitute the Feynman rules in 𝑡-space into the Feynman diagrams. FORM [38] was used to evaluate the traces of Dirac and color matrices.
Regarding topology classification, we first used Apart [39] and MultivariateApart [40] (see also [41][42][43]) to eliminate the linear dependence among Feynman propagators. Then we classified all resulting topologies into integral families with 18 propagators by an in-house code. The code searches for the possible loop momentum transformations to see if one topology can be mapped into another one or not. To reduce the size of the unreduced amplitude, we first employed Reduze2 [44] to eliminate integrals from zero sectors, and then applied shift relations to relate integrals between different sectors.
The IBP reductions for the amplitude were done by the private code Finred, which employs finite field sampling and rational reconstruction techniques [45][46][47]. It is well-established that optimizing the input IBP system can significantly enhance the efficiency of the reduction process. In our study, we achieved this optimization by utilizing the linear algebra method [48] in order to exert control over the generation of squared propagators.
We used the differential equation (DE) method [49] to determine the solutions of the master integrals. The derivation of the system of differential equations for the master integrals follows the same approach as used for the amplitude reduction. In the first step, we select master integrals according to our generic integral ordering, and chose to ignore IBP relations generated from seeds in supersectors even if this leads to missing linear relations between the master integrals. The rational functions in the differential equations are therefore somewhat complicated and involve rational numbers with many digits. Consequently, for their reconstruction, we employed a significant number of samples for the variables and several finite fields, each with a cardinality of order (2 63 ). We chose to reveal these missing "hidden" relations by exploring the so-called scaling relations (see e.g. [50]) of the master integrals. In our case, the scaling relations read
where 𝐼 𝑖 represents the 𝑖-th master integral and [𝐼 𝑖 ] denotes mass dimension of 𝐼 𝑖 , and we set [𝑡] = 0, [Δ] = -1 such that both Δ ⋅ 𝑝 and 𝑡 are dimensionless. In practice, through IBP reductions, the left-hand side of ( 13) can not always be reduced to the simple form on the right-hand side. By enforcing the above scaling relations, we obtained several extra relations among the master integrals. Those additional relations greatly simplified the DE system. In the current case, the 𝑁 𝑓 𝐶 3
𝐹 contribution, we found 658 remaining master integrals, and the corresponding DE system with respect to 𝑡 can be cast into 𝜖-form [51] by applying the codes CANONICA [52,53] and Libra [54,55]. We obtained
where we had set 𝑝 2 = -1 and Δ ⋅ 𝑝 = 1. ⃗ 𝐼 is the vector of the new canonical master integrals, 𝐀 (𝑖) are matrices involving rational numbers only, and 𝑡 𝑖 = 0, ±1, 2. Interestingly, in addition to the letters appearing in harmonic polylogarithms (HPLs) [56], a new letter 𝑡 -2 appears. This new letter enters the solutions of canonical master integrals starting from transcendental weight 7 only and does not contribute to the 𝑁 𝑓 𝐶 3 𝐹 contributions to the non-singlet splitting functions. It would be curious to see if the new letter contributes, e.g., to the finite part of corresponding off-shell OMEs or not. We leave it to future study.
The boundary conditions can be conveniently chosen in the limit 𝑡 → 0, where the linear propagators trivialize, and additional relations between the master integrals allow for their further reduction. The resulting master integrals are four-loop self-energy integrals [57,58], which in the present case were mapped to the master integrals for two-point functions in [59,60]. By mapping to self-energy master integrals, we were able to determine the boundary conditions for all 658 master integrals. In practice, it is easier to first apply the regularity conditions: no branch cuts can be generated
in the Taylor series of equation (12). We thus solved the canonical differential equations in the limit 𝑡 → 0 by keeping 𝜖 to all orders, and we set 𝑐 𝑎 to zero for terms 𝑐 𝑎 𝑡 ±𝑎𝜖 (𝑎 is a positive integer) in the resulting solutions. In this way, we expanded the canonical solutions to transcendental weight 7 in terms of HPLs and generalized polylogarithms (GPLs) with the letter 𝑡 -2.
The amplitude reduction was performed directly in terms of the canonical basis, where we can use the anticipated factorization of the 𝜖 and 𝑡 dependence in the denominators, and construct the denominators first [61,40]. This helps to reduce the number of numerical samples required to reconstruct the functional dependence in 𝜖 and 𝑡. We subsequently inserted the solutions of the canonical basis into the amplitude and expanded the resulting amplitude order by order in 𝜖. We observed the emergence of Harmonic Polylogarithms (HPLs) with weights up to 6 in the single pole of the amplitude. For this 𝜖-expanded amplitude, we reconstructed also the rational numbers from their images in various finite fields. It is expected that the rational numbers appearing in the 𝜖-expanded amplitude are simpler, thus fewer finite fields are required for their reconstruction.
In this manner, we expressed the result for the 𝑁 𝑓 𝐶 3 𝐹 contribution to the bare four-loop OME, denoted as ⟨𝑞|𝑂 B ns |𝑞⟩ B , in terms of Harmonic Polylogarithms (HPLs) in parameter-𝑡 space. Subsequently, we transformed this expression to 𝑛-space using harmonic sums [62,63], aided by the
HarmonicSums package [64][65][66][67][68][69]. This transformation yielded the 𝑛-space representation for the bare OME.
In the previous section, we obtained the 𝑁 𝑓 𝐶 3 𝐹 contribution to four-loop bare OME to the single pole in 𝜖 in 𝑛-space. The constant 𝑍 ns can be readily extracted from the renormalization equation (11). Notice that we also need the 3-loop corrections to the bare OME to order 𝜖 0 , which was obtained previously for all color structures in [22]. We checked explicitly that the obtained 𝑍 ns has the same form as in equation ( 8), and thus allows the determination of the 𝑁 𝑓 𝐶 3 𝐹 contributions to 𝛾 (3) ns from the single pole of 𝑍 ns . The flavor singlet but valence non-singlet contribution 𝛾 s ns in (10) vanishes for the color structure 𝑁 𝑓 𝐶 3 𝐹 , which allows us to write down our result in the following unified form,
with )
)
)
) ) 2 -192 (𝑛 + 1) 3 + 640 3(𝑛 + 1) 4 -320 3𝑛 + 112
)
)
) )
)
)
)
T. Gehrmann, A. von Manteuffel, V. Sotnikov et al.
-
)
)
)
)
)
)
)
)
) -1280 3(𝑛 + 1) 5 + 224 3𝑛 6 + 1184 3(𝑛 + 1) 6 -32 .
(
Here, we follow the notation used in [14] and omit the argument 𝑛 of the harmonic sums defined by
Our result contains harmonic sums up to weight 6. Unlike for the singlet anomalous dimensions we computed before, we notice that only two kinds of denominators 1∕𝑛 and 1∕(𝑛 + 1) appear in the above equation. Furthermore, the coefficients of the first power of 1∕𝑛 and 1∕(𝑛 + 1) differ by a minus sign only, thus we can always write them as a single term, for example,
The above all-𝑛 result in ( 16) is new. Evaluating the result for fixed 𝑛, we found full agreement with the fixed 𝑛 ≤ 16 results derived in [13].
Moreover, the all-𝑛 results for 𝜁 4 and 𝜁 5 terms have been derived in [70] and [13], respectively. For these results, we also found full agreement.
Applying an inverse Mellin transformation to equation (1), we obtained the corresponding splitting functions in momentum fraction 𝑥-space.
This is achieved with the help of the function InvMellin in the package HarmonicSums; alternatively, the method proposed in [71] could be used. The explicit expressions for the corresponding splitting functions are provided in the ancillary files.
It is interesting to study the various limits of the splitting functions. In the limit 𝑥 → 0, the result is free of power divergences and reads
) log(𝑥)
For 𝑃 (3),+ ns , the double logarithmically enhanced terms proportional to log(𝑥) 𝑘 with 𝑘 = 5, 4 have been derived in reference [72], and we find full agreement.
In the limit 𝑥 → 1, 𝑃
where the plus distribution is defined as
for a continuous test function 𝑓 (𝑥).
(
We note that in this notation, the perturbative expansions of the coefficients are again defined in powers of 𝑎 𝑠 , that is,
and similarly for the other coefficients. Interestingly, all results shown in equation ( 23) have been derived before, numerically for 𝐵
analytically for the others [74][75][76][77][78]. We find perfect agreement with the literature, thus providing another strong check of our all-𝑛 result in (16).
The coefficient of
denoted by 𝐴(𝑎 𝑠 ) determines the cusp anomalous dimension [79], and we find that our result above agrees with the 𝑁 𝑓 𝐶 3 𝐹 contribution to four-loop cusp anomalous dimension in [74][75][76]. The coefficient of 𝛿(1 -𝑥) denoted by 𝐵(𝑎 𝑠 ) is called the virtual anomalous dimension. In [73], the numerical result 𝐵 4 | | |𝑁 𝑓 𝐶 3 𝐹 = 80.780 ± 0.005 has been obtained, which agrees well with our analytic result shown above. The numeric result of [73], the analytic four-loop collinear anomalous dimensions of [80,76], together with the soft-rapidity correspondence derived in [81][82][83], allowed the numerical determination of the four-loop rapidity anomalous dimensions [84,85]
where for 𝛽(𝑎 𝑠 ) the limit 𝜖 → 0 of equation ( 7) is implied. Performing expansions for all-order results in the equation ( 26) to 𝑎 4 𝑠 , it reads (see also [13])
where the four-loop quantities 𝐶 4 and 𝐷 4 depend on 𝐴, 𝐵 and the 𝛽 function from lower-loop orders only. Our results verify the above conjecture for the color factor 𝑁 𝑓 𝐶 3 𝐹 explicitly.
We analytically computed the 𝑁 𝑓 𝐶 3 𝐹 contribution to the four-loop, non-singlet anomalous dimension for arbitrary Mellin moments 𝑛 for the first time. The method is based on the framework of the operator product expansion, through the computations of off-shell operator matrix elements. In contrast to the singlet case, the renormalization of the non-singlet contributions computed here is conceptually straight-forward. We introduced a tracing parameter to replace symbolic exponents depending on 𝑛. In this way, we were able to employ standard integration-by-parts reductions and