Polarised Amplitudes and Soft-Virtual Cross Sections for bb̄ → ZH at NNLO in QCD

Production of the Higgs boson, H in association with a massive vector boson, V , i.e., the V H process, plays an important role in the explorations of Higgs physics at the Large Hadron Collider, both for a precise study of Higgs’ Standard Model couplings and for probing New Physics. In this publication we present the two-loop corrections in massless quantum chromodynamics (QCD) to the amplitude of the Higgs production associated with a Z boson via the bottom quark-antiquark annihilation channel with a non-vanishing bottom-quark Yukawa coupling, which is a necessary ingredient of the full next-to-nextto-leading-order QCD corrections to the V H process in the five-flavour scheme. The computation is performed by projecting the D-dimensional scattering amplitude directly onto an appropriate set of Lorentz structures related to the linear polarisation states of the Z boson. We provide analytic expressions of the complete set of renormalised polarised amplitudes in terms of polylogarithms of maximum weight four. To give an estimation of the size of contributions from amplitudes considered in this work, we compute numerically the resulting cross sections under the soft-virtual approximation. We also take the opportunity to make a dedicated discussion regarding an interesting subtlety appearing in the conventional form factor decomposition of amplitudes involving axial currents regularised in D dimensions. ar X iv :1 91 0. 06 34 7v 1 [ he pph ] 1 4 O ct 2 01 9


Introduction
Ever since the discovery of the 125 GeV Higgs boson at the Large Hadron Collider [1,2] (LHC), the detailed investigation of its dynamical properties, i.e., how it interacts with (other) known fundamental particles, remains among the major research topics of the current and future particle physics programs. The interactions of this Higgs boson explored so far are in accord with the predictions of the Standard Model (SM), and considerable improvement of the experimental precision is expected with the future high luminosity LHC program. (See, for instance, ref. [3] and references therein for a brief overview.) One of the recent achievements regarding the investigation of the 125 GeV Higgs boson's couplings to fermions was the direct observation of its decay to a pair of bottom quarks by the ATLAS and CMS experiments [4,5]. The main contributing channels to this important experimental result are from production processes in which the Higgs boson is produced in association with a W or Z boson, known as the V H process, and the associated electroweak vector boson is typically chosen to be reconstructed via its leptonic decays. The presence of the vector boson in the final state in addition to the Higgs boson is crucial to substantially reduce the SM backgrounds, for instance by requiring a large transverse momentum of the associated vector boson [6]. Additional selection criteria are also imposed to enrich the signal V H events over backgrounds eventually to a manageable level to allow for such an experimental observation [4,5]. Accordingly, on the theoretical side, it is thus very desirable to have a precise knowledge about the V H process at hadron colliders, especially to meet the foreseeable precision requirements from future experiments for studies of the 125 GeV Higgs boson (as well as potentially non-standard Higgs bosons) with ever more details.
Given the aforementioned phenomenological importance of V H productions, there have been many computations available in the literature on this subject aiming to improve theoretical predictions as much as possible. Main production channels at the LHC include the quark-induced and gluon-induced V H processes. The focus of this article is a part of the b-quark-induced ZH process that involves a non-vanishing Yukawa coupling λ b between the b quark and the Higgs boson. At the tree level there are three contributing diagrams for the b-quark-induced ZH process 1 , as shown in figure 1. The diagram (A) gives an Figure 1. Feynman diagrams at leading order s-channel contribution with the same structure as that of the Drell-Yan production, which has been studied extensively in QCD up to order O(α 2 s ) in refs. [7][8][9][10][11][12][13][14][15] and to O(α 3 s ) in refs. [16][17][18][19]. The presence of this contribution is independent of λ b as the Higgs boson is radiated by bremsstrahlung of the Z boson, which makes this channel particularly valuable for studying V V H vertex. In diagrams (B) and (C) the Higgs boson is coupled directly to the b quark via its Yukawa interaction. Since the Z boson appears in this non-Drell-Yan type diagrams only as an external on-shell leg, there is no explicit gauge-fixing parameter involved and this part of the contributions is manifestly gauge invariant. If one keeps the b quark kinematically massless (m b = 0), these two types of contributions do not interfere and hence can be treated separately. This is because the chirality of the b quark line is preserved in diagram (A) but flipped exactly once in diagrams (B) and (C).
Starting from O(α 2 s ) in QCD a class of two-loop diagrams appears contributing to the b-quark initiated ZH process in which the Higgs couples directly to a closed top-quark loop. This was studied in refs. [12,12,20] by making use of asymptotic expansions in the heavy-top limit. In addition, at order α 2 s the gluon-fusion induced ZH production also opens up, which actually contributes considerably due to the rather large gluon luminosity at the LHC [20][21][22]. While these Drell-Yan type contributions and top-quark loop induced corrections (in the heavy-top limit) have already been studied extensively in the literature, all of which are independent of λ b , the work presented in this article focuses on the order α 2 s QCD corrections to the leading-order diagrams (B) and (C). We keep the b quark kinematically massless in our analytic computation of two-loop amplitudes involved. Keeping λ b = 0 while m b = 0 amounts to retaining just the leading contribution of the result with full b quark mass dependence, an approximation similar to what was adopted in a series of impressive works on H → bb [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37] 2 . Recently, by keeping the non-zero bottom Yukawa coupling, the di-Higgs production cross section at next-to-next-to-leading order (NNLO) in QCD at threshold is also achieved in ref. [39].
In the SM, quarks acquire their masses through Yukawa interactions with one Higgs doublet which develops a non-vanishing vacuum expectation value (vev ) leading to electroweak symmetry breaking. Consequently the quark Yukawa coupling strengths are simply proportional to their respective masses divided by vev. However, in certain beyond SM scenarios with more than one Higgs doublet, e.g., the Minimal Supersymmetric SM [40] (MSSM), the bottom-quark Yukawa coupling can be enhanced w.r.t. the top-quark Yukawa coupling in the large tanβ region, where tanβ is the ratio of vev s of up-and down-type Higgs fields in the Higgs sector of the MSSM. Such scenarios motivate the detailed investigation of the b-quark-induced ZH process that is proportional to λ b . To estimate the size of contributions from the channels considered in this work we compute numerically their cross sections under the soft-virtual approximation [28,41] in the five-flavour scheme.
Apart from being motivated by its potential phenomenological relevance in Higgs physics, there are also a few theoretical questions in the computation of the non-Drell-Yan type bb → ZH scattering amplitude, which we think are of interest. We address subtleties appearing in the conventional form factor decomposition of loop amplitudes involving an axial current in D dimensions (with a non-anticommuting prescription for γ 5 ): whether we need to include all evanescent Lorentz structures 3 to end up with correct results in computations made in D-dimensions; whether the particular regularisation prescription implied by projectors prescribed recently in ref. [42] remains unitary at higher orders once applied to this scattering process, etc. The two-loop correction to the non-Drell-Yan type diagrams in figure 1 provides a non-trivial case to investigate these issues. We will address these questions in detail through the work presented in this article.
The article is organised as follows. In the next section, we set up the kinematics and specify a few basic aspects of the computation. In section 3 prescriptions used for constructing projectors are exposed in detail. In particular, the subsection 3.1 is devoted to a self-contained discussion of a set of projectors that correspond directly to amplitudes in the linear polarisation basis. The Lorentz structures needed in the conventional form factor decomposition and the corresponding projectors are given in subsection 3.2 with a dedicated discussion of an interesting subtlety appearing in the conventional form factor decomposition of amplitudes involving axial currents regularised in D dimensions. Section 4 gives details of the ultraviolet renormalisation prescription we adopted as well as that of the verification of the axial quantum anomaly relation in our computational setup. In section 5, we outline the flowchart, especially the tool chain, employed to accomplish the computation of amplitudes presented through this work. We verify in section 6 the universal infrared divergences appearing in our ultraviolet renormalised amplitudes, with emphasis on the regularisation-scheme independence of the four-dimensional finite remainders obtained in this way. In particular, we address the aforementioned questions. To estimate the size of the contributions arising from channels considered in this work we compute numerically the cross sections in the soft-virtual approximation in section 7. We conclude in section 8.

Preliminaries
We consider the production of a scalar Higgs boson, H, in association with a massive vector boson, Z, through bottom quark anti-quark annihilation b(p 1 ) +b(p 2 ) → Z(q 1 ) + H(q 2 ) . (2.1) Here b(b) denotes the bottom quark (anti-quark). The quantity within the parenthesis represents the momentum of the corresponding particle satisfying on-shell conditions p 2 i = 0, q 2 1 = m 2 z , q 2 2 = m 2 h , where m z and m h are the mass of the Z and Higgs boson, respectively. The Mandelstam variables are defined as The physical region of the phase space is bounded by t u = q 2 1 q 2 2 such that it satisfies where κ is the Källén function defined by κ(s, q 2 1 , q 2 2 ) ≡ s 2 + q 4 1 + q 4 2 − 2(sq 2 1 + q 2 1 q 2 2 + sq 2 2 ) .

(2.4)
As already mentioned in the introduction, we keep a non-zero Yukawa coupling but only for the b quark, which otherwise is set to be kinematically massless. We compute through this article only QCD corrections to the non-Drell-Yan type diagrams in figure 1 that depends on λ b , in n f = 5 flavour massless QCD. The Z-boson interacts with all massless quarks through the respective vector and axial couplings. Any multi-loop calculation involving axial coupling in dimensional regularisation [43,44] (DR) faces the problem of defining the inherently 4-dimensional objects, Dirac's γ 5 (and Levi-Civita symbol µνρσ ), properly in D-dimensions. In this article, we follow the definition of γ 5 in dimensional regularisation which was introduced by 't Hooft-Veltman [43] and Breitenlohner-Maison [45] However, the γ 5 defined through the above equation no longer fully anti-commutes with the D-dimensional γ µ , which has profound consequences in computations involving axial currents in D-dimensions. In particular, the aforementioned definition of γ 5 affects the ultraviolet renormalisation non-trivially [46][47][48], to be addressed in details in section 4. In order to define the Hermitian axial current correctly we need to symmetrise it [47,49] before using the definition (2.5) By combining (2.5) and (2.6), we obtain [47,49] which is used in D-dimensions for our calculation. The contraction of pairs of Levi-Civita symbols appearing in the calculation is made through where all the indices carried by space-time metric tensors on the right hand side are (by definition) considered in D dimensions [50]. The symbol [ ] around the indices represents the anti-symmetric combination. In the next section, we discuss the projector method that is adopted to compute helicity amplitudes of the scattering process (2.1).

The Prescription of Projectors
In this work helicity amplitudes of the process (2.1) are obtained by projecting the Ddimensional amplitude in its defining Feynman-diagrammatic representation directly onto a minimal set of D-dimensional projection operators following the approach proposed in [42]. Since there has not been much discussion of this approach in the literature, we will first provide a quick recap of the essentials involved, while the reader is referred to ref. [42] for more details.

Projectors for Linearly Polarised bbZH Amplitudes
Polarised amplitudes carry bookkeeping external polarisation state vectors. What typically stops one from viewing the product of these state vectors as a projector defined in the usual sense is simply the fact that these external state vectors are not explicitly given solely in terms of external momenta and/or algebraic constants. These are essentially the defining criteria of the usual projectors in order to end up with projections that are Lorentz invariant and dependent only on external kinematics. The polarisation projectors as prescribed in [42] are based on the momentum basis representations of external state vectors, and all their open Lorentz indices are by definition taken to be D-dimensional to facilitate a uniform projection with just one dimensionality D=g µ µ . Essentially, the momentum basis representations of polarisation state vectors allow us to find a Lorentz covariant representation of the tensor products of external particles' state vectors (for both bosons and fermions, massless or massive) solely in terms of external momenta and algebraic constants, such as the metric tensor, Levi-Civita symbol and Dirac matrices. As all their open Lorentz indices are promoted to be D-dimensional, no dimensional splitting is ever introduced for loop momenta and/or Lorentz indices of inner vector fields. In this way, the contraction of Lorentz indices is made commutable with loop integration. It is also owing to this fact that despite being different from the the conventional dimensional regularisation scheme [51] (CDR), all UV renormalisation constants and integrated IR subtraction coefficients needed to complete a full fixed order computation of polarised physical observables can be directly recycled from results obtained within CDR. The bare (helicity) amplitudes resulting from this computational scheme are in general different from those defined in main-stream regularisation schemes, such as CDR [51], 't Hooft-Veltman (HV) [43], Dimensional-Reduction (DRED) [52][53][54], and Four-Dimensional-Helicity (FDH) schemes [55,56], while the properly defined finite remainders will of course be the same 4 . All these points will be verified explicitly for the scattering amplitudes of (2.1) up to two-loop order presented in this publication, in particular in section 6.
There are a number of reasons responsible for generating differences in the bare amplitudes as regularised in this way from their counterparts defined in the aforementioned main-stream dimensional regularisation schemes, as explained in ref. [42]. In particular, in our scheme the number of polarisation degrees-of-freedom (d.o.f) of an external particle is always equal to the number of its physical d.o.f in 4 dimensions while, however, all open Lorentz indices are by definition set to be D-dimensional to facilitate technically a uniform projection in D dimensions. Notably, all bookkeeping Levi-Civita symbols ubiquitously appearing in the momentum basis representations of certain transversal polarisation states are by definition manipulated according to (2.8) with all resulting space-time metric tensors set to be D-dimensional. Notice that in both construction and application of these D-dimensional projectors onto amplitudes, there is no need to appeal to their Lorentz tensor decomposition representations first (and hence there is no question of whether or not the amplitude reconstructed from form factor decomposition is faithful in D dimensions).
The amplitude of (2.1), being multi-linear in the state vectors of the external particles to all loop orders in perturbative calculations, can be schematically parameterised as The symbol Γ µ ≡ Γ µ vec + Γ µ axi denotes a matrix in the Dirac spinor space with one open Lorentz index µ which may be carried by either the elementary Dirac matrix γ µ or one of the external momenta involved (such as q µ 1 ). It consists of contributions from the vector and axial coupling of the Z boson, denoted respectively by Γ µ vec and Γ µ axi . As emphasized in the introduction, this work deals with the non-Drell-Yan type contributions to (2.1) in massless QCD where a non-vanishing Yukawa coupling λ b is retained. Because of the Yukawa coupling vertex on the external b quark line, all non-vanishing amplitudes of this gauge invariant class of diagrams involve two external massless spinors, u(p 1 ) and v(p 2 ), with opposite chirality, because there is one flip of helicity between these two external b quark spinors. As a consequence, the power of elementary Dirac matrices in Γ µ sandwiched between the two external massless b quark spinors must be even in order to have a nonvanishing matrix element between such a pair of spinors (with opposite chirality).
According to ref. [42], we construct the following list of projectors, which in their bookkeeping forms read as where the N s = 1, N p = γ 5 , and ε µ j with j = T, Y, L denote the three linear polarisation eigenstates of Z(q 1 ) identified in the center-of-mass reference frame of the collision. To be more specific about this, we choose ε µ T to be the transversal polarisation within the scattering plane determined by the three linearly independent external momenta, p 1 , p 2 , q 1 . The other transversal polarisation ε µ Y is orthogonal to p 1 , p 2 , and q 1 , and is constructed using the Levi-Civita symbol. The third physical polarisation state of the Z boson, its longitudinal polarisation denoted by the vector ε µ L , has its spatial part aligned with its own momentum q 1 .
The momentum basis representations of ε µ j thus defined in terms of p 1 , p 2 , q 1 can be determined in the following way. We first write down a Lorentz covariant parameterisation ansatz for the ε µ j and then solve the orthogonality and normalisation conditions of linear polarisation state vectors for the linear decomposition coefficients. Once we have established a definite Lorentz covariant decomposition form in 4 dimensions solely in terms of external momenta and kinematic invariants, this form will be used as the definition of the corresponding polarisation state vector in D dimensions. Following this line, it is rather straightforward to arrive at the following result: where the squares of the normalisation factors are We start with the primitive bookkeeping form of the projectors in (3.2), and then substitute (3.3) and subsequently simplify the expressions using 4-dimensional Lorentz and Dirac algebra (and 4-dimensional equations of motion) as much as possible. In this way one ends up with one definite form of projectors to be used in D-dimensional calculations (see discussion in the beginning of section 3.2.2). Eventually, after substituting the nonanticommutating γ 5 prescription, i.e. (2.5) and (2.7), we choose the following reduced definite form of the 6 projectors to be used in the projections of (3.1) according to the D-dimensional algebra: where we have pulled out the normalisation factors as defined in (3.4) and the symbol γγγγ ≡ − i 24 µνρσ γ µ γ ν γ ρ γ σ is introduced. Upon pulling out the respective normalisation factors, all 6 projectors thus constructed have only a polynomial dependence on the kinematic variables and external momenta, which is very advantageous in practice. Computations of the projection of the amplitude (3.1) onto each of these 6 structures proceed in the usual way: the product of Dirac spinors will be cast into a trace of products of Dirac matrices with the aid of (unpolarised) Landau density matrices of external spinors. The contraction of pairs of Levi-Civita symbols appearing in the projection is made in accordance with (2.8) where all indices of the resulting metric tensors are set to be D-dimensional by definition. Furthermore, since by construction the momentum basis representations of polarisation state vectors in (3.3) fulfill all the defining physical constraints, the index contraction between these projectors (3.5) and the amplitude (3.1) is therefore always simply done with the space-time metric tensor g µν . Since the process (2.1) in question has only 3 linearly independent external momenta, upon using total momentum conservation, there will be no term that involves explicitly the Levi-Civita tensor after the projections have been made using (3.5). It is thus straightforward to see that the projections with P µ 1 , P µ 3 , and P µ 5 receive contributions only from the vector part of the Z boson coupling to quark, while contributions from the Z boson's axial coupling to quark survives only in the projections with P µ 2 , P µ 4 , and P µ 6 . In this way, the axial and vector part of the amplitude (3.1) get separated naturally.
Polarised amplitudes of this scattering process are thus first extracted in the linear polarisation basis using the 6 projectors as described above, from which helicity amplitudes can be readily composed afterwards. Helicity eigenstates of a massless spinor coincide with its chiral eigenstates, and according to the usual convention, a left chiral massless u(v)type spinor has a negative (positive) helicity. Therefore in the end we compose 6 polarised amplitudes by the following linear combinations: where the superscript, e.g. ++T , denotes the respective polarisations of the b quark, anti-b quark and the Z boson. Note that, each polarised amplitude can be decomposed into its vector and axial part, where the superscript [j] runs over all the six polarisation configurations shown in (3.6). The modulus squares of the normalisation factors involved are The complex phase factor of a helicity amplitude is unphysical by itself, and they can be freely set among those associated with independent external polarisation states. For definiteness, we choose the positive realvalued roots of the expressions in (3.8) when evaluating (3.6) in the physical regions.
The scattering amplitude in (3.1) can be written in terms of the polarised amplitudes in (3.6) (albeit, with a bit of abuse of notation) as follows: where [j] runs over all the six polarised states listed in (3.6) and φ j are the aforementioned corresponding unphysical phase factors. In computing any physical quantity where the squared modulus of the matrix element enters, this phase factor does not contribute owing to the orthogonality of the six polarised amplitudes: Regarding the polarisation states of the (massive) Z boson, its helicity states can be constructed optionally, as circular polarisation states from the linear ones, e.g.

Projectors for Conventional Form Factors
For the purpose of cross checking the results obtained with the projectors (3.5) as prescribed in the preceding section, we performed also the conventional form factor decomposition of (3.1). The following two subsections are thus devoted to a description of Lorentz structures included in the decomposition and the derivation of corresponding form factor projectors. Special attention is given to discussing a subtlety in the form factor decomposition of (3.1) with axial current vertices in D dimensions, which is similar to that known in the four-quark scattering amplitude qq → QQ discussed in [57,58].

The vector part
By Lorentz covariance the amplitude (3.1) can be expressed as a linear combination of a finite Lorentz structure basis at any finite order. These structures are constrained by physical requirements such as on-shell kinematics and symmetries of the dynamics. Here because of the attachment of a Yukawa interaction on the massless b quark line, the power of elementary Dirac matrices in Γ µ should be even to have a non-vanishing matrix element betweenv(p 2 ) and u(p 1 ) with opposite chirality. Under the condition of even powers in elementary Dirac matrices, combined with the P-even constraint from the vector coupling, we see that there are only 4 possible linearly independent Lorentz structures. Thus one can write down for the vector part of the amplitude (3.1): after taking into account also the equations of motion for the on-shell massless spinors v(p 2 ) and u(p 1 ). Notice that linear completeness of the Lorentz structure basis employed in (3.12) holds in D dimensions without demanding the transversality of the Z boson's physical polarisation states 5 , because we prefer to use the simple metric tensor g µν instead of the physical polarisation sum rule for index contraction in projections. This choice of Lorentz structures also leads to simple form factors clearly seen at the tree level. The linear decomposition coefficients, F vec,1 , · · · , F vec,4 , are Lorentz invariant functions of external kinematics, which are often called form factors. In addition these functions depend also on the value of dimensionality in dimensional regularisation, and their concrete expressions depend on the perturbative order at which they are computed. Once given the linear basis as well as the index contraction rule, one can then compute the Gram matrix of this linear basis, and its inverse gives us the projectors with which one obtains the Lorentz-invariant form factors. The projectors for the pure vector form factors defined in (3.12), derived under the aforementioned conditions, are still compact enough to be documented explicitly: With the tensor amplitude reconstructed according to (3.12) with form factors projected using (3.13), we can then compute the polarised amplitudes as defined in (3.6) and compare with the direct computation using (3.5), albeit only for contributions from the vector coupling of Z boson. In fact it is sufficient, and more convenient, to perform the comparison for the projections P µ iv (p 2 ) Γ µ u(p 1 ) with i = 1, 2, · · · , 6. To this end, one can first calculate matrix elements of each of the four Lorentz basis structures in (3.12) between the external state vectors as defined in (3.3). Namely one projects each of these Lorentz basis structures onto the list of projectors in (3.5), computed using the D dimensional algebra in exactly the same way as how the six projections P µ iv (p 2 ) Γ µ u(p 1 ) are computed. The resulting matrix provides us with the linear transformation needed to combine the vector form factors in (3.12) to get the projections P µ iv (p 2 ) Γ µ u(p 1 ) (keeping only vector coupling contributions). These matrix elements are purely rational functions in the Mandelstam variables of external kinematics (and also the dimensionality D in general). Multiplying these transformation matrix elements with the respective vector form factors accordingly gives us the projections P µ iv (p 2 ) Γ µ u(p 1 ), and subsequently polarised amplitudes in (3.6). We will come back to this in section 6. In view of this discussion, it is also worthy of mentioning that to alleviate some of the related difficulties within the conventional form factor decomposition approach, it is advocated in ref. [59] to combine the step of projecting out form factors with the step of composing helicity amplitudes in one go.
Because the four basis structures in (3.12) are linearly complete in D dimensions for M vec in (3.1) (regardless of the QCD loop order), the tensor amplitude reconstructed according to (3.12) with form factors projected using (3.13) (retaining their full D dependencies) is thus a faithful representation of M vec in D dimensions. Owing to this, the aforementioned comparison can be made directly for the un-renormalised un-subtracted bare virtual amplitudes, and the agreement should be exact with full D dependence. This is indeed what we saw in our calculation as will be discussed further in section 6, which serves as a strong check of the computational setup we have established so far.
In contrast, this kind of agreement is no longer observed regarding the axial part of the amplitude (3.1), i.e.,v(p 2 ) Γ µ axi u(p 1 ), to be discussed in detail in the following subsection. This is because in this case it is a non-trivial task to exhibit, beyond the tree level amplitude, the full linearly-complete basis for the (conventional) form factor decomposition in D dimensions (with a non-anticommuting γ 5 prescription). However, we will show through this work that as long as one is only concerned with physical quantities (or quantities that actually contribute to physical observables) it is not necessary to find out and use such an "ultimate" basis for form factor decomposition in a calculation done with dimensional regularisation.

The axial part
Since we are employing a constructive prescription of γ 5 in DR, i.e. (2.5), where the full anti-commutativity is sacrificed (in exchange of the cyclicity of trace), algebraic equivalences between expressions in 4 dimensions under the condition of a fully anti-commuting γ 5 may no longer hold true for their D-dimensional counterparts. Indeed the manual symmetrisation needed in the proper definition of an Hermitian axial current in DR, as reflected in (2.6), is already a prominent example of such a point. However, full anti-commutativity of γ 5 is routinely exploited when establishing a linearly independent (and complete) basis of Lorentz structures for amplitudes involving axial vertices in 4 dimensions (in addition to the basic Dirac algebra and also on-shell equations of motion). This implies that a form factor decomposition basis that is claimed to be linearly independent and complete in 4 dimensions for an amplitude involving axial vertices may not necessarily remain truly linearly complete for this amplitude regularised in D dimensions. As a matter of fact, such an interesting point can be verified straightforwardly for M axi in (3.1) already at the tree level.
Since this subtle point has not yet been discussed much in the literature, let us be more specific about it here. First of all, following the reasoning similar to what was done for the vector part in the previous subsection, we could write down the following four linearly independent Lorentz structures for the axial part of the amplitude ( which is linearly complete in 4 dimensions under the condition of even powers in elementary Dirac matrices combined with the P-odd constraint from the axial coupling. One then substitutes the 't Hooft-Veltman definitions, i.e., (2.5) and (2.7), for the bookkeeping γ 5 in (3.15). This yields a list of structures to be used in decomposing M axi . Once given the linear basis (3.15), one can then compute the Gram matrix of this linear basis whose inverse gives us the projectors with which one obtains the corresponding Lorentz-invariant form factors, similar to the vector part. In particular, pairs of Levi-Civita tensors will be contracted in accordance with (2.8) with the resulting space-time metric tensors assumed to be D dimensional. The projectors for the form factors corresponding to axial Lorentz structures with their bookkeeping form given in (3.15), derived under the aforementioned conditions, read as: The Gram matrix associated with (3.15) used in deriving the above projectors has a matrix rank 4, confirming the fact that these four structures (3.15) are linearly independent (in D dimensions). Now, to address the question of whether M axi lives in a space linearly spanned by (3.15) in D dimensions, we enlarge this list by appending the tree level expression M axi directly given by its Feynman-diagrammatic representation with the axial current regularised as described in section 2. We then compute the Gram matrix of this enlarged list of five elements in exactly the same way as before, and we find out that the matrix rank is increased to 5 rather than staying at 4. This linear dependence test clearly implies that the tree-level axial amplitude M axi (in D dimensions) is not linearly dependent on the four structures (3.15) according to the usual D-dimensional computational prescription. However, the similar kind of linear dependence test can be performed for the vector part as well, and there we do see that the matrix rank of the accordingly enlarged list is not increased but stays at 4, confirming the fact that the M vec does live in the space linearly spanned by (3.12) in a D-dimensional computation. At least for M axi at the tree level, it is not too difficult to find out a linearly complete Lorentz structure basis in D dimensions for it (e.g. simply enlarging the basis (3.15) by collecting additional tensor structures until enough 6 ). This leads to a larger linear decomposition basis, which in turn leads to projectors more complicated than (3.16). However, by the argument made above, a priori, a D-dimensional linearly-complete basis at the tree level is again not guaranteed to be linearly complete and faithful for M axi at loop orders. Such an issue is very similar to what was already observed in the D-dimensional form factor decomposition of four-quark scattering amplitude qq → QQ in [57,58], albeit of a different technical origin.
There is also another manifestation of this important point in the case at hand, which one can easily examine. One can first apply the thus-derived projectors (3.16) onto M axi at the tree level to get the corresponding axial form factors. Then the axial part of the amplitude (3.1) can be reconstructed similarly to what was done for the vector part in accordance with the formula (3.12). However, one can verify explicitly that with the amplitude reconstructed in this way one can not recover the unpolarised squared modulus of the axial part of amplitude (3.1) (at the tree level) computed using the physical polarisation sum rule in CDR directly from its original Feynman diagrammatic representation. This tells us again that the computation of the tensor amplitude with (3.15) is not algebraically identical to the original M axi (even at the tree level) in D dimensions. Because of this we can not write down an exact decomposition identity for M axi just in terms of the basis (3.15) in analogy to its vector counterpart (3.12).
Moreover, the projectors in (3.16) are not very convenient to use in practice, especially if one considers some more complicated versions with additional D-dimensional linearly independent Lorentz structures incorporated in the form factor decomposition basis. However, a priori, it is not clear whether omitting additional (evanescent) Lorentz structures in the form factor decomposition of a D-dimensional amplitude, i.e., using a D-dimensional linearly-incomplete decomposition basis like (3.15) forv(p 2 ) Γ µ axi u(p 1 ) at loop orders, would still always lead to correct results, because for sure the so reconstructed amplitude is not algebraically identical to the original defining form (given by Feynman diagrams). Bold moves of this kind have already been tried out, e.g. in the computation of qq → tt amplitudes to two-loop order in ref. [60]. On top of this undecided situation, there is also no clear statement in the literature about whether one can always set the dimension variable D = 4, i.e. in the curly brackets of the expressions in (3.16), and still expect with confidence that the amplitude reconstructed in this way yields correct results.
Let us now summarize the issues we have just discussed regarding the form factor decomposition of an amplitude involving axial vertices in D dimensions, in particular the expressionv(p 2 ) Γ µ axi u(p 1 ) at hand.
1. It is not easy to construct the full D-dimensional linearly complete basis for the form factor decomposition ofv(p 2 ) Γ µ axi u(p 1 ), especially at increasing loop orders.
Concerns like these were among the motives that drove us to take the approach as prescribed in ref. [42]. However, it is still very interesting to investigate whether the aforementioned bold moves, especially combining the second and third point above, would work all the way through. This is then what we did for the axial part of amplitude (3.1), in addition to the computation using the projectors discussed in section 3.1, in order to answer this question. Namely, not only did we choose a set of form factor decomposition structures that are known to be linearly incomplete in D dimensions (albeit, linearly independent and complete in 4 dimensions), but we also set manually all explicit D in the resulting form factor projectors to be 4. In other words, for the form factor decomposition of M axi we used (3.16) with D=4 in the curly brackets, i.e., a set of form factor projectors that are essentially identical to the respective expressions that would be used in a four-dimensional form factor decomposition 7 . Namely we define our "axial form factors" (throughout this article) as where the [4] in the superscript denotes the fact that we set D = 4 in the original P µ i,axi . Note that, however, when applying the so determined "axial form factor projectors" to M axi , any pair of Levi-Civita symbols will always be contracted in accordance with (2.8) where all indices of the resulting metric tensors are D-dimensional [50]. Subsequently we build up an intermediate axial amplitudeM µ axi defined as M µ axi ≡ F 1,axiv (p 2 ) γ 5 u(p 1 ) p µ 1 + F 2,axiv (p 2 ) γ 5 u(p 1 ) p µ 2 + F 3,axiv (p 2 ) γ 5 u(p 1 ) q µ 1 + F 4,axiv (p 2 ) γ µ γ 5 / q 1 u(p 1 ) .

(3.19)
This quantityM µ axi is not algebraically identical (or faithful) to the original Feynman amplitude M µ axi =v(p 2 ) Γ µ axi u(p 1 ) in D dimensions, due to issues discussed above, while we expect that the two should eventually lead to the same properly defined finite remainder in four dimensions.
Indeed, with well established results for polarised amplitudes obtained using physical projectors defined in section 3.1, not relying on any explicit Lorentz tensor decomposition of the original Feynman amplitude (both conceptually and technically), we have verified that eventually the same finite remainders for M axi were obtained with both approaches, and hence we confirm a positive answer to the question raised above for the process (2.1). We will come back to this later in section 6. The comparison regarding these axial amplitudes can proceed as described at the end of the last subsection, but with one crucial difference compared to their vector counterparts: the agreement can only be expected at the level of 7 Let us stress that in this acrobatic axial form factor decomposition there is no imposition of any explicit dimensional splitting, which would have allowed one to keep strictly only 4 dimensional d.o.f. in the external projectors (clearly separated from all orthogonal evanescent ones) such as typically done in the HV scheme.
properly defined finite remainders in four dimensions, due to the subtle points explained above.

UV renormalisation
The bare scattering amplitude of the process (2.1) beyond the leading order (LO) is not finite, and it contains poles in dimensional regulator (≡ 1/2(4 − D)) arising from ultraviolet, soft and collinear regions of the loop momenta. The UV renormalisation of the amplitude (3.1) is done in the modified minimal subtraction (MS) scheme. To be more specific, the UV divergences in the vector part of the amplitude (3.1) in massless QCD are handled by the QCD coupling constant renormalisation througĥ and the renormalisation of the Yukawa coupling througĥ where S ≡ exp [ (ln 4π − γ E )] with Euler constant γ E = 0.5772.... In (4.1), (4.2) and throughout this article, hat (ˆ) represents the bare quantity, a s ≡ α s /4π is the strong coupling constant, µ is an auxiliary mass-dimensionful parameter introduced in dimensional regularisation to keep the coupling constants dimensionless and µ R is the usual renormalisation scale. These two QCD renormalisation constants are given to two-loop order in the MS scheme by [61]: The quantities C A = N c and C F = (N 2 c − 1)/(2N c ) are the eigenvalues of the quadratic Casimir operators in the adjoint and the fundamental representations of the SU(N c ) gauge group, respectively. However, (4.1) and (4.2) are not sufficient to properly UV renormalise the axial part of (3.1) with γ 5 regularised in DR as described in section 2. The particular prescription for the γ 5 -related quantities, such as (2.7) for the axial current which violates the anti-commutativity of γ 5 with γ µ , gives rise to some additional spurious UV poles in in dimensional regularisation which has to be removed manually [46]. Manifestation of this issue shows up as the loss of correct chiral Ward identities for axial currents in the original prescriptions by 't Hooft-Veltman [43] and Breitenlohner-Maison [45] (at higher loop orders). This amendment is typically realized in form of some additional UV renormalisation constants [46][47][48]62] introduced specifically for the axial part of the amplitude (on top of the aforementioned renormalisations common to the vector part). The concrete expressions of these additional finite renormalisation constants depend on the treatment of the Levi-Civita tensor appearing in the constructive expression eq.(2.5) for the nonanticommutating γ 5 . In practice, it is convenient to have this rectification done with the so-called flavour singlet and non-singlet Feynman diagrams separated, which we now turn to in the following subsections.

UV Renormalisation of the Flavour Non-Singlet Current
Samples of the flavour non-singlet, or non-anomalous, set of Feynman diagrams are shown in figure 2 at one-loop order and in figure 3 at two-loop order, where there is no appearance of the closed triangle fermion loop with axial Z coupling vertex (related to the quantum anomaly to be addressed in the next subsection). For these diagrams, the axial part of the coupling between the Z boson and the b quark takes place through the flavour non-singlet axial current where I 3 denotes the third component of the (weak) isospin operator. The appearance where, as already mentioned, the hat (ˆ) represents the bare operator. The aforementioned renormalisation constants up to two-loop order are given by [47]: The properly UV renormalised non-singlet axial current in (4.6) is non-anomalous and should be conserved for massless quarks, i.e. the standard Ward identity holds and the renormalisation invariance of the current is restored (just like its vector counterpart). While, due to the presence of the non-zero Yukawa coupling λ b in the amplitude (3.1) under consideration, the currentv(p 2 ) Γ µ axi u(p 1 ) (which is a Green correlation function involving the axial current operator J ns µ,A ) is no longer conserved even though the b quark mass is set to be zero. This can be verified explicitly by replacing the polarisation vector of the Z boson by its momentum in the amplitude M axi which yields q 1,µv (p 2 ) Γ µ axi u(p 1 ) = 0, and this expression is proportional to λ b . This non-conservation is "classically expected" and hence not regarded as a "quantum anomaly", which we will turn to in the next subsection. While on the other hand q 1,µv (p 2 ) Γ µ vec u(p 1 ) = 0 is still true.

UV Renormalisation of the Flavour Singlet Current and ABJ Anomaly
The flavour-singlet (or anomalous) Feynman diagrams, characterised by featuring a triangle fermion loop with axial Z coupling, start to appear at two-loop order for the process in question. They are shown in To UV renormalise the singlet axial current properly, as in the previous case, in addition to performing an overall MS renormalisation, we also need to multiply with an additional finite renormalisation constant [48,51,62]: The renormalisation constants up to two-loop order read as [48]: Unlike the flavour non-singlet axial current, the singlet axial current does not satisfy the standard Ward identity even in the massless quark limit. This current exhibits anomalous properties which are known as axial or Adler-Bell-Jackiw (ABJ) anomaly [63,64]. The operator relation for the ABJ anomaly of massless axial current reads [65]: where GG ≡ µνρσ G a µν G a ρσ and G a µν is the gluonic field strength tensor. The subscript R of the composite local operators on both sides of (4.11) represents that these composite local operators need to be properly renormalised so that this operator relation holds. The operator on the left-hand side is renormalised multiplicatively in the same way as the current (4.9) itself, while the renormalisation of the GG can be found in [48,66] in detail. The concrete expression of Z s 5,A in (4.10) is determined precisely such that the properly UV renormalised singlet axial current in (4.9) has an anomaly that does obey (4.11), i.e. , that preserves the one-loop character of this relation. Technically speaking, (4.11) is irrelevant regarding the contributions to M axi from this flavour singlet set of 6 Feynman diagrams all depicted in figure 4. This is because the polarisation states of the external on-shell Z boson in the physical amplitude from these diagrams, either transversal or longitudinal as listed in (3.3), are always orthogonal to its own momentum. However, we verified this operator level relation explicitly within our computational setup for the amplitude M axi in order to perform an indirect check of our treatment of these anomalous diagrams. To be a bit more specific about this, let us denote the contribution to M axi from these 6 anomalous Feynman diagrams byv(p 2 ) Γ µ ABJ u(p 1 ) ε * µ . The Green correlation function of the external fields with the lefthand side of (4.11) (or rather, the matrix element of the external operator ∂ µ J s µ,A R between the vacuum and the external on-shell states of our process (2.1) excluding the Z boson) reads asv(p 2 ) Γ µ ABJ u(p 1 ) q 1,µ . Despite the non-vanishing anomaly ∂ µ J s µ,A R at the operator level,v(p 2 ) Γ µ ABJ u(p 1 ) q 1,µ would have completely vanished (since b quark is taken massless) if it were not due to the presence of the b quark Yukawa coupling in our calculation. Replacing ∂ µ J s µ,A R by the operator on the right-hand side of (4.11) leads to 3 one-loop diagrams, as depicted in figure 5, which are connected to those in figure 4 by shrinking the fermion triangle loop into an effective vertex between two gluons and one imaginary or auxiliary pseudo-scalar. These one-loop diagrams would have also been vanishing if there would be no Higgs radiated from the (massless) b quark, simply due to the tension between the angular momentum conservation and chirality conservation along the massless b quark line (without Yukawa coupling). In the end we computed these one-loop diagrams andv(p 2 ) Γ µ ABJ u(p 1 ) q 1,µ , both of which are proportional to the Yukawa coupling λ b , and found an exact agreement between the two and hence verified (4.11) in our computational setup.
In view of the classification of contributing diagrams as discussed above, the scattering amplitude (3.1) as well as the polarised ones defined in (3.6) are thus conveniently separated into pieces with their respective renormalisations as follows: where the superscript [j] runs over all the six polarisation configurations as in (3.6). Upon applying the projectors (defined in section 3) onto the Feynman amplitudes, we compute the bare polarised amplitudesM

Computation of the Amplitudes
Despite the presence of modern techniques based on unitarity (or on-shell cuts), the Feynman diagrammatic approach of computing loop amplitudes still remains a popular choice which is employed in this article, especially in view of many available efficient tools for Feynman diagram generations and manipulations for a large variety of theory models. 9 .
The technical aspects of the computation of bare Feynman amplitudes in (4.12) closely follow the steps used in the recent calculation of doubly massive form factors [67] in maximally supersymmetric Yang-Mills theory. Feynman diagrams are generated symbolically using QGRAF [68]. For the non-Drell-Yan type contributions to the process (2.1), there are 2 and 10 diagrams at the tree and one-loop level, respectively. At two-loops, there are 153 flavour non-singlet and 6 flavour singlet diagrams generated in our setup. The symbolically generated diagrams are passed through a series of in-house codes based on FORM [69] in order to apply the Feynman rules, perform SU(N c ) colour and D-dimensional Lorentz and Dirac algebras. Upon multiplying the properly constructed projectors, every projected Feynman amplitude is expressed as a linear combination of a large number of scalar Feynman integrals belonging to the family of four-point amplitudes with two offshell legs of different virtualities. Using the liberty of transforming the loop momenta, all the scalar Feynman integrals are categorised into three and six different integral families at one-and two-loop order, respectively, with the help of REDUZE2 [70,71]. These scalar integrals can be reduced to to linear combinations of a much smaller number of loop integrals, called master integrals (MI), using integration-by-parts (IBP) [72,73]. To perform the IBP reduction, we use LiteRed [74] along with Mint [75] at one-loop and Kira [76,77] at two-loop level. Upon performing the IBP reduction, we get at the two-loop level 134 MI which are further reduced to a set of 84 independent integrals by taking into account additional relations and crossings of the external momenta. In a parallel setup, we first obtain the unreduced symbolic form of polarised amplitudes using an extension of the program GoSam [78][79][80] at one-loop and two-loop order. The list of unreduced loop integrals appearing is then extracted and fed to Kira [76,77] to obtain a table of generally usable IBP rules. Insertion of the IBP table and subsequent simplification of rational coefficients in front of MI are performed with an in-house routine based on a parallelised usage of Mathematica and fermat [81]. We have checked that the final reduced forms of all projected amplitudes in terms of MI are identical for these two parallel setups.
The set of independent MI involved in our amplitudes agrees with those in the fourpoint amplitudes with two massive legs with different virtualities which were first computed in refs. [82][83][84]. A subset of these MI were computed in refs. [85,86]. Later, an independent computation of these MI was performed in ref. [87] where the solutions are optimised for efficient numerical evaluation. For our present calculation, we use the optimised solutions of the MI computed in ref. [87] which are available in HepForge [88] in computer readable format. In order to use the optimised solutions, we introduce the same set of dimensionless quantities defined in ref. [87] constructed out of the Mandelstam variables (2.2) This choice of variables also rationalises the root of κ. The physical region is bounded by the constraints The symbol alphabet of the MI involves 19 letters (see ref. [87] for a detailed description). The results of the MI are computed as Laurent series in up to transcendental weight 4 which enables us to get the one-loop polarised amplitudes, M [j], (1) , to O( 2 ) and the two-loop ones, M [j], (2) in (4.13), to O( 0 ) in dimensional regularisation. While the analytic results of the UV renormalised polarised amplitudes defined in (3.6) and (4.12) are too lengthy to be presented explicitly here, we thus attach these results as ancillary files in Mathematica format along with the arXiv submission 10 . Details about conventions and variables of the saved analytic expressions can be found in the ReadMe.txt submitted along with these files. In the following section, we will discuss the universal infrared divergences present in these UV renormalised amplitudes and the (fourdimensional) finite remainders subsequently defined.

IR factorisation and RS Independent Finite Remainders
In a gauge theory like QCD, the UV-finite amplitudes beyond leading order typically contain divergences arising from the soft and collinear configurations of the loop momenta, which appear as poles in dimensional regulator . Fortunately these IR divergences systematically factor out from the amplitudes to all orders in perturbation theory [89,90], which demonstrates a kind of simple universal structure in these additional divergences. The explicit form of these factorised universal IR poles in massless QCD, being only dependent on the nature of the external coloured particles, was first determined up to two-loop order in terms of universal IR subtraction operators in ref. [91]. In ref. [92], a detailed derivation was presented by exploiting the factorisation and resummation properties of scattering amplitudes which was later generalised to all orders in terms of a soft anomalous dimension matrix in refs. [93,94].

The RS Independent Finite Remainders of M [j]
As briefly mentioned at the beginning of section 3.1 and also argued in detail in ref. [42], the bare and UV renormalised polarised amplitudes resulting from the projectors constructed in section 3.1 are in general different from their counterparts defined in dimensional regularisation schemes such as CDR, HV, DRED and FDH (and hence they should not be compared blindly without conversion). However, crucially, the properly defined finite remainders in four dimensions are guaranteed to be the same, known as the regularisation-scheme (RS) independence of these objects. In other words, the projection prescription adopted here implies a specific regularisation prescription that, despite being different from its companions, remains also unitary in the sense as defined in refs. [95,96]. As long as one is only concerned with quantities that actually contribute to physical observables, identical results should always be obtained.
To appreciate this crucial point, recall that in this particular prescription, all open Lorentz indices of the (physical) polarisation projectors are set to be D-dimensional and no dimensional splitting is ever introduced, just like in CDR, which consequently ensures the commutation between Lorentz index contraction in the projection and loop integration. This means that applying the projectors defined in (3.5) directly to the original Feynman-diagrammatic representation of a loop amplitude should lead to the same polarised amplitudes as would be obtained by applying these projectors to a conceivable (faithful) D-dimensional form factor decomposition representation of that amplitude. The difference between this particular regularisation prescription and the CDR, with γ 5 regularised and treated in accordance with refs. [47,48], concerns merely the external particles' state vectors, which are decoupled from the UV and/or IR singularities contained in the loop integration, as is evident from the multiplicative UV renormalisation and IR factorisation as discussed in the preceding sections. Regarding the latter part, captured by certain "universal" factors, both regularisation prescriptions are exactly the same. Thus the regularisation convention implied by projectors defined in (3.5) shares the identical set of renormalisation constants (as given in section 4) and IR anomalous dimensions as the CDR (with γ 5 regularised and treated technically in accordance with refs. [47,48]). From this point of view one can readily expect to end up with the same (four-dimensional) finite remainder in this regularisation prescription as one would obtain from a computation purely within the CDR (and also the same as in any other variant of the dimensional regularisation).
The IR pole structures in the UV renormalised polarised amplitudes (4.13) can be exhibited through a parameterisation in terms of the universal IR subtraction operators I (i) ( ), depicted in [91] as (6.1) The explicit expressions of the IR subtraction operators for the current process are given in the CDR scheme by with the cusp anomalous dimension K and the quantity H (2) ( ) [93]: Based on the discussion given above, we thus expect that the explicit expressions of the IR subtraction operators I (i) ( ) in CDR given by (6.2) can be directly used to predict the IR poles contained in our UV renormalised polarised amplitudes according to Catani's IR factorisation formula (6.1). Indeed this is precisely what we have observed in all our polarised UV renormalised amplitudes, defined by (4.12) and power expanded in (4.13). This thus serves as a strong check on the correctness of our computation. Upon subtracting all these predicted IR singular pieces from the UV-finite polarised amplitudes, we obtain the finite remainders M . The analytic expressions of these finite remainders are too lengthy to be presented explicitly here, but they can be extracted from our UV renormalised polarised amplitudes, attached as ancillary file in Mathematica format, according to (6.1).

Same Finite Remainders Recovered from Form Factor Decomposition
As alluded to in section 3.2, we like to cross-check these finite remainders, M vec , where the comparison is really straightforward. The vector amplitude reconstructed according to (3.12) with form factors F i,vec projected using (3.13) is truly a faithful representation of M vec in D dimensions (defined by its original Feynman diagrammatic representation), because the basis structures in (3.12) are linearly complete w.r.tv(p 2 ) Γ µ vec u(p 1 ) in D dimensions regardless of the QCD loop order. Owing to this, the comparison can actually be made directly for the unrenormalised un-subtracted bare virtual amplitudes, and the agreement should be exact with full D dependence (i.e. no expansion and truncation in ). Furthermore, since the linear composition with the projections P µ iv (p 2 ) Γ µ u(p 1 ) to the polarised amplitudes M [j] , given in (3.6), is fixed independent of howv(p 2 ) Γ µ u(p 1 ) is computed, it is thus sufficient, and technically more convenient, to compare directly at the level of these projections. As already sketched at the end of section 3.2.1, the transformation connecting form factors F i,vec to P µ iv (p 2 ) Γ vec,µ u(p 1 ) reads as where i runs from 1 to 6 as in (3.5) and the contraction involved in the big parentheses should be done according to D dimensional Lorentz/Dirac algebra. We have checked explicitly that there is an exact agreement for these six bare quantities (with full D dependence) between the direct projection calculation and the one with a detour of conventional D-dimensional form factor decomposition. From this, it follows straightforwardly that the 4-dimensional finite remainders M vec,fin will be the same because the UV renormalisations and IR subtractions proceed identically in both computations.
Concerning the vector part of the amplitude (3.1) done with a strictly D-dimensional faithful form factor decomposition, there is really no surprise that the same finite remainders are obtained. Regarding the axial part, the projectors for projecting our "axial form factors" defined in (3.18) are a bit un-conventional due to two points as discussed in section 3.2.2: (1) the basis set (3.15) is known to be linearly incomplete in D dimensions for M axi in question; (2) all explicit D appearing in the corresponding projectors (3.16) are set manually to be 4 even thought it is not an overall D dependence. This makes our form factor decomposition for the axial part not qualified as being called D-dimensional faithful. Still we defined an intermediate axial amplitudeM axi in (3.19) from the so determined "axial form factors". As mentioned before it would be very interesting to see whether the correct 4-dimensional finite remainders M [j] axi,fin could still be obtained in a computation based on such an acrobatic version of axial form factor decomposition.
The comparison can proceed following a similar line as the vector part just discussed above, with one important exception: there is now no more any reason for expecting that the amplitudes before UV renormalisation and IR subtraction agree, and indeed they do not in our computations. Consequently, we first properly renormalise our axial form factors F i,axi in (3.18) according to section 4, and then apply the previously discussed IR subtractions onto the resulting UV-finite axial form factors to end up with their finite remainders in four dimensions, denoted as F i,axi,fin respectively. Afterwards these fourdimensional regular objects are combined together according to F 1,axi,fin P µ iv (p 2 ) γ 5 u(p 1 ) p 1,µ + F 2,axi,fin P µ iv (p 2 ) γ 5 u(p 1 ) p 2,µ + F 3,axi,fin P µ iv (p 2 ) γ 5 u(p 1 ) q 1,µ + F 4,axi,fin P µ iv (p 2 ) γ µ γ 5 / q 1 u(p 1 ) , (6.5) where i runs again from 1 to 6 as in (3.5) while the contractions involved in the big square bracket can be done according to four-dimensional Lorentz/Dirac algebra. We have verified explicitly that these four-dimensional finite remainders composed in (6.5) agree exactly with those from the renormalised and subtracted P µ iv (p 2 ) Γ axi,µ u(p 1 ) computed from direct projections as described in section 3.1. Therefore in this work we confirm that, regarding the process (2.1), it is not necessary to construct and use the "ultimate" D-dimensional axial decomposition basis, as long as one is only concerned with physical quantities (or quantities that actually contribute to physical observables in four dimensions). An acrobatic version of axial form factor decomposition such as described in section 3.2.2 is sufficient even in a calculation done with dimensional regularisation. Of course, a similar statement applies also to the vector part of the amplitude (3.1), which we also checked explicitly in our computational setup.

Cross Section to NNLO in the Soft-Virtual Approximation
As we discussed in the introduction 1, the dominant mode for the ZH production at LO is the quark-antiquark annihilation which has three channels viz. s, t and u, as shown in figure 1. The t, u channels result only from b quark annihilation due to the presence of b quark Yukawa coupling. In the SM, the b quark Yukawa coupling is small owing to its proportionality to the b quark mass. Moreover, the distributions of b quark inside the proton is much smaller than those of other light quarks. Consequently, the contributions arising from the t-and u-channels are expected to be much smaller than that of the s-channel. In previous analyses as reviewed in the introduction, contributions from t, u-channels are not taken into account due to their expected small size. Nevertheless, it is important to have a definite estimate of the size of these t, u-channel contributions for the current precision studies. In this section, we study the numerical impact of the t, u-channel contributions. Since these channels do not interfere with the s-channel, they can be treated separately. We obtain their contributions up to NLO in QCD using Madgraph [97] and to NNLO in the soft-virtual approximation [28,39,41], which is also known as the threshold limit. The latter requires the knowledge of the two-loop amplitudes, as well as the tree-level and one-loop amplitudes to high powers in , which are presented for the first time in this article.
The hadronic cross section of the ZH production that results from the t, u-channels of b quark initiated partonic sub-processes is given by Here, x 1,2 are fractions of momenta of incoming hadrons carried away by bottom quarks b andb, and f a is the parton distribution function (PDF) normalised at the factorisation scale µ F . The mass factorised partonic cross section is given by σ ZH bb,tu . Both the LO and NLO contributions to σ ZH bb,tu are straightforward to obtain. We use Madgraph [97] for this. For the NNLO cross section, in addition to the virtual contributions computed in this article, the double-real and real-virtual contributions are also needed. The inclusion of the exact real radiations is beyond the scope of the present work. While, in the absence of this, owing to the universality of soft and collinear contributions, an approximated result of these real radiation corrections can be computed following [28,39,41], which combined with the virtual part gives a finite quantity, called the soft-virtual (SV) cross section. At the partonic level, the dominant contribution to the SV cross section comes from the terms proportional to distributions of the kind δ(1 − z) and D j (z) which is defined as The variable z is defined as Q 2 /s where Q 2 = (q 1 + q 2 ) 2 is the invariant mass square of the final state ZH. For the current process, such contributions can arise only from the b-quark initiated sub-processes. Hence, we write the hadronic cross section resulting from t, u channels as Thanks to the factorisation of the matrix elements as well as the phase space into soft and hard parts in the soft limit, the SV part (σ ZH,SV bb,tu ) of the inclusive cross section is found to be σ ZH,SV bb,tu where ∆ ZH,SV bb is the mass factorised partonic cross section in the threshold limit which is calculable order by order in strong coupling constant a s . Suppressing the obvious arguments, we expand ∆ ZH,SV bb as Using the one and two loop amplitudes computed in this article, the universal soft distribution functions and the Altarelli-Parisi kernels, ∆ SV,(j) bb are obtained in terms of cusp A q , soft f q and collinear B q anomalous dimensions. Setting the renormalisation and factorisation scales at Q, i.e., µ R = µ F = Q, we find Here denotes the real part of the quantity within the respective parenthesis. The quantity M (j) k is the coefficient of a j s k in the expansion of the matrix element M in (3.9) through: Using the relation between the unpolarised squared matrix element squares and the squared modulus of the polarised amplitudes given in (3.10), we relate these quantities order by order in a s which read where the polarised amplitudes are expanded in powers of a s in (4.13). Using the results of the polarised amplitudes which are computed in this article to O( 4 ), O( 2 ) and O( 0 ) at tree, one-and two-loop level, respectively, we have obtained the square of the matrix elements up to two-loop level. Moreover, by expanding the square of the matrix elements that appeared on the left-hand side of the aforementioned equation in powers of following (7.7) and comparing these with the explicit results obtained in this article, we have extracted the relevant quantities appearing in the results of the SV cross section in (7.6).
The cusp [98][99][100][101][102][103][104], soft, and collinear anomalous dimensions are given by The universal constants G q,(j) k arising from the soft-collinear distribution are given by [28] G q,(1) 1 where, ζ 2 = 1.64493407 · · · , ζ 3 = 1.20205690 · · · . Now using the matrix elements along with other universal quantities in (7.6), we obtain the SV cross section up to NNLO which is finite in the limit → 0. The cancellation of soft, collinear singularities among virtual, real emissions and mass factorisation kernels in the SV cross section serves as a check of the correctness of our renormalised polarised amplitudes, despite that they are given for a specific regularisation prescription implied by projectors devised in section 3.1. The size of the partonic SV cross section at NLO and NNLO are 16 KB and 29 MB respectively, hence we skip presenting them explicitly. Instead the results are included in Mathematica format as ancillary files along with the arXiv submission.
In the rest of this section, we discuss the numerical result of the total cross section of ZH production at the LHC. The center-of-mass energy is taken to be √ S = 13 TeV. We use MMHT2014 PDF sets obtained through LHAPDF interface [110], and use the five-flavour scheme throughout. The renormalisation and factorisation scales are set at µ R = µ F = m h + m z . The Fermi constant is G F = 1.16637 × 10 −5 GeV −2 . The width of Z is given by Γ Z = 2.4952. We use running b quark mass in the Yukawa coupling. We choose m b (m b ) = 4.18 and evolve it to appropriate scales using the renormalisation group equation. In table 1, we present QCD (and electroweak) corrections to the ZH production at the LHC@13TeV at various orders, with contributions from t, u-channels separated from that of the s-channel. The s-channel contributions include only those processes where the Higgs boson couples to Z, whereas the (t + u)-channel contributions contain those where the Higgs boson is radiated from the b quark. The s-channel contributions are obtained using the program vh@nnlo [105] where the electroweak (EW) corrections [106,107] are also implemented. The NLO QCD correction increases the total inclusive cross section of the ---- Table 1. QCD and electroweak (EW) contributions to ZH production at the LHC with center of mass energy 13 TeV. The σ ZH gg refers to the contribution coming from the gluon initiated subprocesses. The top quark loop contribution is denoted by σ ZH qq (top). All the numbers of the cross-sections are in unit pb.
ZH production at the LHC@13TeV by 31%, and the NNLO Drell-Yan-like correction alone by another 4%, and the threshold contribution at N 3 LO gives less than 1%. Apart from the Drell-Yan like contributions, starting from NNLO in QCD there appear also sub-processes induced by the top quark loop from which the Higgs is radiated. In addition, there is also the gluon initiated sub-processes at NNLO which give rise to about 6% correction for ZH production. (See table 2 in ref. [19] for the sizes of various contributions to ZH production at the LHC at different center-of-mass energies.) Regarding the contributions from the t, u-channels, we obtain the result up to NLO in QCD using Madgraph [97], and compute the hadronic cross section (7.4) at the NNLO using an in-house fortran routine that employs GiNaC [108] for the numerical evaluation of polylogarithms [109]. From the numbers in the last column of table 1, we see that these t, u-channels give a three orders of magnitude smaller contribution compared to that of the s-channel one.

Conclusions
A precise knowledge of the ZH process is important for fully exploring Higgs physics at hadron colliders, especially to meet the foreseeable precision requirements from future collider experiments for detailed studies of the 125 GeV Higgs boson (as well as potentially non-standard Higgs bosons). Through the work presented in this article we provide the analytic results of the two-loop massless QCD corrections to the b-quark-induced ZH process that involves a non-vanishing b-quark Yukawa coupling λ b , which is a necessary ingredient of the complete O(α 2 s ) QCD corrections to this process in the five-flavour scheme. The computation was performed by projecting the D-dimensional scattering amplitudes (in their Feynman-diagrammatic representations) directly onto an appropriate set of Lorentz structures related to the linear polarisation states of the Z boson, as explained in detail in section 3.1. We emphasize that the renormalised polarised amplitudes, attached as ancillary files, are given for a specific regularisation prescription implied by these special projectors in use, and hence they should not be compared blindly with their counterparts defined in various other DR schemes (e.g. CDR, HV, DRED and FDH, etc). However, the crucial point is that this special dimensional regularisation prescription is also unitary as argued in ref. [42], and more importantly it shares the identical set of renormalisation constants and IR anomalous dimensions as the CDR (with γ 5 -related objects treated tech-nically in accordance with refs. [47,48]). Because of this, one will always end up with the same properly defined RS-independent finite remainders, which is eventually required for computing physical observables. We have verified explicitly that the IR poles contained in our renormalised amplitudes match with those predicted by Catani's IR factorisation formula in section 6 as well as those according to the soft-virtual factorisation formula in section 7. This serves as a strong check of the correctness of our results.
In addition, regarding the vector part of the amplitude (3.1) considered in this publication, we have cross-checked the finite remainders of these renormalised polarised amplitudes, defined according to (6.1), with those obtained from a strictly D-dimensional conventional form factor decomposition approach and found exact agreement. Concerning the conventional form factor decomposition of amplitudes involving axial currents regularised in D dimensions (with a non-anticommuting γ 5 prescription like (2.5)), we have dedicated a few specialised discussions to an interesting issue resembling a similar one that appears in qq → QQ [57,58], albeit of a different technical origin. Despite all concerns discussed in section 3.2.2, through computations accomplished in this article we confirm that regarding the process (2.1) it is not necessary to construct and use the "ultimate" D-dimensional axial decomposition basis. An acrobatic version of axial form factor decomposition such as described in section 3.2.2 is sufficient even for a calculation done with dimensional regularisation (as long as one is only concerned with physical quantities in four dimensions). And a similar statement applies also to the vector part of the amplitude (3.1) covered in our computations.
For a quantitative estimate of the size of the non-Drell-Yan contributions, we have studied their numerical impact on the total cross section of ZH production in section 7.1. We find that these non-Drell-Yan processes give about three orders of magnitude smaller contributions compared to that of the s-channel Drell-Yan-like processes. Nevertheless, the availability of these polarised amplitudes for the process (2.1) would allow us to combine the b-quark-induced ZH production with the subsequent decay of the Z boson with full spin correlations accounted for.