Three-loop vertex integrals at symmetric point

This paper provides details of the massless three-loop three-point integrals calculation at the symmetric point. Our work aimed to extend known two-loop results for such integrals to the three-loop level. Obtained results can find their application in regularization-invariant symmetric point momentum-subtraction (RI/SMOM) scheme QCD calculations of renormalization group functions and various composite operator matrix elements. To calculate integrals, we solve differential equations for auxiliary integrals by transforming the system to the $\varepsilon$-form. Calculated integrals are expressed through the basis of functions with uniform transcendental weight. We provide expansion up to the transcendental weight six for the basis functions in terms of harmonic polylogarithms with six-root of unity argument.


Introduction
Many physically significant quantities in QCD can be extracted from the three-point Green functions. There is remarkable progress [1][2][3][4] in calculating QCD renormalization group functions in the minimal subtraction(MS) renormalization scheme [5]. We can choose exceptional momenta routing(one of the external momenta set to zero) and restrict ourselves to considering only the divergent parts of diagrams for these calculations. Due to the unphysical nature of the minimal subtraction scheme, a group of momentum subtraction schemes(MOM) [6] is widely used in calculations requairing finite parts of the three-point functions. Such physical schemes are crucial for the Lattice calculations, where one has access to vertex functions directly. Calculation of various vertices in regularization invariant momentum subtraction (RI/MOM) schemes provides a connection between the lattice calculations and the MS scheme results. By choosing exceptional momenta routing, putting one of the external momenta to zero, we define RI/MOM scheme. With symmetric point kinematics configuration with all squares of external momenta equal, we define RI/SMOM scheme.
The main difficulty in the RI/MOM and RI/SMOM scheme calculation is the necessity to know finite parts of vertex functions in the chosen kinematics. If we require only finite parts of the propagator type integrals for the exceptional momenta routing, then for the vertices in the RI/SMOM case, we need to know three-point integrals at the symmetric point. Calculation of the symmetric point integrals is the main difficulty of the RI/SMOM scheme. The present paper focuses on how to solve the problem at the three-loop level.
There are many results based on the calculation of the two-loop three-point functions at the symmetric point. Namely three-loop RI/SMOM beta-functions [7], two-loop correction to the relation between RI/SMOM and MS quark mass [8,9], and renormalization of Wilson operator matrix elements [7][8][9][10][11]. We recently extended some of these calculations to the three-loop level analytically [12,13], using results of the present paper. Independently, numerical results for renormalization of different operators matrix elements appeared in the series of papers [14,15].
Integrals of our interest are given on the left of the Fig. 1 and condition p 2 1 = p 2 2 = q 2 = −1 defines the symmetric subtraction point. Before providing details of the three-loop calculation, we need to review techniques used at two loops. The first two-loop calculation at the symmetric point [16] used a large momentum expansion procedure. This technique does not rely on any symmetric point integrals knowledge and requires only results for massless propagator integrals. In subsequent two-loop calculations [7][8][9][10][11], authors adopted Integration-By-Parts(IBP) [17,18] reduction to the minimal set of two-loop master integrals known for a long time, even in general kinematics not restricted to the symmetric point [19][20][21]. This way, the problem of the master integrals calculation had not appeared explicitly in those works.
On the other hand, a very restricted subset of the three-loop master integrals needed for analytical calculations at the symmetric point [12,13] is known [22].
Our starting point for three-loop master integrals calculation is a paper [23], where linear reducibility of the two-loop three-point integrals with off-shell legs is demonstrated after an appropriate variable change. Linear reducibility of three-loop integrals with offshell kinematics was indicated in the E.Panzer Ph.D. thesis [24]. Since symmetric point integrals are a particular case of considered integrals, we expect them to be expressed in terms of generalized polylogarithms(GPL) as a direct consequence of the linear reducibility. The same property will have auxiliary integrals with p 2 1 = p 2 2 , but q 2 arbitrary, and we will be able to reduce differential equations(DE) system for such integrals to the so-called εform [25].
As the primary calculation method, we solve differential equations system for auxiliary integrals. Due to performance issues with the most complicated integrals, we use direct integration in terms of GPLs for partial checks only. Also, the method of direct integration requires construction of a new set of master integrals without sub-divergencies, which is possible to achieve by shifting space-time dimension and increasing some propagators powers [26]. Conversion between the set of finite integrals and the original basis is an additional complication of the method. However, an alternative based on the DE system solution is also highly nontrivial. We need to consider more complicated integrals depending Table 1. Two-loop auxiliary topologies on at least one dimensionless variable and provide boundary conditions(BC) sufficient to fix all needed integration constants in the DE solution. The paper organized as follows: we define auxiliary topologies containing all integrals of our interest in section 2, we describe the solution of differential equations for auxiliary integrals in section 3, we describe our approach to fixing boundary conditions with a detailed two-loop example in section 4, construction of the uniform transcendentality weight basis of integrals can be found in section 5, and in section 6 we provide the list of calculated master integrals. Actual results for the integrals are in the supplementary materials to the paper.

Notations for integrals topologies
To reduce many integrals in the course of calculations [12,13], we apply IBP reduction to the small set of master integrals considered in the present paper. To perform the reduction efficiently, we use the Laporta algorithm [27] implemented in the package Reduze2 [28]. We define auxiliary topologies containing the complete set of propagators to express all the appearing scalar products for the algorithm's work. To uniquely identify an integral inside topology, we use a vector of integer numbers corresponding to propagators' powers.
For a vertex integral at the symmetric point, we assign external momenta according to the left part of the Fig. 1 and set p 2 1 = p 2 2 = q 2 = −1. This definition makes all integrals real, and self-energy integrals are in one-to-one correspondence with integrals from the MINCER package [29,30]. Later, this property will be helpful for the boundary conditions fixing procedure described in section 4. At the one-loop level, we have a single topology A 1 with three propagators (2.1) At the two-loop level, we have two seven-propagator topologies, A 2 and B 2 where X = {A, B} and propagators' momenta are defined explicitly in table 1. At three-  Table 2 provides explicit expressions for momenta P i for each of the three-loop topologies with X = {A, B, C}. Table 3. Distribution of the number of integrals over topologies.
The number of unique master integrals after IBP reduction are given in table 3. We present all master integrals graphs in tables 4, 5 and 6.

Differential equations for auxiliary integrals
As stated in the introduction, our main calculation method is a solution of the DE system for specially constructed auxiliary integrals. We consider integrals similar to original three-point ones, but now with q 2 not fixed, and introduce scaleless variable x with q 2 = x p 2 1 . At first glance, this makes us calculate more complex integrals than the initial set. However, the advantage is that now we have access to the DE system, connecting point x = 1 with any other point, where boundary conditions can be constructed easily. In our calculation, we consider the limit x → 0 for fixing boundary conditions and the limit x → ∞ for checks.
We analyze the same topologies as before, but now with q 2 not fixed for auxiliary integrals. After performing IBP reduction and identification of master integrals for each of topologies, we have the following number of master integrals: 3 for A 1 , 12 for A 2 , 11 for B 2 , 88 for A 3 and 91 for B 3 . Since for topology C 3 at the symmetric point, we have only one integral, which factorizes into lower-loop integrals from topologies A 1 and B 2 , we do not consider it here.
To differentiate master integrals in variable x and reduce the result back to master integrals, we make use of the package Reduze 2 [28]. To convert the obtained DE system to the ε-form with first apply a variable change x = − (z−1) 2 z . In a new variable z, singular points of the DE system correspond to the following set of singular points in an old variable x: z The point z = e ± iπ 3 is a solution we are interested in, limit z = 1 we use in section 4 as boundary conditions to fix integration constants in the DE solution. The solution in the point z = 0 is used to derive known results for the there-loop massless form-factor type integrals [31][32][33], which is necessary for checking obtained results. We do not consider the point z = −1 in the present paper.
After the variable change, the system of differential equations for the vector of integrals f (z, ε) for each of the considered topologies has the following form With matrix M having only a limited set of ε-independent singularities of the finite order in the set of points This form of the DE system is ideally suited for further conversion to the canonical form [25]. Only one complication is the appearance of sixth-roots of unity residues at the three-loop order. For one and two-loop integrals, the DE systems are singular only in z = {0, 1, −1}.
Following the strategy from [25] and with algorithm [34] implemented in the publically available package epsilon [35], we convert the original DE system to another DE system for a new basis of integrals g(z, ε) = T −1 (z, ε) f (z, ε) (3.4) Obtained DE system matrix has only Fuchsian singularities, and ε dependence completely factorizes. Applying epsilon for reduction of the original system to the form (3.4), we perform all the steps of the algorithm [34] except the last one. Due to the performance issues, we were forced to find a constant matrix transforming the system to the form with ε factorized manually. Another reason to make the last step manually is that the found transformation matrix is not unique, as not unique the DE system matrix in ε-form. This freedom in the transformation matrix and related freedom in the ε-form matrix we use in section 5, where we construct a new basis of integrals with uniform transcendentality(UT) weight.
The differential equations system (3.4) has an excellent property that, after expanding all the master integrals in ε, the differential equations for the series coefficients g i (z) with g(z, ε) = g i (z)ε i decouple entirely, and the solution for the particular coefficient has the form: To perform integration, we utilize GPLs [36] integration rules, and starting with some order n of ε-expansion where all g n ≡ 0 we proceed by induction. At the next step we associate integration constants with zero weight GPLs, since G(; z) ≡ 1. Singularities of the system (3.4) define the alphabet of GPLs. In the examined case, we have a subset of the full six-root-of-unity alphabet considered in [37], and for lower loop orders, the alphabet of HPLs containing only {0, 1, −1} is enough. GPLs are directly integrated with G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t). (3.6) Each order of ε-expansion of the constructed solution is built from the GPLs of the uniform transcendental weight multiplied by some unknown constants. In the general case, actual values of integration constants r i can spoil the uniform transcendentality(UT) structure of the solution, and we address the question of restoring UT expansion in section 5. Our next step is to fix integration constants r i by providing appropriate boundary conditions.

Fixing boundary conditions
To fix integration constants in the solution of DE from the previous section, we consider its behavior near the point z = 1. Since the point z = 1 is a singular point of the DE system (3.4), we can not take the limit directly due to the logarithmic singularities log(1−z) appearing in the solution. Thanks to the Fuchsian form of the DE system (3.4), the leading order term of expansion in a small variablez = 1 − z can be constructed directly from the matrix residue A 1 [38] 1 : lim z→0 g(ε, z) =z εA 1 c(ε) + O(z). Momenta routing used for the asymptotic expansion coefficients calculation is given in Fig. 1(right). With this routing, our kinematics constraints are as follows P 2 = x Q 2 and P · Q = x 2 Q 2 , and results of the expansion in the limit x → 0 we can obtain from the large momentum procedure(LMP) considering Q large. For more details on the LMP, see [42] and references therein. Expansion of the arbitrary three-point vertex integral with external momenta P, Q in large momentum Q has the form: To fix coefficients a S i,j of the expansion, we use setup based on the package EXP. First, using EXP, we identify all needed subgraphs to expand in large external momenta, then we expand to sufficiently high order in the large momentum Q, further we perform tensor reduction to separate powers of the scalar products Q 2 and P · Q from loop integrals dependent on a single external momenta Q. The remaining two-point integrals are calculated with MINCER, which is fast enough to provide us with expansion terms to high order in variablē z. We determine coefficient a i,j as a linear combination of massless propagator-type master integrals, these integrals we can expand later up to the required order in ε.
After conversion of obtained expansions back with T −1 , we truncate the series to get the leading term for g and fix c(ε) from Eq. (4.1). Contributions of different subgraphs from different integrals to the same set of functions c provide us with a strong check on the validity of the procedure. Now expanding both solution of the system (3.4) in z → 1 and leading order solution (4.1) with c(ε) substituted, we can fix all integration constants up to required order in ε simply expanding MINCER master integrals deep enough in ε.
As a real example, we consider the DE system for non-trivial two-loop integrals J 1 = I A2 0,1,0,0,1,1,0 , J 2 = I A2 0,1,1,0,1,0,0 , J 3 = I A2 0,1,1,0,1,1,0 , J 4 = I A2 0,2,1,0,1,1,0 . (4.4) DE system for integrals j(z) now is in the ε-form with numeric matrices Due to the Fuchsian form of the DE system, the leading order solution in the limit z → 1 is determined by the matrix exponent: Withz = 1 − z it reads when we keep leading terms inz only Here all c i are functions of ε, and we see that some of the functions c i can be fixed independently from different integrals expansions stemming as an additional cross check.
In provided example c 1 can be fixed both from the j 1 and j 3 . With EXP we obtain the following expansions for original integrals: Where T 1 is the two-loop massless sunset integral with Q 2 = −1. Above we provide only orders needed to fix all c i (ε) entering (4.8). After multiplication with T −1 we find the required functions: After fixing the required integration constants, the obtained solution can be checked by considering the limit z → 0, corresponding to infinitely large (q 2 → ∞) external momentum squared. The hard subgraph's contribution in this limit allows us to extract results for the massless form-factor integrals and compare them with expressions known up to the threeloop order [32,33].
As before due to the Fuchsian form of the DE system, we can follow procedure described in [39] and construct general form of the expansion in z → 0 Expanding integrals g(z) in ε, we can find several first terms of a Laurent expansion of c j (ε) around ε = 0. With the help of the matrix T , we can construct expansion for the original basis of integrals f , which we utilize to extract terms with specific noninteger powers of z. We can extract the value I FF for the massless form-factor integral with L-loops and E internal edges from the hard subgraph contribution of the form: From the present paper results, we found agreement with the results of the paper [33] for all master integrals up to the transcendental weight six.

Basis of integrals with uniform transcendentality weight
The main problem with the already obtained solution is that its expansion coefficients lack the UT property. For practical applications like [12,13], we need to expand some of the integrals to high order in ε corresponding to the transcendental weight seven. Only part of the integrals with GPLs up to the transcendental weight six enters the final result, and all GPLs with transcendental weight seven cancel in the sum. At intermediate steps, manipulations with a weight seven GPLs are highly complicated, since reduction rules are available up to the weight six only [37].
Possible solution to the problem is a new basis of UT weight functions, which need to be known up to the transcendentality weight six to express results containing GPLs up to the weight six. Explicitly we want to find a basis of pure functions U = T −1 UT f , with ε-expansion in the form U = 6 j=0 u j ε j + O(ε 7 ), where each coefficient u j has uniform transcendental weight j.
The UT integral basis construction is a complicated task even with special tools designed to attack the problem [43,44]. Fortunately, with already available results we can easily avoid all these difficulties.
First, we notice that if the integrals with arbitrary z have UT weight, this property also holds for the integrals at the symmetric point(z → e iπ 3 ). The weight of GPLs and all z-dependent prefactors do not change transcendentality after taking the limit, and we can focus on constructing UT basis for auxiliary integrals.
Second, due to the UT property, the DE system for UT integrals is also in ε-form. Our goal is to find a transformation matrix between canonical master integrals from section 3 and the new UT basis. As mentioned in section 3, the last step of reduction to the ε-form has afreedom, and we are looking for the transformation to UT basis by varying these parameters.
Another important observation is that for the z-dependent integrals with UT weight, this property also holds when considering their expansion near singular points. For example leading terms of expansion in z → 1 used for fixing boundary conditions in section 4.
Since ε-expansion of the prefactorz a i ε in (5.1) has transcendental weight zero, all coefficients C i for expansion of the UT integral J UT also have UT weight. In considered limit, all coefficients C i are built from massless propagator integrals only. After reduction to master integrals, most of them are known for arbitrary D, and for the remaining, transformation to the UT basis is known from [45]. In this way, we find transformation matrix T UT by analyzing expansion coefficients in the limit z → 1.
All UT integrals constructed in this way are regular in the limit z → e iπ 3 and going back to the original integrals and then reducing to the basis of symmetric point integrals, we obtain representation for integrals of our interest in terms of pure functions.

Results and conclusion
In addition to the described method, we use IBP reduction to the new basis of finite integrals to follow the strategy described in the paper [26]. With package HyperInt [46], we find an analytical solution for several integrals, except the most complicated, which is in complete agreement with results of the previuos sections. All integrals calculated in the paper were checked numerically with the sector decomposition approach implemented in the package pySecDec [47].
In tables 4, 5 and 6, we present all master integrals needed for calculations at the symmetric point up to the three-loop order. Analytical results for integrals can be found in supplementary files to the paper and contain a transformation matrix to the basis of UT weight functions and expansion of latter up to the transcendental weight six.
Up to the three-loop level, expansion of the basis functions is expressed in terms of HPLs with the sixth-root of unity argument. Instead of the basis constructed in paper [48], which has an attractive feature, that reduction rules for real(imaginary) part of HPLs contain only real(imaginary) parts of basis functions. We construct a new basis allowing us to present results in a more compact form. Also, our basis contains as basis elements products of lower weight functions, which simplifies products of lower-loop integrals during renormalization. The total number of elements N w at weight w is the same N w = 2F 2w , where F n is the n-th Fibonacci number and explicitly is given by N w = {2, 6, 16, 42, 110, 288} for the weights considered in the paper.
In the paper, we calculated all massless three-point master integrals at the symmetric point. As the solution method, we apply differential equations for specially constructed auxiliary integrals and fix boundary conditions from the large momentum expansion. Results for integrals are expressed through the basis of functions with uniform transcendental weight. We provide expansion in ε for these functions in terms of the HPLs with the sixth root of unity argument up to the transcendental weight six. The obtained result can find its application in future calculations in the RI/SMOM scheme and as the boundary conditions for three-point integrals in more general kinematics.