Atomic structure calculations of helium with correlated exponential functions

The technique of quantum electrodynamics (QED) calculations of energy levels in the helium atom is reviewed. The calculations start with the solution of the Schr\"odinger equation and account for relativistic and QED effects by perturbation expansion in the fine-structure constant $\alpha$. The nonrelativistic wave function is represented as a linear combination of basis functions depending on all three interparticle radial distances, $r_1$, $r_2$ and $r = |\vec{r}_1-\vec{r}_2|$. The choice of the exponential basis functions of the form $\exp(-\alpha r_1 -\beta r_2 -\gamma r)$ allows us to construct an accurate and compact representation of the nonrelativistic wave function and to efficiently compute matrix elements of numerous singular operators representing relativistic and QED effects. Calculations of the leading QED effects of order $\alpha^5m$ (where $m$ is the electron mass) are complemented with the systematic treatment of higher-order $\alpha^6m$ and $\alpha^7m$ QED effects.


I. INTRODUCTION
The helium atom is the simplest many-body atomic system in the nature. Since the advent of quantum mechanics, helium was used as a benchmark case for developing and testing various calculational approaches of many-body atomic theory. Today, the nonrelativistic energy of various helium electronic states can be computed with an essentially arbitrary numerical accuracy [1,2]. The same holds also for the leading-order relativistic correction. Subsequently, the quantum electrodynamics (QED) effects in the atomic structure of helium can be clearly identified and studied by comparison of theoretical predictions with the large body of available experimental data. Experimental investigations of helium spectra have progressed rapidly over the years, recently reaching the precision of a few tens of Hertz [3].
For light atomic systems such as helium, relativistic and QED corrections to energy levels can be systematically accounted for by the perturbation expansion in the fine-structure constant α. The starting point of the expansion is the nonrelativistic energy of order α 2 m (= 2 Ry, where m is the electron mass and Ry is the Rydberg energy). The leading relativistic correction is of order α 4 m, whereas QED effects enter first in order α 5 m. A large body of work has been done in recent years to calculate QED effects in helium spectra. Extensive calculations of helium energies were accomplished by Gordon Drake et al. [4][5][6]. Their calculations are complete through order α 5 m and approximately include some higher-order QED effects. The next-order α 6 m QED correction was for a long time known only for the fine-structure intervals [7,8]. For individual energy levels, these effects were derived and calculated numerically by one of us (K.P.) [9][10][11]. The higher-order α 7 m QED effects were evaluated by us first for the fine structure [12][13][14] and just recently for the triplet n = 2 states of helium [15][16][17].
The purpose of this article is to review and systematize the technique of calculations of the helium atomic structure, developed in numerous investigations over the last three decades. The starting point of the calculations is the Schrödinger equation, which is solved variationally after expanding the wave function into a finite set of explicitly-correlated basis functions depending on all three interparticle radial distances. It has been known for a long time [18] that inclusion of the interelectronic distance explicitly into the basis set is crucially important for constructing an accurate representation of the two-electron wave function. Moreover, it has also long been recognized [19] that an accurate wave-function representation should satisfy the so-called cusp conditions at the two-particle coalescence points | r i − r j | = 0. The cusp condition is expressed [20], after averaging over angles and for the singlet states, as where r is an interparticle distance and the parameter λ = 1/2 for the electron-electron and λ = −Z for the electron-nucleus cusp (where Z is the nuclear charge number).
The two most succesful basis sets used in the literature for high-precision calculations of the atomic structure of helium are: the Hylleraas basis set adopted by Drake et al. [4][5][6] and the exponential basis set put forward by Korobov [21,22] and used in numerous calculations of our group. Both these basis sets are explicitly correlated and are able to reproduce the cusp conditions with great accuracy. In the present work we will concentrate on the exponential basis set, because only this basis has been successfully used in calculations of higher-order QED effects so far.

II. WAVE FUNCTIONS
The spatial wave function ψ LML with a specified total angular momentum L and its momentum projection M L for a two-electron atom is standardly represented as where f l1l2 is the radial part of the wave function, r = r 1 − r 2 , andr = r/r. Furthermore, Y l1l2 LML are the bipolar spherical harmonics, where l 1 m 1 l 2 m 2 |lm is the Clebsch-Gordan coefficient and Y lm are the spherical harmonics. We stress that the radial part of the wave function is assumed to be explicitly correlated, i.e., the function f depends on all interparticle distances, r 1 , r 2 , and r. In this case, the sum over l 1 and l 2 in Eq. (2) is restricted [23] by two conditions (A) : l 1 + l 2 = L , or (B) : l 1 + l 2 = L + 1 , (4) which lead to wave functions of different parities (−1) l1+l2 . The bipolar spherical harmonics are usually handled in the spherical coordinates using the apparatus of Racah algebra, see, e.g., Ref. [24]. We find, however, that calculations with explicitly correlated functions are more conveniently performed in Cartesian coordinates. One of the reasons is that the action of numerous momentum operators encountered in calculations is most easily evaluated in the Cartesian coordinate system. The corresponding calculations can easily be automatized and performed with the help of systems of symbolic computations.
For this purpose, the expansion of the wave function is more conveniently made in terms of the bipolar solid harmonics. In order to define them, we start with the solid harmonics, where the normalization coefficient A L is fixed below. The solid harmonics obey the following summation rule, where (r i r j r k . . .) (L) is a traceless and symmetric tensor of the order L constructed from components of the vector r with Cartesian indices i, j, k . . .. and the summation over these Cartesian indices is implicit. The last equation determines A L , which is related to the coefficient of x L in the Legendre polynomial P L (x), specifically, We now define the bipolar solid harmonics Y l1l2 LML as where R ≡ r 1 × r 2 , ξ is an arbitrary vector, and the righthand-side of the above equations does not depend on ξ after the L-fold differentiation. The bipolar solid harmonics are proportional to the corresponding bipolar spherical harmonics with a prefactor that does not depend on angles, so their angular parts are exactly the same. Now, using Eq. (6), we obtain that the bipolar solid harmonics Y l1l2 LML obey the analogous summation rule where Y l1l2 i1..iL are the symmetric and traceless tensors of rank L with Cartesian indices i 1 . . . i L , The summation formula (10) shows that the matrix elements with the spatial wave function can be represented in terms of matrix elements with the Cartesian wave function as follows where Q is an arbitrary spatial operator. Eq. (14) is the Cartesian representation of the spatial wave function used in the present work.
We now present explicit formulas for the Cartesian wave functions for different values of the angular momentum and parity. For L = 0 we have l 1 = l 2 = 0 and only even parity. The wave function is just a scalar, where the upper sign in ± corresponds to the singlet and the lower sign to the triplet state. For L = 1, we have (l 1 , l 2 ) = (0, 1), (1, 0) for the odd parity and (l 1 , l 2 ) = (1, 1) for the even parity. The corresponding wave functions are vectors, The L = 2 odd and even wave functions are second-rank tensors, where we suppressed arguments of the radial functions F and G and the elementary second-rank tensors are defined as Explicit expressions for the L = 3 and L = 4 functions can be found in Appendix A of Ref. [25]. The spatial wave functions are normalized by

III. EVALUATION OF MATRIX ELEMENTS
The spin-dependent wave function with definite values of the total momentum J, its projection M , the angular momentum L, and the spin S is given by where M S is the spin projection, χ SMS is the spin function, and ψ LML is the spatial wave function. As described in the previous section, in our calculations we evaluate all matrix elements in Cartesian coordinates. The spatial wave function with the angular momentum L is represented in the form (14); namely, as a traceless tensor of rank L symmetric in all Cartesian indices carried by r 1 , r 2 , and r 1 × r 2 . In addition, it is assumed that the wave function has a definite symmetry with respect to r 1 ↔ r 2 .
The norm and the expectation value of any spinindependent operator are immediately reduced to the spatial matrix element, where the summation over Cartesian indices is implicit. This equation is sufficient for determining the nonrelativistic wave function and the nonrelativistic energy. The relativistic and QED corrections involve operators depending on the electron spin. The expectation value of an arbitrary operator Q on a state with definite J, for the singlet S = 0 states, is expressed as where I is the unity matrix, J = L + S, S = ( σ 1 + σ 2 )/2, and the trace is performed in the 4-dimensional space of two spins. Further evaluation of the matrix element proceeds by performing the trace of the operators in the spin space, with help of the following trace rules, In the case of a spin-independent operator Q, Eq. (25) is reduced to Eq. (24). For the triplet states one considers three values of J = L−1, L, L+1. The expectation value then takes the form For the spin-independent operators, this equation is equivalent to Eq. (24). The coefficients A JL and B JL are obtained by considering two particular cases, (2) . The left-hand-side of Eq. (31) is then immediately expressed in terms of J and L, whereas the right-hand-side is evaluated by using This consideration gives for J = L+1, L, and L−1, correspondingly. These are all the formulas needed to factorize out the spin dependence of matrix elements and to express them in terms of spatial integrals.
The expectation values of an arbitrary operator Q for the singlet and triplet wave functions are obtained from Eqs. (25) and (31). We now write explicitly the corresponding expressions. The results for the S states are For the P states, we obtain The results for the D states are

IV. INTEGRALS WITH EXPONENTIAL BASIS FUNCTIONS
The radial parts of the wave function (14) are represented as linear combinations of the exponential basis functions, where c k are linear coefficients, N is the size of the basis, and α k , β k , and γ k are nonlinear parameters obtained in the process of the basis optimization. One of the great features of the exponential basis functions is that the evaluation of radial integrals is very simple. A calculation of radial matrix elements of various operators with wave functions of the form (46) is reduced to evaluation of the integrals I(i, j, k), For matrix elements of the nonrelativistic Hamiltonian, only integrals with non-negative values of i, j, and k are required. All such integrals can be obtained by differentiation of the master integral I(0, 0, 0) over the nonlinear parameters, for n i , n j , n k ≥ 0. The expression for the master integral I(0, 0, 0) is very simple, Matrix elements of relativistic corrections involve integrals with additional inverse powers of r 1 , r 2 , and r, whose evaluation requires two additional master integrals. Their expression can be obtained by integrating Eq. (49) with respect to the corresponding nonlinear pa- where Li 2 is the dilogarithm function [26]. Other integrals for relativistic corrections are obtained by differentiating the above formulas for master integrals. We note that Eq. (50) contains a spurious singularity at α = β. The zero in the denominator is compensated by the vanishing logarithm function and thus is not a real singularity but can lead to numerical instabilities. In order to transform Eq. (50) to an explicitly regular form, we introduce a regularized logarithm function ln 1 (x) by separating out the first term of the Taylor expansion, Introducing ln 1 (x) with x = (α − β)/(β + γ) in Eq. (50), we obtain a regular representation of this formula. In practical calculations we encounter more spurious singularities of this kind. They are eliminated with the help of functions ln n (x), which are introduced analogously to ln 1 (x) by separating n first terms of the Taylor expansion of ln(1 + x). Matrix elements of QED corrections involve several integrals with large negative powers of radial distances, like 1/r 3 , 1/r 4 , and even 1/r 5 . Such integrals are singular and need proper definitions. With the exponential functions, it is possible to obtain simple and numerically stable representations for such integrals. The corresponding procedure is described in Appendix A. Numerical results for basic singular integrals for the 2 3 S and 2 3 P states of helium are presented in Table I.
In our calculations of the α 7 m QED effects [27], integrals with ln r were encountered for the first time, where γ E is the Euler gamma constant. Such integrals are evaluated with the help of the following set of master integrals [27]: where Li 3 is the trilogarithm function [26]. Eq. (56) is valid for α > β. The corresponding result for α < β is obtained by the analytic continuation with help of the following identities [26] The result for the case of α = β is straightforwardly obtained from Eq. (56).

V. NONRELATIVISTIC ENERGY AND WAVE FUNCTION
The nonrelativistic Hamiltonian of the helium atom for the infinitely heavy nucleus is where p a = −i ∇ a is the momentum operator of the electron a and Z is the nuclear charge number (Z = 2 for helium). The Schrödinger equation is A direct solution of the Schrödinger equation is standardly substituted by the problem of finding the minimum or, generally, a stationary point of the variational functional The variational eigenvalues obtained from this functional are the upper bounds to the true eigenvalues, and the corresponding eigenvectors provide the linear coefficients c k of the basis-set expansion (46). It is important that the variational principle works equally well for the ground and for the excited states. The finite nuclear mass correction to the nonrelativistic energy is induced by the nuclear kinetic energy operator where M is the nuclear mass and P = − p 1 − p 2 is the nuclear momentum. There are two ways to incorporate the nuclear mass effect to the nonrelativistic energy: (i) to include the operator δ M H into the nonrelativistic Hamiltonian H 0 and solve the nuclear-mass dependent Schrödinger equation and (ii) to solve the Schrödinger equation for the infinitely heavy nucleus and to account for the nuclear mass effects by perturbation theory. In our calculations with the exponential basis we found that the inclusion of δ M H into the nonrelativistic Hamiltonian leads to numerical instabilities for S states (but not for P and higher-L states). So, for S states we account for the nuclear mass effects by perturbation theory (up to the third order in 1/M [28]), whereas for the P and D states we usually include δ M H into the solution of the Schrödinger equation. We checked that for the P and D states both methods yield equivalent results.
It should be mentioned that in the literature it is customary to split the operator δ M H into the mass-scaling and mass-polarization parts, The effect of the mass scaling (caused by the first term in Eq. (63)) can be incorporated into the nonrelativistic Hamiltonian (59) by switching to the reduced-mass atomic units r → µ r, where µ = 1/(1 + m/M ) is the reduced mass. As a result, the mass-scaling term leads to the appearance of the reduced mass prefactor in the nonrelativistic energy E 0 → µ E 0 and only the mass polarization term needs to be accounted for separately. We find it more convenient to keep the nuclear kinetic energy operator in the closed form (62), because this greatly simplifies consideration of higher-order recoil QED effects. Because the nonrelativistic Hamiltonian H 0 does not depend on spin, its matrix elements are immediately reduced to radial integrals with the spatial wave functions according to Eq. (24). Computing the action of gradients ∇ 1,2 on the wave functions (14), we express the matrix elements ψ|H 0 |ψ as a linear combination of integrals I(i, j, k) with i, j, k ≥ 0, which are rational functions of the nonlinear parameters α n , β n , and γ n .
The choice of the nonlinear basis parameters α n , β n , and γ n is crucially important for obtaining an accurate and compact representation of the wave function and the energy E 0 . The general approach is to perform the variational optimization of the basis parameters, by searching for a minimum of the eigenvalue of the Hamiltonian matrix corresponding to the desired reference state. Because the optimization of each individual nonlinear parameter is not effective from the computational point of view, we use the approach introduced by Vladimir Korobov [21]. In this method, the (real) nonlinear parameters α, β, and γ are quasirandomly distributed in the intervals and the parameters A 1,2 , B 1,2 , and C 1,2 are determined by the variational optimization. We note that the nonlinear parameters as well as A 1,2 , B 1,2 , and C 1,2 can be both positive and negative. However, in order to ensure the normalizability of the wave function and its physical behavior at large r 1 , r 2 , and r, we require that where ǫ ∼ √ 2 E io , with E io being the ionization energy. The performance of the basis set can be significantly improved if one introduces several sets of intervals A 1,2 , B 1,2 , and C 1,2 which are optimised variationally. In our calculations we use typically two or three sets of intervals. This can be considered as an analogue of several different exponential scales in the Hylleraas-type calculations by Drake at al. [6,29].
We also note that in calculations for excited 1snl states it is advantageous to include several screened hydrogenic wave functions of the type φ Z 1s ( r 1 ) φ Z−1 n ′ l ( r 2 ) with n ′ ≤ n in the basis, whose parameters are excluded from optimization. This ensures that the variational optimization is localized at the local minimum with the desired principal quantum number n and does not collapse to lower n's.
Our procedure for determination of the nonrelativistic wave function and energy looks as follows. For a given size of the basis N the nonlinear parameters α n , β n , and γ n with n = 1, . . . , N are distributed quasirandomly within the initial set of intervals with parameters A i , B i , and C i . Then, the N × N matrix of the nonrelativistic Hamiltonian H 0 is computed. The linear coefficients c n and the desired reference-state eigenvalue E 0 are determined by the inverse iteration method. The inversion of the Hamiltonian matrix is performed by the LDU decomposition method. This procedure is repeated for different sets of the parameters A i , B i , and C i , searching for the minimum value of the energy eigenvalue.
A disadvantage of working with the exponential basis is that the basis quickly degenerates as N is increased (i.e. the determinant of the Hamiltonian matrix becomes very small), which leads to numerical instabilities in linear algebra routines. Because of this the usage of an extended-precision arithmetics is mandatory. In our calculations we used the Fortran 95 libraries for the octupleprecision (about 64 digits) arithmetics written by V. Korobov [30], quad-double routine by D. H. Bailey, and MP-FUN/MPFR library by D. H. Bailey [31]. Table II shows an example of the convergence of numerical results with the exponential basis with increase of the basis size. We observe that with just N = 200 basis functions one obtains the nonrelativistic energy with about 10-digit accuracy.

VI. RELATIVISTIC CORRECTION
The relativistic correction splits the nonrelativistic energy levels with quantum numbers L > 0 and S > 0 into sublevels according to the value of the total momentum J. This effect is known as the fine structure. It is often convenient to consider separately the centroid energy levels obtained by averaging over all J sublevels, and the fine-structure intervals between individual J sublevels. The centroid energy is defined as The relativistic correction is induced by the Breit Hamiltonian, which is conveniently separated into the spin-independent and the spin-dependent parts, In the leading order of perturbation theory, the spinindependent part H A contributes only to the centroid energy, whereas the spin-dependent part H fs causes the fine structure splitting.

A. Centroid energy
The spin-independent part of the Breit Hamiltonian is given by where P = − p 1 − p 2 is the nuclear momentum. In order to account for the finite nuclear mass effects, the expectation value of the operator H A should be evaluated with the eigenfunctions ψ M of the Schrödinger Hamiltonian with the finite nuclear mass (i.e. the sum of Eqs. (59) and (62)). Alternatively, the wave function ψ M can be constructed by perturbation theory in 1/M . In our calculations, we include the nuclear recoil effect for the relativistic correction perturbatively for the S states, and nonperturbatively for the L > 0 states. The matrix element of H A is reduced to the radial integral with the spatial wave functions according to Eq. (24) and can be evaluated numerically. However, the expectation values of the operators p 4 a and δ 3 (r a ) are slowly converging with respect to the size of the basis because these operators are nearly singular. It is possible to significantly improve the speed of convergence if one transforms these operators to a more regular form [32]. Specifically, for a given nearly singular operator H X we search for another, more regular operator H XR and an additional operator Q X , which satisfy the following equation where {. , .} denotes the anticommutator. It is obvious that H X = H XR , as long as the expectation value is evaluated with the eigenfunctions of the Hamiltonian H 0 . In practice, it is usually possible to find such a pair of operators H XR , Q X that the most singular part of H X is absorbed in the anticommutator. The additional operator Q X is generally a combination of Z/r 1 , Z/r 2 , and 1/r, with the coefficients in front of these terms determined by requiring the cancellation of all Dirac-δ-like contributions. Specifically, we find the following regularized form of the operator H A (without the nuclear recoil) [10] where V = −Z/r 1 − Z/r 2 + 1/r. The operator ∇ 2 1 ∇ 2 2 in the above formula is not self-adjoint and requires an explicit definition. Its action on a trial function φ on the right should be understood as plain differentiation (omitting δ 3 (r)); no differentiation by parts is allowed in the matrix element. It can be checked that the operators H A and H AR satisfy the following equation where Formulas with the finite nuclear mass are analogous but more lengthy; they are given by Eqs. (62)-(67) of Ref. [33]. Table III presents numerical results for the leading relativistic correction to the 2 3 P centroid energy, performed with different basis sets. We observe that, for the same basis size, the number of correct digits for the matrix element is half as much as for the nonrelativistic energy.

B. Fine structure
The fine structure of energy levels is induced by spindependent operators. The spin-dependent part of the Breit Hamiltonian is conveniently written as a sum of three operators with different spin structure, with where κ = α /2π + O(α 2 ) is the anomalous magnetic moment correction and σ a is the vector of Pauli matrices acting on a'th electron. We note that the operators H B , H C , and H D contain radiative corrections in form of the electron anomalous magnetic moment. In this way we account for the complete QED effects of order α 5 m to the fine structure. It should be mentioned that the matrix element of H C is nonzero only if the operator is sandwiched between wave functions with different spin values. Therefore, any symmetrical matrix element of H C vanishes, and this operator does not contribute in the leading order of perturbation theory. We note, however, that H C contributes to the second-order perturbation corrections (in the order α 6 m).
In order to perform the spin-angular reduction in the matrix elements of H fs , it is convenient to introduce spatial operators Q B , Q C , and Q D , explicitly separating the spatial and the spin degrees of freedom, Using Eqs. (31), (38)- (41) and performing traces of the spin operators, we express all matrix elements in terms of spatial radial integrals. For the 3 P states, we obtain where for J = 0, 1, and 2, respectively. For the 3 D states, an analogous calculation yields where for J = 1, 2, and 3, respectively.

VII. LEADING QED CORRECTION
The leading QED contribution is of the order α 5 m. For the fine structure, this contribution is already accounted for by the electron anomalous magnetic moment terms in the Breit Hamiltonian, as given by Eqs. (74)-(76). So, we need to examine only the centroid energy.
The spin-independent mα 5 Hamiltonian representing the leading QED effects was derived in the 1950s by Araki and Sucher [34,35] where ln k 0 is the so-called Bethe logarithm defined as and 1/r 3 ǫ is the regularized 1/r 3 operator (distribution) defined by its matrix elements with an arbitrary smooth function f ( r) as +4 π δ 3 (r) (γ E + ln ǫ) . (90) The nuclear recoil correction to the leading QED contribution consists of two parts, where δ M H is defined by Eq. (62), and δ M H (5) is the recoil addition to the α 5 m Hamiltonian given by [36] Here, δ M ln k 0 is the correction to the Bethe logarithm ln k 0 induced by the nonrelativistic kinetic energy operator P 2 /2, and 1/r 3 a,ǫ is the regularized 1/r 3 a operator defined analogously to Eq. (90).
The recoil correction to the Bethe logarithm δ M ln k 0 is often separated into the mass-scaling and masspolarization parts, where δ p1p2 denotes the perturbation due to the mass polarization operator p 1 · p 2 . The corresponding separation for the 1/r 3 ǫ matrix element reads From the computational point of view, the numerical evaluation of the QED effects involves two new features, as compared to the relativistic correction: matrix elements of the singular operators 1/r 3 and 1/r 3 a and the Bethe logarithm. Calculation of expectation values of singular operators with exponential basis functions is examined in Appendix A; it does not present any computational difficulties. On the contrary, the computation of the Bethe logarithm is rather nontrivial; it is examined in the next section.

A. Bethe logarithm
There are two different approaches developed for the calculation of the Bethe logarithm in few-electron atoms. The first one starts with the definition (89) and uses the basis-set representation of the Hamiltonian as a sum of the spectrum of the eigenfunctions. The difficulty is that the sum in the numerator is nearly diverging because the dominant contribution comes from the high-energy continuum states of the spectrum. This problem is solved by using a basis set whose spectrum of pseudostates spans a huge range of energies [37].
An alternative approach was first introduced by C. Schwartz [23] and further developed by V. Korobov [38][39][40]. Within this method, the Bethe logarithm ln k 0 is represented as an integral over the momentum of the virtual photon, with subtracting the ultraviolet asymptotics and performing the limit, where D = 2πZ δ 3 (r 1 ) + δ 3 (r 2 ) , ∇ ≡ ∇ 1 + ∇ 2 , and The asymptotic expansion of J(k) for large k reads Splitting the integration interval (0, Λ) into two parts (0, K) and (K, Λ), where K is an arbitrary cutoff parameter, we can rewrite Eq. (95) as The above expression is finite, does not depend on K, and is suitable for a numerical evaluation. We now address the angular reduction in the secondorder matrix element J(k) given by Eq. (96). It is performed in several steps. First, we represent the gradient acting on the reference-state wave function ∇ j ψ i1..iL as a sum of irreducible Cartesian tensors, as described in Appendix B. For example, the gradient acting on a Pstate wave function ∇ j ψ i is represented as a sum of the L = 0, L = 1, and L = 2 irreducible Cartesian tensors, which induce, correspondingly, the L = 0, L = 1, and L = 2 angular-momentum contributions from the resolvent. The second-order matrix element of an irreducible tensor Φ i1..iL is transformed as where Φ i1..iJ is the solution of the inhomogeneous Schrödinger equation Inserting the explicit representation of Φ as a sum over the spectrum, we obtain An alternative way to arrive at this expression is to observe that the scalar product Φ|ψ includes an integration over the continuous and a summation over the discreet variables, namely Φ|ψ ≡ Φ i1..iL |ψ i1..iL = i1..iL d 3n r Φ i1..iL * (r) ψ i1..iL (r). The advantage of the integral representation of the Bethe logarithm is that J(k) has a form of the symmetric second-order perturbation correction and thus obeys the variational principle. We therefore can variationally optimize the basis-set representation of the resolvent For lower values of k, the basis can be variationally optimized if one fixes pre-optimized parameters of the more deeply bound states with E n < E 0 .
Our numerical procedure was performed in two steps. First, we optimized the basis for several different scales of the photon momentum, k = 10 i , with typical values of i = 1, .., 4. After that, the computation of the function J(k) was performed with a basis obtained by merging together the optimized sets for the two closest k i points, thus essentially doubling the size of the basis. In the second step, we perform the integration over k. The integral over (0, K) (with the typical choice of K = 10) was calculated analytically, after the full diagonalization of the Hamiltonian matrix. The remaining interval was split into two parts, (K, K 2 ) and (K 2 , ∞), with the typical choice of K 2 = 10 4 . The integral over the former was performed with help of Gauss-Legendre quadratures, after the change of variables t = 1/k 2 . The remaining part of the integral was calculated analytically, after fitting numerical values of J(k) to the known form of the asymptotic expansion, where pol(x) denotes a polynomial of x. The first terms of this expansion are given by Eq. (97), whereas the higher-order coefficients are obtained by fitting. Calculations of the Bethe logarithm for the finite nuclear mass can be performed analogously to the above, or by perturbation theory. The numerical procedure for evaluation of the recoil correction to the Bethe logarithm by perturbation theory is described in Appendix A of Ref. [41]. Table IV presents a comparison of different calculations of the Bethe logarithm for the 2 3 P state of helium. The most accurate results for the ground and excited states of helium are obtained by Korobov in Ref. [40]. Results for He-like ions can be found in Refs. [37,41].

VIII. α 6 m QED EFFECTS
The α 6 m QED corrections to energy levels in atoms are represented by the sum of the expectation value of the effective α 6 m Hamiltonian H (6) and the second-order perturbation correction induced by the Breit Hamiltonian, where H We note that in order to avoid admixture of higher-order contributions in E (6) , we have to retain only the α 4 m part in the definition of the Breit Hamiltonian, i.e., to set the magnetic moment anomaly κ → 0 in the definitions (74)-(76). This is indicated by the superscript "4" in the corresponding operators. Formulas for the effective α 6 m Hamiltonian H (6) are rather lengthy and will not be reproduced here. In the case of fine structure, they were first obtained by Douglas and Kroll in 1974 [42] and later re-derived in Refs. [43,44]. For the energy centroid, the situation is greatly complicated because of the appearance of numerous diverging operators. The corresponding derivation was accomplished by one of us (K.P.), in Ref. [9] for the triplet states and in Ref. [10] for the singlet states of helium. The complete formulas suitable for numerical evaluation can be found in Ref. [25].
The nuclear recoil α 6 m correction has the same structure as the non-recoil one, but the expressions for the operators are much more complicated. This correction was calculated in Ref. [33] for the triplet states and in Ref. [45] for the singlet states of helium.

A. Second-order terms
We now discuss the evaluation of the second-order contributions, represented by the second term in Eq. (103). Such corrections were first calculated for the fine structure by Hambro [46] and by Lewis and Serafino [7]. Later, the fine-structure calculations were greatly improved in Refs. [47][48][49]. For the centroid energies, the second-order corrections were calculated in Refs. [10,11] for the 2S and 2P states and in Refs. [25,50] for the nD states of helium.
It is convenient to rewrite Eq. (103), expressing the second-order perturbation correction more explicitly, (105) We note that the non-symmetrical second-order corrections (the last two terms in the above equation) vanish for the centroid energy, but contribute to the fine structure. The second-order perturbative corrections are calculated as follows. In the first step, we perform traces over the spin degrees of freedom in the matrix elements. Then we decompose the product of a tensor operator Q and the reference-state wave function ψ into the irreducible tensor partsψ, as described in Appendix B. In the last step we calculate the second-order matrix elements induced by the irreducible partsψ as (see Eq. (101)) The numerical evaluation of symmetrical second-order contributions was carried out with the variational optimization of the nonlinear parameters of the basis set for the resolvent 1/(E 0 − H 0 ). Convergence of numerical results is often rather slow, especially for contributions with H AR . This is associated with the fact that the effective wave function |δψ = 1/(E 0 − H 0 ) ′ |H AR has an integrable singularity at r a → 0. In order to represent such wave functions with the exponential basis, very large (both positive and negative) exponents are required. In order to effectively span large regions of parameters, we used non-uniform distributions of the nonlinear parameters. E.g., for the nonlinear parameters α i we used the distributions of the kind [9] with a = 2 and 3, where the variable t i has a uniform quasirandom distribution over the interval (0, 1) and the variables A 1,2 are subjects of variational optimization. An example of the convergence study of the second-order correction H AR 1 (E0−H0) ′ H AR is given in Table V. Numerical evaluation of non-symmetrical second-order contributions was carried out with basis sets, optimized for the corresponding symmetrical corrections.

IX. α 7 m QED EFFECTS
The α 7 m QED correction to energy levels in atoms is given [12] by the sum of the relativistic correction to the Bethe logarithm E L , the expectation value of the effective α 7 m Hamiltonian H (7) , and the perturbation of the α 5 m QED operator by the Breit Hamiltonian, The regularized effective α 5 m Hamiltonian is [17] H (5) where H R is non-Hermitian and is assumed to act on a ket trial function φ on the right.
The relativistic correction to the Bethe logarithm is rather complicated. We will not discuss its calculation here, but direct the reader to original studies. This correction was first calculated for the fine structure of the 2 3 P state; the corresponding calculations for helium and helium-like ions were performed in Refs. [12][13][14]. In our recent investigation [15] we performed a calculation for the energy centroid of the 2 3 S and 2 3 P states. For singlet states of helium, this correction has never been calculated so far.
The derivation of the effective α 7 m Hamiltonian H (7) for helium is an extremely difficult problem. It was first accomplished by one of us (K.P.) for the fine structure in Refs. [12,13]. Recently, we performed [16,17] the derivation of H (7) for triplet states of helium and calculated [27] the corresponding correction to the energies of the 2 3 S and 2 3 P states. For singlet states, the effective α 7 m Hamiltonian is unknown.
From the computational point of view, the main difficulty of the evaluation of the α 7 m correction is the calculation of the Bethe-logarithm contribution E L . The computational scheme is similar to that for the plain Bethe logarithm and is described in Ref. [15]. Conversely, the computation of the expectation value of H (7) and the second-order corrections is very similar to the calculation of the α 6 m corrections.

X. OTHER EFFECTS
The finite nuclear size correction is given by (in relativistic units) where R is the root-mean-square nuclear charge radius, and the expectation value of the Dirac δ functions is assumed to include the finite-nuclear-mass correction induced by δ M H. The higher-order QED effects are approximated on the basis of known results for hydrogenic atoms. Specifically, the hydrogenic one-loop and two-loop corrections for the 2s state of He + are given by [51] An approximation for the higher-order α 8 m QED correction to the ionization energies of the helium atom is obtained from the corresponding hydrogenic 2s contribution by

XI. COMPARISON OF THEORY AND EXPERIMENT
In this section we summarize numerical results of QED calculations of energy levels in 4 He and compare theoretical predictions with available experimental results. Table VI presents such a comparison for transitions between states with the principal quantum number n = 2. We note that our present theoretical uncertainty for the 2 3 S -2 1 S transition is increased as compared to our previous work [28]. The reason is an accidental cancelation of the estimated α 7 m term between the 2 3 S and 2 1 S states in Ref. [28]. Now the α 7 m correction is calculated for the 2 3 S state and the theoretical uncertainty is defined by the 2 1 S state only. Table VI shows good agreement of theory and experiment for the singlet-singlet and triplettriplet transitions but some tension for the singlet-triplet transitions. Specifically, we note a 2.3 σ deviation from the experimental result [53] for the 2 3 S-2 1 P transition (with σ denoting the standard deviation).
Of particular importance is the agreement observed for the 2 3 P -2 3 S transition, because in this case two triplet states are involved, for which the theoretical accuracy is the highest. Theoretical calculations of energies for the 2 3 S and 2 3 P states [17] are complete through order α 7 m, with resulting theoretical uncertainty below 100 kHz, whereas for the 2 1 S and 2 1 P states the theory [28] is complete up to order α 6 m only and the theoretical accuracy is on the level of 1 MHz. For the D states, theoretical calculations [25,50] are also complete at order α 6 m, but the absolute theoretical precision is much higher since the QED effects are smaller. In general, we conclude that for the intrashell n = 2 transitions there is good agreement for transitions between the states with the same spin multiplicity and some tension for the states of different spin multiplicity.
The situation becomes even more strained when we consider ionization energies and transitions involving states with different n's. The corresponding comparison is presented in Table VII. We immediately notice that all differences between theory and experiment are of the same sign and that most of them are outside of the theoretical error bars. The largest discrepancies are found for the 2 3 S 1 -3 3 D 1 and the 2 3 P 0 -3 3 D 1 transition, of 6 and 12 σ, correspondingly. These transitions involve the triplet states, for which theoretical uncertainties are the smallest, so that 0.5 MHz differences from the experimental values lead to large relative deviations.
The comparison in Tables VI and VII suggests that there might be a contribution missing in theoretical calculations of energy levels, which weakly depends on L but strongly depends on the principal quantum number n (the latter is natural because the 1/n 3 scaling is typical for all QED effects). This conjecture was put forward in Ref. [50] and since then strengthened by subsequent calculations and measurements. Such a missing contribution most likely originates from the α 6 m or α 7 m QED corrections because all other theoretical effects are crosschecked against independent calculations [5]. Table VIII presents the comparison of theoretical and experimental results for the fine-structure intervals of the 2 3 P state in 4 He. Theoretical predictions for these intervals are of greater accuracy than for other intervals of the n = 2 manifold. This is both due to the fact that the theory of these intervals [14,57] is complete at the order α 7 m and due to the smallness of QED effects. We observe a generally good agreement between theory and experiment for the fine-structure intervals. The only tension is a 1.4 σ deviation for the P 1,2 interval measured in Ref. [3]. We note that all pre-2010 experimental results were to a greater or lesser degree influenced by unaccounted quantum-interference effects and were reevaluated in Refs. [74,75].
Summarizing, we have reviewed a large amount of work accomplished during the last decades in calculations of QED effects in the atomic structure of the helium atom. The leading-order α 5 m QED effects are nowadays well established by independent calculations and tested by comparison with numerous experiments. However, recent calculations of higher-order α 6 m and α 7 m QED effects revealed some small but systematic deviations from high-precision experimental transition energies. Having in mind the importance of the helium spectroscopy for determination of nuclear properties and fundamental constants, we conclude that further theoretical and experimental efforts are needed in order to find the reasons behind the observed discrepancies. In calculations of the Bethe logarithm and the secondorder perturbation corrections, we encounter a problem of decomposition of products of irreducible Cartesian tensors into the irreducible parts. In this section we collect formulas required for such decompositions. The product of two vectors is represented as a sum of a symmetric and traceless second-rank tensor, a vector, and a scalar, The product of a vector and a symmetric and traceless second-rank tensor is decomposed as where This identity can be verified by contracting Eq. (B2) with δ ij and ǫ ijk . It can be easily extended to the higher-rank tensors Q. Finally, we present the decomposition of the product of two symmetric and traceless tensors P ij and Q kl , required for calculations of second-order corrections involving D-states, P ij Q kl = (P ij Q kl ) (4) + ǫ ika T jal + ǫ jka T ial + ǫ ila T jak + ǫ jla T iak + δ ik T jl + δ il T jk + δ jk T il + δ jl T ik