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

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.


I. INTRODUCTION
The discovery of the Standard Model (SM) Higgs boson by ATLAS [1] and CMS [2] 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 [3]. 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) [4] for the Higgs boson production. Later on, in [5,6], 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 [4]. Thanks to the continued efforts beyond the NLO [7][8][9][10][11][12][13], the most precise prediction to date, namely next-to-next-to-next-to leading order (N 3 LO) prediction [14,15] for the inclusive production of the Higgs boson in the gluon fusion process, is now available (see [16][17][18] for rapidity distributions). In addition, at NNLO accuracy, the tiny effects due to finite top quark mass have already been computed in [19,20]. Electroweak (EW) corrections [21,22] and mixed QCD-EW corrections [23] 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 subdominant 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 [24] taking into account radiative corrections from QCD and EW are known for many of these processes.
Among these subdominant 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 [25] and CMS [26] collaborations have made dedicated efforts to measure them precisely. Among them, bottom Yukawa is one of the most sought 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 [27].
In the SM, bottom Yukawa coupling is less significant with respect to top Yukawa coupling while in the minimal supersymmetric SM (MSSM) [28] the coupling is proportional to 1= cos β which can increase the cross section in some parametric region. The angle β is related to the ratio, denoted by tan β, 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 [29][30][31], called 4FS and 5FS, respectively. In the former, one assumes that proton sea does not contain bottom quarks, 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 nonzero 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 [32][33][34] are known. On the other hand, in 5FS, NLO [35,36], NNLO [37] and the threshold effects at N 3 LO [13,38] (see [16,39] for rapidity distributions) are known for some time. Also, the 5FS cross section providing the dominant cross section in a matched prediction [40,41] is very well known. In [42], resummation of timelike logarithms in the Soft Collinear Effective Theory framework has been performed. Recently, for the bottom quark annihilation, complete N 3 LO corrections [43][44][45] have become available. In addition, the resummation of threshold contributions [46] at N 3 LO þ N 3 LL 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 subdominant 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 [45] and resummed [46] results solely from third order in perturbative QCD.
In [47], pure QED and mixed QCD × QED contributions have been obtained for the Drell-Yan (DY) process through Abelianization [48,49] at orders Oðα 2 Þ and Oðαα s Þ, respectively. In [47], a suitable algorithm is obtained by studying the group theory structure of QCD and QED amplitudes that contribute to the partonic subprocesses 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 subprocesses to the corresponding QED results. This way both pure QED as well as QCD × QED contributions to inclusive production cross section for the Z boson in the DY process have been obtained in [47] at NNLO level. Following this approach, we can in principle proceed to obtain pure QED and mixed QCD × QED contributions to the bottom quark annihilation process from the QCD results. Although the QCD results [37,50] to NNLO are presented for N ¼ 3 of SU(N) and hence Abelianization cannot be used, however, in [51], 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 [37,50] for N ¼ 3. Hence, we could use the results given in [51] and method of Abelianization to obtain pure QED as well as QCD × 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 × QED corrections to inclusive production of the Higgs boson in bottom quark annihilation up to NNLO level in U(1) and SUðNÞ × Uð1Þ. In addition, we reproduce the same for the production of Z boson in the 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 ε ¼ d − 4, d being the space-time dimension. The UV divergences are removed in 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 counterterms using the 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 × QED. We find that they are process independent. Using the universal IR anomalous dimension and following [52], we compute the renormalization constant for the Yukawa coupling in QED as well as in QCD × QED from the form factors (FF) of the Higgs bottom antibottom operator and vector current of DY process.
The paper is organized as follows. In Sec. II, after discussing the theoretical framework, we briefly describe in Sec. II A how we compute higher order QCD and QED radiative corrections to various partonic and photonic channels that contribute to the inclusive cross section. In Sec. II B, we discuss the UV and IR structure of the form factors and cross sections using the K þ G equation and obtain the mass factorized cross sections. In the following subsection, we discuss the Abelianization procedure. The phenomenological impact of our theoretical predictions are presented in Sec. III. Finally we summarize in Sec. IV. The universal constants that appear in the soft distribution function, FFs of vector current and bottom quarks and the mass factorized partonic and photonic cross sections are presented in the Appendices 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 where λ b is the Yukawa coupling which, after the EW symmetry breaking, is found to be m b =v. ψ b ðxÞ and m b denote the bottom quark field and mass, respectively. v is the vacuum expectation value (VEV) of the Higgs field ϕðxÞ. In the SM, the Higgs boson production through bottom quark annihilation is subdominant compared to gluon fusion through the 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 [53], tan β, the ratio of the VEVs of Higgs doublets can increase the contributions resulting from the bottom quark annihilation channel. At LO, where h is the SM-like light Higgs boson, H and A are the heavy and the pseudoscalar Higgs bosons, respectively. The parameter χ 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 [54][55][56] 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 where σ 0 is the Born cross section and f a ðx i ; μ 2 F Þ are parton distribution functions (PDFs) for a ¼ q;q; g and photon distribution function (PHDF) if a ¼ γ. The scaling variables x i are their momentum fractions. Δ cd are the partonic subprocess contributions normalized by the Born cross section. The scales μ R and μ 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. Δ cd can be expanded in powers of the QCD coupling constant a s ¼ g 2 s ðμ 2 R Þ=16π 2 and QED coupling constant a e ¼ e 2 ðμ 2 R Þ= 16π 2 , g s and e being the strong and electromagnetic coupling constants, respectively. That is, after suppressing μ R and μ F dependence, with Δ ð0;0Þ cd ¼ δð1 − zÞ and z ¼ q 2 =s. In the following, we describe the methodology to compute Δ ði;jÞ cd up to second order in the couplings.

A. Methodology
In this section, we briefly describe how higher order perturbative corrections Δ ði;jÞ cd [Eq. (5)] are computed. The details of computational procedure can be found in [57,58]. Beyond the leading order (LO), the partonic channels consists of one-and two-loop virtual subprocesses, realvirtual and single and double-real emissions, some of which are presented in Figs. 1-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.
Subprocesses 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 subprocesses. 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. We have used the program QGRAF [59] to generate virtual as well as real emission Feynman diagrams that contribute to the relevant subprocesses. An in-house FORM [60] code is used to perform all the symbolic manipulations, e.g. performing Dirac, SU(N) color and Lorentz algebra. A large number of loop integrals show up in the virtual diagrams. The integration-by-parts identities are used through a Mathematica based package, LITERED [61] to reduce them to a minimum set of master integrals (MI). For the virtual processes, at two-loop level, the form factors in QCD, QED and mixed QCD × QED require four MIs. For those processes that involve pure real emissions with or without virtual diagrams, we use the method of reverse unitarity [8] that allows one to use integration-by-parts identities to reduce the resulting phase-space integrals to a set of few MIs, the later can be found in [62]. For the real-virtual type of processes at NNLO level we need nine MIs for QCD, eight MIs for both QED and mixed QCD × QED. For the pure real emissions at NNLO level, we need 24 MIs for each of QCD, QED and mixed QCD × QED processes. Finally we obtain contributions to each subprocess, containing UV and IR singularities as poles in ε ¼ 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 × QED. In order to explore the IR structure, we study the production of a Z boson in hadron colliders, namely the DY process to the same accuracy in QCD, QED and QCD × 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 [13,52,63], 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 λ b up to two-loop level in QCD, QED and QCD × 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 × QED.

B. UV and IR structures in QED and QCD × QED
Having computed all the partonic channels that contribute to the hadronic cross sections in QED and QCD × 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 the Sudakov K þ G equation and for the later, following [12,13] we exploited the universal structure of the soft distribution function resulting from the soft gluon/photon emissions.
In order to remove the UV divergences that result from virtual subprocesses, we use the renormalization constants Z a c , c ¼ s, e for the QCD and QED coupling constants, respectively and Z λ b for the Yukawa coupling. The Yukawa coupling Z λ b receives contributions from both QCD and QED. Z a s and Z a e relate the bare couplingsâ s ¼ĝ 2 s =16π 2 of QCD andâ e ¼ê 2 =16π 2 of QED to the renormalized ones a s ðμ 2 R Þ and a e ðμ 2 R Þ, respectively, at the renormalization scale μ R in the following way: where a c ¼ fa s ; a e g. Here, S ε ≡ exp½ðγ E − ln 4πÞ ε 2 is the phase-space factor in d-dimensions, γ E ¼ 0.5772… is the Euler-Mascheroni constant and μ is an arbitrary mass scale introduced to makeâ s andâ e dimensionless in d-dimensions. The renormalization constant Z a c up to two loops is given by where β ij and β 0 ij are QCD and QED beta functions, respectively. In the present case, only one loop β, i.e. β 00 and β 0 00 [64], appear. They are given by 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 , n l is the number of active quark and lepton flavors and e q , e l refers to electric charge for quark q and lepton l, respectively. The renormalization constant Z b λ ða s ; a e Þ satisfies the renormalization group equation: whose solution in terms of the anomalous dimensions γ ði;jÞ b up to two loops is found to be Note that while the UV singularities factorize through Z λ b , singularities from QCD and QED mix from two loops onward. For QCD, γ ði;0Þ b is known to four loops [65]. In this paper, using the universal IR structure of the amplitudes and cross sections in QED, we determine γ ði;jÞ b up to two loops in QED i.e. for ði; jÞ ¼ ð0; 1Þ; ð0; 2Þ and in QCD × QED i.e. for ði; jÞ ¼ ð1; 1Þ.
We begin with the bare form factorsF I ðâ s ;â e ; Q 2 ; μ 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 the Sudakov integro-differential equation for these FFs, analogous to the QCD one. In dimensional regularization, they take the following form: where fa c g ¼ fa s ; a e g and Q 2 ¼ −q 2 is the invariant mass of the final state particle (dilepton 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 the one-loop level.
In other words, if we factorize IR singularities from the FFs, the resulting IR singular function cannot be written as a product of pure QCD and pure QED functions. More specifically, there will be terms proportional to a i s a j e , where i; j > 0, which will not allow factorization of QCD and QED ones. Hence, K I will have IR poles in ε from pure QED and pure QCD in every order in perturbation theory and in addition, from QCD × QED starting from 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 × 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 ε → 0 and the later contain only logarithms in q 2 . Note thatF I are renormalization group (RG) invariant so does the sum K I þ G I . Thus, the RG invariance ofF I implies 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 ðμ 2 R Þ and a e ðμ 2 R Þ as Unlike K I , G I do not contain any IR singularities but depend only on Q 2 and hence we expand them as where the first term is the boundary condition on each G I at μ 2 R ¼ Q 2 . Expanding A I in powers of a s and a e and using RG equations for QCD and QED couplings, we obtain The form factorsF I that we computed in this paper in QCD, QED and QCD × QED up to two-loop level can be used to extract the cusp anomalous dimensions (A ði;jÞ I ) by comparing them against Eq. (21). The explicit expressions for the unrenormalizedF I s are presented in Appendix B. We find A ði;jÞ I up to two loops as To obtain the process independent part of soft gluon/ photon contributions in the real emission subprocesses, we follow the method described in [12,13], where the soft distribution function for the inclusive cross section for producing a colorless state was obtained from the form factors and partonic subprocess cross sections involving real emissions of gluons. The soft distribution functions denoted by Φ 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 Φ b ¼ Φ q ¼ C F =C A Φ g holds up to the three-loop level [12,13,63]. We can use the partonic subprocesses of either DY process or the Higgs boson production in bottom quark annihilation namelyσ qq orσ bb normalized by the square of the bare form factorF q orF b to obtain Φ I . In general Φ I , which is a function of the scaling variable z ¼ q 2 =s, is defined as with Z q ¼ 1 and Z b ¼ Z λ b being the overall renormalization constant. The symbol C refers to "ordered exponential" which has the following expansion: Here ⊗ is the Mellin convolution and fðzÞ is a distribution of the kind δð1 − zÞ and D i . The plus distribution D i is defined as We can compute the UV finiteσ IĪ every order in renormalized perturbation theory. Since, we have not determined Z λ b , we can only compute the unrenormalized partonic cross sectionσ IĪ ¼σ IĪ =Z 2 I . From the explicit results forσ IĪ and the form factorsF I , using Eq. (25) we obtain Φ I up to second order in a s , a e and a s a e . We find Φ q ¼ Φ b up to second order in the couplings demonstrating the universality. In [12,13], it was shown that the soft distribution function Φ I satisfies the Sukakov K þ G equation analogous to the form factorF I due to similar IR structures that both of them have, order by order in perturbation theory. That is, Φ I satisfies where the IR singularities are contained inK and the finite part inḠ. RG invariance of Φ I implies Note that the same anomalous dimensions govern the evolution of bothK I andḠ I . This ensures that the soft distribution function contains right soft singularities to cancel those from the form factor leaving the bare partonic cross section to contain only initial state collinear singularities. The later will be removed by mass factorization by appropriate Altarelli-Parisi kernels. ExpandingK I ðfa c gÞ andḠ I ðfa c ðq 2 Þg; 1; ε; zÞ in powers of fa c g as has been done for K I ðfa c gÞ and G I ðfa c gÞ [see Eqs. (15) and (19) where, for up to two loops, G Assuming B [Eq. (22)] which is known to second order. They are found to be Alternatively, assuming B given in Eq. (22). Substituting the above UV anomalous dimensions in Eq. (10), we obtain Z λ b to second order in the couplings.
Using the renormalization constants Z a s , Z a e and Z λ b for the coupling constants α s , α e and the Yukawa coupling λ b , we obtain UV finite partonic cross sections. The soft and collinear singularities arising from gluons/photons/fermions in the virtual subprocesses cancel against those from the real subprocesses when all the degenerate states are summed up, thanks to the Kinoshita-Lee-Nauenberg theorem [66,67]. 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 Γ qq and Γ qg up to the two-loop level for the Uð1Þ and SUðNÞ × Uð1Þ cases. Since Γ qq and Γ qg are governed by the splitting functions P qq and P qg , we extract them to second order in couplings. In [48], these splitting functions up to NNLO level, both in QED and QCD × QED, were obtained using the Abelianization procedure. The splitting functions that we have obtained by demanding finiteness of the mass factorized cross section, agree with those in [48]. The mass factorized partonic cross section for each partonic subprocess up to NNLO in QED and in QCD × QED are presented in Appendix C along with the known NNLO QCD results [37]. In the next section, we use them to study their numerical impact at the LHC energies.

C. Abelianization procedure
In [47], QCD × 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 × QED results. Unlike in [47], 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 counterterms. Using these results at NNLO in QCD, QCD × 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 1 . This is absent for the processes where a photon is present instead of a 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 [47].
Rule 1: quark-quark initiated cases. QCD QCD × QED QED when both initial quarks are bottom quarks.
Rule 2: quark-gluon initiated cases: (after multiplying 2C A C F for the initial state gluon). QCD QCD × QED QED Rule 3: gluon-gluon initiated cases: (after multiplying 2C A C F for each initial state gluon).

III. RESULTS AND PHENOMENOLOGY
In this section, we study the numerical impact of pure QED and mixed QCD × 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 (c.m.) energy of ffiffiffi S p ¼ 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 [68], MRST [69], CT14 [70] and PDF4LHC17. The PDFs, PHDFs and the strong coupling constant a s can be obtained, using the LHAPDF-6 [71] interface. We have used the following input parameters for the masses and the couplings: Both a s ðμ R Þ and m b ðμ R Þ are evolved using appropriate QCD β-function coefficients and quark mass anomalous dimensions respectively. However, we have considered fixed α e ¼ 4πa e throughout the computation.
The Higgs boson production cross section from bottom quark annihilation at the present energy of LHC is not substantial. For example, at 13 TeV, the third order QCD corrections contribute almost 1% to the NNLO, while the mixed QCD × QED corrections contribute around 0.1% on top of the NNLO contributions. However, for the high luminosity LHC, measuring them at higher center of mass energy (c.m.) would give larger contributions and it will improve the precision. Hence, we have first studied how the cross section varies with the c.m. of LHC. In Fig. 4, we plot the inclusive production cross sections at various orders in perturbative QCD and QED for the range of c.m. energies between ffiffiffi S p ¼ 6 to 22 TeV. In the inset, the index "ij" indicates that QCD at ith order and QED at jth order in perturbative theory are included (e.g. "NNLO 11" indicates NNLO mixed QCD × 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 (μ R ) and factorization (μ 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 a larger Born contribution compared to pure QCD. However, if we consider the same running of Yukawa coupling, the NLO QCD effects are 50-500 times larger than NLO QED effects, depending on the scale choice. 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 × QED at various orders in perturbation theory which we have The Δ i;j indicates sole ith order QCD and jth order QED corrections to the total contribution. For example, Δ 11 means Δ 0;0 þ Δ 1;0 þ Δ 0;1 þ Δ 1;1 .
The Δ ij indicates the total contribution, in other figures, which is denoted by either LO, NLO or NNLO; e.g. Δ 11 means NNLO 11 .
In Table II, a similar study has been performed for ffiffiffi S p ¼ 13 TeV and the scales μ Fixed order predictions depend on the renormalization (μ R ) and factorization (μ F ) scales. The uncertainty resulting from the choice of the scales quantifies the missing higher order contributions. Hence, we have studied their dependence by varying them independently around a central scale. Figure 5 shows the dependence of the cross section on the renormalization scale (μ R ) for the fixed choice of the factorization scale μ F ¼ m h =4. It clearly demonstrates the importance of higher order corrections as the μ R variation is much more stable at NNLO 20 compared to the lower orders. In Fig. 6, we present the dependence on the factorization scale (μ F ) keeping the renormalization scale (μ R ) fixed at m h . Similar to the μ R variation, μ F variation improves after adding higher order corrections. To illustrate their dependence when both of the scales are changed simultaneously, we present the cross section by performing seven-point scale variation and the results are listed in Table III. 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 IV    In Table V, we repeat the study for ffiffiffi S p ¼ 13 TeV and μ R ¼ m h and μ F ¼ m h =4. We have also studied the uncertainties resulting from the choice of PDF set [71]. Using NNPDF31, in Fig 7, we plot the variation of the cross section with respect to different choices of PDF and PHDF replica. The central value and PDF uncertainties are given by the average and standard deviation over the replica sample, and are denoted in Fig. 7 by the thick line and shaded region, respectively.

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 subdominant 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 subdominant. Note that this is known to third order in QCD [45]. While this is a subdominant 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 nonfactorizable or mixed QCD × 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 the Sudakov K þ G equation to understand the IR structure in terms of cusp, collinear and soft anomalous dimensions. We demonstrate that the IR singularities from QCD, QED and QCD × 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/nonfactorizable QCD × 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 the QCD one, we find that it contains soft terms from mixed QCD × 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 × 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 [47]. Having obtained the complete NNLO results from QED and QCD × 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 × QED improve the reliability of the predictions.

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 the same as Fig. 4.