The discrepancy between the CDF measurement and the Standard Model theoretical predictions for the $W$-boson mass underscores the importance of conducting high-precision studies on the $W$ boson, which is one of the predominant objectives of proposed future $e^+e^-$ colliders. We investigate in detail the production of $W$-boson pair at $e^+e^-$ colliders and calculate the mixed QCD-EW corrections at the next-to-next-to-leading order. By employing the method of differential equations, we analytically compute the two-loop master integrals for the mixed QCD-EW corrections to $e^+e^- \rightarrow W^+W^-$. By utilizing the Magnus transformation, we derive a canonical set of master integrals for each integral family. This canonical basis fulfills a system of differential equations where the dependence on the dimensional regulator, $\epsilon$, is linearly factorized from the kinematics. Finally, these canonical master integrals are given as Taylor series in $\epsilon$ up to $\epsilon^4$, with coefficients written as combinations of Goncharov polylogarithms up to weight four. Upon applying our analytic expressions of these master integrals to the phenomenological analysis on $W$-pair production, we find that the $\mathcal{O}(\alpha\alpha_s)$ corrections hold substantial significance in the $\alpha(0)$ scheme, especially in the vicinity of the top-pair resonance ($\sqrt{s} = 2\, m_t$) induced by top-loop integrals. However, these corrections can be heavily suppressed by adopting the $G_{\mu}$ scheme.
Comment: 35 pages, 7 figures