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.


INTRODUCTION: A HISTORICAL LINEUP
Dirac's one-electron space−time equation was an ingenious departure from Schrodinger'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. 1,2 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. 3 Breit attempted a two-particle generalization of Dirac's oneelectron theory in a series of papers between 1928 and 1931, 4−7 by adopting Darwin's 1920 calculation of the classical electromagnetic interaction for two moving charges 8 and the quantum mechanical velocity operator obtained with Dirac's formalism. 7 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". 5,6 This procedure was later explained based on Dirac's hole theory by Brown and Ravenhall. 9 The problem, called Brown−Ravenhall (BR) disease, related to the artificial coupling of the positive-and negative-energy states of Dirac's theory when naively 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 demon-strated that bound states of helium-like two-electron systems represented by the Breit equation have (unphysical) finite lifetimes (on the order of α 3 E h , where α is the fine-structure constant). 10,11 It has been also discussed 12 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 13 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 α 2 E h ". 11 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 14 and the development of the space−time approach to quantum electrodynamics, 15 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). 16 It is interesting to note that the same equation was written at the end of a paper without explanation by Nambu already in 1950, 17 and it was formulated by Schwinger 18 and also by Gell-Mann and Low 19 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. 20 Salpeter reported results for the hydrogenic case (of one heavy and one light particle, M ≫ m) up to the α 3 (m/M)E h 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 α 3 E h order. 21 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. 22 Sucher's final α 3 E h -order correction formulas were identical with those reported by Araki 23 a year earlier, but we can build on the fundamental ideas explained in Sucher's work for further developments.
In 1974, Douglas and Kroll 24 started their paper on the α 4 E h -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. 25 During the 1990s, Zhang worked on higher-order corrections to the fine-structure splitting (α 5 E h ) and energy levels (α 4 E h ) of helium.
Pachucki initiated a different approach starting from the late 1990s. 26−30 This approach is based on performing a Foldy− Wouthuysen transformation 31 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 32 to the required α order. In 2006, Pachucki reported the complete α 4 E h -order corrections to the energy levels of singlet helium, 33 thereby extending the 1974 work of Douglas and Kroll valid only for triplet states, as well as work from Yelkhovsky 34 and computations from Korobov and Yelkhovsky 35 in 2001 for α 4 E h -order corrections of singlet helium. In 2016, the complete α 4 E h -order corrections derived by Pachucki were used for the ground electronic state of the H 2 molecule with fixed protons. 36 Most recently, the 'Foldy−Wouthuysen−Pachucki' approach has been used to derive α 5 E h -order contributions for triplet states of helium. 37,38 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. 39 −43 For atomic and molecular  computations, the exact equal-time form of the BS equation,  as introduced by Salpeter and Sucher, 20,22 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 12,44,45 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. 48 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, 22 Douglas and Kroll, 24 as well as Salpeter. 20 The second part of the article highlights our recent work, 49−54 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 12,44,45 from the 1980s implying a possible generalization and Broyles' work from 1987 55 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, 56,57 and the corresponding two-time (equal-time) Green function idea developed by Shabaev. 32 Comparison of the Salpeter− Sucher approach with the quasi-potential method is left for future work.

Introductory Ideas and Propagators
The Dirac equation for a particle of mass m 1 and with the free-particle Hamiltonian and the α 1 and β 1 Dirac matrices. Feynman pointed out in 1949 14 that instead of working with the Hamiltonian equation, it is often more convenient to use the corresponding Green function or propagator For a Dirac particle in an external scalar field, Φ 1 z 1 stands for the charge number of the active particle and e is the elementary charge. Simple calculation 14 shows that the G 1 propagator can be obtained from the G 1 (0) free-particle propagator through the integral equation (corresponding to subsequent interaction events of the particle with the external field) as According to Feynman's combination of the electronic and positronic theory in a consistent manner, 14 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", 58 and the propagator is where ϕ̅ n (1) = ϕ n (1) *β 1 is the Dirac adjoint. The G 1 (0) (x 1 , x 1 ′) freeparticle 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 is where α = e 2 /(4π) is the fine-structure constant in natural units To describe a two-particle system, we can consider the G(x 1 , x 2 ; x 1 ′, x 2 ′) two-particle propagator or amplitude which describes that particles 1 and 2 get from x 1 ′, x 2 ′ to x 1 , x 2 space− time points. For noninteracting particles, the two-particle propagator is the simple product of the one-particle propagators, G 1 (x 1 , x 1 ′)G 2 (x 2 , x 2 ′). For interacting two-particle systems Salpeter and Bethe, 16 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 D F μν 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, corresponds 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 by 16,22,24,59 where Σ 1 is the sum of the one-electron self-energy contributions. Douglas and Kroll, 24 following the last chapter of Sucher's work, 22 formulated the two-electron equations by formally including all radiative corrections. This formulation also relies on the work of Mathews and Salam 60 who explained that the Bethe−Salpeter equation can be renormalized with the replacement of G 1 ′, G 2 ′, γ 1μ , γ 2μ , and the D F photon propagator by 22,24 We note that Mathews and Salam mostly formulated their renormalization approach based on series expansion, whereas Kallen 61 and Lehmann 62 defined renormalization terms without the use of power series expansion in the interaction constant. Karplus and Kroll 63 and Jauch and Rohrlich 64 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 22,24 (following the field theoretical derivation of the BS equation by Gell−Mann and Low 19 ), 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. which 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 as

Coordinate and Fourier Transformation: The Total and Relative Time and Energy
where the average ("total") time and relative time was introduced as It is important to note that E is the total energy of the twoparticle system and it corresponds to the T total time, eq 26. Similarly to T and t, we define Then, we obtain the following equation for the Ψ(r 1 , r 2 , t) space-and relative-time wave function + = r r r r r r r r r r with the interaction kernel depending only on the relative time variables where it is exploited that the external field is time independent, i.e., K′(x 1 , x 2 ; x 1 ′, x 2 ′) 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 22 and Douglas and Kroll 24 continued the calculation in momentum space, and we follow this line of thought. The r 1 , r 2 space coordinates of the two particles and the t relative time are replaced with the p 1 , p 2 momenta and the ε relative energy. The relative-time and relative-energy wave functions are connected by the seven-dimensional Fourier transformation while the interaction kernel in momentum space is defined as 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 I interparticle kernel, which is obtained as the sum of j I ( ) operators corresponding to irreducible "diagrams". 16,22 Action of the interaction kernel for a single-photon exchange (written in the Coulomb gauge), eqs 12−14, on some f(p 1 , p 2 , ε) two-particle function depending also on the ε relative energy can be written in the momentum-space representation as A more compact operator form of C and T is obtained by using k = (k, ω) momentum-energy translation operators. The one-particle translation operators are the two-particle translation operator is and for later convenience, we also define the notation Then, the Coulomb and transverse parts of the one-photon exchange, eqs 43 and 44, can be written as It 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: The remainder, i.e., difference of the transverse and the Breit interactions, is the retarded part, which we label as while the instantaneous contributions (Coulomb−Breit) can be handled "together": To write down the mathematical expression for more complicated j I ( ) interactions including multiple (Coulomb and/or transverse) photons (e.g., Figure 1), Sucher 22 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 p 1 = (p 1 , ε) and p 2 = (p 2 , −ε); (c) label all remaining electron lines with conserving the four-momentum.
Second, for a fully labeled diagram, j I ( ) can be formulated by writing 1. α̃1 i or α̃2 i for a transverse interaction vertex; 2. a factor S 1 (p 1 − k) for an intermediate electron line labeled with p 1 − k on the path of 1 and a factor S 2 (p 2 + k′) for an intermediate electron line labeled with p 2 + 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, is where the k-translation of the one-particle Hamiltonians gives Sucher calculated corrections to the energy up to order α 3 E h , Douglas and Kroll calculated the fine-structure splitting up to order α 4 E h , 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. 16

A Practical Wave Equation: the Exact Equal-Time Bethe−Salpeter Equation and Emergence of the No-Pair Dirac−Coulomb(−Breit) Hamiltonian
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 i instantaneous part and the "rest" The instantaneous part, i , induces only a trivial shift for the ε relative energy, and hence the effect of the relative energy can be integrated out where i is only a momentum shift integral operator In 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 as By integrating both sides with respect to the relative energy and using eqs 63 and 66, we obtain where the operator identity was used in the last step Next, we define the one-particle positive-and negative-energy projection operators for particles i = 1 and 2 by which, at this point, contains a purely formal definition for the one-particle Hamiltonian absolute value operator which also means that In short, we can also write If there is no external field, e.g., pre-BO description of two spin-1/2 particles, then, ± i reduces to the free-particle projector: 67 since for U i = 0, and the eigenvalues of the Hamiltonian absolute value operator are Using the 1 and 2 notation, we can write the 1 1 and 2 1 propagators in 1 as and similarly according to Feynman's prescription 14 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 i and Φ are independent of ε, so we only need to calculate with the two-particle projectors, To evaluate these ε integrals, we use Cauchy's residue theorem: (78) 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 identities 22 where the second identity holds in general, the first is valid only for commuting A and B operators. All in all, we obtain where the short notation is introduced, Using this result, we can rewrite eq 68 and finally obtain, the exact, equal-time Bethe−Salpeter (eBS) equation Equation 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 positiveenergy projected or no-pair Dirac−Coulomb(−Breit) equation It is important to note that in the present derivation, 20,22,24 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 12,45 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, 20 Sucher, 22 and Douglas and Kroll, 24 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, 16 i.e., the full Coulomb(−Breit) ladder.

Phenomenology: Why and When the Equal-Time Equation Is Useful?
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 singlephoton 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.

PROSPECTS REGARDING
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 22 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 energydependent perturbation (here E ( )), where the no-pair Hamiltonian is and Γ(E, Φ i ) stands for the (reduced) resolvent The Φ i,n functions are eigenfunctions of the i no-pair Hamiltonian with instantaneous (i) interactions. A useful relation for the quasi-Green function is obtained as follows. Using the = + + + Thus, the quasi-Green function can be written as Aiming for a given order result, simplifications are possible. In Sucher's α 3 E h calculation, 22 it was sufficient to consider only the first two terms in the expansion of eq 91, furthermore, the exact energy could be approximated by . These approximations essentially led to first-and second-order Rayleigh−Schrodinger-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 corrections The technically more involved part includes an integral for the ε relative energy and carries retardation and crossed-photon contributions (e.g., We can start by considering the first term of the expansion where the ε integral is within the square brackets, since it is surrounded by equal-time quantities (without any relativeenergy 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 oneparticle propagators: Then, using eq 102, we can write where in the last step, we used the short notation for the oneelectron propagators, eqs 75 and 76. It is convenient to consider the propagators as the sum of electronic and positronic contributions, S 1 = S 1+ + S 1− and S 2 = S 2+ + S 2− . Even for complicated kernels with multiple photon exchanges, remembering the sign of the imaginary component of the ε pole (positive for S 1− , S 2+ and negative for S 1+ , S 2− ) is useful for the identification of nonvanishing contributions. Furthermore, depending on the actual interaction, one can make arguments (following Sucher 22   , and retained only the positive-energy space contribution between 1 and i ; furthermore, we inserted the relationship of eq 105 and exploited the fact that Φ i is the solution of the no-pair equation.

Transverse and Retardation Correction
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 ≈ E C , then we can proceed for = T (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, + 1 and + 2 can be suppressed next to the Φ C no-pair wave function (and we assume z 1 z 2 = +1 for simplicity),

ACS Physical Chemistry Au
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 obtain 2 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 is With further manipulation, Sucher obtained the Coulomb ladder correction to T ++ resulting in the appearance of the interacting no-pair DC Hamiltonian in the resolvent The 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).

Initial Thoughts about Crossed Photon Corrections
× C C 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 ≈ E i approximation): which can be simplified by repeated use of the residue theorem. 22 Direct evaluation (or possible approximation) of the remaining integrals is a future task for precise no-pair wave ACS Physical Chemistry Au pubs.acs.org/physchemau Review functions computed by numerical solution of the no-pair wave eq (sections 4 and 5).

NUMERICAL SOLUTION OF THE NO-PAIR DIRAC−COULOMB−BREIT EIGENVALUE EQUATION
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. 49−52 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 ) = m e = 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).

No-Pair Hamiltonian
In the coordinate-space representation, the no-pair Dirac− Coulomb−Breit (DCB) Hamiltonian, eq 90, is where 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. H i [4] (i = 1, 2) is the single-particle Dirac Hamiltonian of eq 5 shifted by the m i c 2  (117) and the U i external potential is due to the nuclei with Q A = Z A electric charge, fixed at position R A and z i refers to the electric charge of ith active particle.
The third term of H [16] stands for the Coulomb interaction of the particles (with r 12 = |r 1 − r 2 |) The symbol ⊡ stands for a block-wise direct product (also called Tracy−Singh product 68−71 ), which allows us to retain in the many-particle quantities the block structure of the oneparticle Dirac matrix expressed with the σ i (i = 1, 2, 3) Pauli matrices: The simple generalization of this one-particle balance to the two-particle case is where σ 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" 49

Spatial and Spinor Basis
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, 71 with λ = l, s (corresponding to the large and small components) and m = ± 1/2 (corresponding to the spin projection, s z ), 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, M S = 0 singlet and S = 1, M S = 0, +1, −1 triplet), i.e.
Spin-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. 49−53

Symmetry Adaptation of the Basis
For identical spin-1/2 fermions, antisymmetrized basis functions must be used where 12 exchanges coordinate space labels and 12 . In two-particle computations with different spin-1/2 particles, the antisymmetrization step is, of course, omitted. 54 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 P G [16] operation projecting onto an irreducible representation of G can be realized by accounting for both the spatial and spin part of the problem. 50,53

A Variational Procedure
We approximate eigenfunctions of H [16] in the {|ψ iq (16) ⟩} basis by the linear combination This is a linear variational problem for the c iq coefficients and a nonlinear variational problem for the basis function parameters, {A i [2] , s i }. 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. 49−52,75−82 For calculating the matrix elements of eq 138, the positiveenergy projection of the Hamiltonian must be carried out. The matrix representation of the ++ [ ] 16 projector is constructed by using the positive-energy eigenstates of the noninteracting twoparticle Hamiltonian, H 1 [4] ⊡ I [4] + I [4] ⊡ H 2 [4] , represented as a matrix over the actual basis space. [49][50][51][52]54 Selection of the "positive-energy" two-electron states can be realized approximately by "cutting" the noninteracting spectrum based on some energetic condition, 49,50 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. 83 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. 49,50 4.6. Perturbative Inclusion of B [16] 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] [16] X 12 [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 (secondorder) perturbative and variational inclusion of the Breit interaction resulted in very small energy differences (on the order of a few ppb relative difference), 51,52 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 50 (Sec. 5.1), which indicates the relevance of the development of a variational relativistic procedure.

Extension to a Pre-Born−Oppenheimer Relativistic Framework
We have originally formulated and implemented the equations for two-electron systems with fixed external charges, i.e., for Born−Oppenheimer-like relativistic computations. 49−52 Most recently, it became possible to generalize these ideas to twoparticle systems without external charges, i.e., pre-Born− Oppenheimer-like 77,84−86 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.

AN OVERVIEW OF NUMERICAL RESULTS
Before our work, a few "high-precision" Dirac−Coulomb computations have been reported in the literature for heliumlike ions, 10,71,83,87 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, 10,11,83,88−90 the role and the correct form of the kinetic balance condition, 71,89,90 the "choice" of a good projector for correlated computations. 83,87,91 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 ACS Physical Chemistry Au pubs.acs.org/physchemau Review 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. 46,47 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. 78,92−97 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.

Variational vs Perturbation Theory: Perturbative Benchmark for the No-Pair Energies through α Scaling
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, 49−52 as well as for twoparticle positronium-like systems without external charges. 54 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, 4−7 the complete ε 3 correction was first reported by Araki 23 in 1957 and Sucher 22 in 1958, the ε 4 correction to triplet states of helium was derived by Douglas and Kroll 24 in 1974 and also for singlet states by Yelkhovsky 34 (and computations with Korobov 35 ) in 2001 and by Pachucki 33 in 2006. There are currently ongoing efforts 37,38 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. 98 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 as In 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. 50,51 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 22 (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 nopair corrections and compare with our α 3 -order coefficient from variational results of helium-like ions (and two-electron molecules). 50−52 Similar α 3 E h -order results are available for hydrogen-and positronium-like two-particle systems from Fulton and Martin. 21 All implementation details and extensive comparison with the perturbative results have been reported, 49−52 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 49−54 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.

ACS Physical Chemistry Au
pubs.acs.org/physchemau Review