Magnetic dipole operator from chiral effective field theory for many-body expansion methods

Many-body approaches for atomic nuclei generally rely on a basis expansion of the nuclear states, interactions, and current operators. In this work, we derive the representation of the magnetic dipole operator in plane-wave and harmonic-oscillator basis states, as needed for Faddeev calculations of few-body systems or many-body calculations within, e.g., the no-core shell model, the in-medium renormalization group, coupled-cluster theory, or the nuclear shell model. We focus in particular on the next-to-leading-order two-body contributions derived from chiral effective field theory. We provide detailed benchmarks and also comparisons with quantum Monte Carlo results for three-body systems. The derived operator matrix elements represent the basic input for studying magnetic properties of atomic nuclei based on chiral effective field theory.


I. INTRODUCTION
Calculating the electromagnetic structure of nuclei is a powerful tool to explore and test nuclear theory.The weak electromagnetic coupling compared to the strong interaction allows for a perturbative treatment of these processes, so that the nuclear structure content can be separated with great control.The electromagnetic interaction between the nucleus and external photons can in general be described by a currentcurrent interaction.While quantum electrodynamics (QED) describes the current of the external probe, nuclear theory deals with the nuclear current.To first approximation, the interaction between the photon and an atomic nucleus can be expressed in terms of the sum of photon interactions with all the individual nucleons.This approximation is equivalent to retaining only one-body contributions in the nuclear current, while all possible higher-body operators are neglected.Even though these leading terms provide the dominant contributions, higher-order contributions, especially from two-body operators are crucial for precise predictions of electromagnetic observables.
The modern approach to quantitatively understanding lowenergy nuclear physics in terms of ab initio calculations is based on effective field theory (EFT), most notably chiral EFT.It provides a systematic expansion of the strong interaction between nucleons as well as electroweak interactions with a direct connection to the fundamental theory of quantum chromodynamics (QCD) and its symmetries [1][2][3].A power-counting scheme orders the expansion terms according to decreasing importance in powers of (Q/Λ b ) ν , with Q the typical momentum scale governing processes in the nu-cleus, which is of the order of the pion mass m π , and Λ b the breakdown scale Λ b = 500 − 600 MeV.Leading order (LO) terms, i.e., ν = −2 for electromagnetic currents, include the dominant one-body contributions mentioned earlier, while next-to-leading order (NLO) and next-to-next-to-leading order (N 2 LO) terms, etc., add contributions of decreasing importance.The systematic expansion provides a way to improve calculations and to determine uncertainties arising from neglected higher orders [4,5].Furthermore, EFT provides a consistent derivation of nuclear forces and currents.To date, there have been several efforts to derive electromagnetic nuclear currents within the framework of chiral EFT.In Refs.[6][7][8] time-ordered perturbation theory was used to obtain current operator expressions up to next-to-next-to-next-to-leading order (N 3 LO) in the chiral expansion, while Refs.[9][10][11] used the method of unitary transformation.Both methods agree on the current operators at the order we employ in this work.However, at higher orders disagreements occur; for a detailed discussion see Ref. [11].
Calculating the electromagnetic structure of nuclei involves evaluating the electromagnetic nuclear current operator J µ = (ρ, j), with charge operator ρ and three-vector current operator j, between initial and final states of the nuclear system |i⟩ and | f ⟩.The Fourier transform of the current operator contains information about the charge and magnetization densities inside the nucleus.Because the nuclear states have a definite angular momentum, it is useful to decompose the nuclear current into its multipole components.For example, the current operator j can be expressed in terms of electric and magnetic multipole operators, the long-wavelength limits of which correspond to the electric and magnetic moment operators, where the magnetic dipole contribution is the focus of this work.With the magnetic dipole operator, one can calculate groundstate properties like the nuclear magnetic moment, defined by where J and M are the nuclear spin and its projection, respectively, and ξ represents all other quantum numbers relevant to describe the state.In addition, one can calculate magnetic arXiv:2308.00136v2[nucl-th] 4 Dec 2023 transitions between nuclear states.The probability of an initial state of the nucleus to emit or absorb a photon and transition to a final state is given by Fermi's golden rule [12]: where λ is the angular momentum of the photon with energy E γ .We use units with ℏ = c = ε 0 = 1, and the last term above is the transition strength Here, O mag λ represents the magnetic multipole operator in the long-wavelength limit (momentum transfer Q → 0) [13,14], with the reduced matrix element ⟨ξ f J f || O mag λ ||ξ i J i ⟩.In the case λ = 1 the multipole operator is the magnetic dipole operator O mag 1z = 3 4π µ z , and the corresponding transition is referred to as the dipole transition or M1 transition.
Studies of electromagnetic properties calculated with chiral EFT currents combined with nuclear states obtained from chiral EFT interactions or phenomenological potentials have to date been focused on few-nucleon systems and light nuclei.Deuteron and trinucleon electromagnetic form factors, radii, and moments have been studied up to N 3 LO in Refs.[15][16][17][18][19][20][21][22][23][24], with the most recent result for the charge and quadrupole form factors of the deuteron pushing the calculation to fifth order in the chiral expansion [25,26].In Ref. [27], magnetic moments and electromagnetic transitions of light nuclei up to A ≤ 9 have been calculated with a hybrid method, combining phenomenological wave functions with chiral magnetic dipole operators up to N 3 LO, based on quantum Monte Carlo methods to solve the many-body problem.More recently, the first full chiral EFT calculation of the ground-state magnetic moment and the lowest magnetic transition in 6 Li has been presented [28].All studies identified that current operator contributions beyond LO are important to improve agreement with experimental magnetic properties, with two-body contributions entering at NLO having the largest impact.
Higher-order corrections to the current operator are clearly necessary for improving the agreement with experimental results and can provide an explanation to long-standing discrepancies between theory and experiment.For example, the systematically smaller beta-decay rate in nuclei compared to free neutrons can, in part, be explained by the coupling of the weak force to two nucleons [29].In spite of this evidence for weak processes, the magnetic structure of heavier nuclei has so far only been studied without two-body currents (2BCs).Most ab initio many-body methods that calculate medium-mass nuclei are based upon basis expansion methods [30][31][32][33][34]. To perform computations, these frameworks require operators expanded in a computational basis that is commonly constructed based on harmonic-oscillator (HO) states.In this work, we provide partial-wave matrix elements for the LO one-body-current and NLO two-body current operators in a two-body momentumspace basis as well as partial-wave matrix elements for the corresponding LO and NLO magnetic dipole operators in HO bases.A straightforward implementation of these matrix elements can be used to calculate magnetic properties of medium-mass nuclei.We validate our expressions by comparing the trinucleon magnetic moments obtained from the the magnetic form factors with Faddeev calculations against the magnetic dipole operator used in Jacobi no-core shell-model (NCSM) calculations.Figure 1 displays the strategy of our trinucleon magnetic moments calculations.
This paper is organized as follows.In Sec.II we introduce the current operators that are employed and show how to obtain the magnetic dipole operators from them.Section III provides the expressions for the matrix elements of the various operators with respect to the different bases.The results and comparison of the trinucleon magnetic moments are presented in Sec.IV.Finally, we conclude in Sec.V.

II. NUCLEAR MAGNETIC MOMENTS
Nuclear magnetic moments can be calculated using two related methods.The first uses the magnetic form factor at zero momentum transfer.The form factor is the Fourier transform of the magnetization density of the nucleus and is obtained by calculating the expectation value of the nuclear current operator.The second method obtains the nuclear magnetic moment by directly evaluating the magnetic moment operator, which is the long-wavelength limit of the dipole term of the multipole expansion of the current operator.Below we specify the current and magnetic moment operators.

A. Magnetic form factor normalization
The one-body current operator at LO (∼ eQ −2 ) is given in momentum space by [21] where e is the elementary charge, m N the nucleon mass, k and k ′ are the initial and final nucleon momenta, K = (k + k ′ )/2, σ is the vector of Pauli spin matrices, and Q is the spatial part of the momentum transfer associated with the photon The functions e N Q 2 and µ N Q 2 are given by with G S /V E (G S /V M ) the isoscalar S and isovector V nucleon electric (magnetic) form factors, respectively.At zero momentum transfer, the form factors are known to be M (0) = 0.880 µ N , and G V M (0) = 4.706 µ N , where µ N = eℏ/2m proton is the nuclear magneton.For all form-factor calculations in this work we employ the nucleon parametrization FIG. 1. Schematic of the two methods used in this work to obtain the magnetic moments of the triton (t) and helion (h).The left part of the figure shows the steps for the current operator evaluated in a momentum-space basis to calculate the magnetic form factor, while the right part demonstrates the equivalent steps for the magnetic dipole operator evaluated between harmonic-oscillator basis states.derived by Ye et al. [35].This parametrization includes twophoton exchange corrections as well as information from new high-precision electron-nucleon scattering data, including uncertainties.
At NLO (∼ eQ −1 ), the leading 2BC operators enter.They are connected to the one-pion-exchange interaction and their momentum-space expressions are given by [21] with the axial coupling g A = 1.27, the pion decay constant F π = 92.3MeV, the averaged pion mass m π = 138.039MeV, the momentum transfers q i = k ′ i − k i , the two-body center of mass momenta 2 )/2, and the Pauli isospin matrices τ i , operating on nucleon i.Here momentum conservation implies P ′ = P + Q. Figure 2 shows the LO and NLO diagrams for the current operator.The first and second term in the square bracket in Eq. ( 7) correspond to diagram (b) and (c), respectively, which are commonly referred to as the "seagull" and "pion-in-flight terms."Higher orders of the current oper-ator have been derived see, e.g., Refs.[7,[9][10][11]21].Also, we note that the 2BC operator is not regularized in this work.In that sense, the operator is not fully consistent with employed nuclear interactions.Very recently, a consistent implementation was achieved by using semilocal coordinate-space regularization [36].
For the triton (t) and helion (h), the magnetic form factor is given by [27] where Q = |Q|, J = 1/2 and M J are the total three-body angular momentum and its projection, T = 1/2 and M T are the three-body isospin and its projection, |Ψ⟩ represents the three-body state, and we have suppressed the other quantum numbers of the triton or helion.As mentioned previously, the magnetic moments of triton and helion are given by the form factors at zero momentum transfer:

B. Magnetic moment operator
The magnetic moment operator is determined from the nuclear current operator in momentum space by [7] The current operator can be expanded as a sum of one-and many-body operators, resulting in a similar expansion for the magnetic moment operator where µ 1b,i is the single-nucleon contribution and µ 2b,i j the two-body part.The two-body magnetic dipole operator at NLO includes contributions from diagrams (b) and (c) in Fig. 2. The one-body magnetic dipole operator is given by [37] with the orbital angular momentum ℓ i .Because 2BC operators are translationally invariant with respect to the two-body center of mass R i j = (r i + r j )/2, the center-of-mass motion can be factored out as Accordingly, Eq. ( 10) splits into two parts, where one term depends only on the intrinsic coordinates and the other also on the center of mass.The intrinsic magnetic moment operator and is then given by [7] whereas the center-of-mass dependence is contained in the socalled "Sachs" term [38] µ Sachs 2b,i j = This division into two parts can be made for 2BC operators at any order.Summing the contributions of the seagull and the pion-inflight terms yields for the total NLO intrinsic operator (20) with f (r i j ) = 1+1/(m π r i j ) and the unit vector ri j of r i j = r i −r j .The result for the NLO Sachs term is given by where V 1π (r i j ) is the coordinate-space one-pion-exchange potential without isospin dependence: Here, S i j (r i j ) = 3(r i j • σ i )(r i j • σ j ) − σ i • σ j and h(r i j ) = 1 + 3/(m π r i j ) + 3/(m π r i j ) 2 .At NLO, the Sachs term is determined by the one-pion-exchange potential only.

III. OPERATOR MATRIX ELEMENTS
Matrix elements of operators expanded with respect to a chosen computational basis are essential components for calculating observables in few-and many-body methods.This section focuses on expanding the one-and two-body current and magnetic dipole operators into matrix elements with respect to a specific basis.First, we show the expansion of the current operators with respect to one-, two-, and three-body momentum-space Jacobi bases, where the three-body result is expressed in terms of the one-and two-body matrix elements.Next, we expand the magnetic dipole operator contributions with respect to one-, two-, and three-body relative HO bases.As the Sachs term explicitly depends on the two-body centerof-mass coordinate, embedding it into three-body Jacobi basis needs additional steps.We will show this embedding for the Sachs term in detail.
The single-particle partial-wave momentum-space basis states we employ are given by with the absolute value of the single-particle momentum k i = |k i |, orbital angular momentum l i , spin s i = 1 2 , total angular momentum j i and its projection m j i , and isospin t i = 1 2 along with its projection m t i for nucleon i.Relative two-body quantum numbers are denoted by capital letters and relative two-body momentum-space basis states are defined by with the relative momentum p = 1 2 (k 1 − k 2 ) of two nucleons, p = |p|, the relative orbital angular momentum L, twobody spin S , total angular momentum J and its projection M J , and total isospin T and its projection M T .The collective index α 2b defines the set of two-body quantum numbers α 2b = {L, S , J, T }.Three-body basis states are constructed by coupling a third nucleon to the two-body system and defining it relative to the center of mass of the nucleon pair with momentum p: Here is the second Jacobi momentum, whereas ℓ, s, j, and t are the corresponding spin, isospin and angular momentum quantum numbers [39].The total three-body angular momentum and isospin are denoted by J and T , respectively.Here again, the collective index α = {L, S , J, T, l, s, j, J, T } contains the partial-wave quantum numbers that define the state.
Partial-wave HO states are constructed in a similar manner, with the only difference being that the momentum is exchanged by the principle HO quantum number.Accordingly, the single-particle HO basis states are given by while the two-body basis states become and the three-body HO basis states are specified by Note that the Jacobi coordinates used in our NCSM calculations are not exactly the same as those used in the Faddeev calculations.All definitions are provided in Appendix A.

A. Partial-wave expanded current operator
Generally, the matrix elements of one-body and two-body current operators, as defined in Eqs. ( 4) and (7), can be expressed in the following form within the three-body partialwave basis defined in Eq. ( 25) [39,40]: with and the spinor spherical harmonics This factorized representation is very useful in practice as onebody operator can be represented in a natural way in terms of the quantity Q (Q, q, q ′ ), while for two-body operators the dynamics of the interaction with the external probe can be parametrized naturally via P (Q, p, p ′ ).

One-body operators
For the representation of one-body currents it is convenient to choose the coordinates such that the external probe interacts with the "last" particle, which for the 3N system is the nucleon with Jacobi momentum q resp.q ′ describing its motion relative to the subsystem of the other two nucleons.Specifically, we have: which implies for the Jacobi momenta: Consequently, the matrix elements of P(Q, p, p ′ ) are trivial since both particles with the relative momentum p = p ′ are just spectator particles: The quantity Q(Q, q, q ′ ) on the other hand contains the dynamics of the probe described by the one-body current and is given by: The evaluation is straightforward and can be performed numerically or partially analytically, depending on the specific form of the current operator.

Two-body operators
In the case of 2BC operators it is most convenient for the practical evaluation of the matrix elements to choose the two-body subsystem characterized by the Jacobi momenta p resp.p ′ to interact with the external probe, i.e., k ′ 3 = k 3 .We first express the single-particle momenta in terms of the Jacobi and center-of-mass momenta [39]: As a result of the interaction process, the center-of-mass momentum changes, P ′ 3N = P 3N + Q, and hence: while the Jacobi momentum p = (k 2 − k 1 )/2 and p ′ = (k ′ 2 − k ′ 1 )/2 are again unaffected by the interaction with the external probe.Overall, we obtain: The quantity Q is a current-independent function that in this case depends only on the kinematics specified by the momentum Q, while the two-body quantity P contains all the information about the current operator.

B. Harmonic-oscillator expanded magnetic dipole operator
In this section, we show matrix elements expressed in the HO basis, which enter the three-body Jacobi NCSM calculations.Without loss of generality we choose the external mo-mentum Q to be along the z direction.In the Jacobi NCSM, a wave function |Ψ 3b ⟩ is obtained through diagonalization of the Hamiltonian and expressed as a superposition of antisymmetrized HO basis states: Note that the antisymmetrized HO basis |i⟩ is not the same as the state defined in Eq. ( 28) and computed as with the coefficient of fractional parentage ⟨Nnα|i⟩ [41,42].
Through the antisymmetrization, it is clear that expectation values do not depend on the choice of the three-body Jacobi coordinate, i.e., one can choose the spectator particle.For example, the basis definitions ( 25) and ( 28) take spectator particle as the third one.Exploiting this, one can find for one-body operators.A similar expression can be found for two-body operators: The main tasks are to find expressions for and Here, we introduced doubly reduced matrix element with respect to spin and isospin, where κ is the isospin rank of the magnetic moment operator.This does not lose any information, and one can always restore normal matrix elements by means of the Wigner-Eckart theorem.

Harmonic-oscillator basis
To compute matrix elements within the HO basis, we first define the momentum-space representation of radial oscillator wave functions for a single-particle Rnℓ (k) = ⟨k ℓ|n ℓ⟩, defined by the overlap between momentum and HO eigenstates, given by the slightly modified definition from Ref. [43] Rnℓ with b ≡ 1/ √ m N ω the oscillator length in terms of the oscillator frequency ω and the nucleon mass m N , and L ℓ n (x) are generalized Laguerre polynomials.Similarly, coordinate-space radial wave functions R nℓ (r) = ⟨r ℓ|n ℓ⟩ are given by which are connected to the momentum-space functions through a Fourier-Bessel transform where the overlap is described by spherical Bessel functions ⟨r ℓ|k ℓ⟩ = √ 2/π j ℓ (kr).

One-body operator
Matrix elements of the one-body magnetic dipole operator defined in Eqs. ( 13)-( 16) are given by For the spin term µ κ spin , the required reduced matrix element for Eq. ( 47) is A similar expression can be found for the orbital contribution µ κ orb : ) Notice that there is a factor 2/3, coming from the transformation from single-particle to three-body Jacobi coordinates.Also, we have used that we can choose the orbital angular momentum of the three-body center-of-mass coordinate as zero since the intrinsic and center-of-mass motions do not couple.

Intrinsic magnetic dipole operator
After angular momentum recoupling, the reduced matrix element of the intrinsic magnetic dipole operator (20) with respect to a relative two-body HO basis (27) can be computed through with Also, the matrix element of the intrinsic two-body magnetic moment operator for Eq. ( 48) can be written as

Sachs operator
In the previous section, we described how to embed a twobody operator that only depends on the relative coordinate between two nucleons with respect to a relative two-body basis and a three-body Jacobi basis.Here, we consider a two-body operator depending on the two-nucleon center-of-mass, in addition to the relative coordinate.The Sachs contribution to the NLO magnetic dipole operator is of this type.
We start by constructing basis states that explicitly include the two-body center-of-mass motion.Such a basis is denoted by where the subscript 12 expresses relative quantities between nucleon 1 and 2, while NN indicates quantities related to the two-body center of mass.Coupling the relative and center-ofmass angular momenta generates the total angular momentum of the two-body system J rc with projection M J rc .A schematic representation of this basis is displayed in the right half of Fig. 3, which shows the momenta associated to a two-body system with respect to the origin by black dots labeled 1 and 2.
FIG. 3. Schematic of two different coordinate systems representing a three-nucleon system.The left part represents the Jacobi coordinate system characterized by (P 3N , p, q), while the right part shows the coordinate system described by (P NN , p, k 3 ).
The matrix element of the Sachs operator with respect to the states defined in Eq. ( 61) is given by In addition to the integration over r, the matrix elements of the Sachs operator are integrated over the center-or-mass coordi-nate R. Same as the intrinsic contribution, the isospin rank of the operator κ is 1.
Finding an expression for the matrix element of the Sachs operator in terms of three-body Jacobi basis states requires more work compared to the matrix elements we evaluated in previous sections.First, the two-body basis we defined in Eq. (61) has to be extended to include a third nucleon.This is achieved by coupling the total two-body angular momentum and isospin to the total angular momentum and isospin of the third nucleon to obtain three-body quantities.Second, this basis allows to evaluate the Sachs operator in terms of three-body states which include the two-body center-of-mass motion, so that the result in Eq. ( 62) can be used to express the matrix element.This expression, however, is unsuitable to calculate expectation values of the operator because the NCSM wave functions are expressed with three-body Jacobi states.Therefore, in a third step, we determine the overlap between the two different three-body bases.
To consider three particles in a basis which includes the two-body center-of-mass motion, we add a third nucleon represented by n 3 (ℓ 3 1 2 ) j 3 , which is defined with respect to the origin, and couple its angular momentum j 3 with the total twobody angular momentum J rc to the total three-body angular momentum J tot : The right part of Fig. 3 shows the coordinate system that corresponds to this basis state.
The matrix elements of the Sachs operator using the states defined in Eq. ( 64) can be represented in terms of reduced two-body matrix elements from Eq. ( 62), angular momentum and isospin coupling factors, and a third particle which is diagonal in all its quantum numbers: (65) This result shows that the majority of the work to calculate three-body matrix elements consists in determining the twobody matrix elements.
To calculate the overlap between the three-body basis states defined in Eq. ( 64) and the three-body Jacobi states from Eq. ( 28), the three-body center-of-mass motion has to be included.This is done by coupling the three-body angular momentum J with the center-of-mass orbital angular momentum L 3N to the total angular momentum J tot , so that the basis from Eq. ( 28) is extended to where the subscript 3N denotes quantities related to the threebody center of mass.This state corresponds to the coordinate representation in the left part of Fig. 3.The matrix elements of the Sachs operator in the three-body Jacobi basis can be obtained by carrying out the following transformation: tems.Also note that the Kronecker deltas indicate the relative two-body quantum numbers of the two different bases to be the same, which is expected as they essentially represent the same subsystem.A detailed derivation of this result can be found in Appendix B. Note that the object in the left-hand side in Eq. (67) takes the required form in Eq. ( 48) as we set N 3N = L 3N = 0.

IV. RESULTS
In this section, we examine the magnetic form factors and the magnetic moments of the trinucleons using different nuclear interactions based on chiral EFT.First, we present results for the magnetic form factors by evaluating the current operator in terms of the matrix elements presented in Section III A with corresponding partial-wave expanded wave functions, which are obtained by solving the three-body Faddeev equations.Magnetic moments are calculated as the zeromomentum-transfer limits of these form factors.We furthermore show results for the magnetic moments obtained from the expanded magnetic dipole operator expressions discussed in Sec.III B, based on NCSM wave functions.Finally, we compare Faddeev and NCSM calculations against each other to benchmark our results for the magnetic dipole operator.
We use the non-local chiral NN interactions by Entem, Machleidt, and Nosyk (EMN) [44] from LO to N 3 LO with cutoffs Λ = 420, 450, and 500 MeV.These are supplemented with 3N interactions at the same orders, with 3N low-energy constants (LECs) c D and c E determined by fits to the triton binding energy and nuclear matter saturation properties, and with a nonlocal three-body regulator with cutoff Λ 3N identical to the NN cutoff [45].A systematic study of the dependence of the magnetic observables on the three-body (LECs) c D and c E is, however, not pursued, because the magnetic form factors turn out to be nearly independent of the 3N interaction, and the magnetic moments even less so [46].In addition, we also consider the Entem and Machleidt (EM) [47] interaction at N 3 LO, with a cutoff Λ = 500 MeV.
The availability of the EMN potentials at each order of the expansion makes it possible to calculate theoretical uncertainty estimates of neglected higher-order terms based on the convergence pattern of observables.We use a Bayesian model as outlined in [5,48,49] to provide a statistical approach to calculate these uncertainties for the form-factor results.This method determines a posterior distribution which captures all the information about the neglected higher-order terms, from which degree-of-belief (DoB) intervals are calculated.For the evaluation, we employ a prior set C 0.25−10 with Λ b = 650 MeV, specified and publicly made available as a code in Ref. [48].A characteristic momentum scale of p = 2/3Q is used to calculate the 68% and 95% DoBs of the form factors. Overall, the characteristic scale for momentum Q transferred to a nucleus with mass number A is set by (A − 1)/A Q according to Ref. [50].
The trinucleon magnetic form factors and dipole moments are computed as and respectively.Here, |Ψ J,T M J , (F)⟩ and |Ψ J,T M J , (NCSM)⟩ are the Faddeev and NCSM wave functions, respectively, with total three-body spin and isospin specified by J = T = 1/2 and maximally projected total angular momentum states, i.e., M J = ±J.The isospin projection M T = −1/2 or M T = 1/2 determines whether the wave function represents a triton or a helion, respectively.
In order to perform Faddeev calculations, we truncate the basis by choosing a maximal value for the relative total twobody angular momentum J. Our calculations include partial waves up to J ≤ 6, which generates 42 distinct combinations of one-and two-body quantum numbers.This truncation proves to be sufficient to obtain converged results with respect to the basis states, so that any variation observed in the results are attributed to the interactions.

A. Magnetic form factor
We use the EMN interactions to calculate the magnetic form factors of the trinucleons with LO and NLO current operators.Figure 4 shows the triton (left column) and the helion (right column) magnetic form factors, in units of µ N , for calculations with the one-body current operator only (top row) and including 2BC corrections at NLO (bottom row) as a function of the momentum transfer Q, in units of fm −1 , and are compared to experimental results which are summarized by the hatched band [51].Results are given for cutoffs Λ = 420, 450, and 500 MeV at N 3 LO by the dashed-dotted, dashed, and solid lines, respectively, together with the 68% (light band) and 95% (dark band) DoBs for the 500 MeV result.These bands terminate around Q ∼ 4.6 fm −1 , because the Bayesian method has a limited range of validity.Lower orders of the interaction are used to calculate the truncation uncertainty, displayed by the colored bands, but not explicitly shown otherwise.
The magnetic form-factor results for 3 He obtained with only the one-body current disagree with experiment for all cutoffs over the entire momentum-transfer region and underestimate the data until the minimum, even when the truncation uncertainty is considered.For 3 H, the results disagree below 3.8 fm −1 when the truncation uncertainty is taken into account, yet fall within the 95% DoB interval at larger momentum transfer.Although the truncation uncertainty clearly increases as Q grows, note that the logarithmic scale overemphasizes the uncertainty at large Q in the figures.As expected, the cutoff dependence increases as well with momentum transfer, and the bands for the different cutoff values separate around Q ∼ 3 fm −1 ; nevertheless they remain within the 68% DoB intervals.At Q = 0, the results for Λ = 500 MeV, FIG. 4. Triton (left column) and helion (right column) magnetic form factors, in units of µ N , as a function of the momentum transfer, in units of fm −1 .The hatched band represents a parametrization of the elastic scattering data [51].The top row shows the result with the one-body current operator only, while the bottom row includes the 2BC contributions.The solid, dashed, and dashed-dotted lines represent the results for interactions at N 3 LO with cutoff Λ = 500, 450, and 420 MeV, respectively.Light and dark shaded bands represent the 95% and 68% DoBs.
F M, t (0) = 2.583 µ N and F M, h (0) = −1.767µ N , deviate from the experimental values for the triton (2.9789624659(59) µ N ) and helion (−2.127625307(25) µ N [52]) magnetic moments, and within the low-momentum transfer regime an approximate constant offset from experiment is found.In the following, we will examine the form factor normalization, i.e., the magnetic moment, in more detail.It is well known that higher-order two-body current operator corrections are sizable at Q = 0 [21,53], and thus they will impact the offset observed at low momentum transfers.The bottom row of Fig. 4 shows results with 2BC corrections included.Values at low momentum transfer are shifted up for both nuclei, but still disagree with experiment.Note that the DoB intervals of the one-body and two-body results do not overlap at low momentum transfers, which is a conse-quence of the inconsistent inclusion of current operators compared to the order of the interaction.At higher Q, the minimum is shifted to higher momentum transfers, so that the central value of the helion band is slightly too low and the central value of the triton band reproduces the minimum, confirming once again that 2BCs provide essential corrections to the onebody current operator.Because the chiral truncation uncertainty is expected to grow for increasing momentum transfers, the goal of exactly reproducing the minimum is too strong.We observe that within the truncation uncertainty the higher Q region fully overlaps with experimental data for both nuclei.The cutoff variation is slightly reduced with respect to the one-body result and falls well within the 68% DoBs.Stronger conclusions can only be made by including higher-order operators in the calculation.
FIG. 5. Same as Fig. 4, but compared to results from Ref. [21] and without DoB intervals.The green and yellow bands correspond to a variation of the cutoff scale from Λ = 500 MeV to Λ = 600 MeV.
In Fig. 5, we compare our results to calculations from Piarulli et al. [21], which are given by the yellow and green bands.The upper row displays the comparison of the onebody current operator (yellow band), while the bottom row shows results with 2BCs included (green band).Their results are obtained by calculating the expectation value of the operators with wave functions generated by the hyperspherical harmonics framework.The bands in this case represent the variation of the cutoff from Λ = 500 MeV to Λ = 600 MeV of the employed chiral interaction and therefore have a different interpretation compared to our bands.
Our one-body current results for 3 He with cutoffs Λ = 450 and 500 MeV fall within the yellow band over the entire momentum range, while our results for 3 H overestimate the minimum and the high momentum region.Considering that the operator is identical, the differences could only be explained by the use of different chiral interactions.This shows that apart from the corrections to the operator, also the interaction strongly influences the high momentum transfer region.
The green bands in the bottom row include 2BC operators up to N 3 LO, whereas our results only include the leading NLO contributions.At zero and low momentum transfers a small difference is observed, we will clarify its origin when discussing the magnetic moment below.For the minimum and the high momentum region, the observation is similar to the one-body comparison.Assuming that the truncation uncertainty for the results from Ref. [21] would be comparable to our findings, all results at high momentum transfers would agree and even be consistent with experiment.
In order to systematically study the low momentum transfer region, we display the trinucleon magnetic moments obtained from F M (0) in Fig. 6 as a function of increasing chiral order, for cutoffs 450 MeV (indicated by downward triangles and connected by dashed lines) and 500 MeV (indicated by cir-

LO
Comp. - 3 He FIG. 6. Triton and helion magnetic moments in units of µ N as a function of increasing order of the chiral expansion for the EMN NN + 3N interaction at cutoffs 450 MeV (downward triangles and dashed lines) and 500 MeV (circles and solid lines).Results from calculations with LO, NLO, NLO + NLO 2BC, N 2 LO + NLO 2BC, and N 3 LO + NLO 2BC interactions are shown in orange, yellow, green, blue, and red, respectively.Blue plus and cross symbols show the magnetic moment including NLO and N 2 LO current operator corrections, respectively, from Ref. [21].The black diamonds represent the experimental values for both nuclei [52].
cles and solid lines).They are compared to the experimental values of the triton and the helion, as well as to results from Piarulli et al. [21] (blue pluses and crosses), which include NLO 2BC corrections and the relativistic one-body correction to the magnetic moment at N 2 LO.The magnetic moments can be understood as one-nucleon hole with respect to 4 He, and the single-particle limit for triton (helion) of 2.793 (-1.913) µ N is reasonably close to the computed result.We observe that increasing the order of the chiral interaction used for the bound-state calculation has almost no effect on the magnetic moment of the trinucleons.However, adding NLO 2BC corrections (shown as "NLO + NLO 2BC" in the figure) changes the values by ∼10% and improves agreement with experiment.Our final result, labeled "N 3 LO + NLO 2BC," agrees well with both results from Piarulli et al., which implies that the relativistic correction to the one-body operator is very small.On the other hand, comparing the result to the experimental values suggests that important corrections to the operator are still missing to explain the remaining 5−7% discrepancy.Two-body corrections to the current operator at N 3 LO introduce new LECs that have to be fixed before predictions can be made [7,9].Different strategies exist for this procedure, and commonly a combination of observables is chosen which includes the isoscalar µ S and the isovector combination µ V of the trinucleon magnetic moments to constrain the new LECs [7,21]. 1 As a result, the experimental 2 The isoscalar (µ S ) and isovector (µ V ) combinations of the trinucleon mag-TABLE I. Triton and helion magnetic moments and their (cumulative) contributions, in µ N , from the form factor normalization, the magnetic dipole operator, as well as experimental values [52].The NCSM wave function is computed at N max = 40 and ℏω = 20 MeV.magnetic moments of the trinucleons are reproduced exactly if these higher-order corrections to the operator are taken into account.Therefore, tests of higher order 2BC require finite momentum transfer or nuclei beyond A = 3.

B. Test of magnetic dipole operator
In this section we present the trinucleon magnetic moments obtained from the magnetic dipole operator with the Jacobi NCSM, as discussed in Section III B. The three-body NCSM calculations are done with NuHamil code [54].We examine their convergence behavior and benchmark them to the magnetic moments obtained from the form factors in momentum space.Because our main goal is to benchmark the magnetic dipole operator matrix elements, we only consider the EM N 3 LO interaction with cutoff Λ = 500 MeV, without 3N interactions.
Figure 7 shows the convergence of the triton (left) and helion (right) magnetic moments as a function of N max for four different HO frequencies ℏω = 10, 20, 30, and 40 MeV.The top row displays the one-body magnetic dipole operator results relative to the one-body magnetic form factor normalization, while the bottom row shows the NLO-corrected result, represented by the label "µ HO, 1b + µ HO, 2b ", with the corresponding form-factor normalization in absolute terms.The two-body contribution consists of the intrinsic and Sachs terms, i.e., µ HO, 2b = µ NLO, intrinsic 2b + µ NLO, Sachs   2b   , where the majority of the correction is contributed by the intrinsic component.
Below we present numerical values and discuss the impact of both contributions in more detail.The results show a systematic convergence towards the desired results obtained from the form factor normalizations.Remarkably, very high values of N max are required in order to obtain converged results.This pattern follows the slow convergence of the three-body energy, which is a consequence of describing a loosely bound system in a HO basis.
Table I displays the magnetic moments obtained from both methods.The left column contains the results for the triton netic moments are defined by µ S ≡ µ t + µ h and µ V ≡ µ t − µ h .
FIG. 7. Convergence of the triton (left column) and helion (right column) magnetic moments as a function of N max for four different HO frequencies ℏω = 10, 20, 30, and 40 MeV obtained with the EM NN interaction with cutoff Λ = 500 MeV.The upper row shows the relative deviation, in %, of the magnetic moment calculated from the magnetic dipole operator compared to the form factor normalization, while the lower row shows the convergence of the magnetic moment including the NLO contributions, consisting of the intrinsic and Sachs components, towards the NLO form factor result (dashed line), in µ N .
magnetic moment µ t from the form-factor normalization and the magnetic dipole operator and the right column gives the same results for the helion µ h .The first row presents results from calculations with the one-body operator only ("LO"), while the second row shows results with the NLO 2BCs included ("NLO").Contributions from the latter to the magnetic moment are shown separately by the rows indicated with "intrinsic" and "Sachs."As noted in Section III B, this separation cannot be made for the form factor calculation, hence only the total values can be compared.At the bottom, the experimental results for both nuclei are given.The effect on the ground-state magnetic moment of the μNLO, intrinsic 2b operator accounts for the bulk of the correction and its influence amounts to around 10%, while the Sachs operator, which requires much more resources to calculate, has a minor effect of 0.5-1%.Total results from both methods agree with each other within ≪ 1%.
The excellent agreement between the magnetic moments of the triton and helion obtained from both methods gives strong confidence that the partial-wave decomposition of the LO and in particular the NLO magnetic dipole operator has been carried out correctly.After a transformation to single-particle coordinates (given in Appendix C), the matrix elements presented in this work can be used in many-body basis-expansion frameworks that are capable of calculating observables for heavier nuclei.For example, a recent NCSM calculation of the magnetic moment and a magnetic transition of 6 Li based on the developments presented in this work showed that the NLO corrections to the magnetic dipole operator are essential and improve the agreement with experiment [28].

V. SUMMARY AND CONCLUSIONS
In this paper, we studied the nuclear magnetic dipole operator obtained from chiral EFT current operators with a particular focus on the two-body NLO contribution.We discussed the general connection between the current operator and the magnetic dipole operator, and presented the coordinate space expressions for the NLO magnetic dipole operator.The magnetic dipole operator from 2BCs can be split into two terms: the intrinsic and Sachs terms, where the Sachs term depends explicitly on the two-body center-of-mass, in addition to the relative coordinate.We derived in detail the partial-wave decomposed matrix elements of the operators in the corresponding single-particle, two-body, and three-body bases.For the current operator we employed momentum-space basis states, to easily accommodate for the momentum dependence of the operator, while the magnetic dipole operator was evaluated with respect to HO basis states through Eqs. ( 58) and (62).
The decomposition of the latter has been checked by calculating the trinucleon magnetic moments and benchmarking them to form factor normalization results.
We provided results for the trinucleon magnetic form factors based on the EMN potentials, which allowed us to also estimate the uncertainty arising from truncating the chiral expansion.This uncertainty estimate, which by default grows for increasing momentum transfer, indicates that the precise reproduction of the minimum is not an important discriminator for EFTs.Our results agree well with values previously obtained in the literature using different chiral interactions and are consistent with experiment if the uncertainty in the EFT truncation is taken into account.In addition, our results for the magnetic dipole operator show a very good convergence and agree well with the form factor normalization results, demonstrating that the coordinate-space expression and the partialwave decomposition of the dipole operator are correct.
This work establishes a starting point for many-body expansion methods to incorporate the HO partial-wave decomposed NLO dipole matrix elements for calculations of electromagnetic observables.Such studies could validate already obtained results for A ≤ 9 systems [7,28] and will extend the ab initio analysis of NLO corrections to magnetic observables to medium-mass nuclei based on chiral EFT interactions and consistent current operators.As a next step, the partialwave decomposition of the magnetic dipole operator could be pushed to higher orders in the chiral expansion.This will test the chiral expansion and reduce the truncation uncertainty.
For our Faddeev calculations, the Jacobi momenta are defined as with the single-particle momenta k i , i = 1, 2, 3.The corresponding conjugate coordinates are given by with the single-particle coordinates r i , i = 1, 2, 3.For our NCSM calculations, we use a symmetric choice as the Talmi-Moshinsky bracket [55,56] is defined with them.The momenta and coordinates are then defined as and

(A4)
Appendix B: Three-body overlap for Sachs term evaluation Here we give a derivation of Eq. (68).To this end, we consider the following recoupling so that we can factor out the orbital part: (B1) To factorize the three-body center-of-mass part, one can use the coordinate transformation The above coordinate transformation ensures the following transformation using the Talmi-Moshinsky bracket [55,56], with the notation given in Ref. [57]: , n ′ ℓ ′ : λ|N NN L NN , n 3 ℓ 3 : λ⟩ d=2 .

Appendix C: Transformation to single-particle basis
For applications to medium-mass and heavier systems, current developments will need to be combined with basisexpansion methods such as coupled-cluster theory [30] or the in-medium-similarity renormalization group [31,33].Then, the matrix elements need to be expressed in terms of the single-particle coordinates rather than relative and centerof-mass coordinates.They are related by Talmi-Moshinsky transformations, already mentioned in Appendex B.
Our goal here is to show the transformation for Here the state |n 1 ℓ 1 j 1 , n 2 ℓ 2 j 2 : J tot ⟩ is the antisymmetrized product of {n 1 , ℓ 1 , j 1 } and {n 2 , ℓ 2 , j 2 } states coupling to the total angular momentum J tot in the proton-neutron (p-n) basis.
One can find the following transformation

FIG. 2 .
FIG.2.Diagrams for the LO (top row) and NLO (bottom row) contributions to the electromagnetic current operator, indicated by their scaling according to eQ ν .Solid lines represent nucleons, while dashed and wiggly lines represent pions and photons.Diagrams (b) and (c) are the leading 2BCs given by the seagull and pion-in-flight contribution, respectively.Note that the one-body charge operator is represented by diagram (a) too, but with order eQ −3 .