NNLO QCD$\oplus$QED corrections to Higgs production in bottom quark annihilation

TL;DR

This paper calculates the NNLO QED and mixed QCD-QED corrections to Higgs production via bottom quark annihilation at the LHC, improving precision in theoretical predictions for collider experiments.

Contribution

It provides the first systematic inclusion of NNLO QED and mixed QCD-QED corrections to this process, analyzing infrared structure and deriving relevant contributions.

Findings

QED and QCD-QED corrections significantly affect Higgs production rates.

Infrared poles are governed by universal anomalous dimensions.

Numerical impact of NNLO corrections is quantified at LHC energies.

Abstract

We present next-to-next-to leading order (NNLO) quantum electrodynamics (QED) corrections to the production of the Higgs boson in bottom quark annihilation at the Large Hadron Collider (LHC) in the five flavor scheme. We have systematically included the NNLO corrections resulting from the interference of quantum chromodynamics (QCD) and QED interactions. We have investigated the infrared (IR) structure of the bottom quark form factor up to two loop level in QED and in QCD$\times$QED using K+G equation. We find that the IR poles in the form factor are controlled by the universal cusp, collinear and soft anomalous dimensions. In addition, we derive the QED as well as QCD$\times$QED contributions to soft distribution function as well as to the ultraviolet renormalization constant of the bottom Yukawa coupling up to second order in strong coupling and fine structure constant. Finally, we report our findings on the numerical impact of the NNLO results from QED and QCD$\times$QED at the LHC energies taking into account the dominant NNLO QCD corrections.

Figures

Tables

Table 1. Table 1: Individual contributions in (pb) to various perturbative orders at S 𝑆 \sqrt{S} =14 TeV.

	${\Delta }^{0,0}$	${\Delta }^{1,0}$	${\Delta }^{0,1}$	${\Delta }^{2,0}$	${\Delta }^{1,1}$	${\Delta }^{0,2}$
LO00	1.0181
NLO10	1.1362	-0.1810
NLO01	1.2219		0.0030
NNLO20	1.1433	-0.1683		-0.1935
NNLO11	1.1542	-0.1699	0.0029		-0.0005
NNLO02	1.2422		0.0031			-4 ${10}^{-6}$

Table 2. Table 2: Individual contributions in (pb) to various perturbative orders at S 𝑆 \sqrt{S} =13 TeV.

	${\Delta }^{0,0}$	${\Delta }^{1,0}$	${\Delta }^{0,1}$	${\Delta }^{2,0}$	${\Delta }^{1,1}$	${\Delta }^{0,2}$
LO00	0.3911
NLO10	0.4588	0.1557
NLO01	0.4935		0.0003
NNLO20	0.4726	0.1614		0.0220
NNLO11	0.4771	0.1630	0.0003		1.5 ${10}^{-4}$
NNLO02	0.5135		0.0003			6 ${10}^{-6}$

Table 3. Table 3: 7-point scale variation at S 𝑆 \sqrt{S} =13 TeV.

($\frac{{\mu }_R}{m_h},\frac{{\mu }_F}{m_h}$)	(2,$\frac{1}{2}$)	(2,$\frac{1}{4}$)	(1,$\frac{1}{2}$)	(1,$\frac{1}{4}$)	(1,$\frac{1}{8}$)	($\frac{1}{2}$,$\frac{1}{4}$)	($\frac{1}{2}$,$\frac{1}{8}$)
NNLO20 (pb)	0.707	0.643	0.690	0.656	0.562	0.661	0.606
NNLO11 (pb)	0.759	0.602	0.780	0.641	0.445	0.682	0.498
NNLO02 (pb)	0.728	0.465	0.804	0.514	0.250	0.574	0.279

Table 4. Table 4: Result using different PDFs at S 𝑆 \sqrt{S} =14 TeV.

	MRST	NNPDF	CT14	PDF4LHC
NNLO20 (pb)	0.7805	0.7816	0.7574	0.8546
NNLO11 (pb)	0.9691	0.9867	0.9644	1.0625
NNLO02 (pb)	1.2020	1.2453	1.2288	1.3123

Table 5. Table 5: Result using different PDFs at S 𝑆 \sqrt{S} =13 TeV.

	MRST	NNPDF	CT14	PDF4LHC
NNLO20 (pb)	0.6610	0.6561	0.6398	0.7178
NNLO11 (pb)	0.6451	0.6406	0.6259	0.6996
NNLO02 (pb)	0.5252	0.5139	0.5030	0.5605

Equations

$${\mathcal{L}}_b=-{\lambda }_b\varphi (x){\bar{\psi }}_b(x){\psi }_b(x),$$

$$\frac{{\lambda }_t^{\mathrm{MSSM}}}{{\lambda }_b^{\mathrm{MSSM}}}=f_{\varphi }(\chi )\frac{m_t}{m_b}\frac{1}{\tan \beta },$$

$$f_{\varphi }(\chi )=\{\begin{array}{l}-\cot \chi \text{for}\varphi =h,\\ \tan \chi \text{for}\varphi =H,\\ \cot \beta \text{for}\varphi =A,\end{array}$$

$$\sigma (S,q^2)={\sigma }_0({\mu }_R^2)$$

$$\times {\Delta }_{cd}(s,q^2,{\mu }_F^2,{\mu }_R^2),$$

$${\Delta }_{cd}(z,q^2)=\sum _{i,j=0}^{\infty }{a}_s^i{a}_e^j{\Delta }_{cd}^{(i,j)}(z,q^2),$$

$$\displaystyle\frac{\hat{a}_{c}}{\big{(}\mu^{2}\big{)}^{\frac{\varepsilon}{2}}}S_{\varepsilon}=\frac{a_{c}(\mu_{R}^{2})}{\big{(}\mu_{R}^{2}\big{)}^{\frac{\varepsilon}{2}}}~{}Z_{a_{c}}\left(a_{s}(\mu_{R}^{2}),a_{e}(\mu_{R}^{2}),\varepsilon\right)\,,$$

$$Z_{a_s}$$

$$Z_{a_e}$$

$$\beta_{00}=\frac{11}{3}C_{A}-\frac{4}{3}n_{f}T_{F}\,,\quad\beta_{00}^{\prime}=-\frac{4}{3}\Big{(}N\sum_{q}e_{q}^{2}\Big{)}\,.$$

$${\mu }_R^2\frac{d}{d{\mu }_R^2}\ln Z_{\lambda }^b=\frac{\epsilon }{4}+{\gamma }_b^{(i,j)}(a_s({\mu }_R^2),a_e({\mu }_R^2))$$

$$Z_{{\lambda }_b}($$

$$\displaystyle+a_{s}^{2}\Big{\{}\frac{1}{\varepsilon^{2}}\Big{(}2\big{(}\gamma^{(1,0)}_{b}\big{)}^{2}+2\beta_{00}\gamma^{(1,0)}_{b}\Big{)}+\frac{1}{\varepsilon}\gamma^{(2,0)}_{b}\Big{\}}$$

$$\displaystyle+a_{e}^{2}\Big{\{}\frac{1}{\varepsilon^{2}}\Big{(}2\big{(}\gamma^{(0,1)}_{b}\big{)}^{2}+2\beta_{00}^{\prime}\gamma^{(0,1)}_{b}\Big{)}+\frac{1}{\varepsilon}\gamma^{(0,2)}_{b}\Big{\}}$$

$$\displaystyle+a_{s}a_{e}\Big{\{}\frac{1}{\varepsilon^{2}}\big{(}4\gamma^{(1,0)}_{b}\gamma^{(0,1)}_{b}\big{)}+\frac{1}{\varepsilon}\big{(}\gamma^{(1,1)}_{b}\big{)}\Big{\}}\,.$$

$$\displaystyle Q^{2}\frac{d}{dQ^{2}}\ln\hat{F}_{I}=\frac{1}{2}\Big{[}$$

$$\displaystyle+G_{I}\Big{(}\{\hat{a}_{c}\},\frac{Q^{2}}{\mu_{R}^{2}},\frac{\mu_{R}^{2}}{\mu^{2}},\varepsilon\Big{)}\Big{]}\,,$$

$$\displaystyle\mu_{R}^{2}\frac{d}{d\mu_{R}^{2}}K_{I}\Big{(}\{\hat{a}_{c}\},\frac{\mu_{R}^{2}}{\mu^{2}},\varepsilon\Big{)}=-A_{I}(\{a_{c}(\mu_{R}^{2})\})\,,$$

$$\displaystyle\mu_{R}^{2}\frac{d}{d\mu_{R}^{2}}G_{I}\Big{(}\{\hat{a}_{c}\},\frac{Q^{2}}{\mu_{R}^{2}},\frac{\mu_{R}^{2}}{\mu^{2}},\varepsilon\Big{)}=A_{I}(\{a_{c}(\mu_{R}^{2})\})\,,$$

$$A_I(\{a_c({\mu }_R^2)\})=\sum _{i,j}{a}_s^i({\mu }_R^2)a_e^j({\mu }_R^2)A_I^{(i,j)}$$

$$\displaystyle K_{I}(\mu_{R}^{2},\varepsilon)=\sum_{i,j}\hat{a}_{s}^{i}\hat{a}_{e}^{j}\Big{(}\frac{\mu_{R}^{2}}{\mu^{2}}\Big{)}^{(i+j)\frac{\varepsilon}{2}}S_{\varepsilon}^{(i+j)}K_{I}^{(i,j)}(\varepsilon)\,,$$

$$K_I^{(1,0)}$$

$$K_I^{(0,1)}$$

$$K_I^{(2,0)}$$

$$K_I^{(0,2)}$$

$$K_I^{(1,1)}$$

$$\displaystyle G_{I}\Big{(}\{\hat{a_{c}}\},\frac{Q^{2}}{\mu_{R}^{2}},\frac{\mu_{R}^{2}}{\mu^{2}},\varepsilon\Big{)}$$

$$+{\int }_{\frac{Q^2}{{\mu }_R^2}}^1\frac{d{\lambda }^2}{{\lambda }^2}{A}_I(\{a_c({\lambda }^2{\mu }_R^2)\})$$

$${\int }_{\frac{Q^2}{{\mu }_R^2}}^1\frac{d{\lambda }^2}{{\lambda }^2}$$

$$\displaystyle\times S_{\varepsilon}^{(i+j)}\Big{[}\Big{(}\frac{Q^{2}}{\mu_{R}^{2}}\Big{)}^{(i+j)\frac{\varepsilon}{2}}-1\Big{]}K^{(i,j)}(\varepsilon)\,.$$

$$G_I(\{a_c(Q^2)\},1,\epsilon )=\sum _{i,j}{a}_s^i(Q^2)a_e^j(Q^2)G_I^{(i,j)}(\epsilon ),$$

$$\displaystyle\ln\hat{F}_{I}=\sum_{i,j}\hat{a}_{s}^{i}\hat{a}_{e}^{j}\Big{(}\frac{Q^{2}}{\mu^{2}}\Big{)}^{(i+j)\frac{\varepsilon}{2}}S_{\varepsilon}^{(i+j)}\hat{\mathcal{L}}_{F_{I}}^{(i,j)}(\varepsilon)\,,$$

$${\hat{\mathcal{L}}}_{F_I}^{(1,0)}$$

$${\hat{\mathcal{L}}}_{F_I}^{(0,1)}$$

$${\hat{\mathcal{L}}}_{F_I}^{(2,0)}$$

$$\displaystyle-\beta_{00}G_{I}^{(1,0)}(\varepsilon)\Big{)}+\frac{1}{2\varepsilon}\Big{(}G_{I}^{(2,0)}(\varepsilon)\Big{)}\,.$$

$${\hat{\mathcal{L}}}_{F_I}^{(0,2)}$$

$$\displaystyle-\beta_{00}^{\prime}G_{I}^{(0,1)}(\varepsilon)\Big{)}+\frac{1}{2\varepsilon}\Big{(}G_{I}^{(0,2)}(\varepsilon)\Big{)}\,.$$

$${\hat{\mathcal{L}}}_{F_I}^{(1,1)}$$

$$G_I^{(i,j)}(\epsilon )$$

Full text

††preprint:  IMSc/2019/05/05, PSI-PR-19-12, TIF-UNIMI-2019-7

I Introduction

The discovery of the Standard Model (SM) Higgs boson by ATLAS Aad:2012tfa and CMS Chatrchyan:2012xdj collaborations at the Large Hadron Collider (LHC) has not only put the SM in a strong footing but also opened up a plethora of physics programs that can probe physics beyond the SM (BSM). Since the Higgs boson couples dominantly to heavy fermions and massive vector bosons, the corresponding observables are expected to be sensitive to new physics. In order to make definitive claims in the context of BSMs, it is henceforth extremely important to understand the Higgs sector of the SM. This is possible thanks to dedicated efforts from the LHC collaborations to measure the properties of the SM Higgs boson to unprecedented accuracy. Both ATLAS and CMS collaborations have already measured very precisely the partial width of the Higgs bosons both in the SM as well as in several BSM scenarios. This is the beginning of an era of precision physics at the LHC. These studies will be incomplete without the precise theoretical predictions in the SM as well as BSM.

At hadron colliders, where underlying scattering events are dominated by strong interaction, quantum effects are unavoidable. Quantum Chromodynamics (QCD), the theory of strong interaction, plays an important role at the LHC. Often, one finds that the leading order predictions from perturbative QCD are unreliable due to unphysical scales such as renormalization and factorization scales and also due to missing higher order radiative corrections. Radiative corrections from QCD are also large. Inclusion of such corrections improves the reliability of the predictions not only by making them more precise, but also by reducing the dependency on the unphysical scales.

At the LHC, the dominant production channel for the Higgs boson is gluon fusion through top quark loops Graudenz:1992pv . Owing to the complexities involved with the two loop massive Feynman integrals, an effective theory where top quark is integrated out, was proposed to obtain the first result at next-to-leading order (NLO) Dawson:1990zj for the Higgs boson production. Later on, in Djouadi:1991tka ; Spira:1995rr , the NLO corrections, taking into account the mass of the top quark, were shown to be very close to the prediction from the effective theory approach Dawson:1990zj . Thanks to the continued efforts beyond the NLO Harlander:2002wh ; Anastasiou:2002yz ; Ravindran:2003um ; Ravindran:2004mb ; Moch:2005ky ; Ravindran:2005vv ; Ravindran:2006cg , the most precise prediction till date, namely next-to-next-to-next-to leading order (N3LO) prediction Anastasiou:2015ema ; Mistlberger:2018etf for the inclusive production of the Higgs boson in the gluon fusion process, is now available (See Ravindran:2006bu ; Ahmed:2014uya ; Dulat:2018bfe for rapidity distributions). In addition, at NNLO accuracy, the tiny effects due to finite top quark mass have already been computed in Harlander:2009mq ; Pak:2009dg . Electroweak (EW) corrections Aglietti:2004nj ; Actis:2008ug and mixed QCD-EW corrections Anastasiou:2008tj are shown to improve the predictions.

The predictions from the perturbative QCD for the dominant production channel have reached the level of precision which now requires inclusion of the contributions from the sub-dominant channels. For example, one includes production channels such as vector boson fusion, associated production with a vector boson, bottom quark annihilation etc. In addition, the precise predictions deFlorian:2016spz taking into account radiative corrections from QCD and EW, are known for many of these processes.

Among these sub-dominant processes, production of the Higgs boson in bottom quark annihilation has been a topic of interest both in the SM as well as BSM contexts. In the SM, Yukawa couplings of the Higgs boson to the quarks and leptons are free parameters and precise determination of the couplings is possible at the LHC. These couplings are highly sensitive to scales of new physics as the mass ($m_{h}$) of the Higgs boson is close to the EW scale. Hence, both ATLAS ATLAS-CONF-2019-005 and CMS Sirunyan:2018koj collaborations have made dedicated efforts to measure them precisely. Among them, bottom Yukawa is one of the most sought one and it is a challenging task for experimentalists. Associated production of the Higgs boson with vector bosons or with top quarks and its subsequent decay to bottom quarks have been studied to achieve this. In addition, some interesting proposals can be found in Englert:2015dlp .

In the SM, bottom Yukawa coupling is less significant with respect to top Yukawa coupling while in the Minimal Supersymmetric SM (MSSM) Gunion:1992hs the coupling is proportional to $1/\cos\beta$ which can increase the cross section in some parametric region. The angle $\beta$ is related to the ratio, denoted by $\tan\beta$, of the vacuum expectation values of two Higgs doublets. The production of Higgs boson(s) in perturbative QCD is studied in four flavor and five flavor schemes Buza:1996wv ; Bierenbaum:2009mv ; Blumlein:2018jfm , called 4FS and 5FS, respectively. In the former, one assumes that proton sea does not contain bottom quark, and they are radiatively generated from gluons in the proton. These bottom quarks can annihilate to produce the Higgs boson. Their contributions are enhanced by logarithms of bottom quark mass spoiling the perturbation theory. Hence they need to be resummed to obtain reliable predictions. In the 5FS, one can avoid these logarithms by introducing non-zero bottom quark distributions in the proton. They are present due to pair production of bottom quarks from the gluons in the proton sea. Since, the leading order contribution in 5FS is two to one, while in 4FS, it is two to three, computations beyond the leading order are relatively easier in 5FS. In 4FS, only NLO QCD effects Dittmaier:2003ej ; Dawson:2003kb ; Wiesemann:2014ioa are known. On the other hand, in 5FS, NLO Dicus:1998hs ; Balazs:1998sb , NNLO Harlander:2003ai and the threshold effects at N3LO Ravindran:2006cg ; Ahmed:2014cha (see Ravindran:2006bu ; Ahmed:2014era for rapidity distributions) are known for sometime. Also, 5FS cross-section providing the dominant cross-section in a matched prediction Bonvini:2016fgf ; Forte:2015hba is very well known. In Ebert:2017uel , resummation of time-like logarithms in SCET framework has been performed. Recently, for the bottom quark annihilation, complete N3LO corrections Gehrmann:2014vha ; Ahmed:2014pka ; Duhr:2019kwi have become available. In addition, the resummation of threshold contributions H:2019dcl at N3LO+N3LL accuracy have also been included.

Unlike the dominant channel, gluon fusion to the Higgs boson, bottom quark annihilation has not received much attention in the context of EW corrections, presumably because it is already sub-dominant at the LHC. In this paper, we make the first attempt to include the QED corrections to the inclusive production to this channel. We expect that these corrections could be comparable to the fixed Duhr:2019kwi and resummed H:2019dcl results solely from third order in perturbative QCD.

In deFlorian:2018wcj , pure QED and mixed QCD$\times$QED contributions have been obtained for the Drell-Yan (DY) process through Abelianization deFlorian:2015ujt ; deFlorian:2016gvk at orders ${\cal O}(\alpha^{2})$ and ${\cal O}(\alpha\alpha_{s})$, respectively. In deFlorian:2018wcj , a suitable algorithm is obtained by studying the group theory structure of QCD and QED amplitudes that contribute to the partonic sub-processes of DY production. The algorithm contains a set of transformations on the color factors/Casimirs of SU(N) that transforms QCD results for the partonic sub-processes to the corresponding QED results. This way both pure QED as well as QCD$\times$QED contributions to inclusive production cross section for the Z boson in DY process have been obtained in deFlorian:2018wcj at NNLO level. Following this approach, we can in principle proceed to obtain pure QED and mixed QCD$\times$QED contributions to the bottom quark annihilation process from the QCD results. Although the QCD results Maltoni:2003pn ; Harlander:2003ai to NNLO are presented for $\text{N}=3$ of SU(N) and hence Abelianization can not be used, however, in Majhi:2010zg , resonant production of sleptons in a R-parity violating supersymmetric model was studied where radiative corrections from SU(N) gauge fields with $n_{f}$ fermions were included to NNLO level. Since, sleptons couple only to fermions in this model through Yukawa coupling, these NNLO corrections coincide with the results of Maltoni:2003pn ; Harlander:2003ai for $\text{N}=3$. Hence, we could use the results given in Majhi:2010zg and method of Abelianization to obtain pure QED as well as QCD$\times$ QED results for bottom quark annihilation to the Higgs boson. However, in order to scrutinize the very approach of Abelianization, we explicitly compute pure QED and QCD$\times$QED corrections to inclusive production of the Higgs boson in bottom quark annihilation up to NNLO level in U(1) and SU(N)$\times$U(1). In addition, we reproduce the same for the production of Z boson in DY process. The computation beyond the leading order involves evaluation of virtual and real emission processes. The contributions from them are sensitive to ultraviolet (UV) and infrared (IR) divergences. We compute them in dimensional regularization, hence divergences appear as poles in dimensional parameter $\varepsilon=d-4$, $d$ being the space-time dimension. The UV divergences are removed in ${\overline{\rm MS}}$ scheme. The IR divergences result from soft gluons and massless collinear partons. The former is called the soft divergence and later collinear divergence. While soft divergences cancel between virtual and real emission processes in the inclusive cross section, the collinear divergences are removed by mass factorization. We determine, both UV as well as mass factorization counter terms using factorization property of the inclusive cross section and obtain collinear finite contributions to the Higgs boson production in bottom quark annihilation and Z boson production in DY. We determine IR anomalous dimensions up to two-loop level in both QED and QCD$\times$QED. We find that they are process independent. Using the universal IR anomalous dimension and following Ahmed:2015qpa , we compute the renormalization constant for the Yukawa coupling in QED as well as in QCD$\times$QED from the form factors (FF) of Higgs bottom anti-bottom operator and vector current of DY process.

The paper is organized as follows. In section II, after discussing the theoretical frame work, we briefly describe in sub-section II.1 how we compute higher order QCD and QED radiative corrections to various partonic and photonic channels that contribute to the inclusive cross section. In sub-section II.2, we discuss the UV and IR structure of the form factors and cross sections using K+G equation and obtain the mass factorized cross sections. In the following sub-section, we discuss about the Abelianization procedure. The phenomenological impact of our theoretical predictions are presented in section III. Finally we summarize in section IV. The universal constants that appear in soft distribution function, FFs of vector current and bottom quarks and the mass factorized partonic and photonic cross sections are presented in the Appendix A, B and C, respectively.

II Theoretical framework

The Lagrangian that describes the interaction of the Higgs boson with the bottom quarks is described by the Yukawa interaction and is given by

[TABLE]

where $\lambda_{b}$ is the Yukawa coupling which, after the EW symmetry breaking, is found to be ${m_{b}}/{v}$. $\psi_{b}(x)$ and $m_{b}$ denote the bottom quark field and mass, respectively. $v$ is the vacuum expectation value (vev) of the Higgs field $\phi(x)$. In the SM, the Higgs boson production through bottom quark annihilation is sub-dominant compared to gluon fusion through top quark loop. One finds that the bottom Yukawa coupling is 35 times smaller than top quark Yukawa coupling and in addition, the bottom quark flux in the proton-proton collision is much smaller than the gluon flux. However, in the MSSM Nilles:1983ge , $\tan\beta$, the ratio of the vev s of Higgs doublets can increase the contributions resulting from the bottom quark annihilation channel. At LO,

[TABLE]

with

[TABLE]

where $h$ is the SM like light Higgs boson, $H$ and $A$ are the heavy and the pseudo-scalar Higgs bosons, respectively. The parameter $\chi$ is the mixing angle between weak and mass eigenstates of the neutral Higgs bosons $h$ and $H$. We set $m_{b}=0$ except in the Yukawa coupling  Aivazis:1993pi ; Collins:1998rz ; Kramer:2000hn as it is much smaller than the other energy scales in the process. The number of active flavors is taken to be $n_{f}=5$ and we work in the Feynman gauge.

The inclusive production of a colorless state in hadronic collisions is given by

[TABLE]

where $\sigma_{0}$ is the Born cross section and $f_{a}(x_{i},\mu_{F}^{2})$ are parton distribution functions (PDFs) for $a=q,\overline{q},g$ and photon distribution function (PHDF) if $a=\gamma$. The scaling variables $x_{i}$ is their momentum fractions. $\Delta_{cd}$ are the partonic sub-process contributions normalized by the Born cross section. The scales $\mu_{R}$ and $\mu_{F}$ are renormalization and factorization scales. $S$ and $s=x_{1}x_{2}S$ are hadronic and partonic center of mass energy, respectively. $q^{2}$ is the invariant mass of the final colorless state. $\Delta_{cd}$ can be expanded in powers of the QCD coupling constant $a_{s}=g_{s}^{2}(\mu_{R}^{2})/16\pi^{2}$ and QED coupling constant $a_{e}=e^{2}(\mu_{R}^{2})/16\pi^{2}$, $g_{s}$ and $e$ being the strong and electromagnetic coupling constants, respectively. That is, after suppressing $\mu_{R}$ and $\mu_{F}$ dependence,

[TABLE]

with $\Delta_{cd}^{(0,0)}=\delta(1-z)$ and $z=q^{2}/s$. In the following, we describe the methodology to compute $\Delta_{cd}^{(i,j)}$ up to second order in the couplings.

II.1 Methodology

In this section, we briefly describe how higher order perturbative corrections $\Delta_{cd}^{(i,j)}$ (Eq. 5) are computed. The details of computational procedure can be found in Ahmed:2016qhu . Beyond the leading order (LO), the partonic channels consists of one and two loop virtual sub processes, real-virtual and single and double-real emissions, some of which are presented in Fig. 1, Fig. 2 and Fig. 3. The black line with an arrow indicates the bottom quark, the wavy line the photon, the curly line the gluon and the Higgs boson is indicated by the dashed line.

Sub-processes involving virtual diagrams are sensitive to UV singularities. Due to the presence of massless gluons and photons, we encounter soft singularities in both virtual and real emission sub-processes. In addition, we encounter collinear singularities, as we treat all the quarks including the bottom quark massless. We use dimensional regularization to regulate all these singularities.

$\times$

\times$$\times

\times$$\times

\times$$\times$$\times

We have used the program QGRAF Nogueira:1991ex to generate virtual as well as real emission Feynman diagrams that contribute to the relevant sub-processes. An in-house FORM Vermaseren:2000nd code is used to perform all the symbolic manipulations, e.g. performing Dirac, SU(N) color and Lorentz algebra. Large number of loop integrals show up in the virtual diagrams. The integration-by-parts identities are used through a Mathematica based package, LiteRed Lee:2013mka to reduce them to a minimum set of master integrals. For those processes that involve pure real emissions with or without virtual diagrams, we use the method of reverse unitarity Anastasiou:2002yz that allows one to use IBP identities to reduce the resulting phase-space integrals to a set of few master integrals, the later can be found in Anastasiou:2012kq . Finally we obtain contributions to each sub-process, containing UV and IR singularities as poles in $\varepsilon=d-4$.

In the next section, we study the UV and IR structure of FF and the soft distribution function that contribute to NNLO level in QCD, QED and QCD$\times$QED. In order to explore the IR structure, we study the production of a $Z$ boson in hadron colliders, namely DY process to the same accuracy in QCD, QED and QCD$\times$QED. In particular, we focus our attention to the FF and the soft distribution function that contribute to the inclusive DY production cross section. Following Ravindran:2006cg ; Ahmed:2015qpa ; Ahmed:2014cla , we demonstrate the factorization of IR singularities in both the FFs and show how to extract the process independent cusp ($A$), collinear ($B$) and soft ($f$) anomalous dimensions from them. Using the FF of the bottom quark and the process independent soft distribution function we can extract the UV anomalous dimension of the Yukawa coupling $\lambda_{b}$ up to two loop level in QCD, QED and QCD$\times$QED. Finally, we demonstrate the factorization of collinear singularities and the mass factorization leading to IR finite partonic contributions to inclusive hadronic cross sections for both the Higgs boson and DY productions up to NNLO level in QCD, QED and QCD$\times$QED.

II.2 UV and IR structures in QED and QCD$\times$QED

Having computed all the partonic channels that contribute to the hadronic cross sections in QED and QCD$\times$QED, we use them to study the UV and IR structure of the FFs and soft gluon/photon emissions in the Higgs boson and DY productions. For the former, we have used Sudakov K+G equation and for the later, following Ravindran:2005vv ; Ravindran:2006cg we exploited the universal structure of soft distribution function resulting from the soft gluon/photon emissions.

In order to remove the UV divergences that result from virtual sub-processes, we use the renormalization constants $Z_{a_{c}},c=s,e$ for the QCD and QED coupling constants, respectively and $Z_{\lambda_{b}}$ for the Yukawa coupling. The Yukawa coupling $Z_{\lambda_{b}}$ receives contributions from both QCD and QED. $Z_{a_{s}}$ and $Z_{a_{e}}$ relate the bare couplings $\hat{a}_{s}=\hat{g}_{s}^{2}/16\pi^{2}$ of QCD and $\hat{a}_{e}=\hat{e}^{2}/16\pi^{2}$ of QED to the renormalized ones $a_{s}(\mu_{R}^{2})$ and $a_{e}(\mu_{R}^{2})$, respectively, at the renormalization scale $\mu_{R}$ in the following way,

[TABLE]

where $a_{c}=\{a_{s},a_{e}\}$. Here, $S_{\varepsilon}\equiv\exp[\left(\gamma_{E}-\ln 4\pi\right)\frac{\varepsilon}{2}]$ is the phase-space factor in $d$-dimensions, $\gamma_{E}=0.5772...$ is the Euler-Mascheroni constant and $\mu$ is an arbitrary mass scale introduced to make $\hat{a}_{s}$ and $\hat{a}_{e}$ dimensionless in $d$-dimensions. The renormalization constant $Z_{a_{c}}$ up to two-loops are given by

[TABLE]

where $\beta_{ij}$ and $\beta^{\prime}_{ij}$ are QCD and QED beta functions, respectively. In the present case, only one loop $\beta$ i.e. $\beta_{00}$ and $\beta_{00}^{\prime}$Cieri:2018sfk appear. They are given by

[TABLE]

Here, $C_{A}=N$ is the adjoint Casimir of $SU(N)$. We also denote the fundamental Casimir $C_{F}=(N^{2}-1)/2N$ for later use. $n_{f}$ is the number of active quark flavors and $e_{q}$ refers to electric charge for quark $q$. The renormalization constant $Z_{\lambda}^{b}(a_{s},a_{e})$ satisfies the renormalization group equation:

[TABLE]

whose solution in terms of the anomalous dimensions $\gamma^{(i,j)}_{b}$ up to two loops is found to be

[TABLE]

Note that while the UV singularities factorize through $Z_{\lambda_{b}}$, singularities from QCD and QED mix from two loop onward. For QCD, $\gamma^{(i,0)}_{b}$ is known to four loops Vermaseren:1997fq . In this paper, using the universal IR structure of the amplitudes and cross sections in QED, we determine $\gamma^{(i,j)}_{b}$ up to two loops in QED i.e. for $(i,j)=(0,1),(0,2)$ and in QCD$\times$QED i.e. for $(i,j)=(1,1)$.

We begin with the bare form factors $\hat{F}_{I}(\hat{a}_{s},\hat{a}_{e},Q^{2},\mu^{2})$, $I=q,b$ where $q$ denotes the DY process and $b$ denotes the Higgs boson production in bottom quark annihilation. Note that, these FFs are computed in the perturbative framework where both QCD as well as QED interactions are taken into account simultaneously and hence they depend on both QCD and QED coupling constants. In addition, we find that the UV renormalized FFs demonstrate the factorization of IR singularities. Using gauge and renormalization group invariance, we propose Sudakov integro-differential equation for these FFs, analogous to the QCD one. In dimensional regularization, they take the following form

[TABLE]

where $\{a_{c}\}=\{a_{s},a_{e}\}$ and $Q^{2}=-q^{2}$ is the invariant mass of the final state particle (di-lepton pair in the case of DY and single Higgs boson for the case of Higgs production). Explicit computation of the form factors shows that IR singularities, resulting from QCD and QED interactions not only factorize but also mix beyond one loop level. In other words, if we factorize IR singularities from the FFs, the resulting IR singular function can not be written as a product of pure QCD and pure QED functions. More specifically, there will be terms proportional to $a_{s}^{i}a_{e}^{j}$, where $i,j>0$, which will not allow factorization of QCD and QED ones. Hence, $K_{I}$ will have IR poles in $\varepsilon$ from pure QED and pure QCD in every order in perturbation theory and in addition, from QCD$\times$QED starting from ${\cal O}(a_{s}a_{e})$. On the other hand, overall factorization of IR singularities implies that the constants $K_{I}$ contain the IR singularities from QCD, QED and QCD$\times$QED, while the $G_{I}$s will have IR finite contributions. Since, the IR singularities of FFs have dipole structure, $K_{I}$ will be independent of $q^{2}$ while $G_{I}$s will be finite in $\varepsilon\rightarrow 0$ and the later contain only logarithms in $q^{2}$. Note that $\hat{F}_{I}$ are renormalization group (RG) invariant so does the sum $K_{I}+G_{I}$. Thus, the RG invariance of $\hat{F}_{I}$ implies

[TABLE]

where $A_{I}$ are the cusp anomalous dimensions. The solutions to the above RG equations for $K_{I}$ can be obtained by expanding the cusp anomalous dimensions ($A_{I}$) in powers of renormalized coupling constants $a_{s}(\mu_{R}^{2})$ and $a_{e}(\mu_{R}^{2})$ as

[TABLE]

and $K_{I}$ as

[TABLE]

where $A^{(i,0)}$ and $A^{(0,i)}$ result from pure QCD and pure QED interactions and $A^{(i,j)},i,j>0$ from QCD$\times$QED. Using RG equations for the couplings $a_{s}$ and $a_{e}$, the perturbative solutions to Eq.(12) are found to be,

[TABLE]

Unlike $K_{I}$, $G_{I}$ do not contain any IR singularities but depend only on $Q^{2}$ and hence we expand them as

[TABLE]

where the first term is the boundary condition on each $G_{I}$ at $\mu_{R}^{2}=Q^{2}$. Expanding $A_{I}$ in powers of $a_{s}$ and $a_{e}$ and using RG equations for QCD and QED couplings, we obtain

[TABLE]

Expanding the finite function $G_{I}(a_{s}(Q^{2}),a_{e}(Q^{2}),1,\varepsilon)$ as,

[TABLE]

substituting the solutions of $K_{I}$ and $G_{I}$ in Eq.(II.2) and performing the integration over $Q^{2}$ we get

[TABLE]

where,

[TABLE]

Following Ravindran:2004mb , we expand $G_{I}^{(i,j)}(\varepsilon)$ around $\varepsilon=0$ in terms of collinear ($B_{I}^{(i,j)}$), soft ($f_{I}^{(i,j)}$) and UV ($\gamma_{I}^{(i,j)}$) anomalous dimensions as

[TABLE]

with

[TABLE]

The form factors $\hat{F}_{I}$ that we computed in this paper in QCD, QED and QCD$\times$QED up to two loop level can be used to extract the cusp anomalous dimensions ($A_{I}^{(i,j)}$) by comparing them against Eq.(II.2). We find $A_{I}^{(i,j)}$ up to two loops as

[TABLE]

Unlike $A_{I}^{(i,j)}$, the other anomalous dimensions $B_{I}^{(i,j)}$, $f_{I}^{(i,j)}$ and $\gamma^{(i,j)}_{I}$ ($\gamma^{(i,j)}_{q}$ is zero) can not be disentangled either from $\hat{F}_{q}$ or $\hat{F}_{b}$ alone. In order to disentangle $B_{I}^{(i,j)}$ and $f_{I}^{(i,j)}$, we study the partonic cross sections resulting from soft gluon and soft photon emissions as they are only sensitive to $f_{I}^{(i,j)}$.

To obtain the process independent part of soft gluon/photon contributions in the real emission sub-processes, we follow the method described in Ravindran:2005vv ; Ravindran:2006cg , where the soft distribution function for the inclusive cross section for producing a colorless state was obtained from the form factors and partonic sub-process cross sections involving real emissions of gluons. The soft distribution functions denoted by $\Phi_{J}$, are governed by cusp ($A_{J}$) and soft anomalous dimensions $f_{J}$, where $J=q,b,g$. It is also known that the identity $\Phi_{b}=\Phi_{q}=C_{F}/C_{A}\Phi_{g}$ holds up to three loop level Ravindran:2005vv ; Ravindran:2006cg ; Ahmed:2014cla . We can use the partonic sub-processes of either DY process or the Higgs boson production in bottom quark annihilation namely $\hat{\sigma}_{q\overline{q}}$ or $\hat{\sigma}_{b\overline{b}}$ normalized by the square of the bare form factor $\hat{F}_{q}$ or $\hat{F}_{b}$ to obtain $\Phi_{I}$. In general $\Phi_{I}$, which is function of the scaling variable $z=q^{2}/s$, is defined as,

[TABLE]

with $Z_{q}=1$ and $Z_{b}=Z_{\lambda_{b}}$ being the overall renormalization constant. The symbol $\mathcal{C}$ refers to “ordered exponential” which has the following expansion:

[TABLE]

Here $\otimes$ is the Mellin convolution and $f(z)$ is a distribution of the kind $\delta(1-z)$ and $\mathcal{D}_{i}$. The plus distribution $\mathcal{D}_{i}$ is defined as,

[TABLE]

We can compute the UV finite $\hat{\sigma}_{I\overline{I}}$ every order in renormalized perturbation theory. Since, we have not determined $Z_{\lambda_{b}}$, we can only compute the unrenormalized partonic cross section $\tilde{\sigma}_{I\overline{I}}=\hat{\sigma}_{I\overline{I}}/Z_{I}^{2}$. From the explicit results for $\tilde{\sigma}_{I\overline{I}}$ and the form factors $\hat{F}_{I}$, using Eq. (25) we obtain $\Phi_{I}$ up to second order in $a_{s}$, $a_{e}$ and $a_{s}a_{e}$. We find $\Phi_{q}=\Phi_{b}$ up to second order in the couplings demonstrating the universality. In Ravindran:2005vv ; Ravindran:2006cg , it was shown that the soft distribution function $\Phi_{I}$ satisfies Sukakov K+G equation analogous to the form factor $\hat{F}_{I}$ due to similar IR structures that both of them have, order by order in perturbation theory. That is, $\Phi_{I}$ satisfies

[TABLE]

where, the IR singularities are contained in $\overline{K}$ and the finite part in $\overline{G}$. RG invariance of $\Phi_{I}$ implies

[TABLE]

Note that, the same anomalous dimensions govern the evolution of both $\overline{K}_{I}$ and $\overline{G}_{I}$. This ensures that the soft distribution function contains right soft singularities to cancel those from the form factor leaving bare partonic cross section to contain only initial state collinear singularities. The later will be removed by mass factorization by appropriate Altarelli-Parisi kernels. Expanding $\overline{K}_{I}(\{a_{c}\})$ and $\overline{G}_{I}(\{a_{c}(q^{2})\},1,\varepsilon,z)$ in powers of $\{a_{c}\}$ as has been done for $K_{I}(\{a_{c}\})$ and $G_{I}(\{a_{c}\})$ (see Eqs.(15,19)), with the replacements of $K_{I}^{(i,j)}$ by $\overline{K}_{I}^{(i,j)}$ and

[TABLE]

the solution to Eq.(II.2) is found to be

[TABLE]

where,

[TABLE]

$\overline{G}_{I}^{(i,j)}(\varepsilon)$ is related to finite function $\overline{G}_{I}(\{a_{c}(q^{2})\},1,\varepsilon,z)$ defined in Eq.(30) through the distributions $\delta(1-z)$ and $\mathcal{D}_{j}$. Thus expanding $\overline{G}^{(i,j)}(\varepsilon)$ in terms of the $a_{s}(q^{2}(1-z)^{2})$ and $a_{e}(q^{2}(1-z)^{2})$ we write,

[TABLE]

where $q_{z}^{2}=q^{2}(1-z)^{2}$. Following, Ravindran:2005vv ; Ravindran:2006cg , the IR finite $\overline{\mathcal{G}}^{(i,j)}_{I}(\varepsilon)$ can be expanded as

[TABLE]

where, for up to two loops

[TABLE]

Comparing the soft distribution functions $\Phi_{I}$, $I=q,b$, obtained from the explicit computation up to second order in coupling constants against the formal solution given in Eqs.(II.2), we can obtain $A_{I}^{(i,j)}$ and $f_{I}^{(i,j)}$ for $(i,j)=(1,0),(0,1),(1,1),(2,0),(0,2)$. Finally, we obtain $f_{I}^{(1,0)}=f_{I}^{(0,1)}=f_{I}^{(1,1)}=0$ and

[TABLE]

Now that we have $f_{I}^{(i,j)}$, it is now straightforward to obtain $B_{q}^{(i,j)}$ in Eq.(22) from the explicit results on $G^{(i,j)}_{q}$ as $\gamma^{(i,j)}_{q}=0$ for DY. This way we obtain,

[TABLE]

Assuming $B_{b}^{(i,j)}=B_{q}^{(i,j)}$, we determine the UV anomalous dimension, $\gamma_{b}^{(i,j)}$ from $G^{(i,j)}_{b}$ (Eq.(22)) which is known to second order. They are found to be

[TABLE]

Alternatively, assuming $B_{b}^{(i,j)}=B_{q}^{(i,j)}$ and $f_{b}^{(i,j)}=f_{q}^{(i,j)}$, we can determine $\gamma^{(i,j)}_{b}$ by comparing the difference $G_{b}^{(i,j)}-G_{q}^{(i,j)}$ obtained using DY and Higgs boson form factors $\hat{F}_{q}$ and $\hat{F}_{b}$ at $\varepsilon=0$ against the formal decomposition of $G_{I}^{(i,j)}$ given in Eqs.(22). Substituting the above UV anomalous dimensions in Eq.(II.2), we obtain $Z_{\lambda_{b}}$ to second order in the couplings.

Using the renormalization constants $Z_{a_{s}}$, $Z_{a_{e}}$ and $Z_{\lambda_{b}}$ for the coupling constants $\alpha_{s}$, $\alpha_{e}$ and the Yukawa coupling $\lambda_{b}$, we obtain UV finite partonic cross sections. The soft and collinear singularities arising from gluons/photons/fermions in the virtual sub-processes cancel against those from the real sub-processes when all the degenerate states are summed up, thanks to the KLN theorem Kinoshita:1962ur ; Lee:1964is . What remains at the end, is the initial state collinear singularity, which can be removed by mass factorization. Collinear factorization allows us to determine the mass factorization kernels $\Gamma_{qq}$ and $\Gamma_{qg}$ up to two-loop level for $U(1)$ and $SU(N)$ $\times$ $U(1)$ cases. Since $\Gamma_{qq}$ and $\Gamma_{qg}$ are governed by the splitting functions $P_{qq}$ and $P_{qg}$, we extract them to second order in couplings. In deFlorian:2015ujt , these splitting functions up to NNLO level, both in QED and QCD$\times$QED, were obtained using the Abelianization procedure. The splitting functions that we have obtained by demanding finite-ness of the mass factorised cross section, agree with those in deFlorian:2015ujt . The mass factorized partonic cross section for each partonic sub-process up to NNLO in QED and in QCD$\times$QED are presented in the Appendix along with the known NNLO QCD results Harlander:2003ai . In the next section, we use them to study their numerical impact at the LHC energies.

II.3 Abelianization procedure

In deFlorian:2018wcj , QCD$\times$QED corrections to the DY process were obtained by studying the $SU(N)$ color factors in Feynman diagrams that contribute to QCD corrections. This led to an algorithm namely Abelianization procedure which provides a set of rules that transform QCD results into pure QED and mixed QCD$\times$QED results. Unlike in deFlorian:2018wcj , without resorting to Abelianization rules, we have performed explicit calculation to obtain the contributions resulting from all the partonic and photonic channels taking into account both UV and mass factorization counter terms. Using these results at NNLO in QCD, QCD$\times$QED and in QED, we find a set of rules that can relate QCD and QED results. Note that if there is a gluon in the initial state, averaging over its color factor gives a factor $\frac{1}{N^{2}-1}$. This is absent for the processes where photon is present instead of gluon in the initial state. Also, for pure QCD or QED, the gluons or photons are degenerate and hence one needs to account for a factor of 2. Keeping these in mind, we arrive at a set of relations among QCD and QED results. We have listed them in the following tables for various scattering channels. They are found to be consistent with the procedure used in deFlorian:2018wcj .

Rule 1 : quark-quark initiated cases

[TABLE]

${}^{*}e_{q}^{2}=e_{b}^{2}$ when both initial quarks are bottom quarks.

Rule 2 : *quark-gluon initiated cases

*(After multiplying $2C_{A}C_{F}$ for the initial state gluon)

[TABLE]

Rule 3 : *gluon-gluon initiated cases

*(After multiplying $2C_{A}C_{F}$ for each initial state gluon)

[TABLE]

III Results and phenomenology

In this section, we study the numerical impact of pure QED and mixed QCD$\times$QED corrections over the dominant QCD corrections up to NNLO level to the production of the Higgs boson in bottom quark annihilation at the LHC, mainly for the center of mass (CM) energy of $\sqrt{S}=13$ TeV. Since we include QED effects, we need PHDF inside the proton in addition to the standard PDFs. For this purpose, we use NNPDF 3.1 LUXqed set Bertone:2017bme , MRST Martin:2004dh , CT14 Schmidt:2015zda and PDF4LHC17. The PDFs, PHDFs and the strong coupling constant $a_{s}$ can be obtained, using the LHAPDF-6 Buckley:2014ana interface. We have used the following input parameters for the masses and the couplings:

[TABLE]

Both $a_{s}(\mu_{R})$ and $m_{b}(\mu_{R})$ are evolved using appropriate QCD $\beta$-function coefficients and quark mass anomalous dimensions respectively. However, we have considered fixed $\alpha_{e}=4\pi a_{e}$ throughout the computation.

The Higgs boson production cross section from bottom quark annihilation at the present energy of LHC is not substantial. However, for the high luminosity LHC, measuring them at higher center of mass energy (CM) would give larger contributions and it will improve the precision. Hence, we have first studied how the cross section varies with the CM of LHC. In Fig. 4, we plot the inclusive production cross sections at various orders in perturbative QCD and QED for the range of CM energies between $\sqrt{S}=6$ to $22$ TeV. In the inset, the index ‘ij’ indicates that QCD at ‘i’-$th$ order and QED at ‘j’-$th$ order in perturbative theory are included (e.g. ‘NNLO 11’ indicates NNLO mixed QCD$\times$QED). In Fig. 4, we have used NNPDF31_lo_as_0118, NNPDF31_nlo_as_0118_luxqed and NNPDF31_nnlo_as_0118_luxqed for LO, NLO and NNLO, respectively. The renormalization ($\mu_{R}$) and factorization ($\mu_{F}$) scales are kept fixed at $m_{h}$ and $m_{h}/4$, respectively. We note that in Fig. 4, the pure QED contributions are large. This is due to the fact that we consider leading order QCD running of Yukawa coupling which gives larger Born contribution compared to pure QCD. In order to understand this in more detail, we study the impact of different contributions to the cross sections resulting from QCD, QED and mixed QCD$\times$QED at various orders in perturbation theory which we have tabulated in Table 1, for $\sqrt{S}=14$ TeV and for the scale choice $\mu_{R}=\mu_{F}=m_{h}$.

The $\Delta^{i,j}$ indicates sole $i$-th order QCD and $j$-th order QED corrections to the total contribution. For example, NNLO11 means $\Delta^{0,0}+\Delta^{1,0}+\Delta^{0,1}+\Delta^{1,1}$.

In Table 2, a similar study has been performed for $\sqrt{S}=13$ TeV and the scales $\mu_{R}=m_{h}\,,~{}\mu_{F}=m_{h}/4$.

Fixed order predictions depend on the renormalization ($\mu_{R}$) and factorization ($\mu_{F}$) scales. The uncertainty resulting from the choice of the scales quantify the missing higher order contributions. Hence, we have studied their dependence by varying them independently around a central scale.

Fig. 5 shows the dependence of the cross section on the renormalization scale ($\mu_{R}$) for the fixed choice of the factorization scale $\mu_{F}=m_{h}/4$. It clearly demonstrates the importance of higher order corrections as the $\mu_{R}$ variation is much more stable at NNLO20 compared to the lower orders.

In Fig. 6, we present the dependence on the factorization scale ($\mu_{F}$) keeping the renormalization scale ($\mu_{R}$) fixed at $m_{h}$. Similar to the $\mu_{R}$ variation, $\mu_{F}$ variation improves after adding higher order corrections. To illustrate their dependence when both the scales are changed simultaneously, we present the cross section by performing 7-point scale variation and the results are listed in Table 3. We have used NNPDF31_nnlo_as_0118_luxqed for this study.

The perturbative predictions also depend on the choice of PDFs and PHDFs. There are several groups which fit them and are widely used in the literature for the phenomenological studies. In order to estimate the uncertainty resulting from the choice of PDFs and PHDFs, in Table 4, we present the NNLO results from various PDF sets, for $\sqrt{S}$=14 TeV and $\mu_{R}=\mu_{F}=m_{h}$.

In Table 5, we repeat the study for $\sqrt{S}$=13 TeV and $\mu_{R}=m_{h}$ and $\mu_{F}=m_{h}/4$.

We have also studied the uncertainties resulting from the choice of PDF set Buckley:2014ana . Using NNPDF31, in Fig 7, we plot the variation of the cross section with respect to different choices of PDF and PHDF templates keeping the central set as the reference. The thick line is obtained using the central set. The shaded region resulting from other sets quantifies the uncertainty.

IV Discussion and Conclusion

Precision studies is one of the prime areas at the LHC. Measuring the parameters of the SM to unprecedented accuracy can help us to improve our understanding of the dynamics that governs the particle interactions at high energies. This is possible only if the accuracy of theoretical predictions is comparable to that of the measurements. Thanks to the on-going efforts from experimentalists and theorists, there are stringent constraints on various physics scenarios in the pursuit of searching for the physics beyond the SM. The efforts to compute the observables that are related to top quarks and Higgs bosons have been going on for a while as these observables are sensitive to high scale physics. Since the dominant contributions to these processes are known to unprecedented accuracy, inclusion of sub-dominant contributions along with radiative corrections is essential for any consistent study. In this context, the present article explores the possibility of including EW corrections to Higgs boson production in bottom quark annihilation which is sub-dominant. Note that, this is known to third order in QCD Duhr:2019kwi . While this is a sub-dominant process at the LHC, in certain BSM contexts, the rates are significantly appreciable leading to interesting phenomenological studies. Since, the computation of full EW corrections is more involved, as a first step towards this, we compute all the QED corrections, in particular, to the inclusive Higgs boson production in bottom quark annihilation up to second order in QED coupling constant $a_{e}$, taking into account the non-factorizable or mixed QCD$\times$QED effects through $a_{s}a_{e}$ corrections. The computation involves dealing with QED soft and collinear singularities resulting from photons and the massless partons along with the corresponding QCD ones. Understanding the structure of these QED IR singularities in the presence of QCD ones, is a challenging task. We have systematically investigated both QCD and QED IR singularities up to second order in their couplings taking into account the interference effects. We use Sudakov K+G equation to understand the IR structure in terms of cusp, colliner and soft anomalous dimensions. We demonstrate that the IR singularities from QCD, QED and QCD$\times$QED interactions factorize both at the FF, as well as at the cross section level. While the IR singularities factorize as a whole, the IR singularities from QCD do not factorize from that of QED leading to mixed/non-factorizable QCD$\times$QED IR singularities. In addition, by computing the real emission processes in the limit when the photons/gluons become soft, we have studied the structure of soft distribution function. While the later demonstrates the universal structure analogous to QCD one, we find that it contains soft terms from mixed QCD$\times$QED that do not factorize either as a product of those from QCD and QED separately. Using the universal IR structure of the observable, we have determined the mass anomalous dimension of the bottom quark and hence the renormalization constant for the bottom Yukawa. We also discussed the relation between the results from pure QED and pure QCD as well as between QCD$\times$QED and QCD through Abelianization. We have determined a set of rules that relate them and they are found to be consistent with those observed in the context of DY deFlorian:2018wcj . Having obtained the complete NNLO results from QED and QCD$\times$QED, we have systematically included them in the NNLO QCD study to understand their impact at the LHC energy. We find that the corrections are mild as expected. However, we show that the higher order corrections from QED and QCD$\times$QED improve the reliability of the predictions.

Acknowledgment

A.A.H, A.C, P.M and V.R would like to acknowledge the support of the CNRS LIA (Laboratoire International Associé) THEP (Theoretical High Energy Physics) and the INFRE-HEPNET (IndoFrench Network on High Energy Physics) of CEFIPRA/IFCPAR (Indo-French Centre for the Promotion of Advanced Research). We thank G. Ferrera and A. Vicini for useful discussions.

Appendix A $\mathcal{\overline{G}}_{I,ij}^{(k)}$s of the soft distribution function

The constants $\mathcal{\overline{G}}_{I,ij}^{(k)}$ in the soft distribution function are given by,

[TABLE]

Appendix B Form factors

We present the analytic expressions of the form factors and the finite partonic cross sections for all the partonic channels. The labeling is same as Fig. 4

The unrenormalized form factor ($\hat{F}_{I}$) can be written as follows in the perturbative expansion of unrenormalized strong coupling constant ($\hat{a}_{s}$) and unrenormalized fine structure constant ($\hat{a}_{e}$)

[TABLE]

$I=q,b$ denotes the Drell-Yan pair production and the Higgs boson production in bottom quark annihilation, respectively. The coefficients ${\cal F}_{1}^{q},{\cal F}_{2,0}^{q},{\cal F}_{2,1}^{q}$ and ${\cal F}_{2,2}^{q}$ are

[TABLE]

The coefficients ${\cal F}_{1}^{b},{\cal F}_{2,0}^{b},{\cal F}_{2,1}^{b}$ and ${\cal F}_{2,2}^{b}$ are

[TABLE]

Appendix C $\Delta_{cd}^{(i,j)}$ for bottom quark annihilation from QCD, QED and QCD$\times$QED up to NNLO

In the following, we present finite partonic cross sections $\Delta_{cd}^{(i,j)}$ as defined in Eq. 5, up to NNLO level in the strong and electro-magnetic coupling constants. In QCD, $\Delta_{cd}^{i,0}$ for bottom quark annihilation is already known Harlander:2003ai ; Majhi:2010zg . In the following, $\Delta_{cd}^{i,0},i=1,2$ is in $SU(N)$ gauge theory, while $\Delta_{cd}^{0,j},j=1,2$ is in $U(1)$ gauge theory.

[TABLE]

The corresponding results from the QED and QCD$\times$QED are found to be

[TABLE]

Partonic cross sections contributing to pure NLO and NNLO QED corrections:

[TABLE]

The constants $\zeta_{i}=\sum_{k=1}^{\infty}\frac{1}{k^{i}}$, $k\in{\mathbb{N}}$ denote the Riemann’s $\zeta$-functions, e.g., $\zeta_{2}=1.64493406684822643647\ldots\,$ and $\zeta_{3}=1.20205690315959428540\ldots$. The Spence functions Li${}_{2}(x)$ and Li${}_{3}(x)$ are defined by

[TABLE]

and the Nielson function S${}_{1,2}(x)$ is given by

[TABLE]