Dispersion Approach in Two-Loop Calculations

The higher-order corrections become increasingly important with experiments reaching sub-percent level of uncertainty as they look for physics beyond the Standard Model. Our goal is to address the full set of two-loop electroweak corrections to M{\o}ller or electron-proton scattering. It is a demanding task which requires an application of various approaches where two-loop calculations can be automatized. We choose to employ dispersive sub-loop insertion approach and develop two-loop integrals using two-point functions basis. In that basis, we introduce a partial tensor reduction for many-point Passarino-Veltman functions, which later could be used in computer algebra packages. In this paper, we have considered self-energy, triangle and box sub-loop insertions into self-energy, vertex and box topology.


I. INTRODUCTION
The electroweak precision searches for the physics beyond the standard model (BSM) frequently demand a subpercent level of accuracy from both experiment and theory. From the theory perspective, this can be achieved by extending the perturbation expansion of the scattering matrix element to the two-loop level. However, since the electroweak (EW) interaction usually introduce differentmass propagators and higher-order tensor Feynman integrals, the two-loop EW calculations can become increasingly complicated. In general, in electroweak case, it is very challenging or even not possible to find analytic results beyond the one-loop level, so one would either have to resort to various approximations or purely numerical methods. See, for example, a comprehensive overview of a variety of numerical loop integration techniques in [1], a general case of the two-loop two-point function for arbitrary masses in [2], and a method of calculating scalar propagator and vertex functions based on a double integral representation in [3,4]. The more recent developments on analytical evaluation of two-loop self-energies can be found in [5][6][7][8][9][10], and on numerical evaluation of general n-point two-loop integrals using sector decomposition in [11,12]. In addition, we normally have to evaluate several thousands Feynman graphs. Obviously, such a voluminous task should be delegated to the computer-based calculations. Although there has been a strong progress, we still have an ongoing problem of having to deal with cumbersome expressions and consequently be forced to introduce approximations. In many existing techniques, the tensor integrals have to be reduced to master scalar integrals which increases size of the final expressions dramatically. In order to address this, we employ a dispersive approach to subloop diagrams, and introduce a partial tensor reduction of the two-loop graphs. In general, a subloop can be represented through a dispersion tensor integral operator with a relatively simple propagatorlike structure. Dispersion tensor integral numerator could be then absorbed into the effective Feynman vertices or propagators, and the secondloop integration will acquire an additional propagator. The idea of the subloop insertions with the help of the dispersive approach was implemented for the self-energies [13,14] and partially for the vertex graphs with the help of Feynman parametrization [15]. We extend this for self-energy, vertex and box subloop insertions of the general tensor structure. In addition, we apply the reduction of the three-, four-, and five-point tensor coefficient functions insertions to the derivative representation of the two-point tensor function. Of course, as in previous works, the treatment of the UV and IR divergences requires subtractions derived from the two-loop EW counterterms and introduction of the photon mass regulator. In Sec. II, we start with basic definitions and ideas of dispersive treatment of subloop insertions. After that, we consider self-energies, vertex and box insertions and provide a partial tensor reduction. Section III considers specific example of two-loop self energy with vertex-type insertion and provide numerical comparison with [13] for the kinematical region below and above threshold.

II. SUBLOOP APPROACH
Generally, a dispersion relation allows to express a loop integral through the known imaginary part: Here, q 2 is the external momentum squared and s 0 is the branch-point position on the real axis. The imaginary part ℑLðq 2 Þ can be calculated from discontinuities of the loop integral using Cutkosky rules. If we consider the subloop insertion represented by the self-energy, triangle, or box, we can extract an imaginary part of two-, three-, and four-point tensor coefficient functions from the routines such as FF [16] and LOOPTOOLS [17], which are already implemented numerically. This leaves us with a problem of expressing the two-loop tensor integrals in terms of many-point tensor coefficient function. We start with definition of a general two-loop tensor integral in the dimensional regularization: where q 1;2 are the integration momenta in the first and second loops, respectively. The momenta k i;j;l represent various combinations of the external momenta p i;j;l from a two-loop graph. The masses of internal particles are defined as m i;j;l . For the processes specifically related to the parityviolating scattering, a subloop topology would be defined by either self-energy, triangle or box insertions.

A. Self-energy subloop
The self-energy subloop could be inserted into another self-energy, triangle or box topology (see Fig. 1). After replacing self-energy subloop by the dispersion integral, graphs from Fig. 1 could be reduced to graphs shown in Fig. 2. More specifically, for fermions or vector bosons, the self-energy subloop can be defined in form of the Lorentz covariant terms: Here, in Eq. (3), Σ V−V T;L ðq 2 Þ represents transverse and longitudinal parts of diagonal and mixed vector boson self-energies. In Eq. (4), Σ f L;R;S ðq 2 Þ represents left, right and scalar parts of the fermion self-energy graph. The ω AE ¼ where ℑB i;ij;ijk ðs; m 2 α ; m 2 β Þ can be easily computed using LOOPTOOLS or FF libraries. As a result, Eqs. (3) and (4) can be rewritten in the following form: The summation in Eq. (6) is done over all possible internal particles in the self-energy subloop carrying masses m α and m β . Using Eq. (6), we can now write general expressions for the two-loop topologies in Fig. 1. In the case of a subloop represented by a vector boson self-energy, we get where F μ 1 μ 2 ðq 2 ; s; m α ; m β Þ is defined as For the subloop insertion stemming from fermion selfenergy, we can write and In both Eq. (7) and Eq. (9), the integration over the second loop momenta can be written in the form of the Passarino-Veltman many-point tensor coefficient functions. For selfenergy subloop insertion specifically, we can simplify the second-loop integral by means of partial fraction decomposition. Since the momentum running in the adjacent propagators and dispersion term ðs − q 2 − iϵÞ −1 is the same, we can write This gives the following result for a case of the vector boson self-energy insertion Functions L a;μ 1 μ 2 ;ν 1 …ν R and L b;μ 1 μ 2 ;ν 1 …ν R depend on two-, three-, and four-point tensor coefficient functions for the topologies in Fig. 1, defined from left to right, respectively. Each of the many-point functions are dependent on the integration parameter s, masses and combinations of external momenta. For the two-loop result with fermion self-energy insertions, we can write: As in the previous case, the functions N 1;M;1 a;b;c;ν 1 …ν R are defined through functions B, C and D in similar fashion. Integration in Eqs. (12) and (13) can be completed numerically after subtraction of the UV divergences. In case where it is possible to perform subloop subtraction (i.e., there is no global UV divergence), we can use the self-energy subloop which has already subtracted terms. For example, γ − γ self-energy is: Clearly, Σ γ−γ ð0Þ ¼ 0, but we keep this term in Eq. (14) anyway for the dispersive subtraction. This way, the second-loop integral in Eq. (6) will acquire term q 4 s 2 , which will cause the cancellation of two γ propagators around γ − γ insertion, and hence we can omit use of partial fraction expansion in Eq. (11) and still write the final second-loop integration in the form of B, C and D functions for the topologies defined from left to right in Fig. 1. The same can be done for Z − Z, γ − Z and W − W mixings insertions. For Z − Z and W − W insertions, we can write, in the on-shell renormalization scheme, the following [18]: Equations (15) and (16) would cancel Z − Z or W − W and γ − Z propagators in the second loop integration, and, as a result, functions L a;b in Eq. (12) will again depend on B, C and D Passarino-Veltman tensor coefficient functions. In the case of f − f subloop, we can apply the following on-shell subtractions: Here, the functions I L;R;S have the following integral representation: Substitution of Eq. (17) into the second-loop integration will result in the cancelation of ðq 2 − m 2 f Þ in I L;R;S with one of the fermion propagators. Terms a L;R;S do not result in cancellations, but they also do not introduce the dispersion denominator ðs − q 2 − iϵÞ −1 into the second-loop integration. After applying partial fraction expansion of the denominators with the same momentum, functions N a;b;c would depend on B, C and D tensor coefficient functions only.
The structure of the insertions in Eq. (6) suggests that we can introduce V − V and f − f effective mixing propagators with a dispersion integral removed. All the functions of parameter s could be left un-evaluated during the secondloop integration. This gives us the possibility to employ a computer-algebra approach, where the second-loop integral could be evaluated analytically, and after subtractions the dispersion integration can be carried out numerically. For V − V effective mixing we can write that as a combination of the transverse and longitudinal propagators: When evaluating the second-loop integral we can leave imaginary parts of Σ T;L un-evaluated and get analytical structure for the two-loop graph. In case subtraction is possible at the subloop level, V − V effective propagators would have the following structure: Here, the functions T V−V ðs; m 2 V Þ and ℑΣ T;L are independent of the second-loop momenta and could be left unevaluated until dispersion integration is performed. For the specific V − V mixings, functions in Eq. (20) are defined in Table I. For fermion mixing, we can also introduce following effective propagator For the subtracted f − f subloop insertions, we can replace Eq. (21) and introduce, with the help of Eqs. (17) and (18), the following set of effective fermion propagators derived from the first-loop integration:Π where The UV-finite functions d L;R;S ðs; m 2 α ; m 2 β Þ can be derived from Eq. (18) using a dispersive representation of the constants a L;R;S : The implementation of the effective propagators in Eqs. (19), (20), (21), and (22) is straightforward in the computer-algebra packages such as FormCalc or Form ( [17,19]), and as a result it is possible to construct the two-loop self-energies matrix elements in the analytical form. If additional subtractions are needed, it can be done later using the second-order EW counterterms ( [15,20]). Finally, the last step of the calculations would be numerical integration. This can be done with the help of the numerical libraries from LOOPTOOLS or FF and integration routines such as VEGAS or QUADPACK. Since the dispersion integration would involve many-point tensor coefficient functions (above the two-point functions), a numerical stability could become a concern. Out of all the Passarino-Veltman functions, only two-point tensor coefficient functions have well-defined analytical structure and therefore the most stable numerically. It would be most desirable if we could write three-, four-, and five-point functions (for the cases of triangle and box type of the subloop insertion) which enter functions L a;b and N a;b;c in Eqs. (12) and (13) as some representation of the two-point tensor coefficient functions. This can be achieved if we combine the Feynman trick [20,21] with derivative representation of the many-point functions. Let us now consider three-, four-, and five-point functions separately.
We start with the scalar three-point function and later consider results for the C i;ij;ijk Passarino-Veltman functions. The general expression for the three-point scalar function is given by With the help of Feynman trick, we can join the first two propagators in Eq. (23), and after shifting momentum wherex ¼ 1 − x. Term ½ðτ − ðp 1x þ p 2 ÞÞ 2 − m 2 12 −2 can be replaced after shifting mass m 2 12 by a small parameter λ.
As a result, the expression for the three-point function can be reduced to the derivative representation of two-point function [see Fig. 3]: On the contrary, in [21], the differentiation in Eq. (26) would be done with respect to m 2 12 , but this would require the analytical differentiation first (since m 2 12 is a function of the Feynman parameter), and then integration over the Feynman parameter. In our case, the integration over the Feynman parameter, and then differentiation with respect to λ can be all done numerically. Using the approach from Eq. (26), we can replace the many-point tensor functions in the second-loop integration by the derivatives of the twopoint function. Taking into account that the LOOPTOOLS libraries have numerical implementation for the regularization of UV and IR divergences, we can perform the twoloop calculations using the two-point functions basis and, later, after the appropriate subtractions, perform the dispersive and the Feynman integration numerically. The final step would be a numerical evaluation of the derivative of the integrated result. Since LOOPTOOLS and FF libraries are capable of computing tensor coefficient functions of the higher-rank tensors, we can derive the partial tensor reduction for the three-point function in terms of the two-point basis. For the function C μ 1 …μ N , we can write For the vector C μ , we can use the tensor decomposition and in parallel apply the Feynman trick on the right-hand side of Eq. (27). After using the derivative approach from Eq. (26), and shifting the momenta as before, we get where B 0;μ ≡ B 0;μ ððp 1x þ p 2 Þ 2 ; m 2 3 ; m 2 12 þ λÞ. Using the tensor decomposition for B μ ¼ −ðp 1μx þ p 2μ ÞB 1 , and matching terms in front of p 1μ and p 2μ , we can solve for C 1;2 in terms of B 0;1 functions: The same idea can be extrapolated to the higher orders of the three-point functions. The results for the C μν and C μνα partial reduction are given in Appendix A. An important advantage of the partial tensor reduction is that we can substantially reduce the size of the final expressions in the two-loop integrals. The reduction of the three-point functions to the two-point basis can be also employed in the dispersive representation of C 0;μ;μν;μνα which could prove helpful if we use the triangle subloop insertion. For the four-point function, we will use an analogous approach. We start with the general structure of scalar D 0 function: In this case, we join the first three propagators, and after shifting momentum Obviously, the reduction to the two-point B 0 function is achieved by the second-order differentiation with respect to mass shift parameter λ: Again, the partial reduction for D μ;μν;μνρ;μνρσ can be done in the similar way as in Eq. (28) (Appendix B). The five-point functions are reduced as before with the help of the Feynman trick, the third-order differentiation with respect to the mass shift parameter λ, and shift of momentum q ¼ τ − Σ 4 i¼1 p i . Here, we can write where m 2 x − p 2 12ȳ y − p 2 123z z þ 2xyðp 1 p 12 Þ þ 2xzðp 1 p 123 Þ þ 2yzðp 12 p 123 Þ and p ij ¼ p i þ p j , p ijk ¼ p i þ p j þ p k . Results of the reduction for E functions are given in Appendix C for the tensor coefficients functions up to the fourth rank. As one can see, we can construct the twoloop integrals with the self-energy subloop insertions using only the two-point basis. Employing the effective propagators, we can now construct model files for the FeynArts [22] package and obtain the two-loop results in FormCalc or FORM. The final UV-and IR-regularized expressions would have to be renormalized by means of subtractions. If subloop subtraction is not possible, we can regularize the insertion by a cutoff of the dispersion integral. In this case, the renormalization constants from counterterms responsible for the cancellations of subloop UV divergences, would have to be calculated using the dispersive representation with the same cutoff. In the last stages, numerical dispersion integration and differentiation would have to be done. The same ideas can be extrapolated in the case of the triangle-type of insertions in the self-energy, vertex or box diagrams.

B. Triangle subloop
Examples of the triangle subloop insertion in two-loops topology are shown in Fig. 4. Our starting point here would be to construct the dispersive representation of the threepoint function, which later could be used in the second-loop integration. To simplify, we will consider the case in which one of the external legs of the triangle insertion is put onshell (see Fig. 5). This could be a case shown in Fig. 4, for the triangle insertion in the box acting as the second loop. Considering that all particles in the loop are scalars, the graph in Fig. 5 is a three-point scalar function, and using the notation in Eq. (23), we can write where we used p 1 ¼ −k, p 2 ¼ q 2 and p 3 ¼ k − q 2 . In order to replace the three-point function in Eq. (33) by the dispersive representation [23], we can again use the Feynman trick and join the first two propagators. It is important to apply the The two-point function can be easily written as a dispersion integral and substituted into Eq. (34): The branch point of the two-point function is on the real axes (m 2 12 > 0), and hence the dispersion integral in Eq. (35) is well defined. Since B 0 function is UV-divergent, the integral in Eq. (35) would diverge, which is addressed by introducing a cutoff Λ. After differentiating numerically, Λ dependence will cancel. We checked numerically that the final result in Eq. (35) in fact does not depend on the cutoff parameter Λ. Using LOOPTOOLS libraries, we can compare left-and right-hand parts of Eq. (35) and see good agreement (see Fig. 6). In general, there will be cases when m 2 12 ¼ m 2 1x þ m 2 2 x − p 2 1 xx could become negative for the specific values of Feynman parameters or external momenta. In this case, we would choose an integration contour over the upper half-part of the complex plane, and as a result perform dispersion integration from −Λ 2 to Λ 2 cut-off. To demonstrate, we consider a scalar two-point function with arbitrary imaginary mass. Since B 0 function is UV-divergent, we will take derivative of B 0 with respect to mass shift parameter λ, and for lim λ→0 ∂ ∂λ B 0 ðp 2 ; m 2 1 ; −jm 2 2 j þ λÞ we can write: If we take, for example, Λ 2 ¼ 10 10 GeV 2 , m 2 1 ¼ 1 GeV 2 and jm 2 2 j ¼ 9 GeV 2 , we observe that results obtained by means of Eq. (36) and LOOPTOOLS shows up discrepancy only after 6th digit (see Table II). Using the idea outlined in Eq. (35), we can also derive dispersive representation for the higher order three-point function. With the help of Eq. (29) and Appendix A, we show in Fig. 7 results for the C 1;2 and C 00;11;12;22 functions. Since C 00 is UV divergent, we obtained finite result after subtraction at q 2 2 ¼ 0. The results are in good agreement with LOOPTOOLS (see Fig. 8). We also have tested third-rank three-point functions and found an excellent agreement. In general, we can write that triangle subloop can be replaced by the following effective coupling OperatorD is defined asD ¼ lim λ→0 ∂ ∂λ R 1 0 dx R Λ 2 rðx;λÞ ds…, and function rðx; λÞ has the following structure: rðx; λÞ ¼ ðm 3 þ ðm 2 12 þ λÞ 1=2 Þ 2 θðm 2 12 Þ − Λ 2 θð−m 2 12 Þ. Momentum p 1 can be a combination of external momenta only. Momentum at the subloop level, then Eq. (37) has to be modified as follows: The general algorithm of calculations with the triangle subloop insertions could be summarized as follows. First, calculate one-loop triangle insertion in Passarino-Veltman basis. Next, replace all three-point functions by dispersive representation using rotation of the arguments, so that momentum p 1 does not depend on the second-loop integration momentum. After that, add the term ðs − ðp 2 þ p 1x Þ 2 − iϵÞ −1 from Eq. (37), or (37), or from Eq. (38) to the second-loop integration. The next stage is to get the second-loop integral, again in Passarino-Veltman basis, and rewrite three-, four-, and five-point functions in the two-point basis. After that, apply the subtractions for the second loop using counterterms in a given renormalization scheme. The final stage is to evaluate the dispersion and Feynman parameters integrals, and do numerical differentiation at least once.

C. Box subloop
Box subloops can be incorporated in double boxes, ladder boxes and crossed-vertex graphs (see Fig. 9). As in the case of triangle subloop, our starting point would be an example of the double-box diagram. Here, we would define dispersive representation of the four-point insertion. In our case, all the external particles are on-shell, and masses on the top and bottom lines of Fig. 10 are equal to m 1 and m 2 , respectively. Following the momenta notation in Fig. 10, we can write for the left box subloop: After joining the first three propagators, using the mass-shift parameter approach, shifting momenta τ ¼ q 1 þ q 2 þ k 3 − k 1 and applying dispersive representation of the two-point function, we get the following result for the box subloop insertion: Mass m 2 123 ¼ m 2 3 ðx − yÞ þ x 2 m 2 2 þ y 2 m 2 1 − 2xyðk 1 k 2 Þ could become negative for specific values of x and y, so the function rðx; y; λÞ is defined as rðx; y; λÞ ¼ ðm 4 þ ðm 2 123 þ λÞ 1=2 Þ 2 θðm 2 123 Þ − Λ 2 θð−m 2 123 Þ. Here, θðxÞ is a usual step function. The second-loop integral will receive, from Eq. (40), an additional propagator and would become a four-point function with dispersion parameter s playing a role as the mass: OperatorÎ λxys has the following definition:Î λxys ¼ lim λ→0 dy R Λ 2 rðx;y;λÞ ds…. For the second-loop integration, we would also apply approach outlined in Eq. (31), and after joining the first three propagators, using mass-shift parameter ξ, final two-loop result can be expressed completely in the two-point basis: OperatorÎ ξzω is defined asÎ ξzω ¼ lim ξ→0 For the generalized box subloop, we can replace it by four-particle coupling where operatorD defined asD ¼ lim λ→0

III. NUMERICAL EXAMPLE
Lets consider a case with the triangle (three-point function) insertion applied as an example. Namely, we will consider the diagram in Fig. 11 with all the masses are different and external particles are off-shell. Following the notation of [13], a general expression for that graph is: After joining the first two propagators, shifting momenta as τ ¼ q 1 þ q 2 , and applying differentiation with mass shift parameter λ, we get where m 2 12 ¼ m 2 1x þ m 2 2 x − p 2x x and rðx; λÞ ¼ ðm 3 þ ðm 2 12 þ λÞ 1=2 Þ 2 θðm 2 12 Þ − Λ 2 θð−m 2 12 Þ. For the timelike process when p 2 > 0, a value of m 2 12 could become negative for the conditions above the threshold. That is addressed by using Eq. (36). In spacelike cases of p 2 < 0, the dispersion representation of the insertion is well defined. The second-loop integral, with the help of Eq. (26), can be written in two-point function basis. The final two-loop result can now be given in a compact form: where m 2 45 ¼ m 2 4ȳ þ m 2 5 y − p 2ȳ y. The integration and differentiation in Eq. (46) can be done numerically. Numerical integration, in the case of below threshold conditions, have used the global adaptive algorithm implemented in the integration package of Mathematica program. For the above threshold conditions, we have implemented numerical CUBA libraries from [24,25]. These libraries are specifically designed for the multidimensional integration and employ Monte Carlo or cubature rules of polynomial degree algorithms. They showed the best convergence for the integration above threshold conditions. After that, we can compare our results with [13] (see Table III). We get a good agreement of the results and extend our calculations to the spacelike regime as well. In Table III, we also provide computing time and compare it to [13] for the case of twodimensional integration using VEGAS routine. As we can see, computing time from [13] is nearly constant, but in our case it is highly dependent on the value of jp 2 j. Further analysis shows that in our case we have to apply numerical differentiation with respect to mass-shift parameters, which is sensitive to the variations of λ and ξ in Eq. (46). Sensitivity grows when jp 2 j is getting large, and in order to achieve stability in numerical differentiation, precision of the integration over Feynman parameters should increase substantially. This is primarily the cause of the increase in integration time. The computing time of the dispersion integral is usually in the order of fraction of a second, and numerically very stable. In the case of the box-type insertion, we would need, on top of the multidimensional integration, to deal with the second-order differentiation. That is a much more challenging task and will require twoloop graphs, with box insertions, to be evaluated using C++, FORTRAN, or PYTHON languages. This will be addressed in follow-up publications dedicated to the numerical evaluations of the dispersive representations of two-loop graphs with the box-type insertions.

IV. CONCLUSION
In this work, we have applied the dispersive treatment approach of the subloop insertion and represented the twoloop results in the two-point function basis. The secondloop integration was reduced to the two-point basis with the use of the partial tensor reduction. The partial tensor reduction simplifies analytical expressions considerably to the point that it is possible to employ computer algebra evaluating the two-loop calculations analytically and carry out integration and differentiation numerically after that. As an example, we have compared our results for the double self-energy shown in Fig. 11 with [13] and found an excellent agreement.
APPENDIX B: D μ , D μν , D μνρ , AND D μνρσ dy…, we list only the final results D μ :