R ( κ ) = κ 〈 ⃗ ∇ 2 〉 D + ln(2 κ ) + $_{4}$√$_{2}$$_{Z}$ κ 1 / 2 − 2 Z $^{2}$(ln κ + 1) κ . (A5) w ( k ) = pol ( 1 √$_{k}$ ) + ln k k pol ( 1 k ) , (A6) ln( k$_{0}$ ) = 1 D lim K →∞ [ 〈 ⃗ ∇ 2 〉 K + D ln(2 K ) + ∫ K 0 dk kJ ( k ) ] , (A1) where D = 2 πZ 〈 δ $^{3}$( r$_{1}$ ) + δ $^{3}$( r$_{2}$ ) 〉 , ⃗ ∇ ≡ ⃗ ∇ 1 + ⃗ ∇ $_{2}$, and J ( k ) = 〈 ⃗ ∇ 1 E$_{0}$ − H$_{0}$ − k ⃗ ∇ 〉 . (A2) ln( k$_{0}$ )$_{M}$ = m M [ R$_{M}$ ( κ ) + 1 D ∫ κ 0 dk k J$_{M}$ ( k ) + ∫ ∞ κ dk w$_{M}$ ( k ) k 2 ] , (A7) ln( k$_{0}$ ) = R ( κ ) + 1 D ∫ κ 0 dk k J ( k ) + ∫ ∞ κ dk w ( k ) k 2 , (A3) J$_{M}$ ( k ) = 2 〈 φ ∣ ⃗ ∇ 1 E$_{0}$ − H$_{0}$ − k ⃗ ∇ ∣ δφ 〉 + 〈 φ ∣ ⃗ ∇ 1 E$_{0}$ − H$_{0}$ − k [ δE − ⃗ p$_{1}$ · ⃗ p$_{2}$ ] 1 E$_{0}$ − H$_{0}$ − k ⃗ ∇ ∣ φ 〉 , (A8) w ( k ) = k 3 D J ( k ) + k 2 D 〈 ⃗ ∇ 2 〉 + k − $_{2}$√$_{2}$$_{Zk}$ 1 / 2 + 2 Z $^{2}$ln k , (A4) Acknowledgments Appendix A: Bethe logarithm 8 and R ( κ ) is a simple function obtained by integrating out the separated asymptotic expansion terms, theoretical calculations of helium-like atoms. In the present investigation, we performed a calculation of the energy levels of the n = 1 and n = 2 states of light helium-like atoms with the nuclear charge Z = 2 . . . 12, within the approach complete up to orders mα 6 and m $^{2}$/Mα$^{5}$. An extensive analysis of the 1 /Z expansion of individual corrections was carried out and comparison with results of the complementary approach was made whenever possible. Our general conclusion is that the results obtained within the approaches based on the Zα and the 1 /Z expansion are consistent with each other up to a high level of precision. However, further im- provement of numerical accuracy of the all-order, 1 /Z - expansion results and their extension into the lower- Z region is needed in order to safely merge the two comple- mentary approaches. The calculational scheme employed for the evaluation of Eq. (A3) is similar to that previously used [15] for the relativistic corrections to the Bethe logarithm. At the first step, a careful optimization of nonlinear basis-set parameters was carried out for a sequence of scales of the photon momentum: k$_{i}$ = 10 i and i = 1 , . . . , i$_{max}$ , with i$_{max}$ = 5 for the S states and i$_{max}$ = 4 for the P states. The optimization was performed with incrementing the size of the basis until the prescribed level of convergence is achieved for the function w ( k ). At the second step, the integrations of the photon momentum k were performed. For a given value of k , the function J ( k ) was calculated with a basis obtained by merging together the optimized bases for the two closest k$_{i}$ points, thus essentially dou- bling the number of the basis functions. The function w ( k ) was obtained from J ( k ) according to Eq. (A4). We wish to thank Tomasz Komorek for participation in the beginning of the project. Computations reported in this work were performed on the computer clusters of St. Petersburg State Polytechnical University and of Institute of Theoretical Physics, University of Warsaw. The research was supported by NIST through Precision Measurement Grant PMG 60NANB7D6153. V.A.Y. ac- knowledges additional support from the “Dynasty” foun- dation and from RFBR (grant No. 10-02-00150-a). The integral over k ∈ [0 , κ ] was calculated analytically, by diagonalizing the Hamiltonian matrix and using the spectral representation of the propagator. The value of the auxiliar parameter κ was set to κ = 100. The integral over k ∈ [ κ, ∞ ) was separated into two parts, k < 10 i max and k > 10 i max . The first part was evaluated numeri- cally by using the Gauss-Legendre quadratures, after the change of variables t = 1 /k $^{2}$. The second part was ob- tained by integrating the asymptotic expansion of the function w ( k ). The coefficients of this expansion were obtained by fitting the numerical data for w ( k ) to the form Following the approach of Refs. [34, 35], the Bethe log- arithm (13) is expressed in terms of an integral over the momentum of the virtual photon, where pol( x ) denotes a polynomial of x . The total num- ber of fitting parameters was about 9 − 11. The range of k to be fitted and the exact form of the fitting function were optimized so as to yield the best possible results for the known asymptotic expansion terms of J ( k ). The first-order perturbation of the Bethe logarithm by the mass-polarization operator can be represented as [10] For the purpose of numerical evaluation, the integration over the photon momentum k is divided into two regions by introducing the auxiliar parameter κ , where where the function w ( k ) represents the residual obtained from J ( k ) by removing all known terms of the large- k asymptotics,