# Quantum electrodynamics of two-body systems with arbitrary masses

Jacek Zatorski  Faculty of Physics, University of Warsaw, Pasteura 5, 02-093
Warsaw, Poland  Vojtěch Patkóš  Faculty of Mathematics and Physics, Charles
University, Ke Karlovu 3, 121 16 Prague 2, Czech Republic  Krzysztof Pachucki
Faculty of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland

###### Abstract

We perform a calculation of quantum electrodynamics effects in excited states
with $l>1$ of arbitrary two-body systems up to $\alpha^{6}\,\mu$ order. The
obtained results are valid for hadronic atoms, as long as the strong
interaction effects are negligible. We demonstrate that for circular states
with $l\sim n$, the effective expansion parameter is $Z\,\alpha/n$, which
extends the applicability of the derived formulas to heavy ions. Moreover,
inclusion of higher-order terms is feasible, which indicates that accurate
measurements of excited states of, for example, muonic or antiprotonic
hydrogenic ions can be used to determine the fundamental constants and to
search for the existence of yet unknown long-range interactions.

††preprint: Version 1.0

##  I Introduction

Finite nuclear mass corrections to atomic energy levels are often the main
limiting factor of theoretical predictions, for example in the isotope shift
of regular atoms, or transition frequencies in exotic atoms such as muonium
$e^{-}\,\mu^{+}$ [1]. It is because there is no fundamental equation for a
two-body system with arbitrary masses which accounts for relativity. This is
in contrast with the Dirac equation, which holds when one of the particles is
infinitely heavy. Although there is an exact formula [2], which gives the
leading $m/M$ correction valid for an arbitrary nuclear charge, the more
general result for an arbitrary mass ratio is not known. In the alternative
approach one employs QED and expansion in the fine structure constant $\alpha$
to derive an analytic formula for nuclear recoil corrections in consequent
orders in $\alpha$. Here we demonstrate that one can derive exact formulas for
$\alpha$ expansion coefficients for excited rotational $l>1$ states, which are
valid for arbitrary masses of its constituents. We perform calculations up to
the order $\alpha^{6}\,\mu$ and indicate the possibility of further
improvements. This would allow for extremely accurate theoretical predictions
for excited hydrogenic states as long as the contact interactions, e.g. from
strong forces, are negligible. A good example are excited states of
antiprotonic atoms, like $\bar{p}\,p$ or $\bar{p}\,\alpha$. If corresponding
measurements are available, this would allow for more accurate determination
of fundamental constants or tests of the existence of unknown long-range
interactions, similarly to precision tests performed with antiprotonic helium
[3, 4, 5], but much more accurate.

Several accurate results have been obtained for excited states of hydrogenic
systems including two-loop QED in the infinite nuclear mass limit, and leading
$m/M$ corrections. This, in combination with $\mu$H Lamb shift measurement,
allowed for the determination of the Rydberg constant and the proton charge
radius. The measurement of magnetic interactions in hydrogen-like carbon is
currently the best determination of the electron mass in atomic mass units.
The measurement of muonium $1S-2S$ will in the future provide the $\mu$ mass.
If accurate measurements are available for hydrogenic rotational states [6],
this would allow not only for more accurate determination of fundamental
constants, but also will allow for the search of long-range interactions. For
all these examples, accurate theoretical predictions of energy levels are
necessary, and in this work we demonstrate the availability of highly accurate
results for rotational levels.

The rest of this paper is organized as follows. In Section II we outline the
method of calculation. Section III contains details of the calculation of the
leading relativistic correction to energy. In Section IV we present the
leading QED correction to energy levels, and Section V is devoted to higher-
order QED corrections to the energy levels, followed by a summary in Section
VI.

##  II Expansion of energy in powers of the fine structure constant $\alpha$

In this work we consider two-body systems consisting of a negatively charged
particle of mass $m_{1}$ and a positively charged nucleus of mass $m_{2}$.
From the theoretical perspective the relevant atomic transitions are those
between states of the angular momentum number $l\approx n$. For that reason,
during calculations we will omit contributions coming from operators
containing short-range interactions, namely those of the nuclear size. We will
neglect also (electronic) vacuum polarization, which, however, can be included
separately, and its significance depends on $n$, $l$, and constituent masses.

As a consequence, the energy of a two-body system with masses $m_{1},m_{2}$,
charges $e_{1},e_{2}$, spins $s_{1},s_{2}$, and g-factors $g_{1},g_{2}$ can be
expressed as a series in powers of the fine structure constant $\alpha$,

$\displaystyle
E(\alpha)=E^{(0)}+E^{(2)}+E^{(4)}+E^{(5)}+E^{(6)}+o(\alpha^{7})\,,$ (1)

where each individual term $E^{(j)}$ is of the order $\alpha^{j}$. In
particular,

$\displaystyle E^{(0)}=m_{1}+m_{2}\,,$ (2)

and $E^{(2)}$ is the eigenvalue of the nonrelativistic two-body Hamiltonian
$H_{0}=H^{(2)}$ in the center of mass frame,

$\displaystyle
H_{0}=\frac{p^{2}}{2\,\mu}+\frac{e_{1}\,e_{2}}{4\,\pi}\,\frac{1}{r}\,.$ (3)

When we set $e_{1}=-e,e_{2}=Z\,e$, it is equal to

$E^{(2)}=E_{0}=-\frac{(Z\,\alpha)^{2}\,\mu}{2\,n^{2}}\,,$ (4)

with $\mu=m_{1}\,m_{2}/(m_{1}+m_{2})$ being the reduced mass.

##  III Leading relativistic correction $E^{(4)}$

$E^{(4)}$ is the leading relativistic correction, namely the expectation value
of the Breit Hamiltonian $H^{(4)}$ with the nonrelativistic wave function,

$\displaystyle H^{(4)}=$ $\displaystyle\ -\frac{\vec{p}^{\,4}}{8\,m_{1}^{3}}\
-\frac{\vec{p}^{\,4}}{8\,m_{2}^{3}}+\frac{e_{1}\,e_{2}}{4\,\pi}\,\biggl{\\{}\frac{1}{2\,m_{1}\,m_{2}}\,p^{i}\,\biggl{(}\frac{\delta^{ij}}{r}+\frac{r^{i}\,r^{j}}{r^{3}}\biggr{)}\,p^{j}+\frac{g_{1}\,g_{2}}{4\,m_{1}\,m_{2}}\,\biggl{[}\frac{s_{1}^{i}\,s_{2}^{j}}{r^{3}}\,\biggl{(}\delta^{ij}-3\,\frac{r^{i}\,r^{j}}{r^{2}}\biggr{)}-\frac{8\,\pi}{3}\,\vec{s}_{1}\vec{s}_{2}\,\delta^{(3)}(\vec{r})\biggr{]}$
$\displaystyle\
-\frac{\vec{r}\times\vec{p}}{2\,r^{3}}\cdot\biggl{[}\frac{g_{1}}{m_{1}\,m_{2}}\,\vec{s}_{1}+\frac{g_{2}}{m_{1}\,m_{2}}\,\vec{s}_{2}+\frac{(g_{2}-1)}{m_{2}^{2}}\,\vec{s}_{2}+\frac{(g_{1}-1)}{m_{1}^{2}}\,\vec{s}_{1}\biggr{]}\biggr{\\}}-\frac{e_{1}e_{2}}{6}\bigl{(}\left\langle
r_{E}^{2}\right\rangle_{1}+\left\langle
r_{E}^{2}\right\rangle_{2}\bigr{)}\,\delta^{(3)}(\vec{r})\,,$ (5)

where $\vec{r}=\vec{r}_{12}$, $\vec{p}=\vec{p}_{1}=-\vec{p}_{2}$, and we
neglect the quadrupole moment of any of particles. The g-factor is related to
the magnetic moment anomaly by $g=2(1+\kappa)$. For the proton it amounts to
$g_{p}=5.585\,694\,689\,3(16)$, and for the antiproton it is the same. For the
case of a spin $1/2$ point-like particle the charge radius is

$\displaystyle\left\langle r_{E}^{2}\right\rangle=$ $\displaystyle\
\frac{3}{4\,m^{2}}\,.$ (6)

For a spinless particle it is $\left\langle r_{E}^{2}\right\rangle=g=0$. Let
us note that Hamiltonian (5) does not account for any potential annihilation
effects, which would be present, for example, in positronium. It also does not
account for strong interaction effects, which are present for hadronic
particles.

From the Hamiltonian in Eq. (5) one obtains $E^{(4)}$ by taking the
expectation value with the nonrelativistic eigenfunction of $H_{0}$. The
result for a state with the principal quantum number $n$ and the angular
momentum $l$ is

$\displaystyle E^{(4)}=$
$\displaystyle\mu^{3}(Z\alpha)^{4}\bigg{\\{}\frac{1}{8\,n^{4}}\,\biggl{(}\frac{3}{\mu^{2}}-\frac{1}{m_{1}\,m_{2}}\biggr{)}-\frac{1}{\mu^{2}(2l+1)\,n^{3}}+\frac{2\,\delta_{l0}}{3\,n^{3}}\bigl{(}\left\langle
r_{E}^{2}\right\rangle_{1}+\left\langle
r_{E}^{2}\right\rangle_{2}\bigr{)}+\frac{\delta_{l0}}{m_{1}\,m_{2}\,n^{3}}$
$\displaystyle+\frac{2}{l(l+1)(2l+1)n^{3}}\bigg{[}\vec{L}\cdot\vec{s}_{1}\bigg{(}\frac{1+2\kappa_{1}}{2m_{1}^{2}}+\frac{1+\kappa_{1}}{m_{1}m_{2}}\bigg{)}+\vec{L}\cdot\vec{s}_{2}\bigg{(}\frac{1+2\kappa_{2}}{2m_{2}^{2}}+\frac{1+\kappa_{2}}{m_{1}m_{2}}\bigg{)}$
$\displaystyle-\frac{6(1+\kappa_{1})(1+\kappa_{2})}{m_{1}m_{2}\,(2\,l-1)(2\,l+3)}s_{1}^{i}s_{2}^{j}(L^{i}L^{j})^{(2)}\bigg{]}+\frac{8\,\delta_{l0}}{3m_{1}m_{2}\,n^{3}}(1+\kappa_{1})(1+\kappa_{2})\,\vec{s}_{1}\cdot\vec{s}_{2}\bigg{\\}}\,.$
(7)

In the last equation we have introduced a symmetric traceless tensor
$(L^{i}L^{j})^{(2)}$, which is defined as

$(L^{i}L^{j})^{(2)}=\frac{1}{2}\big{(}L^{i}L^{j}+L^{j}L^{i})-\frac{\delta^{ij}}{3}\vec{L}^{2}\,.$
(8)

Expression (7) is valid for arbitrary spin of both particles, point-like or
hadronic. We will now explore particular cases of the general result in Eq.
(7).

###  III.1 Dirac limit of an infinitely heavy and point-like nucleus

For an infinitely heavy $M=m_{2}$ and point-like nucleus we set $m=m_{1}$,
$\vec{s}=\vec{s}_{1}$, and $\kappa=\kappa_{1}$ to obtain

$\displaystyle
E^{(4,0)}=m\,(Z\alpha)^{4}\bigg{\\{}\frac{3}{8\,n^{4}}-\frac{1}{(2\,l+1)\,n^{3}}$
$\displaystyle+\frac{1+2\kappa}{l\,(l+1)\,(2\,l+1)\,n^{3}}\,\vec{L}\cdot\vec{s}+\frac{2\,\delta_{l0}}{3\,n^{3}}\,m^{2}\left\langle
r_{E}^{2}\right\rangle\bigg{\\}}\,.$ (9)

If the first particle is also point-like, then

$\displaystyle E^{(4,0)}=$ $\displaystyle
m\,(Z\alpha)^{4}\bigg{\\{}\frac{3}{8\,n^{4}}-\frac{1}{(2\,l+1)\,n^{3}}$
$\displaystyle+\frac{1}{l\,(l+1)\,(2\,l+1)n^{3}}\,\vec{L}\cdot\vec{s}+\frac{\delta_{l0}}{2\,n^{3}}\bigg{\\}}\,.$
(10)

Using the total angular momentum $\vec{J}=\vec{L}+\vec{s}$, it can be
rewritten in the form

$E^{(4,0)}=m\,(Z\alpha)^{4}\bigg{\\{}\frac{3}{8\,n^{4}}-\frac{1}{(2j+1)\,n^{3}}\bigg{\\}}\,.$
(11)

This can be compared to the result obtained from the Dirac equation, which is

$\displaystyle E_{D}=$ $\displaystyle\ m\,f(n,j)\,,$ (12) $\displaystyle
f(n,j)=$ $\displaystyle\
\Bigg{[}1+\frac{(Z\alpha)^{2}}{\big{(}n-\big{(}j+\frac{1}{2}\big{)}+\sqrt{(j+\frac{1}{2})^{2}-(Z\alpha)^{2}}\big{)}^{2}}\Bigg{]}^{-\frac{1}{2}}\,.$
(13)

Expanding this up to the order $(Z\alpha)^{4}$ we get

$E_{D}=m\bigg{[}1-\frac{(Z\alpha)^{2}}{2n^{2}}+(Z\alpha)^{4}\bigg{(}\frac{3}{8\,n^{4}}-\frac{1}{(2j+1)\,n^{3}}\bigg{)}\bigg{]}\,.$
(14)

The second term is nonrelativistic energy, and the third term is the
relativistic correction, in agreement with Eq. (11).

###  III.2 Atom with a spinless, point nucleus and point, spin-$1/2$ particle

For the case of a spinless point nucleus of finite mass $M=m_{2}$ and spin
$1/2$ point particle of mass $m=m_{1}$ with $\vec{s}=\vec{s}_{1}$ we get

$\displaystyle
E^{(4)}=\frac{\mu\,(Z\alpha)^{4}}{n^{3}}\bigg{\\{}\frac{1}{8\,n}\,\biggl{(}3-\frac{\mu^{2}}{m\,M}\biggr{)}-\frac{1}{(2l+1)}$
$\displaystyle+\bigg{(}\frac{\mu^{2}}{m^{2}}+\frac{2\,\mu^{2}}{mM}\bigg{)}\bigg{[}\frac{\delta_{l0}}{2}+\frac{1}{l(l+1)(2\,l+1)}\,\vec{L}\cdot\vec{s}\bigg{]}\bigg{\\}}\,.$
(15)

If the spinless nucleus is much heavier than the point spin-$1/2$ particle,
then the leading recoil correction from Eq. (15) equals

$E^{(4,1)}=\frac{m^{2}}{M}\,(Z\alpha)^{4}\bigg{\\{}-\frac{1}{2\,n^{4}}+\frac{1}{(2j+1)\,n^{3}}\bigg{\\}}\,.$
(16)

This is in agreement with the known recoil correction to the Dirac energy [2],

$\displaystyle E_{\rm rec}(n,j)=\frac{m}{2\,M}\,\bigl{(}1-f(n,j)^{2}\bigr{)}$
(17)

Expanding this to the fourth order in $Z\alpha$ we obtain a result in
agreement with Eq. (16).

###  III.3 Two point spin-1/2 particles with equal mass

For two point particles with equal mass $m_{1}=m_{2}=m$ and spin 1/2 we can
define the total spin $\vec{S}=\vec{s}_{1}+\vec{s}_{2}$ and get

$\displaystyle\
E^{(4)}=\frac{m}{8}\frac{\alpha^{4}}{n^{3}}\,\bigg{\\{}\frac{11}{8\,n}-\frac{4}{(2l+1)}+\frac{4}{3}\delta_{l0}\,\vec{S}^{2}$
$\displaystyle\
+\frac{3}{l(l+1)(2\,l+1)}\bigg{[}\vec{L}\cdot\vec{S}-\frac{2\,S^{i}S^{j}(L^{i}L^{j})^{(2)}}{(2\,l-1)(2\,l+3)}\bigg{]}\bigg{\\}}$
(18)

where

$\displaystyle\vec{L}\cdot\vec{S}=$ $\displaystyle\
\frac{1}{2}\big{[}j(j+1)-l(l+1)-s(s+1)\big{]}\,,$ (19) $\displaystyle
S^{i}S^{j}(L^{i}L^{j})^{(2)}=$ $\displaystyle\
\big{(}\vec{L}\cdot\vec{S})^{2}+\frac{1}{2}\vec{L}\cdot\vec{S}-\frac{1}{3}l(l+1)s(s+1)\,.$
(20)

This relativistic energy does not account for annihilation terms.

###  III.4 Both particles spinless and point-like

Finally, the case of both particles being spinless and point-like is

$\displaystyle
E^{(4)}=\frac{\mu(Z\alpha)^{4}}{n^{3}}\bigg{[}\frac{1}{8\,n}\biggl{(}3-\frac{\mu^{2}}{m_{1}\,m_{2}}\biggr{)}-\frac{1}{(2l+1)}+\delta_{l0}\frac{\mu^{2}}{m_{1}\,m_{2}}\bigg{]}\,.$
(21)

In the infinite nuclear mass $M=m_{2}$ limit

$\displaystyle E^{(4,0)}=$ $\displaystyle\
m(Z\alpha)^{4}\bigg{\\{}\frac{3}{8\,n^{4}}-\frac{1}{(2l+1)\,n^{3}}\bigg{\\}}\,.$
(22)

It has the same form as result (11) from the Dirac equation with $j$ replaced
by $l$, while the leading recoil correction is

$\displaystyle E^{(4,1)}=$ $\displaystyle\
\frac{m^{2}}{M}(Z\alpha)^{4}\bigg{\\{}-\frac{1}{2n^{4}}+\frac{1+\delta_{l0}}{(2l+1)\,n^{3}}\bigg{\\}}\,,$
(23)

which is different from the case of the Dirac particle in Eq. (16).

##  IV Leading QED correction $E^{(5)}$

QED effects for $l>0$ are partially accounted for in the $g$-factor, which is
present in the Breit-Pauli Hamiltonian $H^{(4)}$. Additional QED corrections
are represented as $E^{(5)}$, and for $l>0$ they have the form

$\displaystyle E^{(5)}=$ $\displaystyle\
-\frac{14\,(Z\,\alpha)^{2}}{3\,m_{1}\,m_{2}}\,\left\langle\frac{1}{4\,\pi}\,\frac{1}{r^{3}}\right\rangle-\frac{2\,\alpha}{3\,\pi}\,\biggl{(}\frac{1}{m_{1}}+\frac{Z}{m_{2}}\biggr{)}^{2}$
$\displaystyle\
\times\left\langle\,\vec{p}\,(H_{0}-E_{0})\,\ln\bigg{[}\frac{2\,(H_{0}-E_{0})}{\mu(Z\,\alpha)^{2}}\bigg{]}\,\vec{p}\right\rangle\,.$
(24)

The second term is related to the so-called Bethe logarithm $\ln[k_{0}(n,l)]$
by

$\displaystyle\ln[k_{0}(n,l)]\equiv\frac{n^{3}}{2\mu^{3}(Z\alpha)^{4}}$
$\displaystyle\times\left\langle\phi\left|\vec{p}\,(H_{0}-E_{0})\ln\biggl{[}\frac{2(H_{0}-E_{0})}{\mu(Z\alpha)^{2}}\biggr{]}\,\vec{p}\right|\phi\right\rangle$
(25)

and has been tabulated for many hydrogenic states. Let us also write $E^{(5)}$
as a function of quantum numbers $n$ and $l$. Making use of definition (IV)
and Eq. (138) we obtain

$\displaystyle E^{(5)}=$ $\displaystyle\
-\frac{7}{3\pi}\frac{(Z\alpha)^{5}\mu^{3}}{m_{1}m_{2}}\,\frac{1}{l(l+1)(2l+1)\,n^{3}}\,$
$\displaystyle\
-\frac{4}{3\pi}\biggl{(}\frac{1}{m_{1}}+\frac{Z}{m_{2}}\biggr{)}^{2}\frac{\alpha(Z\alpha)^{4}\mu^{3}}{n^{3}}\,\ln
k_{0}(n,l)$ (26)

in agreement with Eq. (19) from Ref. [7].

##  V Higher-order QED correction $E^{(6)}$

The total contribution of the order $m\,\alpha^{6}$ to energy can be split
into two parts: i) the expectation value of the first-order operators with
nonrelativistic wave function, and ii) the second-order contribution induced
by the Breit Hamiltonian $H^{(4)}$,

$E^{(6)}=\langle H^{(6)}\rangle+\langle
H^{(4)}\,\frac{1}{(E_{0}-H_{0})^{\prime}}\,H^{(4)}\rangle\,.$ (27)

Here the prime in $(E_{0}-H_{0})^{\prime}$ means that we exclude the reference
state from the resolvent. The effective operator $H^{(6)}$ can be derived
within the framework of NRQED. In general, both the first-order and the
second-order contributions separately contain singular operators that cancel
each other when combined together. However, because we are interested in
states with $l>1$ then all singularities vanish. Omitting contact
interactions, the NRQED Hamiltonian for a spin $s=0,1/2$ and finite size
particle, from Ref. [8] is

$\displaystyle\ H_{\rm
NRQED}=e\,A_{0}+\frac{1}{2\,m}\bigl{[}\vec{\pi}^{\,2}-2\,e\,(1+\kappa)\,\vec{s}\cdot\vec{B}\bigr{]}$
$\displaystyle\
-\frac{1}{8m^{3}}\bigl{[}\vec{\pi}^{\,4}-2e\,\\{p^{2}\,,\,\vec{s}\cdot\vec{B}\\}-2e\kappa\\{\vec{\pi}\cdot\vec{B},\vec{\pi}\cdot\vec{s}\\}\bigr{]}$
$\displaystyle\
+\frac{p^{6}}{16m^{5}}-\frac{1+2\kappa}{4m^{2}}\,e\,\vec{s}\cdot(\vec{E}\times\vec{\pi}-\vec{\pi}\times\vec{E})$
$\displaystyle\
-\Bigl{(}1+\frac{4\,\kappa}{3}\Bigr{)}\;\frac{3\,e}{16\,m^{4}}\,\vec{s}\,\\{p^{2},\vec{\nabla}A_{0}\times\vec{p}\\}$
$\displaystyle\
+\frac{\,e}{32\,m^{4}}\bigg{(}1+\frac{s_{a}(s_{a}+1)}{3}\bigg{)}[p^{2},[p^{2},A_{0}]]$
$\displaystyle\
-\frac{1}{2}\,\Bigl{(}\alpha_{E}-\frac{s_{a}(s_{a}+1)}{3\,m^{3}}\Bigr{)}\,e^{2}\,\vec{E}^{2}\,,$
(28)

where $[X\,,Y]\equiv X\,Y-Y\,X$ is the commutator of two operators,
$\\{X\,,\,Y\\}\equiv X\,Y+Y\,X$ is the anticommutator,
$\vec{\pi}=\vec{p}-e\,\vec{A}$, $e$ is the charge, $\kappa$ is the magnetic
moment anomaly, which is related to the $g$-factor by $g=2\,(1+\kappa)$, and
$e^{2}\,\alpha_{E}$ is the static electric dipole polarizability. For point-
like particles $\alpha_{E}=0$.

Using the NRQED Hamiltonian in Eq. (28) we can derive the effective operator
$H^{(6)}$ for the case of two spinless particles, for one spinless and one
spin-1/2 particle, and for two spin-1/2 particles. The derivation of effective
operators following Ref. [9] is presented in Appendix A. The general form of
$H^{(6)}$ is

$H^{(6)}=\sum_{i=0...9}\delta H_{i}\,.$ (29)

The individual effective operators are the following

$\displaystyle\delta H_{0}=$
$\displaystyle\frac{p^{6}}{16m_{1}^{5}}+\frac{p^{6}}{16m_{2}^{5}}\,,$ (30)
$\displaystyle\delta H_{1}=$
$\displaystyle\sum_{a}-\frac{3Z\alpha}{16m_{a}^{4}}\biggl{(}1+\frac{4\kappa_{a}}{3}\biggr{)}\vec{L}\cdot\vec{s}_{a}\,\\{p^{2},\frac{1}{r^{3}}\\}$
$\displaystyle+\frac{1}{32m_{a}^{4}}\bigg{(}1+\frac{s_{a}(s_{a}+1)}{3}\bigg{)}\,[p^{2},[p^{2},V]]\,,$
(31) $\displaystyle\delta H_{2}=$
$\displaystyle\frac{Z\alpha}{4m_{1}^{2}m_{2}^{2}}(1+2\kappa_{1})(1+2\kappa_{2})$
$\displaystyle\,(\vec{s}_{2}\times\vec{p})^{i}\left(\frac{\delta^{ij}}{r^{3}}-3\frac{r^{i}r^{j}}{r^{5}}\right)(\vec{s}_{1}\times\vec{p})^{j}\,,$
(32) $\displaystyle\delta H_{3}=$
$\displaystyle\sum_{a}\frac{1}{4m_{a}^{3}}\bigg{(}\\{p^{2},\,\vec{p}_{a}\cdot
e_{a}\vec{\cal{A}}_{a}\\}+\\{p^{2},\,\vec{s}_{a}\cdot\vec{\nabla}_{a}\times
e_{a}\vec{\cal{A}}_{a}\\}\bigg{)},$ (33) $\displaystyle\delta H_{4}=$
$\displaystyle\sum_{a}\frac{e_{a}^{2}}{2m_{a}}\,\vec{\cal{A}}_{a}^{2}\,,$ (34)
$\displaystyle\delta H_{5}=$
$\displaystyle\sum_{a}\,\frac{e_{a}^{2}\,(1+2\kappa_{a})}{2m_{a}^{2}}\,\vec{s}_{a}\cdot\vec{\cal{E}}_{a}\times\vec{\cal{A}}_{a}\,,$
(35) $\displaystyle\delta H_{6}=$ $\displaystyle\,\delta H_{6}^{A}+\delta
H_{6}^{B}+\delta H_{6}^{C}\,,$ (36)

where

$\displaystyle\delta H_{6}^{A}=$
$\displaystyle\frac{Z\alpha}{16m_{1}m_{2}}\bigg{\\{}\frac{2Z^{2}\alpha^{2}}{r^{3}}+\frac{iZ\alpha
r^{i}}{r^{3}}\bigg{[}\frac{p^{2}}{2m_{2}},\frac{r^{i}r^{j}-3\delta^{ij}\,r^{2}}{r}\bigg{]}p^{j}$
$\displaystyle-\,p^{i}\bigg{[}\frac{r^{i}r^{j}-3\delta^{ij}\,r^{2}}{r},\frac{p^{2}}{2m_{1}}\bigg{]}\frac{iZ\alpha
r^{j}}{r^{3}}-p^{i}\bigg{[}\frac{p^{2}}{2m_{2}},\bigg{[}\frac{r^{i}r^{j}-3\delta^{ij}\,r^{2}}{r},\frac{p^{2}}{2m_{1}}\bigg{]}\bigg{]}p^{j}\bigg{\\}}+(1\leftrightarrow
2)\,,$ (37) $\displaystyle\delta H_{6}^{B}=$
$\displaystyle\frac{Z\alpha}{4m_{1}m_{2}}\bigg{\\{}(1+\kappa_{1})\bigg{[}\bigg{(}\vec{s}_{1}\times\frac{\vec{r}}{r}\bigg{)}^{i},\,\frac{p^{2}}{2m_{1}}\bigg{]}\frac{iZ\alpha
r^{i}}{r^{3}}-(1+\kappa_{2})\frac{iZ\alpha
r^{i}}{r^{3}}\bigg{[}\frac{p^{2}}{2m_{2}},\bigg{(}\vec{s}_{2}\times\frac{\vec{r}}{r}\bigg{)}^{i}\bigg{]}$
$\displaystyle+\,(1+\kappa_{1})\bigg{[}\frac{p^{2}}{2m_{2}},\bigg{[}\bigg{(}\vec{s}_{1}\times\frac{\vec{r}}{r}\bigg{)}^{i},\,\frac{p^{2}}{2m_{1}}\bigg{]}\bigg{]}p^{i}+(1+\kappa_{2})p^{i}\bigg{[}\frac{p^{2}}{2m_{2}},\bigg{[}\bigg{(}\vec{s}_{2}\times\frac{\vec{r}}{r}\bigg{)}^{i},\,\frac{p^{2}}{2m_{1}}\bigg{]}\bigg{]}\bigg{\\}}+(1\leftrightarrow
2)\,,$ (38) $\displaystyle\delta H_{6}^{C}=$
$\displaystyle-\frac{Z\alpha(1+\kappa_{1})(1+\kappa_{2})}{4m_{1}^{2}m_{2}^{2}}\bigg{[}p^{2},\bigg{[}p^{2},\,\vec{s}_{1}\vec{s}_{2}\frac{2}{3r}+s_{1}^{i}s_{2}^{j}\frac{1}{2r}\bigg{(}\frac{r^{i}r^{j}}{r^{2}}-\frac{\delta^{ij}}{3}\bigg{)}\bigg{]}\bigg{]}\,,$
(39) $\displaystyle\delta
H_{7}=\sum_{a}\frac{e_{a}^{2}\,(1+2\kappa_{a})}{2m_{a}^{2}}\;\vec{s}_{a}\cdot\vec{\cal{E}}_{a}\times\vec{\cal{A}}_{a}$
$\displaystyle+\frac{ie_{a}\,(1+2\kappa_{a})}{8m_{a}^{3}}\,\big{[}\,\vec{\cal{A}}_{a}\cdot(\vec{p}_{a}\times\vec{s}_{a})+(\vec{p}_{a}\times\vec{s}_{a})\cdot\vec{\cal{A}}_{a},p_{a}^{2}\,\big{]}\,,$
$\displaystyle\delta
H_{8}=-\frac{1}{2}\sum_{a}\bigg{(}\alpha_{Ea}-\frac{s_{a}(s_{a}+1)}{3\,m_{a}^{3}}\bigg{)}\,\frac{Z^{2}\alpha^{2}}{r^{4}}\,,$
(41) $\displaystyle\delta
H_{9}=\sum_{a}\frac{e_{a}\kappa_{a}}{4m_{a}^{3}}\\{\vec{p}_{a}\cdot\vec{\nabla}_{a}\times\vec{\cal{A}}_{a},\vec{p}_{a}\cdot\vec{s}_{a}\\}\,.$
(42)

The static fields ${\cal A}_{i}$ and ${\cal E}_{i}$ are defined as

$\displaystyle e_{1}{\cal{A}}^{i}_{1}=$
$\displaystyle-\frac{Z\alpha}{2r}\left(\delta^{ij}+\frac{r^{i}r^{j}}{r^{2}}\right)\frac{p_{2}^{j}}{m_{2}}-\frac{Z\alpha}{m_{2}}\frac{\left(\vec{s}_{2}\times\vec{r}\right)^{i}}{r^{3}}(1+\kappa_{2})\;,$
(43) $\displaystyle e_{2}{\cal{A}}^{i}_{2}=$
$\displaystyle-\frac{Z\alpha}{2r}\left(\delta^{ij}+\frac{r^{i}r^{j}}{r^{2}}\right)\frac{p_{1}^{j}}{m_{1}}+\frac{Z\alpha}{m_{1}}\frac{\left(\vec{s}_{1}\times\vec{r}\right)^{i}}{r^{3}}(1+\kappa_{1})\;,$
(44) $\displaystyle e_{1}\vec{\cal{E}}_{1}=$
$\displaystyle-Z\alpha\frac{\vec{r}}{r^{3}},\,\,e_{2}\vec{\cal{E}}_{2}=Z\alpha\frac{\vec{r}}{r^{3}}\,.$
(45)

In the case when one of the particles is spinless and the other one has spin
$1/2$, the effective operators (30-42) are changed in such a way that
$\vec{s}_{1}=s_{1}=0$ in $\delta H_{i}$ and also in $e_{2}\vec{\cal{A}}_{2}$.
While, in the case of two spinless particles, the effective operators (30-42)
are modified by setting $\vec{s}_{1}=\vec{s}_{2}=s_{1}=s_{2}=0$ and the same
for static fields $e_{a}\vec{\cal{A}}_{a}$.

The evaluation of the expectation value of $H^{(6)}$ with nonrelativistic wave
function and the second-order contribution with the Breit Hamiltonian for
individual cases of particles is then presented in Appendices B, C, and D. The
results for energies are presented in the following subsections.

###  V.1 Two particles with spins $s_{1}=s_{2}=0$

We will begin with the case when both the positively charged nucleus and
negatively charged particle are spinless. The final result for the energy
$E^{(6)}$ is the sum of the first-order and second-order contributions
evaluated in Appendix B. Defining

$\displaystyle f^{(6)}(n,j)=$
$\displaystyle-\frac{5}{16n^{6}}+\frac{3}{2(2j+1)n^{5}}-\frac{3}{2(2j+1)^{2}n^{4}}$
$\displaystyle-\frac{1}{(2j+1)^{3}n^{3}}$ (46)

we obtain the correction of the order $m\,\alpha^{6}$

$\displaystyle E^{(6)}=$
$\displaystyle\mu(Z\alpha)^{6}\biggl{[}f^{(6)}(n,l)+\frac{\mu^{2}}{m_{1}\,m_{2}}\bigg{(}\frac{3}{16\,n^{6}}-\frac{8\,l(l+1)-3}{2\,(2l-1)(2l+1)(2l+3)\,n^{5}}+\frac{6}{(2l-1)(2l+1)(2l+3)\,n^{3}}\bigg{)}$
$\displaystyle-\frac{\mu^{4}}{16\,m_{1}^{2}\,m_{2}^{2}\,n^{6}}+\frac{2\,\mu^{3}\,(\alpha_{E1}+\alpha_{E2})}{(2l-1)(2l+1)(2l+3)}\bigg{(}\frac{1}{n^{5}}-\frac{3}{l(l+1)\,n^{3}}\bigg{)}\biggr{]}\,.$
(47)

Similarly to $E^{(4)}$, we can explore particular limiting cases of the
general result in Eq. (47). Firstly we perform the limit where mass $M=m_{2}$
of the nucleus is infinitely heavy, to obtain the result in agreement with the
Klein-Gordon equation:

$E^{(6,0)}=m(Z\alpha)^{6}\,f^{(6)}(n,l)\,.$ (48)

We may also expand Eq. (47) up to the order $m^{2}/M$ to obtain the first-
order recoil correction

$\displaystyle E^{(6,1)}=$ $\displaystyle\
(Z\alpha)^{6}\,\frac{m^{2}}{M}\,f^{(6,1)}_{0}(n,l)\,,$ (49)

where we introduce the function

$\displaystyle
f^{(6,1)}_{0}(n,l)=\frac{1}{2n^{6}}+\frac{6-10l(l+1)}{(2l-1)(2l+1)(2l+3)n^{5}}$
$\displaystyle+\frac{3}{2(2l+1)^{2}n^{4}}+\frac{3+28l(l+1)}{(2l-1)(2l+1)^{3}(2l+3)n^{3}}\,.$
(50)

This concludes the evaluation of the correction to energy for a purely
spinless two-body system.

###  V.2 Two particles with spins $s_{1}=0,\,s_{2}=1/2$

For the case when one of the particles is spinless while the other one has
spin 1/2, the derivation of the first-order and second-order contributions is
presented in Appendix C. In the resulting expression for energy $E^{(6)}$ we
can separate the part which is already contained in the final result for two
spinless particles, Eq. (47). If we denote the spinless energy (47) as
$E^{(6)}_{s_{1}=s_{2}=0}$, then

$\displaystyle E^{(6)}=E^{(6)}_{s_{1}=s_{2}=0}$
$\displaystyle+\frac{\mu(Z\alpha)^{6}}{l(l+1)(2l-1)(2l+1)(2l+3)}\big{[}A+B\,\vec{L}\cdot\vec{s}_{2}\big{]}\,.$
(51)

With the help of substitutions

$\displaystyle\lambda_{0}=$ $\displaystyle\
\frac{(2\,l-1)(2\,l+3)}{2\,l+1}=2\,l+1-\frac{4}{2\,l+1}\,,$ (52)
$\displaystyle\lambda_{1}=$ $\displaystyle\
\frac{3}{2l(l+1)}+\frac{4}{(2l+1)^{2}}\,,$ (53) $\displaystyle\lambda_{2}=$
$\displaystyle\
-\frac{3}{4l^{2}}+\frac{13}{4l(l+1)}-\frac{3}{4(l+1)^{2}}+\frac{8}{(2l+1)^{2}}\,,$
(54)

and the definition

$X=\frac{X_{5}}{n^{5}}+\frac{X_{4}}{n^{4}}+\frac{X_{3}}{n^{3}}$ (55)

for $X=A,B$ the coefficients $A$ and $B$ in Eq. (V.2) can be written as

$\displaystyle A_{5}=$ $\displaystyle\
\frac{l(l+1)}{2}\bigg{[}g_{2}^{2}\frac{\mu^{2}}{4m_{2}^{2}}+(-2-3g_{2}+g_{2}^{2})\frac{\mu^{3}}{2m_{2}^{3}}+\frac{\mu^{4}}{m_{2}^{4}}\bigg{]}\,,$
$\displaystyle A_{4}=$ $\displaystyle\
\frac{3\lambda_{0}}{4}\bigg{[}-g_{2}^{2}\frac{\mu^{2}}{2m_{2}^{2}}+g_{2}\frac{\mu^{3}}{m_{2}^{3}}-\frac{\mu^{4}}{2m_{2}^{4}}\bigg{]}\,,$
$\displaystyle A_{3}=$ $\displaystyle\
g_{2}^{2}\bigg{(}\lambda_{1}-\frac{9}{2}\bigg{)}\frac{\mu^{2}}{4m_{2}^{2}}+\big{(}3(2+5g_{2}-g_{2}^{2})-2g_{2}\lambda_{1}\big{)}\frac{\mu^{3}}{4m_{2}^{3}}$
$\displaystyle\ +(\lambda_{1}-9)\frac{\mu^{4}}{4m_{2}^{4}}\,$ (56)

and

$\displaystyle B_{5}=$ $\displaystyle\
(2l+1)\lambda_{0}\bigg{[}-2g_{2}\frac{\mu}{m_{2}}-3(g_{2}+1)\frac{\mu^{3}}{2m_{2}^{2}}+\frac{3\mu^{4}}{2m_{2}^{4}}\bigg{]}+\big{(}8(5+3g_{2})l(l+1)-3(10+4g_{2}+g_{2}^{2})\big{)}\frac{\mu^{2}}{4m_{2}^{2}}\,,$
$\displaystyle B_{4}=$ $\displaystyle\
\frac{3\lambda_{0}}{l(l+1)}\bigg{[}g_{2}l(l+1)\frac{\mu}{m_{2}}+\frac{g_{2}^{2}-4l(l+1)}{4}\frac{\mu^{2}}{m_{2}^{2}}-\frac{g_{2}}{2}\frac{\mu^{3}}{m_{2}^{3}}+\frac{\mu^{4}}{4m_{2}^{4}}\bigg{]}\,,$
$\displaystyle B_{3}=$ $\displaystyle\
2g_{2}(3-\lambda_{1})\frac{\mu}{m_{2}}+\big{(}g_{2}^{2}\lambda_{2}-6(1+g_{2})+2\lambda_{1}\big{)}\frac{\mu^{2}}{m_{2}^{2}}+(3-\lambda_{2})\bigg{[}2g_{2}\frac{\mu^{3}}{m_{2}^{3}}-\frac{\mu^{4}}{m_{2}^{4}}\bigg{]}\,.$
(57)

In order to verify the correctness of our results we will compare them with
known results for limiting cases. We restrict now to point-like particles by
setting $\alpha_{Ea}=0$. Firstly, we perform the limit of
$m_{1}\rightarrow\infty$ in Eq. (V.2), i.e., the limit of an infinitely heavy
spinless particle. Defining

$k=(l-j)(2j+1)$ (58)

then Eq. (V.2) leads to the nonrecoil result

$\displaystyle\lim_{m_{1}\rightarrow\infty}E^{(6)}=m_{2}(Z\alpha)^{6}\bigg{[}f^{(6)}(n,j)$
$\displaystyle+\kappa_{2}\bigg{(}\frac{-9+19k+16k^{2}}{2|k|(2k-1)(2k+1)(2k+3)\,n^{5}}-\frac{3}{2k^{2}(2k+1)\,n^{4}}-\frac{-3-5k+37k^{2}+66k^{3}+24k^{4}}{2|k|k^{2}(k+1)(2k-1)(2k+1)^{2}(2k+3)\,n^{3}}\bigg{)}$
$\displaystyle+\kappa_{2}^{2}\bigg{(}\frac{3(k+1)}{2|k|(2k-1)(2k+1)(2k+3)\,n^{5}}-\frac{3}{2k^{2}(2k+1)^{2}\,n^{4}}-\frac{-3-5k+55k^{2}+120k^{3}+60k^{4}}{2|k|k^{2}(k+1)(2k-1)(2k+1)^{3}(2k+3)\,n^{3}}\bigg{)}\bigg{]}\,.$
(59)

For $\kappa_{2}=0$ and both $l=j\pm 1/2$ the expression (V.2) coincides with
the Dirac result acquired by expanding the expression on the right-hand side
of Eq. (12) up to the order $(Z\alpha)^{6}$. The leading recoil correction of
the order $m_{2}^{2}/m_{1}$ can be compared with the result from the
literature. Defining the function

$\displaystyle
f^{(6,1)}_{1/2}(n,k)=\frac{1}{2n^{6}}-\frac{-3-2k+14k^{2}+10k^{3}}{|k|(2k-1)(2k+1)(2k+3)n^{5}}+\frac{3}{8k^{2}n^{4}}+\frac{-3-2k+12k^{2}+56k^{3}}{8|k|k^{2}(2k-1)(2k+1)(2k+3)n^{3}}$
(60)

then by expanding Eq. (V.2) for large $m_{1}$ up to the order
$m_{2}^{2}/m_{1}$, we get the leading recoil correction

$\displaystyle
E^{(6,1)}=(Z\alpha)^{6}\,\frac{m_{2}^{2}}{m_{1}}\,\biggl{[}f^{(6,1)}_{1/2}(n,k)$
$\displaystyle+\kappa_{2}\bigg{(}-\frac{-12+27k+22k^{2}}{|k|(2k-1)(2k+1)(2k+3)n^{5}}+\frac{3}{k^{2}(2k+1)n^{4}}+\frac{-3-5k+49k^{2}+96k^{3}+36k^{4}}{|k|k^{2}(k+1)(2k-1)(2k+1)^{2}(2k+3)n^{3}}\bigg{)}$
$\displaystyle+\kappa_{2}^{2}\bigg{(}-\frac{9+11k}{2|k|(2k-1)(2k+1)(2k+3)n^{5}}+\frac{9}{2k^{2}(2k+1)^{2}n^{4}}+\frac{3\big{(}-3-5k+57k^{2}+128k^{3}+68k^{4}\big{)}}{2|k|k^{2}(k+1)(2k-1)(2k+1)^{3}(2k+3)n^{3}}\bigg{)}\bigg{]}\,.$
(61)

This is in agreement with the known result from [10] for $l=1$ and
$\kappa_{2}=0$. For the more general case of $l>0$ the energy correction can
be obtained by means of combining the expansion of Eq. (17) up to the order
$(Z\alpha)^{6}$ with Eq. (75) from Ref. [11]. The resulting sum can be written
as follows

$\displaystyle
E^{(6,1)}(n,j,l)=\frac{m_{2}^{2}}{m_{1}}(Z\alpha)^{6}\biggl{[}\frac{1}{2n^{6}}-\frac{2}{(2j+1)n^{5}}+\frac{3}{2(2j+1)^{2}n^{4}}+\frac{1}{(2j+1)^{3}n^{3}}\biggr{]}+\frac{m_{2}^{2}}{m_{1}}\frac{(Z\alpha)^{2}}{2\mu^{4}}\left\langle\frac{\vec{L}^{2}}{r^{4}}\right\rangle\,,$
(62)

where $\langle 1/r^{4}\rangle$ is given in Eq. (139). Comparing equation (61)
for $\kappa_{2}=0$ with (62) for $l=j\pm 1/2$ confirms that they agree.

The next term in the mass ratio expansion, of the order $m_{2}^{3}/m_{1}^{2}$
for $\kappa_{2}=0$, is

$\displaystyle E^{(6,2)}=$
$\displaystyle(Z\alpha)^{6}\,\frac{m_{2}^{3}}{m_{1}^{2}}\,\biggl{[}-\frac{15}{16n^{6}}+\frac{-9-22k+84k^{2}+72k^{3}}{4|k|(2k-1)(2k+1)(2k+3)n^{5}}-\frac{3(2k-1)}{8k^{2}(2k+1)n^{4}}$
$\displaystyle-\frac{3-k-20k^{2}+184k^{3}+480k^{4}+304k^{5}}{8|k|k^{2}(k+1)(2k-1)(2k+1)^{2}(2k+3)n^{3}}\bigg{]}\,,$
(63)

which can be considered as a new result.

Secondly, we can verify that in the limiting case of Eq. (V.2), when the
spin-$1/2$ particle becomes infinitely heavy, we obtain the result in
agreement with that obtained from the Klein-Gordon equation, i.e.,

$\lim_{m_{2}\rightarrow\infty}E^{(6)}=m_{1}(Z\alpha)^{6}f^{(6)}(n,l)\,.$ (64)

Similarly as in the leading relativistic correction $E^{(4)}$ this coincides
with the Dirac result in Eq. (V.2) by replacing $j$ by $l$ and setting
$\kappa_{2}=0$. Lastly, we also present the next term in the expansion of Eq.
(V.2) in terms of the $m_{1}^{2}/m_{2}$ ratio. Namely,

$\displaystyle E^{(6,1)}=$
$\displaystyle(Z\alpha)^{6}\,\frac{m_{1}^{2}}{m_{2}}\,\biggl{[}f^{(6,1)}_{0}(n,l)+g_{2}\,\vec{L}\cdot\vec{s}_{2}\biggl{(}-\frac{2}{l(l+1)(2l+1)n^{5}}+\frac{3}{l(l+1)(2l+1)^{2}n^{4}}+\frac{1+6l(l+1)}{l^{2}(l+1)^{2}(2l+1)^{3}n^{3}}\biggr{)}\biggr{]}\,.$

where $f^{(6,1)}_{0}(n,l)$ is defined in Eq. (V.1).

###  V.3 Two particles with spins $s_{1}=s_{2}=1/2$

Finally, we present the result for the energy $E^{(6)}$ for the case of two
spin-1/2 particles. The corresponding derivation of the first-order and
second-order contributions is presented in Appendix D. We express the result
using a set of basic spin-dependent operators. To do that we use the
identities presented in Appendix G. Moreover, to keep the formulas compact, we
introduce the following symbols

$\displaystyle\lambda_{3}=-\frac{3}{(2l-1)(2l+3)}+\frac{1}{2(2l+1)^{2}}\,,$
(66)
$\displaystyle\lambda_{4}=-\frac{3}{4l(l+1)}-\frac{3}{(2l-1)(2l+3)}-\frac{3}{2(2l+1)^{2}}\,,$
(67)
$\displaystyle\lambda_{5}=\frac{3}{2l^{2}}-\frac{1}{2l(l+1)}+\frac{3}{2(l+1)^{2}}\,,$
(68)
$\displaystyle\lambda_{6}=\frac{3}{4l^{2}}-\frac{9}{4l(l+1)}+\frac{3}{4(l+1)^{2}}+\frac{4}{(2l-1)(2l+3)}$
$\displaystyle-\frac{6}{(2l+1)^{2}}\,,$ (69)
$\displaystyle\lambda_{7}=-\frac{4}{(2l-1)(2l+3)}+\frac{2}{3(2l+1)^{2}}\,,$
(70)
$\displaystyle\lambda_{8}=-\frac{3}{l^{2}}+\frac{7}{l(l+1)}-\frac{3}{(l+1)^{2}}-\frac{72}{(2l-1)(2l+3)}$
$\displaystyle+\frac{20}{(2l+1)^{2}}\,,$ (71)
$\displaystyle\lambda_{9}=-\frac{1}{l^{2}}+\frac{3}{l(l+1)}-\frac{1}{(l+1)^{2}}-\frac{12}{(2l-1)(2l+3)}$
$\displaystyle+\frac{8}{(2l+1)^{2}}\,,$ (72)
$\displaystyle\lambda_{10}=-\frac{3}{l^{2}}+\frac{13}{l(l+1)}-\frac{3}{(l+1)^{2}}+\frac{36}{(2l-1)(2l+3)}$
$\displaystyle+\frac{32}{(2l+1)^{2}}\,,$ (73)
$\displaystyle\lambda_{11}=\frac{15}{l^{2}}-\frac{29}{l(l+1)}+\frac{15}{(l+1)^{2}}-\frac{12}{(2l-1)^{2}}$
$\displaystyle+\frac{116}{(2l-1)(2l+3)}-\frac{96}{(2l+1)^{2}}-\frac{12}{(2l+3)^{2}}\,,$
(74)

in addition to $\lambda_{0}$, $\lambda_{1}$, and $\lambda_{2}$ defined in the
previous section. The final result for $E^{(6)}$ of two spin-1/2 particles is
then

$E^{(6)}=E^{(6)}_{s_{1}=s_{2}=0}+\frac{\mu(Z\alpha)^{6}}{l(l+1)(2l-1)(2l+1)(2l+3)}\bigg{[}A+B\,\vec{L}\cdot\vec{s}_{1}+C\,\vec{L}\cdot\vec{s}_{2}+D\,\vec{s}_{1}\cdot\vec{s}_{2}+F\,(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}\bigg{]}\,,$
(75)

where we have again subtracted the spinless result $E^{(6)}_{s_{1}=s_{2}=0}$
given by Eq. (47). Using once more the definition (55) for $X=A,B,C,D,F$ we
get the following coefficients

$\displaystyle A_{5}=$ $\displaystyle\
\frac{l(l+1)\,(g_{1}-2)}{4}\bigg{(}\frac{g_{1}+2}{2}\frac{\mu^{2}}{m_{1}^{2}}+(g_{1}-1)\frac{\mu^{3}}{m_{1}^{3}}\bigg{)}+\frac{\mu^{4}}{m_{1}^{2}m_{2}^{2}}\frac{l(l+1)}{2}+(1\leftrightarrow
2)\,,$ $\displaystyle A_{4}=$ $\displaystyle\
-\frac{3\lambda_{0}}{16}-\frac{3\lambda_{0}}{4}(g_{1}-2)\bigg{(}\frac{g_{1}+2}{2}\frac{\mu^{2}}{m_{1}^{2}}-\frac{\mu^{3}}{m_{1}^{3}}\bigg{)}-\frac{\mu^{4}}{m_{1}^{2}m_{2}^{2}}\frac{3}{8(2l+1)}\bigg{(}4l(l+1)+\frac{3g_{1}^{2}g_{2}^{2}}{16}\,-3\bigg{)}+(1\leftrightarrow
2)\,,$ $\displaystyle A_{3}=$ $\displaystyle\
\frac{\lambda_{1}-3}{8}+\frac{(g_{1}-2)}{4}\bigg{(}\frac{(g_{1}+2)}{2}(-9+2\lambda_{1})\frac{\mu^{2}}{m_{1}^{2}}-(-9+3g_{1}+2\lambda_{1})\frac{\mu^{3}}{m_{1}^{3}}\bigg{)}$
$\displaystyle\
+\frac{\mu^{4}}{m_{1}^{2}m_{2}^{2}}\frac{1}{2}\bigg{(}-\frac{9}{2}+\lambda_{3}+\frac{g_{1}^{2}g_{2}^{2}\,\lambda_{4}}{16}-\lambda_{4}\bigg{)}+(1\leftrightarrow
2)\,,$ (76)

where the symbol $(1\leftrightarrow 2)$ stands for the exchange of masses
$m_{1}\leftrightarrow m_{2}$ and $g$-factors $g_{1}\leftrightarrow g_{2}$ of
both particles. For the spin-orbit contribution coefficient $B$ we get the
following result,

$\displaystyle B_{5}=$ $\displaystyle\
\frac{\mu^{2}}{8m^{2}_{1}}\big{(}24(g_{1}-1)-6g_{1}^{2}-3g_{1}g_{2}(1-g_{2})+16l(l+1)(2-g_{1})\big{)}+\frac{\mu^{2}}{8m_{1}m_{2}}g_{1}\big{(}48+3g_{2}(1-g_{2})-64l(l+1)\big{)}$
$\displaystyle\
-\frac{3\mu^{3}}{8m_{1}^{3}}\big{(}12(1-g_{1})-g_{1}g_{2}(1-g_{2})-16l(l+1)(1-g_{1})\big{)}+\frac{\mu^{4}}{m_{1}^{2}m_{2}^{2}}\bigg{(}-\frac{9}{2}-\frac{3}{8}g_{1}g_{2}+\frac{3}{32}g_{1}^{2}g_{2}^{2}+6l(l+1)\bigg{)}\,,$
$\displaystyle B_{4}=$ $\displaystyle\
\frac{3}{8l(l+1)}\bigg{[}\frac{\mu^{2}}{m^{2}_{1}}\lambda_{0}\big{(}2(g_{1}^{2}-1)+g_{1}g_{2}(1-g_{2})+8l(l+1)(g_{1}-1)\big{)}+\frac{\mu^{2}}{m_{1}m_{2}}g_{1}\lambda_{0}\big{(}g_{2}(g_{2}-1)+8l(l+1)\big{)}$
$\displaystyle+\frac{\mu^{3}}{m_{1}^{3}}\lambda_{0}(4(1-g_{1})-g_{1}g_{2}(1-g_{2}))+\frac{\mu^{4}}{m_{1}^{2}m_{2}^{2}}\frac{9g_{1}^{2}g_{2}^{2}+(g_{1}g_{2}+2)\big{(}32l(l+1)-24\big{)}}{8(2l+1)}\bigg{]}\,,$
$\displaystyle B_{3}=$ $\displaystyle\
\frac{\mu^{2}}{2m^{2}_{1}}\big{(}-6+\big{(}2g_{1}(g_{1}-1)+g_{1}g_{2}(1-g_{2})\big{)}\lambda_{2}+(1-g_{1})\lambda_{5}\big{)}+\frac{\mu^{2}}{2m_{1}m_{2}}g_{1}\big{(}12-(2+g_{2}(1-g_{2}))\lambda_{2}-\lambda_{5}\big{)}$
$\displaystyle\
+\frac{\mu^{3}}{2m_{1}^{3}}\big{(}12(g_{1}-1)+\big{(}4(1-g_{1})-g_{1}g_{2}(1-g_{2})\big{)}\lambda_{2}\bigg{)}+\frac{\mu^{4}}{2m_{1}^{2}m_{2}^{2}}\bigg{(}-6+\big{(}2+g_{1}g_{2}\big{)}\lambda_{2}+\frac{3}{8}g_{1}^{2}g_{2}^{2}\lambda_{6}\bigg{)}\,.$
(77)

The coefficient $C$ is obtained from $B$ simply by the exchange
$(1\leftrightarrow 2)$. The spin-spin coefficient $D$ reads

$\displaystyle D_{5}=$ $\displaystyle\
\frac{\mu^{2}}{m_{1}m_{2}}\frac{l(l+1)\,g_{1}g_{2}}{6}+\frac{\mu^{4}}{m_{1}^{2}m_{2}^{2}}\frac{4l(l+1)}{3}\,,$
$\displaystyle D_{4}=$ $\displaystyle\
-\frac{\mu^{3}}{m_{1}^{2}m_{2}}(g_{1}-2)\lambda_{0}-\frac{\mu^{2}}{2m_{1}m_{2}}\lambda_{0}\big{(}2(1-g_{1})+g_{1}g_{2}\big{)}-\frac{\mu^{4}}{m_{1}^{2}m_{2}^{2}}\frac{3\,\big{(}g_{1}^{2}g_{2}^{2}-16\big{)}+64l(l+1)}{32(2l+1)}+(1\leftrightarrow
2)\,,$ $\displaystyle D_{3}=$ $\displaystyle\
\frac{\mu^{3}}{m_{1}^{2}m_{2}}\frac{2(g_{1}-2)(\lambda_{1}-3)}{3}+\frac{\mu^{2}}{12m_{1}m_{2}}\big{(}24(g_{1}-1)-15g_{1}g_{2}+4\big{(}2(1-g_{1})+g_{1}g_{2}\big{)}\lambda_{1}\big{)}$
$\displaystyle+\frac{\mu^{4}}{2m_{1}^{2}m_{2}^{2}}\bigg{(}\frac{1}{12}\big{(}g_{1}^{2}g_{2}^{2}-16\big{)}\lambda_{4}-6+\lambda_{7}\bigg{)}+(1\leftrightarrow
2)\,.$ (78)

Finally, the coefficient $F$ of the tensor spin-spin interaction is the most
complicated and equals

$\displaystyle F_{5}=$ $\displaystyle\
\frac{1}{(2l-1)(2l+3)}\bigg{[}\frac{\mu^{3}}{m^{2}_{1}m_{2}}\frac{9}{4}(g_{1}-2)\big{(}4-6g_{2}+3g_{1}g_{2}+16\,l(l+1)\big{)}$
$\displaystyle+\frac{\mu^{2}}{2m_{1}m_{2}}\bigg{(}18(1-g_{1})-\frac{45}{4}g_{1}g_{2}+l(l+1)\big{(}47g_{1}g_{2}+72(1-g_{1})\big{)}\bigg{)}$
$\displaystyle+\frac{\mu^{4}}{m_{1}^{2}m_{2}^{2}}\bigg{(}\lambda_{0}(5-3g_{1}g_{2})(2l+1)+6+\frac{9}{2}g_{1}g_{2}-\frac{9}{8}g_{1}^{2}g_{2}^{2}\bigg{)}\bigg{]}+(1\leftrightarrow
2)\,,$ (79) $\displaystyle F_{4}=$ $\displaystyle\
-\frac{\mu^{3}}{m_{1}^{2}m_{2}}\frac{3(g_{1}-2)}{4l(l+1)(2l+1)}\big{(}24-18g_{2}+9g_{1}g_{2}+16l(l+1)\big{)}$
$\displaystyle+\frac{\mu^{2}}{m_{1}m_{2}}\frac{3}{8l(l+1)(2l+1)}\big{(}48(g_{2}-1)-15g_{1}g_{2}+l(l+1)(32(g_{2}-1)-28g_{1}g_{2})\big{)}$
$\displaystyle+\frac{3\mu^{4}}{2m_{1}^{2}m_{2}^{2}}\bigg{(}\frac{48-72g_{1}g_{2}+15g_{1}^{2}g_{2}^{2}+l(l+1)(12g_{1}^{2}g_{2}^{2}-64)}{16l(l+1)(2l+1)}-\frac{3\,g_{1}^{2}g_{2}^{2}}{4\lambda_{0}}\bigg{)}+(1\leftrightarrow
2)\,,$ (80) $\displaystyle F_{3}=$ $\displaystyle\
\frac{\mu^{3}}{4m_{1}^{2}m_{2}}(g_{1}-2)\big{(}8\lambda_{8}+9(g_{1}-2)g_{2}\lambda_{9}\big{)}+\frac{\mu^{2}}{8m_{1}m_{2}}\bigg{(}\frac{300g_{1}g_{2}}{\lambda_{0}(2l+1)}-g_{1}(7g_{2}-4)\lambda_{10}-8(g_{1}-2)\lambda_{8}+36(g_{1}-1)g_{2}\lambda_{9}\bigg{)}$
$\displaystyle+\frac{\mu^{4}}{2m_{1}^{2}m_{2}^{2}}\bigg{(}-\lambda_{10}-\frac{12}{\lambda_{0}(2l+1)}+\frac{9}{2}\,g_{1}g_{2}\bigg{(}\lambda_{9}+\frac{4}{\lambda_{0}(2l+1)}\bigg{)}+\frac{1}{16}g_{1}^{2}g_{2}^{2}\lambda_{11}\bigg{)}+(1\leftrightarrow
2)\,.$ (81)

The above coefficients can be verified against the positronium limit of (75)
by setting $m_{1}=m_{2}$, $Z=1$, $\kappa_{1}=\kappa_{2}=0$ and both particles
being point-like. The obtained result is in agreement with the former one from
Ref. [12]. It is worth mentioning that although we derived the result for the
case of $l>1$ the final formula gives a correct result also for $l=1$.
Moreover, we would like to point out that in Ref. [13] during the evaluation
of Eq. (125) for $\delta E_{5}$ an error was made for the case of $l>1$.
Namely, in the final result for $\delta E^{(6)}(n,j,l,s)$ instead of the value
$l$, the number $1$ was inserted in the part coming from $\delta E_{5}$
affecting results in Eqs. (205), (209), (213), and (217). After correcting
this error, the results of Ref. [13] are in accordance with those of Ref. [14]
as well as with this work.

Exploring now the limiting case of Eq. (75) for two point-like particles when
the particle of mass $m_{1}$ is infinitely heavy, we set $\kappa_{2}=0$ and
obtain the Dirac result

$E^{(6,0)}=m_{2}(Z\alpha)^{6}\,f^{(6)}(n,j)\,,$ (82)

where we used the total angular momentum $\vec{J}=\vec{L}+\vec{s}_{2}$. The
first-order recoil correction given by expanding Eq. (75) for large $m_{1}$ up
to the order $m_{2}^{2}/m_{1}$, using Eq. (60), reads

$\displaystyle
E^{(6,1)}=(Z\alpha)^{6}\,\frac{m_{2}^{2}}{m_{1}}\,\biggl{[}f^{(6,1)}_{1/2}(n,k)$
$\displaystyle+\frac{g_{1}}{2}\bigg{\\{}\vec{s}_{1}\cdot\vec{s}_{2}\bigg{(}\frac{2}{3(2l-1)(2l+1)(2l+3)\,n^{5}}-\frac{2}{l(l+1)(2l+1)^{2}\,n^{4}}-\frac{2(-3-11l+25l^{2}+72l^{3}+36l^{4})}{3l^{2}(l+1)^{2}(2l-1)(2l+1)^{3}(2l+3)\,n^{3}}\bigg{)}$
$\displaystyle+\vec{L}\cdot\vec{s}_{1}\bigg{(}-\frac{-21+32l+32l^{2}}{2l(l+1)(2l-1)(2l+1)(2l+3)\,n^{5}}+\frac{3}{2l^{2}(l+1)^{2}\,n^{4}}+\frac{-3-5l+19l^{2}+48l^{3}+24l^{4}}{2l^{3}(l+1)^{3}(2l-1)(2l+1)(2l+3)\,n^{3}}\bigg{)}$
$\displaystyle+(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}\bigg{(}\frac{-63+116l+116l^{2}}{l(l+1)(2l-1)^{2}(2l+1)(2l+3)^{2}\,n^{5}}-\frac{3(3+20l+20l^{2})}{l^{2}(l+1)^{2}(2l-1)(2l+1)^{2}(2l+3)\,n^{4}}$
$\displaystyle-\frac{-9-75l-115l^{2}+640l^{3}+2120l^{4}+2160l^{5}+720l^{6}}{l^{3}(l+1)^{3}(2l-1)^{2}(2l+1)^{3}(2l+3)^{2}\,n^{3}}\bigg{)}\bigg{\\}}\bigg{]}$
(83)

and can be considered as a new result.

##  VI Summary and conclusions

We have studied two-body systems with constituents having arbitrary mass,
arbitrary magnetic moment, and spin $s=0$ or $s=1/2$, in states with angular
momentum $l>1$. The expansion of energy levels up to $\alpha^{6}\,\mu$ has
been obtained in an analytic form and verified against special cases where one
of the particles is much heavier. The first two expansion terms in the mass
ratio have a universal character. The leading term coincides with the solution
of the Klein-Gordon or Dirac equation, while the term proportional to the mass
ratio does not depend on the spin of the heavier particle. The next expansion
terms are not universal and depend on both particle spins. In addition, the
special case of positronium $m_{1}=m_{2}$, and $\kappa_{1}=\kappa_{2}=0$ was
used to verify the correctness of our results. Although in our calculations we
assumed $l>1$, our results were in agreement also for the case of $l=1$. It is
not proof, but it may indicate that the presented results are valid in the
more general case of $l>0$ for point particles.

The obtained general formulas are valid for many two-body systems, including
muonic and antiprotonic atoms. This demonstrates that precise measurements of
energy levels can serve not only for the determination of fundamental
constants, but also for the search of unknown interactions in the range not
accessible in normal atoms, namely from 1 MeV up to 100 MeV. Although there
are no definite plans to study rotational states of, for example, protonium or
$\bar{p}\,\alpha$, the availability of a high-precision result is crucial for
planning the corresponding measurements.

One such measurement is considered at NIST [6], which involves rotational
states of heavy hydrogen-like ions, with the purpose of an independent
determination of the Rydberg constant. The vacuum polarization is completely
negligible, and formulas obtained in this paper will be valid also for
circular states of heavy ions because the effective expansion parameter for
these states is $(Z\alpha/n)^{6}$. If necessary, inclusion of the $\alpha^{7}$
contribution for an arbitrary mass ratio is feasible along the lines of Ref.
[15]. This means that one may explore the potential of high accuracy results
for an arbitrary two-body system in search for physics beyond the Standard
Model.

###### Acknowledgements.

This work was supported by the National Science Center (Poland) Grant No.
2017/27/B/ST2/02459. JZ would like to thank Gregory Adkins for spotting the
problem with $\delta E_{5}$.

##  Appendix A Derivation of the effective operator $H^{(6)}$

Using the NRQED Hamiltonian in Eq. (28) we derive individual first-order
effective operators of the order $\alpha^{6}$ for the case of two spin-$1/2$
particles with omitting all the local terms. The derivation proceeds along the
same lines as in Ref. [9]. For spinless particles the effective operators are
obtained by omitting the corresponding spin terms.

The first effective operator comes from the fourth term in the NRQED
Hamiltonian $H_{\textrm{NRQED}}$ in Eq. (28) and represents the higher-order
relativistic correction to the kinetic energy,

$\delta H_{0}=\frac{p_{1}^{6}}{16m_{1}^{5}}+\frac{p_{2}^{6}}{16m_{2}^{5}}\,.$
(84)

$\delta H_{1}$ comes from the interaction between particles when one vertex is
given by the sixth and seventh terms in $H_{\textrm{NRQED}}$, and the other
one is $eA_{0}$. It leads to

$\displaystyle\delta H_{1}=$ $\displaystyle\
\sum_{a}\frac{3}{16m_{a}^{4}}\bigl{(}1+\frac{4\kappa_{a}}{3}\bigr{)}\vec{s}_{a}\cdot\\{p_{a}^{2},e_{a}\vec{\cal{E}}_{a}\times\vec{p}_{a}\\}$
$\displaystyle+\frac{1}{32m_{a}^{4}}\bigg{(}1+\frac{s_{a}(s_{a}+1)}{3}\bigg{)}\,[p_{a}^{2},[p_{a}^{2},V]]\,,$
(85)

where the potential $V$ is

$V=\frac{e_{1}e_{2}}{r}=-\frac{Z\alpha}{r}$ (86)

and the electrostatic fields $\vec{\cal{E}}_{a}$ are defined in Eq. (45).

$\delta H_{2}$ arises when both particles interact via the fifth term in
$H_{\textrm{NRQED}}$. It leads to the effective operator

$\displaystyle\delta H_{2}=$ $\displaystyle\
\frac{e_{1}e_{2}}{4m_{1}^{2}m_{2}^{2}}(1+2\kappa_{1})(1+2\kappa_{2})\,$
$\displaystyle(\vec{s}_{2}\times\vec{p}_{2})^{i}\frac{1}{r^{3}}\left(\delta^{ij}-3\frac{r^{i}r^{j}}{r^{2}}\right)(\vec{s}_{1}\times\vec{p}_{1})^{j}\,.$
(87)

$\delta H_{3}$ arises when one of the particles interacts through part of the
second term in $H_{\textrm{NRQED}}$,

$-\frac{e}{m}\,\vec{p}\cdot\vec{A}-\frac{e}{m}\;(1+\kappa)\vec{s}\cdot\vec{B}\,,$
(88)

while the other particle through part of the third term in
$H_{\textrm{NRQED}}$,

$\displaystyle-\frac{1}{8m^{3}}\left(\pi^{4}-2e\vec{s}\cdot\vec{B}\,\pi^{2}-2\pi^{2}e\vec{s}\cdot\vec{B}\right)\,.$
(89)

If we define static vector potentials by Eqs. (43) and (44), then we can write
the effective operator $\delta H_{3}$ as

$\displaystyle\delta
H_{3}=\frac{1}{4}\sum_{a=1,2}\frac{1}{m_{a}^{3}}\big{\\{}p^{2}_{a},\vec{p}_{a}\cdot
e_{a}\vec{\cal{A}}_{a}+\vec{s}_{a}\cdot(\vec{\nabla}_{a}\times
e_{a}\vec{\cal{A}}_{a})\big{\\}}.$ (90)

$\delta H_{4}$ comes from the exchange of two transverse photons between
particles. One of the particles interacts twice by the term (88), and the
other one by $\frac{e^{2}}{2m}\,\vec{A}^{2}$. The resulting effective operator
can then be obtained simply by replacing the vector potential $\vec{A}$ by the
static field $\vec{\cal{A}}$,

$\delta H_{4}=\sum_{a=1,2}\frac{e_{a}^{2}}{2m_{a}}\,\vec{\cal{A}}_{a}^{2}\,.$
(91)

$\delta H_{5}$ corresponds to the case when one of the particles interacts
through the fifth term in $H_{\textrm{NRQED}}$, and the other one through (88)
and $eA_{0}$,

$\delta
H_{5}=\sum_{a=1,2}\,\frac{e_{a}^{2}\,(1+2\kappa_{a})}{4m_{a}^{2}}\,\vec{s}_{a}\cdot(\vec{\cal{E}}_{a}\times\vec{\cal{A}}_{a}-\vec{\cal{A}}_{a}\times\vec{\cal{E}}_{a}).$
(92)

$\delta H_{6}$ represents the retardation correction to the single transverse
photon exchange between particles when both particles interact through (88).
The corresponding $\alpha^{6}$ contribution is

$\displaystyle\langle\delta H_{6}\rangle=\sum_{a\neq
b}\sum_{b}-\frac{e_{a}e_{b}}{m_{a}m_{b}}\int\frac{d^{3}k}{(2\pi)^{3}\,2\,k^{4}}\bigg{(}\delta^{ij}-\frac{k^{i}k^{j}}{k^{2}}\bigg{)}$
$\displaystyle\bigg{\langle}\big{(}\vec{p}_{a}+\vec{s}_{a}\times\vec{\nabla}_{a}\big{)}^{i}\,e^{i\vec{k}\cdot\vec{r}_{a}}(H_{0}-E_{0})^{2}\big{(}\vec{p}_{b}+\vec{s}_{b}\times\vec{\nabla}_{b}\big{)}^{j}\,e^{-i\vec{k}\cdot\vec{r}_{b}}\bigg{\rangle}\,.$
(93)

Commuting out the expression $(H_{0}-E_{0})^{2}$ and performing the
integration we obtain the effective operator $\delta H_{6}$ in Eq. (36).

$\delta H_{7}$ is a retardation correction in a single transverse photon
exchange where one vertex is (88), and the other one comes from the fifth term
in $H_{\textrm{NRQED}}$,

$-\frac{e(1+2\kappa)}{4m^{2}}\;\vec{s}\cdot(\vec{E}\times\vec{p}-\vec{p}\times\vec{E})\,.$
(94)

The resulting effective operator is then

$\displaystyle\delta
H_{7}=\sum_{a=1,2}\frac{e_{a}^{2}\,(1+2\kappa_{a})}{4m_{a}^{2}}\;\vec{s}_{a}\cdot\big{(}\vec{\cal{E}}_{a}\times\vec{\cal{A}}_{a}-\vec{\cal{A}}_{a}\times\vec{\cal{E}}_{a}\big{)}$
$\displaystyle+\frac{ie_{a}\,(1+2\kappa_{a})}{8m_{a}^{3}}\,\big{[}\,\vec{\cal{A}}_{a}\cdot(\vec{p}_{a}\times\vec{s}_{a})+(\vec{p}_{a}\times\vec{s}_{a})\cdot\vec{\cal{A}}_{a},p_{a}^{2}\,\big{]}\,.$
(95)

$\delta H_{8}$ comes as a correction to the Coulomb interaction between
particles. One particle interacts through the eighth term in
$H_{\textrm{NRQED}}$, and the other one by $eA_{0}$. Then

$\displaystyle\delta H_{8}=$ $\displaystyle\
-\frac{1}{2}\sum_{a}\bigg{(}\alpha_{Ea}-\frac{s_{a}(s_{a}+1)}{3\,m_{a}^{3}}\bigg{)}\,\frac{Z^{2}\alpha^{2}}{r^{4}}\,.$
(96)

Finally, the effective operator $\delta H_{9}$ arises when one of the
particles interacts through part of the third term in $H_{\textrm{NRQED}}$,

$\frac{e\,\kappa}{4m^{3}}\\{\vec{\pi}\cdot\vec{B},\vec{\pi}\cdot\vec{s}\\}\,,$
(97)

and the second particle through (88), giving rise to the single transverse
photon exchange without retardation. The resulting effective operator can be
obtained simply by replacing $\vec{B}$ in the last equation by its static
form. Namely, we can write

$\delta
H_{9}=\sum_{a}\frac{e_{a}\,\kappa_{a}}{4m_{a}^{3}}\\{\vec{p}_{a}\cdot\vec{\nabla}_{a}\times\vec{\cal{A}}_{a},\vec{p}_{a}\cdot\vec{s}_{a}\\}\,.$
(98)

This concludes the derivation of effective operators of the order $\alpha^{6}$
for two spin-1/2 particles for states with $l>1$.

##  Appendix B Expectation values for particles with $s_{1}=s_{2}=0$

To simplify the expectation values we set $Z\alpha=1$. The final result then
needs to be multiplied by a factor $(Z\alpha)^{6}$. Let us denote the
individual contributions in Eq. (27) as

$E^{(6)}=\delta_{1}E^{(6)}+\delta_{2}E^{(6)}\,.$ (99)

We start by evaluating the first-order part $\delta E^{(6)}_{1}$ and again
stress that we omitted all the local Dirac delta-like contributions because we
are focusing on states with $l>1$.

#### First-order contribution $\delta_{1}E^{(6)}$.

The first-order correction to the energy is the expectation value of operator
$H^{(6)}$, and for two spinless particles we use the effective operators (42)
omitting all spins. We separate $\delta_{1}E^{(6)}$ into expectation values of
individual operators appearing in $H^{(6)}$; thus, in general we have

$\delta_{1}E^{(6)}=\sum_{i=0}^{9}\delta E_{i}\,,$ (100)

where

$\delta E_{i}=\left\langle\delta H_{i}\right\rangle\,.$ (101)

However, for spinless particles the effective operators $\delta H_{2}$,
$\delta H_{5}$, $\delta H_{7}$, and $\delta H_{9}$ vanish, so the
corresponding contributions to energy are zero,

$\delta E_{2}=\delta E_{5}=\delta E_{7}=\delta E_{9}=0\,.$ (102)

The first nonzero contribution is the higher-order relativistic correction to
kinetic energy, which can be evaluated to

$\displaystyle\delta
E_{0}=\frac{\mu^{2}}{2}\biggl{(}\frac{1}{m_{1}^{5}}+\frac{1}{m_{2}^{5}}\biggr{)}\biggl{(}\mu
E_{0}^{3}+3\mu E_{0}^{2}\left\langle\frac{1}{r}\right\rangle$
$\displaystyle+3\mu
E_{0}\left\langle\frac{1}{r^{2}}\right\rangle+\mu\left\langle\frac{1}{r^{3}}\right\rangle+\frac{1}{2}\left\langle\frac{1}{r^{4}}\right\rangle\biggr{)}\,,$
(103)

where $E_{0}=-\frac{\mu}{2n^{2}}$. The expectation value of operators
$1/r^{k}$ will be resolved with the help of the formulas in Appendix E. The
next corrections are

$\displaystyle\delta E_{1}=$ $\displaystyle\
-\frac{\mu}{8}\biggl{(}\frac{1}{m_{1}^{4}}+\frac{1}{m_{2}^{4}}\biggr{)}\left\langle\frac{1}{r^{4}}\right\rangle\,,$
(104) $\displaystyle\delta E_{3}=$ $\displaystyle\
\frac{\mu}{m_{1}m_{2}}\biggl{(}\frac{1}{m_{1}^{2}}+\frac{1}{m_{2}^{2}}\biggr{)}\bigg{\\{}2\mu
E_{0}^{2}\left\langle\frac{1}{r}\right\rangle+4\mu
E_{0}\left\langle\frac{1}{r^{2}}\right\rangle$ $\displaystyle+\bigl{(}2\mu-
E_{0}l(l+1)\bigr{)}\left\langle\frac{1}{r^{3}}\right\rangle+\frac{1}{2}\bigl{(}1-l(l+1)\bigr{)}\left\langle\frac{1}{r^{4}}\right\rangle\bigg{\\}}\,,$
(105) $\displaystyle\delta E_{4}=$ $\displaystyle\
\frac{E_{0}}{m_{1}m_{2}}\left\langle\frac{1}{r^{2}}\right\rangle+\frac{1}{m_{1}m_{2}}\left\langle\frac{1}{r^{3}}\right\rangle$
$\displaystyle+\frac{1}{2\mu\,m_{1}m_{2}}\biggl{[}1-\frac{3l(l+1)}{4}\biggr{]}\left\langle\frac{1}{r^{4}}\right\rangle\,.$
(106)

The retardation correction due to single transverse photon exchange is

$\displaystyle\delta
E_{6}=-\frac{E_{0}}{2m_{1}m_{2}}\left\langle\frac{1}{r^{2}}\right\rangle+\frac{1}{m_{1}m_{2}}\biggl{[}-\frac{1}{4}+\frac{\mu
E_{0}(1+l(l+1))}{m_{1}m_{2}}\biggr{]}\left\langle\frac{1}{r^{3}}\right\rangle$
$\displaystyle+\frac{1}{m_{1}m_{2}}\biggl{[}-\frac{1}{8\mu}+\frac{\mu}{m_{1}m_{2}}+l(l+1)\biggl{(}\frac{3}{8\mu}+\frac{\mu}{m_{1}m_{2}}\biggr{)}\biggr{]}\left\langle\frac{1}{r^{4}}\right\rangle+\frac{3}{2m_{1}^{2}m_{2}^{2}}\biggl{(}1-\frac{l^{2}(l+1)^{2}}{4}\biggr{)}\left\langle\frac{1}{r^{5}}\right\rangle\,.$
(107)

Finally, correction $\delta E_{8}$ is

$\delta
E_{8}=-\frac{1}{2}\,(\alpha_{E1}+\alpha_{E2})\,\left\langle\frac{1}{r^{4}}\right\rangle\,,$
(108)

which concludes the evaluation of the first-order expectation values for two
spinless particles.

#### Second-order contribution $\delta_{2}E^{(6)}$.

The second-order contribution of the order $(Z\alpha)^{6}$ is induced by the
Breit Hamiltonian $H^{(4)}$,

$\delta_{2}E^{(6)}=\langle
H^{(4)}\frac{1}{(E_{0}-H_{0})^{\prime}}\,H^{(4)}\rangle\,.$ (109)

In general, such a contribution contains singularities that must be isolated
by means of regularization and then cancelled when combined with the
corresponding first-order contribution. All such singularities are, however,
proportional to the Dirac delta function and thus vanish for $l>0$. Therefore,
we can ignore these singularities and write the Breit Hamiltonian for spinless
particles in the form (omitting local terms)

$H^{(4)}=-A\bigg{[}\frac{1}{r}\,H_{0}+\,H_{0}\,\frac{1}{r}\bigg{]}-\frac{B}{r^{2}}+\frac{C}{r^{3}}\,,$
(110)

where

$\displaystyle A=$ $\displaystyle\
\frac{\mu}{m_{1}m_{2}}+\frac{\mu^{2}}{2}\biggl{(}\frac{1}{m_{1}^{3}}+\frac{1}{m_{2}^{3}}\biggr{)}\,,$
(111) $\displaystyle B=$ $\displaystyle\
\frac{2\mu}{m_{1}m_{2}}+\frac{\mu^{2}}{2}\biggl{(}\frac{1}{m_{1}^{3}}+\frac{1}{m_{2}^{3}}\biggr{)}\,,$
(112) $\displaystyle C=$ $\displaystyle\ \frac{l(l+1)}{2m_{1}m_{2}}\,.$ (113)

After elementary calculations we can express the second-order contribution in
terms of the first-order operators $\left\langle 1/r^{k}\right\rangle$ listed
in Appendix E and $\mathcal{E}_{ij}$ listed in Appendix F. Then,

$\displaystyle\delta_{2}E^{(6)}=4A^{2}E_{0}\left\langle\frac{1}{r}\right\rangle^{2}-4A^{2}E_{0}\left\langle\frac{1}{r^{2}}\right\rangle+2AB\left\langle\frac{1}{r}\right\rangle\left\langle\frac{1}{r^{2}}\right\rangle$
$\displaystyle-2AB\left\langle\frac{1}{r^{3}}\right\rangle-2AC\left\langle\frac{1}{r}\right\rangle\left\langle\frac{1}{r^{3}}\right\rangle+\,\biggl{(}2AC-\frac{A^{2}}{2\mu}\biggr{)}\left\langle\frac{1}{r^{4}}\right\rangle$
$\displaystyle+4A^{2}E_{0}^{2}\,\mathcal{E}_{11}+4ABE_{0}\,\mathcal{E}_{12}+B^{2}\,\mathcal{E}_{22}-4ACE_{0}\,\mathcal{E}_{13}$
$\displaystyle-2BC\,\mathcal{E}_{23}+\,C^{2}\,\mathcal{E}_{33}\,.$ (114)

Combining the first-order contribution $\delta_{1}E^{(6)}$ with the second-
order term $\delta_{2}E^{(6)}$, and using the expectation values from
Appendices E and F we obtain the total correction to the energy $E^{(6)}$ for
two spinless particles, presented in Eq. (47).

##  Appendix C Expectation values for particles with $s_{1}=0,\,s_{2}=1/2$

The next case is when one particle is spinless and the other one has spin
$s_{2}=1/2$. In contrast to the previous section the effective operators
$\delta H_{5}$ and $\delta H_{7}$ do not vanish, and operators $\delta H_{1}$,
$\delta H_{3}$, $\delta H_{4}$, and $\delta H_{6}$, and also the Breit
Hamiltonian $H^{(4)}$ now have spin-dependent parts.

#### First-order contribution $\delta_{1}E^{(6)}$.

The first-order contribution is given by the sum (100) with the individual
contributions resolved as follows. Contribution $\delta E_{0}=\langle\delta
H_{0}\rangle$ is the same as in the spinless case, given by Eq. (B). The next
term differs from the spinless case and is

$\displaystyle\delta E_{1}=$ $\displaystyle\ \delta
E_{1}^{s_{1}=s_{2}=0}-\mu\biggl{[}\frac{1}{32m_{2}^{4}}+\frac{\vec{L}\cdot\vec{s}_{2}}{m_{2}^{4}}\biggl{(}\frac{3}{4}+\kappa_{2}\biggr{)}\biggr{]}\left\langle\frac{1}{r^{4}}\right\rangle$
$\displaystyle-\frac{\mu
E_{0}}{m_{2}^{4}}\biggl{(}\frac{3}{4}+\kappa_{2}\biggr{)}\vec{L}\cdot\vec{s}_{2}\left\langle\frac{1}{r^{3}}\right\rangle\,.$
(115)

Here $\delta E_{1}^{s_{1}=s_{2}=0}$ stands for the spinless result (104). We
will use this notation throughout the following sections.

The next term $\delta E_{2}=0$ again vanishes. The remaining nonvanishing
contributions are

$\displaystyle\delta E_{3}=$ $\displaystyle\ \delta E_{3}^{s_{1}=s_{2}=0}$
$\displaystyle-\frac{\mu}{m_{1}m_{2}}E_{0}\biggl{(}\frac{(1+\kappa_{2})}{m_{1}^{2}}+\frac{1}{m_{2}^{2}}\biggr{)}\vec{L}\cdot\vec{s}_{2}\left\langle\frac{1}{r^{3}}\right\rangle$
$\displaystyle-\frac{\mu}{m_{1}m_{2}}\biggl{(}\frac{(1+\kappa_{2})}{m_{1}^{2}}+\frac{1}{m_{2}^{2}}\biggr{)}\vec{L}\cdot\vec{s}_{2}\left\langle\frac{1}{r^{4}}\right\rangle\,,$
(116) $\displaystyle\delta E_{4}=$ $\displaystyle\ \delta
E_{4}^{s_{1}=s_{2}=0}+\frac{1}{m_{1}m_{2}}\biggl{[}\frac{1}{4m_{2}}\,(1+\kappa_{2})^{2}$
$\displaystyle-\frac{(1+\kappa_{2})}{2m_{2}}\,\vec{L}\cdot\vec{s}_{2}\biggr{]}\left\langle\frac{1}{r^{4}}\right\rangle\,,$
(117) $\displaystyle\delta E_{5}=$ $\displaystyle\
-\frac{(1+2\kappa_{2})}{4m_{1}m_{2}^{2}}\,\vec{L}\cdot\vec{s}_{2}\left\langle\frac{1}{r^{4}}\right\rangle\,,$
(118) $\displaystyle\delta E_{6}=$ $\displaystyle\
E_{6}^{s_{1}=s_{2}=0}+\frac{\mu}{2m_{1}m_{2}^{3}}\,(1+\kappa_{2})\,\vec{L}\cdot\vec{s}_{2}\left\langle\frac{1}{r^{4}}\right\rangle\,,$
(119) $\displaystyle\delta E_{7}=$ $\displaystyle\
-\frac{(1+2\kappa_{2})}{4m_{1}m_{2}^{2}}\,\biggl{(}1+\frac{\mu}{m_{2}}\biggr{)}\,\vec{L}\cdot\vec{s}_{2}\left\langle\frac{1}{r^{4}}\right\rangle\,.$
(120)

Contribution $\delta E_{8}$ can be obtained from the spinless case (108) by
transition $\alpha_{E2}\rightarrow\alpha_{E2}-(4m_{2}^{3})^{-1}$, and $\delta
E_{9}=0$.

#### Second-order contribution $\delta_{2}E^{(6)}$.

For the second-order contribution we repeat the derivation from the spinless
case with one difference. The Breit Hamiltonian is again expressed in the form
(110) with the coefficients $A$ and $B$ being the same as in the spinless
case. However, coefficient $C$ now contains a spin-dependent part. The
derivation proceeds along the same line leading to the result (B) where for
the coefficient $C$ we obtain

$C=\frac{\vec{L}^{2}}{2m_{1}m_{2}}+\biggl{(}\frac{1+2\kappa_{2}}{2m_{2}^{2}}+\frac{1+\kappa_{2}}{m_{1}m_{2}}\biggr{)}\vec{L}\cdot\vec{s}_{2}$
(121)

instead of (113). To evaluate the term $(\vec{L}\cdot\vec{s}_{2})^{2}$ we use
Eq. (162).

##  Appendix D Expectation values for particles with $s_{1}=s_{2}=1/2$

Finally we evaluate the most complicated case, i.e., two spin-1/2 particles,
such as in antiprotonic hydrogen or positronium. For the first-order
contribution we will use the effective operators (30-42) with both spins
nonzero.

#### First-order contribution $\delta_{1}E^{(6)}$.

The contribution to the energy of the order $(Z\alpha)^{6}$ for states with
$l>1$ coming from the first-order operators is again expressed in the form of
series (100). The effective operators $\delta H_{i}$ may contain spin-
dependent terms from both particles. The first contribution $\delta E_{0}$ is
the same as in the spinless case, given by Eq. (B). The next contribution is

$\displaystyle\delta E_{1}=\delta
E_{1}^{s_{1}=s_{2}=0}-\mu\biggl{[}\frac{1}{32m_{1}^{4}}+\frac{1}{32m_{2}^{4}}$
$\displaystyle+\frac{\vec{L}\cdot\vec{s}_{1}}{m_{1}^{4}}\biggl{(}\frac{3}{4}+\kappa_{1}\biggr{)}+\frac{\vec{L}\cdot\vec{s}_{2}}{m_{2}^{4}}\biggl{(}\frac{3}{4}+\kappa_{2}\biggr{)}\biggr{]}\left\langle\frac{1}{r^{4}}\right\rangle$
$\displaystyle-\mu
E_{0}\bigg{[}\frac{1}{m_{1}^{4}}\biggl{(}\frac{3}{4}+\kappa_{1}\biggr{)}\vec{L}\cdot\vec{s}_{1}+\frac{1}{m_{2}^{4}}\biggl{(}\frac{3}{4}+\kappa_{2}\biggr{)}\vec{L}\cdot\vec{s}_{2}\bigg{]}\left\langle\frac{1}{r^{3}}\right\rangle\,.$

Correction $\delta E_{2}$ does not vanish now and is

$\displaystyle\delta
E_{2}=\frac{1}{4m_{1}^{2}m_{2}^{2}}(1+2\kappa_{1})(1+2\kappa_{2})\bigg{[}\bigg{(}4\mu
E_{0}\left\langle\frac{1}{r^{3}}\right\rangle$
$\displaystyle+4\mu\left\langle\frac{1}{r^{4}}\right\rangle-12(l-1)(l+2)\left\langle\frac{1}{r^{5}}\right\rangle\bigg{)}\frac{(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}}{(2l-1)(2l+3)}$
$\displaystyle+\bigg{(}4\mu
E_{0}\left\langle\frac{1}{r^{3}}\right\rangle+4\mu\left\langle\frac{1}{r^{4}}\right\rangle-3(l-1)(l+2)\left\langle\frac{1}{r^{5}}\right\rangle\bigg{)}$
$\displaystyle\times\frac{\vec{s}_{1}\cdot\vec{s}_{2}}{3}\bigg{]}\,,$ (123)

where we have used the expectation value identity

$\bigg{\langle}\delta^{ij}-3\frac{r^{i}r^{j}}{r^{2}}\bigg{\rangle}=\frac{6}{(2l-1)(2l+3)}\big{\langle}(L^{i}L^{j})^{(2)}\big{\rangle}\,.$
(124)

The next terms are

$\displaystyle\delta E_{3}=$ $\displaystyle\ \delta
E_{3}^{s_{1}=s_{2}=0}+\frac{\mu}{m_{1}m_{2}}\bigg{[}\bigg{(}-E_{0}\,\left\langle\frac{1}{r^{3}}\right\rangle-\left\langle\frac{1}{r^{4}}\right\rangle\bigg{)}\bigg{(}\bigg{(}\frac{1+\kappa_{2}}{m_{1}^{2}}+\frac{1}{m_{2}^{2}}\bigg{)}\vec{L}\cdot\vec{s}_{2}+\bigg{(}\frac{1+\kappa_{1}}{m_{2}^{2}}+\frac{1}{m_{1}^{2}}\bigg{)}\vec{L}\cdot\vec{s}_{1}\bigg{)}$
$\displaystyle+6\bigg{(}\frac{1+\kappa_{2}}{m_{1}^{2}}+\frac{1+\kappa_{1}}{m_{2}^{2}}\bigg{)}\left\langle\frac{1}{r^{3}}\bigg{(}E_{0}+\frac{1}{r}\bigg{)}\frac{(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}}{(2l-1)(2l+3)}\right\rangle\bigg{]}$
(125) $\displaystyle\delta E_{4}=$ $\displaystyle\ \delta
E_{4}^{s_{1}=s_{2}=0}+\frac{1}{m_{1}m_{2}}\biggl{[}\frac{1}{4m_{1}}\,(1+\kappa_{1})^{2}+\frac{1}{4m_{2}}\,(1+\kappa_{2})^{2}-\frac{(1+\kappa_{1})}{2m_{1}}\,\vec{L}\cdot\vec{s}_{1}-\frac{(1+\kappa_{2})}{2m_{2}}\,\vec{L}\cdot\vec{s}_{2}\biggr{]}\left\langle\frac{1}{r^{4}}\right\rangle\,,$
(126) $\displaystyle\delta E_{5}=$ $\displaystyle\
-\frac{(1+2\kappa_{1})}{4m_{1}^{2}m_{2}}\,\bigg{[}\vec{L}\cdot\vec{s}_{1}-4(1+\kappa_{2})\bigg{(}\frac{(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}}{(2l-1)(2l+3)}+\frac{\vec{s}_{1}\cdot\vec{s}_{2}}{3}\bigg{)}\bigg{]}\left\langle\frac{1}{r^{4}}\right\rangle+(1\leftrightarrow
2)\,,$ (127) $\displaystyle\delta E_{6}=$ $\displaystyle\ \delta
E_{6}^{s_{1}=s_{2}=0}+\frac{\mu}{m_{1}m_{2}}\biggl{[}\frac{1}{2m^{2}_{1}}\,(1+\kappa_{1})\,\vec{L}\cdot\vec{s}_{1}+\frac{1}{2m^{2}_{2}}\,(1+\kappa_{2})\,\vec{L}\cdot\vec{s}_{2}\biggr{]}\left\langle\frac{1}{r^{4}}\right\rangle$
$\displaystyle-\frac{\mu}{m_{1}^{2}m_{2}^{2}}(1+\kappa_{1})(1+\kappa_{2})\bigg{(}\frac{2}{3}\vec{s}_{1}\cdot\vec{s}_{2}-\frac{(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}}{(2l-1)(2l+3)}\bigg{)}\left\langle\frac{1}{r^{4}}\right\rangle\,,$
(128) $\displaystyle\delta E_{7}=$ $\displaystyle\
-\frac{\mu\,(1+2\kappa_{1})}{4m_{1}^{3}m_{2}}\bigg{[}\vec{L}\cdot\vec{s}_{1}+4(1+\kappa_{2})\bigg{(}\frac{(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}}{(2l-1)(2l+3)}+\frac{\vec{s}_{1}\cdot\vec{s}_{2}}{3}\bigg{)}\bigg{]}\left\langle\frac{1}{r^{4}}\right\rangle+(1\leftrightarrow
2)+\delta E_{5}\,.$ (129)

Contribution $\delta E_{8}$ is obtained from the spinless case (108) by
transition $\alpha_{Ej}\rightarrow\alpha_{Ej}-(4m_{j}^{3})^{-1}$. Finally,

$\displaystyle\delta E_{9}=$ $\displaystyle\
\frac{\kappa_{1}(1+\kappa_{2})}{2m_{1}^{3}m_{2}}\bigg{[}\bigg{(}8\mu
E_{0}\left\langle\frac{1}{r^{3}}\right\rangle+8\mu\left\langle\frac{1}{r^{4}}\right\rangle-6(l-1)(l+2)\left\langle\frac{1}{r^{5}}\right\rangle\bigg{)}\frac{(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}}{(2l-1)(2l+3)}$
$\displaystyle+\bigg{(}-\frac{4}{3}\mu
E_{0}\left\langle\frac{1}{r^{3}}\right\rangle-\frac{4}{3}\mu\left\langle\frac{1}{r^{4}}\right\rangle+(l-1)(l+2)\left\langle\frac{1}{r^{5}}\right\rangle\bigg{)}\vec{s}_{1}\cdot\vec{s}_{2}\bigg{]}+(1\leftrightarrow
2)\,.$ (130)

The total first-order correction $\delta_{1}E^{(6)}$ is a sum
$\sum_{i=0}^{9}\delta E_{i}$.

#### Second-order contribution $\delta_{2}E^{(6)}$.

The second-order contribution in the case of two particles with spin 1/2 is
more complicated than in the previous cases due the to presence of a tensor-
like term in the Breit Hamiltonian. We write the Hamiltonian $H^{(4)}$ in the
form

$\displaystyle
H^{(4)}=-\frac{\mu^{2}}{2}\bigg{(}\frac{1}{m_{1}^{3}}+\frac{1}{m_{2}^{3}}\bigg{)}H_{0}^{2}-\frac{\mu}{r^{2}}\bigg{[}\frac{\mu}{2}\bigg{(}\frac{1}{m_{1}^{3}}+\frac{1}{m_{2}^{3}}\bigg{)}+\frac{2}{m_{1}m_{2}}\bigg{]}-\mu\bigg{[}\frac{\mu}{2}\bigg{(}\frac{1}{m_{1}^{3}}+\frac{1}{m_{2}^{3}}\bigg{)}+\frac{1}{m_{1}m_{2}}\bigg{]}\bigg{(}\frac{1}{r}H_{0}+H_{0}\frac{1}{r}\bigg{)}$
$\displaystyle+\frac{1}{r^{3}}\bigg{[}\bigg{(}\frac{1+\kappa_{1}}{m_{1}m_{2}}+\frac{1+2\kappa_{1}}{2m_{1}^{2}}\bigg{)}\vec{L}\cdot\vec{s}_{1}+\bigg{(}\frac{1+\kappa_{2}}{m_{1}m_{2}}+\frac{1+2\kappa_{2}}{2m_{2}^{2}}\bigg{)}\vec{L}\cdot\vec{s}_{2}+\frac{\vec{L}^{2}}{2m_{1}m_{2}}\bigg{]}+\frac{3}{m_{1}m_{2}}(1+\kappa_{1})(1+\kappa_{2})\frac{(n^{i}n^{j})^{(2)}s^{i}_{1}s^{j}_{2}}{r^{3}}\,,$

with $(n^{i}n^{j})^{(2)}=n^{i}n^{j}-\delta^{ij}/3$ being a symmetric,
traceless tensor. The evaluation of the second-order contribution with the
Breit Hamiltonian in the form (D) proceeds along the same line as in the
previous cases using identities in Appendices F and G with the exception of
the term

$\bigg{\langle}\frac{(n^{i}n^{j})^{(2)}s^{i}_{1}s^{j}_{2}}{r^{3}}\frac{1}{(E_{0}-H_{0})^{\prime}}\frac{(n^{k}n^{l})^{(2)}s^{k}_{1}s^{l}_{2}}{r^{3}}\bigg{\rangle}\,.$
(132)

There are three possible contributions that correspond to three possible
values of the angular momentum of the intermediate state $l^{\prime}=l,l\pm
2$. Let us denote by $Q=(n^{i}n^{j})^{(2)}s_{1}^{i}s_{2}^{j}$ and by
$\hat{P}_{l}$ the projection operator into the subspace with angular momentum
$l$. Assuming that the outer states have the angular momentum $l$, we get for
the diagonal term

$\displaystyle Q\,\hat{P}_{l}\,Q=$ $\displaystyle\
\frac{l(l+1)}{24(2l-1)(2l+3)}\bigg{(}1+\frac{4}{3}\vec{s}_{1}\cdot\vec{s}_{2}\bigg{)}$
$\displaystyle\
-\frac{1}{8(2l-1)(2l+3)}\vec{L}\cdot(\vec{s}_{1}+\vec{s}_{2}\big{)}$
$\displaystyle-\frac{(2l-3)(2l+5)}{6(2l-1)^{2}(2l+3)^{2}}(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}\,.$
(133)

For off-diagonal terms we get

$\displaystyle Q\,\hat{P}_{l+2}\,Q=$ $\displaystyle\
\frac{(l+1)(l+2)}{16(2l+1)(2l+3)}\bigg{(}1+\frac{4}{3}\vec{s}_{1}\cdot\vec{s}_{2}\bigg{)}$
$\displaystyle+\frac{(l+2)}{8(2l+1)(2l+3)}\vec{L}\cdot\big{(}\vec{s}_{1}+\vec{s}_{2}\big{)}$
$\displaystyle+\frac{(l+2)}{2(2l+1)(2l+3)^{2}}(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}\,,$
(134)

and

$\displaystyle Q\,\hat{P}_{l-2}\,Q=$ $\displaystyle\
\frac{l(l-1)}{16(2l-1)(2l+1)}\bigg{(}1+\frac{4}{3}\vec{s}_{1}\cdot\vec{s}_{2}\bigg{)}$
$\displaystyle-\frac{(l-1)}{8(2l-1)(2l+1)}\vec{L}\cdot\big{(}\vec{s}_{1}+\vec{s}_{2}\big{)}$
$\displaystyle+\frac{(l-1)}{2(2l-1)^{2}(2l+1)}(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}\,.$
(135)

The radial part of the second-order contribution (132) is, in the case of the
diagonal term $l^{\prime}=l$, equal to $\mathcal{E}_{33}$, Eq. (F). For off-
diagonal terms the Eqs. (F) and (F) are used.

##  Appendix E Expectation values of the basic first-order operators

Here we list the expectation values used in this work. Derivation of the
expectation values for operators $1/r^{k}$ can be found in the third section
of the first chapter of book [16], and setting $Z\alpha=1$ they are

$\displaystyle\left\langle\frac{1}{r}\right\rangle=\frac{\mu}{n^{2}},$ (136)
$\displaystyle\left\langle\frac{1}{r^{2}}\right\rangle=\frac{2\mu^{2}}{\left(2l+1\right)n^{3}},$
(137)
$\displaystyle\left\langle\frac{1}{r^{3}}\right\rangle=\frac{2\mu^{3}}{l(l+1)(2l+1)n^{3}},$
(138)
$\displaystyle\left\langle\frac{1}{r^{4}}\right\rangle=4\mu^{4}\frac{3n^{2}-l(1+l)}{l(l+1)(2l-1)(2l+1)(2l+3)n^{5}},$
(139)
$\displaystyle\left\langle\frac{1}{r^{5}}\right\rangle=\frac{4\mu^{5}\big{(}1-3l(1+l)+5n^{2}\big{)}}{(l-1)l(l+1)(l+2)(2l-1)(2l+1)(3+2l)n^{5}}.$

These equations are valid for $l>1$.

##  Appendix F The second-order matrix elements

In this section we present the expectation values of the second-order
operators of the form

$\displaystyle\mathcal{E}_{nm}=\left\langle\frac{1}{r^{n}}\,\frac{1}{(E_{o}-H_{0})^{\prime}}\,\frac{1}{r^{m}}\right\rangle\,.$
(141)

We again set $Z\alpha=1$ for simplification of the results. When the angular
momentum of the intermediate states is the same as for the reference state
then we obtain the following results:

$\displaystyle\mathcal{E}_{11}=-\frac{\mu}{2n^{2}},$ (142)
$\displaystyle\mathcal{E}_{12}=-\frac{2\mu^{2}}{(2l+1)n^{3}},$ (143)
$\displaystyle\mathcal{E}_{22}=-\mu^{3}\left(\frac{6}{(2l+1)^{2}\,n^{4}}+\frac{4}{(2l+1)^{3}\,n^{3}}\right),$
(144)
$\displaystyle\mathcal{E}_{13}=\frac{\mu}{2l(l+1)}\left(2\,\mathcal{E}_{12}-\left\langle\frac{1}{r^{2}}\right\rangle\right)\,,$
(145)
$\displaystyle\mathcal{E}_{23}=\frac{\mu}{l(l+1)}\left(\,\mathcal{E}_{22}-\left\langle\frac{1}{r^{3}}\right\rangle\right)\,,$
(146)
$\displaystyle\mathcal{E}_{33}=\frac{\mu^{2}}{l^{2}(l+1)^{2}}\,\mathcal{E}_{22}-\frac{\mu^{2}}{l^{2}(l+1)^{2}}\left\langle\frac{1}{r^{3}}\right\rangle$
$\displaystyle-\frac{3\mu}{2l(l+1)}\left\langle\frac{1}{r^{4}}\right\rangle\,.$
(147)

The second-order expectation values involving the $1/r^{3}$ operator were
derived with the help of the identity

$\frac{1}{r^{3}}=\frac{\mu}{l(l+1)}\left(\frac{1}{r^{2}}-\left[D_{r},H_{0}\right]\right)\,,$
(148)

where $D_{r}=\partial_{r}+\frac{1}{r}$. For off-diagonal terms with $1/r^{3}$
the evaluation is more complicated. In the case of $l^{\prime}=l+2$ we need to
calculate the expression

$\bigg{\langle}\frac{1}{r^{3}}\frac{1}{(E_{0}-H_{l+2})^{\prime}}\frac{1}{r^{3}}\bigg{\rangle}$
(149)

in which we defined the radial part of the Hamiltonian as

$H_{l}=-\frac{1}{2\mu}\bigg{(}\partial_{r}^{2}+\frac{2}{r}\partial_{r}-\frac{l(l+1)}{r^{2}}\bigg{)}-\frac{1}{r}\,.$
(150)

Using the identities

$\displaystyle\frac{1}{r^{3}}=$ $\displaystyle\
H_{l}\hat{\alpha}-\hat{\alpha}H_{l+2}\,,$ (151)
$\displaystyle\frac{1}{r^{3}}=$ $\displaystyle\
\hat{\beta}H_{l}-H_{l+2}\hat{\beta}\,,$ (152)

where operators $\hat{\alpha}$ and $\hat{\beta}$ are defined by

$\displaystyle\hat{\alpha}=$ $\displaystyle\
-\frac{\mu}{3(l+1)(l+2)}D_{r}-\frac{\mu(2l+3)}{3(l+1)(l+2)}\frac{1}{r}\,$
$\displaystyle+\frac{\mu^{2}}{3(2l+3)(l+1)(l+2)}$ (153)
$\displaystyle\hat{\beta}=$ $\displaystyle\
\frac{\mu}{3(l+1)(l+2)}D_{r}-\frac{\mu(2l+3)}{3(l+1)(l+2)}\frac{1}{r}\,$
$\displaystyle+\frac{\mu^{2}}{3(2l+3)(l+1)(l+2)}\,,$ (154)

we can simplify the radial part to

$\displaystyle\bigg{\langle}\frac{1}{r^{3}}\frac{1}{(E_{0}-H_{l+2})^{\prime}}\frac{1}{r^{3}}\bigg{\rangle}=\langle\hat{\alpha}\,(E_{0}-H_{l+2})\,\hat{\beta}\rangle$
$\displaystyle=-\frac{2\mu^{2}(2l+3)}{9(l+1)^{2}(l+2)^{2}}E_{0}\bigg{\langle}\frac{1}{r^{2}}\bigg{\rangle}$
$\displaystyle-\frac{\mu^{2}(7+9l+3l^{2})}{3(l+1)^{2}(l+2)^{2}(2l+3)}\bigg{\langle}\frac{1}{r^{3}}\bigg{\rangle}$
$\displaystyle-\frac{\mu(27+45l+25l^{2}+4l^{3})}{18(l+1)^{2}(l+2)^{2}}\bigg{\langle}\frac{1}{r^{4}}\bigg{\rangle}\,$
(155)

For the radial part with the angular momentum of the intermediate states equal
to $l^{\prime}=l-2$ we repeat the calculation using

$\displaystyle\frac{1}{r^{3}}=$ $\displaystyle\
H_{l}\hat{\delta}-\hat{\delta}H_{l-2}\,,$ (156)
$\displaystyle\frac{1}{r^{3}}=$ $\displaystyle\
\hat{\gamma}H_{l}-H_{l-2}\hat{\gamma}\,,$ (157)

with

$\displaystyle\hat{\delta}=$ $\displaystyle\
-\frac{\mu}{3l(l-1)}D_{r}-\frac{\mu(1-2l)}{3l(l-1)}\frac{1}{r}-\frac{\mu^{2}}{3l(l-1)(2l-1)}\,,$
$\displaystyle\hat{\gamma}=$ $\displaystyle\
\frac{\mu}{3l(l-1)}D_{r}-\frac{\mu(1-2l)}{3l(l-1)}\frac{1}{r}-\frac{\mu^{2}}{3l(l-1)(2l-1)}\,.$

Then the radial part is

$\displaystyle\bigg{\langle}\frac{1}{r^{3}}\frac{1}{(E_{0}-H_{l-2})^{\prime}}\frac{1}{r^{3}}\bigg{\rangle}=\langle\hat{\delta}\,(E_{0}-H_{l-2})\,\hat{\gamma}\rangle$
$\displaystyle=\frac{2\mu^{2}(2l-1)}{9(l-1)^{2}l^{2}}E_{0}\bigg{\langle}\frac{1}{r^{2}}\bigg{\rangle}+\frac{\mu^{2}(1-3l+3l^{2})}{3(l-1)^{2}l^{2}(2l-1)}\bigg{\langle}\frac{1}{r^{3}}\bigg{\rangle}$
$\displaystyle+\frac{\mu(-3+7l-13l^{2}+4l^{3})}{18(l-1)^{2}l^{2}}\bigg{\langle}\frac{1}{r^{4}}\bigg{\rangle}\,.$
(160)

Eqs. (F) and (F) are used for evaluation of the second-order contribution for
two spin-1/2 particles.

##  Appendix G Reduction of spin-angular operators

To calculate the second-order contribution in $E^{(6)}$ we use the following
identities to evaluate products of spin-angular operators:

$\displaystyle(\vec{L}\cdot\vec{s}_{1})^{2}=$
$\displaystyle\frac{l(l+1)}{4}-\frac{1}{2}\vec{L}\cdot\vec{s}_{1}\,,$ (161)
$\displaystyle(\vec{L}\cdot\vec{s}_{2})^{2}=$
$\displaystyle\frac{l(l+1)}{4}-\frac{1}{2}\vec{L}\cdot\vec{s}_{2}\,,$ (162)
$\displaystyle\frac{1}{2}\big{\\{}\vec{L}\cdot\vec{s}_{1},\vec{L}\cdot\vec{s}_{2}\big{\\}}=$
$\displaystyle(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}+\frac{l(l+1)}{3}\vec{s}_{1}\cdot\vec{s}_{2}\,,$
$\displaystyle\frac{1}{2}\big{\\{}(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j},\vec{L}\cdot\vec{s}_{1}\big{\\}}=$
$\displaystyle\frac{(2l-1)(2l+3)}{24}\vec{L}\cdot\vec{s}_{2}$
$\displaystyle-\frac{3}{4}(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j},$ (164)
$\displaystyle\frac{1}{2}\big{\\{}(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j},\vec{s}_{1}\cdot\vec{s}_{2}\big{\\}}=$
$\displaystyle\frac{1}{4}(L^{i}L^{j})^{(2)}s_{1}^{i}s_{2}^{j}\,$ (165)
$\displaystyle(\vec{s}_{1}\cdot\vec{s}_{2})^{2}=$
$\displaystyle\frac{3}{16}-\frac{1}{2}\vec{s}_{1}\cdot\vec{s}_{2}\,.$ (166)

The second-order contribution is then expressed as a linear combination of the
basic spin-angular operators in the same way as the first-order contribution.

## References

  * [1] B. Ohayon, G. Janka, I. Cortinovis, Z. Burkley, L. d. S. Borges, E. Depero, A. Golovizin, X. Ni, Z. Salman, A. Suter, C. Vigo, T. Prokscha, and P. Crivelli (Mu-MASS Collaboration), Phys. Rev. Lett. 128, 011802 (2022).
  * [2] H. Grotch and D. R. Yennie, Z. Phys 202, 425 (1967); Rev. Mod. Phys. 41, 350 (1969).
  * [3] Eberhard Widmann, Few-Body Systems 45, 165 (2009).
  * [4] R. S. Hayano, M. Hori, D. Horvath and E. Widmann, Rep. Prog. Phys. 70, 1995 (2007).
  * [5] M. Hori et al., Phys. Rev. Lett. 94, 063401 (2005).
  * [6] J. N. Tan, S. M. Brewer, and N. D. Guise, Phys. Scr. 2011, 014009 (2011).
  * [7] A. Veitia and K. Pachucki, Phys. Rev. A 69, 042501 (2004).
  * [8] J. Zatorski and K. Pachucki, Phys. Rev. A 82, 052520 (2010).
  * [9] K. Pachucki, Phys. Rev. A 71, 012503 (2005).
  * [10] E. A. Golosov, I. B. Khriplovich, A. I. Milstein and A. S. Yelkhovsky, Zh. Eksp. Teor. Fiz. 107, 393 (1995) [Sov. Phys. JETP 80, 208 (1995)].
  * [11] U. Jentschura and K. Pachucki, Phys. Rev. A 54, 1853 (1996).
  * [12] G. S. Adkins, B. Akers, M. F. Alam, L. M. Tram, X. Zhang, ”Calculation of higher order corrections to positronium energy levels, Proc. Sci. 353 (2019).
  * [13] J. Zatorski, Phys. Rev. A 78, 032103 (2008).
  * [14] G. S. Adkins, D. B. Cassidy, and J. Pérez-Ríos. ”Precision spectroscopy of positronium: Testing bound-state QED theory and the search for physics beyond the Standard Model.” Physics Reports 975 (2022): 1-61.
  * [15] U.D. Jentschura, A. Czarnecki, and K. Pachucki, Phys. Rev. A 72, 062102 (2005).
  * [16] H.A. Bethe and E.E. Salpeter, Quantum Mechanics Of One- And Two-Electron Atoms (Plenum Publishing Corporation, New York, 1977).