Two-loop Vertices with Vacuum Polarization Insertion

We present the analytic evaluation of the second-order corrections to the massive form factors, due to two-loop vertex diagrams with a vacuum polarization insertion, with exact dependence on the external and internal fermion masses, and on the squared momentum transfer. We consider vector, axial-vector, scalar and pseudoscalar interactions between the external fermion and the external field. After renormalization, the finite expressions of the form factors are expressed in terms of polylogarithms up to weight three.


Introduction
The analytic evaluation of the two-loop corrections to the electron form factors in Quantum Electrodynamics (QED) can be considered among the pioneering projects that triggered the developments of mathematical techniques and concepts for the evaluation of multi-loop Feynman integrals in Quantum Field Theory, which is still ongoing nowadays.Originally addressed by means of dispersion relations and giving a finite mass to the photon, in order to regulate the otherwise divergent integrals [1][2][3], the contributing vertex graphs were also later evaluated [4,5] within the dimensional regularization scheme, by using the differential equations method [6][7][8].The results of [3][4][5] were among the first studies showing that classical polylogarithms could not represent an exhaustive set of functions for Feynman integrals beyond one-loop, therefore, pointing to the need of introducing an extended set of polylogarithms, such as the Harmonic Polylogarithms (HPLs) [9] -later embedded in the wider class of Generalised Polylogarithms (GPLs) [10].HPLs turned out to be useful both for the direct integration of the dispersive integrals [3], and for solving the differential equations of the master integrals (MIs) in terms of iterated integrals [4], in combination with suitable change of variables, needed to rationalize the integration kernels−a procedure that would have been later dubbed as alphabet rationalization.
The feasibility of the analytic evaluation of the electron form factors at two-loop and beyond in QED was the basis of successive studies involving the heavy-quark form factors at two-loop in Quantum Chromodynamics (QCD), in the case of a vector interaction, and later extended, to account for the scalar, pseudoscalar, and pseudo-vector interactions [11][12][13][14][15][16][17][18][19][20][21][22][23].From a more formal point of view, form factors turn out to be also important for investigating the singular behavior of massive amplitudes in gauge theories, through their relation to the corresponding massless approximation [24].
In the context of scattering processes, the same two-loop vertex diagrams considered in the above studies appear, on the one side, among the (factorised) diagrams contributing to processes with the same number of external particles, yet with a larger number of loops [25][26][27][28][29][30], and, on the other side, to processes having the same number of loops, yet with more legs.In this respect, the two-loop threepoint integrals contributing to the vector form factors also appear in the evaluation of the two-loop QED corrections to the amplitude of the four-fermion scattering [31][32][33] and in the two-loop QCD corrections to heavy-quark pair production in the light-quark fusion channel [34].
In this work, we present, for the first time, the contribution to the vertex form factors of a heavy lepton with mass m e , coupled to a generic external particle, i.e. through a vector, axial-vector, scalar and axial couplings, coming from a two-loop graph with the insertion of a (vector) gauge boson vacuum polarization, with a closed-loop of a different type of heavy lepton, with mass m i ̸ = m e .
We hereby address the evaluation of the renormalised form factors by decomposing them, via integrationby-parts identities (IBPs) [35,36] and Laporta's algorithm [37], to a linear combination of seven master integrals, and by evaluating the latter using the Magnus method for differential equations [38].
After UV-renormalization, the renormalized form factors are finite and carry the complete dependence on the squared transferred-momentum q 2 = s, as well as on both the internal and the external lepton masses, respectively m i and m e .They are expressed in terms of GPLs up to weight three; equivalent expressions in terms of classical polylogarithms are given as well.
In the case of vector coupling, the evaluation of the diagram considered in the current work was previously discussed in [39,40], where semi-analytic expressions of the form factors were given as onefold integrals of kernels that are computed analytically by means of the hyper-spherical integration method.The numerical evaluation of our analytic results, by means of [41], is in perfect agreement with the results of [39] implemented in [42].An independent evaluation of the same diagram was also considered in [31], where the MIs were also computed using the method of differential equations.Our calculation is performed using a different change of variables, yielding a simple structure of the system of differential equations.The set of MIs presented in this work have been successfully checked against the numerical values provided by SecDec [43] and AMFLow [44] and are in full agreement with those presented in [31].We also verify that the vector form factor F 2 at zero momentum transfer agrees with the known expression of the anomalous magnetic moment given in [45,46].
The analytic expression of the vector form factors presented in this work can be directly applied in updating the analyses of the next-to-next-to-leading order (NNLO) QED corrections to the four fermion scattering with two massive lepton species.Additionally, they can be used also to complete the analytic evaluation of the QCD corrections to the heavy-quark form factors [13][14][15][16][17], as well as the QCD corrections to the Higgs boson decay in a pair of b b-quarks [22].
The paper is organized as follows.In Section 2, the definition of the vertex function is introduced, and the corresponding form factors for the vector, axial vector, scalar and pseudoscalar are discussed in Section 3. In Section 4, we discuss the integral decomposition and the evaluation of the master integrals.In Section 5, we discuss the renormalization procedure.Finally, the results and conclusions are presented in Sections 6 and 7 respectively.Appendix A contains further details on the structure of the matrices appearing in the system of differential equations obeyed by the master integrals.In Appendix B, we elaborate on the renormalization and on the evaluation of the one-and two-loop counterterm diagrams and of the relative renormalization constants.

Vertex diagrams with vacuum polarisation insertion
We consider the two-loop vertex diagrams V (k) (the index (k) refers to a labelling introduced earlier in the literature [4,5,[13][14][15][16][17]) pertaining to a generic (axial-)vector or (pseudo-)scalar boson of momentum q µ with virtuality q 2 = s that couples to an external on-shell fermion-antifermion pair, respectively carrying momenta p 1 and p 2 with p 2 1 = p 2 2 = m 2 e .The diagrams are subject to correction through the insertion of a vacuum polarization involving a massless vector boson coupled to a fermion-antifermion pair of mass m i ̸ = m e , as depicted in Fig. 1.The corresponding expressions are given by with where: α is the fine structure constant, µ 2 is the mass-scale introduced to keep the fine structure constant dimensionless in dimensional regularization, with d = 4 − 2ϵ, and k i (i = 1, 2) are the loop momenta.The fermion and the gauge boson propagators (in Feynman gauge) are respectively defined as S(q, m) = / q + m q 2 − m 2 + iε and P µν (q) = − g µν q 2 + iε (2.3) with the Minkowski metric g µν .The symbol J depends on the type of external boson, which is given by -)vector coupling, −i (g S + g P γ 5 ) , (pseudo-)scalar coupling. (2.4) The quantities g V , g A , g S , g P represent the vector, axial vector, scalar and pseudo scalar coupling constants respectively.

Form Factors
Following [5,13,14], a vertex Γ µ describing the (axial-)vector coupling of fermions to a gauge boson admits a decomposition in terms of four scalar form factors F 1,2 and G 1,2 (valid at any number of loop and any topology).Depending on whether we consider (axial-)vector or (pseudo-)scalar coupling, the number of form factors changes.For the former case, although the general decomposition involves six form factors, two do not survive due to gauge invariance.Hence, the decomposition in terms of four form factors F 1,2 and G 1,2 reads [5, 13, 14], with The form factors F i and G i (with i = 1, 2) can be extracted from Γ µ by means of suitable projectors, P µ F i and P µ Gi as The explicit expressions for the projectors in d = 4 − 2ϵ dimensions are given by with and In the derivation of the aforementioned projectors P µ Gi (s, p 1 , p 2 ), we assume an anticommuting γ 5 in d dimensions.However, we employ a non-anticommuting γ 5 , as introduced by 't Hooft-Veltman [47] and Breitenlohner-Maison [48] for the computation of form factors in eq.(3.2).The latter is defined as We treat the Levi-Civita symbol following the prescription in [49,50].In particular, the contraction of the ϵ µνρσ with the one from the projector is done according to the usual mathematical identity in four dimensions, but with the Lorentz indices of the resulting spacetime metric tensors all taken as d-dimensional.This is often known as Larin's prescription.As the Feynman diagrams considered in this article do not exhibit any anomalous behaviour−or, in other words, they fulfil non-anomalous Ward identities−the finite remainder is guaranteed to be independent of the prescription (commuting or anticommuting) adopted for γ 5 [51,52] 1 .
The vertex function for the (pseudo-)scalar admits a decomposition in terms of two scalar form factors F S and F P [17] as Γ = −i g S F S (s, m i , m e ) + g P γ 5 F P (s, m i , m e ) (3.9) with g S , g P the scalar and pseudoscalar couplings, respectively.The form factors F S and F P can be extracted from Γ by means of suitable projectors, P S and P P , as where the expressions for the projectors are given by The form factors F i and G i can be computed perturbatively, as series expansions in powers of (α/π), as where the superscripts (1) and ( 2) indicate the number of loops of the contributing diagrams.The analytic evaluation of the contributions of the two-loop diagrams in Figure 1 to the form factors F (2) i and G (2) i , keeping complete dependence on the masses of the internal and of the external fermions is the main result of this work.We denote these contributions as Figure 2: Integral family associated to eq. (4.1).The dashed lines denotes the external leg with momentum q µ , and q 2 = s.Straight thin lines (resp.straight thick lines) denote denominators with mass m e (resp.m i ).Wavy lines denote massless denominators.

Computation of Form Factors
The evaluation of the form factors proceeds by applying the relevant projectors, defined in eq.(2.2) to the vertices Γ (k) µ and Γ (k) .To evaluate the projectors, the Lorentz and Dirac algebra is performed in d dimensions and implemented in the Mathematica packages Package-X [53] and FeynCalc [54] independently.The result of this operation is a linear combination of scalar Feynman integrals, all members of the same integral family given by where , represented in Figure 2. We observe that D 1 , . . ., D 5 carry the momentum flowing through the diagram propagators, while D 6 and D 7 are auxiliary denominators related to irreducible scalar products.
For computational convenience, we define the integration measure of the scalar integrals in eq.(4.1) as, such that the two tadpole MIs read, The relation between the loop integral measure and the scalar integral measure, therefore is, where C ϵ is defined as with the limiting value lim ϵ→0 C ϵ = i/4 .

System of Differential Equations
The basis T i , up to ϵ rescalings, obeys a system of differential equations (sDEQ) whose matrix differential has a linear dependence on ϵ.By using the method of the Magnus/Dyson exponential matrix [38], this set is transformed into a canonical basis I with elements: which obey a sDEQ having the following canonical structure [58]: In this way, the dependence on the ϵ parameter is factorized, and the entries of the total differential matrix dA are rational in the variables x and y.The latter depend on the original variables s, m 2 i and m 2 e , through the relations with inverse whose expressions are obtained with the help of the package RationalizeRoots [59].The special form of eq.(4.8) implies that the solution can be written as a Taylor series expansion in power of ϵ, as with where γ is some regular path in the (x, y)−plane, and I (j−i) (x 0 , y 0 ) is a vector of boundary constants.
In terms of the original variables, the boundary vector corresponds to the (canonical) MIs evaluated in the limit s → 0 and m i = m e .Explicitly, the matrix in eq. ( 4.8) is given by: with and where the coefficient matrices M i , shown in Appendix A, have rational numbers as entries.
The solution to the differential equation is valid in the region where all the letters in eq.(4.14) are real and positive.In terms of the variables s, m i and m e , eq. (4.15) corresponds to the (unphysical) region Given eqs.(4.13, 4.14), the solution of eq.(4.12) can be written in terms of GPLs, defined recursively as

Boundary Conditions
The determination of the boundary constants proceeds by combining both quantitative and qualitative properties of the considered integrals: • the boundary values of I 1,2 are determined via direct integration, taking into account eq.(4.4); • the boundary constants of I 3,4 are determined by considering the subsystem formed by the first four MIs and using a different variable, z = m e /m i .The boundary constants for this subsystem are fixed at z = 0.The boundary constants for I 3,4 , expressed in the original variables x and y, are determined by matching the general solution against the abovementioned subsystem, in the equal mass limit.
• the boundary values of I 5,6,7 are computed in the limit s → 0 with m i = m e , where they are expected to vanish, due to the prefactors appearing in the definition of the canonical bases and the regularity of the MIs T 5,6,7 .
As a useful consistency check, the boundary constants were also evaluated numerically with the package AMFlow [44] at the point s = 0, m i = m e and reconstructed using the PSLQ algorithm [60].
The boundary vector takes a simple form: with We used PolyLogTools [61] for the algebraic manipulation of the GPLs, and GiNaC [41] for their numerical evaluation.The MIs were successfully compared against the numerical values provided by pySecDec [43] and AMFlow, as well as against the set of MIs presented in [31].The analytic expression of the MIs, written in terms of GPLs up to weight w = 4 are provided in the ancillary file <results.m>accompanying this article, as well as the corresponding arXiv version.

Renormalization
The diagrams depicted in Figure 1 constitute a gauge invariant subset of vertex diagrams that depend on both m i and m e , and thus the UV renormalization can be addressed independently of other contributions.The renormalization of other divergent graphs depending solely on a single mass scale m e has been performed separately in [5].The renormalized vertex functions Γ (k)ren µ and Γ (k)ren are defined by the following combination of diagrams, namely by adding two conuterterm diagrams to the unrenormalized vertices, and are computed in a two steps procedure.The first type of counterterm diagrams represents the subtraction of the one-loop sub-divergence, achieved by renormalizing the vacuum polarisation insertion in the on-shell scheme, with = Z (1) where ℓ µ = k µ 1 − p µ 1 is the momentum flowing through the insertion.The one-loop renormalization constant Z (1) 3 is implicitly defined by requiring lim The second type of counterterm diagrams, defined as, S g S 1 + Z (2) cancel the genuine two-loop residual divergences of the vertex.The two-loop renormalization constants with Z (2) j , j ∈ {V, A, S, P }, are implicitly defined by Using IBPs, the renormalization constants admit the following expressions in terms of MIs, Z = a j,1 + a j,2 + a j,3 + a j,4 , where the coefficients a j,k are not shown explicitly.After inserting the expression of the MIs, they can be expressed as Laurent series in ϵ as where v = m i /m e .We note that the difference between Z V and Z A , as well as between Z S and Z P , begins from the finite term onwards.More details on the evaluation of the renormalization constants can be found in Appendix B, where in particular the coefficients for j = V are derived as an illustrative example.

Renormalized Form Factors
For computational convenience, we introduce the rescaled form factors F (k)ren i and G (k)ren i , defined as The renormalized form factors F (k)ren i , with i ∈ {1, 2, S, P }, and G (k)ren i , with i ∈ {1, 2}, are expressed in terms of 74 GPLs up to weight w = 3, out of which 50 depending on x, with weights in the set and 24, on y, with weights in the set {−1, 0, 1, −i, i} .( Alternatively, we also provide their expression in terms of logarithms and classical polylogarithms, which turns out to be convenient for the numerical evaluation and for the series expansion of the form factors.The conversion of the GPLs into classical polylogarithms can be handled with the algorithm developed in [62] 2 .The renormalized form factors are in addition verified to be in numerical agreement with [39].Additionally, all the renormalized form factors are independently re-calculated using the variable choices and MIs presented in [31], and are found in complete agreement.
The expressions of the renormalized form factors F in terms of GPLs, and in terms of classical polylogarithms constitute the main results of this communication.Their expressions, too long to be shown here, as well as an implementation for their numerical evaluation, can be found in the ancillary file <results.m>, and <evaluator.m>,respectively accompanying this article, as well as the corresponding arXiv version.

Anomalous magnetic moment
The limit s → 0 of F (k)ren 2 corresponds to the two-loop, mass dependent contributions to the leptonic g − 2. The limit can be considered at the diagrammatic level−following the discussion in appendix B−and the form factor can be expressed in terms of T 1,2,3,4 .
For generic ϵ, its expression in terms of MIs reads with (6.6) Inserting the explicit expressions for the MIs, and considering just the finite term in the ϵ → 0 limit, eq.(6.6) yields, where z = m e /m i .Eq. (6.7) is an alternative, yet equivalent expression to the one presented in [45] 3 and revisited in [46].It can be written in terms of logarithms and dilogarithms, upon the substitutions

Conclusions
We presented the analytic evaluation of the second-order corrections to the massive form factors, coming from two-loop vertex diagrams with a vacuum polarization insertion, with exact dependence on the external and internal fermion masses, and on the squared momentum transfer.We considered vector, axial-vector, scalar, and pseudoscalar interactions in the coupling between the external fermion and the external field.The calculation was performed within the dimensional regularization scheme.Using integration-by-parts identities, the form factors were decomposed in terms of a basis of seven master integrals.The latter were evaluated by means of the differential equation method, making use of Magnus exponential matrix.The renormalized form factors were expressed in terms of generalised polylogarithms up to weight three, and in addition converted to classical polylogarithms.
The presented results can be considered as the last, missing contributions to the problem of the analytic evaluation of the second-order corrections to the massive form factors in QED and QCD, within the dimensional regularisation scheme, a problem which began to be addressed about two decades ago.The expressions of the form factors evaluated in this work can be straightforwardly applied in the context of the evaluation of the next-to-next-to-leading order virtual QED and QCD corrections to the decay of a massive neutral boson into heavy particles, or to the four (massive) fermion scattering amplitudes.

B Evaluating the Renormalization Constants
In this appendix, we describe in detail the evaluation of the renormalization constants for the vector form factors.A similar procedure has been followed in the other cases.

Renormalization Constant
The renormalization constant Z is defined implicitly by the following equation lim In order to derive its explicit expression, we expand the one-loop two point function in terms of scalar master integrals To evaluate the limit ℓ 2 → 0 it is necessary to expand the two point function, as a power series in ℓ , up to the first order as: so that, by using Taylor series expansion, we can identify By direct inspection, the ℓ 2 → 0 limit can be taken diagrammatically, (as in this example the integral is finite), which implies the value j 0 = 1.
Using IBPs, we can also evaluate the first derivative of the 2-point integral, (which corresponds to the differential equation), and take the ℓ 2 → 0 limit as follows, , (B.7) which simplifies to The latter can be read as an equation in j 1 , whose solution gives j 1 = ϵ/6m 2 i , hence fixing our Taylor series to be .
The one-loop MIs, although simple, are computed with the method of differential equations, using the change of variables as in eq.(4.9) to ensure compatibility with the unrenormalized form factors.

Renormalization Constant
We define Z (2) V through the requirement that The first contribution on the r.h.s.can be evaluated by using eq.(B.12) and eq.(B.9), giving By substituting in the relevant expansions for the master integrals, at leading order, the expression for Z (2) V takes the simple form 20)

Figure 1 :
Figure 1: Two-loop vertex diagrams with vacuum polarization insertion.The left panel corresponds to the (axial-)vector case.The right panel corresponds to the (pseudo-)scalar case.The internal fermion with mass m i is represented with a thick line, while the external fermion with mass m e is represented with a thin line.