Zemach moments and radii of 2,3H and 3,4He
N. Nevo Dinur O. J. Hernandez S. Bacca N. Barnea C. Ji S. Pastore M. Piarulli R. B. Wiringa


































α



    TRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T 2A3, CanadaTRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T 2A3, CanadaDepartment of Physics and Astronomy, University of British Columbia, Vancouver, BC, V6T 1Z4, CanadaInstitut für Kernphysik and PRISMA Cluster of Excellence, Johannes-Gutenberg Universität, Mainz, DE-55128, GermanyInstitut für Kernphysik and PRISMA Cluster of Excellence, Johannes-Gutenberg Universität, Mainz, DE-55128, GermanyTRIUMF, 4004 Wesbrook Mall, Vancouver, BC, V6T 2A3, CanadaDepartment of Physics and Astronomy, University of Manitoba, Winnipeg, MB, R3T 2N2, CanadaRacah Institute of Physics, The Hebrew University, Jerusalem 91904, IsraelKey Laboratory of Quark and Lepton Physics (MOE) and Institute of Particle Physics, Central China Normal University, Wuhan 430079, China Physics Department, Washington University, St Louis, MO 63130, USA Physics Department, Washington University, St Louis, MO 63130, USAPhysics Division, Argonne National Laboratory, Argonne, IL 60439, USA
        We present  benchmark  calculations of Zemach moments and radii of ^2,3H and ^3,4He 
        using various few-body methods. Zemach moments are required to interpret muonic atom data 
        measured by the CREMA collaboration at the Paul Scherrer Institute. Conversely,  radii extracted
        from spectroscopic measurements can be compared with ab initio computations, posing stringent
        constraints on the nuclear model. For a given few-body method, 
        different numerical  procedures can be applied to compute these quantities. 
        A detailed analysis of the numerical uncertainties entering the  
        total theoretical error is presented. Uncertainties from the few-body method and the 
        calculational procedure are found to be smaller than the dependencies on the 
        dynamical modeling and the single nucleon inputs, which are found to be ≲ 2%. 
        When relativistic corrections and two-body currents are accounted for, the 
        calculated moments and radii are in very good agreement with the available experimental data. 
        21.45.+v, 21.10.Ky, 23.20.JsZemach moments and radii of ^2,3H and ^3,4He
            R.B. Wiringa
            December 30, 2023
        ============================================
        
        
        
        
        § INTRODUCTION
        
        
        Recent spectroscopic measurements on muonic atoms have enabled an extraction 
        of the charge radii of the proton <cit.> and 
        deuteron <cit.> with unprecedented precision, exposing 
        inconsistencies with measurements performed on electronic systems, 
        see, e.g., Refs. <cit.>. 
        The emergence of the so-called “proton-radius” and “deuteron-radius” 
        puzzles has attracted the attention of both the experimental and theoretical 
        physics communities. Regardless of the nature of these puzzles,
        it became clear that the precise determination of any nuclear charge radius 
        from spectroscopic measurements on its muonic atom/ion heavily relies 
        on an accurate knowledge of nuclear structure corrections to the 
        muonic spectrum <cit.>. 
        The CREMA collaboration has began investigating other light systems, 
        such as ^3,4He <cit.>, therefore,
        detailed studies on light nuclei are called for, and  demand a careful investigation 
        of all sources of uncertainty.
        
        In a hydrogen-like muonic atom or ion, 
        the energy difference 2S–2P, also called the Lamb shift (LS), is a sensitive 
        probe of the charge distribution of the nucleus (see, e.g., 
        Refs. <cit.> for reviews and 
        Ref. <cit.> and references therein for the most recent calculations). 
        In a Zα expansion up to 5^ th order, with α being the fine-structure
        constant and Z the proton number, this energy shift is related to the rms electric charge radius of the nucleus 
        R_E≡√(⟨R|_E^2|⟩) 
        by
        
            Δ E_ LS  =
            δ_ QED+
            
            𝒜_ QED·
            R_E^2 + δ_ TPE ,
        
        where 
        δ_ QED 
        and 𝒜_ QED 
        are independent of nuclear structure  and 
        are known to a very high accuracy from quantum electro-dynamics (QED).
        The precision of the radius extracted from 
        these measurements is  driven by the uncertainty in the  δ_ TPE term <cit.>. 
        The latter describes the two-photon 
        exchange (TPE) process where two virtual photons transfer energy and momentum to and from the nucleus.
        We note in passing that an analogous expression to Eq. (<ref>) 
        allows the extraction of the Zemach radius ⟨R|_Z |$⟩ (defined in Section <ref>) 
        from the measured hyperfine splitting (HFS) of
        a muonicnSstates <cit.>. Also in this case, accurate 
        nuclear structure calculations of the TPE contribution play a crucial 
        role <cit.>.
        
        
        Theδ_ TPEterm can be separated into elastic and inelastic 
        contributions. In the second case, the nucleus is excited to intermediate states.
        The elastic contributionδ_ Zemis  related to the third Zemach 
        moment of the electric form factorR^3_E_(2), 
        while the inelastic termδ_ polis related to the nuclear polarizability, 
        so thatδ_ TPE= δ_ Zem+δ_ pol.
        Notably, ab initio calculations reported in Ref. <cit.> currently provide the most 
        precise determinations ofδ_ polandδ_ Zemvalues and
        include nucleons' finite sizes but neglect 
        the contributions from two-body currents.
        One of the goals of this paper is to study the effect of two-body currents, 
        nuclear models, and different treatments of single-nucleon finite-sizes.
        
        The puzzles exposed by muonic laser spectroscopy have contributed 
        to the evolution of nuclear theory into a new era of precision,  where the 
        various sources of theoretical uncertainty need to be addressed adequately. 
        It is worth noting that while there has been considerable activity recently 
        devoted to the theoretical evaluation ofδ_ TPEin light muonic 
        atoms <cit.>, 
        the variety of few-body methods used in ab initio nuclear physics have yet to confront the computation 
        of nuclear Zemach moments and similar observables. A famous benchmark of different few-body methods for 
        computing the binding energy and radius of^4He dates back almost two decades <cit.> and 
        thus did not utilize state-of-the-art nuclear forces. More recently, other four-body and even five-body 
        benchmarks were performed, e.g., in Refs. <cit.>, which focused on 
        hadronic scattering rather than on electromagnetic observables.
        Filling this gap is among the goals of this work. To this end, we benchmark
        different ab initio methods on electromagnetic radii, Zemach moments, and 
        other ground-state observables for light nuclei in the mass range of2≤ A≤ 4, 
        which are of interest to the ongoing experimental efforts mentioned above.
        
        
        We focus on ground-state observables that can be readily calculated by 
        the few-body methods adopted here. We  neglectδ_ pol—which requires an additional computational development, 
        as described in Refs. <cit.>.
        
        
        In particular,  we solve theA=2problem using either the Numerov algorithm <cit.> 
        or the harmonic oscillator expansion used in Refs. <cit.>.
        ForA=3and 4, we use  Variational Monte Carlo (VMC) <cit.> and Green's function Monte Carlo 
        (GFMC) <cit.> methods, along with two different implementations of
        the  hyperspherical harmonics (HH) expansions, namely its momentum-space formulation 
        (HH-p) <cit.> and the effective interaction scheme in coordinate-space (EIHH) <cit.>. 
        These are all well-established methods and we do not provide further details on them here, 
        but rather refer the interested reader to the following 
        articles <cit.>.
        
        The paper is structured as follows. In Section <ref>, we define the 
        various electromagnetic observables under study and present the numerical procedures 
        implemented for their computations.
        In Section <ref>, we perform a benchmark in the impulse approximation (IA) using 
        wave functions from different few-body methods forA=2,3and4systems,  and compare the results 
        to experimental data. The agreement with data is reached by including  relativistic corrections 
        and two-body currents, whose contributions are studied only for theA=3nuclei. 
        
        Finally, we probe the sensitivity of our results to variations in both nuclear and nucleonic inputs,
        and in Section <ref> we draw our conclusions.
        
        
        
        
        § NUMERICAL PROCEDURES
        
        
        For a given few-body method that can provide nuclear wave functions, 
        different procedures can be used to calculate both Zemach and regular electromagnetic moments.
        We present a momentum-space, a coordinate-space, and a mixed (momentum & coordinate-spaces)
        formulation.  The latter exploits the respective advantages of the previous two methods.
        
        
        
        
         §.§ Definitions and momentum-space formulation
        
        The electric (E) and magnetic (M) form factors are 
        defined in momentum-space <cit.>
        as expectation values of the  ground state wave function 
        of the A-body nucleus. In particular, the deuteron electric 
        and magnetic form factors are defined, respectively, as <cit.>
            F_E(q^2)   =   1/3∑_M=± 1,0⟨ d;M|ρ(q ẑ)
            | d;M⟩ ,   
            
            F_M(q^2)   =   √(2) m_d/q  Im[ ⟨ d;1| j_y(q ẑ)
            |  d;0⟩ ] ,
        
        where| d;M⟩is the deuteron state with spin projectionJ_z=M,ρandj_ydenote, respectively, the charge
        operator andy-component of the current operator; the momentum transfer qis
        taken along thez-axis (the spin quantization axis), andm_dis the deuteron mass. 
        Form factors are normalized atq^2=0to1and(m_d/m_N) μ_d, 
        respectively, whereμ_dis the deuteron magnetic moment (in units of nucleon Bohr magnetonμ_N).
        
        The charge and magnetic form factors of the trinucleons are 
        derived from <cit.>
            F_E(q^2)    =   1/Z ⟨ +|ρ(q ẑ) | + ⟩ , 
            
            F_M(q^2)    =    -2  m_N/q  Im [  ⟨ -| j_y(q ẑ) |  + ⟩ ]  ,
        
        normalized to1andμ,
        whereμis the magnetic moment  of the three-body system (in units ofμ_N), and|±⟩represent either the^3He state or^3H state in spin projectionsJ_z=± 1/2.
        
        The charge and current operators are  expanded in many-body terms as
        
            ρ(q)   =   ∑_i^A ρ_i (q) + ∑_i<j^A ρ_ij(q) + … , 
             j(q)    =   ∑_i^A  j_i(q) + ∑_i<j^A  j_ij(q) +… .
        
        Calculations that retain only one-body terms in Eq. (<ref>) are typically called impulse approximation (IA)
        computations.
        In this paper, instead, we denote with IA those calculations that make use of only the leading-order (LO) one-body term
        in the chiral expansion of the electromagnetic operator <cit.>, basically
        excluding the relativistic one-body corrections.
        These operators are the standard charge and current one-body operators obtained from the
        non-relativistic reduction of the covariant nucleonic electromagnetic currents.
        In this work, we use two-body currents derived from a chiral effective field theory
        with pions and nucleons up to and including one-loop corrections <cit.>. 
        Note that, contributions from two-body terms enter at next-to-leading order (NLO) and at N^4LO
        in the chiral expansion of the current and charge operators, 
        respectively. Thus, two-body terms are expected to be sizable 
        in observables induced by the current operator and small in those induced by the charge operator. 
        
        The finite size of the nucleon is accounted for by including 
        the proton (p) and neutron (n) electric (E) and magnetic (M) form 
        factors,G^p/n_E/M(q^2). For example, the IA charge operator in the point-nucleon limit reads
        
            ρ( q)= ∑_i^Z exp(  q· r_i) ,
        
        where r_iis the coordinate of thei-th nucleon, and it becomes
        
            ρ( q)=∑_i^A exp( q· r_i) [ ( 1+τ^3_i/2) G^p_E(q^2)+ ( 1-τ^3_i/2) G^n_E(q^2)]  ,
        
         when the nucleonic form factors are included,
        withτ^3_ibeing the third isospin component of thei-th nucleon.G^p/n_E/M(q^2)are typically represented by parameterizations of 
        electron-scattering data, and here we will test the sensitivity of our results
        to different nucleonic inputs. 
        
        The nuclear electromagnetic form factorsF_EandF_Mcan be regarded as distributions, and thus one can define the corresponding 
        momenta at different orders in theq^2expansion.
        The 2^ ndand 4^ thelectric (magnetic) moments can be 
        derived from an expansion near momentum transferq^2→ 0of the charge 
        (magnetic) form factor as 
        
            F_x(q^2)= 1 -1/3!⟨ R_x^2 ⟩ q^2 +1/5!⟨ R_x^4 ⟩ q^4 + … ,
        
        where
        
            ⟨ R_x^2 ⟩   =    - 6  ∂ F_x(q^2)/∂ q^2|_q=0 , 
            ⟨ R_x^4 ⟩   =    60  ∂^2 F_x(q^2)/∂^2 q^2|_q=0 ,
        
        withx=E (M).
        Given the calculated nuclear form factorsF_x(q^2)at small values ofq^2, the⟨R|_x^2,4|$⟩ are then obtained from a quadratic fit 
        as indicated in Eq. (<ref>). From these, of course, follow
        estimates of, e.g.,  the rms  charge radius R_E=√(⟨R|_E^2|⟩), 
        which is measured by Lamb shift experiments in muonic atoms using Eq. (<ref>).
        
        The elastic component of δ_ TPE in Eq. (<ref>), 
        namely δ_ Zem, is directly proportional to the third Zemach 
        moment <cit.>, defined as
        
            R^3_E_(2)  = 48/π∫_0^∞dq/q^4[ F_E^2(q^2) -1 + q^2 ⟨R|_E^2|⟩/3] .
        
        
        The Zemach radius (traditionally called the first Zemach moment) 
        is a quantity of mixed electric and magnetic nature, 
        defined as
        
            ⟨R|_Z |=⟩ -4/π∫_0^∞dq/q^2[ F_E(q^2) F_M(q^2)/F_M(0) -1  ]   .
        
        It was first developed by Zemach in Ref. <cit.> in the context of hyperfine 
        splitting in hydrogen S-states, where the leading correction due to the proton's finite 
        size was shown to be proportional to  R_Z.
        Consequently, R_Z of a spin-half nucleus can be experimentally determined, e.g., 
        from the hyperfine splitting in its muonic hydrogen-like atom/ion, with precision that 
        could rival determinations from electron scattering <cit.>.
        
        
        R_E=R_E^2^1/2, R^3_E_(2), R_E^4, R_Z, R_M=R_M^2^1/2
        and μ are the observables we study in this paper. Since they are all essentially moments
        of electromagnetic distributions, we refer to them cumulatively as “electromagnetic moments”.
        
        We would like to comment on the q integration that enter in the above definitions. Clearly for a certain large value of q, denoted as Q_max, the form factors F_E/M(q^2) 
        are too small to contribute to the integrals in Eqs. (<ref>) and (<ref>). Therefore from q=Q_max up to 
        q=∞ the tail of the integrand is given  by the analytical expression in Eqs. (<ref>) and (<ref>) where F_E/M(q^2) are set to 0. 
        On the other hand, the integrands of the above equations are numerically unstable near q=0.
        Therefore at 0⩽ q ⩽ Q_min, where Q_min is a small value, they are replaced with their low-q^2 approximations 
        
            lim_q^2→ 01/q^4[ F_E^2(q^2)-1+q^2  R_E^2/3]
            = R^4_E/36+R^4_E/60 , 
                  lim_q^2→ 01/q^2[ F_E(q^2) F_M(q^2) -1  ]= - ⟨ R_E^2⟩/6-⟨ R_M^2⟩/6    .
        
        
        Using Eqs. (<ref>)–(<ref>) to calculate electromagnetic moments 
        is hereafter referred to as the momentum-space numerical procedure and denoted with “q-space”.
        
        
        
         §.§ Coordinate-space formulation
        
        The rms charge radius, as well as other even moments, can be readily obtained from  point-nucleon computations in coordinate-space.
        In the IA and in the non-relativistic limit, the  2^ nd and 4^ th moments of the electric charge distribution can be obtained  as
        
            ⟨ R^2_E⟩    =   ⟨ R^2_p ⟩  +  r_p^2+ N/Z r^2_n, 
            R^4_E    =   ⟨ R^4_p ⟩ + r^4_p + N/Z r^4_n
            + 10/3( r^2_p⟨ R^2_p⟩ + N/Z r^2_n⟨ R^2_n⟩) ,
        
        where, the point-proton mean-square radius is calculated as an expectation value on the ground-state wave-function
        
            ⟨ R^2_p ⟩ = ⟨Ψ_0 |1/Z∑_i^Z r^2_i | Ψ_0 ⟩ .
        
        Analogous
        expressions exist for the point-neutron radius ⟨ R^2_n ⟩ and for ⟨ R^4_p ⟩. 
        We perform our benchmark calculations with the
        Kelly parameterization of the nucleon form factors <cit.>.
        Accordingly, the 2^ nd and 4^ th moments of the intrinsic nucleon electric form factors are taken to be 
        r^2_p=0.744(7) fm^2, r^2_n= -0.112(3) fm^2, 
        r^4_p =1.6(1) fm^4, and r^4_n= -0.33(2) fm^4.
        
        
        Using Eqs. (<ref>) to calculate 
        ⟨ R^2_E⟩ and R^4_E will be 
        referred to as the coordinate-space numerical procedure and denoted with “r-space”. 
        
        
        
         §.§ Mixed momentum & coordinate-space formulation
        
        Given F_E(q^2), one can obtain the charge density in coordinate-space, in the non-relativistic limit, as its Fourier transform
        
            ρ_E ( r) = ∫ d^3 q  F_E(q^2)  e^-i q·r .
        
        The n-th electric Zemach moment is defined as
        
            ⟨R|^n_E|_⟩(2) = ∫ d^3 r∫ d^3r' ρ_E (r') ρ_E (|r'-r|) r^n .
        
        By inserting Eq. (<ref>) into Eq. (<ref>)
        one obtains
        
            ⟨R|^n_E|_⟩(2) = 
            ∫_0^∞ dr r^n+1[
            2/π∫_0^∞ dq   q F^2_E(q^2) sin(qr)] ,
        
        which contains integrals on both q and r. 
        In the above we used explicitly only the contribution from the “spherical” 
        part of the charge distribution, which is an approximation for the deuteron, but is  exact for A= 3 and 4. 
        
        This algorithm was found to be very robust, 
        and does not suffer the numerical uncertainty associated with the regularization of Eq. (<ref>) at q=0.
        
        Obviously, F^2_E(q^2) in Eq. (<ref>) may be replaced with either 
        F_E(q^2)· F_M(q^2) or F^2_M(q^2), leading to the calculation of other moments, 
        such as, e.g., R_Z of Eq. (<ref>).
        Additionally, there exist relations between various Zemach and regular moments, 
        e.g.,
        
            2 ⟨R|^2_E|   ⟩=   R^2_E_(2) , 
            
            2 R^4_E   =   R^4_E_(2) - 10/3⟨R|^2_E|^⟩2 ,
        
        which enable the consistent calculation of essentially all the regular and Zemach moments, and particularly all the observables targeted here, via this procedure.
        
        Using Eqs. (<ref>)–(<ref>) to
         calculate electromagnetic moments 
        will be referred to as the mixed-space numerical procedure, which we label with “qr-space”. 
        
        
        
         §.§ Numerical Procedures: Comparison
        
        We apply the numerical procedures detailed above to study the ^3He electric moments. 
        In particular, the r- and qr-space procedures are used in combination with EIHH
        few-body computational method, while the q-space procedure is implemented within the 
        HH-p method.
