We present a fully differential calculation of lepton pair production, taking into account the dominant next-to-next-to-leading order QED corrections as well as next-to-leading order electroweak and polarisation effects. We include all lepton masses, hard photon emission, as well as non-perturbative hadronic corrections. The corresponding matrix elements are implemented in the Monte Carlo framework MCMULE. In order to obtain a numerically stable implementation, we extend next-to-soft stabilisation, a universal technique based on a next-to-leading-power expansion, to calculations with polarised leptons. As an example, we show results tailored to the Belle II detector with the current setup as well as a potential future configuration that includes polarised beams.