All-multiplicity amplitudes with four massive quarks and identical-helicity gluons

We explore the on-shell recursion for tree-level scattering amplitudes with massive spinning particles. Based on the factorization structure encoded in the same way by two different recursion relations, we conjecture an all-multiplicity formula for two gauged massive particles of arbitrary spin and any number of identical-helicity gluons. Specializing to quantum chromodynamics (QCD), we solve the on-shell recursion relations in the presence of two pairs of massive quarks and an arbitrary number of identical-helicity gluons. We find closed-form expressions for the two distinct families of color-ordered four-quark amplitudes, in which all gluons comprise a single color-adjacent set. We compare the efficiency of the numerical evaluation of the two resulting analytic formulae against a numerical implementation of the off-shell Berends-Giele recursion. We find the formulae for both amplitude families to be faster for large multiplicities, while the simpler of the two is actually faster for any number of external legs. Our analytic results are provided in a computer-readable format as two ancillary files.


Introduction
Scattering amplitudes in gauge theory have been long known to allow for drastic analytic simplifications, in particular when expressed using the spinor-helicity formalism [1][2][3][4][5][6]. The prime example is the tree amplitude for n gluons, of which two carry negative helicity [7]: A(1 + , . . . , j − , . . . , l − , . . . , n + ) = i j l 4 1 2 2 3 · · · n−1|n n 1 , (1.1) written here after color ordering (see e. g. [8]). This single-term amplitude, also known as maximally helicity-violating (MHV), is the epitome of simplicity hidden inside a factorially growing avalanche of Feynman diagrams. The all-plus and one-minus helicity configurations are even simpler in that their amplitudes vanish. On the other hand, configurations with three or more minus-helicity gluons (NMHV, N 2 MHV, etc.) involve multiple terms, the number of which goes roughly like n k for an N k MHV amplitude [9,10]. Note that this is still significantly tamer than the naively expected factorial growth.
The possibility for such on-shell analytic simplifications can be attributed to • gauge redundancies, which tie together different Feynman vertices; • fields with more indices than needed to describe their on-shell spin degrees of freedom.
A natural way to sidestep such field-theoretic complications and directly target streamlined amplitude expressions is the purely on-shell method of Britto-Cachazo-Feng-Witten (BCFW) recursion [11,12], which relates higher-point amplitudes to those at lower points. Scattering amplitudes involving massive spinning particles are typically more complicated than those in the purely massless sector. However, the on-shell strategy outlined above applies to them in its entirety and can be implemented in the massive spinor-helicity formalism of Arkani-Hamed, Huang and Huang [13] (for earlier iterations see [14][15][16][17][18]).
The first all-multiplicity results in this formalism were obtained for tree-level QCD amplitudes with two massive quarks by one of the current authors [19]. In this paper, we 1. generalize an n-point result of [19] to arbitrary spinning matter, see eq. (2.4); 2. find new all-multiplicity formulae (3.10) and (3.33) for QCD amplitudes with four massive quarks and (n − 4) positive-helicity gluons.
The latter results exhibit a pattern of increasing analytic complexity, which is reminiscent of pattern seen in the massless N k MHV amplitudes. These observations lead us to 3. explore if our analytic formulae are advantageous to other methods in the context of purely numerical amplitude evaluation.
Indeed, a widely used method for computing tree amplitudes is the off-shell Berends-Giele (BG) recursion [20], which is known to be very robust and efficient in purely numerical calculations [21][22][23][24][25][26][27]. In fact, in massless QCD it was shown [27] to outperform evaluation of closed formulae [9,10] starting at the N 2 MHV level of analytic complexity. For QCD amplitudes with massive quarks, we find that the evaluation of the simple two-quark formulae, quoted from [19] in eqs. (2.3) and (4.1), is increasingly faster than the BG recursion, as the total number of particles n grows. In the four-quark case, we observe a similar situation for one of the two families of color-ordered amplitudes admitting closed-form expressions. For the other such family, however, our analytic results become numerically advantageous only at multiplicities higher than n = 15. This suggests that the complexity of the considered, relatively simple, massive external-particle configurations is already approaching the edge of what it makes sense to tackle analytically, assuming that improving the numerical efficiency of tree-amplitude generation is what one is after.
In order to facilitate the reuse of our results in analytical or numerical studies, we provide the ancillary file TreeAmpResults.txt with flexible Wolfram-friendly implementations of the two-quark amplitude expressions of [19], as well as the new four-quark formulae. Our results involve composite auxiliary spinors, a sample implementation of which can be found in the second ancillary file TreeAmpAuxSpinors.txt.

From 3 to n points for general spin
In this section we describe the on-shell approach to tree-level amplitudes with massive particles and demonstrate its effectiveness by deriving an all-multiplicity all-spin formula for gauged spinning matter.

General gauged-matter amplitude
Let us use the formalism outlined above to write the scattering amplitude for a massive quark-antiquark pair and (n − 2) positive-helicity gluons [19]: A(1 a , 2 + , 3 + , . . . , (n−1) + , n b ) = im 1 a n b [2| n−3 j=2 P 12...j p j+1 + (s 12...j − m 2 ) |n−1] where we label spinors by particle numbers instead of momenta, i. e. |p j α = |j α . By convention, we consider all momenta outgoing and denote momentum sums as P 12...j = p 1 + p 2 + . . . + p j and their Lorentz squares as s 12...j . Slashed matrices P mean either P αα or Pα α = αβ αβ P ββ depending on the spinors surrounding them. Hence the numerator in eq. (2.3) contains (n − 4) factors of (P 12...j )α γ (p j+1 ) γβ + (s 12...j − m 2 )δα β . Their order of matrix multiplication is such that j increases from left to right. Similarly to the MHV amplitude (1.1), the two-quark amplitude above consists of a single term for any number of plus-helicity gluons. Moreover, it depends on the quark spin labels a and b only through the single spinor product 1 a n b = βα |1 a α |n b β . In fact, the amplitude with two massive scalars instead of quarks, derived in [28][29][30], can be obtained by mere replacement 1 a n b → m. We can therefore make a good educated guess for the amplitude with (n − 2) plus-helicity gluons and two massive spin-s particles: 2 .
For a comprehensive exposition of the formalism, we refer the reader to [13]. Our conventions are consistent with the latest arXiv version of [19]. For completeness, in Appendix A we provide explicit spinor parametrizations in terms of momentum components, which can be used for numerical evaluation. 2 The general-spin formula (2.4) first appeared in an explicit form in the lectures [31] and implicitly already in [32] in the context of heavy-particle effective theory. Apart from matching the s = 0, 1/2 results of [19,29,30], it has also been recently derived for s = 1 in [33].
Here and below 1 a n b 2s is a shorthand for the tensor product symmetrized in the (separate) SU(2) index sets {a i } and {b i } 1 a n b 2s = 1 (a 1 n (b 1 1 a 2 n b 2 · · · 1 a 2s ) n b 2s ) . (2.5) The sign prefactors (−1) s are included in eq. (2.4) for the sake of consistency with [34]. The massive particles are understood to transform in some representation of an arbitrary gauge group. As long as we are allowed to pick the adjoint representation of SU(N c ) or U(N c ), we have a well-defined notion of color ordering. Then the permutations of eq. (2.4) may be color-dressed as required by any other group and representation [35]. For instance, in the case of U(1) we have the following amplitude with (n − 2) photons: where Q is the charge of the two charged matter particles. The spin-1/2 case simply corresponds to quantum electrodynamics. Perhaps more interestingly, the spin-1 version of eqs. (2.4) and (2.6) gives a closed-form all-multiplicity expression for the amplitude involving two W bosons emitting (n−2) plus-helicity photons, as follows from the recursive arguments below. Furthermore, the gravitational analogue of the above amplitudes, in which the massless gauge bosons are replaced by positive-helicity gravitons, can be obtained via the Kawai-Lewellen-Tye (KLT) double copy [36,37]. Indeed, in the case of two matter particles, the massive-matter extension of the double copy [34,[38][39][40] is no more complicated than in the purely massless case and therefore has a simple closed form [37]. The spin-1/2 formula (2.3) was proven [19] recursively in QCD with the four-point amplitude as the starting point. The induction step was set up using a massless BCFW shift | n−2] = |n−2] − z|n−1], | n−1 = |n−1 + z|n−2 , (2.7) which by construction preserves the on-shell conditionsp 2 n−2 =p 2 n−1 = 0 and momentum conservationp n−2 +p n−1 = p n−2 +p n−1 . Proving the spin-s formula (2.4) can be approached in the same way, starting from the corresponding four-point result of [34] and assuming the so-called "minimal-coupling" three-point amplitude [13] for the induction step: where |q is an arbitrary reference spinor. We can then effortlessly generalize the spin-1/2 derivation to the spin-s case. Instead of showing here essentially the same calculations as in [19], let us rather discuss an important remaining issue regarding higher spins, which will motivate the need for an independent rederivation of eq. (2.4).
Recall that BCFW recursion [11,12] treats a complex-shifted tree amplitude as a rational function of z and relates the residues of finite poles in z to lower-point amplitudes. Now the spin-1/2 derivation of [19] relied on the vanishing behavior of the n-point amplitude at z → ∞, which is guaranteed for QCD with massive or massless quarks [41,42]. It is also known to be true for gauge theory coupled to scalar matter [43], as well as for a spontaneously broken gauge theory [41]. In the latter case, it follows from the irrelevance of the lower-energy symmetry-breaking effects for the large-momentum behavior of the tree amplitude. so the classical arguments in pure gauge theory (see e. g. [44]) still hold for suitably chosen shifts. In the case where the unbroken part of the gauge symmetry is U(1), the sum over the orderings in eq. (2.6) must result in an amplitude for a W -boson pair emitting (n − 2) photons.
The problematic point is that tame boundary behavior is not a given for higher spins. In fact, for generic higher-spin amplitudes it must fail, because a naive BCFW derivation of the Compton scattering amplitude A(1 {a} , 2 + , 3 − , 4 {b} ) with opposite helicities is already known [34] to produce an expression with unphysical pole structure unless s ≤ 1. One could then proceed to search for additional amplitude contributions which vanish on the physical poles but subtract out the unphysical pole from the naive BCFW result [13,[45][46][47].
In the case of the all-plus formula (2.4), however, we have reasons to suspect that it constitutes a valid spin-s amplitude by itself: • First of all, the BCFW-shifted form of the all-plus amplitude behaves as O(1/z) for z → ∞. Indeed, the denominator contains a factor ofŝ 12...(n−2) − m 2 = O(z), whereas the numerator dependence on z is canceled inˆ p n−2 |n−1] = |n − 2 [n−2|n−1]. Therefore, the proposed formula is consistent with the vanishing boundary behavior.
• Moreover, the only lower-point input needed to build up the all-plus formula is the three-point amplitude (2.8) with the positive-helicity gluon. In other words, the amplitude (2.4) is sensitive exclusively to how the matter is coupled to the self-dual sector of gauge theory. Therefore, if one imagines an effective field theory, from which this formula follows, it may seem plausible that the higher-dimension effective operators in this theory (e. g. of the schematic formΦF n−2 Φ/m 2n−6 ) would behave in a particular way when evaluated on self-dual gauge fields F = i * F , which would lead to a vanishing boundary behavior of the all-plus amplitude. Constructing a fieldtheoretic argument along these lines would be very interesting but would go outside of the on-shell scope of this paper.
• Finally, we can apply a different BCFW shift and see that it is consistent with the same formula (2.4). Let us construct such an argument below.

Massive-massless BCFW shift
Although the BCFW recursion was originally established via shifting two massless momenta [12], it was soon extended to shifts involving massive momenta [30,43]. At first the corresponding massive spinors used to be defined using the massless spinor-helicity formalism with the help of an additional massless reference momentum [14][15][16]. The new massive spinor-helicity formalism [13] allows for a more elegant formulation of such shifts, as explored in [46,[48][49][50][51]. In the rest of this paper, we will be picking one massive particle j and one massless particle k and shift their on-shell spinors as follows [43,49]: where the mass m of particle j is used to ensure that z is dimensionless. In terms of momenta, this [j a , k shift translates tô (2.10) One can immediately see that the momentum conservation and on-shell conditions are preserved,p 2 j = m 2 andp 2 k = 0.
Alternative derivation of spin-s amplitude. Let us immediately employ this shift to construct an inductive argument for the spin-s amplitude (2.4). First of all, note that if we choose to apply it to the last two particles, the factor of 1 a n b 2s stays unaffected by the [n b , n − 1 shift, as does the rest of the formula except for the denominator factor n−2| n−1 . Therefore, the conjectured expression behaves as O(1/z) for z → ∞, similarly to the situation with the massless shift (2.7). This lets us continue to assume a vanishing boundary behavior in our construction of the induction step.
Considering the [n b , n − 1 shift more closely, we see that every color-ordered factorization channel which separates the shifted particles necessarily involves a purely gluonic amplitude with at most one negative helicity. The only case where such an amplitude does not vanish is when it is trivalent and evaluated on complex kinematics. Therefore, the only pole with a nonzero residue is produced in theŝ (n−2)(n−1) -channel: (2.11) To compute the left-hand side of this recursion, we only need the three-gluon MHV, given for completeness in Appendix B, as well as assume that eq. (2.4) holds for the (n − 1)-point spin-s amplitude, which constitutes our induction hypothesis. We find (2.12) On this pole, the shift variable takes the value (2.13) (Here and below we use the notation like n−2|n|n−1] = n−2| p n |n−1] interchangeably.) Hence we may identify the denominator factor n−3|P A = n−3|n−2 , and furthermore (2.14) Finally, rewriting the right-handed spinor for the complex intermediate momentumP A as we conclude that the n-point residue (2.12) coincides with the conjectured expression (2.4). We now have two different derivations of the conjectured all-multiplicity all-spin formula for gauged spinning matter, which is governed by a single gauge coupling and exhibits a vanishing boundary behavior for the used BCFW shifts. All the poles in this amplitude are physical, with correct factorization channels all the way down to the basic three-point amplitudes (2.8). In these respects, it is unique and must belong to an interesting family of effective field theories of massive higher-spin matter coupled to a gauge field. For lower spins s ≤ 1, the validity of eq. (2.4) is guaranteed by the previously established results for scalars [29] and quarks [19,30], as well as by the boundary-behavior arguments of [44,51] for vector particles in the context of a spontaneously broken gauge theory.

QCD amplitudes with two quark pairs
In this section we proceed to the main computations of this paper: four-quark amplitudes with any number of gluons of positive helicity. We consider the more general case of two distinctly flavored quark pairs. The identical-flavor amplitude may always be obtained by setting the two masses to be equal and subtracting two relabelings: where we have color-coded the two quark flavors by red and blue. Moreover, we concentrate on the color-ordered versions of the tree amplitudes in question. There are multiple possible color decompositions [38,52] which can be used to assemble the full color-dressed amplitude from different color orderings. Let us first specialize to orderings of the form A(1, 2, {3, 4} ¡ σ(5, . . . , n)) where σ ∈ S n−4 is an arbitrary permutation and ¡ is the shuffle product. In this setting, {3, 4} ¡ σ(5, . . . , n) amounts to all permutations of {3, 4, 5, . . . , n}, in which the relative ordering of the quark labels 3 and 4 stays fixed. These orderings constitute Melia's amplitude basis [53,54] with respect to the Kleiss-Kuijf (KK) relations [55]. This basis and the corresponding color decomposition [38,56] allows us to fix the relative position of the four quarks and avoid gluon insertions between quarks 1 and 2. Furthermore, the color-kinematic information additionally provided by (the mass-and flavor-extended version of) the Bern-Carrasco-Johansson relations [38,57,58] lets us consider only orderings of the form A(1, 2, 3, σ(4, 5, . . . , n)).
Color-ordered amplitudes contain poles only in subsets of consecutively ordered momenta. This is very helpful for solving on-shell recursions. Indeed, the BCFW relations receive contributions from z-dependent poles, with the z-shifted particles on the different sides of each factorization channel. Consequently, if two adjacent particles are shifted in an n-point ordered amplitude, there are only n − 3 factorization channels to consider.
Note that some ordered amplitudes involving massive and massless particles, such as A(1, 2, 3, 4, 5), lack a pair of adjacent massless particles, so the purely massless shift of the type (2.7) cannot be applied to them. This is why we find the massive-massless shift (2.9) especially useful for computing massive four-quark amplitudes. Since we only consider gluons of positive helicity, we shift the left-handed massless spinors, e. g.
Note that in the limit where the quark 4 is taken to be massless its SU(2) labels d = 1, 2 may be associated to the positive and negative helicities, respectively. Therefore, the above [4 d , 5 shift actually corresponds to two massless BCFW shifts: [4 + , 5 + and [4 − , 5 + . As regards the [4 − , 5 + shift, in QCD it is easy to see [59] that every Feynman diagram vanishes at infinity no worse than the gluon's polarization vectorε + 5 = O(1/z). Showing the vanishing boundary behavior under the massless [4 + , 5 + shift is already more involved as the naive power counting estimates the amplitude to be O(1) [30,43]. For completeness, in Appendix C we provide an argument based on the QCD Feynman rules showing that for the general massive-massless shift [j a , k the leading O(1) terms vanish, thus ensuring the applicability of on-shell recursion.

Gluon insertions between like-flavored quarks
We have found the simplest family of color orderings to be the amplitudes of the form A(1 a , 2 b , 3 c , 5 + , . . . , n + , 4 d ), where quarks 1 and 2 carry mass M , while quarks 3 and 4 have another flavor and mass m. For concreteness, we set 1 and 3 to be outgoing quarks and 2 and 4 to be antiquarks. This can be easily tweaked, since reversing a fermionic line in the color-ordered amplitude (without changing the ordering) simply results in an overall minus sign, as shown in Appendix D. We will be using this way of flipping quark lines to our convenience in the rest of the paper. Moreover, we will henceforth underscore the quarks and bar the antiquarks only when their choice is not otherwise obvious, which will unclutter the amplitude notation.

(3.8)
This formula is explicitly local because we have obtained it from a single pole contribution.
n-point amplitude. We have chosen the above form of the five-point amplitude because it can be seen as a special case of an all-multiplicity formula that we have found. To write it in a compact form, let us introduce a shorthand notation for the propagator denominators, such as D 3i...j = P 2 3i...j − m 2 , as well as the following auxiliary massless spinors: where the matrix factors are ordered by k increasing from left to right. These objects are slight generalizations of the spinors appearing in the numerator of the two-quark formula (2.3), although now they involve quark momenta p 3 and p 4 . Our all-multiplicity can then be written as The last term in the sum (3.10b) contains the auxiliary spinor |d n n+1 ], which is not defined in by eq. (3.9b), because its lower index is larger than the upper index. We supplement our definition by the following additional requirement: The only other spinor contraction appearing in the i = n contribution, n|P 35...n |d n n+1 ], is irrelevant, as it comes with a vanishing prefactor D 4 = p 2 4 − m 2 . In this way, we were able to integrate this contribution into the general sum instead of spelling it out separately.
Note that both |c j i ] and |d j i ] transform under the massless little group of the spinor |j]. For instance, at five points the only such spinor is |d 5 5 ] = |5], because the sum (3.10b) does not yet contribute. It is easy to see that the remaining term (3.10a) produces precisely eq. (3.8). Having thus ensured the base case, we may proceed to prove the complete n-point formula (3.10) by induction.
Inductive proof. We choose to use the BCFW shift [4 d , n + , for which each recursion step has exactly two contributions: (3.12) The first diagram B n factorizes into the four-quark seed amplitude (B.3) and the previously obtained [19] n-point amplitude (2.3) with two external quarks. Thus we immediately get (3.13) where we have taken care to separate the shifted factor n − 1|n from the rest of the denominator. We can also recognize the auxiliary spinor [c 5 n−1 | appearing in the unshifted numerator factor. Similarly to the five-point calculation, we need to insert the relevant pole values for the complex kinematics (determined byD 35...n = 0), for which we obtain (3.14) We also need to eliminate the internal momentumP B :  After slight rearrangement, this can be seen to match the i = n contribution in the sum (3.10b) of the n-point formula.
The second residue C n is where we really need an inductive argument, as it factorizes into an (n − 1)-point amplitude of the same type that we aim to compute: By the induction hypothesis, we may express the right-hand side of the factorization using the formula (3.10), in which the complex kinematics is determined by the poleŝ (n−1)n = 0: , The internal momentumP C may then be decomposed into , (3.19) up to a rescaling, where we have recognized the first appearance of the auxiliary spinor |d n n−1 ], defined in eq. (3.9b). The non-recursive part of the residue can be identified with the factor −i Combining it together with the prefactor in the (n − 1)-point amplitude, we can already observe the formation of the n-point prefactor Let us now inspect how the recursion affects the terms of the formula (3.10) in the curly brackets. We immediately notice that almost all of the featured momentum sums involve bothp 4 andP C and so can be easily reduced to an unshifted form. Namely, where the arrows correspond to going from the n-point expressions to the shifted (n − 1)point expressions on the right-hand side of the residue (3.17). From a similar examination of the auxiliary spinors, we see that [c 5 i−1 | stay entirely unshifted, whereas |d n i ] become Applying these replacements to the first contribution (3.10a) in the curly brackets, we find where we have used the identity After combining eq. (3.24) with the prefactor (3.21) and slightly rearranging the last line, we retrieve the precise form of the first term (3.10a) in the n-point expression. In other words, this term is converted into itself under the recursion (3.12).
We now turn to the part the shifted (n − 1)-point amplitude which involves the sum (3.10b). First, we isolate the last term from the rest of the sum, in which we use n|4|d n n+1 ] → P C |4|d P C P C +1 ] = 1 in line with the supplementary definition (3.11), which simplifies this contribution. Together with the main sum, it gives where we have also taken into account the transition rules (3.22) and (3.23). The only quantity that still needs to be evaluated on the shift kinematics iŝ Here we have replaced |n] = |d n n ] so as to expose the similarity to the summands in the main sum above. We can therefore integrate the last three lines of eq. (3.26) into this sum as its i = n − 1 contribution. If we multiply this sum by the prefactor (3.21) and add the i = n term arising from by the first residue B n , the result will coincide with the corresponding sum in the all-multiplicity amplitude (3.10). This concludes our proof.

Gluon insertions between distinctly flavored quarks
Let us now proceed to presenting results for the other family of color orderings of the form A(1 a , 2 b , 3 c , 4 d , 5 + , . . . , n + ), where all gluons are inserted between quarks 1 and 4 that carry different flavors. We will follow a similar path as the previous subsection.

5-point amplitude.
At five points, the [4 d , 5 shift (3.2) now implies two BCFW residues, the second of which contains two helicity configurations: From this recursion, the amplitude can be derived in the following form: Note that the contributions B n , C n , and E n have already appeared at five points, whereas F n is a new residue appearing for n > 5, and it is where the main recursion happens. Let us proceed directly to the final expression, which we again write in terms of auxiliary spinors that originate from the numerator of the two-quark formula (2.3): Here the momentum sums like P 3...i are always understood in the "clockwise" fashion, namely P 3...i = p 3 + p 4 + . . . + p i−1 + p i . The proof of eq. (3.33) is very similar to that in the previous subsection, and we leave it to Appendix E.2.

Gluons in between both distinctly and like-flavored quarks
Recall that in order to construct the full color-dressed amplitude with four quarks one needs [38] color-ordered amplitudes of the general form A (1, 2, 3, σ(4, 5, . . . , n)), This of course includes the cases where gluons are inserted on both sides of quark 4. Due to the complexity of the recursive structure, we refrain from computing the closed formula in this ordering. Alternatively, here we present a BCFW shift to extract the numerical result of the amplitude. In this ordering, we found it most efficient to use the shift [3 c , 5 + . The factorization channels are   The last term in eq. (3.34) involves a lower-point version of the amplitude on the righthand side, with gluons 5 and 6 replaced by a single gluon with a complex momentum. This is where the recursion step takes place, unless k = 5, in which case this term vanishes, and the amplitude is obtained from the previously computed amplitudes.

Numerical evaluation
In this section we analyze our analytic results from the point of view of the time needed to evaluate them numerically on a computer.
The main benchmark for us will be our own implementation of the off-shell Berends-Giele (BG) recursion [20], also in C++. This method essentially provides an inductive realization of conventional Feynman diagrams and is therefore very flexible and relatively easy to implement. In particular, its validity does not depend on such properties of the theory under consideration as the boundary behavior for large complex-shifted momenta. Furthermore, in the case of massless gauge theory it was found [21,27] to be more efficient than the on-shell BCFW recursion, in the case of generic helicities of the external particles. For these reasons, the BG recursion is the method of choice for numerical evaluation of tree-level amplitudes for phenomenological purposes [22][23][24][25][26].
It is perhaps somewhat counter-intuitive that closed analytic formulae even for massless QCD amplitudes are not always able to outperform the BG recursion. Indeed, [27] observed that the evaluation time of analytic formulae for massless QCD amplitudes [10] (which were obtained by solving [9] the on-shell BCFW recursion in N = 4 supersymmetric Yang-Mills theory) grows faster than the time needed to run a numerically efficient implementation of the BG recursion already at the N 2 MHV level of analytic complexity. Recall that massless amplitudes possess a "complexity-peeling" property, which is reflected by the fact that the color-ordered N k MHV amplitude formulae involve k nested sums.
This "complexity-peeling" property is scrambled in the massive case. Indeed, in the massive spinor-helicity formalism a single amplitude e. g. with 4 massive quarks constitutes a 2×2×2×2 matrix, each element of which can be associated with a massless-quark helicity amplitude in the high-energy limit. So it is natural that the complexity of the analytic formula in the massive case can be no simpler than that of the most complicated massless helicity configuration that it incorporates. In fact, we observe that the massive formulae are generally even more complicated than this naive expectation, as most easily seen in the two-quark case. For instance, the "one-minus" amplitude [19] A(1 a , 3 − , 4 + , . . . , n + , ..k | n−2 j=k P 13...j p j+1 + D 13...j |n] s 3...k n−1 j=k D 13...j 3|1|P 3...k |k 3|1|P 3...k |k+1 × 1 a 2 b 3|1|P 3...k |3 + 1 a 3 2 b 3 s 3...k (4.1) involves a linearly growing number of terms. However, all but one term contain the quark mass, so only the first term survives the massless limit. It gives rise to 2 × 2 = 4 helicity amplitudes, two out of which vanish and the other two constitute the quark counterparts of the Parke-Taylor formula (1.1).

Evaluation timings.
In view of such a non-trivial analytic-complexity behavior of our massive analytic results, it is interesting how long it takes to evaluate them numerically. Figures 1 and 2 show our timing results for QCD amplitudes with two and four quarks, respectively. The depicted timing data were obtained on a 2017 laptop with a Intel Core i7-7500U processor and Ubuntu Linux 20.04. Each amplitude was computed using its explicit BG for qgq All-Plus (α ≈ 3.7) Solved qgq All-Plus formula obtained by solving an on-shell BCFW recursion and via an off-shell BG recursion. During each measurement, the amplitudes were evaluated repeatedly using double-precision floating-point numbers on random kinematic phase-space points for approximately 10 seconds. The time periods needed for each evaluation were averaged, while their fluctuations resulted in error bars that are too small to be visible in the figures behind the data-point markers. As a sanity check, we have also verified that three other laptops (equipped with macOS 11, 12 and Windows 11) produce qualitatively identical results.
The two-quark amplitudes evaluated in figure 1 correspond to the "all-plus" formula (2.3) involving only positive-helicity gluons and to the "one-minus" formula (4.1). Both four-quark formulae (3.10) and (3.33) in figure 2 correspond to all gluons having positive helicity but in different positions with respect to the two quark flavors. We have implemented these formulae using C++. Both the analytic formulae and the BG recursion were implemented for cross-checking purposes and were not fully optimized. The BG recursion is implemented in a bottom-up approach, used in [60], and is similar to what is described in [27], without the need for the cashing of off-shell currents that was employed in [21]. The analytic formulae were implemented as given, without low-level optimization. As only the compiler optimization was used (albeit with the "aggressive" flag -Ofast), it is likely that some additional gain of perhaps up to 50% in performance is achievable for both approaches. However, we do not expect this to affect neither the large-n scaling of the computational complexity nor the comparison between algorithms to a significant degree.
There is an important difference in the way the analytic formulae and the BG recursion are implemented. That is, although the latter uses the same massive spinor parametrization as the former, all the external spin labels are fixed whenever the BG recursion routine is called, since it is natural for an off-shell method to treat massive and massless particles on the same footing. This is in contrast to the on-shell analytic formulae, which fix the  massless helicities but incorporate all massive spin degrees of freedom in one go -and which were implemented accordingly. The formulae were therefore evaluated once for each kinematic phase-space point, whereas the BG recursion was called 4 and 16 times per point in the two-and four-quark cases, respectively.
Our observations based on the evaluation timing data are the following.
• Evaluating the one-term all-plus expression (2.3) is much faster than the BG recursion for any number of particles, as clearly seen from figure 1. The observed dependence on the total number of particles n is roughly linear, which is consistent with the growing complexity of the numerator and denominator in the formula.
• The one-minus formula (4.1) is also faster than the BG recursion, but not by as large a margin as the all-plus expression, see figure 1. The observed dependence on the number of particles is roughly quadratic, with different machines giving the slope α to be in the 1.9-2.1 range. This is consistent with the growing number of terms in combination with their increasing complexity.
• The behavior of our simplest four-quark formula (3.10) with respect to the corresponding four-quark BG recursion, as seen in figure 2, is qualitatively similar to that of the two-quark one-minus amplitude. Namely, the analytic formula is faster and has a milder growth curve than the BG recursion.
• The more complicated four-quark formula (3.33) starts by taking roughly the same time to evaluate as the BG recursion but due to a milder growth curve takes a lead from 9 external gluons (n = 13), see figure 2. On other machines, we have observed that this intersection point may move (as high as n = 22) but continues to exist.  • The BG timings have a very mild dependence on the particle configuration, as is perhaps best seen in the unified figure 3. Indeed, the power-law slopes α for all four recursive amplitude evaluations are in the region 3.4-3.7. In fact, once we accounted for the factor of four difference in the number of evaluations in the two-and fourquark cases, we could observe that the four-quark amplitudes with gluon insertions between like-flavored quarks (labeled in the figures as "QQqgq") are measurably faster to evaluate for the BG recursion than the two-quark amplitudes, for the same number of particles. This is because replacing two gluons by a quark pair may reduce the number of possible factorization channels in an amplitude and hence accelerate the recursion. In the "QQqgq" configuration, this indeed happens because the color ordering entirely prevents the external gluons from coupling to one of the quark lines.
• On the contrary, the analytic formulae depend very strongly on the external particle and helicity configuration, see figure 3. In particular, replacing gluons with quarks increases the evaluation time by a factor increasingly greater than four, at least in the case of identical-helicity gluons.
• However, the slopes of the four-quark analytics seem to be gentler than expected. Formula (3.10) contains an increasing number of terms, each of which grows in complexity with the number of particles n, so it might already seem slightly surprising that its dependence on n appears to be milder than quadratic in view of the observed α = 1.9. This effect, however, is not very robust, measuring this slope on other machines produces results in the range 1.7-2.0, which does include the expected behavior. The more complicated amplitude (3.33), however, even contains a double sum, which suggests a cubic dependence for large n, as opposed to the observed α = 2.1. This is definitely significantly tamer than 3 and is corroborated by our alternative datasets, which place α in the range 1.7-2.1. Our suspicion is that either our amplitude formula (3.33) responds particularly well to compiler optimization or 30 particles is simply not enough for these amplitude formulae to exhibit their true large-n behavior .

Conclusions
In this paper, we have obtained new analytic all-multiplicity results for gauge-theory amplitudes with massive matter. First, we have substantiated the conjecture [31] for the amplitude (2.4) involving two massive spin-s particles and an arbitrary number of positive-helicity gluons by showing that it is consistent with two different types of BCFW shifts. This meshes particularly well with similar results [32] in heavy-particle effective theory [32,61]. Then we have derived the closed formulae (3.10) and (3.33) for two families of QCD amplitudes involving two quark pairs and an arbitrary number of positive-helicity gluons. The corresponding amplitudes with all negative-helicity gluons are of course a collateral result. Namely, they may be obtained simply by flipping the chirality of all SL(2, C) spinors, i. e. by exchanging angle bra and ket spinors with square ket and bra spinors, respectively. Finally, we have explored how much time our analytic formulae take to be evaluated on a computer as opposed to the widely used method of numerical off-shell recursion [20]. We do not claim our implementations of either side of such a comparison to have been sufficiently optimized to be production-ready, but we still believe them to be adequate for our purposes of qualitative comparison. We find the large-multiplicity behavior of our amplitude formulae to still be significantly milder than that of the BG recursion, although the more involved of our two four-quark formulae needs moderately high multiplicities (13-22 for our implementations) to overtake the BG evaluation. This means that solving the BCFW recursion in search of analytic formulae in more complicated cases than considered here may hardly be motivated by the need for numerically faster tree amplitudes. Indeed, on-shell recursion relations with more terms (or terms containing sums in themselves) are likely to produce analytic results with three or more nested summations, which will perform worse than the numerical off-shell recursion, as previously observed in massless QCD [27]. 3 We hope, however, that the all-multiplicity formulae discussed in this paper will provide new analytic data for the exploration of new structures of QCD. For instance, it was the investigation of the geometric origin of the spurious poles [62] in massless gauge-theory amplitudes that has led to breakthroughs in our understanding of the all-loop integrand structure of N = 4 supersymmetric Yang-Mills theory [63][64][65][66], such as its Grassmannian/Amplituhedron geometry. In the massive-quark case at hand, one could observe the presence of spurious poles already in the one-minus two-quark amplitude formula [19], given in eq. (4.1). These poles appear in the form of composite spinor products 3|P |Q|k , where P and Q are massive momenta, which is very similar to what happens in massless NMHV amplitudes. Now the four-quark formulae (3.10) and (3.33) exhibit spurious poles of the new schematic form D P k|Q|c] + D Q k|P |c] , where |c] is a composite spinor. These poles seem less reminiscent of their massless counterparts, and we wonder whether they could be understood e. g. in terms of twistor [62], ambitwistor geometry [67,68] and scattering equations [69][70][71][72][73][74][75][76]. We hope that finding such an understanding of QCD amplitudes involving massive particles will be facilitated by our new explicit analytic results.
to counter the uncertainty of treating the square root. At a glance, the exponents may seem reducible to θ(i √ −1), however, this would spoil the behavior under momentum reversal. Massive spinors can be generically parametrized as The Moreover, the identification det |p a α = det [p a |α = m > 0 can be spoiled by the determinants being negative due to the square-root uncertainty again. In that case, we restore their positivity by introducing an additional factor of i into |p a and −i into [p a |.
The above parametrization is invalid if p + = 0, in which case we use As before, additional factors of ±i should be introduced to ensure positivity of the determinants. The above two parametrizations fail when P = 0, in which case we have p ± = ±p 3 . If it is not zero, we parametrize the massive spinors as In the case where E = p + , another option one can adopt is (A.10) A singular kinematics that invalidates the above parametrizations is P = p ± = 0. This means that either p ⊥ = 0 or p ⊥ = 0, and the momentum reduces to p µ = {p 0 , p 1 , ±ip 1 , 0}.
Then we use the following parametrization Finally, for a massive particle at rest, the momentum is p µ = {p 0 , 0, 0, 0}, in which case we may adopt the parametrization that aligns the spin vectors with the z-axis, (A.12)

B Seed amplitudes for recursion
Here we provide all the basic ingredients sufficient for the QCD computations in this paper. The three-gluon MHV amplitude is given by the Parke-Taylor formula The amplitudes with two quarks and a gluon of either helicity are [19] A(1 a , (B. 2) The last building block involves two pairs of distinctly flavored massive quarks. It is most easily computed from a single color-ordered Feynman diagram: -25 - Here we have used the following construction of the Dirac spinors from the Weyl spinors: We have also abused the notation by writingū p = − p| + [p| and v p = −|p + |p], which is allowed as long as one only contracts indices of the same chirality.

C Boundary behavior argument
Here we prove the vanishing boundary behavior of QCD amplitudes under the massive quark-gluonic shift [j a , k , as given by eq. (2.9), provided that the gluon k has positive helicity. The argument is similar the one in [42] for the massless shift. In an arbitrary tree-level Feynman diagram, we consider the unique path that the shift momentum zr µ follows from the antiquark j to the gluon k: 4 Here K + L free Lorentz indices and the implicit fermion index are understood to be contracted with other parts of the diagram, which only depend on the sump j +p k = p j +p k . Let us examine how the leading term for large z is constructed. First of all, there are K + L different propagator denominators affected by the shift, which altogether provide an O(1/z K+L ) behavior. The two shifted external wavefunctions compensate each other: Moreover, on the left-hand side of the diagram, the term leading in z is obtained from K copies of the shift momentum zr µ from each momentum-dependent gluonic vertex, while the gluon propagator numerators are taken to be constant in the Feynman gauge. Finally, on the right-hand side of the diagram the vertices are constant, whereas the highest power of z is obtained by picking up L copies of zr from each fermion propagator numerator. Therefore, it seems like the diagram (C.1) should behave as O(1).
Before we can see that the leading O(1) term actually vanishes, we note that it involves an odd Dirac-algebra element γ ν zr γ ν 1 zr · · · zr γ ν L , where the index ν is contracted with the z-dependent gluonic tree on the left-hand side of the diagram. This expression may be rewritten in terms of eight basic matrices {γ µ , γ µ γ 5 } µ=0,1,2,3 via the property Thus the fermionic part of the O(1) contribution must involve γ µ |k] = γ µ γ 5 |k] = σ µ |k]. Now we are ready to consider the faith of the Lorentz vector indices. Altogether, there are K + L free indices, which may be distributed among K + L copies of zr µ , one polarization vectorε µ k+ and one Dirac-algebra structure γ µ |k]. We see that at least some index contractions are required to make this happen. Since we have both the metric and Levi-Civita tensors in the mix, the contractions allowed in the leading contribution are up to K+L copies ,ε µ k+ ∝ q|σ µ |k], γ µ |k] , (C. 4) both of which reduce the number of free indices by two. However, all such combinations are exactly zero. Indeed, r 2 = r ·ε + k = r |k] = ε + k |k] = 0, and the most involved check is for where we have reversed the property (C.3) to remove the Levi-Civita contraction. We may thus conclude that individual Feynman diagrams behave as O(1/z) under the shift (2.9).

D Fermionic line reversal
In this appendix we prove that reversing an arrow of a fermionic line in a color-ordered amplitude results in an overall factor of −1. In an arbitrary Feynman diagram, consider the part that directly involves a fermionic line: (D.1) Here the Lorentz indices on the gluonic lines (which for concreteness we put to the right of the fermion line) should be contracted with other parts of the diagram, each carrying a momentum flow P k . According to the color-ordered Feynman rules and the Dirac spinors in eq. (B.4), the relevant part of the above diagram is then where P f,k = p f + k j=1 P j . We omit the factors of i and the propagator denominators. We now wish to transpose the diagram by putting the spinorsf in the front and f in the end. Note that the number of gamma matrices between these spinors vary between n − 1 and 2n − 1, the latter being the only term surviving the massless limit. First, let us for concreteness consider this term, which has an odd number of gamma (and therefore sigma) matrices. Then we can transpose it directly by using the properties Here transposition is understood to switch between the two kinds of sigma matrices, as well as change the order of multiplication, e. g.
So it is clear that the massless term of eq. (D.2) can be rewritten as For an even number of gamma matrices, the internal Weyl index structure is different, and an additional sign appears due to internal contractions of the type αγ βγ = −δ β α : Note that such combinations of gamma matrices occur in the Feynman diagram along with an odd power of mass m, whereas the odd combinations are multiplied by an even power of m. Therefore, we may absorb all the minuses induced by transposition into the mass terms of the propagators. Hence we rewrite the diagram (D.2) as Here in the second line we have written the propagator numerators in the usual form, as well as have accounted for the flipped color-ordered Feynman rules, which amounted to pulling out (−1) 2n−1 = −1 to the front. In this way, we may conclude that we have examined all the consequences of flipping a fermion line and arrived at precisely minus the fermion-flipped version of the diagram (D.1), as promised. Finally, note that the same sign flip would occur even if the gluonic lines in the diagram were arranged differently.

E.1 5-point amplitude with gluon between distinctly flavored quarks
Here we present the derivation of the five-point amplitude (3.29 where the latter identity is true for any choice of |q due to gauge invariance but is perhaps most obvious for |q = |4|5]. Similarly, we will be fixing the reference spinors to our convenience without further comments in the rest of the appendix. Plugging eqs. (E.2) and (E.3) into eq. (E.1), we can show that the residue B 5 coincides with the part of (3.29) which is inside the big square brackets. The contributions C 5 and E 5 are evaluated at the same pole, whereŝ 34 = 0 and hence The contribution C 5 involves an amplitude with two quarks and with two positive-helicity gluons, which is a special case of (2.3). Combining with the three-point amplitude in eq. (B.2), we find The necessary shifted quantities here are . (E. 6) We emphasize that the mass m in the denominator of eq. (E.5) is canceled by the numerator [3 c4d ], so there is no pole in the massless limit. Using that 5 |P 12 |5] = −s 12 on this pole, one can now easily verify that C 5 corresponds to the first line of eq. (3.29).
Lastly, E 5 needs the four-point amplitude with two opposite-helicity gluons, which can be found e. g. in [19]: With this we can write the contribution E 5 as  5 4 ] to s 34 5|4|d 5 4 ] = s 34 and using a couple of Schouten identities. For general n-point amplitudes, we are going to compute every diagram in (3.30). We will show that B n , C n and E n correspond to the i = 5 terms in the sums (3.33b), (3.33c) and (3.33d), respectively. The rest of terms will be recursively generated by F n .
The computation of B n is similar to that of the five-point case. The crucial new step is to recognize n−2 j=5 P (−P B )5...j p j+1 + D (−P B )5...j |n] = |a n 6 ] (E. 9) in the numerator of the all-plus two-quark amplitude (2.4). Gluing the latter with the four-quark amplitude (B.3), we obtain The next two contributions C n , E n come from the same factorization limit in which the two quark pairs are entirely separated. Gluing one of the three-point amplitudes in eq. (B.2) and the (n − 1)-point two-quark amplitude (2.4) with all gluon helicities positive, we obtain This corresponds to the i = 5 term in the sum (3.33c), in which 5|P 34 |3|P 34 |d 5 4 ] is identified with s 34 , as follows from the definition (3.32).
The contribution E n comes from the same factorization channel as C n , but the intermediate gluon has opposite helicity. This means that E n involves the two-quark amplitude with one negative-helicity gluon and the rest of the gluon helicities being positive. Such an amplitude was computed in [19] and is given by eq. (4.1). In it, the minus-helicity gluon is color-adjacent to an outgoing quark, whereas in E n it must stand next to an outgoing antiquark. We flip the overall sign to reverse the arrow of the fermionic line, as described in Appendix D, and combine it with the first three-point amplitude in eq. (B.2). This gives M k|k+1 P C |2|P 3...k |P C P C |P 3...k |a n k+1 ] 12 P C |2|P 3...k |P C − 1P C 2P C s 3...k s 3...k n−1 j=k D 2...j 5 6 P C |2|P 3...k |k P C |2|P 3 Recalling that 5|P 34 |3|P 34 |d 5 4 ] = s 34 , we can match it to the i = 5 term in the sum (3.33d). Finally, the last contribution F n involves a lower-point amplitude of the same type that we are calculating: so we need the inductive hypothesis. Note that the following steps are largely analogous to those following eq. (3.17) in the proof of the ordered amplitude (3.10). This is because its structure is very similar to eq. (E.24). Now the shifted kinematics in this channel is determined byŝ 56  where we see the formation of the n-point prefactor.
Now we consider the effects of the recursion on the terms of the amplitude formula (3.33) in the curly brackets, which we assume to hold at (n − 1) points by the induction hypothesis. Almost all the momentum sums appearing in the formula involve bothp 4 andP F and so remain unshifted: where the arrows indicate the transition from the general n-point expression (3.33) to the specific shifted (n−1)-point amplitude A(1 a , 2  so this term can be written as the i = 6 term in the remaining sum. In comparison with the terms (3.33c) in the n-point formula, we seem still to be off by a factor of 6|4|5]/D 45 , but it is exactly canceled by the prefactor (E.28). We have thus verified that the inductive residue F n converts the sum (3.33c) into itself but lacking the i = 5 term, which we have already retrieved from C n . Similar inductive arguments can be seen to hold for the contribution (3.33a), as well as for the sums (3.33b) and (3.33d), which are complemented by B n and E n , respectively. This concludes the proof of the n-point formula (3.33).