Pionless Effective Field Theory Evaluation of Nuclear Polarizability in Muonic Deuterium
PPionless Effective Field Theory Evaluation of NuclearPolarizability in Muonic Deuterium
Samuel B. Emmons,
1, 2
Chen Ji, ∗ and Lucas Platter
2, 41
Department of Mathematics, Physics, and Computer Science,Carson-Newman University, Jefferson City, TN 37760, USA Department of Physics and Astronomy,University of Tennessee, Knoxville, TN 37996, USA Key Laboratory of Quark and Lepton Physics, Institute of Particle Physics,Central China Normal University, Wuhan 430079, China Physics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA (Dated: September 18, 2020)
Abstract
We calculate the longitudinal structure function of the deuteron up through next-to-next-to-leading order in the framework of pionless effective field theory. We use these results to computethe two-photon polarizability contribution to Lamb shift in muonic deuterium, which can be utilizedto extract the nuclear charge radius of the deuteron. We present analytical expressions order-by-order for the relevant transition matrix elements and the longitudinal structure function, and wegive numerical results for the corresponding contributions to the Lamb shift. We also discuss theimpact of relativistic and other higher-order effects. We find agreement with previous calculationsand explain the accuracy of our calculation. ∗ [email protected] a r X i v : . [ nu c l - t h ] S e p . INTRODUCTION The deuteron is the simplest bound nuclear system and is made up of only two nucle-ons. It is the perfect testing ground for new ideas in nuclear theory since calculations arerelatively simple for this system and can be compared with decades of experimental data.The electromagnetic properties of a nucleus provide insights into its size, shape, and itscontinuum properties. These are therefore sensitive observables that test our understandingof short- and long-range charactistics of the nucleon-nucleon interaction.The charge radius is one of the most elementary electromagnetic observables. It can bemeasured using either elastic electron-nucleus scattering or laser spectroscopy. An analy-sis of the world-averaged electron-deuteron scattering data determined the deuteron root-mean-square (rms) charge radius to be r d = 2 . r d and an analysis based on the world-averageddeuteron spectroscopy data yields r d = 2 . µ -d) found a smaller radius, r d = 2 . . σ from the electronic deuterium ( e -d) spectroscopic result. This apparent differenceof r d in e -d and µ -d coincides with the original proton radius puzzle that spurred muchtheoretical and experimental work after a 7 σ deviation was discovered between the protoncharge radius extracted from the µ -H Lamb shift, r p = 0 . r p = 0 . µ -d spectroscopy is sensitive tothe two-photon exchange (TPE) contribution to the atomic 2 S -2 P level spacing [11]. Theelectromagnetic polarization of the nucleus caused by the muon leads to a distortion ofthe muon-nucleus wave function and affects thereby the atomic spectrum. In the TPEprocess, the nucleus is virtually excited and de-excited by the exchange of two photonswith the muon. Therefore, TPE depends not only on the bound-state properties of thedeuteron, but also on the nucleon-nucleon continuum scattering state and the form of the2lectromagnetic current. TPE plays a crucial role in connecting the physics of nuclearstructure to photonuclear reactions, and the accuracy of related calculations depends howwell known the nuclear Hamiltonian is.The effects of TPE on µ -d observables were originally calculated in Refs. [12–15] andwere recently revisited with improved accuracy using different nuclear models [16–21]. Thecalculations done with Argonne V18, chiral effective field theory ( χ EFT), and zero-rangeapproximated (ZRA) nucleon-nucleon interactions show good agreement with each other,demonstrating the high predictive power and accuracy of the state-of-the-art nuclear models.By analyzing statistical and systematic uncertainties on TPE calculations performed withinthe χ EFT framework, the uncertainty due to nuclear model dependence in these calculationswas probed [18, 19]. Furthermore, TPE in muonic deuterium was also considered using adispersion relation analysis of the scattering data [22]. This approach was also shown to agreewell with the aforementioned nuclear model calculations. Additional work has extended theevaluation of TPE to other light muonic atoms and ions, i.e. , µ H, µ He + , and µ He + [23–27].In this work, we make use of the work of Rosenfelder and Leidemann [13, 15] and therebyalso more recent calculations that use a state-of-the-art nuclear Hamiltonian [16, 18, 20].However, instead of χ EFT, we will use pionless effective field theory ( /π EFT). Calculations inthis framework can be expanded on an order-by-order basis in an expansion parameter that isproportional to the range of the nuclear interaction R over the two-nucleon scattering length a . The momentum scale of the processes considered within this approach are assumed to beof order 1 /a . /π EFT is suited for processes whose momentum scales are well below the pionmass, as in µ -d. For such processes, this approach offers a systematic expansion that is orderby order renormalizable, results whose regulator dependence is transparent and understood,and for the two-nucleon systems frequently also analytic results that reveal the dependenceon physical parameters directly [28]. Additionally, since the muon is approximately 200times heavier than the electron, it orbits closer to the nucleus and may be considered asapproximately non-relativistic, enabling us to neglect relativistic effects at the EFT orderwe consider. Compared with χ EFT, the order-by-order renormalizability and regulatorindependence in /π EFT provides a rigorous systematic uncertainty estimation which is model-independent. Further, the smaller number of parameters in /π EFT make it a powerful toolto explore few-body universality in few-nucleon systems [29].3he organization of this work is as follows. We introduce the basic equations that relatethe inelastic structure functions and electric form factors to the TPE effect in muonic deu-terium in Sec. II. In Sec. III, we explain the pionless EFT Lagrangian that will be utilized inthis work. We then show in Sec. IV how diagrammatic calculations can be used to obtain theinelastic structure function. In Sec. V, we present results for the TPE energy shift obtainedfrom our calculations and compare it to previous calculations. We conclude with a summaryand a discussion of possible extensions of this work.
II. THEORY OF TWO-PHOTON EXCHANGE CONTRIBUTIONS
The two-photon exchange contributions in muonic deuterium can be separated into apart depending on the structure of the atomic nucleus and another that depends on theinternal dynamics of the single nucleon. In this paper, we focus only on the nuclear two-photon exchange contribution to the former, labeled δ A TPE , by considering single nucleons aspoint-like particles. δ A TPE consists of the elastic and the inelastic parts δ A TPE = δ A Zem + δ A pol , (1)where the elastic part δ A Zem corresponds to the nuclear third Zemach moment contributionfirst derived for light muonic atoms by Friar [30] as a nuclear finite-size contribution of order α , where α = 1 / .
036 denotes the fine structure constant. It is given by [31, 32] δ A Zem = − m r α (cid:10) R (cid:11) (2) = − m r α π (cid:90) ∞ dqq (cid:20) F ( q ) − q (cid:104) R (cid:105) (cid:21) , (2)where m r is the muon-deuteron reduced mass, F E ( q ) is the deuteron electric form factor,and (cid:104) R (cid:105) = − ∂F E ( q ) /∂q ] q =0 , where the derivative is taken with respect to a low- q expansion of F E ( q ).The inelastic contribution δ A pol in Eq. (1) is due to the electric polarization of the nucleusin which the deuteron is virtually excited by exchanging two photons with the muon. Thisis related to the integral of the forward virtual Compton amplitude that can be written interms of the nuclear inelastic structure functions [13]. The polarizability term δ A pol may befurther separated into longitudinal and transverse parts as [13, 15] δ A pol = δ A pol ,L + δ A pol ,T , (3)4here δ A pol ,L and δ A pol ,T are defined respectively as [15] δ A pol ,L = − α | φ (0) | (cid:90) ∞ dq (cid:90) ∞ ω th dωK L ( ω, q ) S L ( ω, q ) , (4) δ A pol ,T = − α | φ (0) | (cid:90) ∞ dq (cid:90) ∞ ω th dω [ K T ( ω, q ) S T ( ω, q ) + K S ( ω, q ) S T ( ω, . (5)The variables ( ω, q ) are the four-momentum carried by the exchanged photon, where q = | q | , φ (0) = (cid:112) α m r / π is the atomic 2 S -state wave function at origin. To ensure that only theinelastic regime is considered, we have ω th ≥ B d + q / m N . B d = 2 . m N = 938 .
92 MeV is two times the proton-neutron reduced mass,and q / m N is the recoil energy of the nucleus. S L and S T are the longitudinal and transversedeuteron inelastic structure functions, respectively.The longitudinal integration kernel in Eq. (4) is given by K L ( ω, q ) = 12 E q (cid:20) E q − m µ ) ( ω + E q − m µ ) − E q + m µ ) ( ω + E q + m µ ) (cid:21) , (6)where m µ denotes the muon mass, and E q = (cid:112) q + m µ retains the relativistic kinematicsof the muon. The transverse and seagull kernels of Eq. (5) are provided in Ref. [15]. The seagull term is required to ensure the gauge invariance and to cancel the infrared singularitynear q = 0 in the transverse term. In the Coulomb gauge utilized in this paper, the seagullterm contributes only to δ A pol ,T [13]. In the non-relativistic limit q (cid:28) m µ , the longitudinalkernel in Eq. (6) is approximated in q/m µ expansion by K NRL = 1 q ( ω + q / m µ ) . (7)This kernel’s higher-order terms emerge as relativistic corrections.We note that the expressions for δ A pol ,L and δ A pol ,T in Eqs. (4) and (5) were multipliedby an additional factor R ( µ ) = 0 . α effect, and is thus neglected in this paper for consistency sincethe evaluation of δ A pol and δ A Zem is of order α .In this work, we compute only relevant contributions up through next-to-next-to-leadingorder in /π EFT. Following Rosenfelder [13], we relate the longitudinal part of the structurefunction S L to the transition matrix element M by S L ( ω, q ) = (cid:90) d p (2 π ) δ ( ω − B d − q m N − p m N ) |M| , (8)5here |M| is the squared transition matrix element for the electric density operator, be-tween the deuteron ground state and all intermediate excited states.In the framework of pionless effective field theory /π EFT, we will calculate the squaredmatrix element |M| relevant to the longitudinal deuteron structure function by consideringthe coupling of a single A Coulomb photon to the deuteron. We will discuss below thatcontributions arising from the transverse structure function do not contribute to the orderconsidered here. The deuteron charge form factor F E ( q ) has been evaluated using /π EFTin Ref. [33].
III. PIONLESS EFFECTIVE FIELD THEORY
In this section, we provide a brief overview of /π EFT and the partial-divergence subtraction(PDS) renormalization scheme. We include some details about both the on-shell and off-shell nucleon-nucleon scattering amplitudes that are needed in the transition matrix elementcalculations required for the structure function calculation of Eq. (8).
A. Lagrangian and Feynman rules
The nucleonic part of the /π EFT Lagrangian is given by [33] L = N † (cid:18) i∂ + ∇ m N (cid:19) N − C (cid:0) N T P i N (cid:1) † (cid:0) N T P i N (cid:1) + C (cid:104) ( N T P i N ) † (cid:16) N T P i ←→∇ N (cid:17) + h.c (cid:105) − C (cid:16) N T P i ←→∇ N (cid:17) † (cid:16) N T P i ←→∇ N (cid:17) . (9)where we included the EFT nucleon-nucleon contact interactions in the S -channel up tonext-to-next-to-leading order (NNLO) and P i = σ σ i ⊗ τ / √ S channel. Additionally, P i ←→∇ = P i −→ ∂ + ←− ∂ P i − ←− ∂P i −→ ∂ , and the low-energyconstants (LECs) C i are determined through renormalization by reproducing parameters inthe effective range expansion around the deuteron pole. The neutron-proton spin-tripletscattering phase shift is expanded as k cot δ t = − γ + ρ d k + γ ) + · · · , (10)where γ = √ m N B d denotes the deuteron binding momentum and ρ d = 1 .
764 fm is theeffective range. Assuming that the momentum scale of processes considered here is compa-rable to the deuteron binding momentum, the expansion parameter in /π EFT is γρ d ≈ . m N ≈ ( γρ d ) − ρ − d . Feynman rules corresponding to the La-grangian given in Eq. (9) give a nucleon propagator as S ( p , p ) = [ p − p / m N + iε ] − . Thefirst relativistic correction to a two-nucleon matrix element is the kinetic energy multipliedwith a term proportional to k /m N . Using k ∼ γ and the aforementioned counting for thenucleon mass we see that this correction enters at O (( γρ d ) ).Under the /π EFT expansion, the LECs are expanded analogously by C = C , − + C , + C , ,C = C , − + C , − ,C = C , − , (11)where for a coefficient C n,m , n denotes the power of momentum in the contact term and n + m + 1 indicates the /π EFT order at which C n,m emerges. In the power-divergence-subtraction (PDS) scheme, the expanded LECs of Eq. (11) are given by [33, 34] C , − = − πm N µ − γ ) , C , − = 2 πm N ρ d ( µ − γ ) ,C , = 2 πm N ρ d γ ( µ − γ ) , C , − = − πm N ρ d γ ( µ − γ ) ,C , = − πm N ρ d γ ( µ − γ ) , C , − = − πm N ρ d ( µ − γ ) , (12)where µ is the PDS renormalization scale.The electromagnetic interaction with the nucleon field is included by replacing the four-gradient with the minimally-coupled gauge covariant derivative D µ = ∂ µ + ie Q A µ , where Q =(1 + τ ) / τ acting in isospin spaceand A µ is the electromagnetic gauge field. The Lagrangian for the induced electromagneticinteraction is thus given by L EM = − eN † Q N A + ie m N (cid:104) N † Q ( −←−∇ + −→∇ ) N (cid:105) · A − e m N N † Q A · A QN , (13)where the last term is a two-photon coupling and yields the seagull term in the two-photonexchange, which contributes only to δ pol ,T when Coulomb gauge is adopted [13]. We alsodrop the two-nucleon current from the minimal coupling within the contact interactions,because it only enters at higher orders than the NNLO considered in this paper.The nucleon current from one-photon coupling is given by J µ ( p , p (cid:48) ) = e Q (cid:18) , p + p (cid:48) m N (cid:19) . (14)7 + A LO A LO FIG. 1. Diagrammatic representation of the Lippmann-Schwinger equation for the LO two-bodyscattering amplitude A LO . The round vertex represents insertion of the LO contact term. We see that the transverse current enters the transition matrix element with a factor of k/m N . This implies that the squared matrix element entering the calculation of the trans-verse structure factor is also O (( γρ d ) ). B. Leading order nucleon-nucleon amplitude
We denote the leading order amplitude in the triplet channel as A LO , which is showndiagrammatically in Fig. 1. This is the Lippmann-Schwinger equation and is solved byidentifying the expression in the figure as the iterative sum i A LO ( E ) = − iC , − [1 − I ( E ) A LO ( E )] , (15)where E is the two-nucleon energy in the center of mass frame and I indicates the loopintegral defined in Eq. (A1). I is solved using the PDS scheme [34] and its dependence onthe renormalization scale µ is given in Eq. (A1). The leading-order amplitude A LO ( E ) = − πm N γ + ip , (16)is obtained using Eq. (15) and taking the expression of C , − in Eq. (12), where p = √ m N E is the nucleon-nucleon on-shell relative momentum.8 + + + A LO A LO A LO A LO A NLO
FIG. 2. Diagrammatic representation of the NLO amplitude calculated perturbatively. Thesquare vertex indicates the insertion of an NLO vertex rule.
C. Next-to-leading order amplitude
The diagrams required to evaluate the next-to-leading-order (NLO) are shown in Fig. 2and lead to the NLO half-off-shell amplitude i A NLO ( k, p ; E ) = − iC , [1 + i I i A LO ] − iC , − (1 + i I i A LO ) (cid:20) k + p i I i A LO (cid:21) = − i πm N ρ d γ + ip (cid:20) γ − ip + 12( γ − µ ) (cid:0) k − p (cid:1)(cid:21) . (17)where the loop integrals I , I and the amplitude A LO are evaluated at the center-of-massenergy E = p /m N , and k is the incoming-state momentum. We arrive at the second line ofEq. (17) using the Lippmann-Schwinger equation (15) and using I n = p n I from Eq. (A1).On the energy shell, the incoming momentum equals the on-shell momenta k = p . Thisyields the on-shell NLO amplitude A NLO ( p, p, E ) = − πρ d m N γ − ipγ + ip . (18) D. Next-to-next-to-leading order amplitude
To obtain the NNLO amplitude, we make use of a diagrammatic expression including allNNLO diagrams. After replacing the couplings in the resulting expression using Eq. (12),we find the half-off-shell NNLO amplitude to be A NNLO ( k, p ; E ) = − πm N ρ d γ + ip (cid:20) ( γ − ip ) + γ − ipγ − µ (cid:18) γ + ipγ − µ (cid:19) k − p (cid:21) . (19)The NNLO on-shell amplitude is found after setting k = p and is A NNLO ( p, p, E ) = − πρ d m N ( γ − ip ) γ + ip . (20)The µ -dependence of Eq. (19) is removed in the on-shell amplitude of Eq. (20).9 . Deuteron electric form factor In the vicinity of the bound-state pole p = iγ , the on-shell amplitude is given by A d ( E ) = − Z d E + γ /m N , (21)where Z d is the wave function renormalization factor related to ρ d and γ as Z d = 8 πγm N (1 − ρ d γ ) . (22)Eqs. (21) and (22) are needed in the calculation of the deuteron electric form factor thathas been calculated up to NNLO in /π EFT as [33] F E ( q ) = 11 − ρ d γ (cid:20) γq arctan q γ − ρ d γ (cid:21) . (23)The form factor of Eq. (23) is needed in evaluating the third Zemach moment contributionto the Lamb Shift in muonic deuterium given by Eq. (2). IV. DIAGRAMMATIC CALCULATION OF TRANSITION MATRIX ELEMENTS
In this section, we use /π EFT Feynman diagrams to evaluate the longitudinal transitionmatrix elements up to its NNLO contribution, which is needed for the evaluation of thelongitudinal structure function and the polarizability effect.
A. Transition matrix element at LO
At LO, the transition matrix element from the A -photon excitation of the deuteron isdepicted by the Feynman diagrams in Fig. 3. The first two diagrams, (a) and (b), arecontributions without final state interactions, while the third diagram contains final stateinteractions as indicated by a shaded grey oval symbolizing an insertion of the LO scatteringamplitude. Figure 3 also displays the kinematics we use in the calculation of this matrixelement. Note that the incoming deuteron is at rest, so the total energy of the initialtwo-nucleon state is therefore − B d . The final two-nucleon state has energy E = p /m N .We evaluate the plane-wave contribution by calculating the amplitude without projectingon specific spin-orbit coupled final states and are thereby summing up all electric multipolescontributing to the transition induced by the A photon. The plane-wave final states are10 − B d , ) p + q − p + q p + q − p + q ( − B d , ) A LO ( − B d , ) p + q − p + q (a) (b) (c) FIG. 3. Contribution to the leading order transition matrix element needed for the inelasticlongitudinal structure function. The wavy line denotes the coupling of the electromagnetic currentto the nucleon. represented by two nucleon spinors N T and N for each outgoing nucleon leg. N holds a tensorproduct of the two-component spinor for spin and the two-component spinor for isospin. Theinitial deuteron S state is projected out by the operator P i . The antisymmetrization istreated by including both diagrams where the photon either couples to the top or the bottomleg as shown in diagrams (a) and (b) of Fig. 3.For diagram (a) in Fig. 3, we obtain( i M a LO ) iαβ = i (cid:112) Z d iS ( − B d − ( − p + q / m N , p − q ) N Tα ( p + q Q P † i N β ( − p + q i (cid:112) Z d i ˜ M a LO N Tα ( p + q Q P † i N β ( − p + q . (24)where we defined the amplitude ˜ M a LO = − m N / [ γ + ( p − q / ]. The indices α and β onthe nucleon spinors denote spin and isospin basis indices, and i is the deuteron spin index.Note that a factor of 2 is included in Eq. (24) to account for both possible contractionsleading to that expression. Similarly, diagram (b) in Fig. 3 yields( i M b LO ) iαβ = i (cid:112) Z d i ˜ M b LO N Tα ( p + q P † i Q N β ( − p + q , (25)where the amplitude ˜ M b LO = − m N / [ γ + ( p + q / ].Furthermore, we evaluate the contributions with S final-state interactions at LO, usingthe diagram (c) shown in Fig. 3. The excitation from S bound state to the S scat-tering state is forbidden by the A photon. The iterative sum of final-state interactions isrepresented by the off-shell scattering amplitude P † i i A LO P i .11 LO + + A NLO (a) (b) (c)
FIG. 4. NLO diagrammatic contribution to the matrix element M . Final state S-wave interactionsare indicated by the grey blob. We have omitted one diagram without final state interactions wherethe photon couples to the lower outgoing leg. Using the Feynman rules we obtain for this diagram the transition matrix element( i M c LO ) iαβ = i (cid:112) Z d Tr (cid:104) P † i Q P j (cid:105) i A LO ( E ) N Tα ( p + q P † j N β ( − p + q × (cid:90) d l (2 π ) iS ( − B d + l , l − q iS ( − B d + l + ω, l + q iS ( − l , l − q . (26)Note that a factor of 8 is included for the possible number of contractions in the matrixelement evaluation. The loop integral arising in the calculation of diagram (c) in Fig. 3 welabel as J , and after evaluating the trace Tr[ P † i Q P j ] = δ ij / i M c LO ) iαβ = i (cid:112) Z d i ˜ M c LO N Tα ( − p + q P † i N β ( p + q , (27)where ˜ M c LO = 2 J A LO ( E ). The definition and the derivation of the J integral result aregiven in App. A 2.Below, at NLO and NNLO, we will use the short-hand notation introduced above andwork with the amplitudes ˜ M with √ Z d and the spinor pieces factored out. Amplitudes withsuperscripts a or b will continue to correspond to diagrams without final-state interactions,while those with superscript c correspond to diagrams at a given order that contain final-state interactions. B. Transition matrix element at NLO
Diagrams that contribute to the transition matrix element at NLO are shown in Fig. 4.The diagram with one NLO contact that has no final state interactions is given by diagram12a) of Fig. 4. After antisymmetrization, it leads to the expressions i ˜ M a NLO =2 (cid:20) ( − i C , − (cid:16) i I ( − B d ) (cid:0) p − q (cid:1) + i I ( − B d ) (cid:17) + ( − iC , ) i I ( − B d ) (cid:21) × iS ( − B d − ( − p + q / m N , p − q )= i ρ d m N ( µ − γ ) − (28)and i ˜ M b NLO = i ρ d m N ( µ − γ ) − , (29)where we used that −I ( − B d ) = γ I ( − B d ) and C , = γ C , − .Diagram (b) of Fig. 4 gives the contribution with no explicit contact term insertion, butwith one NLO final-state scattering amplitude. It leads to i ˜ M c, =2 ρ d i A LO ( E ) (cid:90) d k (2 π ) iS ( − B d − k m N , k ) iS ( − B d + ω − k m N , k + q ) × (cid:34) ( γ − ip ) + 12( γ − µ ) ( v − p ) (cid:35) , (30)where v = | k + q | and the half off-shell NLO amplitude A NLO ( v, p, E ) of Eq. (17) dependingon the loop momentum k was inserted. Using the definitions from the appendix for theloops that couple to the photon, we can write i ˜ M c, = 2 i A LO ( E ) ρ d (cid:20) ( γ − ip ) J + 12( γ − µ ) (cid:16) ˜ J − p J (cid:17)(cid:21) = ρ d i A LO ( E ) (cid:20) γ − ip ) J + m N π (cid:21) , (31)where ˜ J is defined in Eq. (A8) and we used the result of Eq. (A10) in App. A 3 to replace (cid:16) ˜ J − p J (cid:17) . Diagram (c) of Fig. 4 with one NLO operator and one LO scattering amplitudeyields an amplitude we label i ˜ M c, . Combining the sum of loop integrals and using theidentities C , = γ C , − and I ( − B d ) = − γ I ( − B d ) allows for this amplitude to be written˜ M c, = C , − I ( − B d )( J + γ J ) A LO ( E )= − ρ d m N µ − γ I ( E ) A LO ( E ) , (32) The diagram label (b) in Fig. 4 is not connected with the superscript c indicating final-state interactions. NLO A NNLO A LO A LO ( a ) ( b ) ( c )( d ) ( e ) ( f ) FIG. 5. NNLO diagrammatic contribution to the matrix element M that includes one insertionof NNLO operator or two insertions of NLO operator. The grey square and triangle denote thecollection of NLO and NNLO operators, respectively. These diagrams make zero contributionsafter renormalization. where the expressions for loop integrals from Eq. (A7) are used. The summation of Eqs. (31)and (32) yields ˜ M c NLO = ˜ M c, + ˜ M c, , = ρ d A LO ( E ) (cid:20) γ − ip ) J + m N π (cid:18) µ + ipµ − γ (cid:19)(cid:21) . (33) C. Transition matrix element at NNLO
In Fig. 5, we show the diagrams contributing at NNLO. Diagrams (a) and (b) in Fig. 5generate the NNLO contributions to M a and M b . The sum of these diagrams give nocontribution because ˜ M a NNLO = ˜ M b NNLO = 0. Diagrams (c) and (d) in Fig. 5 evaluate to zero,too. Diagrams (e) and (f) in Fig. 5 represent the non-zero contributions to the transitionmatrix element at NNLO. Diagram (e) gives the contributions with one insertion of oneNLO operator and one NLO amplitude. Its contribution to M c yields i ˜ M c, = C , − I ( − B d )( J + γ J ) ρ d i A LO ( E )( γ − ip ) , = i A LO ( E ) (cid:16) ρ d (cid:17) m N π ( γ − ip ) µ + ipµ − γ , (34)14here the expression of loop integral in Eq. (A7) is used. Inserting the half-off-shell NNLOamplitude in diagram (f), we obtain i ˜ M c, = 2 i A LO ( E ) (cid:16) ρ d (cid:17) (cid:20) ( γ − ip ) J + γ − ip γ − µ ) ( ˜ J − p J )(1 + γ + ipγ − µ ) (cid:21) = 2 i A LO ( E ) (cid:16) ρ d (cid:17) (cid:20) ( γ − ip ) J + m N π ( γ − ip )(1 + γ + ipγ − µ ) (cid:21) . (35)The regulator-dependence in the NNLO contribution is removed when Eqs. (34) and (35)are summed together as˜ M c NNLO = 2 A LO ( E ) (cid:16) ρ d (cid:17) ( γ − ip ) (cid:20) ( γ − ip ) J + m N π (cid:21) . (36)We note that the S-D wave mixing operator enters at NNLO in the /π EFT Lagrangian,and in principle gives a contribution to the transition matrix element at the same order.However, due to the orthogonality of S-wave and D-wave component, contributions fromthe D-wave projection to the squared matrix element does not interfere with the S-waveprojection. Therefore, the S-D mixing contributions enter at N4LO in the squared matrixelement, and we therefore do not consider this higher-order effect.
D. Matrix element squared
The calculation of the inelastic longitudinal structure function in Eq. (8) requires thesquared matrix elements and a sum over the outgoing nucleon-nucleon spin and isospinstates. This sum leads to traces of products of projection and charge operators. We carrythis out without projecting the outgoing state on a specific spin or isospin coupling. Wewrite |M| = 12 (cid:88) αβ |M aαβ + M bαβ + M cαβ | , (37)where the factor of 1 / |M| = Z d (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ M a − ˜ M b (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ M a + ˜ M b M c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (38)15here the amplitudes on the right hand side can be expanded order-by-order. The first termin Eq. (38) receives no NLO or NNLO corrections since ˜ M a NLO − ˜ M b NLO = 0, as can be seenfrom Eqs.(28) and (29), and ˜ M a NNLO = ˜ M b NNLO = 0.The second term in Eq. (38) receives LO, NLO, and NNLO contributions, but the NNLOcontribution arises only from ˜ M c NNLO , given in Eq. (36). We can therefore write the squaredamplitude up through NNLO as |M| = Z d (cid:40)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ M a LO − ˜ M b LO (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ M a LO + ˜ M b LO
2+ 2 A LO ( E ) (cid:34) J (cid:32) (cid:88) n =0 (cid:16) ρ d (cid:17) n ( γ − ip ) n (cid:33) + ρ d m N π (cid:16) ρ d γ − ip ) (cid:17)(cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:41) . (39) V. RESULTS
Analytical results for the inelastic longitudinal structure function can be calculated byinserting Eq. (39) into the integrand in Eq. (8). Once the structure function is obtained ateach order we consider in this work, we compute δ A pol ,L order-by-order. We extract the electricdipole contribution for comparison with other works that carry out an explicit multipoledecomposition [15, 20]. Additionally, we calculate the third inelastic Zemach moment term δ A Zem . A. Results for the longitudinal structure function
To demonstrate the order-by-order convergence of the /π EFT calculation of the longitudi-nal structure function S L , we extract the ρ d γ power dependence at each order. Besides theexplicit ρ d γ dependence in Eq. (39), we expand the the deuteron renormalization constant Z d in powers of ρ d γ as Z d = 8 πγm N (cid:0) ρ d γ + ρ d γ + · · · (cid:1) , (40)to include the order-by-order correction from the wave function renormalization.At leading order we find the inelastic structure function result S LO L ( q, ω ) = γpπm N (cid:20) m N m N ω − q p + 32 π m N ( γ + p ) × (cid:16)(cid:0) Re (cid:2) J (cid:3)(cid:1) − γp Re (cid:2) J (cid:3) Im (cid:2) J (cid:3) − (cid:0) Im (cid:2) J (cid:3)(cid:1) (cid:17)(cid:21) , (41)16 .5 5.0 7.5 10.0 12.5 15.0 17.5 20.0 22.5 [MeV] S L [ M e V ] q = 20MeV LONLONNLO Z d [MeV] S L [ M e V ] q = 50MeV LONLONNLO Z d FIG. 6. Longitudinal structure function S L as a function of ω at fixed momentum exchange q = 20 MeV (left panel) and q = 50 MeV (right panel). The dotted, dashed, dot-dashed lines givethe result of the LO, NLO and NNLO structure function, respectively. The solid line gives NNLO Z d -improved result. where p = (cid:112) m N ω − γ − q /
4. At NLO, the inelastic longitudinal structure function wecalculate as S NLO L ( q, ω ) = ρ d γS LO L ( q, ω ) + 8 ρ d γm N ( γ + p ) (cid:26) p Re [ J ] − γ Im [ J ]+ 4 πm N (cid:104) γp (cid:0) (Re [ J ]) − (Im [ J ]) (cid:1) + ( p − γ )Re [ J ] Im [ J ] (cid:105)(cid:27) . (42)Lastly, the N2LO part of the inelastic structure function is S NNLO L = ρ d γS NLO L ( q, ω ) + m N γρ d π ( γ + p ) (cid:26) p + 8 πm N (cid:104) γp Re [ J ] + ( p − γ )Im [ J ] (cid:105) + 16 π m N (cid:104) ( p − pγ ) (cid:0) Im [ J ] − Re [ J ] (cid:1) + (6 p γ − γ )Re [ J ] Im [ J ] (cid:105)(cid:27) . (43)The EFT convergence of the structure function is shown in Fig. 6, where S L ( ω, q ) isplotted as a function of ω by fixing q at 20 MeV and 50 MeV in the two different plots in thefigure. Calculations of S L ( ω, q ) at LO, NLO, and NNLO are compared in the plots, whichshow an order-by-order convergence in /π EFT.
B. Benchmark with dipole approximation results
As mentioned previously, we do not perform our calculation using a multipole decomposi-tion as in Refs. [15, 20]. Instead, the integration to determine δ A pol ,L in Eq. (4) was taken using17he complete inelastic longitudinal structure function of the nucleus, given order-by-orderin Eqs. (41), (42), and (43), that implicitly contains all inelastic multipole contributions.However, we extract the inelastic dipole excitation for comparison with previous literature.This is motivated by the fact that the dipole excitation gives the largest contribution toTPE.The electric-dipole excitation of the nucleus arises when the outgoing unbound nucleonsare in a spin-singlet isospin-triplet state. The antisymmetry of the wave function is preservedby the fact that the outgoing N N state has odd orbital angular momentum. The first termin the low-momentum approximation of the squared matrix element in Eq. (38) correspondsto a P wave between the outgoing nucleons and is |M| ≈ Z d (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˜ M a LO − ˜ M b LO (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≈ m N Z d ( p · q ) ( γ + p ) . (44)Inserting the matrix element |M| into a low- q truncated version of Eq. (8) by dropping therecoil energy q / m N gives the dipole part of the inelastic longitudinal structure function as S D ( | q | , ω ) = (cid:90) d p (2 π ) δ ( ω − γ m N − p m N ) |M| = Z d √ m N π ( ω − γ /m N ) / ω q . (45)Plugging the dipole structure function S D ( | q | , ω ) of Eq. (45) and the non-relativistic kernel K NRL from Eq. (7) into the TPE energy-shift equation given by Eq. (4) yields δ A pol ,D = − π α φ (0) Z d m / N √ m r π γ , (46)which matches the expression of the same contribution in Ref. [21]. In the above expression,we replaced m µ with the µ -d reduced mass m r to adjust the truncation of recoil correctionin the low-q approximation. The numerical evaluation of the dipole term given by Eq. (46)yields δ A pol ,D = − .
925 meV.
C. Numerical results for TPE contribution to Lamb shift
The TPE contribution to the Lamb shift in muonic deuterium can be calculated order-by-order in /π EFT to the desired precision. We do so up to NNLO in an expansion in theparameter ρ d γ . To evaluate the elastic TPE, Zemach term δ A Zem , we insert the electric form18 A pol ,L NR limit R limit /π EFT LO -0.962 -0.943 /π EFT NLO -1.346 -1.320 /π EFT NNLO -1.499 -1.470 Z d improved -1.605 -1.574dipole term -1.925 –ZRA ( η -less) -1.590 -1.553ZRA ( η -expansion) -1.590 -1.564 χ EFT ( η -less) -1.588 -1.562 χ EFT ( η -expansion) -1.590 -1.560TABLE I. The longitudinal polarizability δ A pol ,L is calculated in both non-relativistic and relativisticlimit. Results in /π EFT are calculated in this work at LO, NLO, NNLO and with the NNLO Z d -improvement approach. Results from other works in dipole approximation, ZRA and χ EFT areextracted from information in Refs. [20, 21, 27]. factor F E from Eq. (23) into the integral equation (2). δ A Zem is evaluated to be − .
362 meV,which is consistent with the calculation in Ref. [20].The longitudinal polarizability δ A pol ,L is firstly evaluated in the non-relativistic approxima-tion by using the non-relativistic kernel K NRL from Eq. (7) in the TPE sum rule. At leadingorder in /π EFT expansion, we obtain[ δ A pol ,L ] LO = − .
962 meV . (47)The NLO correction to δ A pol ,L consists of two parts, whose sum gives∆[ δ A pol ,L ] NLO = ( − .
393 + 0 . . (48)The dominant contribution arising from ρ d γ expansion in the Z d factor corresponding tothe first term in Eq. (48) is given by ρ d γ [ δ A pol ,L ] LO , and the rest is due to the NLO diagramsin Fig. 4. The dominant contribution to the NNLO correction on δ A pol ,L arises from the ρ d γ expansion in Z d , and is given by ρ d γ ∆[ δ A pol ,L ] NLO . The NNLO diagrams in Fig. 5 lead to acontribution that is signicantly smaller. The sum of the two NNLO contributions yields∆[ δ A pol ,L ] NNLO = ( − .
157 + 0 . . (49)19y summing up the contributions at each order, we obtain δ A pol ,L to be − .
962 meV, − . − .
499 meV at LO, NLO and NNLO respectively.The ( ρ d γ )-expansion of Z d clearly dominates the effective range corrections to δ A pol ,L andeffects of diagrams that include final state interactions beyond NNLO are suppressed. Theaccuracy of the calculation can therefore be improved above NNLO by simply using theresummed wave function renormalization Z d in the evaluation . We first separate out thediagrammatic contributions at LO, NLO and NNLO respectively by setting Z d equaling toits LO value 8 πγ/m N . The sum of the results are then multiplied with an additional factor1 / (1 − ρ d γ ) to match to the full wave function renormalization. By doing so, we have δ A pol ,L = 11 − ρ d γ ( − .
962 + 0 .
009 + 0 . − .
605 meV ± .
066 meV . (50) Z d -improvement can also be applied to the evaluation of the structure function, the resultof which is shown in Fig. 6. The uncertainty presented in Eq. (50) is ∼ ( ρ d γ ) , where new /π EFT parameters at higher orders are expected to enter.In Table I, we show the /π EFT results of δ A pol ,L computed at LO, NLO, NNLO and withthe NNLO Z d -improvement approach. We also shown the comparison with Ref. [20] inthe table, where the calculations were done using two different expansion methods, named η -less method and η -expansion method respectively. Using the η -less method, Ref. [20]evaluated δ A pol ,L in the multipole expansion of the charge operator and summing up multipolecontributions to high orders. In the η -expansion approach, δ A pol ,L in the same non-relativisticand point-proton limit is equivalent to δ A pol ,L = δ (0) D + δ (1) Z + δ (2) R + δ (2) Q + δ (2) D D with notationsgiven in Ref. [27]. Using nuclear potentials in the ZRA, η -less and η -expansion methods bothobtained δ A pol ,L = − .
590 meV [20], which is different by 0 .
8% from this work. Both /π EFTand ZRA are based on the effective range expansion. However, in ZRA, only the bound-state wave function renormalization is range corrected as in Z d . In the work by Hernandez et al. , the intermediate scattering states in the two-photon processes were treated as planewaves subtracted by the bound ground state when using the ZRA potential. Using χ EFTpotential, where final-state interactions are embedded in the diagonalization of the nuclearHamiltonian, δ A pol ,L obtained in η -less and η -expansion methods are respectively − .
588 and − .
590 meV [20]. The results from /π EFT agree with both ZRA and χ EFT calculationswithin the expected uncertainty at NNLO. Note that this is not Z-matching as introduced in Ref. [35] since we do not change the low-energy coeffi-cients used in the evaluation of the transition matrix element. η -expansion results of therelativistic δ A pol ,L in ZRA (or χ EFT), we need to add an additional relativistic correction0 .
037 meV from Ref. [21] (or 0 .
030 meV from Ref. [27]) to the non-relativistic value. Theresulting δ A pol ,L is − .
553 meV in ZRA and − .
560 meV in χ EFT. The relativistic δ A pol ,L calculated using η -less method (noted by ∆ L in Ref. [20]) is − .
562 meV in χ EFT. Wecalculate the relativistic η -less δ A pol ,L in ZRA by using the analytic matrix element fromRef. [20] and obtain − .
564 meV. Results with /π EFT in relativistic limit and comparisonare also shown in Table I. It indicates that agreement within a 1% discrepancy in δ A pol ,L isalso achieved in the relativistic limit. VI. SUMMARY
In this work, we have calculated the longitudinal structure function to NNLO in /π EFT.We have given analytic expressions for the squared matrix element required to calculate thestructure function. At NNLO, only two parameters, i.e. , the deuteron binding momentumand S-wave spin-triplet effective range, are required as experimental input. We furthermoreincluded final and initial state interactions consistently and showed explicitly that final stateinteractions are strongly suppressed at higher orders. We separated out the dipole contri-bution to quantify how much S-wave final state interactions contribute to this process. The /π EFT approach brings the advantage that calculations can be done largely analytical in thetwo-nucleon sector, and very few parameters are directly related to two-nucleon scatteringenter the calculation. At NNLO, we expect for an accuracy of approximately 5%.We compared our results for the energy shift δ A pol ,L in µ -d with recent calculations us-ing χ EFT and ZRA nuclear potentials. We found good agreement with these calculationsand confirm thereby our uncertainty estimate. We limited ourselves to NNLO because thedeuteron channel shape parameter, two-body current counterterm and S-D mixing opera-tor must be included at higher orders, making the analysis significantly more complicated.However, while much more involved, such a higher order calculation might also lead tointeresting results, specifically if the unknown counterterm could be used to describe uni-versal correlations between electromagnetic observables such as the deuteron radius and thenuclear polarizability correction. 21nother extension of this work is to calculate the TPE contribution to the Lamb shift inthe triton and Helium-3. Recent calculations show that the /π EFT is well-suited to describethe electromagnetic properties of the three-nucleon system [36, 37]. In the three-nucleonsystem, /π EFT loses some of its advantages. For example, the wave function has to becalculated numerically, too. We also expect significantly slower convergence and thereforelarger uncertainties as the binding momentum of the three-nucleon states is much largerthan the one of the deuteron. However, it is an important and interesting question by itselfhow well such details of the three-nucleon state can be described in /π EFT.
ACKNOWLEDGMENTS
We thank Richard Hill for discussions in the early stages of this project and Sonia Baccafor comments on this manuscript. This work has been funded by the National Science Foun-dation under Grant No. PHY-1555030, by the Office of Nuclear Physics, U.S. Departmentof Energy under Contract No. DE-AC05-00OR22725 and by the National Natural ScienceFoundation of China under Grant No. 11805078.
Appendix A: Relevant loop integrals in the power-divergence subtraction scheme1. Two-point loop integrals
We define the two-point loop integrals that are required in the calculation up throughNNLO as in Ref. [34], such that I PDS2 n ( E ) = i (cid:90) d q (2 π ) q n S ( q + E, q ) S ( − q , − q )= − m N π p n ( µ + ip ) , (A1)where p = √ m N E is the relative momentum within the two nucleon pair. For deuteronbound state we have p = iγ , with γ = √ m N B d denoting the deuteron binding momentum.22 . Loop integral J In diagram (c) of Fig. 3, the three-point loop integral J needed in Eq. (27) is given as J n = (cid:90) d l (2 π ) l n iS ( − B d + l , l ) iS ( − B d + l + ω, l + q ) iS ( − l , l ) , (A2)which has three nucleon propagators in the integrand. This integral can be solved in variousways. One way that relates it to the calculation of a quantum mechanical matrix elementof the charge operator is to solve it in coordinate space. We can reexpress the diagramas an integral over to coordinate space wave functions using Fourier transform of the twopropagators 1 − B d − k m N + i(cid:15) = − m N (cid:90) d re i k · r e − γr πr , E − ( k + q / m N + i(cid:15) = − m N (cid:90) d re i ( k + q / · r e ipr πr , (A3)where E = p /m N = ω − B d − q / (4 m N ). For J , placing the results of Eq. (A3) in Eq. (A2)gives J = − π (cid:16) m N π (cid:17) (cid:90) d re − γr e ipr sin( qr/ qr . (A4)Evaluating the integral in Eq. (A4), we find J = − m N πq (cid:34) cot − (cid:18) m N ω − q / qγ (cid:19) + i tanh − (cid:32) q (cid:112) m N ω − γ − q / m N ω (cid:33)(cid:35) , (A5)where we define the range of cot − as (0 , π ).In the evaluation of the matrix elements in this work we encounter the sum J + γ J .We can express this sum as J + γ J = m N (cid:90) d k (2 π ) m N ω − γ − q / − k + i(cid:15) . (A6)This integral can be solved using the PDS formula from KSW (cid:0) J + γ J (cid:1) PDS = − m N π (cid:32) µ + i (cid:114) m N ω − γ − q (cid:33) = m N I ( E ) , (A7)with m N E = m N ω − γ − q / p due to energy conservation.23 . Loop integral (cid:101) J At NLO, we encounter the loop diagram (cid:101) J in the loop that has NLO final state interac-tions. We define it as˜ J = (cid:90) d k (2 π ) ( k + q iS ( − B d + k , k ) iS ( − B d + k + ω, k + q ) iS ( − k , − k ) . (A8)Carrying out the countour integration gives˜ J = (cid:90) d k (2 π ) ( k + q i (cid:20) − B d − k m N + i(cid:15) (cid:21) − i (cid:20) − B d + ω − q m N − m N ( k + q + i(cid:15) (cid:21) − . (A9)We can now evaluate the relevant sum that involves (cid:101) J ˜ J − p J = (cid:90) d k (2 π ) (cid:20) ( k + q − p (cid:21) i (cid:20) − B d − k m N + i(cid:15) (cid:21) − × i (cid:20) − B d + ω − q m N − m N ( k + q + i(cid:15) (cid:21) − = m N π ( γ − µ ) , (A10)where we used p = m N ( ω − B d ) − q / J + γ J = m N I ( E ) = − m N π ( µ + ip ) . (A11) [1] I. Sick, “Precise radii of light nuclei from electron scattering,” inPrecision Physics of Simple Atoms and Molecules (Springer Berlin Heidelberg, Berlin,Heidelberg, 2008) pp. 57–77.[2] R. Pohl et al., Metrologia , L1 (2017).[3] R. Pohl et al., Science , 669 (2016).[4] R. Pohl, Nature , 213 (2010).[5] A. Antognini, Science , 417 (2013).[6] P. J. Mohr, D. B. Newell, and B. N. Taylor, Rev. Mod. Phys. , 035009 (2016).[7] A. Beyer et al., Science , 79 (2017).
8] N. Bezginov, T. Valdez, M. Horbatsch, A. Marsman, A. Vutha, and E. Hessels, Science ,1007 (2019).[9] H. Fleurbaey, S. Galtier, S. Thomas, M. Bonnaud, L. Julien, F. Biraben, F. Nez, M. Abgrall,and J. Guna, Phys. Rev. Lett. , 183001 (2018), arXiv:1801.08816 [physics.atom-ph].[10] W. Xiong et al., Nature , 147 (2019).[11] E. Borie, Annals of Physics , 733 (2012).[12] K. Pachucki, D. Leibfried, and T. W. H¨ansch, Phys. Rev. A , R1 (1993).[13] R. Rosenfelder, Nucl. Phys. A , 301 (1983).[14] Y. Lu and R. Rosenfelder, Physics Letters B , 7 (1993).[15] W. Leidemann and R. Rosenfelder, Phys. Rev. C , 427 (1995).[16] K. Pachucki, Phys. Rev. Lett. , 193007 (2011).[17] K. Pachucki and A. Wienczek, Phys. Rev. A , 040503 (2015).[18] O. Hernandez, C. Ji, S. Bacca, N. Nevo-Dinur, and N. Barnea, Physics Letters B , 344(2014).[19] O. Hernandez, A. Ekstrm, N. Nevo-Dinur, C. Ji, S. Bacca, and N. Barnea, Phys. Lett. B , 377 (2018).[20] O. J. Hernandez, C. Ji, S. Bacca, and N. Barnea, Phys. Rev. C , 064315 (2019),arXiv:1909.05717 [nucl-th].[21] J. L. Friar, Phys. Rev. C , 034003 (2013).[22] C. E. Carlson, M. Gorchtein, and M. Vanderhaeghen, Phys. Rev. A , 022504 (2014).[23] J. L. Friar, Phys. Rev. C , 1540 (1977).[24] C. Ji, N. Nevo-Dinur, S. Bacca, and N. Barnea, Phys. Rev. Lett. , 143402 (2013).[25] N. Nevo-Dinur, C. Ji, S. Bacca, and N. Barnea, Physics Letters B , 380 (2016).[26] C. E. Carlson, M. Gorchtein, and M. Vanderhaeghen, Phys. Rev. A , 012506 (2017).[27] C. Ji, S. Bacca, N. Barnea, O. J. Hernandez, and N. Nevo-Dinur, J. Phys. G , 093002(2018), arXiv:1806.03101 [nucl-th].[28] H.-W. Hammer, S. K¨onig, and U. van Kolck, Rev. Mod. Phys. , 025004 (2020),arXiv:1906.12122 [nucl-th].[29] H.-W. Hammer and L. Platter, Annu. Rev. Nucl. Part. Sci. , 207 (2010).[30] J. L. Friar, Annals of Physics , 151 (1979).[31] I. Sick, Phys. Rev. C , 064002 (2014).
32] N. Nevo-Dinur, O. J. Hernandez, S. Bacca, N. Barnea, C. Ji, S. Pastore, M. Piarulli, andR. B. Wiringa, Phys. Rev. C , 034004 (2019).[33] J.-W. Chen, G. Rupak, and M. J. Savage, Nucl. Phys. A653 , 386 (1999), arXiv:nucl-th/9902056 [nucl-th].[34] D. B. Kaplan, M. J. Savage, and M. B. Wise, Phys. Lett.
B424 , 390 (1998), arXiv:nucl-th/9801034.[35] D. R. Phillips, G. Rupak, and M. J. Savage, Phys. Lett.
B473 , 209 (2000), arXiv:nucl-th/9908054.[36] J. Vanasse, Phys. Rev. C , 024002 (2017), arXiv:1512.03805 [nucl-th].[37] J. Vanasse, Phys. Rev. C , 034003 (2018), arXiv:1706.02665 [nucl-th]., 034003 (2018), arXiv:1706.02665 [nucl-th].