  We present a novel next-to-next-to-leading order (NNLO) QCD calculation
matched to parton shower for the production of a pair of $Z$ bosons decaying to
four massless leptons, $p p \to \ell^+ \ell^- \ell'^+ \ell'^- + X$, at the LHC.
Spin correlations, interferences and off-shell effects are included throughout.
Our result is based on the resummed beam-thrust spectrum, which we evaluate at
next-to-next-to-leading-logarithmic (NNLL$'_{\mathcal{T}_0}$) accuracy for the
first time for this process, and makes use of the GENEVA Monte Carlo framework
for the matching to PYTHIA8 shower and hadronisation models. We compare our
predictions with data from the ATLAS and CMS experiments at 13 TeV, finding a
good agreement.
