  At $e^+ \, e^-$ colliders the QED--initial state radiation forms a large part
of the radiative corrections. Their precise and fast evaluation is an essential
asset for the experiments at LEP, the ILC and the FCC-ee, operating at high
luminosity. A long standing problem in the analytic calculation of the
$O(\alpha^2)$ initial state corrections concerns a discrepancy which has been
observed between the result of Berends et al. (1988) \cite{Berends:1987ab} in
the limit $m_e^2 \ll s$ and the result by Bl{\"u}mlein et al. (2011)
\cite{Blumlein:2011mi} using massive operator matrix elements deriving this
limit directly. In order to resolve this important issue we recalculated this
process by integrating directly over the phase space without any approximation.
For parts of the corrections we find exact solutions of the cross section in
terms of iterated integrals over square root valued letters representing
incomplete elliptic integrals and iterations over them. The expansion in the
limit $m_e^2 \ll s$ reveals errors in the constant $O(\alpha^2)$ term of the
former calculation and yields agreement with the calculation based on massive
operator matrix elements, which has impact on the experimental analysis
programs. This finding also explicitly proofs the factorization of massive
initial state particles in the high energy limit including the terms of
$O(\alpha^2)$ for this process.
