We consider the production of a pair of heavy quarks \(Q\overline{Q }\) in association with a generic colour singlet system V at lepton colliders, and present the first analytic calculation of the two-loop soft function differential in the total momentum of the real radiation. The calculation is performed by reducing the relevant Feynman integrals into a canonical basis of master integrals by means of integration-by-parts identities. The resulting integrals are then evaluated by solving a system of differential equations in the kinematic invariants, whose boundary conditions are determined analytically with some care due to the presence of Coulomb singularities. The fully differential soft function is expressed in terms of Goncharov polylogarithms. This result is an essential ingredient for a range of N3LL resummations for key collider observables at lepton colliders, such as the \(Q\overline{Q }V\) production cross section at threshold and observables sensitive to the total transverse momentum of the radiation in heavy-quark final states. Moreover, it constitutes the complete final-final dipole contribution to the fully differential soft function needed for the description of \(Q\overline{Q }V\) production at hadron colliders, which plays an important role in the LHC physics programme.
