Polarized $q \bar{q} \rightarrow Z +$Higgs amplitudes at two loops in QCD: the interplay between vector and axial vector form factors and a pitfall in applying a non-anticommuting $\gamma_5$

We consider QCD corrections to two loops for the polarized amplitudes of $q{\bar q}\to Z +$ Higgs boson. First we show how the polarized amplitudes of $b \bar{b} \rightarrow Z h$ associated with a non-vanishing $b$-quark Yukawa coupling and a scalar or pseudoscalar Higgs boson $h$ can be built up solely from vector form factors (FF) of properly grouped classes of diagrams, bypassing completely the need of explicitly manipulating $\gamma_5$ in dimensional regularization (up to a few"anomalous", i.e., triangle diagrams). We determine the contributions of the triangle diagrams in the heavy top limit. We present the analytic results of the vector FF and the triangle-diagram contributions to the axial vector FF, which are sufficient for deriving the two-loop QCD amplitudes for $b \bar{b} \rightarrow Z h$ with a CP-even and CP-odd Higgs boson $h$. We derive the respective Ward identity for these amplitudes, which are subsequently verified to two-loop order in QCD using these FF. In addition, the FF of a class of corrections to $q \bar{q} \rightarrow ZH$ proportional to the top-Yukawa coupling are obtained analytically to two-loop order in QCD in the heavy-top limit using the Higgs-gluon effective Lagrangian where the top quark is integrated out. We address a pitfall that occurs when applying the non-anticommutating $\gamma_5$ prescription to this class of contributions that has been overlooked so far in the literature. We attribute this issue to the fact that the absence of certain heavy-mass expanded diagrams in the infinite-mass limit of a scattering amplitude with an axial vector current depends on the particular $\gamma_5$ prescription in use.


Introduction
A detailed investigation of the kinematical and dynamical properties of the 125 GeV Higgs boson discovered at the Large Hadron Collider (LHC), i.e., its kinematic profile as well as how it interacts with (other) known fundamental particles, remains among the major research topics of the current and future physics programs.In particular, the production of the Higgs boson in association with a massive electroweak vector boson, known as the VH process, plays an important role in the exploration of Higgs physics at the LHC, both for a precise study of the Higgs boson's couplings to Standard Model particles and for probing new physics.For instance, it supplies the main production channels behind a recent experimental triumph, the direct observation of the Higgs-boson decay to a pair of bottom quarks by the ATLAS and CMS experiments [1,2].This was largely possible owing to the fact that the presence of the associated vector boson offers means to substantially reduce the Standard Model backgrounds, for instance by requiring a large transverse momentum of this vector boson [3].Given foreseeable upgrades in experimental precision at future collider experiments, e.g. the high luminosity LHC program, it is very desirable to have a precise knowledge about the VH process at hadron colliders on the theoretical side as well.
In ref. [13] two of the authors of this article have presented the analytic results of the two-loop massless QCD corrections to the b-quark-induced ZH process pertaining to a non-vanishing b-quark Yukawa coupling λ b .The amplitudes for the polarized Z-boson states were constructed following the prescription of ref. [14].For the treatment of the axial vector current vertex in dimensional regularization the prescription of refs.[15,16] was used, where the γ 5 no longer anticommutes with the Dirac matrices in D dimensions.
We consider, in this article, the QCD corrections to two loops for the polarized amplitudes of q q → Z+ Higgs boson.First we determine the amplitudes for b b → Zh at two loops associated with a non-zero b-quark Yukawa coupling in analytic fashion, both for a scalar (h = H) and a pseudoscalar (h = A) Higgs boson and for a polarized Z boson.We use the well-known fact that an anticommuting γ 5 , denoted by γ AC 5 in this article, can be used in D-dimensional computations [17][18][19][20][21][22] that do not involve the Adler-Bell-Jackiw (ABJ) anomaly [23,24].We show how the respective "non-anomalous" contributions that correspond to diagrams where the Z boson couples to an open quark line, can be built up solely from vector form factors of properly grouped classes of diagrams whose computation does not involve the axial vector current from the outset.Furthermore, we derive and verify explicitly the Ward identities for the QCD corrections to the b-quark Yukawa coupling dependent contributions to b b → Zh, h = H, A, using these vector form factors.In addition, we determine the two-loop "anomalous" contributions to b b → Zh corresponding to diagrams that involve b-and t-quark triangles and the axial vector current.
The QCD corrections to the quark-annihilation induced process q q → Zh include a class of diagrams where the Higgs boson couples directly to a closed top-quark loop that start to appear at the two-loop order in QCD.As the second topic of this work, we calculate a subset of these contributions, the so-called class-I diagrams, to O(α 3 s ) at the amplitude level in the heavy-top limit using the Higgs-gluon effective Lagrangian of Higgs effective field theory (HEFT) where the top quark is integrated out.Our motivation for presenting the discussion of these contributions here is that their computation within HEFT led us to uncover a pitfall in the application of a non-anticommuting γ 5 to this class of contributions: we found that extra counterterms are needed in addition to those that are known from the prescription of [15,16] in order to get correct results that obey respective Ward identities when the axial vector current is dimensionally regularized with a non-anticommuting γ 5 .We stress that this is not to be regarded as a contradiction to the prescription of ref. [15,16], because after all this is a phenomenon that is associated with the use of HEFT in the calculation of this class of diagrams.We attribute the need for such additional counterterms to the fact that the absence of certain heavy-mass expanded diagrams in the infinite-mass limit of a scattering amplitude with an axial vector current actually depends on the particular γ 5 prescription in use.
The article is organized as follows.In section 2, we consider contributions to the polarized amplitude of b b → Zh (h = H, A) that are proportional to the respective b-quark Yukawa coupling.First, we determine the form factor (FF) decomposition of the amplitude of b b → ZH and show how in the case of non-anomalous contributions the axial vector FF can be obtained from the FF associated with the vector current.The latter FF are computed to two-loop order in QCD.Moreover, we determine the b-and t-quark contributions to the axial FF from the "anomalous" diagrams with quark triangles in the limit m t → ∞.Then we show how the vector and axial vector FF of the respective amplitude associated with a pseudoscalar Higgs boson can be obtained from the FF associated with a scalar Higgs.In section 3, we derive the respective Ward identity that the amplitude of section 2 for h = H, A must satisfy and check these identities using the FF computed before.In section 4, we consider a class of contributions, the so-called class-I contributions, to q q → ZH proportional to the top-Yukawa coupling within HEFT where the top quark is integrated out.First, we determine the vector FF and non-anomalous axial vector FF to two loops using an anticommuting γ 5 in dimensional regularization and show that the respective Ward identities are satisfied.In section 5, we recompute the non-anomalous class-I axial FF and also the contributions from the anomalous diagrams in HEFT using a non-anticommuting γ 5 in D dimensions.We find that both for the non-anomalous and anomalous axial vector FF, an additional counterterm is required in order to satisfy the respective (non)anomalous Ward identity.Finally, we investigate the axial vector part of the non-anomalous class-I diagrams with a non-anticommuting γ 5 in the full theory without taking the limit m t → ∞.We conclude in section 6.
2 Axial vector form factors from vector counterparts

Production of a scalar Higgs boson H
In this section, we consider first, for definiteness, the production of a scalar Higgs boson, H, in association with a massive vector boson, Z, through bottom quark anti-quark annihilation: The four-momenta of the particles in (2.1) satisfy the on-shell conditions , where m Z and m H are the masses of the Z and Higgs boson, respectively.The Mandelstam variables are We keep a non-zero Yukawa coupling only for the b quark, which otherwise is taken to be massless, and for the top quark whose contributions are considered in the infinite mass limit.By employing an elegant methodology, we recompute in this section the two-loop QCD corrections to the non-Drell-Yan type diagrams of the process (2.1), shown at tree level in figure 1, that depend on the b-quark Yukawa coupling λ b , in QCD for 5 massless flavors and determine also the contribution of the top quark in the infinite mass limit.These contributions form a gauge-invariant set.The amplitude of these contributions can be parameterized, including QCD corrections to two loops, by ( The symbol Γ µ axi represents a matrix in the spinor space with one open Lorentz index µ that may be carried by either the Dirac matrix γ µ or one of the external momenta involved.It is the sum of the contributions from the vector and axial vector couplings of the Z boson.For purposes discussed below we have factored out λ b and the vector and axial vector couplings of the b quark to the Z boson, denoted by g V,b and g A,b , respectively. At the end of this section, we consider the production of a pseudoscalar Higgs boson A analogous to (2.1) and discuss how the respective scattering amplitude to two-loop order in QCD analogous to (2.3) can be obtained from the vector form factors that determine the amplitude (2.3) and which will be computed next.

The interplay between axial vector and vector form factors
Unlike ref. [13], where a non-anticommuting γ 5 was used for the Z boson axial vector coupling, we compute here all non-anomalous contributions to the amplitude (2.3) to two loops using an anticommuting γ AC 5 in D dimensions.First we consider all two-loop diagrams where i) only the Higgs boson and ii) the Higgs and Z boson are radiated from a closed massless quark loop, such as those shown in figure 2. These contributions vanish because they involve a trace with an odd number of Dirac γ matrices.Thus, in all non-vanishing non-anomalous contributions to two loop order in QCD, the Z boson couples to the open b-quark line and γ AC 5 contained in Γ µ axi can be anticommuted next to an external b-quark spinor.These non-anomalous diagrams can be further divided into two classes, denoted in the following by class-ZH and class-HZ, which correspond to the QCD corrections to the tree-level diagrams (A) and (B), respectively, of figure 1.The reason for this separation is simply the fact that due to the presence of the chirality-flipping Yukawa interaction on the b-quark line a relative minus sign is generated between these two contributions when the anticommuting γ AC 5 is pushed next to the same external b-quark spinor.This is also the reason why here the non-anomalous axial vector form factors are not identical to the vector ones, see below.
Let us turn off for a moment the axial vector coupling of the Z boson and consider only M vec in (2.3).We can decompose it in terms of form factors as follows: In this form factor decomposition, which reflects the chirality flip along the massless bquark line, we have taken into account the equations of motion for the on-shell massless spinors v(p 2 ) and u(p 1 ), but have not confined ourselves to the physical polarization states of the Z boson.The projectors for obtaining the vector form factors defined in (2.4) from Γ µ X (X = ZH, HZ) are derived and explicitly given in ref. [13].In the absence of the axial current, the four basis structures in (2.4) are linearly complete in D dimensions for M vec in (2.3), regardless of the QCD loop order.
The projection of the two sets of vector form factors F i,ZH and F i,HZ encounters no subtlety at all, and their renormalization is standard with the details given in ref. [13].The complete vector form factors, F i,vec , defined by are given by F i,vec = F i,ZH + F i,HZ for i = 1, 2, 3, 4 . (2.6) Restoring the axial vector couplings of the Z boson, the amplitude M axi defined in eq.(2.3) consists of a "non-anomalous" and "anomalous" contribution: These contributions can be decomposed into form factors, in analogy to eq. (2.5).Using γ AC 5 in D dimensions, the non-anomalous axial form factors F i,axi(ns) , defined by are obtained by (2.9) The appearance of the relative minus sign in (2.9) is explained above eq.(2.4).
We have checked that the unrenormalized, unsubtracted vector form factors F i,vec defined in (2.6) agree with those obtained in ref. [13] to two-loop order in D dimensions.Without surprise, the axial form factors F i,axi composed in (2.9) are indeed different from those defined in ref. [13] in their bare form.However, after carrying out the ultraviolet (UV) renormalizations and infrared (IR) subtractions, the finite remainders of these form factors are identical to those given in ref. [13]. 1 Here we remark that the UV renormalization needed for F i,axi(ns) in (2.9) is identical to that of the vector form factors, which is different from what was done in ref. [13] regarding the axial form factors.We have thus crosschecked the previous computation of the non-anomalous axial part of ref. [13] (where a non-anticommuting γ 5 [15] was used).The UV renormalized results of the partial form factors (FF), F i,ZH and F i,HZ , which are the building blocks for composing the complete non-anomalous vector and axial matrix elements, are provided as ancillary files with the arXiv submission.Details about conventions and variables of the analytic expressions can be found in the ReadMe.txtsubmitted along with these files.We present the results to two-loop order, more specifically, the first three coefficients of the partial FF in the expansion where X ∈ {ZH, HZ}, a s = α s /(4π) and µ R is the renormalization scale.

Contributions from diagrams involving quark triangles
The two-loop "anomalous" contributions to the process (2.1) involving quark triangles can be represented by six Feynman diagrams, shown in figure 3. Furry's theorem tells us 1 There is a relative minus sign between the fourth axial form factor defined by (2.9) and the corresponding one defined by equation (3.19) in ref. [13], because the latter one corresponds to the Lorentz structure v(p2) γ µ γ5 / q 1 u(p1).
where m t is the top quark mass.The relative minus sign is due to the definition of M axi in eq.(2.3) (where we pulled out an overall coupling factor g A,b from the contribution of the axial vector) and the opposite weak isospins of the b and t quark.The amplitude M µ axi(s) can be decomposed into form factors F i,axi(s) , in complete analogy to eq. (2.8).
The top-loop induced diagrams that contribute to the process (2.1), all starting at O(α 2 s ), can be divided into a set which is independent of t-quark Yukawa coupling λ t , i.e. the diagrams of figure 3 with the t quark circulating in the triangle loop, and another set dependent on the top Yukawa coupling that will be addressed in section 4. Concerning the contributions to M t,µ axi(s) (m t ), we show below that the last two pentagon-triangle diagrams of figure 3 vanish in the limit m t → ∞.The contributions of the other four box-triangle diagrams are not zero in this limit, but they can be captured by a decoupling matching constant introduced for the effective b bZ interaction, resulting from integrating out the virtual top quark in these diagrams.
Applying the expansion-by-graph procedure (see, e.g.ref. [25] and references therein) to one representative pentagon-triangle diagram, as shown on the l.h.s. of figure 4, in the heavy-top limit m t → ∞, one ends up with the sum of two classes of heavy-mass expanded contributions, indicated by the two terms on the r.h.s. of this figure.They correspond to the hard-hard loop momentum region with |k 2 | ∼ |k 1 | ∼ m t and the soft-hard loop The dotted lines in the asymptotic irreducible graphs [25] on the left of "⊗" indicate external momenta which can all be set to zero if only the leading contribution in the limit m t → ∞ is to be kept.The graphs on the right of "⊗" denote the corresponding co-graphs in the heavy-mass expanded result.momentum region with2 |k 2 | < |k 1 | ∼ m t , respectively.Let us first look at the contribution from the hard-hard momentum region, corresponding to the first term on the r.h.s. of figure 4. The asymptotically irreducible graph [25] here is a two-loop "vacuum" graph that has the mass-dimension -1, and after expansion in small ratios it depends only on m t in the infinite m t limit and polynomially on the external momenta which can be set to zero.This is sufficient to ensure that there is no non-vanishing contribution from this graph in the limit m t → ∞.
We then move on to the second term, originating from the soft-hard momentum region.Here the subgraph that is to be expanded is a triangle t-quark loop that has the same topology as the VVA triangle diagrams drawn in figure 5. Here, counting of the mass dimension alone is no longer sufficient to tell us whether or not this subgraph has a nonvanishing limit when m t goes to infinity.Therefore we compute this subgraph explicitly and take the analytic expressions of all one-loop integrals involved from their implementations in PackageX [26].
Leaving all Lorentz indices from the three gauge vertices open in the diagrams of figure 5, a rank-3 Lorentz tensor amplitude Γ µ 1 µ 2 µ (l 1 , l 2 , m t ) is thus introduced for the sum of these two one-loop diagrams.The tensor Γ µ 1 µ 2 µ (l 1 , l 2 , m t ) can be further split into a vector (spin-1) part and a scalar (spin-0) part with respect to its axial-vector current index µ: where and the physical polarisation projector g µν − qµqν q 2 projects Γ µ 1 µ 2 µ (l 1 , l 2 , m t ) onto the space of vector polarisations (indicated by the subscript v).We refrain from going into the technical details of the one-loop computation involved, e.g.via the form factor decomposition approach, but merely point out the following explicitly verified fact: keeping the momenta l 1 and l 2 fixed, the rank-3 Lorentz tensor amplitude Γ µ 1 µ 2 µ (l 1 , l 2 , m t ) vanishes in the limit m t → ∞.This holds true for , respectively.This implies that the second term on the r.h.s. of figure 4 associated with the soft-hard region vanishes in the limit m t → ∞.Therefore the two pentagon-triangle diagrams do not contribute in this limit.
As an aside we briefly comment on an interesting implication related to the above statement about the triangular one-loop subgraph, which, albeit not new, does not seem to be common knowledge.At the on-shell kinematic configuration l 2 1 = l 2 2 = 0, the term Γ µ 1 µ 2 µ v (l 1 , l 2 , m t ) vanishes completely in four dimensions, owing to the Landau-Yang theorem (because the color factor here is trivial), regardless of the mass of the top quark.The only non-vanishing piece at this configuration is Γ µ 1 µ 2 q s (l 1 , l 2 , m t ), associated with the scalar polarisation state.For m t = 0, Γ µ 1 µ 2 q s (l 1 , l 2 , m t ) is given precisely by the ABJ anomaly.What we would like to emphasize here is that Γ µ 1 µ 2 q s (l 1 , l 2 , m t ) from the oneloop diagrams of figure 5 vanishes in the heavy top limit m t → ∞, as a consequence of a non-trivial cancellation between the pure m t -independent quantum anomalous contribution (i.e. the ABJ anomaly) and the non-vanishing limit of the m t -dependent "classical" contribution at m t → ∞.In other words, the non-vanishing m t → ∞ limit of the non-decoupling m t -dependent "classical" contribution to Γ µ 1 µ 2 q s (l 1 , l 2 , m t ) is exactly opposite to the m t -independent ABJ anomaly contribution.Indeed, this point can be checked straightforwardly by a direct projection, because Γ µ 1 µ 2 q s (l 1 , l 2 , m t ) contains just one Lorentz covariant structure.The corresponding unique projector is indifferent to whether the contribution comes from the m t -independent quantum anomalous part or the m t -dependent "classical" part, and both will be projected out on the same footing.
Let us come back to the discussion of the heavy-top limit of the remaining four boxtriangle diagrams of figure 3 where the results about the triangle subgraphs given above will be used again.Applying the expansion-by-graph procedure to one representative diagram of this topology, as shown on the l.h.s. of figure 6, in the heavy-top limit, one ends up with the sum of two classes of heavy-mass expanded contributions, similar to the case of the pentagon-triangle diagrams.Note that here the b-quark tree-propagator is amputated before applying the heavy-mass expansion procedure, resulting in an off-shell b-quark leg (indicated by the double-thick line in figure 6) whose momentum is considered small compared to m t .The second term on the r.h.s. of figure 6 associated with the soft-hard loop-momentum region vanishes in the limit of infinitely heavy top quark due to the same reason as mentioned above.However, the first term that originates from the hard-hard momentum region has an asymptotic irreducible graph that is a two-loop boxtriangle "vacuum" graph with zero mass-dimension and turns out to have a non-vanishing limit at m t → ∞.As should be clear from the diagrammatic illustration of the heavymass expansion result given in figure 6, this contribution can be presented by a set of local composite operators determined by the co-graph, an effective bbZ-vertex, and the "vacuum" graph in front of it.One key feature of the expanded "vacuum" graph (associated with the hard-hard momentum region) is that its dependence on the external momenta {p 1 , p 2 , q 1 , q 2 } is purely polynomial.In other words, this special vacuum loop amplitude has a regular limit for vanishing external momenta.Since we are only interested in the leading contribution in 1/m t at m t → ∞ where all external momenta in this expanded "vacuum" graph can be put to zero, it is then not hard to see that this expanded box-triangle "vacuum" graph will lead to the same expression regardless of whether or not the b-quark line is on-shell.This non-vanishing infinite-m t limit of the (properly renormalized) r.h.s. of figure 6 can thus be captured by introducing an effective b bZ interaction with a decoupling matching constant C (s) bbZ (being independent of the kinematics), denoted by accounting for integrating out the top-loop from these diagrams at the leading power in the heavy-top expansion. 3Based on the information above, C (s) bbZ (m t ) can be extracted from the heavy-top limit of the two-loop anomalous QCD corrections to the renormalized (on-shell) quark form factors computed in ref. [27], and is given by where g A,t denotes the axial coupling between the top quark and Z boson, C F = (N 2 c − 1)/(2N c ) and T R = 1/2.We note in passing that this expression is equal to twice the order a 2 s Wilson coefficient in front of the effective interaction q qA between a pair of light quarks and a pseudoscalar Higgs boson A resulting from integrating out the heavy top quark, computed in ref. [28].
In summary, the non-vanishing infinite-m t limit of M t axi(s) (m t ) of the diagrams of figure 3 with a t-quark triangle is given by the axial part of the tree-level amplitude of the process (2.1) multiplied by the coefficient (2.13).We attach the analytic results of the properly renormalized axial form factors of the complete contribution M axi(s) of (2.11), defined in complete analogy to (2.8), as ancillary files in Mathematica format along with the arXiv submission.(See ref. [13] for the renormalization constants used.)Finally, we remark that after incorporating the non-vanishing contribution M t axi(s) (m t ), the explicit µ R dependence cancels in the amplitude M axi(s) = M b axi(s) − M t axi(s) (m t ) at the two-loop order.(2.14) We remark that in this case the two diagrams of figure 1 represent the complete tree-level amplitude (because a tree-level ZZA vertex does not exist).Using γ AC in D dimensions the amplitude

Production of a pseudoscalar
proportional to the pseudoscalar b-quark Yukawa coupling λb can be constructed to twoloop order QCD with the above vector form factors F i,HZ and F i,ZH .Notice that the subscripts, vec and axi, of the amplitudes in (2.15) refer to the respective coupling between the Z boson and the b quark.It is straightforward to see that the set of Lorentz basis structures, and consequently the corresponding form factor projectors previously used, in the form factor decomposition for the amplitude involving the CP-even Higgs boson H exchange their roles in the current process (2.14).For instance, the set of Lorentz structures for decomposing M vec appearing in (2.5) is now the one needed for the form factor decomposition of M A axi defined in (2.15).Consequently the corresponding "vector" projectors will project out the axial form factors of the production amplitude of a pseudoscalar Higgs boson A in association with a Z boson.It can be verified to two-loop order, conveniently using γ AC 5 , that the following relationships hold among the form factors: (up to an overall phase factor depending on the parameterization of the general Yukawa couplings.For our choice (3.1) it is an overall factor i. It is suppressed here for simplicity.).The results of F i,ZH and F i,HZ were obtained in section 2.2.Furthermore, for the contribution from the triangle diagrams analogous to figure 3, we have and the computation of F i,axi(s) of M axi was discussed in section 2.3.
3 The Ward identities for b b → ZH, ZA In this section we derive and subsequently check, using the form factors F i,HZ and F i,HZ obtained above, the Ward identity for the QCD virtual corrections to two loops to the process (2.1) and (2.14) with a CP-even and CP-odd Higgs boson, respectively, keeping a non-vanishing Yukawa coupling only for the b quark.Let us emphasize again that no Higgs bremsstrahlung from the Z boson is considered here.

Derivation of the Ward identity
The classical Lagrangian that encodes all the aforementioned information reads as where G a µν denotes the gluon field strength tensor, gives rise, at the classical level, to the following equation: If one had kept the b quark massive, then there will be also a b-quark mass-dependent term appearing on the right-hand side of (3.5).At the quantum level the current J µ 5 suffers from the ABJ anomaly [23,24], i.e., the term −a s µνρσ G a µν G a ρσ /2 appears4 in addition on the right-hand side of (3.5).Proper renormalization of fields and interactions involved are understood wherever needed.We will come back to this later.
Next we derive a relation between the S-matrix element of the process (2.1), and likewise for (2.14), and a respective S-matrix element where the Z boson has been removed.We call this relation a Ward identity, although the axial vector current involved in this relation is not even partially conserved (cf.eq.(3.5)).
CP-even Higgs boson: The corresponding S-matrix element is where µ Z denotes the polarization vector of the Z boson.We consider here only the socalled non-anomalous QCD contributions to the S-matrix element (3.6), namely eq.(3.5) will be used without the ABJ anomaly term.
Applying the LSZ reduction formalism [29] to the S-matrix element (3.6) one obtains5 where disc.denotes disconnected terms that do not contribute here.Let us further denote Contracting M µ with the four-momentum of the Z boson we obtain where we used the conservation of the vector current, the divergence of the axial vector current, eq.(3.5) (in the absence of ABJ anomaly), and translation invariance.From eqs. (3.6) -(3.9) we obtain the following relation which we call a Ward identity: This relation holds to all orders in the QCD coupling (in the absence of the anomalous diagrams) with the dynamics specified by the Lagrangian (3.1) in its renormalised form.
Notice that the kinematics of the matrix element on the right-hand side of this equation obeys p 1 + p 2 − q 2 = q 1 , where q 2 1 = m 2 Z .This is due to the external momentum insertion q µ 1 introduced by the local composite operator b(0)iγ 5 b(0)H(0) (which is understood to be normal-ordered).To lowest order in perturbation theory one gets We remark that the overall sign of the right-hand side follows from the definition of the initial two-fermion state.The three-point pseudoscalar vertex (3.11) represents an incoming b and b quark with four-momenta p 1 and p 2 , respectively, and an outgoing Higgs boson with four-momentum q 2 .It is depicted in figure 7.

CP-odd Higgs boson:
The corresponding S-matrix element is (3.12) The respective Ward identity for (3.12) can be derived in completely analogous fashion using (3.5).We obtain Evaluating the right-hand side of (3.13) to lowest order we get The remarks made below (3.10) and (3.11) apply also here.

Checking the Ward identity
Next we verify the Ward identity (3.10) for a CP-even Higgs boson H at the level of the UV renormalized and IR subtracted finite remainders in four dimensions.From the discussion of section 2.2 we see that the virtual (non-anomalous) two-loop QCD corrections to the amplitude of the process (2.1) are given by where the previously calculated analytic expressions of the vector form factors F i,ZH and F i,HZ will be inserted and all spinor products are to be evaluated in four dimensions.The left-hand side of the Ward identity (3.10) is obtained by contracting the amplitudes (3.15) with the four-momentum of the Z boson.We get ) The form factors F i,ZH and F i,HZ can be UV-renormalized and IR-subtracted as outlined in detail in ref. [13], which subsequently leads to their respective finite remainders in four dimensions.Inserting these finite remainders of F i,ZH and F i,HZ into (3.17) and (3.18), respectively, we obtain the finite remainders of the left-hand side (3.10) in four dimensions.
In particular we have verified that q 1 • M vec = 0. (This holds, of course, already before renormalization and subtraction.) The right-hand side of (3.10) consists, up to the factor −2g A,b λ b , of the tree-level threepoint pseudoscalar vertex (3.11) and its one-loop and two-loop QCD virtual corrections, albeit with a special kinematic configuration as explained above.Example diagrams are shown in figure 7.All contributions are proportional to the Lorentz structure v(p 2 ) γ 5 u(p 1 ).We used here an anticommuting γ AC 5 for the pseudoscalar vertex.After carrying out the UV renormalization and IR subtraction of these QCD virtual corrections, we obtained an expression that agrees analytically with the aforementioned finite remainder of the lefthand side of (3.18).
In the case of a CP-odd Higgs boson, where the vector and axial form factors, i.e.F A i,vec and F A i,axi(ns) , are given by (2.16) in terms of F i,ZH and F i,HZ , it is straightforward to verify the Ward identity (3.13) in completely analogous fashion.To be more specific, with the parameterization of the Yukawa couplings as in (3.1), we have F A j,vec = i F j,vec and F A j,axi(ns) = i F j,axi(ns) .The two sets of form factor-factor decomposition bases of the amplitudes involving H and A differ just by an additional γ 5 sandwiched between spinors.Likewise, the r.h.s. of the Ward identities (3.10) and (3.13) differ by the factor iγ 5 sandwiched between spinors, as far as the Lorentz structure is concerned.Thus the form factors (2.16) fulfill the Ward identity (3.13).
4 Top-Yukawa coupling dependent top-loop contributions to q q → ZH in Higgs effective theory Starting from O(α 2 s ) in QCD, a new class of two-loop diagrams contributes to the quark antiquark annihilation initiated Zh process where the Higgs boson couples directly to a closed top-quark loop, and hence the top-quark Yukawa coupling gets involved.This was studied for h = H for instance in ref. [4,5] by making use of asymptotic expansions in the heavy-top limit.In ref. [5], the two-loop virtual top-quark contributions to q q → ZH proportional to λ t were classified into two sets: class-I (with examples shown in figure 8) and class-II, depending on whether the Z boson couples to the external light quark in the initial state or to the virtual top-quark loop (from which the Higgs boson is radiated), giving rise to different electroweak coupling factors.These contributions were not covered in ref. [13] where the computations were made in n f = 5 massless QCD with a non-vanishing Yukawa coupling for the b quark only (as discussed in the preceding sections except for the section 2.3 where the top-induced triangle diagrams are included in addition).As the second part of the work presented in this article, we compute the contributions of the class-I diagrams to O(α 3 s ) in the heavy-top limit using the Higgs effective field theory (HEFT) to two-loop order.We confine ourselves in the following to a scalar Higgs boson h = H.Our motivation for presenting these contributions here in some detail is a problem in applying a non-anticommuting γ 5 that we encountered in the computation of this class of contributions.
We parameterize in the following the general Yukawa coupling λ t of a CP-even Higgs boson to the top quark by where −m t /v with v = 246 GeV is the SM top-Yukawa coupling and the dimensionless parameter c t depends on the specific Higgs model.

HEFT and UV renormalization
In ref. [5] it was pointed out that applying the heavy-mass expansion procedure to the class-I diagrams of figure 8 that leads to terms similar to those depicted in figure 4, albeit with Z and H exchanged, the expanded terms featuring the effective q qZH or q qH vertex vanish to leading power in 1/m t .Only the terms that involve the effective Higgs-gluon-gluon (Hgg) vertex contribute in the infinite m t limit.Thus, the leading approximation in powers of 1/m t of the diagrams figure 8 can be described with the Higgs effective Lagrangian (see e.g.ref. [30]) where the top quark is integrated out.The diagrams of figure 8 are then reduced to the one-loop diagrams shown in figure 9.
Figure 9.The class-I diagrams of the top-quark contributions to q q → ZH in the limit m t → ∞ in the leading power approximation.The black blob indicates the effective Hgg vertex.
However, as will become clear at the end of the next section, it turns out that the validity of this point depends on the particular γ 5 prescription in use.To be more specific, it is not true for the axial current regularized using a non-anticommuting γ 5 , which can be deduced from our computations described below.This is one of the key results conveyed through the remaining sections.
In the Higgs effective theory, where the top-quark degrees of freedom are integrated out, the Lagrangian density that encapsulates the interaction between the scalar Higgs boson and gluons is given by (neglecting terms that are not relevant here) where c t and v are defined in and below eq.(4.1), respectively, and C H denotes the Wilson coefficient that is determined for a Standard Model Higgs boson by matching the effective n f = 5 flavor theory to the full (n f + 1)-flavor theory order-by-order in the QCD coupling.
To second order in a s (µ 2 R ) it is given by [31,32] where µ R is the renormalization scale and C A , C F denote the quadratic Casimir operators of the SU(N c ) color gauge group in the adjoint and fundamental representations, respectively.The effective operator (4.2) must be renormalized, in addition to performing the QCD coupling renormalization (done in the MS scheme), in order to get rid of all the UV poles appearing in the scattering amplitudes.This is achieved by with the operator renormalization constant [33][34][35]] The coefficients of the QCD β-function are given by Figure 10.Examples of the class-I two-loop diagrams proportional to λ t involving only an open massless quark line q for q q → ZH in HEFT.

Form factors of the class-I contributions using an anticommuting γ 5
We consider now the class-I λ t -dependent contributions to q q → ZH in HEFT.The leadingorder contributions are depicted in figure 9. Examples of the two-loop non-anomalous QCD corrections that involve only an open massless quark line q where the Z boson couples to are shown in figure 10, and the two-loop QCD corrections that involve in addition a closed quark loop where the Z boson couples to are displayed in figure 11.We denote the corresponding contributions to the amplitude by A in order to distinguish it from the contributions discussed in section 2: where We note that, by construction, the index contraction between these projectors (4.9) and the ε * µ -stripped amplitude is to be done with the D-dimensional space-time metric tensor g µν .
The UV renormalization of the amplitude (4.8) proceeds as explained in section 4.1.The technical aspects of the computation of the vector form factors closely follow the steps as explained in detail in ref. [13].We remark that the two-loop amplitudes (or form factors) (4.8) involves 117 master integrals, for which we take the analytic expressions computed in ref. [36] available in HepForge [37] in computer readable format.The UV renormalized vector form factors F i,vec at one as well as two loops were checked for exhibiting the universal infrared structures [38][39][40][41][42].This serves as a strong check of our computations.
The expansion of the UV renormalized form factors in powers of a s is defined by The analytic results of these UV renormalized vector form factors are too lengthy to be presented here, but they can be provided upon demand from the authors 6 .We decompose the amplitude A axi associated with the axial vector current as follows: where the contribution A axi(ns) covers the one-and two-loop diagrams of figures 9 and 10, respectively, while the term A axi(s) results from the non-vanishing two-loop diagrams of figure 11.The dependence of this term on g A,b only will be discussed in section 5.2 below.The axial components of the form factors are defined by decomposing A axi in analogy to (4.8) as follows (we drop here the additional subscript 'ns' or 's' for ease of notation): and these form factors are expanded in powers of a s in analogy to the vector counterparts in (4.11).
For the computation of the non-anomalous part A axi(ns) , respectively the form factors F 2,axi(ns) , we use an anticommuting γ AC in D dimensions.Then A axi(ns) respects chiral invariance which implies that In order to check eq. ( 4.14) we derive the projectors that correspond to the decomposition (4.13) (with an anticommuting γ AC 5 ).With these projectors we computed the F (l) i,axi(ns) at one and two loops (l = 1, 2) and confirm eq.(4.14).
For the anomalous two-loop diagrams, where only axial vector part survives, we employ a non-anticommuting γ 5 in dimensional regularization [15,16,[43][44][45].Our results will be discussed in section 5.2.First we cross-check whether the correct non-anomalous axial form factors (i.e., those given by eq.(4.14)) are obtained using a non-anticommuting γ 5 as employed in the literature.Here we found something quite intriguing, which we now turn to in the following section.
5 A pitfall in applying a non-anticommuting γ 5 to λ t -dependent contributions to q q → ZH in HEFT In our attempt to compute the two-loop QCD contributions to the class-I diagrams 7 in HEFT, shown in figure 9, we noticed a technical pitfall in applying a non-anticommuting γ 5 .To our surprise, this γ 5 issue manifests itself already in the LO (one-loop) contributions.

5.1
The class-I axial form factors computed using a non-anticommuting γ 5 Our computation of the axial vector form factors, F i,axi , as introduced in (4.13), using a non-anticommuting γ 5 follows closely ref. [13], especially with regard to the axial form factor decomposition (although a new set of chirality-preserving form factor decomposition basis is needed here).The non-anticommuting γ 5 in dimensional regularisation is defined by [43,44] where the Levi-Civita symbol ε µνρσ is treated according to refs.[15,45,46].The usage of this definition has profound implications in higher-order computations involving an axial vector current.One of the main messages conveyed through refs.[13,14] is that even if the loop amplitudes are not defined or regularized strictly in the 't Hooft-Veltman scheme, expressions for projectors derived in four dimensions are still sufficient and lead to correct results (for physical observables), for projectors corresponding to form factors and to (linearly) polarised amplitudes alike.This is particularly helpful in evaluating loop amplitudes that involve axial vector currents (if a non-anticommuting γ 5 is used).If we use an anticommuting γ 5 , then as just discussed in section 4.2, we get for i = 1, 2, 3, 4 to two loops (leaving the anomalous diagrams aside), where the superscript AC indicates the use of an anticommuting γ 5 .Accordingly, conservation of the non-singlet light-quark current implies where A µ axi(ns) is given by (4.13) with F AC i,axi(ns) inserted.Alternatively, if we use the nonanticommuting γ 5 (NAC) (5.1), then the vector part does not change of course, but for the axial FF we find at LO in the four-dimensional limit ( = (4 − D)/2 = 0): and Consequently, this difference in the fourth form factor leads to the violation of the Ward identity q 1,µ A µ,NAC axi(ns) = 0 (5.6) already at LO.With the explicit analytic expressions of these form factors at hand, we find that the expected relations (4.14) can be restored at LO if we introduce an additional amendment term, i.e. subtraction term, that is to be added to A µ,NAC axi(ns) .We denote this term by where the constant factor C ≡ a s (−4C F ) C H v collects the overall perturbative power a 2 s of the LO amplitude (and the effective coupling prefactors) and [γ µ γ 5 ] L denotes the axial vector current matrix renormalized according to the prescription [15,16].(At this order no renormalization of refs.[15,16] gets involved.)The additional renormalization constant Z h 5 (a s ) = 1 + O(a s ) has an expansion in a s which we will determine explicitly to the first order in a s from our NLO computation to be presented later.Diagrammatically, this amendment term can be viewed as introducing a four-point local composite operator corresponding to diagram 12 with a multiplicative factor Z h 5 (a s ) to be determined order by order in a s .With the extra amendment term (5.7), the desired properties are restored: holds for all four form factors and q 1,µ A µ,NAC axi(ns) = 0 is then also fulfilled.
The source of the inequality (5.5) can be traced at LO to be solely due to the one-loop box diagram, i.e., the left-most one in figure 9, while the contributions of the two triangle one-loop LO diagrams respect F NAC i,axi(ns) = F i,vec (i = 1, 2, 3, 4) in the four-dimensional limit.We emphasize that each of the three LO diagrams of figure 9 is finite.However, the Feynman amplitude of the finite one-loop box diagram contains terms that are separately divergent.The expressions obtained using on the one hand an anticommuting γ 5 and on the other hand, a non-anticommuting one, lead to different D-dependent coefficients in front of these divergent terms, with differences being suppressed by at least one power in (D-4).The crucial point is that the (D-4) difference between these two expressions is not an overall prefactor.It is then not hard to conceive that some non-vanishing evanescent anti-commutators are generated when one shifts the non-anticommuting γ 5 from inside the loop correction of the axial vector vertex to the outside of the loop, which then leads to the observed discrepancy.

The NLO QCD corrections UV renormalization of the non-anomalous HEFT diagrams
We move on and discuss the NLO QCD correction to the aforementioned LO diagrams using a non-anticommutating γ 5 in HEFT where a similar phenomenon happens.Let us first consider the NLO QCD corrections of the non-anomalous type to the leading order diagrams in HEFT, shown in figure 9.These corrections correspond to the two-loop diagrams of order α 3 s shown in figure 10.There are also a few non-zero contributions of the anomalous type, i.e., non-vanishing diagrams involving quark triangles, at this perturbative order, shown in figure 11, which we will discuss separately.All the projectors used in computing these non-anomalous NLO QCD corrections are still the same as those used for the LO diagrams, and the only technical complication comes from performing the UV renormalization of the axial form factors regularized with a non-anticommuting γ 5 .The renormalization procedure described in section 4.1 is sufficient to renormalize the vector form factors which is verified by checking the respective Ward identity.
We use now the prescription of refs.[15,16] for the axial vector current, sometimes referred to as Larin's prescription for short, and therefore the corresponding axial vector current renormalization constants in the MS scheme get involved non-trivially at this order of perturbation theory, in addition to the QCD coupling and operator renormalization.For the non-anomalous NLO QCD corrections to the diagrams of figure 9, the normal form of the Ward identity should still hold for the axial vector current, exactly the same as for the vector counterpart.However, to our surprise, if one just incorporates the usual UV counterterms arising from coupling constant, operator renormalization, and the compensation terms dictated by Larin's prescription, the results for the axial form factors are still wrong, which manifests itself in the following checks.
• The remaining pole structures in the axial form factors renormalized in this way do not match with the prediction according to Catani's IR factorisation formula [38].
• The 0 -order terms of the, renormalized and subtracted, vector and respective axial form factors differ, which subsequently implies the violation of the axial Ward identity for the process at hand.The solution we found, which is now not hard to guess based on the above exposition of how to correct the LO result, is that one should also consistently compute the perturbative contributions to the extra local composite operator given in (5.7) that include both i) the perturbative expansion of Z h 5 (a s ) in a s and ii) the NLO (one-loop) corrections to this local composite operator (where the axial vector current involved is again treated by Larin's prescription).Following this line, we determine the expression for Z h 5,ns (a s ) by imposing the equality between the finite remainders of the vector and (non-anomalous) axial form factors.We get To summarize, the non-anomalous axial amplitude A axi(ns) computed at two loops in HEFT using a non-anticommuting γ 5 is renormalized according to8 A µ,NAC axi(ns) (a s ) = Z ns 5,L (a s )Z ns A,L (a s )Z H (a s ) Âµ,NAC axi(ns) (â s ) + J µ,NAC ns . (5.9) In order to distinguish here bare and renormalized quantities we denote the unrenormalized amplitude and QCD coupling with a hat.The counterterm vertex in (5.9) is given by = Z h 5,ns (a s )Z ns 5,L (a s )Z ns A,L (a s ) C v(p 2 ) γ µ γ 5 u(p 1 ) . (5.10) The constant C is given below eq.(5.7).The quantity Z H is the operator renormalization constant (4.5) in the effective Lagrangian in HEFT; Z ns 5,L and Z ns A,L are the renormalization constants for the non-singlet axial vector current in Larin's prescription [15], which we list here for convenience: If one does not invoke the Larin counterterms (5.11) then, of course, only the product Z h,T 5,ns (a s ) ≡ Z h 5,ns (a s )Z ns 5,L (a s )Z ns A,L (a s ) as a whole can be determined in a chosen renormalization scheme (e.g. in the MS scheme used here) by demanding that the correct physical results (see above) are obtained.

UV renormalization of the anomalous diagrams in HEFT
Regarding the anomalous diagrams at this perturbative order, shown by the non-vanishing ones in figure 11, a similar treatment as in the non-anomalous case can be applied, albeit with a bit more complexity due to the ABJ anomaly [23,24].Let us first discuss the UVrenormalized Ward identity that the non-vanishing four diagrams of figure 11, where only the axial vector current contributes, must satisfy.Only the massless b-quark triangles make a non-zero contribution to the last four diagrams of figure 11.The top quark contribution is omitted as we work here in n f = 5 flavor HEFT.The contributions involve the coupling factor g A,b c t .The operator relation for the ABJ anomaly of the massless axial vector b-quark current J 5µ = bγ µ γ 5 b reads: where G G = − µνρσ G a µν G a ρσ .The subscript R indicates that these composite local operators need to be properly renormalized in order that this operator relation holds.Let us denote the renormalized contribution to A axi(s) as defined in (4.12) from the anomalous Feynman diagrams of figure 11 by Taking the divergence of J 5µ amounts to substituting ε * µ → q 1,µ .By restoring the ABJ anomaly term on the right-hand side of (3.5) and repeating steps similar to eqs.(3.7) -(3.9), we obtain for the anomalous contributions to q q → ZH in HEFT the relation: where the kinematics of the matrix element on the r.h.s.obeys p 1 + p 2 − q 2 = q 1 .As in the case of eq.(3.10) there is an external momentum insertion q µ 1 by the composite operator G G(0) R .Note that, at the two-loop order considered here, both sides of (5.14) are finite upon proper UV renormalization to be determined below.
The matrix element on the r.h.s. of (5.14) can be computed in perturbation theory order by order in a s .The first non-vanishing term corresponds to the one-loop diagrams depicted in figure 13.However, the relation (5.14) does not hold with the expressions we get for both sides.At the order we are considering, i.e.NLO in a s w.r.t. the LO diagrams of figure 9, none of the flavor-singlet axial current Z 5 factors of ref. [16] (e.g.listed in eq.(4.10) in ref. [13]) contributes, because their non-vanishing perturbative terms start at order a 2 s relative to the leading terms.In view of our treatment of the non-anomalous contributions discussed previously, we therefore introduce for the renormalization of the anomalous two-loop contributions of figure 11 the local composite operator given in (5.7) as additional counterterm into the game, albeit with a new undetermined coefficient Z h 5,s (a s ).The expression of Z h 5,s (a s ) given in (5.15) was determined such that the properly UV renormalized singlet axial current has a finite anomaly that does obey (5.12) and the Ward identity (5.14).With the concrete analytic expressions at hand, we obtain We remark that the renormalization of the right-hand side of (5.14) at the perturbative order in question involves the following mixed counterterm: s ) is determined by the requirement of minimally subtracting all poles of the Feynman diagrams in figure 13 (under the convention of setting the coupling factors associated with the effective Hgg vertex and the local operator G G to be one).
To summarize, our explicit computation of the class-I contributions to q q → ZH and their QCD corrections in HEFT, as presented in the preceding sections, shows that if one uses a non-anticommuting γ 5 in dimensional regularization one has to introduce an additional counterterm Z h 5 (a s ) C qR [γ µ γ 5 ] L q R Z µ H into the final properly renormalized effective Lagrangian.To be more specific, the complete form of the properly renormalized effective Lagrangian with a non-anticommuting γ 5 that one should use in computing the QCD corrections to O(α 3 s ) to the class-I contributions to q q → ZH in HEFT reads as where L c + L heff R denotes the renormalized form of L c + L heff given by (3.1) and (4.2).
The explicit perturbative expressions for Z h 5,ns (a s ) and Z h 5,s (a s ) were determined in (5.8) and (5.15) to the first order in a s .Furthermore, κ ns = c t g A,q and κ s = c t g A,b .
We conclude this subsection with a short comment on another class of non-vanishing top-loop contributions to the amplitude of q q → ZH proportional to λ t .In these contributions, which are called class-II in ref. [5] and start at two loops, i.e. order a 2 s , the Z boson couples to a virtual top-quark loop, from which the Higgs boson is also radiated.Charge conjugation invariance dictates that only the axial vector current contributes at order a 2 s .
Thus these contributions are proportional to the of couplings λ t g A,t .It was shown in ref. [5], where the axial vector current was regularized according to Larin's prescription, that in the limit m t → ∞ the non-vanishing part of the class-II contributions involves only one structure, namely qγ µ γ 5 qZ µ H.The investigation of the NLO QCD corrections to these two-loop class-II contributions in the limit m t → ∞ is beyond the scope of this work and will be deferred to a future investigation.

5.3
The class-I Feynman diagrams at two loops without taking the heavy-top limit Our investigations above of the class-I type contributions with HEFT show that correct results (i.e., results that obey chiral invariance and the correct Ward identity) are obtained for both the vector part and the (non-anomalous) axial part with an anticommuting γ 5 .But when employing a non-anticommuting γ 5 in the computation of the same set of (non-anomalous) Feynman diagrams in HEFT, it seems that there are some missing pieces, which calls for additional amendment terms as explicitly determined above.In order to have a better understanding of this issue, we would like to know how these class-I Feynman diagrams behave in the full six-flavor theory without taking heavy-top limit.In particular, we want to check whether the equality (5.2) holds in a computation with a non-anticommuting γ 5 , which then implies that the Ward identity (5.3) is satisfied without the need of additional amendment terms (at two-loop order).There are in total 6 Feynman diagrams for the class-I type contributions at two-loop order with a finite top-quark mass, with samples depicted in figure 8.We generate the (unreduced) symbolic expressions using an extension of GoSam [47,48].After applying the integration-by-parts (IBP) [49,50] relations obtained with Kira [51], they are reduced into linear combinations of 55 master loop integrals.The hardest one is a 7-propagator loop integral that depends on 5 scales (2 Mandelstam variables and 3 physical masses), corresponding to the topology of the first diagram in figure 8, for which no analytic results are known yet.For our purpose of a numerical verification of (5.2), all 55 master integrals are computed using pySecdec [52] at one chosen kinematic point.
Both the vector and axial vector form factors of the two-loop class-I diagrams in the full theory, defined in complete analogy to (4.8) and (4.13), respectively, and projected out using the same projectors as employed in the previous computations done in HEFT, are finite without renormalization or infrared subtraction, because they are the Born level contributions of this type.However, the master integrals involved could and indeed contain spurious poles which all cancel in their final linear combinations making up the form factors.In this numerical check, we tried three different sets of master bases: i) loop integrals without irreducible numerators, ii) loop integrals where irreducible numerators are favoured over denominators raised to powers, and also iii) quasi-finite loop integrals [53], each of which is not unique in general and partially depends on the integral-ordering in use.Concerning the particular choice we take, the quasi-finite master basis exhibits the least spurious poles, albeit still starting from 1/ 2 , while the master basis with numerators performs the worst in the numerical evaluation with pySecdec [52].We note that the master basis without irreducible numerators determined by Kira [51] for our integral family exhibits automatically the feature where their coefficients in the IBP table, and hence the reduced amplitudes, have their denominators' D-dependence factorized from the external kinematics [54,55].Furthermore, we observe that vector and axial form factors not only share exactly the same master basis, but also their respective rational coefficients agree, albeit, only to the leading term in their Laurent expansions around D=4 dimensions 9 , just as in their one-loop counterparts in HEFT discussed at the end of section 5.1.Under this condition, if one is only interested in checking the difference between these two sets of form factors, then all pieces needed to this end are in fact those that are used for demonstrating the cancellation of all spurious poles in these finite form factors.This property can be exploited to greatly improve the level of numerical accuracy of the comparison without increasing the computational cost, as the most complicated parts of master integrals required to obtain the form factors or amplitudes to 0 are not needed for this purpose.Eventually, we obtain agreement between the vector and axial form factors of the two-loop class-I diagrams in the full six-flavor theory (i.e., without the effective Higgs-gluon vertex) at the chosen kinematic point within the numerical uncertainty.In particular, taking advantage of the aforementioned insight, the suspicious 4-th vector and axial vector form factor agree with each other to the fourth significant digit, which we deem to be quite sufficient for our purpose 10 .
With the outcome of this critical check we conclude the following.If one naively computes the class-I diagrams in HEFT with a non-anticommuting γ 5 , some terms are missing, namely terms that involve the effective vertex qγ µ γ 5 qZ µ H.We have restored and determined these additional counterterms by enforcing the respective Ward identities as discussed in section 5.1 and 5.2.Computing analytically this particular missing piece directly by applying the heavy top-mass expansion procedure to these two-loop class-I diagrams and confirming the result deduced above would also be an interesting thing to do, which we however defer to future work.Our investigations conducted above show that the presence or absence of diagrams with the effective vertex qγ µ γ 5 qZ µ H in the heavy top-mass expansion of the class-I diagrams depends on the specific γ 5 prescription in use.When the axial vector current is regularized using a non-anticommuting γ 5 these are apparently needed.This observation further strengthens the common lore that one should be very careful with taking claims (or assuming conditions) which were established with an anticommuting γ 5 in a computation where a non-anticommuting γ 5 is employed instead in dimensional regularization.

Conclusions
In the first part of this article (sections 2 and 3), we first the vector current and non-anomalous axial vector current amplitudes b b → Zh proportional to the b-quark Yukawa coupling at two loops for a CP-even and CP-odd Higgs boson h = H, A, respectively.We showed that these polarized amplitudes can be obtained, when the b quark is taken to be kinematically massless, solely from the vector form factors of properly grouped classes of diagrams for ZH production.The computation of these form factors does not involve the axial vector current and hence γ 5 .Subsequently, we have compared with the previous results on b b → ZH of ref. [13] to two-loop order (where different projectors were used with a non-anticommuting γ 5 ).As expected, agreement of the axial part of the amplitudes can only be obtained at the level of properly defined finite remainders in four dimensions.Furthermore, the Ward identities for the QCD corrections to these b-quark Yukawa coupling-dependent contributions to b b → Zh, h = H, A, are derived and verified explicitly in section 3.
In addition, we determined the two-loop contributions to b b → Zh corresponding to diagrams that involve b-and t-quark triangles and the axial vector current in the heavy-top limit (cf.figure 3).For m t → ∞ the top-loop induced triangle diagrams are not vanishing.The explicit dependence on the renormalization scale µ R cancels in the total two-loop triangle contribution.
In the second part of this article (sections 4 and 5), we considered a class of top-Yukawa coupling dependent contributions to the amplitude of q q → ZH, namely the so-called class-I diagrams.Here the Higgs boson is radiated from the internal top-quark loop while the Z boson is emitted from the external light quark line.We computed these contributions to O(α 3 s ) in the heavy-top limit using the Higgs-gluon effective Lagrangian (HEFT).We obtained the analytical expressions of the UV renormalized vector form factors to twoloop order and verified their infrared poles by comparing to Catani's infrared factorization formula.For computing the axial vector form factors of the non-anomalous diagrams, an anticommuting γ 5 can be employed, which results in exactly the same UV renormalized expressions as their vector counterparts.
In an attempt to re-compute the QCD corrections to the same class-I Feynman diagrams in HEFT, but with a non-anticommuting γ 5 , a technical pitfall was noticed and discussed in detail in section 5.One may look at this issue from two perspectives.If one limits the scope to be within HEFT, then there are some additional local composite operators that one has to include when using a non-anticommuting γ 5 in the computation of class-I contributions to q q → ZH, as summarized in (5.17).On the other hand, if one looks at it from the point of view of the original full six-flavor theory, then our analysis in section 5.3 implies the following.If a non-anticommuting γ 5 is used in the axial vector current, then in the infinite top-mass limit certain heavy-mass expanded diagrams survive that are absent in a respective computation where an anticommuting γ 5 is used.Therefore, our results show that the presence or absence of certain heavy-mass expanded diagrams in the infinite-mass limit of a scattering amplitude with an axial vector current actually

Figure 1 .
Figure 1.Leading order Feynman diagrams that involve the bottom-Higgs Yukawa coupling.

Figure 2 .
Figure 2. Examples of diagrams where the Higgs boson or the Higgs and the Z boson are coupled to a closed quark loop.

Figure 3 .
Figure 3.The two-loop diagrams proportional to λ b that involve quark triangles.We denote the first four diagrams by box-triangles and the last two by pentagon-triangles.

Figure 4 .
Figure 4.The heavy-mass expansion of a representative pentagon-triangle diagram with k 1,2 denoting the two loop momenta.The solid thick lines represent the massive t-quark propagators.The dotted lines in the asymptotic irreducible graphs[25] on the left of "⊗" indicate external momenta which can all be set to zero if only the leading contribution in the limit m t → ∞ is to be kept.The graphs on the right of "⊗" denote the corresponding co-graphs in the heavy-mass expanded result.

Figure 5 .
Figure 5.The VVA triangle diagrams appearing as subgraphs in the heavy-mass expansion of figure 4.

Figure 6 .
Figure 6.The heavy-mass expansion of a representative box-triangle diagram with k 1,2 denoting the two loop momenta, drawn in a similar fashion as in figure 5, albeit, with a difference: the double-thick line represents the the amputated b-quark tree propagator, not to be confused with the solid thick triangle representing top quark loop.
Higgs boson A Now let us discuss the production of a CP-odd Higgs boson A in association with a Z boson by b b annihilation, b(p 1 ) + b(p 2 ) → Z(q 1 ) + A(q 2 ) .
and H, A denotes a CP-even and CP-odd Higgs boson, respectively.The kinetic terms of the Higgs bosons and of the Z boson are not listed in (3.1).Performing a continuous global U (1) transformation and applying the Noether theorem to the classical Lagrangian (3.1), b(x) → e iα b(x) and b(x) → e −iα b(x) ,leads to the well known conservation law for the vector current∂ µ J µ = 0 ,(3.3) which holds exactly also at the quantum level.Performing the continuous global U A (1) transformation b(x) → e iαγ 5 b(x) and b(x) → b(x)e iαγ 5 (3.4)

Figure 7 .
Figure 7. Example diagrams at tree level, one-loop and two-loop corresponding to the right hand side of the Ward identity (3.10).The solid arrow indicates the insertion of the external momentum −q 1 into the vertex.

Figure 8 .
Figure 8.The two-loop class-I diagrams of the top contributions to q q → ZH proportional to λ t .The thick solid lines denote the massive top quark.

Figure 13 .
Figure 13.Feynman diagrams contributing to the right-hand side of (5.14).The blob associated with H denotes the effective Hgg vertex while the blob associated with the auxiliary pseudoscalar A results from the local operator G G.