We describe in detail the calculation of the two-loop corrections to the QED Bhabha scattering cross section due to the vacuum polarization by heavy fermions. Our approach eliminates one mass scale from the most challenging part of the calculation and allows us to obtain the corrections in a closed analytical form. The result is valid for arbitrary values of the heavy fermion mass and the Mandelstam invariants, as long as s, t, u ≫ m 2 e .High energy electron-positron or Bhabha scattering [1] is among the most important and carefully studied processes in particle physics. It provides a very efficient tool for luminosity determination at electron-positron colliders and thus mediates the process of extracting physical information from the raw experimental data [2]. The small-angle Bhabha scattering is particularly effective as a luminosity monitor at high-energy colliders. 4 The large-angle Bhabha scattering is used to measure the luminosity at colliders operating at the center-of-mass energy, √ s, of a few GeV, such as BABAR/PEP-II, BELLE/KEKB, BES/BEPC, KLOE/DAΦNE, and CMD, SND/VEPP-2M [5]. 5 Moreover, it will be also used to disentangle the luminosity spectrum at the ILC [6,7]. Bhabha scattering involves stable charged leptons both in the initial and the final states and, therefore, it can be measured experimentally with very high precision. At LEP, the experimental error in the luminosity measurement has been reduced to 0.4 permille [8] and it is expected to be even smaller at the ILC: the goal of the TESLA forward calorimeter collaboration is to reach the experimental accuracy of 0.1 permille in the first year of run [9]. Finally, at the low-energy accelerators DAΦNE and VEPP-2M the cross section of the large-angle scattering is measured with the accuracy of about 1 permille [10,11]. In the phenomenologically most interesting cases of low energy or small angle scattering, the Bhabha cross section is QED dominated, with the electroweak and hadronic effects being strongly suppressed. Therefore, it can be reliably computed in perturbative QED, with the accuracy limited only by uncalculated high order corrections. These properties make in such a way that Bhabha scattering is an ideal "standard candle" for electron-positron colliders.Realistic simulations of the Bhabha events, which take into account the detector geometry and experimental cuts, are performed by means of sophisticated Monte Carlo generators, such as BHLUMI [12], BABAYAGA [5,13], BHAGENF [14], BHWIDE [15], MCGPJ [16], and SABSPV [17]. To match the experimental needs, the two-loop QED corrections must be included into the theoretical analysis and incorporated into the event generators. Since the theoretical accuracy directly affects the luminosity determination and may jeopardize the high-precision physics program at electron-positron colliders, remarkable efforts were devoted to the study of the radiative corrections. The one-loop corrections have been known in the full electroweak theory for a long time [18]. The two-loop electroweak corrections ar...