  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.
