This paper presents the calculation of the total cross-section for direct Higgs boson production in hadron collisions at next-to-next-to-leading order (NNLO) in perturbative QCD. A new algorithmic technique is introduced for evaluating inclusive phase-space integrals, based on the Cutkosky rules, integration by parts, and the differential equation method for computing master integrals. The paper discusses the numerical impact of the O(αs²) QCD corrections to the Higgs boson production cross-section at the LHC and the Tevatron. The Higgs boson is the only missing particle in the Standard Model (SM) of electroweak interactions. Its discovery would confirm the SM and provide insights into mass generation and physics beyond the SM. Direct searches at LEP restrict the Higgs boson mass to be greater than 114.1 GeV, while precision electroweak measurements suggest a value around 90 GeV. The SM remains perturbative up to about 1 TeV, suggesting a relatively light Higgs boson that could be observed at the Tevatron or LHC. Gluon fusion through top-quark loops is expected to be the dominant Higgs production mechanism at these facilities. Theoretical estimates of the cross-section for Higgs boson production via gluon fusion, based on computations up to next-to-leading order (NLO) in perturbative QCD, are insufficient. The leading-order (LO) cross-section is proportional to αs²(μ²), and exhibits a strong dependence on the choice of the scale μ. Including O(αs) corrections decreases the scale dependence but increases the cross-section by about 70%. Therefore, evaluating the next order in the perturbative expansion is important to enhance the credibility of the theoretical predictions. To compute the cross-section to NNLO, the matrix elements for the O(αs²) virtual corrections to gg → H, the O(αs) virtual corrections to gg → Hg, qg → Hq, and q̄q → Hg, and the tree-level matrix elements for the processes gg → Hqg, gg → Hq̄q, qg → Hqg, q̄q → Hgg, and q̄q → Hq̄q must be combined. The inclusive cross-section requires integrating over the loop-momenta in the virtual amplitudes and the phase-space of the real particles in the final state. Both real and virtual corrections are divergent in four dimensions. The amplitudes are regularized using dimensional regularization (d = 4 - 2ε), and the ultraviolet divergences are removed by renormalizing in the MS scheme. The remaining divergences arise from initial state collinear radiation and are absorbed into the parton distribution functions, yielding a finite cross-section. The calculation can be simplified byThis paper presents the calculation of the total cross-section for direct Higgs boson production in hadron collisions at next-to-next-to-leading order (NNLO) in perturbative QCD. A new algorithmic technique is introduced for evaluating inclusive phase-space integrals, based on the Cutkosky rules, integration by parts, and the differential equation method for computing master integrals. The paper discusses the numerical impact of the O(αs²) QCD corrections to the Higgs boson production cross-section at the LHC and the Tevatron. The Higgs boson is the only missing particle in the Standard Model (SM) of electroweak interactions. Its discovery would confirm the SM and provide insights into mass generation and physics beyond the SM. Direct searches at LEP restrict the Higgs boson mass to be greater than 114.1 GeV, while precision electroweak measurements suggest a value around 90 GeV. The SM remains perturbative up to about 1 TeV, suggesting a relatively light Higgs boson that could be observed at the Tevatron or LHC. Gluon fusion through top-quark loops is expected to be the dominant Higgs production mechanism at these facilities. Theoretical estimates of the cross-section for Higgs boson production via gluon fusion, based on computations up to next-to-leading order (NLO) in perturbative QCD, are insufficient. The leading-order (LO) cross-section is proportional to αs²(μ²), and exhibits a strong dependence on the choice of the scale μ. Including O(αs) corrections decreases the scale dependence but increases the cross-section by about 70%. Therefore, evaluating the next order in the perturbative expansion is important to enhance the credibility of the theoretical predictions. To compute the cross-section to NNLO, the matrix elements for the O(αs²) virtual corrections to gg → H, the O(αs) virtual corrections to gg → Hg, qg → Hq, and q̄q → Hg, and the tree-level matrix elements for the processes gg → Hqg, gg → Hq̄q, qg → Hqg, q̄q → Hgg, and q̄q → Hq̄q must be combined. The inclusive cross-section requires integrating over the loop-momenta in the virtual amplitudes and the phase-space of the real particles in the final state. Both real and virtual corrections are divergent in four dimensions. The amplitudes are regularized using dimensional regularization (d = 4 - 2ε), and the ultraviolet divergences are removed by renormalizing in the MS scheme. The remaining divergences arise from initial state collinear radiation and are absorbed into the parton distribution functions, yielding a finite cross-section. The calculation can be simplified by