The Higgs boson decay into bottom quarks is the dominant decay channel contributing to its total decay width, which can be used to measure the bottom quark Yukawa coupling and mass. This decay width has been computed up to \( \mathcal{O}\left({\alpha}_s^4\right) \) for the process induced by the bottom quark Yukawa coupling, assuming massless final states, and the corresponding corrections beyond \( \mathcal{O}\left({\alpha}_s^2\right) \) are found to be less than 0.2%. We present an analytical result for the decay into massive bottom quarks at \( \mathcal{O}\left({\alpha}_s^3\right) \) that includes the contribution from the top quark Yukawa coupling induced process. We have made use of the optical theorem, canonical differential equations and the regular basis in the calculation and expressed the result in terms of multiple polylogarithms and elliptic functions. We propose a systematic and unified procedure to derive the ϵ-factorized differential equation for the three-loop kite integral family, which includes the three-loop banana integrals as a sub-sector. We find that the \( \mathcal{O}\left({\alpha}_s^3\right) \) corrections increase the decay width, relative to the result up to \( \mathcal{O}\left({\alpha}_s^2\right) \) , by 1% due to the large logarithms \( {\log}^i\left({m}_H^2/{m}_b^2\right) \) with 1 ≤ i ≤ 4 in the small bottom quark mass limit. The coefficient of the double logarithm is proportional to CA – CF, which is the typical color structure in the resummation of soft quark contributions at subleading power.
