  We present a calculation of all matching coefficients for $N$-jettiness beam
functions at next-to-next-to-next-to-leading order (N$^3$LO) in perturbative
quantum chromodynamics (QCD). Our computation is performed starting from the
respective collinear splitting kernels, which we integrate using the axial
gauge. We use reverse unitarity to map the relevant phase-space integrals to
loop integrals, which allows us to employ multi-loop techniques including
integration-by-parts identities and differential equations. We find a canonical
basis and use an algorithm to establish non-trivial partial fraction relations
among the resulting master integrals, which allows us to reduce their number
substantially. By use of regularity conditions, we express all necessary
boundary constants in terms of an independent set, which we compute by direct
integration of the corresponding integrals in the soft limit. In this way, we
provide an entirely independent calculation of the matching coefficients which
were previously computed in arXiv:2006.03056.
