Two-loop QCD corrections to <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML" display="inline"> <mml:msup> <mml:mi>e</mml:mi> <mml:mo>+</mml:mo> </mml:msup> <mml:msup> <mml:mi>e</mml:mi> <mml:mo>−</mml:mo> </mml:msup> <mml:mo stretchy="false">

We present a next-to-next-to-leading-order (NNLO) calculation within the nonrelativistic QCD framework for the process <a:math xmlns:a="http://www.w3.org/1998/Math/MathML" display="inline"> <a:msup> <a:mi>e</a:mi> <a:mo>+</a:mo> </a:msup> <a:msup> <a:mi>e</a:mi> <a:mo>−</a:mo> </a:msup> <a:mo stretchy="false">→</a:mo> <a:msup> <a:mi>Z</a:mi> <a:mo>*</a:mo> </a:msup> <a:mo stretchy="false">→</a:mo> <a:mi>J</a:mi> <a:mo>/</a:mo> <a:mi>ψ</a:mi> <a:mo>+</a:mo> <a:mi>J</a:mi> <a:mo>/</a:mo> <a:mi>ψ</a:mi> </a:math> . We find that the NNLO contribution is two to three times larger than the next-to-leading-order contribution, which itself is already three to five times larger than the leading-order result. In the high-energy limit, we provide analytic expressions for the leading-power coefficients in the asymptotic expansion of the two-loop amplitudes. Our results are directly applicable to the bottomonium process <e:math xmlns:e="http://www.w3.org/1998/Math/MathML" display="inline"> <e:mi>Z</e:mi> <e:mo stretchy="false">→</e:mo> <e:mi mathvariant="normal">ϒ</e:mi> <e:mo>+</e:mo> <e:mi mathvariant="normal">ϒ</e:mi> </e:math> . Using the obtained hadronic amplitudes, we predict the decay width of the <j:math xmlns:j="http://www.w3.org/1998/Math/MathML" display="inline"> <j:mi>Z</j:mi> </j:math> boson into these rare dicharmonium and dibottomonium final states.