A Phase space element for VBF Higgs pair production process

# Next-to-next-to-leading order QCD corrections to light Higgs pair production via vector boson fusion in type-II two-Higgs-doublet model

## Abstract

We present the precision predictions on the pair production of light, -even Higgs in weak vector boson fusion (VBF) up to the QCD next-to-next-to-leading-order (NNLO) at hadron colliders within the -conserving type-II two-Higgs-doublet model (2HDM(II)) by adopting the structure function approach. We investigate the model parameter dependence, residual uncertainties from the factorization/renormalization scale, PDFs and on the integrated cross section at the QCD NNLO, and find that the NNLO QCD corrections can reduce the scale uncertainty significantly. By analyzing the kinematic distributions of final Higgs bosons, we can extract the -even Higgs resonance via channel as a means of probing the extension of the Standard Model (SM) Higgs sector.

PACS: 14.80.Ec, 12.38.Bx, 12.60.Fr

## 1 Introduction

Both ATLAS and CMS collaborations at the Large Hadron Collider (LHC) have discovered a neutral boson whose properties are compatible with the Standard Model (SM) Higgs boson [1, 2, 3]. The nature of this particle, including its properties and couplings, is currently being established [3, 4]. So the next important step most probably is the quest for the origin of electroweak symmetry breaking (EWSB). To achieve this goal, measurement of the Higgs self-interactions is necessary, which is the only way to reconstruct the Higgs potential, and determine whether the new particle is the SM Higgs boson or one of an enlarged Higgs sector of new physics. Thus, it is useful to explore the implication of the current Higgs search results on models beyond the SM.

One of the simplest extensions of the SM Higgs sector is the two-Higgs-doublet model (2HDM) [5, 6]. It predicts the existence of two neutral -even Higgs bosons ( and ), one neutral -odd Higgs (), and two charge Higgs bosons (). In addition to their masses, two additional parameters are introduced in the theory: the ratio of the vacuum expectation values (VEVs) of the two Higgs doublets , and the mixing angle between the two -even Higgs fields . There are many types of 2HDMs, each differing in the way that the two Higgs doublets couple to the fermions (for a comprehensive review, see [6]). In this paper, we only consider the 2HDM of type-II (2HDM(II)), which is designed to avoid flavor-changing couplings of the neutral Higgs bosons by one Higgs doublet coupling solely to up-type and the other to down-type fermions. And this model shares many of the features of the Higgs sector of the Minimal Supersymmetric Standard Model (MSSM).

To understand the Higgs self-interactions, the only accessible process is double Higgs production. At the LHC, the most important SM Higgs boson pair production channels have been systematically surveyed in Refs.[8]. The main processes are: (1) gluon-gluon fusion, , through heavy-quark loop, (2) vector boson fusion (VBF) , where vector bosons are radiated off quarks and fusion to Higgs pair, (3) top-quark pair associated Higgs boson pair production , and (4) double Higgs strahlung , where Higgs bosons are radiated off gauge bosons. In the SM, Higgs pair production via VBF has the second largest cross section and offers a clean experimental signature of two centrally produced Higgs bosons with two hard jets in the forward/backward rapidity region. Hence, it is meaningful to investigate the properties of trilinear Higgs self-interactions in this clean reaction. In this paper, we focus on the light -even Higgs pair production via VBF process within the 2HDM(II) to survey the properties of the trilinear Higgs self-couplings and appearing in the Higgs potential [7, 8]. In the previous research works, the VBF Higgs boson pair production process was surveyed in the 2HDM(II) at the QCD NLO [15, 16].

Due to the smallness of the QCD interference between the two inclusive final proton remnants, the VBF single/pair Higgs production at the leading order (LO) can be viewed as a double deep-inelastic scattering (DIS) process in a very good approximation, and the production rate can be computed by adopting the well-known structure function (SF) approach. Apart from the interference effect, the SF approach can still be exactly employed at the QCD next-to-leading order (NLO) [9, 10, 11]. Recently, the SF approach was used to calculate the VBF single/pair Higgs production at hadron colliders in the SM up to the QCD next-to-next-to-leading order (NNLO) [12, 13, 14]. In this paper, we will implement the SF approach to calculate the VBF Higgs pair production in the 2HDM(II) up to the QCD NNLO, and provide not only the total cross section, but also some kinematic distributions of the final Higgs bosons.

The paper is organized as follows. In Sec.2, we give a brief introduction to the 2HDM(II). The description of the SF approach and the strategy of the QCD NNLO calculation are presented in Sec.3. In Sec.4, we give the numerical results and focus on the theoretical uncertainty and some kinematic distributions. A short summary is given in Sec.5. Finally, we present the analytic expressions for the phase space element and matrix elements of the VBF Higgs pair production processes in Appendix.

## 2 Two-Higgs-Doublet Model of Type-II

The 2HDM contains two scalar doublets, and , with weak hypercharge . The most general Higgs potential with , and symmetries has the form as [5]

 V(Φ1,Φ2) = m211Φ†1Φ1+m222Φ†2Φ2+12λ1(Φ†1Φ1)2+12λ2(Φ†2Φ2)2 (2.1) +λ3(Φ†1Φ1)(Φ†2Φ2)+λ4(Φ†1Φ2)(Φ†2Φ1)+12λ5[(Φ†1Φ2)2+h.c.],

where , and are all real parameters, and transform under the discrete symmetry as

 Φ1→Φ1,        Φ2→−Φ2. (2.2)

After EWSB, the neutral components of and acquire VEVs and , respectively, which are determined by the vacuum conditions of

 Missing or unrecognized delimiter for \left Missing or unrecognized delimiter for \left (2.3)

and satisfy 1. We parameterize the two Higgs doublets as

 Φi=⎛⎝ϕ+i1√2(vi+Ri+iIi)⎞⎠,       (i=1,2). (2.4)

The Higgs mass matrices are diagonalized by performing the following rotation transformations:

 (H0h0)=R(α)(R1R2),    (G0A0)=R(β)(I1I2),    (G+H+)=R(β)(ϕ+1ϕ+2), (2.5)

where the rotation matrix is defined as

 R(θ)=(cosθsinθ−sinθcosθ), (2.6)

is the mixing angle between the two -even Higgs fields and , and . The fields and are Nambu-Goldstone bosons and their three degrees of freedom are got “eaten” by the longitudinal components of and bosons, and induce the masses of weak gauge bosons. Therefore, the 2HDM predicts five scalar particles: , , and . We may choose the following seven independent “physical” parameters as the inputs of the Higgs sector:

 mh0,  mH0,  mA0,  mH±,  sinα,  tanβ,  v. (2.7)

Then the quartic couplings can be expressed in terms of these physical parameters as

 λ1=1v2cos2β(m2h0sin2α+m2H0cos2α),    λ2=1v2sin2β(m2h0cos2α+m2H0sin2α), λ3=2m2H±v2+sin2αv2sin2β(m2H0−m2h0),    λ4=1v2(m2A0−2m2H±),    λ5=−1v2m2A0. (2.8)

The tree-level couplings of , and to the SM gauge bosons and fermions with respect to the corresponding couplings of the SM Higgs boson are presented in Table 1. It should be mentioned that the couplings of the -even Higgs bosons and have the same structures as those of the SM Higgs boson, while the Feynman rules for the interactions contain an additional factor since is a pseudoscalar. We can see from the table that when , the couplings of the light -even Higgs to gauge bosons and fermions approach the corresponding SM ones, and the couplings of the heavy -even Higgs to weak gauge bosons approach zero. Therefore, the -even Higgs decouples from the VBF process in the SM limit of . In this work we use as an input parameter of the Higgs sector instead of to manifest the effects on the Higgs couplings to gauge bosons involved in the VBF pair production, considering the fact that the Higgs couplings to weak gauge bosons are proportional to and .

## 3 Calculation Strategy

The SF approach is a very good approximation for studing the VBF processes at hadron colliders, which is accurate at a precision level well above the typical residual scale and parton distribution function (PDF) uncertainties [12]. This approximation is based on the absence or smallness of the QCD interference between the two inclusive final proton remnants. The Higgs boson pair production via VBF is a pure electroweak process at the LO, see Fig.1. There are two types of topological Feynman diagrams (- and -channel) contributing to the VBF Higgs pair production at parton level. The cross section is approximately contributed only by the squared - and -channel amplitudes, while their interference contribution is below . Therefore, the VBF Higgs pair production can be viewed as the double deep-inelastic scattering (DIS) of two (anti)quarks with two virtual weak vector bosons independently emitted from the hadronic initial states fusing into a Higgs boson pair [8]. The cross section can be calculated in terms of the charged-current and neutral-current hadronic structure functions by adopting the SF approach [9]. This method has been implemented to calculate the NNLO QCD corrections to the single Higgs production via VBF [12, 13]. Analogous to the case of the VBF single Higgs production, the nonfactorization contribution to the VBF Higgs pair production can also be safely neglected [13]. In this paper we adopt the SF approach to calculate the total inclusive cross section and differential distributions in the 2HDM(II) at the QCD NNLO for the VBF Higgs pair production .

The differential cross section for the VBF Higgs pair production can be expressed as [13]

 dσ = ∑(V1V2)12S2G2FM2V1M2V21(Q21+M2V1)21(Q22+M2V2)2WV1μν(x1,Q21)MμρV1V2M∗νσV1V2WV2ρσ(x2,Q22) (3.1) ×d3→PX1(2π)32EX1d3→PX2(2π)32EX2ds1ds2dPS2(k1,k2)(2π)4δ4(P1+P2−PX1−PX2−∑j=1,2kj),

where , is the Fermi constant, is the center-of-mass energy of the hadron collider, represents the phase space of the final two Higgs bosons, stands for the matrix element for the VBF subprocess , the physical scale is given by for and are the usual DIS variables, and are the invariant mass of the -th proton remnant. At the end of Eq.(3.1) there includes the four-body final state phase space element for the VBF Higgs pair production process, which is expressed explicitly in Appendix A.

The DIS hadronic tensor can be expressed in terms of the standard DIS structure functions as

 WVμν(xi,Q2i)=(−gμν+qi,μqi,νq2i)FV1(xi,Q2i)+^Pi,μ^Pi,νPi⋅qiFV2(xi,Q2i)+iϵμναβPαiqβi2Pi⋅qiFV3(xi,Q2i),  (V=Z,W±), (3.2)

where is the completely antisymmetric tensor and the momentum is defined as

 ^Pi,μ=Pi,μ−Pi⋅qiq2iqi,μ. (3.3)

Due to the conservation and the identity of the two final Higgs bosons, the matrix element for the process is the same as that for the process, i.e., . Here we depict the Feynman diagrams for the and processes in Fig.2 and Fig.3, respectively, and the explicit expressions for are presented in Appendix B. Then the squared DIS hadronic tensor in Eq.(3.1) can be written in the form as

 WV1μν(x1,Q21)MμρV1V2M∗νσV1V2WV2ρσ(x2,Q22)=3∑i,j=1CV1V2ijFV1i(x1,Q21)FV2j(x2,Q22), (3.4)

where can be automatically generated by using the Mathematica packages FeynArts [19] and FeynCalc [20].

Within the QCD factorization formalism, the structure functions can be expressed as convolutions of the PDFs in proton with the short-distance Wilson coefficient functions. We denote the gluon, quark and antiquark PDFs at the factorization scale by , and , respectively, where the subscript indicates the flavor of the (anti)quark. It is often convenient to write the DIS structure functions in terms of the gluon, and the following singlet and non-singlet quark distributions,

 qs=nf∑i=1(qi+¯qi),                                               (singlet), qvns=nf∑i=1(qi−¯qi),   q±ns,ij=(qi±¯qi)−(qj±¯qj),      (non-singlets). (3.5)

For the -exchange neutral current, the DIS structure functions can be written as follows [13]:

 FZi(x,Q2) = 2fi(x)∫10dy∫10dzδ(x−yz)nf∑j=1(v2j+a2j) ×[q+ns,j(y,μf)C+i,ns(z,Q,μr,μf)+qs(y,μf)Ci,q(z,Q,μr,μf)+g(y,μf)Ci,g(z,Q,μr,μf)], FZ3(x,Q2) = 2f3(x)∫10dy∫10dzδ(x−yz)nf∑j=12vjaj (3.6) ×[q−ns,j(y,μf)C−3,ns(z,Q,μr,μf)+qvns(y,μf)Cv3,ns(z,Q,μr,μf)],

where , , , , and the non-singlet quark densities are obtained from as

 q±ns,i=nf∑j=1q±ns,ij ,     (i=1,...,nf). (3.7)

The vector and axial-vector couplings of quark pair to boson used in Eqs.(3) are given by

 vi=I3i−2Qisin2θW,        ai=I3i, (3.8)

where and are the electric charge and weak isospin of the quark , respectively.

For the -exchange charged current, the DIS structure functions are expressed as follows:

 FW∓i(x,Q2) = fi(x)∫10dy∫10dzδ(x−yz)1nfnf∑j=1(v2j+a2j) ×[±δq−ns(y,μf)C−i,ns(z,Q,μr,μf)+qs(y,μf)Ci,q(z,Q,μr,μf)+g(y,μf)Ci,g(z,Q,μr,μf)], FW∓3(x,Q2) = f3(x)∫10dy∫10dzδ(x−yz)1nfnf∑j=12vjaj (3.9) ×[±δq+ns(y,μf)C+3,ns(z,Q,μr,μf)+qvns(y,μf)Cv3,ns(z,Q,μr,μf)],

where the non-singlet quark densities are defined in terms of as

 δq±ns=∑i∈up, j∈downq±ns,ij , (3.10)

and the vector and axial-vector couplings for charged current are given by

 vi=ai=1√2. (3.11)

We can see from Eqs.(3) and Eqs.(3) that the renormalization and factorization scales for quark densities in each proton (, , and ) enter in Eq.(3.1). The Wilson coefficient functions in Eq.(3) and Eq.(3) parameterize the hard partonic scattering process and can be perturbatively expanded in powers of . Up to the second order in , , and the perturbative expansion of these Wilson coefficient functions reads

 C±i,ns = δ(1−x)+as[c(1),±i,ns+LMP(0),±ns] (3.12) +a2s[c(2),±i,ns+LM(P(1),±ns+c(1),±i,ns⊗(P(0),±ns−β0))+L2M(12P(0),±ns⊗(P(0),±ns−β0)) +β0LR(c(1),±i,ns+LMP(0),±ns)],                                    (i=1,2,3),
 Ci,q = δ(1−x)+as[c(1)i,q+LMP(0)qq] (3.13) +a2s[c(2)i,q+LM(P(1)qq+c(1)i,q⊗(P(0)qq−β0)+c(1)i,g⊗P(0)gq) +L2M(12P(0)qq⊗(P(0)qq−β0)+12P(0)qg⊗P(0)gq) +β0LR(c(1)i,q+LMP(0)qq)],                         (i=1,2),
 Ci,g = as[c(1)i,g+LMP(0)qg] (3.14) +a2s[c(2)i,g+LM(P(1)qg+c(1)i,q⊗P(0)qg+c(1)i,g⊗(P(0)gg−β0)) +L2M(12P(0)qq⊗P(0)qg+12P(0)qg⊗(P(0)gg−β0)) +β0LR(c(1)i,g+LMP(0)qg)],                         (i=1,2),

where , , , is referred to the one-loop beta-function coefficient, and represents the standard Mellin convolution. It should be noted that 2

 P(0),±ns=P(0)qq,                     c(1),±i,ns=c(1)i,q,    (i=1,2,3). (3.15)

The two-loop order quark-quark splitting function and the quark singlet DIS coefficient functions are usually expressed as

 P(1)qq=P(1),+ns+P(1)ps,                     c(2)i,q=c(2),+i,ns+c(2)i,ps,    (i=1,2), (3.16)

where and are the pure-singlet contributions at the second order of . All the DIS coefficient functions , , , , , and the splitting functions , , , , , , used in Eqs.(3.12)-(3.14) are given in Refs.[21, 22, 23, 24, 25]. They can be easily evaluated in terms of harmonic polylogarithms [26]. In this paper we adopt the Fortran program Hplog [27] to implement numerical calculation of harmonic polylogarithms.

## 4 Numerical results and discussion

In this section we present the integrated cross sections and some kinematic distributions for the light -even Higgs pair production via VBF at , and proton-proton colliders up to the QCD NNLO by employing the SF approach. In our numerical calculations we use the following values for the electroweak parameters:

 MW=80.385 GeV,     MZ=91.1876 GeV,     GF=1.1663787×10−5 GeV−2. (4.1)

The Weinberg angle is fixed in the on-shell scheme as . We choose the 2HDM(II) input parameters at two benchmark points, B1 and B2, for demonstration and comparison, whose related parameters are listed in Table 2. The parameters at both the B1 and B2 points survive in the present theoretical and experimental constraints [17]. At the benchmark point B1 we have and there exists resonance effect in the VBF -pair production process. While at the benchmark point B2 there does not exist resonance effect, and the corresponding results should be the same with those in the SM case for the VBF production process. The width of can be calculated by using 2HDMC program [18], and at the benchmark point B1 we get the total decay width of boson being .

We adopt the MSTW2008 PDFs [28] in the convolutions of parton densities with Wilson coefficient functions. In the calculations of the -fusion contributions (see Fig.2), we take the -quark as a massless parton and the number of massless flavors in Eqs.(3). While in the evaluations of the -fusion process (see Fig.3), the initial -quark is not included since it would produce a top-quark in the final state. In the following analysis we take for simplicity and the typical central value of the renormalization/factorization scale is fixed by the corresponding vector-boson momentum transfer for 3, if there is no other statement. Furthermore, we put a lower bound of in order to keep in the perturbative regime, and the independence of the integrated cross section on this technical cut has been checked numerically.

### 4.1 Dependence on 2HDM(II) parameters

The integrated cross section for the VBF light, neutral -even Higgs boson pair production is related to the 2HDM(II) parameters, such as the two -even Higgs boson masses, ratio of the VEVs and the mixing angle between the two -even Higgs bosons. In this subsection we study the dependence of integrated cross section for the VBF production on the related model parameters at the LHC by adopting above event selection scheme.

In Fig.4(a) we depict the LO and NNLO QCD corrected integrated cross sections as functions of with the other related model parameters being the values at the benchmark point B1. We see from the figure that there is a steep increment at the position of due to the on-shell decay of , and resonance effect enhances the production rate obviously in the region of . It shows also that the QCD corrections up to NNLO always increase the LO cross section particularly for the large mass.

Fig.4(b) shows the dependence of the LO and QCD NNLO corrected integrated cross sections on the ratio of the VEVs . There we fix all the 2HDM(II) parameters are the values of the benchmark point B1 except , which varies from to . The figure demonstrates that both the LO and NNLO corrected total cross sections reach their minimal and maximal values at the positions about and , respectively. In the region of , both the LO and the NNLO QCD corrected cross sections exceed .

We plot Fig.4(c) to show the dependence of the LO and NNLO QCD corrected integrated cross sections on the parameter with the other related 2HDM(II) parameters being fixed at the benchmark point B1, i.e., and . It shows obviously that both the LO and the NNLO QCD corrected total cross sections reach their maxima at the position of , and then decrease with the increment of from to .

### 4.2 Theoretical uncertainties of integrated cross section

In order to make a precision comparison between the theoretical predictions and experimental measurements, we should assess thoroughly the theoretical uncertainties affecting the central predictions of the total cross sections. For some production processes at hadron colliders, such as process, the theoretical uncertainty mainly comes from the missing higher order corrections, PDFs and .

#### Scale uncertainty

The uncertainty due to missing higher order radiative corrections can be estimated by varying the factorization/renormalization scale around a central value that is taken close to the physical scale of the process. A conventional range of variation for the VBF process is

 14Q≤μ≤4Q, (4.2)

where the central value of and is the virtuality of the vector bosons which fuse into the Higgs boson pair. In Figs.5(a) and (b) we present the scale dependence of the LO, QCD NLO and NNLO corrected integrated cross sections for the VBF production at the LHC at the benchmark points B1 and B2, respectively. The central values of the integrated cross sections and the corresponding errors due to missing higher order radiative corrections are listed in Table 3. From Figs.5(a, b) and Table 3 we find that the scale uncertainties of integrated cross sections can be significantly reduced by including higher order radiative corrections. For the benchmark B1 (B2), the corresponding relative upper and lower scale relative uncertainties, defined as: , and