  The Higgs boson decay to a massive bottom quark pair provides the dominant
contribution to the Higgs boson width. We present an exact result for such a
decay induced by the bottom quark Yukawa coupling with next-to-next-to-leading
order (NNLO) QCD corrections. We have adopted the canonical differential
equations in the calculation and obtained the result in terms of multiple
polylogarithms. We also compute the contribution from the decay to four bottom
quarks which consists of complete elliptic integrals or their one-fold
integrals. The result in the small bottom quark mass limit coincides with the
previous calculation using the large momentum expansion. The threshold
expansion exhibits power divergent terms in the bottom quark velocity, which
has a structure different from that in $e^+e^-\to t\bar{t}$ but can be
reproduced by computing the corresponding Coulomb Green function. The NNLO
corrections significantly reduce the uncertainties from both the
renormalization scale and the renormalization scheme of the bottom quark Yukawa
coupling. Our result can be applied to a heavy scalar decay to a top quark
pair.
