Three-photon exchange nuclear structure correction in hydrogenic systems
Krzysztof Pachucki Vojtěch Patkóš Vladimir A. Yerokhin

Version 1.0




[]krp@fuw.edu.pl
[]www.fuw.edu.pl/ krp

Faculty of Physics, Warsaw University,
             Pasteura 5, 02-093 Warsaw, Poland


Faculty of Mathematics and Physics, Charles University,  Ke Karlovu 3, 121 16 Prague
2, Czech Republic


Center for Advanced Studies, Peter the Great St. Petersburg Polytechnic University,
Polytekhnicheskaya 29, 195251 St. Petersburg, Russia





The complete relativistic O(α^2) nuclear structure correction to the energy levels of
ordinary (electronic) and muonic hydrogen-like atoms is investigated. The elastic part of the
nuclear structure correction is derived analytically. The resulting formula is valid for an
arbitrary hydrogenic system and is much simpler than analogous expressions previously reported in
the literature. The analytical result is verified by high-precision numerical calculations. The
inelastic O(α^2) nuclear structure correction is derived for the electronic and muonic
deuterium atoms. The correction comes from a three-photon exchange between the nucleus and the
bound lepton and has not been considered in the literature so far. We demonstrate that in the
case of deuterium, the inelastic three-photon exchange contribution is of a similar size and of
the opposite sign to the corresponding elastic part and, moreover, cancels exactly the model
dependence of the elastic part.  The obtained results affect the determination of nuclear charge
radii from the Lamb shift in ordinary and muonic atoms.


31.30.jr, 36.10.Ee, 14.20.Dh

Three-photon exchange nuclear structure correction 
 in hydrogenic systems
    Vladimir A. Yerokhin
    December 30, 2023
==========================================================================



§ INTRODUCTION


The determination of the nuclear charge radii from atomic spectra is a very interesting test of the
Standard Model of fundamental interactions. The lepton universality, namely the identical
interaction strength of all leptons, ensures that the nuclear charge radii derived from the
ordinary (electronic) and the muonic atoms should be exactly the same. However, a series of
experiments on μH <cit.> and μD <cit.> and (still unpublished) measurements
on μ^3He and μ^4He <cit.> revealed significant discrepancies for the determined
nuclear charge radii, as compared to those derived from the corresponding electronic atoms. In
order to verify these discrepancies one should carefully examine all possible sources of
uncertainties in the spectroscopic determinations of the nuclear charge radii.

The main theoretical uncertainty of the Lamb shift in light muonic atoms comes from our
insufficient knowledge of the nuclear internal structure. The nuclear structure corrections are
usually divided into the elastic and the inelastic part. The elastic part (also referred to as the
finite nuclear size correction) is induced by a static distribution of the nuclear charge  and can
be obtained by solving the Dirac equation. The inelastic nuclear correction  is much more
complicated; it encompass the nuclear dipole polarizability and higher-order contributions. To deal
with the nuclear corrections, one performs an expansion of the binding energy in powers of the fine
structure constant α and examines the expansion terms one after another.

The leading nuclear effect is of order α^4 and of pure elastic origin. The first-order
O(α) nuclear-structure correction (often referred to as the two-photon exchange
contribution) has both elastic and inelastic parts and was extensively studied both for the
electronic and the muonic atoms <cit.>. One of the interesting results was a significant cancellation between the elastic
and the inelastic O(α) nuclear contributions.

The next-order O(α^2) nuclear structure correction comes from the three-photon exchange
between the bound lepton and the nucleus. Only the elastic part of this correction has been
addressed in the literature so far <cit.>. In the present work we demonstrate that the
inelastic O(α^2) contribution is significant and partially cancels its elastic counterpart.
We also derive formulas for the complete O(α^2) nuclear correction in deuterium. Our
calculation is performed in the nonrecoil limit and neglects the magnetic dipole and electric
quadrupole moments of the nucleus. The results obtained affect determinations of nuclear charge
radii from the precision spectroscopy of ordinary and muonic atoms. However, they are not able to
explain the previously reported discrepancy between the H-D and μH-μD isotope shift
<cit.>.


We now introduce notations for the nuclear radii that will be extensively used throughout this
paper. r_C denotes the root mean square (rms) charge radius of an arbitrary nucleus, r_C ≡√(⟨ r^2⟩). We will use specific notations for several important nuclei: r_C(
H) ≡ r_p for the rms radius of the proton, r_C( D) ≡ r_d for the rms radius of
the deuteron, and r_s for the deuteron structure radius r^2_s = r^2_d - r^2_p. Since we neglect
the finite nuclear mass effects, there is no  3/(4 m_p^2) term in r_d^2. We define r_CC =
√(⟨ r^4 ⟩) for an arbitrary nucleus, with the specific
cases of r_CC( H) ≡ r_pp for the proton, r_CC( D) ≡ r_dd for the
deuteron, and r_ss for the corresponding structure radius of the deuteron. r_Z is the third
Zemach moment defined below by Eq. (<ref>). We will also introduce two new effective nuclear
radii of arbitrary nuclei, r_C1 and r_C2, defined by Eqs. (<ref>) and (<ref>),
respectively. The corresponding specific notations are r_C1( H) ≡ r_p1 and
r_C2( H) ≡ r_p2 for the proton and r_d1 and r_d2 for the deuteron,
respectively.





§ LEADING FINITE NUCLEAR SIZE CORRECTION

In this section we rederive well-known results for the leading nuclear correction of order
α^4, which is of pure elastic (finite nuclear size) origin and induced by the one-photon
exchange between the bound lepton and the nucleus. This derivation sets the ground for our further
evaluation of higher-order corrections.

Let us assume that the nucleus is a scalar particle with the charge density ρ(q^2) in the
momentum space. The electron-nucleus interaction potential in momentum space is then

    V(q^2) = -ρ(q^2) 4 π Z α/q^2 .

The expansion coefficients of ρ in q^2,

    ρ(q^2) = 1+ρ'(0) q^2 + 1/2 ρ”(0)  q^4 + … ,

can be interpreted in terms of momenta of the nuclear charge distribution ⟨ r^2⟩ and
⟨ r^4⟩,

    ρ'(0)    =    -⟨ r^2⟩/6, 
    ρ”(0)    =   ⟨ r^4⟩/60 .

From the second term in the right-hand side of Eq. (<ref>), one immediately obtains the leading
finite nuclear size correction to the potential,

    δ V = -ρ'(0) 4 π Z α δ^(3)(r) ,

and to the energy level of a hydrogenic system,

    E^(4)_ fns = ⟨ϕ|δ V|ϕ⟩ =
      2 π/3 Z α ⟨ r^2⟩ ϕ^2(0) ,

where, for nS states,

    ϕ^2_nS(0) = < δ^(3)(r)> = (μ Zα)^3/π n^3 ,

and μ = m/(1+m/M) is the reduced mass of the atom.

In order to establish the importance of higher-order effects, we will need numerical values of the
leading finite nuclear size effect in hydrogen and deuterium. The corresponding results, obtained
assuming r_p = 0.840 87 fm and r_d = 2.125 62 fm, are, for the electronic atoms,

    E^(4)_ fns(2S-1S,H)        =    -1 368 396 r_p^2 =   -9 67 541  Hz , 
    
      E^(4)_ fns(2S-1S,D)        =    -1 369 513 r_d^2 = -6 187 818  Hz ,

and for the muonic atoms,

    E^(4)_ fns(2P-2S,μ H)    =    -5.197 45  r_p^2 =   -3.674 92  meV , 
    
      E^(4)_ fns(2P-2S,μ D)    =    -6.073 18  r_d^2 = -27.440 22  meV .


We also observe that one of the relativistic O(α^2) corrections comes from the third term in
Eq. (<ref>),

    δ^(2) V = 1/2 ρ”(0) 4 π Z α ∇^2δ^(3)(r) .

Since its expectation value on nS states is singular, we will use dimensional regularization and
combine this part with other O(α^2) corrections to obtain a finite result.




§ TWO-PHOTON EXCHANGE NUCLEAR STRUCTURE: MUONIC ATOMS

 In this section we address the leading O(α) nuclear structure contribution
E^(5) in muonic atoms, which originates from the two-photon exchange between the bound lepton
and the nucleus.

The elastic part E^(5)_ fns can be obtained from the forward two-photon scattering
amplitude at zero momentum

    E^(5)_ fns = ϕ^2(0) ∫d^3q/(2 π)^3  Tr[γ^0 V 1/p-m γ^0 V (I+γ^0)/4] ,

with V = -4 π Z α/q^2 and p=(m, q⃗). This leads to the so-called Friar correction
<cit.>,

    E^(5)_ fns   =   
    -(4 π Z α)^2 ϕ^2(0) 2 m∫d^3q/(2 π)^3 ρ^2(q^2)-1-2 q^2 ρ'(0)/q^6
       =     -π/3 ϕ^2(0) (Z α)^2 m r_Z^3 ,

where

    r_Z^3 =  ∫ d^3r_1∫ d^3r_2 ρ(r_1) ρ(r_2) |r⃗_1-r⃗_2|^3 .

As pointed out in Refs. <cit.>, it is important to consider the Friar
correction E^(5)_ fns together with the corresponding inelastic part, because of a
cancellation between them, occurring both for the muonic and the ordinary atoms. For this reason,
we do not separate out E^(5)_ fns but absorb it in the total nuclear structure correction
E^(5).



 §.§ Muonic hydrogen


The inelastic two-photon exchange correction in μH has been extensively studied in the
literature (see Ref. <cit.> and references therein). It is also given by the forward
scattering amplitude and can be parameterized in terms of two spin-independent structure functions
of the proton. Using dispersion relations, these functions are usually expressed in terms of the
cross section of the inelastic photon scattering off the proton, which is extracted from
experiment. The main problem of this approach is that one of the dispersion relations involves
subtractions that can only be obtained from theory, and this introduces the dominant uncertainty.

There is good agreement between different calculations of the two-photon exchange correction, with
the final result of E^(5) (2P_1/22S, μ H) = E^(5)_ fns +
E^(5)_ pol = 0.033 2(20) meV assumed by the CREMA collaboration <cit.> in
their determination of the proton charge radius. It is convenient to parameterize this result in
terms of an effective radius r_pF, in analogy to Eq. (<ref>),

    E^(5)(μ H) =  -π/3 ϕ^2(0) (Z α)^2 m r_pF^3 ,

with

    r_pF^3 = 3.270 (197)   fm^3 .

This parametrization will be used below in our calculation of the inelastic contribution in other
muonic atoms, see Eq. (<ref>).



 §.§ Muonic atoms other than hydrogen


For all nuclei other than the proton, the inelastic contribution is dominated by the electric
dipole polarizability. For muonic atoms, one may assume the nonrelativistic approximation, so the
second-order correction due to the electric dipole nuclear excitation is

    E^(5)_ pol0 = α^2 ⟨ϕ ϕ_N|
    d⃗·r⃗/r^3 1/E_N + E_0 - H_N -H_0 d⃗·r⃗/r^3 |ϕ ϕ_N⟩ ,

where d⃗ is the electric dipole operator divided by the elementary charge, and H_0 and
H_N are the nonrelativistic Coulomb Hamiltonian for the muon and the nucleus, respectively. To
the leading order in α, one may neglect the Coulomb interaction and replace
ϕ(r)→ϕ(0) to obtain a compact formula for the leading two-photon exchange
contribution,

    E^(5)_ pol0 = -4 π α^2/3 ϕ^2(0)
      ⟨ϕ_N|d⃗ √(2 m/H_N-E_N)d⃗ |ϕ_N⟩ ,

which contributes 1.910 meV to the 2P-2S transition energy in muonic deuterium
<cit.>.

There are many corrections to the leading contribution <cit.>, the most interesting of them being the one that partially cancels the
Friar correction. To show this, following Ref. <cit.>, we consider the muonic
matrix element P for the nonrelativistic two-photon exchange

    P =∑_i,j⟨ϕ|α/|r⃗-R⃗_i|1/(H_0-E_0+E)α/|r⃗-R⃗_j'||ϕ⟩,

where H_0 is the nonrelativistic Hamiltonian for the muon (electron) in the nonrecoil limit, and
R⃗_i is a position of the ith proton with respect to the nuclear mass center. Using the
on-mass-shell approximation, subtracting  the leading Coulomb interaction, the finite nuclear size,
and the electric dipole polarizability, and expanding in the small parameter √(2 m E) |R⃗_i-R⃗_j'|, we obtain

    P    =   α^2 ϕ^2(0) ∫d^3 q/(2 π)^3 (4 π/q^2)^2 (E+q^2/2 m)^-1 [e^i q⃗·(R⃗-R⃗')-1+q^2/6 (R⃗-R⃗')^2]
    
       ≈   ∑_i,jπ/3 m α^2 ϕ^2(0) |R⃗_i-R⃗_j'|^3
    (1-1/5 √(2 m E)|R⃗_i-R⃗_j'| +… ) .

The corresponding correction to the atomic energy is

    E^(5)_ pol1 = -∑∫ dE∫d^3 R d^3R' ϕ^*_N(R⃗) ϕ_E(R⃗) ϕ^*_E(R⃗') ϕ_N(R⃗') P

Let us consider only the first, E-independent term. When ϕ_E = ϕ_N, it corresponds to the
elastic part, namely, the Friar correction given by Eq. (<ref>). However, the inclusion of all
excited states leads to

    E^(5)_ pol1 = -π/3 m α^2 ϕ^2(0) ∑_i,j=1^Z⟨ϕ_N||R⃗_i-R⃗_j|^3|ϕ_N⟩ ,

which is much different from Eq. (<ref>), in particular it vanishes for deuterium.

There are further nuclear polarizability corrections which were extensively studied in the
literature <cit.>. It is convenient to
write the final result for the two-photon exchange nuclear structure correction separating out the
contribution due to the two-photon exchange with individual nucleons,

    E^(5) = E^(5)_ pol - π/3 m α^2 ϕ^2(0) [Z r^3_pF + (A-Z) r^3_nF] ,

where E^(5)_ pol = E^(5)_ pol0 + E^(5)_ pol1 + ….  Such representation
of the nuclear structure correction is particularly advantageous for calculating isotope shifts,
since the individual nucleon contributions partially cancel each other in the difference, together
with the corresponding uncertainties. Calculating the two-photon exchange with individual nucleons,
we take the effective proton radius r^3_pF from Eq. (<ref>), whereas for the neutron we assume
the corresponding parameter to be four times smaller than that of the proton, r_nF^3 =
r_pF^3/4, with uncertainty of 100%. This choice of r_nF^3 is in agreement with results
summarized in Ref. <cit.> but requires further investigations.

Despite the fact that the literature results for E^(5)_ pol in μD reported by
different groups <cit.> are in good
agreement with each other (see a summary in Ref. <cit.>), one should bear in mind that
a number of higher-order effects exist that have not yet been addressed in any of the previous
studies. Specifically, it has not so far been possible to include nucleon relativistic corrections
to the coupling of the nucleus to the electromagnetic field. We thus believe that all theoretical
predictions of E^(5)_ pol in μD should bear an uncertainty whose relative value is 
approximately the ratio of the average nucleon binding energy to the nucleon mass, which is about 1%.

Summarizing our analysis of the existing literature results, we adopt the sum of entries p_1…
p_12 labelled as “Our choice” in Table 3 in Ref. <cit.> as currently the best value of the
two-photon nuclear polarizability correction to the 2P_1/2 2S transition energy
in μD, and ascribe the uncertainty of 1% to it,

    E^(5)_ pol(2P_1/2 2S,μ D) = 1.6625 (166)  meV .

The above uncertainty of E^(5)_ pol is about 50% larger than the corresponding
estimate of ± 0.0107 meV given in Table 3 of Ref. <cit.>. Finally, we add the
individual nucleon part in Eq. (<ref>) and obtain the total two-photon nuclear structure
correction to the 2P_1/2 2S transition energy in μD,

    E^(5)(2P_1/2 2S,μ D) = 1.7110 (194)  meV ,

which almost coincides with the corresponding result  of 1.7096 (200) meV from
Ref. <cit.>, as given by Eq. (17) of that work.




§ TWO-PHOTON EXCHANGE NUCLEAR STRUCTURE: ELECTRONIC ATOMS


The elastic (finite nuclear size) part of the two-photon exchange nuclear structure correction
for electronic atoms is given by the same formula as for the muonic atoms, Eq. (<ref>).




 §.§ Electronic hydrogen


We calculate the elastic part of the nuclear structure correction for hydrogen according to
Eq. (<ref>) and using the result for the third Zemach moment from Ref. <cit.> obtained
by averaging values measured in scattering experiments,

    r_pZ = 1.587 (26) r_p .

The corresponding result for the 2S–1S transition is

    E^(5)_ fns(2S  1S,e H) = 0.0307 (15)  .


The inelastic part of the two-photon exchange nuclear structure correction E^(5)_ pol was
derived in the logarithmic approximation in Ref. <cit.>,

    E^(5)_ pol(2S  1S,e H) = -m α ϕ^2(0) [5 α_p -β_p] lnE̅_p/m ,

where E̅_p is the average proton excitation energy and α_p and β_p are the
static proton polarizabilities extracted from experiment. Using the same average proton excitation
energy E̅_p = 410 MeV as in Ref. <cit.> and the updated results for the
proton polarizabilities <cit.>,

    α_p    =    10.65 (35)(20)(30) ^3 ,
    β_p    =    3.15 (35)(20)(30) ^3 ,

we obtain the result for the 2S–1S transition of

    E^(5)_ pol(2S  1S,e H) = 0.0567 (85)  ,

where, following Ref. <cit.>, we assumed a 15% uncertainty due to the leading
logarithmic approximation.

The total result for the two-photon nuclear structure correction in electronic hydrogen is

    E^(5)(2S  1S,e H)  = E^(5)_ fns + E^(5)_ pol = 0.0874 (86)  ,

which could be compared with the corresponding result of 0.091 (11) kHz from
Ref. <cit.>.




 §.§ Electronic atoms other than hydrogen


Similarly to the muonic atoms, it is convenient to write the total two-photon exchange nuclear
structure correction separating out the contribution due to the interaction with individual
nucleons,

    E^(5) = E^(5)_ pol - m α ϕ^2(0) [π/3 α Z r^3_pZ
      +Z [5 α_p -β_p] lnE̅_p/m
       +(A-Z) [5 α_n -β_n] lnE̅_n/m] .

In the above formula, the first and the second terms in the brackets represent the elastic and the
inelastic interactions with individual protons, respectively, whereas the third term comes from the
inelastic interaction with individual neutrons. The parameters for the protons are the same as for
hydrogen, whereas for the neutrons we use the experimental polarizabilities <cit.>,

    α_n    =    11.55 (125)(20)(80) ^3 ,
    β_n    =    3.65 (125)(20)(80) ^3 ,

and the same value of E̅_n = 410 MeV as for the proton. We note that the elastic interaction
of the bound electron with the nucleus as a whole is absorbed in E^(5)_ pol, reflecting
the fact that the third Zemach moment correction for a compound nucleus largely cancels out between
the elastic and inelastic parts in the same way as in muonic atoms.

Similarly to the muonic atoms, the nuclear polarizability correction E^(5)_ pol in
Eq. (<ref>) comes from the electric dipole polarizability, which, however, takes a very
different form for the electronic atoms. Since in this case the nonrelativistic approximation is
not valid, one should consider the complete two-photon exchange and keep the relativistic form of
the matrix elements,

    E_ pol^(5)   =   
    i e^2 ϕ^2(0) ∫d ω/2 π ∫d^3k/(2 π)^3 ω^2 (δ^ik-k^i k^k/ω^2)/ω^2-k^2 (δ^jl-k^j k^l/ω^2)/ω^2-k^2
          × Tr[(γ^j 1/p-k-mγ^i +
    γ^i1/p+k-m γ^j) (γ^0+I)/4] 
          ×⟨ϕ_N|d^ k 1/E_N-H_N-ω 
    d^ l|ϕ_N⟩ + …
       =    E^(5)_ pol1 + E^(5)_ pol2 +  E^(5)_ pol3 +… ,

where p=(m,0⃗). Assuming that the electron mass is much smaller than the nuclear excitation
energy, the leading nuclear polarizability correction becomes

    E^(5)_ pol1 = -m α^2 ϕ^2(0) 2/3 ⟨ϕ_N|d⃗ 1/H_N-E_N [19/6 + 5 ln2 (H_N-E_N)/m] d⃗ |ϕ_N⟩ .

The corresponding contribution to the 2S–1S transition in ordinary deuterium is
19.26 (6) kHz <cit.>.

Various small corrections to the electric dipole polarizability for electronic atoms were
considered by Friar in <cit.>. In particular, it was shown there that the Zemach
contribution for deuterium vanishes in the same way as for the muonic deuterium. Furthermore, the
higher-order terms in the m/(H_N-E_N) expansion of Eq. (<ref>) give rise to a correction

    E^(5)_ pol2 =
        -m^3 α^2 ϕ^2(0) 2/3 ⟨ϕ_N|d⃗ 1/(H_N-E_N)^3 [-283/80 + 15/4 ln2 (H_N-E_N)/m] d⃗ |ϕ_N⟩ ,

which contributes 0.106 kHz to the 2S–1S transition in ordinary deuterium <cit.>.
Another important correction is the one due to the magnetic suceptibility <cit.>,

    E^(5)_ pol3 =
        m α^2 ϕ^2(0) 2/3 ⟨ϕ_N|μ⃗ 1/(H_N-E_N)' [-1/6 + ln2 (H_N-E_N)/m] μ⃗ |ϕ_N⟩ ,

where μ is the magnetic moment operator divided by an elementary charge.
It leads to a correction of -0.307(2)(6) kHz to the
2S–1S transition in eD <cit.>.

There were further corrections to the electric dipole polarizability considered in
Ref. <cit.>. However, we are convinced that they were not treated correctly and,
moreover, that there are many more relativistic corrections of the same order. For this reason we
disregard the additional corrections from Ref. <cit.> and assume the total
polarizability correction to be the sum of Eqs. (<ref>), (<ref>), and (<ref>).
Specifically, the result for the nuclear polarizability to the 2S–1S transition in electronic
deuterium is

    E_ pol^(5)(2S 1S,D) = 19.06 (20)  kHz .

Adding the individual nucleon part contribution of 0.149 (22) kHz, we obtain the total
two-photon exchange nuclear structure contribution of

    E^(5)(2S 1S,D) = 19.21 (20)  kHz ,

which could be compared with the sum of the nuclear polarizability correction and the third Zemach
contribution from Ref. <cit.>, 18.70 (7) + 0.51 = 19.21 (7) kHz, perfect
agreement of the numerical values being probably accidental.




§ THREE-PHOTON EXCHANGE ELASTIC CONTRIBUTION

This contribution has been studied by different methods and a number of authors, of note
analytically by Friar in Ref. <cit.>, and numerically by solving the Dirac equation in
the field of finite size nucleus <cit.>. Here we present an alternative analytical
approach, which leads to much simpler analytic formulas. A numerical verification of our formulas
is given in Appendix <ref>.

In the standard analytic approach, one applies the perturbation theory to the Dirac energies with
the perturbing potential δ V = V - V_0 where V_0 = -Z α/r and V is the Coulomb
potential from the finite size nucleus,

    δ^(1)E    =   ⟨ψ̅|δ V|ψ⟩
    δ^(2)E    =   ⟨ψ̅|δ V 1/(p-γ^0 V_0-m)' δ V|ψ⟩
    δ^(3)E    =   ⟨ψ̅|δ V 1/(p-γ^0 V_0-m)' δ V 1/(p-γ^0 V_0-m)' δ V|ψ⟩

One can use the exact Dirac wave function and the reduced Dirac propagators to calculate the
O(α^2) correction to the finite nuclear size <cit.>, which we call the elastic
three-photon exchange correction. However, we will not use the above formulas but employ a
different approach, which we call the scattering amplitude approach. In this approach, the
O(α^2) relativistic correction to the finite nuclear size is induced by the elastic three
photon exchange. The corresponding correction can be divided into the low and the high energy
momentum exchange parts, E^(6)_ fns = E_L + E_H. These parts are calculated as follows.



 §.§ Three-photon exchange: low energy part

The low-energy part E_L is again split into two parts

    E_L    =    E_L1 + E_L2 ,
    
       E_L1   =   ⟨δ V 1/(E-H)' δ V⟩  + ⟨δ^(2) V⟩ , 
    
       E_L2   =   ⟨ϕ|1/8 m^2 ∇^2(δ V) +
       1/4 m^2 σ⃗·∇⃗(δ V)×p⃗|ϕ⟩
          
      + 2 ⟨ϕ|δ V1/(E-H)' [-p^4/8 m^3 +π Z α/2 m^2 δ^3(r)
       ]|ϕ⟩ ,

where E_L1 is the nonrelativistic contribution proportional to r_C^4, and E_L2 is the
relativistic part proportional to r_C^2. All these matrix elements are calculated in
d=3-2 ϵ dimensions. The following results are obtained for the nS states,

    E_L1(nS)    =    [ρ'(0)]^2 16 (Z α)^6/n^3[-1/n -1/2 + γ -lnn/2+Ψ(n)]
           +(Z α)^3 [-1/ϵ + 4 ln(Z α)] 4 [ρ'(0)]^2 ⟨π δ^(d)(r)⟩+
      4 (Z α)^6 ρ”(0)/n^5 ,
    
      E_L2(nS)    =   ρ'(0) 4 (Z α)^6/n^3[9/4n^2-1/n-5/2 +γ -lnn/2+Ψ(n)]
           +(Z α)^3 [-1/ϵ + 4 ln(Z α)] ρ'(0) ⟨π δ^(d)(r)⟩ ,

and for the nP states,

    E_L(nP_1/2)    =   (Z α)^6 ( -9/4 ρ'(0) + 3 ρ”(0)) R'_n1(0)^2  ,
    
      E_L(nP_3/2)    =   (Z α)^6 3 ρ”(0) R'_n1(0)^2  ,

where

    R'_n1(0)^2 = 4/9 n^3 (1-1/n^2) .

For all higher-L states E_L vanishes.



 §.§ Three-photon exchange: high energy part

We start by introducing the two potentials in d-dimensions that will appear in the evaluation of
the high energy-part E_H,

    V_d(r)    =    4 π ∫d^dq/(2 π)^d e^iq⃗·r⃗ ρ(q^2)/q^2 , 
    
      V_d^(2)(r)    =    4 π ∫d^dq/(2 π)^d e^iq⃗·r⃗ ρ(q^2)/q^4 .

Their large r asymptotics are

    V_d(r)     =     V(r) +  , 
    
      V_d^(2)(r)     =     V^(2)(r) + ρ'(0)  V(r) +  ,

and in d=3,

    V(r)     =    1/r +  , 
    
      V^(2)(r)     =     -r/2 + ρ'(0)/r +  ,

where the local terms vanish outside the nucleus.

Now we proceed to the derivation of the high-energy part E_H. It is given by the three-photon
scattering amplitude with momenta p_i = (m,q⃗_i),

    E_H    =    -(4 π Z α)^3 ϕ^2(0) ∫d^dq_1/(2 π)^d ∫d^dq_2/(2 π)^d ρ(q_1^2)/q_1^4 ρ(q_2^2)/q_2^4 ρ((q⃗_1-q⃗_2)^2)/(q⃗_1-q⃗_2)^2
          × Tr[(p_1+m) γ_0  (p_2+m) (γ_0+I)/4] .

The above trace equates to 4 m^2 +q⃗_1 q⃗_2, so we can split E_H into the
nonrelativistic and relativistic parts,

    E_H  = E_H1 + E_H2 .

The nonrelativistic part E_H1 is

    E_H1   =    -(4 π Z α)^3 ϕ^2(0) 4 m^2∫d^dq_1/(2 π)^d∫d^dq_2/(2 π)^d ρ(q_1^2)/q_1^4 ρ(q_2^2)/q_2^4 ρ((q⃗_1-q⃗_2)^2)/(q⃗_1-q⃗_2)^2
       =   
     -ϕ^2(0) (Z α)^3 4 m^2∫ d^dr   V_d(r) [V_d^(2)(r)]^2  .

In order to calculate this integral, we split the integration region into r<Λ and r≥Λ. The first integral is finite in d=3 but diverges at large Λ, and in the second
integral one can use the asymptotic form of potentials,

    E_H1   =   
     -ϕ^2(0) (Z α)^3 4 m^2{4 π ∫^Λ dr r^2 V(r) [V^(2)(r)]^2
      +∫_Λ d^d r  V(r) [ V^(2)(r)  + ρ'(0)  V(r)]^2}
       =   ϕ^2(0) (Z α)^3 4 m^2 4 π ∫^∞ dr ln(r/r_C) d/dr r^3{V(r) [V^(2)(r)]^2 -1/r (r/2-ρ'(0)/r)^2
        }
          
      -ϕ^2(0) (Z α)^3 4 m^2 ρ'(0)^2 {4 π ln(Λ/r_C)
      +∫_Λ d^d r  V(r)^3} .

The expression under the first integral is a local function of r, so this integral is effectively
over the nuclear size, which allows us to introduce an effective nuclear radius r_C1 as

    lnr_C1/r_C +2= 36/r_C^4 ∫^∞ dr ln(r/r_C) d/dr r^3{V(r) [V^(2)(r)]^2 -1/r (r/2-ρ'(0)/r)^2} .

So, the first O(α^2) correction E_H1 is represented in the following form

    E_H1 = 4 π ϕ^2(0) (Z α)^3 4 m^2 r_C^4/36 [1/4 ϵ + 5/2+γ +ln(r_C1 m)] .

The relativistic part E_H2 is
  
    E_H2   =    -(4 π Z α)^3 ϕ^2(0)∫d^dq_1/(2 π)^d∫d^dq_2/(2 π)^d ρ(q_1^2)/q_1^4 ρ(q_2^2)/q_2^4 ρ((q⃗_1-q⃗_2)^2)/(q⃗_1-q⃗_2)^2 q⃗_1q⃗_2
      
       =   ϕ^2(0) (Z α)^3 ∫ d^dr  {2 π ρ(r) V_d^(2)(r) - [V_d(r)]^2} V_d^(2)(r)  ,

  and we proceed in a similar way as in the case of E_H1, namely

    E_H2   =   ϕ^2(0) (Z α)^3 {4 π ∫^Λ dr r^2
      [2 π ρ(r) V^(2)(r) - [V(r)]^2] V^(2)(r)
          
      -∫_Λ d^d r [ V(r)]^2 ( V^(2)(r)  + ρ'(0)  V(r))}
       =   
      -ϕ^2(0) (Z α)^3 4 π ∫^∞ dr ln(r/r_C)
      d/dr r^3[2 π ρ(r) [V^(2)(r)]^2 - [V(r)]^2 V^(2)(r)
        
           - 1/r^2 (r/2-ρ'(0)/r)]
      -ϕ^2(0) (Z α)^3 ρ'(0) [ 4 π ln(Λ/r_C)
      + ∫_Λ d^d r [ V(r)]^3]  .

The expression under the first integral is a local function, so we can introduce the second
effective nuclear radius r_C2 as

    lnr_C2/r_C-1 = 6/r_C^2 ∫^∞ dr ln(r/r_C)
      d/dr r^3[2 π ρ(r) [V^(2)(r)]^2 - [V(r)]^2 V^(2)(r)
        - 1/r^2 (r/2-ρ'(0)/r)] .

So, the second O(α^2) correction is given by

    E_H2 = -4 π (Z α)^3 ϕ^2(0) r_C^2/6 [1/4 ϵ -1/2+γ + ln(r_C2 m)]  .




 §.§ Three-photon elastic exchange: total result

The complete O(α^2) finite nuclear size correction for an arbitrary nucleus is given by the
sum E^(6)_ fns = E_L + E_H = E_L1 + E_L2 + E_H1 + E_H2, with the result

    E^(6)_ fns(nS)    =   -(Z α)^6 m^3 r_C^2 2/3 n^3 [9/4n^2 - 3 -1/n +2 γ -lnn/2+Ψ(n)
            +ln(m r_C2 Z α)]
            
           + (Z α)^6 m^5 r_C^4 4/9 n^3 [-1/n + 2 + 2 γ -lnn/2+Ψ(n)
            + ln(m r_C1 Z α)]
            
          
            +(Z α)^6 m^5 r_CC^4 1/15 n^5 ,
    
      E^(6)_ fns(nP_1/2)    =   (Z α)^6 m (m^2 r_C^2/6 + m^4 r_CC^4/45) 1/n^3 (1-1/n^2)  , 
    
      E^(6)_ fns(nP_3/2)    =   (Z α)^6 m^5 r_CC^4 1/45 n^3 (1-1/n^2)  ,
    
      E^(6)_ fns(nL_J)    =    0   L>1 ,

where r_CC^4 = ⟨ r^4⟩ and the effective nuclear charge radii r_C1 and r_C2
defined by Eqs. (<ref>) and (<ref>) encode the high-momentum contributions and are expected
to be of the order of r_C. Equations (<ref>)-(<ref>) are valid both for electronic and
muonic atoms. However, in the case of the electronic atoms, the terms proportional to r_C^4 and
r_CC^4 in these formulas are smaller than the next-order correction and thus should be
neglected.

Equations (<ref>)-(<ref>) depend on the nuclear model through the effective nuclear charge
radii r_C1 and r_C2. However, as we demonstrate below, the terms with r_C1 and r_C2 exactly cancel in the sum with the corresponding inelastic contribution, so their model
dependence is irrelevant for the finite result.

Table <ref> presents our results for the effective charge radii r_C1 and r_C2,
for two models of the nuclear charge distribution. The ratios listed in the table are proven to be
independent of the nuclear model parameters λ and a, thus making the corresponding
results valid for an arbitrary nucleus and both for the electronic and muonic atoms.

The formulas for E^(6)_ fns have been derived in the nonrecoil limit, i.e., assuming the
infinite nuclear mass. This is different from the approach by Friar in Ref. <cit.>, in
which he replaced the lepton mass by the reduced mass of the system. We do not think such a
replacement is valid. However, recoil corrections for muonic atoms are significant and can be
partially accounted for by the (μ/m)^3 scaling factor that comes from the square of the
nonrelativistic wave function at origin.

Equations (<ref>)-(<ref>) can be compared with the analytical results by Friar derived within
a different approach <cit.>. However, the formulas of Ref. <cit.> are so
complicated that a direct comparison is not possible, except for the state dependence which is in
perfect agreement. A comparison of numerical results presented in Sec. <ref> shows a
reasonable but not perfect agreement. In order to verify our formulas, we performed a
high-precision numerical calculation, by solving the Dirac equation numerically and identifying the
O(α^2) finite nuclear size contribution, as described in Appendix <ref>. Perfect agreement
between analytical and numerical approaches confirms the correctness of Eqs. (<ref>)-(<ref>).






§ INELASTIC THREE-PHOTON EXCHANGE CORRECTION IN MUONIC DEUTERIUM

The inelastic three-photon exchange nuclear structure correction has not yet been studied in the
literature and is the main topic of this work. With momenta of order q∼√(2 m Λ),
which is much lower than the inverse of the nuclear size, the muon kinetic energy becomes
comparable to the characteristic nuclear excitation energy Λ, and thus the muon starts to
probe the nuclear structure and see individual nucleons. It means that the total correction does
not involve contributions coming from muon momenta of the order of the inverse of nuclear size, and
there is no place for the elastic high energy parts encoded in r_C1 and r_C2 effective
nuclear radii.

We represent the total nuclear structure correction E^(6) as a sum of several parts,

    E^(6) = E^(6)_1 + E^(6)_2 + E_C + E^(6)_ np =
      E^(6)_ fns + E^(6)_1,  pol + E^(6)_2,  pol + E_C + E^(6)_ np ,

where the elastic part E^(6)_ fns = E^(6)_1,  fns + E^(6)_2,  fns was
calculated in the previous section, E_1^(6) and E_1,  pol^(6) are proportional to
r_C^4, E_2^(6) and E_2,  pol^(6) are proportional to r_C^2, E_C is the known
Coulomb distortion correction <cit.>, and E^(6)_ np is
the contribution due to the inelastic interaction with individual nucleons.

Below we calculate E^(6)_1 and E^(6)_2. Because the two-photon exchange nuclear structure
correction  E^(5)_ pol was previously shown to be dominated by the electric dipole types
of the nuclear polarizability, we will assume that the same holds for the three-photon exchange
nuclear structure correction.



 §.§ Inelastic contribution ∝ R^2

We represent the total E_2^(6) correction as a sum of several parts

    E_2^(6) = E_L2 + E_C + E_R + E_H2(p)  ,

calculated in the following.
Let us consider the two radiative photon exchange between the muon and the nucleus, taking into
account the Coulomb interaction V. The corresponding energy shift is

    δ E    =   
    i e^2 ∫d ω/2 π ∫d^3k_1/(2 π)^3 ∫d^3k_2/(2 π)^3 ω^2 (δ^ik-k_1^i k_1^k/ω^2)/ω^2-k_1^2 (δ^jl-k_2^j k_2^l/ω^2)/ω^2-k_2^2
          ×[
    ⟨ψ̅|γ^j e^i k⃗_2r⃗ 1/p-γ^0 V - γ^o ω-mγ^i e^i k⃗_1r⃗|ψ⟩
          
    + ⟨ψ̅|γ^i e^i k⃗_1r⃗ 1/p-γ^o V+γ^0 ω-m γ^j e^i k⃗_2r⃗ |ψ⟩]
    
          ×⟨ϕ_N|d^ k 1/E_N-H_N-ω 
    d^ l|ϕ_N⟩,

where d⃗ is the dipole moment operator divided by the elementary charge. In the case of
deuteron d⃗ is equal to the position of the proton with respect to the mass center d⃗ =
R⃗. In the nonrelativistic limit, δ E takes the well known form of Eq. (<ref>).
The corresponding low-energy α^6 contribution is obtained from Eq. (<ref>) by assuming
that the muon momenta are of the order of m α. Then E_0-H_0  can be neglected in
comparison to E_N-H_N and one obtains (with d=3-2 ϵ),

    δ_L E = α^2 ⟨ϕ|1/r^4|ϕ⟩_ϵ 1/d ⟨ϕ_N| R⃗  1/E_N-H_N R⃗ |ϕ_N⟩ ,

and

    ⟨ϕ|1/r^4|ϕ⟩_ϵ
      = ⟨[∇ V(r)]^2⟩
      = ⟨1/r^4⟩ +ϕ^2(0) 4 π (-1/2 ϵ +2 lnα + 2)
       ,

were ⟨ 1/r^4 ⟩ is defined as an integral from a small radius a to infinity with the
1/a and ln a+γ terms subtracted out. The high-energy α^6 part is obtained by
assuming that muon momenta are of the order √(2 m Λ). Then we can use the
explicit Coulomb correction

    δ_H E    =   
    -i e^2 ∫d ω/2 π ∫d^3k_1/(2 π)^3 ∫d^3k_2/(2 π)^3 ω^2 (δ^ik-k_1^i k_1^k/ω^2)/ω^2-k_1^2 (δ^jl-k_2^j k_2^l/ω^2)/ω^2-k_2^2 4 π α/(k⃗_1+k⃗_2)^2
           Tr[(
    γ^j 1/p + k_2-m γ^0 1/p - k_1-m γ^i
    
    +γ^0 1/p -k_1-k_2-m γ^j 1/p - k_1-m γ^i
    
          
    +γ^j 1/p +k_2-m γ^i 1/p + k_1+k_2 -mγ^0
    )  (γ^0+I)/4] ϕ^2(0)
    
          ×[⟨ϕ_N|R^ k 1/E_N-H_N-ω 
    R^ l|ϕ_N⟩ +
    ⟨ϕ_N|R^ l 1/E_N-H_N+ω 
    R^ k|ϕ_N⟩]  ,

where k_1 = (ω, k⃗_1), k_2 = (-ω, k⃗_1), and p = (m,0⃗). Assuming that
(H_N-E_N)/m is small, one performs an expansion and obtains

    δ_H E    =   
    -π α^3 ϕ^2(0) 1/d [⟨ϕ_N|R⃗4/H_N-E_N (1/2 ϵ-1+ln2-ln(H_N-E_N)/m) R⃗ |ϕ_N⟩
          
      +⟨ϕ_N|R⃗(1/2 ϵ+13/2-13 ln2
        -5 ln(H_N-E_N)/m)R⃗ |ϕ_N⟩]  .

The sum of δ_L E in Eq. (<ref>) and δ_H E in Eq. (<ref>) gives the leading
Coulomb distortion correction E_C,

    E_C    =    -1/3 ⟨1/r^4⟩ ⟨ϕ_N|R⃗ 1/H_N-E_N R⃗ |ϕ_N⟩
          
        - 4 π/3 ϕ^2(0) ⟨ϕ_N| R⃗ 1/H_N-E_N [1+ln(2 m α^2/H_N-E_N)]R⃗ |ϕ_N⟩ ,

where

    ⟨1/r^4⟩_nS   =   8/n^3 [-5/3 + 1/2 n +1/6 n^2+γ+Ψ(n)-lnn/2] ,
    ⟨1/r^4⟩_nP   =   2 (3 n^2-2)/15 n^5 ,

and the relativistic correction E_R,

    E_R =-π α^3 ϕ^2(0) r^2_s/6 (1/ϵ+41/3-26 ln2
        -10 ln⟨ E ⟩_2/m)  ,

where r^2_s = ⟨ R^ 2⟩ is the deuteron structure radius and

    ln⟨ E⟩_2/m = 1/r_s^2 ⟨ϕ| R⃗ ln(H_N-E_N)/m R⃗ |ϕ⟩ .

Although there is no elastic high-energy part, the individual proton contributes

    E_H2(p) = -π α^3 ϕ^2(0) r_p^2/6 [1/ϵ -2+4 γ + 4 ln(r_p2 m)] .

The last contribution E_L2 is exactly the same as the one of Eq. (<ref>) with the deuteron
radii.

Finally, the total nuclear structure contribution E^(6)_2 ∝ R^2 is given by the sum of
the elastic and inelastic parts,

    E^(6)_2(nS)    =    E_C(nS)-α^6 m^3 2/3 n^3 [
      r_d^2 (9/4n^2 - 3 -1/n +2 γ -lnn/2+Ψ(n)
      +lnα)
          
    +r^2_s (47/12-13/2 ln2 -5/2 ln⟨ E ⟩_2/m-γ)
    +r_p^2 ln(r_p2 m)] .

It is remarkable that the part of the elastic contribution depending on the effective deuteron
radius r_d2 does not show up in total E^(6)_2. Separately, the expression
for the inelastic contribution is

    E^(6)_ 2,pol = E^(6)_2 - E^(6)_ 2,fns - E_C    =    -2 α^6 m^3/3 n^3 δ_l0 [
    r^2_s (47/12-13/2 ln2 -5/2 ln⟨ E ⟩_2/m) 
          
    +r_p^2 (γ + ln(r_p2 m))
    -r_d^2 (γ+ln(r_d2 m))]  .

The averaged excitation energy ⟨ E⟩_2 has been calculated using the AV18
potential <cit.> with the result

    ⟨ E⟩_2 = 7.374  MeV,

exactly the same as for μD and eD, because the lepton mass cancels out between the left and
the right side of Eq. (<ref>).



 §.§ Inelastic contribution ∝ R^4

We represent the total E_1^(6) correction as a sum of four parts

    E_1^(6) = E_L1 + E_Q + E'_Q + E_H1(p)  .

The middle energy contribution E_Q comes from momenta q∼√(2 m Λ). We derive it by considering
the nonrelativistic three Coulomb photon exchange,

    E_Q = -⟨ϕ,ϕ_N|
      α/|r⃗-R⃗| 1/E-H_0 + E_N - H_N α/|r⃗-R⃗| 1/E-H_0 + E_N - H_N α/|r⃗-R⃗||ϕ,ϕ_N⟩ ,

where H_0 is the lepton kinetic energy operator. In the corresponding electronic matrix element
P_Q, one neglects the lepton binding energy E,

    P_Q  =
    -ϕ^2(0) α^3 4 m^2 (4 π)^3∫d^dq_1/(2 π)^d d^dq_2/(2 π)^d e^i q⃗R⃗_1/q^2 1/q^2+2 m Λ e^-i (q⃗+q⃗')R⃗_2/(q⃗ + q⃗')^2 1/q'^2+2 m Λ' e^i q⃗'R⃗_3/q'^2 ,

expands in R_i, keeping terms ∝ R^4, with the result

    E_Q = -ϕ^2(0) α^3 m^2 π [
      ⟨ R^4⟩ (4/15 ln2-2/5)
      +⟨ R^2⟩^2(-10/27 + 2/3 ln2 +2/9 ln⟨ E ⟩_1/m -8/9 β - 1/9 ϵ)] ,

where

    ln⟨ E⟩_1/m   =   
      -1/⟨ R^2⟩^2 [
        ⟨ 0|R^2 ln(H-E)'/m R^2|0⟩
        -6/5 ⟨ 0|R^i ln(H-E)'/m R^2 R^j|0⟩
          
      +3/10 ⟨ 0|(R^i R^j-δ^ij R^2/3) ln(H-E)'/m (R^i R^j-δ^ij R^2/3)|0⟩]  .

The average energy ⟨ E⟩_1 does not depend on the lepton mass m, since the
dependence on m  cancels out between the left and right side of above equation. We calculate
⟨ E⟩_1 by using the AV18 deuteron potential <cit.>, with the result

    ⟨ E⟩_1 = 2.932  MeV .

Equation (<ref>) involves the dimensionless parameter β defined by

    β =
      -∑_Λ_1,Λ_2'    [   ⟨ 0|R^i R^j+3 δ^ij R^2|Λ_1⟩ ⟨Λ_1|R^i|Λ_2⟩ ⟨Λ_2|R^j|0⟩
       +   ⟨ 0|R^i|Λ_1⟩ ⟨Λ_1|δ^ij R^2-3 R^i R^j|Λ_2⟩ ⟨Λ_2|R^j|0⟩
       +   ⟨ 0|R^i|Λ_1⟩ ⟨Λ_1|R^j|Λ_2⟩ ⟨Λ_2|δ^ij R^2-3 R^i R^j|0⟩] 3/10 ⟨ R^2⟩^2 
      f(Λ_1/Λ_2)  ,

with

    f(x)  = x ln(1+1/√(x)) - √(x) - ln(1+√(x))  .

Since f weakly depends on its argument, one can replace the argument of f in Eq. (<ref>) by
its averaged value to obtain

    β = f(⟨Λ_1/Λ_2⟩)  .

For the estimation of β we will assume that ⟨Λ_1/Λ_2⟩ = 1, 2, 1/2
and thus obtain β = -1.0 (0.2).

There is an additional contribution E'_Q that includes the finite proton size. It is obtained by
inserting the proton electric formfactor in the Coulomb interaction in Eq. (<ref>) and expanding
in R_i up to the second order,

    E'_Q =  ϕ^2(0) α^3 m^2 r_s^2 r^2_p 4 π/9 (13/3 + 1/2 ϵ
         -5 ln 2 -1/r_s^2 ⟨ 0|R⃗ ln(H-E)/m R⃗|0⟩) .

The remaining contributions E_L1 is given by Eq. (<ref>) with the deuteron radii, whereas
E_H1(p) is given by Eq. (<ref>) with the proton charge radii.


Adding all parts together, the total nuclear structure contribution E^(6)_1 ∝ R^4  is

    E^(6)_1(nS)    =   α^6 m^5 4/9 n^3 [r_d^4 (-1/n + γ -lnn/2+Ψ(n)
            + lnα) + r_dd^4 3/20 n^2
          
      -r_ss^4(3/5 ln2-9/10)
      +r_s^4 (1/3 - 3/2 ln2 - 1/2 ln⟨ E ⟩_1/m + 2 β) 
          + r_p^4 (2+γ +ln(r_p1 m))
            + r_s^2 r^2_p (10/3-5 ln 2 - ln⟨ E⟩_2/m)]  .

Again, the part of the elastic contribution depending on the effective deuteron radius r_d1
is not present in  total E^(6)_1. The expression for the separate  inelastic contribution is

    E^(6)_ 1,pol   =      E^(6)_1 - E^(6)_ 1,fns = α^6/n^3 m^5 δ_l0 [
      -r_ss^4 (4/15 ln2-2/5)
      +2/9 r_s^4 (2/3 - 3 ln2 - ln⟨ E ⟩_1/m + 4 β) 
           + 4/9 r_p^4 (2+γ +ln(r_p1 m))
            + 4/9 r_s^2 r^2_p (10/3-5 ln 2 - ln⟨ E⟩_2/m)
            -4/9 r_d^4 (2+γ +ln(r_d1 m))].




 §.§ Total inelastic part


The sum E^(6)_ 1,pol + E^(6)_ 2,pol, as given by Eqs. (<ref>) and (<ref>), is
the total three-photon exchange inelastic nuclear structure contribution, which is the main result
of this work. It should be pointed out that several approximations have been made in our derivation
of this result. First, we ignored the magnetic dipole and the electric quadrupole moments of
deuteron. Second, we neglected the higher orders in (H_N-E_N)/m. These approximations contribute
to the uncertainty of the inelastic part, which we estimate as 10%.

The remaining part of the total three-photon exchange nuclear structure contribution of
Eq. (<ref>) is the contribution due to the interaction with individual nucleons E^(6)_ np.
We have little knowledge about the inelastic three-photon exchange between the muon
and the proton but we expect it could be accounted for in terms of the same effective radii
r_p1 and r_p2 as in the elastic part. We estimate the uncertainty associated with
E^(6)_ np in μH and in μD by applying Eq. (<ref>) to the proton (r_C →
r_p) and making the following substitution,

    ln r_p1   →   ln r_p1 ± 1 , 
    ln r_p2   →   ln r_p2 ± 1 .

It should be mentioned that E^(6)_ np does not contribute to the μD-μH
isotope shift difference.




§ INELASTIC THREE-PHOTON EXCHANGE CORRECTION IN ORDINARY DEUTERIUM

The total inelastic nuclear structure α^6 correction for ordinary deuterium is split into
three parts,

    E^(6) = E_L2 + E_R + E_H2(p) + E^(6)_ np ,

where E_L2 is given by Eq. (<ref>) with the deuteron radii, and E_H2(p) by Eq.
(<ref>) with the proton radii, while E_R is a Coulomb correction to the electric dipole
polarizability, as given by Eq. (<ref>), and E^(6)_ np is the correction due to
the interaction with individual nucleons.

Assuming that H_N-E_N is much larger than the electron mass m, we obtain

    E_R =
    -π α^3 ϕ^2(0) 1/d ⟨ϕ_N|R^ k[
     1/2 ϵ + 5/2 - 2 ln2 (H_N-E_N)/m]
     R^ k|ϕ_N⟩ + O(m/H_N-E_N)  .

We note that the neglected O(m/(H_N-E_N)) terms do not vanish for the l>0 states. Therefore,
the nuclear polarizability correction does not vanish for the l>0 states, but it is additionally
suppressed by the ratio of the electron mass to the nuclear excitation energy.

The correction due to the interaction with individual nucleons E^(6)_ nucleon is expected
to be small and, moreover, it cancels out in the eD-eH isotope shift difference. In order to
estimate the three-photon exchange of the bound electron with the proton, we use the same
argumentation as in Ref. <cit.> to obtain

    E^(6)_ np(p) = 2 π/3 α^3 ϕ^2(0) r_p^2 lnE̅_p /m ,

and assume the uncertainty of 100%. The above correction is proportional to the squared charge
radius, so the corresponding contribution for the neutron is negligible.

Our final result for the three-photon exchange nuclear structure correction in deuterium is given
by

    E^(6)(nS)    =    -α^6 m^3 2/3 n^3 [r_d^2 (9/4n^2 - 3 -1/n + γ -lnn/2+Ψ(n)
        +lnα)
          
       +r_s^2 (23/12 - ln2 ⟨ E ⟩_2/m)
       +r_p^2 (γ + ln(r_p2 m))] + E^(6)_ np .

Separately, the inelastic part E^(6)_ pol = E^(6) - E^(6)_ fns is

    E^(6)_ pol = -2 α^6/3 n^3 m^3 δ_l0 [
        r_s^2 (23/12 - ln2 ⟨ E ⟩_2/m)
       +r_p^2 (γ + ln(r_p2 m))
       -r_d^2 (γ + ln(r_d2 m))]
       + E^(6)_ np .

We note that the fermion mass m cancels exactly in the expression in square brackets in the above equation.



§ RESULTS AND SUMMARY



Our numerical results for the three-photon exchange nuclear structure corrections are presented in
Table <ref>. The elastic part has been calculated with the exponential model of the
nuclear charge distribution. It is displayed in the table separately for the comparison with the
literature results. This part does not bear any uncertainty because its dependence on the charge
distribution model cancels out exactly in the sum with the inelastic part. We observe a reasonable
(although not perfect) agreement with the literature results summarized in Table <ref>.

The inelastic three-photon exchange nuclear structure correction was calculated only for the
electronic and muonic deuterium atoms; the corresponding results are presented in
Table <ref>. We find that the inelastic contribution for deuterium is of opposite sign
as compared to its elastic counterpart and changes significantly the total mα^6 nuclear
structure contribution. In the case of eD, the change is of about 30%, while for μD, the
inelastic part reverses the sign of the overall contribution. For electronic and muonic hydrogen, we
present only estimations for the inelastic three-photon exchange contribution.

Our results for the three-photon exchange nuclear structure corrections affect determinations of
the hydrogen-deuterium nuclear charge radii differences derived from the spectroscopic observations
of the isotope shifts in electronic and muonic hydrogen and deuterium
<cit.>. For the electronic H-D isotope shift of the 1S–2S
transition, our result shifts the total theoretical prediction by 0.8 kHz, which is slightly larger
than the theoretical error of 0.6 kHz assumed in Ref. <cit.>. There are, however,
further corrections to the summary of theoretical contributions presented in
Ref. <cit.>, so we had to update it. Our review of the present status of the theory of
the H-D isotope shift described in Appendix <ref> leads us to the updated result for the
nuclear charge radius difference determined from the measurement of the H-D isotope shift of the
1S–2S transition <cit.>,

    δ r^2 [] ≡ r_d^2-r_p^2  = 3.820 70 (31)  fm^2 ,

which agrees with but is twice as accurate as the previous value of 3.820 07 (65) fm^2
obtained in Ref. <cit.>.

For muonic hydrogen and deuterium, our result for the inelastic three-photon exchange nuclear
structure contribution to the 2P_1/2–2S transition energy of 0.008 75 (88) meV shifts
the deuteron-proton charge radius difference determined in Ref. <cit.> by 0.0014 fm^2,
with the result

    δ r^2 [] ≡ r^2_d-r^2_p  = 3.8126 (34) ^2 .

The results derived from the electronic and muonic atoms disagree by about 2 σ, which
confirms the discrepancy previously observed in Ref. <cit.>.









In summary, we have calculated the complete three-photon exchange O(α^2) nuclear structure
correction to energy levels and the isotope shift of hydrogen-like muonic and electronic atoms. Our
formula for the elastic contribution is valid for an arbitrary hydrogenic system and is much
simpler than corresponding formulas in the literature <cit.>. The inelastic part has
been derived for muonic and electronic deuterium  only. Calculations of the three-photon inelastic
contribution for He^+ and heavier elements are possible but are more complicated. At the same time,
one may expect the inelastic contribution to be as large as the elastic part, which is a sizeable
correction in He^+, about 1% of the total nuclear nuclear size effect.


The authors acknowledge support from the National Science Center of Poland
(Grant No. 2012/04/A/ST2/00105).
V.P. acknowledges support also from the Czech Science Foundation - GAČR (Grant No. P209/18-00918S), and
V.A.Y. acknowledges support also from the Ministry of Education and Science of the Russian Federation
Grant No. 3.5397.2017/6.7.





10

pohl:10 R. Pohl et al.,
Nature (London) 466, 213  (2010).

pohl:16 R. Pohl et al.,
Science 353, 669 (2016).

pohl:priv R. Pohl,
Private communication, 2018.

friar:97:b J. L. Friar and G. L. Payne,
Phys. Rev. C 56, 619 (1997).

pachucki:11:mud K. Pachucki,
Phys. Rev. Lett. 106, 193007 (2011).

pachucki:15:mud K. Pachucki and A. Wienczek,
Phys. Rev. A 91, 040503 (2015).

hernandez:14 O. Hernandez, C. Ji, S. Bacca, N. N. Dinur, and N. Barnea,
Phys. Lett. B 736, 344  (2014).

hernandez:18 O. Hernandez, A. Ekström, N. Nevo Dinur, C. Ji, S. Bacca, and N. Barnea,
Phys. Lett. B 778, 377 (2018).

friar:79:ap J. L. Friar,
Ann. Phys. (NY) 122, 151  (1979).

birse:2012 M. C. Birse and J. A. McGovern,
Eur. Phys. J. A 48, 120 (2012).

antognini:13:ap A. Antognini, F. Kottmann, F. Biraben, P. Indelicato, F. Nez, and
    R. Pohl,
Ann. Phys. (NY) 331, 127  (2013).

krauth:16:ap J. J. Krauth, M. Diepold, B. Franke, A. Antognini, F. Kottmann, and R. Pohl,
Ann. Phys. (NY) 366, 168  (2016).

borie:12 E. Borie,
Ann. Phys. 327, 733  (2012).

khriplovitch:00 I. Khriplovich and R. A. Sen'kov,
Physics Letters B 481, 447 (2000).

demissie:16 B. Demissie and H. W. Grießhammer,
https://arxiv.org/abs/1612.0735 (2016).

mohr:16:codata P. J. Mohr, D. B. Newell, and B. N. Taylor,
Rev. Mod. Phys. 88, 035009 (2016).

friar:97:a J. L. Friar and G. L. Payne,
Phys. Rev. C 55, 2764 (1997).

indelicato:13 P. Indelicato,
Phys. Rev. A 87, 022501 (2013).

wiringa:95 R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla,
Phys. Rev. C 51, 38 (1995).

parthey:10 C. G. Parthey, A. Matveev, J. Alnis, R. Pohl, T. Udem, U. D. Jentschura,
  N. Kolachevsky, and T. W. Hänsch,
Phys. Rev. Lett. 104, 233001 (2010).

jentschura:11 U. D. Jentschura, A. Matveev, C. G. Parthey, J. Alnis, R. Pohl, T. Udem,
  N. Kolachevsky, and T. W. Hänsch,
Phys. Rev. A 83, 042505 (2011).

sick:14 I. Sick,
Phys. Rev. C 90, 064002 (2014).

sick:08 I. Sick,
Phys. Rev. C 77, 041302 (2008).

franke:17 B. Franke, J. J. Krauth, A. Antognini, M. Diepold, F. Kottmann, and R. Pohl,
Eur. Phys. J. D 71, 341 (2017).

diepold:18 M. Diepold et al.,
https://arxiv.org/abs/1606.05231 (2018).

pachucki:06:hesinglet K. Pachucki,
Phys. Rev. A 74, 022512 (2006).

shabaev:04:DKB V. M. Shabaev, I. I. Tupitsyn, V. A. Yerokhin, G. Plunien, and G. Soff,
Phys. Rev. Lett. 93, 130405 (2004).

yerokhin:15:recprl V. A. Yerokhin and V. M. Shabaev,
Phys. Rev. Lett. 115, 233002 (2015).

heisse:17 F. Heiße, F. Köhler-Langes, S. Rau, J. Hou, S. Junck, A. Kracke,
  A. Mooser, W. Quint, S. Ulmer, G. Werth, K. Blaum, and S. Sturm,
Phys. Rev. Lett. 119, 033001 (2017).

yerokhin:16:recoil V. A. Yerokhin and V. M. Shabaev,
Phys. Rev. A 93, 062514 (2016).

davydychev:93 A. I. Davdychev and J. B. Tausk, Nucl. Phys. B 397, 123 (1993).










§ DIMENSIONAL REGULARIZATION FOR BOUND STATES

 The principles of dimensional regularization for bound states have been described in
Ref. <cit.>. Here we only present formulas without derivation which have been
used in the presented calculations. The dimension of space is assumed to be d=3-2 ϵ. The
surface area of the d-dimensional unit sphere is

    Ω_d = 2 π^d/2/Γ(d/2) .

The Coulomb potential in d dimensions is of the form

    V(r) = ∫d^d k/(2 π)^d 4 π/k^2 e^i k⃗·r⃗ = C_1/r^1-2 ϵ ,

where

    C_1 = π^ϵ-1/2 Γ(1/2-ϵ)  .

The elastic contribution involves another potential of the form

    V^(2)(r) = ∫d^d k/(2 π)^d 4 π/k^4  e^i k⃗·r⃗
    = C_2 r^1+2 ϵ,

where

    C_2 = 1/4 π^ϵ-1/2 Γ(-1/2-ϵ).

Futher, we used the following integration formulas

    ∫_Λ d^d r [ V(r)]^3    =   
     -[(4 π)^ϵ Γ(1+ϵ)]^2 4 π  [1/4 ϵ + 1/2+γ+ln(Λ)] ,
    ∫d^dk/(2π)^d1/k^2α1/(k-q)^2β   =   [q^2]^d/2-α-β/[4π]^d/2Γ(α+β-d/2) Γ(d/2-α) Γ(d/2-β)/Γ(d-α-β) Γ(α) Γ(β) ,

and <cit.>

    I    =   ∫d^d k/(2 π)^d ∫d^d q/(2 π)^d 1/[k^2]^n_1 1/[(k-q)^2+m_2^2]^n_2 1/[q^2+m_3^2]^n_3
       =   m_3^2 (d-n_1-n_2-n_3)/(4 π)^d Γ(d/2-n_1) Γ(n_1+n_2-d/2) Γ(n_1+n_3-d/2) Γ(n_1+n_2+n_3-d)/Γ(2 n_1+n_2+n_3-d) Γ(n_2) Γ(n_3) Γ(d/2)
          ×       _2F_1(n_1+n_2+n_3-d,n_1+n_2-d/2,2 n_1+n_2+n_3-d,1- m_2^2/m_3^2) .




§ NUMERICAL VERIFICATION OF THE ELASTIC CONTRIBUTION

 The finite nuclear size (fns) correction can be calculated numerically to all orders
in , by computing the energy eigenvalue of the Dirac equation with the extended-size nuclear
potential and subtracting the analytical point-nucleus result. Knowing the leading α^4 and
α^5 fns corrections analytically, we also can identify the higher-order fns residual from
the numerical all-order results.

The main problem in determining the fns correction numerically is that the corresponding effect is
very small for light electronic atoms. So, for the 2s state of hydrogen, the relativistic
O(α^2) fns correction yields a 1× 10^-13 fraction of the binding energy. In order
to  make an extensive comparison between the numerical and analytical approaches, we performed
numerical calculations for Z as low as Z = 0.25. To make sure that possible numerical
uncertainties do not interfere with the comparison, we determined the binding energies with a
20-digit numerical precision and made sure that the nuclear charge distribution model is exactly
the same as in analytical calculations.

In order to compute the eigenvalues of the Dirac equation, we use the Dual Kinetic Balance method
<cit.> with the finite basis set of B-splines. Because of high accuracy demands,
we implemented this method in quadruple (about 32 digits) arithmetics, similarly as it was done
recently in calculations of the recoil corrections <cit.>. About 200-250 basis
functions were sufficient to reach the required 20-digit numerical accuracy for the binding
energies.

We obtained the relativistic fns correction E_ fns^(6+) that contains contributions of order
()^6 and higher as

    E_ fns^(6+) = E_ fns - E_ fns^(4) - E_ fns^(5) ,

where E_ fns is determined numerically by solving the Dirac equation, E_ fns^(4) is
given by Eq. (<ref>), E_ fns^(5) is given by Eq. (<ref>), and r_Z is evaluated for
the same nuclear model as in the numerical calculation.

The comparison of our all-order numerical results for E_ fns^(6+) with the analytical
()^6 result E^(6) given by Eq. (<ref>) is presented in Fig. <ref>. Both numerical
and analytical results are obtained with the exponential model of the nuclear charge distribution
(see Table <ref>). We plot the scaled function with the leading , r_C, and n
dependence removed,

    F_ fns^(6+) = E_ fns^(6+)/m^3 (Zα)^6 r_C^2/n^3 .

As can be seen from the figure, agreement between the numerical and analytical results is
excellent.







§ HYDROGEN-DEUTERIUM 1S–2S ISOTOPE SHIFT



In this section we update the summary of all available theoretical contributions for the eH-eD
isotope shift of the 1S–2S transition frequency reviewed previously by Jentschura et al.
<cit.>. We use the following values for the fundamental constants, the fine-structure
constant,

    α^-1 = 137.035 999 139 (31) ,

and the Rydberg constant,

    R_∞c = 3.289 841 960 355 (19)× 10^15 z ,

from CODATA 2014 <cit.>. The electron-proton mass ratio we take from the recent
measurement by Heiße et al. <cit.>,

    m_p/m_e = 1 836.152 673 346 (81) .

We note that this value is twice as accurate but 3 σ off from the CODATA 2014 value
<cit.>. For the electron-deutron mass ratio we use the CODATA value
<cit.>,

    m_D/m_e = 3 670.482 967 85 (13) .


In the present section we follow the notations and conventions of Ref. <cit.>. We
will not repeat the full review of the theory but only indicate the entries therein that need to be
updated. The changes are as follows.

(i) The updated result for the leading (Dirac) contribution to the isotope shift (Eq. (28) of
Ref. <cit.>) is

    Δ f_ i = 671 004 071.107 (64)  ,

the change being due to the updated values of the electron-nucleus mass ratios.

(ii) Our present result for the two-photon exchange nuclear structure correction specified by
Eqs. (<ref>) and (<ref>),

    E^(5)(,1S2S) = 19.12 (20)  ,

replaces the sum of Δν_9 given by Eq. (40) of Ref. <cit.> and E_
NS,(b) given by Eq. (45) therein, amounting to 19.11 (2) kHz.

(iii) Our present result for the three-photon exchange nuclear structure correction,

    E^(6)(,1S2S) = -2.029 (41)  ,

replaces the sum of E_ NS,(c) given by Eq. (47) of Ref. <cit.> and Δν_11=± 0.5 kHz given by Eq. (43) therein, amounting to -2.828 ± 0.5 kHz.

(iv) The entry for the higher-order pure recoil ν_5 = -3.41 (32) kHz (Eq. (33) of
Ref. <cit.>) is replaced by the complete all-order (in ) result by Yerokhin and
Shabaev <cit.>. The corresponding correction to the energy is

    δ E(nS) = m^2/M()^5/π n^3[(4 ln 2 - 7/2)π  + ()^2  G_ rec
    + δ_ fns P ] ,

where G_ rec(1S,Z=1) = 9.720 (3), G_ rec(2S,Z=1) = 14.899 (3), δ_ fns
P(nS, H) = -0.000 184 (1) in the case of hydrogen and δ_ fns P(nS, D) =
-0.000 786 (6) for deuteron <cit.>. In the result, the
updated contribution is

    Δν_5 = -3.058  .


(v) For the radiative recoil contribution (Eq. (36) of Ref. <cit.>), we use the
estimation of uncertainty from Ref. <cit.>, which is about three times larger than
the one of Ref. <cit.>,

    Δν_6 = -5.38 (35)  .


The final theoretical value of

    Δ f_ th = 670 999 567.88 (42)  ,

replaces the previous result Δ f_ th() =
670 999 566.90 (89). Combining the theoretical value Δ f_ th with
the experimental result from <cit.>, we obtain the updated result for the
mean-square charge-radii difference

    r^2_d - r^2_p = 3.820 70 (31) ^2 ,

which agrees with but is twice as accurate as the previous value of 3.820 07 (65) fm^2
<cit.>.

