We compute \( \mathcal{O}\left({\alpha}^2Z\right) \) radiative corrections to superallowed β decays with a heavy-particle effective field theory that systematically describes the interactions of low-energy ultrasoft photons with nuclei. We calculate two-loop virtual and one-loop real-virtual amplitudes by reducing the Feynman integrals to a set of master integrals, which we solve analytically using a variety of techniques. These techniques can be applied to other phenomenologically interesting observables. The ultrasoft corrections can then be combined with contributions arising from the exchange of potential photons to obtain the complete \( \mathcal{O}\left({\alpha}^2Z\right) \) correction to the decay rate, with resummation of large logarithms of the electron energy times the nuclear radius. We find that \( \mathcal{O}\left({\alpha}^2Z\right) \) ultrasoft loops induce a relative correction to the decay rate that ranges from 0.7 ∙ 10−3 in the decay of 10C to 3.6 ∙ 10−3 in the decay of 54Co, and will thus impact the extraction of Vud at the permille level. We show that the inclusion of these corrections reduces the residual renormalization scale dependence of the decay rate to a negligible level, making missing ultrasoft perturbative corrections a subdominant source of theoretical uncertainty.
