The $O(\alpha^2)$ Initial State QED Corrections to $e^+e^- \rightarrow \gamma^*/Z_0^*$

We calculate the complete $O(\alpha^2)$ initial state radiation corrections to $e^+ e^-$ annihilation into a neutral vector boson in a direct analytic computation without any approximation. The corrections are represented in terms of iterated incomplete (elliptic) integrals over alphabets of square--root valued letters. Performing the limit $s \gg m_e^2$, we find discrepancies with the earlier results of Ref.~\cite{Berends:1987ab} and confirm results obtained in Ref.~\cite{Blumlein:2011mi} where the effective method of massive operator matrix elements has been used, which works for all but the power corrections in $m^2_e/s$. In this way, we also confirm the validity of the factorization of massive partons in the Drell--Yan process. We add non--logarithmic terms at $O(\alpha^2)$ which have not been considered in previous calculations. The final results in the limit $s \gg m_e^2$ can be given in terms of Nielsen integrals.