Next-to-leading power two-loop soft functions for the Drell-Yan process at threshold

We calculate the generalized soft functions at $\mathcal{O}(\alpha_s^2)$ at next-to-leading power accuracy for the Drell-Yan process at threshold. The operator definitions of these objects contain explicit insertions of soft gauge and matter fields, giving rise to a dependence on additional convolution variables with respect to the leading power result. These soft functions constitute the last missing ingredient for the validation of the bare factorization theorem to NNLO accuracy. We carry out the calculations by reducing the soft squared amplitudes into a set of canonical master integrals and we employ the method of differential equations to evaluate them. We retain the exact $d$-dimensional dependence of the convolution variables at the integration boundaries in order to regulate the fixed-order convolution integrals. After combining our soft functions with the relevant collinear functions, we perform checks of the results at the cross-section level against the literature and expansion-by-regions calculations, at NNLO and partly at N$^3$LO, finding agreement.