Next-to-next-to-leading order $\text{QCD} \otimes \text{EW}$ corrections to $Z$-boson pair production at electron-positron colliders

We present a comprehensive analytic calculation of the next-to-next-to-leading order $\text{QCD} \otimes \text{EW}$ corrections to $Z$-boson pair production at electron-positron colliders. The two-loop master integrals essential to this calculation are evaluated using the differential equation method. In this work, we detail the formulation and solution of the canonical differential equations for the two-loop three-point master integrals with two on-shell $Z$-boson external legs and a massive internal quark in the loops. These canonical master integrals are systematically expanded as Taylor series in the dimensional regulator, $\epsilon = (4-d)/2$, up to the order of $\epsilon^4$, with coefficients expressed in terms of Goncharov polylogarithms up to weight four. Upon applying our analytic expressions of these master integrals to the phenomenological analysis of $Z$-pair production, we observe that the $\mathcal{O}(\alpha \alpha_s)$ corrections manifest at a level of approximately one percent compared to the leading-order predictions, underscoring their significance for comparisons with future high-precision experimental data.