## D. Numerical evaluation 

Evaluation of the two-electron self-energy correction requires momentum-space calculation of the zero- and one-potential terms of the ”modified self-en ergy” contributions ( Δ E SE( A ) $_{SQED}$, Δ E SE( B ) $_{SQED}$, Δ E SE( E ) $_{SQED}$, and Δ E SE( G ) $_{SQED}$) and zero-potential terms of the ”modified vertex” contribu - tions ( Δ E SE( C ) $_{SQED}$, Δ E SE( F ) $_{SQED}$, and Δ E SE( H ) $_{SQED}$). For the zero- and one-potential terms of the ”modified self-energy” contributions we employ the numerical proced ure, developed for the self-energy dia- gram [51, 52]. 

The ”magnetic-vertex” part (with the operator T$_{0}$ in the vertex) of the ”modified vertex” con- tributions is somewhat similar to the vertex part of the one-electron self-energy correction to the g factor or the hyperfine splitting. Hence, our treatment of it s zero-potential term is based on the corresponding calculational procedures developed in Ref. [14] for the g factor and in Refs. [9, 36] for the hyperfine splitting. The angular integration, howev er, required significant generalization, in both cases, due to the interelectronic-interaction matrix elements. The ”interaction-vertex” part (with the interelectronic-interaction operator I (Δ) in the vertex) resembles the vertex part of the two-electron self-energy correction to the energy levels. Therefore, we have developed the nu- merical algorithm for the corresponding zero-potential term on the basis of the one presented in Ref. [53]. The main difference is the structure of the angular integrals which is substantially more complicated, due to the interaction with the external magnetic field. 

The many-potential terms of the ”modified self-energy” and ” modified vertex” parts, as well as the complete contribution of the ”double vertex” part are calculated in the coordinate space. Angular integration and summation over intermediate angular momentum projections is carried out in the standard way. The summation over the complete spectrum of the Dirac equation at fixed angular quantum numbers κ$_{1}$$_{,}$$_{2}$$_{,}$$_{3}$ is performed using the dual-kinetic-balance (DKB) approach [54] with the basis functions constructed from B -splines [55]. The infinite summation over κ$_{1}$$_{,}$$_{2}$$_{,}$$_{3}$ is terminated at | κ | = 10 – 15 and the rest of the sum is evaluated by the least-square inverse- polynomial fitting. In order to perform the integration over ω we employ two different contours, thus performing an additional cross-check. Both of them involve the Wick rotation. The first contour is the same as in Ref. [51]. The integration is performed along the imaginary axis. Besides, the contributions from the poles arising from the terms with ε$_{n}$ ≤ ε$_{a,b}$ must be calculated. The advantage of this contour is that the calculation is less time-consuming, in particular, due to the fact that Im[ I ( ω )] = 0 when Re[ ω ] = 0 . The second contour was proposed in Ref. [52]. The 

