A three-body form factor at sub-leading power in the high-energy limit: planar contributions

The Standard Model (SM) of particle physics contain a couple of massive particles at the electroweak scale, including the Higgs boson, the top quark and the electroweak gauge bosons. An important quest of particle physics is to study the properties of these particles at the energy and luminosity frontiers with high precision, and probe new physics beyond the SM. This requires high precision theoretical predictions for the production processes involving these massive particles. However, the calculation of scattering amplitudes involving massive particles is generally much more difficult than the calculation of massless scattering amplitudes. On the other hand, due to the high energies of the Large Hadron Collider (LHC) and future colliders, the masses of the particles are often small compared to other kinematic invariants in the scattering processes. With such scale hierarchies, the perturbative scattering amplitudes and cross sections can develop large logarithms involving the ratios between the small masses and the large kinematic invariants in the high energy limit.

The standard way to deal with scale hierarchies is the method of factorization. Such factorization can be organized into different orders in an expansion parameter $ \lambda \sim m^2/E^2 $, where m denotes the low scale of the small masses, while E denotes the high scale of other kinematic invariants. The leading power (LP) corresponds to order $ \lambda^0 $. The factorization formula at the LP has been well-understood [1−10]. It has been shown in [10] that a massive scattering amplitude in the high energy limit can be factorized into a massless amplitude, a soft function and several collinear functions (one for each external leg). Such a factorization formula can be used to predict the structure of the large logarithms of the form $ \ln(m/E) $ at higher orders in perturbation theory, and also allows the resummation of these logarithms to all orders.

Given the importance of the factorization approach, it is essential to investigate the behavior of the massive amplitudes at sub-leading powers in λ. On one hand, this allows us to learn about the typical size of power corrections, and therefore provides an uncertainty estimate of calculations based on the LP factorization formula. On the other hand, once we understand the factorization structure at sub-leading powers, we can extrapolate the high-energy approximations to intermediate energy ranges, and eventually combine them with low-energy approximations based on threshold/soft factorization. This will lead to an adequate description of the scattering amplitudes in the whole phase space.

During the past years, studies of small-mass factorization at the next-to-leading power (NLP) in λ have emerged in the literature [11−28], either within the framework of soft-collinear effective theory (SCET) [29−35] or based on diagrammatic analysis. The analysis of sub-leading soft emissions has been carried out in [12]. The construction of an easy-to-use helicity operator basis has been studied in [13−15]. Anomalous dimensions of sub-leading power N-jet operators have been calculated in [17, 19]. Based on power counting and region analysis, the NLP factorization formula in the small-mass limit has been proposed for the Yukawa theory [18] and for quantum electrodynamics (QED) [22, 27]. An important ingredient in the QED factorization formula, the NLP jet function, has recently been computed in [28]. However, the NLP factorization formula has not been generically proven.

An important validation of the factorization formula is to compare its fixed-order expansion with a direct computation of some massive scattering amplitudes or form factors in the small-mass limit. The $ 1 \to 2 $ massive quark form factor has been computed at two loops [36−38], and three-loop calculations are in progress [7, 39−49]. However, $ 1 \to 2 $ kinematics does not capture all possible structure at the two-loop order and beyond. In this paper, we initiate a study of a two-loop $ 1 \to 3 $ massive form factor at sub-leading powers. At two loops, non-trivial correlations can only occur among at most 3 external particles. Therefore, the structure extracted from $ 1 \to 3 $ form factors is generic enough to be applied to other form factors and scattering amplitudes. To this end, our ultimate goal is to study the $ Q\bar{Q}g $ form factor in quantum chromodynamics (QCD), where Q represents a massive quark, and g is the gluon. Due to the complexity of this problem, we will begin with a simpler one in QED, i.e., the $ e^+e^-\gamma $ form factor.

The scattering amplitudes involving massive electrons in QED are phenomenologically interesting on their own. In the future plans of high-energy physics experiments, various high-energy and high-luminosity $ e^+e^- $ colliders have been proposed [50−56]. On these machines, it is important to understand the standard QED processes as precise as possible. These processes provide crucial inputs for the calibration of the detectors, and for a better understanding of beam parameters. The studies in this work are therefore useful for these applications.

The remainder of the paper is organized as follows: we introduce the notation in Section 2, and describe the method to determine independent coefficients in the small-mass expansion in Section 3. In Section 4, we demonstrate the solution of the differential equations satisfied by the independent coefficients, and discuss the final results. We briefly summarize in Section 5 where an outlook for future works is also described.

We consider the $ 1 \to 3 $ QED process

The independent kinematic variables are

where m is the mass of the electron. For convenience, we define the following dimensionless variables

The results for the tree-level amplitude $ {\cal{M}}^{(0)} $ and the one-loop amplitude $ {\cal{M}}^{(1)} $ have been well-established [57, 58]. In this work, we are concerned with the two-loop amplitude $ {\cal{M}}^{(2)} $. We generate the relevant two-loop diagrams and the corresponding amplitudes using $\mathrm{QGRAF}$ [59]. These diagrams can be categorized into eight integral families: five planar and three non-planar. We only consider the planar families in this work. The amplitudes containing Lorentz and Dirac indices can be decomposed into an appropriate basis with scalar coefficients (form factors). For simplicity, in this work we demonstrate our method using the the interference term in the squared-amplitude $ {\cal{F}}^{(2)} = \langle {\cal{M}}^{(0)} | {\cal{M}}^{(2)} \rangle $ where the indices are summed over. We use $\mathrm{FeynCalc}$ [60−63] and $\mathrm{FORM}$ [64, 65] to manipulate the expressions, and express the results in terms of two-loop scalar integrals. These scalar integrals can be expressed as

where $ a = \sum_i a_i $, and we have multiplied an appropriate power of $ (-s_{123}) $ to make the scalar integrals dimensionless. Each set $ \{D_i\} $ defines an integral family. There are 5 planar families involved in this work. The corresponding $ \{D_i\} $ are given by

Note that an integral family may contain several topologies, where a topology is defined according to which $ D_i $'s appear in the denominator of Eq. (4) (i.e., $ a_i > 0 $). The topologies belonging to each of the 5 families are illustrated in Fig. A1 of Appendix A.

Figure A1.  Two-loop planar topologies. Thick lines represent propagators with mass m or external legs with $ p^2 = m^2 $. Thin lines represent massless propagators or external legs. Red lines represent external legs with $ p^2 = s_{123} $.

For each family, we use $ \mathrm{Kira} $ [66, 67] to reduce the scalar integrals to a set of MIs $ {\cal{I}}_i(x,y,z) $ by solving integration-by-parts (IBP) relations. The contribution from this family to the two-loop squared-amplitude can then be decomposed as

where μ is the renormalization scale, while the coefficients $ {\cal{A}}_{i} $ are rational functions, and can be easily expanded in the limit of small x. The remaining task is then to compute the master integrals in the high-energy limit $ x \to 0 $, while keeping the exact dependence on the variables y and z. To illustrate our method, we focus on the family (5c). There are 123 MIs in this family. In the high-energy limit, these MIs admit the asymptotic expansion

The purpose of this work then reduces to the calculation of the coefficients as a function of y and z. Note that the two-loop integrals considered in this work have no soft-collinear overlapped divergences, and the minimal power of $ \epsilon $ is $ -2 $. The integrals also have no power-like singularities in the limit $ m \to 0 $, and therefore the minimal power of x is $ 0 $. The maximal power of $ \log(x) $, $ n_{3,\max} $, depends on the value of $ n_1 $ and $ n_2 $. In practice, we will truncate the series in $ \epsilon $ and x for the purpose of calculating the squared-amplitude $ {\cal{F}}^{(2)} $ up to $ x^1 $ and $ \epsilon^0 $. This corresponds to the next-to-leading power (NLP) in the high-energy expansion. Nevertheless, our method can be used to compute higher power terms as well.

To compute the coefficients, we adopt the method of differential equations [68−70]. From the IBP reduction, we can construct the differential equations of the MIs $ {\cal{I}}_i $ with respect to the variable $ t \in \{x,y,z\} $:

where $ {\cal{P}}^{(t)}_{i,j} $ are rational functions of $ \epsilon $, x, y and z. The differential equations with respect to x lead to relations among the coefficients. Therefore, we only need to compute a set of independent coefficients, which is conceptually similar to the master integrals. After determining the set of independent coefficients, we can employ their differential equations with respect to y and z to obtain their analytic expressions.

We apply the expansion (7) into the differential equations (8) with $ t = x $. We note that

That is, taking derivative with respect to x always decreases the power of x, but may or may not decrease the power of $ \log(x) $. On the other hand, the matrix elements $ {\cal{P}}^{(x)}_{i,j} $ may also contain poles at $ x = 0 $. Therefore, the differential equations with respect to x lead to linear relations among the coefficients $ {\cal{C}}_{i,n_1,n_2,n_3}(y,z) $. From these linear relations, we can determine a set of independent coefficients. We will refer to these independent coefficients as "master coefficients" (MCs), and express the remaining coefficients as linear combinations of the MCs.

Before solving the linear relations (there are an infinite number of them), we need to constrain the highest values of the integers $ n_1 $ and $ n_2 $ for each i, i.e, the maximal powers of $ \epsilon $ and x in the expansion of the master integral $ {\cal{I}}_i $ (note that the maximal power of $ \log(x) $ is naturally determined given the values of $ n_1 $ and $ n_2 $). The first constraint comes from the required orders in the expansion of the squared-amplitude $ {\cal{F}}^{(2)} $. As we have mentioned, we truncate the expansion in x up to the NLP, i.e., order $ x^1 $. In addition, we truncate the expansion in $ \epsilon $ up to order $ \epsilon^0 $. These requirements impose constraints on $ n_1 $ and $ n_2 $ for each master integral $ {\cal{I}}_i $. However, we find that these constraints are too tight for our purpose: the MCs determined from these constrained linear relations are not closed under differentiation. The reason is that there are cancellations among different topologies in a family when their contributions go into Eq. (6).

In order to obtain a closed system of differential equations, we need to slightly loosen the constraints. For that purpose, we expand the coefficients in the differential equations as well:

The lowest values of k and l are then crucial for the determination of the highest expansion orders. As an example, we consider the differential equations with respect to x. The constraints from y and z can be similarly studied. Plugging Eqs. (7) and (10) into (8), we find that

This is the master equation for relations among expansion coefficients. For this system to be closed, it is necessary that all coefficients contributing to a given $ (n_1,n_2,n_3) $ are incorporated in the above equation.

Let's focus on the equation with $ i=11 $, which arises from the differential equation of $ {\cal{I}}_{11} $ with respect to x. In this case, one has $ j \in \{1,8,11\} $ on the right-hand side of Eq. (11). For the NLP accuracy of the squared amplitude, one requires both $ {\cal{I}}_{11} $ and $ {\cal{I}}_{1} $ to be expanded to order $ x^3 $, since the minimal power of x in $ {\cal{A}}_1 $ and $ {\cal{A}}_{11} $ is $ -2 $, which means that we need to consider the equation with $ n_2 = 3 $. However, the minimal value of l for $ i=11 $, $ j=1 $ is $ l=-2 $. Therefore, the coefficient with $ n_2-1-l=4 $ will appear on the right-hand side of the equation. For this reason, we need to expand $ {\cal{I}}_{1} $ up to order $ x^4 $.

Considerations similar to the above also applied to the constrains on $ n_1 $, i.e., the expansion orders in $ \epsilon $. We need to perform such analysis for each sector, and update the constraints until all equations are closed. As a final outcome of such analyses, we provide the number of MCs for each integral family in Table 1.

Integral families	$ F_1 $	$ F_2 $	$ F_3 $	$ F_4 $	$ F_5 $
Number of MCs	247	359	383	158	199

Table 1.  The number of MCs in each integral family.

After obtaining a system of linear equations of the coefficients, we can solve the equations to express all coefficients in terms of the MCs. Before doing that, we need to decide a set of rules for choosing the MCs. This is similar to the rules for choosing the MIs in Feynman integral reduction. In practice, this requires to determine between two coefficients which one is preferred over the other.

The choice of MCs will greatly influence the process of solving their differential equations with respect to y and z. We would like them to have as simple analytic expressions as possible. From experience, this then corresponds to selecting coefficients with smaller $ n_1 $, smaller $ n_2 $ and larger $ n_3 $. Coefficients in lower sectors are also preferred over those in higher sectors.

After deciding the preference, we can solve the linear system using the "user-defined system" functionality of $ \mathrm{Kira} $. We can then construct the system of differential equations of the MCs. For that we need to re-express the derivatives of the MCs in terms of the MCs again. This step can be accomplished by $ \mathrm{Kira} $ as well. For example, we write the differential equation of a MC with respect to y as

The above equation can be provided to $ \mathrm{Kira} $ along with other linear relations. One then obtains the left-hand side expressed in terms of MCs. Finally, we arrive at the differential equations

where $ I=(i,n_1,n_2,n_3) $ and $ J=(j,m_1,m_2,m_3) $ denote two subscript sequences that label the MCs.

After constructing the differential equations, we now want to find the analytic solutions for the MCs. Ideally, the solution procedure can be greatly simplified if the connection matrices $ {\cal{A}}^{(y)}_{I,J}(y,z) $ and $ {\cal{A}}^{(z)}_{I,J}(y,z) $ take a strictly triangular form. In that case the solutions can be easily written in terms of iterated integrals. This is similar in spirit to the idea of canonical differential equations [71]. However, the existing tools [72−76] for finding canonical differential equations do not easily adapt to the current case. Fortunately, the connection matrices are already in a block-triangular form with at most $ 3 \times 3 $ blocks. These can be iteratively converted to a strictly triangular form, as will be demonstrated in the following.

We first consider the simple cases where the derivatives of a MC only depend on itself and the already-solved MCs:

where $ {\cal{G}}_{y} $ and $ {\cal{G}}_{z} $ come from the already-solved MCs. If the coefficients $ {\cal{A}}_{y} $ and $ {\cal{A}}_{z} $ are both zero, the above equations can be easily solved by direct integration. This corresponds to a strictly triangular system discussed earlier. We now consider the case where $ {\cal{A}}_{y} $ and $ {\cal{A}}_{z} $ are nonzero.

For convenience, we introduce the following functions:

Multiplying the MC by the two functions, the differential equations Eq. (14) are transformed into:

One can see that only already-solved MCs appear on the right-hand side of the above two equations. We would like to find a transformation function $ {\cal{T}}(y,z) $ such that

This is an example of a strictly triangular system of differential equations, where the derivatives of the transformed MCs only depend on already-solved ones. Comparing Eq. (17) with Eq. (16), one can see that

where the function $ f_z(z) $ satisfies the differential equation

and similar for $ f_y(y) $. The fact that the right-hand side of the above equation is independent of y follows from the compatibility condition of Eq. (14):

which further follows from the linear-independence of $ {\cal{C}}(y,z) $ with respect to $ {\cal{G}}_y $, $ {\cal{G}}_z $ and their derivatives. From the above analysis, we can write

Building on this result, we get general solution to the differential equaiton (17):

where $ {\cal{C}}_{y}(y) $ is a function that only depends on y, and similarly for $ {\cal{C}}_{z}(z) $. These two functions are related by the second equal sign in the above formula, and are only fixed once the arbitrary constant terms of the indefinite integrals are chosen. Using the fact that $ {\cal{C}}_{y}(y) $ is independent of z, we can derive the differential equation of $ {\cal{C}}_{z}(z) $ with respect to z:

which allows us to determine the analytic expression of $ {\cal{C}}_{z} $ (as well as $ {\cal{C}}_{y} $) up to a constant of integration. This constant should be fixed by a boundary condition, to be discussed later.

The above procedure of transforming the differential equations to a strictly triangular form can be naturally extended to $ 2 \times 2 $ blocks and $ 3 \times 3 $ blocks, as will be discussed in the following subsections.

We now consider the $ 2 \times 2 $ blocks, where the differential equations take the form

where $ t \in \{y, z\} $, $ {\cal{A}}^{(t)}_{ij} $'s are rational functions of y and z, and $ {\cal{G}}^{(t)}_i $'s come from already-solved MCs. To solve the equations, we would like to find a linear combination $ {\cal{C}}' = b_1 {\cal{C}}_1 + b_2 {\cal{C}}_2 $ with $ b_1 $ and $ b_2 $ being rational functions of y and z, such that the derivatives of $ {\cal{C}}' $ only depend on itself. Taking the derivatives, one has

We therefore require

Once we find a suitable pair of $ b_1 $ and $ b_2 $, the differential equation of $ {\cal{C}}' $ becomes

The above equation resembles the one studied in Section 4.1. We can employ the method there to get rid of the term $ {\cal{A}}'_t {\cal{C}}' $, and effectively transform the differential equation to a strictly triangular form. After obtaining the solution of $ {\cal{C}}' $, we can use it to derive analytic expressions of $ {\cal{C}}_1 $ and $ {\cal{C}}_2 $. Taking $ {\cal{C}}_1 $ as an example, we have

Since $ {\cal{C}}' $ is known, the above equation can again be solved using the method of Section 4.1. The expression of $ {\cal{C}}_2 $ can then be obtained from those of $ {\cal{C}}' $ and $ {\cal{C}}_1 $.

From the above discussions, it is clear that the task boils down to finding a particular solution of $ b_1 $ and $ b_2 $ satisfying Eq. (26). Note that Eq. (26) is actually a differential equation for the ratio $ b_2/b_1 $, which can be rewritten as

This is the so-called Riccati equation. There are two ways to find a solution.

Method one: brute-force solution. For the cases we encounter in this work, the Riccati equation can be solved with the help of $\mathrm{Mathematica}$. As an example, we may need to solve the system with

For convenience, we can set $ b_1 = 1 $, and the Riccati equation with respect to y reads

Using $\mathrm{Mathematica}$, we can obtain a general solution of the above equation, which involves arbitrary functions of z. We require that the solution is rational, and also satisfies the corresponding equation with respect to z. This give rises to a particular solution

In practice, this method is sometimes slow, and we may resort to the second method.

Method two: solution by ansatz. Since the equation only constrains $ b_2/b_1 $, we can assume that both $ b_1 $ and $ b_2 $ are polynomials of y and z with integer coefficients. We can then make the ansatz

where the coefficients $ c_{ij}^{(1)} $ and $ c_{ij}^{(2)} $, as well as the degree n, are yet to be determined.

We start with a small value of n, and substitute the ansatz into Eq. (26). This then leads to a system of quadratic equations of the coefficients. If the system of equations only allows a trivial solution (i.e., all coefficients being zero), it means that the value of n is too small. We then increase the value of n until find a non-trivial solution.

For the example considered above, we find that a non-trivial solution can be obtained with $ n=3 $:

which is equivalent to the result of the first method.

Similar to the $ 2\times2 $ cases, the differential equations for $ 3\times3 $ blocks have the form

The method for dealing with this system is analogous. First, we need to construct a linear combination $ {\cal{C}}' = b_1 {\cal{C}}_1 + b_2 {\cal{C}}_2 + b_3 {\cal{C}}_3 $, whose differential equations only involve itself. This boils down to the following equation

This equation can be rewritten as

One can observe that the two unknowns $ b_2/b_1 $ and $ b_3/b_1 $ are still coupled by the two equations, and this system is challenging to solve by brute force. Fortunately, the second method introduced in the previous subsection still works, and offers an algorithmic way to tackle this kind of problems. We make polynomial ansatz for $ b_1 $, $ b_2 $ and $ b_3 $, and solve for the coefficients from Eq. (36). After that, we can solve for $ {\cal{C}}' $ using the method of Section 4.1. We are then left with a $ 2\times2 $ system of $ {\cal{C}}_1 $ and $ {\cal{C}}_2 $, which can be solved following the strategy of Section 4.2.

Let's again demonstrate our approach using an example, where

A particular solution can be found by using the ansatz with $ n=2 $, and is given by

This give rise to a decoupled differential equation

We don't go into the details of the remaining steps here.

In this work, we do not encounter blocks more than $ 3 \times 3 $. But as one can see, the method introduced above can be applied to more complicated cases as well. This proceeds in an iterative way: one first solves for a linear combination whose differential equation is decoupled, which then leads to a simpler system.

To obtain the solution to the system of differential equations, we still need to determine the boundary conditions for each coefficient. For that we use $\mathrm{AMFlow}$ [77] to numerically compute the small-mass expansion of all MIs at a specific kinematic point. During this process, all variables (except $ m^2 $ but including $ \epsilon $), are set to exact values (rational numbers). For each value of $ \epsilon $, the output has the following structure

where $ \alpha_i(\epsilon) $ is a linear function of $ \epsilon $, $ w_{ij}(\epsilon) $ depends on $ \epsilon $ and the kinematic variables (which we have suppressed), and $ j_{\max} $ is specified according the required expansion order (in $ m^2 $) of the MIs. By varying $ \epsilon $ with fixed kinematic variables, we can reconstruct the coefficients in $ \alpha_i(\epsilon) $, as well as the function $ w_{ij}(\epsilon) $ as a power series in $ \epsilon $:

To determine the coefficients $ a_{ijk} $, we need to sample at least $ (k_{\max}-k_{\min}+1) $ different values of $ \epsilon $. Note that the precision of the reconstructed $ a_{ijk} $ also depends on the number of samples. Therefore, we need to choose an appropriate number according to the required expansion order (in $ \epsilon $) of the MIs and the required precision for the coefficients. In practice, this number is about $ 50 $.

With the reconstructed coefficients, we can reexpand Eq. (41) in $ \epsilon $, and obtain the small-mass expansion of the MIs at the chosen kinematic point in the form

The minimal powers of $ \epsilon $ and x, and the maximal powers of $ \log(x) $ for a given pair $ (n_1,n_2) $ can be read off from the above expansion. Note that $ f_{i,n_1,n_2,n_3} $ is just the coefficient $ {\cal{C}}_{i,n_1,n_2,n_3}(y,z) $ at a particular kinematic point, and the general solution of $ {\cal{C}}_{i,n_1,n_2,n_3}(y,z) $ has already been obtained from the differential equations. Such a solution can be written in terms of generalized polylogarithms (GPLs) with the help of $\mathrm{PolyLogTools}$ [78], up to some unknown constants. These constants can then be reconstructed as transcendental numbers by the PSLQ algorithm, using the high-precision numbers (about 100 decimal digits) of $ f_{i,n_1,n_2,n_3} $. The transcendental basis is generated from the elements $ \{\pi,\ln(2),\zeta_2,\zeta_3,\zeta_4,\text{Li}_4(1/2),\zeta_5\} $ up to weight 5. It is worth mentioning that the DEs may demand higher order coefficients in $ \epsilon $, as discussed in Sec. 3.1. That's the reason why weight-5 functions and constants appear, despite the fact that two-loop amplitudes up to order $ \epsilon^0 $ only require weight-4 at most. Note that the form of the analytic expressions can be different in different kinematic regions. We concretely work in the Euclidean region where all planar MIs have no imaginary part. The expressions in other kinematic regions can be obtained via analytic continuation. Since our results are written in terms of GPLs, their analytic structure is very well understood. The branch cuts of GPLs are uniquely fixed by their symbols, and whenever one analytic continues across a branch cut, an extra term proportional to $ \pm i\pi $ is generated, where the sign of $ i\pi $ depends on the direction of going across the cut. See, e.g., [79−81] for more detailed discussions.

With the analytic expressions for the MIs, we can combine them to obtain the planar contributions to the squared-amplitude $ {\cal{F}}^{(2)} $ up to order $ x^1 $ and $ \epsilon^0 $. It can be written in terms of GPLs up to weight $ 5 $, with the following symbol letters:

There are two additional letters $ \{1-y,\,2-y-z\} $ appearing in the MIs. But they cancel out in the final expression for the squared-amplitude. Note that the original massive amplitude involves elliptic integrals, but the expanded amplitude can be expressed by GPLs up to the NLP. The result for the squared-amplitude is very compact, with a size of about 400 KB. We attach the expression as an electronic file to this paper.

We have performed several sanity checks on our calculation. We have applied our method to the one-loop amplitude where the complete result can be easily obtained. Upon expansion in the high-energy limit, the complete result agrees exactly with our calculation. We have also numerically computed the MIs at several different kinematic points other than the chosen boundary point, and find that the results agree with the outcomes from our analytic expressions. In particular, we have verified kinematic points beyond the Euclidean region, and find that the imaginary parts agree as well. These checks demonstrate the reliability of our method.

In this work, we initiate a study of the sub-leading power contributions to the multi-parton massive form factors in the high-energy limit, where the parton masses are much smaller than their energies. These form factors provide crucial information to formulate and validate sub-leading power factorization theorems for generic scattering amplitudes. Such factorization theorems are important to resum the mass logarithms beyond the leading power, and to generate approximate results for scattering amplitudes at higher loops.

While the two-loop massive quark form factors are available, the $ 1 \to 2 $ kinematics is not generic enough to be extended to multi-parton scattering amplitudes. Therefore, we focus on $ 1 \to 3 $ form factors where two-loop multi-parton correlations can be studied. To establish the calculational techniques, we start with the two-loop planar contributions to a QED form factor $ \gamma^* \to e^+e^-\gamma $ in this work. Although only two-particle correlations are present in this case, it is enough to demonstrate that our method is capable to obtain analytic expressions of the relevant loop integrals as a small-mass expansion.

Our calculations are based on the differential equations satisfied by the expansion coefficients of the master integrals. They can be derived up to an arbitrary order in the power expansion using the differential equations satisfied by the master integrals. The differential equations with respect to the mass m are then utilized to derive linear relations among the expansion coefficients. Similar to the IBP reduction, these linear relations can be solve to express all coefficients in terms of a finite set of master coefficients. The differential equations of the master coefficients with respect to momenta invariants are then employed to solve for their analytic expressions. The solutions are expressed by GPLs, which are then combined and lead to a compact analytic result for the planar contributions to the form factor.

There are several obvious follow-ups to this work in sight. A rigorous verification of factorization at NLP necessitates the inclusion of non-planar contributions to the amplitudes. Our method can be readily applied to non-planar contributions to this form factor. The main obstacle in such a calculation lies in the IBP reduction. We have attempted with the state-of-art reduction tools such as $ \mathrm{Kira} $, $\mathrm{NeatIBP}$ [82, 83] and $\mathrm{Blade}$ [84], but have not succeeded within a reasonable amount of time. It will be necessary to develop more efficient reduction methods for such cutting-edge calculations. Another obvious generalization is to consider the $ Q\bar{Q}g $ form factor in QCD, where genuine three-parton correlations can occur. This involves more integral families, but the method established in this work can still be straightforwardly applied. Finally, the differential equations for the master coefficients are solved by brute-force in this work. It will be interesting to investigate whether these differential equations can be casted into a canonical form which can be solved directly as iterated integrals. We leave these investigations to future works.

As mentioned in the main text, we have 5 integral families, which contains 19 planar topologies in total, as shown in Fig. A1. The MIs in each integral family are listed as follows: