Edit Mátyus,Dávid Ferenc,Péter Jeszenszki,Ádám Margócsy
The Bethe–Salpeter QED Wave Equation for Bound-State Computations of Atoms and Molecules
Bethe−Salpeter equation,no-pair Dirac−Coulomb−Breit equation,wave equation,instantaneous interaction,retardation,pair corrections,variational relativistic computations,explicitly correlated Gaussian
Interactions in atomic and molecular systems are dominated by electromagnetic forces and the theoretical framework must be in the quantum regime. The physical theory for the combination of quantum mechanics and electromagnetism, quantum electrodynamics has been “established” by the mid-twentieth century, primarily as a scattering theory. To describe atoms and molecules, it is important to consider bound states. In the nonrelativistic quantum mechanics framework, bound states can be efficiently computed using robust and general methodologies with systematic approximations developed for solving wave equations. With the sight of the development of a computational quantum electrodynamics framework for atomic and molecular matter, the field theoretic Bethe–Salpeter wave equation expressed in space–time coordinates, its exact equal-time variant, and emergence of a relativistic wave equation, is reviewed. A computational framework, with initial applications and future challenges in relation with precision spectroscopy, is also highlighted.
The Bethe–Salpeter QED Wave Equation for Bound-State Computations of Atoms and Molecules

Interactions in atomic and molecular systems are dominated by electromagnetic forces and the theoretical framework must be in the quantum regime. The physical theory for the combination of quantum mechanics and electromagnetism, quantum electrodynamics has been “established” by the mid-twentieth century, primarily as a scattering theory. To describe atoms and molecules, it is important to consider bound states. In the nonrelativistic quantum mechanics framework, bound states can be efficiently computed using robust and general methodologies with systematic approximations developed for solving wave equations. With the sight of the development of a computational quantum electrodynamics framework for atomic and molecular matter, the field theoretic Bethe–Salpeter wave equation expressed in space–time coordinates, its exact equal-time variant, and emergence of a relativistic wave equation, is reviewed. A computational framework, with initial applications and future challenges in relation with precision spectroscopy, is also highlighted.

Dirac’s one-electron space–time equation was an ingenious departure from Schrödinger’s time-dependent wave equation to have a Lorentz covariant description, but it was strange due to the introduction of hole theory that seemed a bit artificial. In a recorded discussion from 1982, Dirac modestly admitted to Hund that for him, it took a year, perhaps two, to understand the role of the negative-energy states. Breit attempted a two-particle generalization of Dirac’s one-electron theory in a series of papers between 1928 and 1931, by adopting Darwin’s 1920 calculation of the classical electromagnetic interaction for two moving charges and the quantum mechanical velocity operator obtained with Dirac’s formalism. Already from the beginning, it was apparent that the Breit equation was not Lorentz covariant. Nevertheless, Breit used this “quasi-relativistic” equation in a perturbation theory approach imposing the Pauli approximation to the four-(16-)component wave function. Good agreement with experiment was obtained after discarding a term from the result “by hand”. This procedure was later explained based on Dirac’s hole theory by Brown and Ravenhall. The problem, called Brown–Ravenhall (BR) disease, related to the artificial coupling of the positive- and negative-energy states of Dirac’s theory when naïvely applied to two-particle systems, survived also in the modern literature and it is commonly used to explain the failure of the two-(many-)particle Breit equation. A recent numerical study demonstrated that bound states of helium-like two-electron systems represented by the Breit equation have (unphysical) finite lifetimes (on the order of α3Eh, where α is the fine-structure constant). It has been also discussed that there was no BR dissolution problem for isolated two-particle systems, like positronium, when modeled with the Breit equation, and it was numerically demonstrated by a finite-element computation that the energy levels of the two-particle Breit equation (with Coulomb interactions, no external fields) are “stable”. Nevertheless, the two-particle Breit equation is still incorrect, or in other words, “correct only up to order α2Eh”. A consistent and Lorentz covariant many-particle theory was put forward by the development of quantum electrodynamics (QED). As a natural continuation of Feynman’s two papers in 1949 on the reinterpretation of the mathematical solutions of the Dirac equation and the development of the space–time approach to quantum electrodynamics, Salpeter and Bethe in 1951 published a Lorentz-covariant wave equation for two interacting particles (with stating that generalization to more than two particles is straightforward). It is interesting to note that the same equation was written at the end of a paper without explanation by Nambu already in 1950, and it was formulated by Schwinger and also by Gell-Mann and Low during 1951. In 1952, Salpeter used this equation for the hydrogen atom in combination with perturbation theory and an instantaneous interaction kernel. Probably, this was the first formulation of the exact equal-time equation for two-particle systems. Salpeter reported results for the hydrogenic case (of one heavy and one light particle, M ≫ m) up to the α3(m/M)Eh order, and he stated that the calculation can be generalized to any masses. In 1954, Fulton and Martin calculated the energy levels for an arbitrary two-fermion system, such as positronium, up to α3Eh order. In 1958, Sucher’s PhD thesis represented another important step forward using the formalism and extending Salpeter’s work to a two-electron system in an external Coulomb field, for the example of the helium atom. Sucher’s final α3Eh-order correction formulas were identical with those reported by Araki a year earlier, but we can build on the fundamental ideas explained in Sucher’s work for further developments. In 1974, Douglas and Kroll started their paper on the α4Eh-order corrections to the fine-structure splitting of helium with a good review of Sucher’s work by extending the formalism with explicitly writing also the radiative terms in the wave equation (Sucher only highlighted the steps at the end of his work). Then, in 1989, Adkins elaborated this direction for positronium, still relying on a perturbative expansion with respect to the nonrelativistic reference for practical calculations. During the 1990s, Zhang worked on higher-order corrections to the fine-structure splitting (α5Eh) and energy levels (α4Eh) of helium. Pachucki initiated a different approach starting from the late 1990s. This approach is based on performing a Foldy–Wouthuysen transformation of the Dirac operator in the Langrangian density—thereby linking the formalism to the nonrelativistic theory from the outset, and then, collecting corrections to the poles of the equal-time Green function to the required α order. In 2006, Pachucki reported the complete α4Eh-order corrections to the energy levels of singlet helium, thereby extending the 1974 work of Douglas and Kroll valid only for triplet states, as well as work from Yelkhovsky and computations from Korobov and Yelkhovsky in 2001 for α4Eh-order corrections of singlet helium. In 2016, the complete α4Eh-order corrections derived by Pachucki were used for the ground electronic state of the H2 molecule with fixed protons. Most recently, the ‘Foldy–Wouthuysen–Pachucki’ approach has been used to derive α5Eh-order contributions for triplet states of helium. In contrast to using a nonrelativistic reference (as in all previous work), we aim for a relativistic QED approach, in which some (well-defined, many-particle) relativistic wave equation is first solved to high precision and used as a reference for computing “QED” (retardation, pair, and radiative) corrections up to a required accuracy. Such an approach appears to be feasible along the lines formally started by Bethe, Salpeter, Sucher, Douglas, and Kroll. These authors performed calculations by hand, so in the end, they had to rely on approximations based on the nonrelativistic formalism. Nowadays, we can use the power of modern computers to first numerically solve a many-particle relativistic wave equation, and then, compute corrections to the relativistic energy. It is also necessary to add that there have been several articles on understanding and solving the original, space–time Bethe–Salpeter (BS) equation. For atomic and molecular computations, the exact equal-time form of the BS equation, as introduced by Salpeter and Sucher, appears to be more promising. In this approach, a two-particle relativistic Hamiltonian and corresponding wave equation emerges, for which numerical strategies for solving wave equations, including the variational method, can be used. During the 1980s, Sucher published review articles about the (formal) connection of the equal-time BS wave equation with the relativistic quantum chemistry framework and computational methodologies, e.g., refs (46 and 47). In this context, it is necessary to mention the excellent book of Lindgren who further developed these ideas for orbital-based many-body applications in chemistry. We consider the renormalized, “mixed gauge”, two-particle Bethe–Salpeter equation as the starting point for a theoretical framework of atoms and molecules and with relevance for spectroscopic applications. This theoretical framework is reviewed in the first part of the paper by relying on work by Sucher, Douglas and Kroll, as well as Salpeter. The second part of the article highlights our recent work, algorithmic details of a computer implementation and numerical results for two spin-1/2 particles with and without a fixed, external Coulomb field, i.e., with relevance for relativistic Born–Oppenheimer (BO) as well as relativistic pre-Born–Oppenheimer (pre-BO) computations. Although in the present review, we focus on the theory and a numerical procedure for two-particle systems, we mention Sucher’s series of papers from the 1980s implying a possible generalization and Broyles’ work from 1987 about presenting a line of thoughts connecting field theory and an N-particle no-pair Dirac–Coulomb–Breit wave equation. Regarding a relativistic QED approach, we also mention the quasi-potential method, which originates from Logunov, Tavkhelidze, and Faustov working during the 1960–70s, and the corresponding two-time (equal-time) Green function idea developed by Shabaev. Comparison of the Salpeter–Sucher approach with the quasi-potential method is left for future work.

The Dirac equation for a particle of mass m1 and x1 = (r1, t1) space–time coordinates iswith the free-particle Hamiltonianand the α1 and β1 Dirac matrices. Feynman pointed out in 1949 that instead of working with the Hamiltonian equation, it is often more convenient to use the corresponding Green function or propagatorFor a Dirac particle in an external scalar field, Φ1with stands for the charge number of the active particle and e is the elementary charge. Simple calculation shows that the G1 propagator can be obtained from the G1(0) free-particle propagator through the integral equation (corresponding to subsequent interaction events of the particle with the external field) asAccording to Feynman’s combination of the electronic and positronic theory in a consistent manner, the propagator is expressed with the eigenvalues and eigenfunctions of the Dirac Hamiltonian as the sum over positive-energy (electronic) states moving forward in time, and the negative sum over negative-energy (positronic) states moving backward in time. Feynman defined the free-particle propagator this way corresponding to eq 3. The arguments can be taken over for a particle in an external field, eq 4, which is known as the “Furry picture”, and the propagator iswhere ϕ̅n(1) = ϕn(1)*β1 is the Dirac adjoint. The G1(0)(x1, x1′) free-particle propagator is recovered for eigenvalues and eigenfunctions of the Dirac equation with Φ1 = 0. Regarding the external field in the present work, only the scalar potential due to the Coulomb field of the fixed nuclei will be relevant, e.g., for helium-like systems with the nucleus fixed at the origin and with Z nuclear charge number, the interaction energy iswhere α = e2/(4π) is the fine-structure constant in natural units (ℏ = c = ϵ0 = 1). To describe a two-particle system, we can consider the G(x1, x2; x1′, x2′) two-particle propagator or amplitude which describes that particles 1 and 2 get from x1′, x2′ to x1, x2 space–time points. For noninteracting particles, the two-particle propagator is the simple product of the one-particle propagators, G1(x1, x1′)G2(x2, x2′). For interacting two-particle systems Salpeter and Bethe, following Feynman, devised an integral equation, called Bethe–Salpeter (BS) equation:where K is the interaction function. In particular, K must contain only the so-called “irreducible” interactions, since the corresponding consecutive, so-called “reducible”, interactions are already included, “iterated” to all orders by the integral equation. The simplest interaction function, K(1), corresponds to the single photon exchange (see also section 2.3) with γiμ = βi(αi, 1)where DFμν is the photon propagator, which takes a simple, manifestly covariant form in Feynman gauge:To describe the interparticle interaction in atoms and molecules, it is more convenient to use the Coulomb gauge, in which the interaction is the sum of the Coulomb (C, the dominant part) and the transverse (T) contributions,withwherecorresponds to the transverse components of α perpendicular to k. (In eq 13, we highlighted the well-known coordinate–space form of the Coulomb interaction, the primarily important interaction term in quantum chemistry.) If radiative corrections are accounted for, the one-electron (one-particle) propagator is replaced byor equivalentlywhere Σ1 is the sum of the one-electron self-energy contributions. Douglas and Kroll, following the last chapter of Sucher’s work, formulated the two-electron equations by formally including all radiative corrections. This formulation also relies on the work of Mathews and Salam who explained that the Bethe–Salpeter equation can be renormalized with the replacement of G1′, G2′, γ1μ, γ2μ, and the DF photon propagator byWe note that Mathews and Salam mostly formulated their renormalization approach based on series expansion, whereas Källen and Lehmann defined renormalization terms without the use of power series expansion in the interaction constant. Karplus and Kroll and Jauch and Rohrlich carried out explicit calculations for the Σ1*, Π*, Λ1μ* renormalized electron self-energy, photon self-energy, and vertex correction operators to order α for the case of no external potentials. For renormalization, it is necessary to work in the Feynman gauge, eq 11. At the same time, binding of the particles in atomic and molecular systems is dominated by the Coulomb interaction, eq 13, which can be identified by writing the interaction operators in the Coulomb gauge. According to Sucher’s arguments (following the field theoretical derivation of the BS equation by Gell–Mann and Low), it is valid to perform renormalization of the radiative terms in the Feynman gauge, and then, use the resulting expressions for the interacting problem, written in the Coulomb gauge. This special procedure is known as the mixed-gauge representation.

Equation 9 can be rewritten for the wave function of a bound state (e.g., Chapter 6 ref (65) or Chapter 12 of ref (66)), formally including now also the radiative effects, asor in shortwhich, using eq 17, can be rearranged to (we note the missing imaginary unit in ref (24))where the full “interaction” kernel, containing also the radiative corrections, was defined asFrom rearrangement of the operator form of eq 4, G1–1 = −iβ1[i∂/∂t1 – H1] and using β1β1 = β2β2 = 1, we obtainwhich is a (space–time) wave equation which accounts for (nonradiative) interactions and radiative corrections on an equal footing. Since eq 25 includes the “own” time for both particles, but the U external interaction (if any) is time independent in our frame of reference, we can write the two-particle wave function aswhere the average (“total”) time and relative time was introduced asIt is important to note that E is the total energy of the two-particle system and it corresponds to the T total time, eq 26. Similarly to T and t, we defineThen, we obtain the following equation for the Ψ(r1, r2, t) space- and relative-time wave functionwith the interaction kernel depending only on the relative time variableswhere it is exploited that the external field is time independent, i.e., K′(x1, x2; x1′, x2′) depends on T and T′ only through the T′ – T difference, and T represents only a constant shift for the T′ integration variable. Both Sucher and Douglas and Kroll continued the calculation in momentum space, and we follow this line of thought. The r1, r2 space coordinates of the two particles and the t relative time are replaced with the p1, p2 momenta and the ε relative energy. The relative-time and relative-energy wave functions are connected by the seven-dimensional Fourier transformationwhile the interaction kernel in momentum space is defined asand it acts as an integral operator,Then, eq 29 can be rewritten asandwithand their inverse define the one-particle propagators (for the E total and ε relative energy), which will be used in later sections:with the four-vector variables p1 = (p1, ε) and p2 = (p2, –ε). (and similarly ) is the momentum-space form of the one-particle Dirac Hamiltonian, eq 5. In this representation, the interaction operators are integral operators, as it was indicated in eq 33 for . (and ) labels the external-field Coulomb operator. For the example of a single nucleus fixed at the origin, eq 8:

According to eq 24, the full interaction kernel contains contributions both from “interparticle” interactions and from radiative contributions:In what follows, we focus on the construction of the interparticle kernel, which is obtained as the sum of operators corresponding to irreducible “diagrams”. Action of the interaction kernel for a single-photon exchange (written in the Coulomb gauge), eqs 12–14, on some f(p1, p2, ε) two-particle function depending also on the ε relative energy can be written in the momentum-space representation as A more compact operator form of and is obtained by using k = (k, ω) momentum-energy translation operators. The one-particle translation operators arethe two-particle translation operator isand for later convenience, we also define the notationThen, the Coulomb and transverse parts of the one-photon exchange, eqs 43 and 44, can be written aswithandwithIt is interesting to note that the Coulomb interaction carries only a trivial shift in the relative-energy dependence, and this corresponds to saying that the Coulomb interaction acts through momentum transfer and the interaction is instantaneous. At the same time, the transverse part has a nontrivial relative-energy dependence, and this is related to the finite propagation speed of the overall interaction (retardation). At the same time, the retardation contribution to the transverse part is small and it is convenient to separate the instantaneous part, which is called the Breit interaction:withThe remainder, i.e., difference of the transverse and the Breit interactions, is the retarded part, which we label aswhile the instantaneous contributions (Coulomb–Breit) can be handled “together”: To write down the mathematical expression for more complicated interactions including multiple (Coulomb and/or transverse) photons (e.g., Figure 1), Sucher derived and summarized the following simple rules, which we call Sucher’s (interaction) rules. First, (a) label each interaction line with four vectors, k, k′, k″, etc. with assigning each line a specific sense (convenient to choose the same for all lines), e.g., from 2 to 1; (b) label the final parts of the world lines of the fermions 1 and 2 with p1 = (p1, ε) and p2 = (p2, –ε); (c) label all remaining electron lines with conserving the four-momentum. Second, for a fully labeled diagram, can be formulated by writing 1. α̃1i or α̃2i for a transverse interaction vertex; 2. a factor S1(p1 – k) for an intermediate electron line labeled with p1 – k on the path of 1 and a factor S2(p2 + k′) for an intermediate electron line labeled with p2 + k′ on the path of 2; while writing down the factors, it is necessary to preserve the order of events along a world line, i.e., factors for “later” events along a world line stand to the left of factors corresponding to “earlier” events; 3. to the right of these expressions a factor κC(k) for a Coulomb interaction line labeled with k = (k, ω) and a factor κT(k′) for a transverse interaction line with k′ = (k′, ω′); 4. in addition to each κC and κT, an η(k)/(−2πi) factor appears, if the interaction is from 2 to 1 (or an η(−k)/(−2πi) factor if the interaction is from 1 to 2). It is also useful to note that the effect of the η(k) = η(k, ω) four-momentum translation on the one-particle propagators, eqs 38 and 39, iswhere the k-translation of the one-particle Hamiltonians gives Sucher calculated corrections to the energy up to order α3Eh, Douglas and Kroll calculated the fine-structure splitting up to order α4Eh, and they have included the following interactions:It is necessary to compile “by hand” only the irreducible interactions, and all reducible diagrams are automatically included in the solution of the BS equation.

Let us exploit the fact that, in atoms and molecules, the dominant part of the interaction is instantaneous (Coulomb or Coulomb–Breit), so, it is convenient to write the kernel as the sum of a instantaneous part and the “rest”The instantaneous part, , induces only a trivial shift for the ε relative energy, and hence the effect of the relative energy can be integrated outwhere is only a momentum shift integral operatorwithIn eq 63, it is also important to note the emergence of the equal-time wave function,Next, we rearrange eq 34 and separate the instantaneous part of the interaction asBy integrating both sides with respect to the relative energy and using eqs 63 and 66, we obtainwhere the operator identity was used in the last stepNext, we define the one-particle positive- and negative-energy projection operators for particles i = 1 and 2 bywhich, at this point, contains a purely formal definition for the one-particle Hamiltonian absolute value operatorwhich also means thatIn short, we can also writeIf there is no external field, e.g., pre-BO description of two spin-1/2 particles, then, reduces to the free-particle projector:since for Ui = 0, and the eigenvalues of the Hamiltonian absolute value operator are . Using the and notation, we can write the and propagators in asand similarlyaccording to Feynman’s prescription of adding a complex number with a small negative imaginary value to the mass and the limit is taken from the positive side for a consistent electron-positron theory (here, the energy replaces Feynman’s mass and δ > 0 with δ → 0+). The first term in eq 68 contains a relative-energy integral, but and Φ are independent of ε, so we only need to calculatewith the two-particle projectors, (σ, σ′ = + or – ). To evaluate these ε integrals, we use Cauchy’s residue theorem:where the summation goes through the poles of f within the domain surrounded by the simple closed curve γ. We can choose the positive γ contour (Oγ counterclockwise, sgnγ = +1), but identical results are obtained from using the negative γ′ contour (Oγ′ clockwise, sgnγ′ = −1). Since this is an important step of the calculation, we proceed term by term with the evaluation of eq 77.The result of this short calculation can be summarized as operator identitieswhere the second identity holds in general, the first is valid only for commuting A and B operators. All in all, we obtainwhere the short notation is introduced,Using this result, we can rewrite eq 68and finally obtain, the exact, equal-time Bethe–Salpeter (eBS) equationwithEquation 88 is the central equation to our work. It is obtained by equivalent mathematical manipulations from the original space-time Bethe–Salpeter equation, eq 9, it is a homogeneous, linear equation for the equal-time wave function, Φ, which depends only on the momenta (or coordinates) of the two fermions. At the same time, the exact equal-time equation is a nonlinear eigenvalue equation for the E energy, since the term also depends on E (through ). We can arrive at a useful initial description of atoms and molecules, by first neglecting , and starting with the solution of the positive-energy projected or no-pair Dirac–Coulomb(−Breit) equationIt is important to note that in the present derivation, the projector is defined according to eqs 70 and 74, and it is connected to the emergence of the no-pair two-particle Dirac Hamiltonian, eqs 77–86. Variants of the no-pair DC(B) equation are commonly used in relativistic quantum chemistry. Sucher analyzed the connection to relativistic quantum chemistry methodologies, in which the Dirac–Hartree–Fock projector is a popular (and natural) choice, and came to the conclusion that the use of that projector is also valid, but then, during the evaluation of the corrections, one has to correct for the difference (which may be complicated). During our work, we stick to the original definition, eq 70 for two particles in an external field and eq 74 for an isolated two-fermion system. Corresponding numerical results (for helium- and for positronium-like systems) are reviewed in section 5. During the present calculation, which follows closely the work by Salpeter, Sucher, and Douglas and Kroll, it was critical to retain the relative energy between the particles. Integration for the relative energy resulted in the emergence of the no-pair, two-electron Dirac Hamiltonian with instantaneous (Coulomb or Coulomb–Breit) interactions. Emergence of the two-particle Hamiltonian naturally occurs for a certain choice of the projector. At the moment, we understand in eq 89 as some “quasi potential” for a DC(B) interacting reference. The DC(B) reference, i.e., numerical solution of eq 90, already contains all reducible interaction diagrams of the instantaneous kernel, i.e., the full Coulomb(−Breit) ladder.

In atoms and molecules, the interaction of electrons and atomic nuclei (considered now as point-like, quasi-elementary particles) are dominated by electromagnetic forces. To capture most of the binding energy in these systems, it is convenient to work in the Coulomb gauge, since the instantaneous Coulomb interaction dominates the binding. Subtle magnetic effects can be accounted for by including also the instantaneous Breit interaction in the treatment. We can define an equal-time equation, a simple, linear HΨ = EΨ-type wave equation, by retaining the instantaneous part (Coulomb or Coulomb–Breit) of the interaction mediated by (subsequent exchanges of) a single photon (at a time) and the positive-energy solutions of matter. The remaining part of the exact equal-time equation can be obtained by integrating through the relative energy (relative time) of the interacting particles (in addition to a simple energy-independent correction term for double-pair instantaneous corrections in the first two terms of eq 89). The exact equal-time equation form is useful, if the correction obtained from the relative-energy integral is small. In atoms and molecules, it can be anticipated that it is small, because the electromagnetic interaction is relatively weak. During the “infinitely” long lifetime of bound systems, infinitely many photon exchanges occur, but these exchanges are mostly consecutive, there are not “many photons” present at the same time. The binding of atoms and molecules is dominated by a single photon exchange at a time, and during the lifetime of the system, there is an infinite ladder of single-photon exchanges, the Coulomb ladder or the Coulomb–Breit ladder with noncrossing steps. The effect of crossing photons can be identified under a high-energy resolution, and as a small effect it can potentially be accounted for as a (low-order) perturbative correction to the interaction ladder. If the interaction was much stronger (α was larger), there were more interaction-mediating particles present at the same time, crossed diagrams would be more important, and the equal-time separation and the no-pair approximation would be less useful.

The equal-time two-particle wave equation with instantaneous interactions could be formulated at the price of the appearance of a complicated potential energy-like term, which contains an integral with respect to the relative energy of the particles, and which can be considered as some effective potential due to the full-fledged description of the photon field for an interacting two-particle reference (no-pair DC or DCB). Sucher formulated low-order perturbative corrections to the no-pair DC(B) energy using Brillouin–Wigner perturbation theory (BWPT). The advantage of the BWPT energy formula is that it remains formally unchanged for an energy-dependent perturbation (here ),where the no-pair Hamiltonian isand Γ(E, Φi) stands for the (reduced) resolventThe Φi,n functions are eigenfunctions of the no-pair Hamiltonian with instantaneous (i) interactions. A useful relation for the quasi-Green function is obtained as follows. Using the completeness relation, we can writeThus, the quasi-Green function can be written as Aiming for a given order result, simplifications are possible. In Sucher’s α3Eh calculation, it was sufficient to consider only the first two terms in the expansion of eq 91, furthermore, the exact energy could be approximated by E ≈ Ei, i.e., and Γ(E, Φi) ≈ Γ(Ei, Φi). These approximations essentially led to first- and second-order Rayleigh–Schrödinger-type correction formulas. For the inclusion of in numerical computations, it is convenient to write it as the sum of two terms: is algebraically straightforward and corresponds to (noncrossing) pair correctionsThe technically more involved part includes an integral for the ε relative energy and carries retardation and crossed-photon contributions (e.g.,, and in Figure 1),For numerical computations, the inverse can be expanded asWe can start by considering the first term of the expansionwhere the ε integral is within the square brackets, since it is surrounded by equal-time quantities (without any relative-energy dependence, cf. eqs 88 and 89). Equation 100 can be considered as the next-order term to the known relation, eq 85: A useful identity for the inverse of the product of the one-particle propagators: Then, using eq 102, we can writewhere in the last step, we used the short notation for the one-electron propagators, eqs 75 and 76. It is convenient to consider the propagators as the sum of electronic and positronic contributions, S1 = S1+ + S1– and S2 = S2+ + S2–. Even for complicated kernels with multiple photon exchanges, remembering the sign of the imaginary component of the ε pole (positive for S1–, S2+ and negative for S1+, S2–) is useful for the identification of nonvanishing contributions. Furthermore, depending on the actual interaction, one can make arguments (following Sucher and Douglas and Kroll) about the relative importance of the contribution from the electronic and positronic subspaces. A useful relation regarding : If Ei and Φi is the eigenvalue and eigenfunction of the no-pair Hamiltonian, eq 90, thenwhere the last step can be made, since . This relation with the approximation is used during the calculations. To proceed, we can first consider a first-order perturbative correction using the no-pair eigenfunction (of eq 90):where we used the approximation , and retained only the positive-energy space contribution between and ; furthermore, we inserted the relationship of eq 105 and exploited the fact that Φi is the solution of the no-pair equation.

If we consider the solution of the no-pair DC equation, Φi = ΦC and approximate the total energy in the correction with the no-pair DC energy, E ≈ EC, then we can proceed for (Figure 1) by using the fact that the ε integral for the “homogeneous” terms (1+1+ and 2+2+) gives zero contributions (the ε-poles are either both in the positive or in the negative imaginary half plane), and only the mixed (1+2+ and 2+1+) terms have a nonvanishing contribution: Furthermore, and can be suppressed next to the ΦC no-pair wave function (and we assume z1z2 = +1 for simplicity), Then we proceed along the counterclockwise integration contour for calculating the ω integral in eq 108, and find one pole at ω = −k (+iΔ) with k = |k|, and thus obtainwhere is the noninteracting two-particle Hamiltonian. A similar calculation can be carried out for exchanged 1 and 2, and thus, the full, positive-energy transverse correction to the no-pair DC wave function iswhich reproduces Sucher’s result, eq 5.26 of ref (22). With further manipulation, Sucher obtained the Coulomb ladder correction to T++ resulting in the appearance of the interacting no-pair DC Hamiltonian in the resolventThe important part of this correction is due to the retardation of the interaction, which can be obtained by separating the instantaneous part according to (page 75 of ref (22)):where the first term gives rise to the Breit operator (cf., eqs 53 and 54)and the second term gives the perturbative retardation correction In 1958, Sucher did not have access to the numerical solution of the no-pair DC equation, so he introduced a series of approximations (including the Pauli approximation) to have final expressions for the nonrelativistic wave function. Nowadays, computer power allows us to compute and converge to “high precision” the numerical solution of the no-pair eigenvalue equation (Sec. 4), so it is a challenge to develop algorithms and computational procedures using an accurate relativistic wave function for the evaluation of perturbative corrections of . For future research, it will be a task to find practical expressions and procedures for the evaluation of the correction terms. Since the correction terms are written in an operator form (without making assumptions about using some special, e.g., one-particle, basis representation), it remains a technical and computational task to evaluate the integrals for a basis representation allowing high-precision numerical results (sections 4 and 5). For general, many-(two-) particle basis functions it may turn out to be convenient to group certain terms together (e.g., retardation and self-energy), which would otherwise be evaluated separately (i.e., in computations with one-particle basis functions).

is the simplest crossed-photon correction (Figure 1). Using Sucher’s rules (section 2.3), we can formulate the correction integral to first-order perturbation theory (and using the E ≈ Ei approximation):which can be simplified by repeated use of the residue theorem. Direct evaluation (or possible approximation) of the remaining integrals is a future task for precise no-pair wave functions computed by numerical solution of the no-pair wave eq (sections 4 and 5).

This section provides a brief overview of the practical aspects of solving the no-pair Dirac–Coulomb or Dirac–Coulomb–Breit equation with explicitly correlated trial functions. Explicitly correlated, i.e., two-particle, basis functions make it possible in practice to converge the energy to a precision where comparison of the 16-component results with precise and accurate perturbative computations (nonrelativistic QED) established in relation with precision spectroscopy is interesting and has been unexplored until recently. For the sake of this comparison, we focus on atoms and molecules of light elements, but in principle, the theoretical and algorithmic framework presented in this work is not limited to low Z systems (unlike finite-order nrQED). Starting from this section, we replace the natural units (ℏ = c = ϵ0 = 1), used in the previous section and common in molecular physics, with Hartree atomic units (ℏ = e = 1/(4πϵ0) = me = 1), convenient for quantum chemistry computations. We also note that it is not assumed that the mass of the spin-1/2 particles equals the electron mass, and so, we continue to explicitly write out the particle mass. Furthermore, the practical solution of the no-pair wave equation, eq 90, is carried out in coordinate space, instead of using the momentum-space representation, which was useful for writing down the interactions (section 2.3) and working with the propagators (section 2.4).

In the coordinate-space representation, the no-pair Dirac–Coulomb–Breit (DCB) Hamiltonian, eq 90, iswhere we wrote the projectors around the entire operator, not only around the interaction, so we can deal with only the positive-energy block, which was decoupled from the Brown–Ravenhall (+ – and − +) and negative-energy (− −) blocks already in eq 90. Hi[4](i = 1, 2) is the single-particle Dirac Hamiltonian of eq 5 shifted by the mic2 rest energyand the Ui external potential is due to the nuclei with QA = ZA electric charge, fixed at position RAand zi refers to the electric charge of ith active particle. The third term of H[16] stands for the Coulomb interaction of the particles (with r12 = |r1 – r2|)while the last term represents the instantaneous Breit interaction in coordinate representation, cf. eqs 53, 54, and 113: The symbol ⊡ stands for a block-wise direct product (also called Tracy–Singh product), which allows us to retain in the many-particle quantities the block structure of the one-particle Dirac matrix expressed with the σi (i = 1, 2, 3) Pauli matrices:Furthermore, we explicitly indicate the (k × k) dimensionality of the matrices by the [k] superscript. For the numerical implementation, we write the Hamiltonian with the σi Pauli matrices:where neglecting the B Breit term (zeroing the antidiagonal blocks) defines the no-pair Dirac–Coulomb (DC) Hamiltonian, HDC[16]. Regarding the projector, it is important to remember that the two-particle Dirac operator with instantaneous (Coulomb or Coulomb–Breit) interactions emerges from the Bethe–Salpeter equation with the operator projecting onto the positive-energy states of the noninteracting problem (section 2.4). In this context, a two-particle operator without this projector appears to be an ad hoc construct without simple connection to quantum electrodynamics.

The no-pair Hamiltonian, eq 116, is bounded from below, and thus, development of (precise) variational procedures to solve its eigenvalue equation is highly relevant for practical application of the theory. To define a good basis set, it is important to ensure a faithful matrix representation of the pi2 = pi·pi identity. During our work, fulfillment of this relation is ensured by using the so-called “restricted kinetic balance” condition, relying on the (σ[2]·a)(σ[2]·b) = (a·b)I[2] + i(a × b)σ[2] property of the Pauli matrices:The simple generalization of this one-particle balance to the two-particle case iswhere σ1[4] = σ[2] ⊗ I[2] and σ2[4] = I[2] ⊗ σ[2] and ⊗ denotes the usual Kronecker product. We have implemented the kinetic balance condition in an operator form, i.e., as a “metric”and detailed operator expressions can be found in previous work.

Finding the spectrum of H[16] (in the X-KB metric) numerically requires a finite set of basis functions, which we first define as the product of a (two-particle) spatial function and an elementary spinor (vector),Regarding the spatial part, two-particle functions can be efficiently represented in the floating explicitly correlated Gaussian (ECG) basis,where . The shift vector and the positive-definite Ai[2] matrix elements are parameters of each ECG to be determined via variational optimization (vide infra). The main advantage of working with ECGs lies in the fact that ECG matrix elements of various operators can be calculated analytically. Regarding the spinor part, 42 = 16 elementary spinors can be constructed, they are of the formwith λ = l, s (corresponding to the large and small components) and m = ± 1/2 (corresponding to the spin projection, sz),In eq 128, the “1” and “2” symbols are shown to highlight the particle index, which is defined by the position of the vector in the Kronecker product. Instead of the elementary spin representation, we can use a spinor basis which is adapted to the two-particle spin eigenstates (S = 0, MS = 0 singlet and S = 1, MS = 0, +1, −1 triplet), i.e.withandSpin-adapted functions are useful both from an interpretational and also from a practical point of view, as they make a direct connection with nonrelativistic results. This direct connection to nonrelativistic computations can be exploited for systems in which relativistic effects are small and the nonrelativistic basis parametrization provides a good starting point for relativistic computations.

For identical spin-1/2 fermions, antisymmetrized basis functions must be usedwhere the antisymmetrizer acts both on the coordinate and spinor spacewhere exchanges coordinate space labels andacts on the spinor components. In particular, = = . In two-particle computations with different spin-1/2 particles, the antisymmetrization step is, of course, omitted. Furthermore, if the system in consideration possesses additional spatial symmetries carried by elements of point group G, then, it is useful to adapt the basis functions also to these symmetries. A PG[16] operation projecting onto an irreducible representation of G can be realized by accounting for both the spatial and spin part of the problem.

We approximate eigenfunctions of H[16] in the {|ψiq(16)⟩} basis by the linear combinationwhich results in the generalized eigenvalue equation,where the Hamiltonian and the overlap matrix elements are calculated asThis is a linear variational problem for the ciq coefficients and a nonlinear variational problem for the basis function parameters, {Ai[2], si}. The coefficients are found by solving eq 137 with a given set of parameters, and the parameters can be refined by minimization of the energy for a selected eigenstate. This optimization procedure, along with (analytic) evaluation of ECG integrals has been implemented in the QUANTEN program package. For calculating the matrix elements of eq 138, the positive-energy projection of the Hamiltonian must be carried out. The matrix representation of the projector is constructed by using the positive-energy eigenstates of the noninteracting two-particle Hamiltonian, H1[4] ⊡ I[4] + I[4] ⊡ H2[4], represented as a matrix over the actual basis space. Selection of the “positive-energy” two-electron states can be realized approximately by “cutting” the noninteracting spectrum based on some energetic condition, or more precisely, by rotating the spectrum to the complex plane via complex rescaling of the electron coordinates. This complex rotation (CR) allows us to distinguish three different “branches” of the noninteracting two-electron system (positive-, Brown–Ravenhall, and negative-energy states), in principle, for any finite rotation angle. In practice, an optimal range for the angle can be found by some numerical experimentation (considering the finite precision arithmetic and the finite basis set size). For the low-Z end of the helium isoelectronic series, the cutting and the CR approaches resulted in practically identical energies, with a relative difference (much) less than 1 ppb.

Since the energetic contribution of the Breit to Coulomb interaction is small, the B[16] term of the DCB Hamiltonian can be treated as a perturbation to the DC problem, which corresponds to the H[16] = HDC[16] + partitioning. The Rayleigh–Schrödinger-type perturbative corrections to the DC energy (up to first- or second-order) are evaluated asandwhere {|Ψn⟩} and EDC,n++ are eigenfunctions and eigenvalues of the no-pair DC Hamiltonian, and BKB[16] = X12[16]†B[16]X12[16]. Brown–Ravenhall states do not require any further caution in the perturbative calculations either, since the Ψn zeroth-order states are within the positive-energy (++) space. For low nuclear charge numbers (low Z), the (second-order) perturbative and variational inclusion of the Breit interaction resulted in very small energy differences (on the order of a few ppb relative difference), which means that the one- and two-Breit photon exchange dominates the “magnetic part” of the interaction. For higher values of the nuclear charge, the difference between the two approaches is anticipated to be larger (to be explored in later work), as higher-order perturbative corrections become more important. These effects are automatically included in the variational solution, which, after all, can be thought of as the infinite-order summation of ladder diagrams. On the other hand, it is interesting to note that higher-order corrections due to the “Coulomb ladder” are significant already beyond Z = 1 (Sec. 5.1), which indicates the relevance of the development of a variational relativistic procedure.

We have originally formulated and implemented the equations for two-electron systems with fixed external charges, i.e., for Born–Oppenheimer-like relativistic computations. Most recently, it became possible to generalize these ideas to two-particle systems without external charges, i.e., pre-Born–Oppenheimer-like relativistic computations, by using a center-of-momentum frame, by considering the operators and definition of the projector according to section 2.4, which results in the emergence of a 16-component no-pair DC(B) Hamiltonian for the relative (internal) motion. The formalism, implementation details, and numerical results, tested with respect to available perturbative corrections according to section 5, are reported in ref (54).

Before our work, a few “high-precision” Dirac–Coulomb computations have been reported in the literature for helium-like ions, but the different computational procedures (with slightly different technical and theoretical details) delivered (slightly) different numerical results. Direct (and useful) comparison of these results with high-precision atomic experiments was not possible due to other, important missing (e.g., radiative and nuclear recoil) corrections. At the same time, many questions and concerns appeared in the literature regarding the Dirac–Coulomb(−Breit) “model” taken as a “starting point” and its use in a variational-type approach, the role and the correct form of the kinetic balance condition, the “choice” of a good projector for correlated computations. There had been even more controversy (and fewer solid data or formal result) regarding the inclusion of the Breit interaction in a variational treatment. Most of the observations have their own right in their own context, but the literature was very fragmented and the proper origins for a Dirac–Coulomb-based variational-type procedure with potential utility for precision spectroscopy had been obscure. At the same time, it is important to add that (various) Dirac–Coulomb(−Breit) Hamiltonian-based computational procedures have already been successfully used for compounds of heavier elements in relativistic quantum chemistry and in relation with (lower) chemical energy resolution. We have in mind (high) precision spectroscopy experiments for “calculable” systems, calculable to an in principle “arbitrary” precision, if the fundamental equations are known. So, we have had anticipated that a precise solution of some (appropriate) variant of a DC(B)-type wave equation is an important step, but it is at most halfway to the solution of the full problem, i.e., for delivering values for direct comparison with precision spectroscopy experiments. For this reason, it was of utmost importance to find good anchors for our work to established results and to the (more) complete theory, i.e., relativistic quantum electrodynamics. The primary and essential “anchor” for our work was, of course, the connection to the field theoretic Bethe–Salpeter equation that was reviewed in section 2. This formal connection clearly defines the form of the operator, the projector, and the (wave) equation which we solve, as well as, in principle, all correction terms due to retardation, pair, and radiative effects. In addition to this formal “benchmark”, it was necessary to establish numerical benchmarks to be able to check “intermediate” numerical results. Extensive testing of numerical results became possible by finding connections to (part of) the already established perturbative relativistic and QED approach based on a nonrelativistic reference. This perturbative route, sometimes called nonrelativistic QED (nrQED), is currently the state of the art for compounds of light elements, which are “calculable” systems to an almost “arbitrary precision”, and has been extensively tested in relation with precision spectroscopy experiments. The fundamental limitation of nrQED is connected with the finite-order of the available corrections in α (including also Zα), which limits the overall accuracy of the results, and this limitation provided the motivation for the present research program.

Using a computer implementation of the algorithmic details summarized in section 4, we computed the no-pair DC and DCB energies for a series of two-electron atomic and molecular systems with fixed nuclei, as well as for two-particle positronium-like systems without external charges. Are these numerical results correct? Do they (with the corresponding wave functions) represent a solid intermediate step for further potential computation of increasingly accurate relativistic QED energies for these systems? Direct comparison with experiment, due to missing corrections carried by (and the nuclear motion for the BO-type computations), is not relevant at the current stage. Numerical results of (more) “complete” nrQED computations have been extensively tested with respect to experiments, and apart from known (and conjectured) limitations of the nrQED framework, these results provide us current numerical benchmarks. At the same time, comparison of our variational no-pair Dirac–Coulomb(−Breit) energies with nrQED is not immediately obvious. In nrQED, the total (electronic) energy is written as the sum of the nonrelativistic (nr) energy and correction terms for increasing orders of the α fine-structure constant:The ε2 correction has been known as the Breit–Pauli Hamiltonian expectation value basically since Breit’s work during 1928–1931, the complete ε3 correction was first reported by Araki in 1957 and Sucher in 1958, the ε4 correction to triplet states of helium was derived by Douglas and Kroll in 1974 and also for singlet states by Yelkhovsky (and computations with Korobov) in 2001 and by Pachucki in 2006. There are currently ongoing efforts for the computation of the ε5 correction to triplet states of helium-like systems. Furthermore, for comparison with experiment, the effect of the nuclear motion is also accounted for in addition to eq 141. A recent review provides an overview of the current status for positronium-like systems. At the same time, a precise variational solution of the no-pair DCB equation provides us with the no-pair or positive-energy projected energy to all orders of α (all orders of Zα), for which the following α series can be formally written asIn eqs 141 and 142, we underlined the quantities that are primarily computed. In a variational computation, we obtain E++ to a certain numerical precision, and we want to check these computations, for testing the (correctness of the result of a complex) implementation and computational work and for gaining more insight and understanding to the numbers. We do not directly have access to the formal α expansion of the no-pair energy (right-hand side of eq 142), but by computing E++(α) for a series of slightly varied α values, we can fit an α-polynomial to the result. Coefficients of this fitted polynomial deliver us values for ε2, ε3++, ε4++, ..., resulting from (a series of) variational computations, and these values can be directly compared (tested) with respect to the relevant (part of the) nrQED corrections (right-hand side of eq 141). The second-order term in eq 142, ε2, is the same as in eq 141. Beyond second-order, the εn++ term contains only part of the εn “complete” nth-order nrQED contribution. Sucher calculated perturbative corrections to the nonrelativistic energy (in this sense, similar in spirit to nrQED), but fortunately, he reported also the no-(and single- and double-)pair part of the contributions. So, we could easily use his no-pair corrections and compare with our α3-order coefficient from variational results of helium-like ions (and two-electron molecules). Similar α3Eh-order results are available for hydrogen- and positronium-like two-particle systems from Fulton and Martin. All implementation details and extensive comparison with the perturbative results have been reported, and in the more recent papers, ref (53) regarding triplet contributions and in ref (54) about a Dirac relativistic pre-Born–Oppenheimer framework for two-particle systems without external charges. In a nutshell, an excellent agreement of the no-pair (BO and pre-BO) variational results is observed through the α scaling procedure for a series of systems, which represents an important milestone for the development of a computational relativistic QED framework for future use in relation with precision spectroscopy.

With relevance for testing and development of the fundamental theory of atomic and molecular matter, a relativistic quantum electrodynamics framework for two-spin-1/2 fermion systems (with or without external fixed nuclei) has been reviewed starting from the field theoretic Bethe–Salpeter (BS) equation. By exploiting the fact that the dominant part of the interaction (Coulomb or Coulomb–Breit) is instantaneous, it is convenient to rewrite the original BS equation to an exact equal-time form, which contains the no-pair Dirac–Coulomb(−Breit) Hamiltonian and a correction term, which carries retardation, pair, and radiative corrections. Since this correction term is anticipated to be small, a perturbative treatment has been considered. Initial ideas have been reviewed for such a perturbative treatment assuming that a highly precise approximation to the DC(B) wave function, solution of the no-pair equation including the instantaneous Coulomb(−Breit) interaction ladder, is available. For computing highly precise approximations to the no-pair DC(B) energy and wave function, implementation details have been reviewed for an explicitly correlated, variational, no-pair DC(B) computational procedure with the Born–Oppenheimer approximation as well as for extension to a pre-Born–Oppenheimer relativistic framework. The computed variational no-pair energies are tested through their α fine-structure constant dependence with respect to the relevant parts of the order-by-order computed nonrelativistic QED (nrQED) corrections. Regarding future work, it is important to realize and implement the evaluation of perturbative corrections for the retardation, pair, and radiative corrections using the variational no-pair DC(B) wave functions already computed for a series of two-particle systems. Further important work will include generalization of the theoretical framework to N-particle systems (including electron, positron, muon, and spin-1/2 nuclei), i.e., which appears to be feasible through the following steps: (a) starting from an N-particle Bethe–Salpeter wave equation; (b) identification of the relevant irreducible interaction kernels; (c) exploitation of the instantaneous character of the dominant part of the interaction; (d) emergence of the N-particle no-pair DCB wave equation for the noninteracting projectors plus a correction term including integral(s) for the relative energies; (e) solution of the no-pair DCB wave equation to high precision using explicitly correlated basis functions and a variational procedure; (f) accounting for the retardation, pair, and radiative corrections by perturbation theory; (g) testing the intermediate results with respect to the relevant terms (if known) from nrQED; (+◊) accounting for the motion of the nuclei (for systems with spin-1/2-nuclei, e.g., H2+, H2, or H3+), by using a many-particle pre-Born–Oppenheimer no-pair DCB approach through generalization of ref (54). At the moment, this endeavor appears to define an extensive research program. The present work reviewed a promising starting point based on the fundamental theory (QED) and outlined necessary practical steps. Various technical and conceptual difficulties can be foreseen. Nevertheless, the success of many-particle Dirac–Coulomb(−Breit) methodologies in relativistic quantum chemistry targeting a much lower, i.e., chemical energy resolution, as well as the limitations due to finite-order nrQED expressions suggest that the development of a computational relativistic QED framework, targeting the spectroscopic energy resolution for testing and further developing the fundamental theory of atomic and molecular matter, is relevant.