  We compute the complete set of two-loop beam functions for the transverse
momentum distribution of the leading jet produced in association with an
arbitrary colour-singlet system. Our results constitute the last missing
ingredient for the calculation of the jet-vetoed cross section at small veto
scales at the next-to-next-to-leading order, as well as an important ingredient
for its resummation to next-to-next-to-next-to-leading logarithmic order. Our
calculation is performed in the soft-collinear effective theory framework with
a suitable regularisation of the rapidity divergences occurring in the
phase-space integrals. We discuss the occurrence of soft-collinear mixing terms
that might violate the factorisation theorem, and demonstrate that they vanish
at two loops in the exponential rapidity regularisation scheme when performing
a multipole expansion of the measurement function. As in our recent computation
of the two-loop soft function, we present the results as a Laurent expansion in
the jet radius $R$. We provide analytic expressions for all flavour channels in
$x$ space with the exception of a set of $R$-independent non-logarithmic terms
that are given as numerical grids. We also perform a fully numerical
calculation with exact $R$ dependence, and find that it agrees with our
analytic expansion at the permyriad level or better. Our calculation allows us
to define a next-to-next-to-leading order slicing method using the leading-jet
$p_T$ as a slicing variable. As a check of our results, we carry out a
calculation of the Higgs and $Z$ boson total production cross sections at the
next-to-next-to-leading order in QCD.
