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

R. Seutin<EMAIL_ADDRESS>Max-Planck-Institut für
Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany  Technische
Universität Darmstadt, Department of Physics, 64289 Darmstadt, Germany
ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für Schwerionenforschung
GmbH, 64291 Darmstadt, Germany  O. J. Hernandez  Institut für Kernphysik and
PRISMA Cluster of Excellence, Johannes Gutenberg-Universität Mainz, 55128
Mainz, Germany  T. Miyagi<EMAIL_ADDRESS>Technische Universität Darmstadt, Department of Physics, 64289 Darmstadt,
Germany  ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für
Schwerionenforschung GmbH, 64291 Darmstadt, Germany  Max-Planck-Institut für
Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany  S. Bacca
<EMAIL_ADDRESS>Institut für Kernphysik and PRISMA Cluster of Excellence,
Johannes Gutenberg-Universität Mainz, 55128 Mainz, Germany  Helmholtz
Institute Mainz, GSI Helmholtzzentrum für Schwerionenforschung GmbH, 64289
Darmstadt, Germany  K. Hebeler<EMAIL_ADDRESS>Technische
Universität Darmstadt, Department of Physics, 64289 Darmstadt, Germany
ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für Schwerionenforschung
GmbH, 64291 Darmstadt, Germany  Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  S. König<EMAIL_ADDRESS>North
Carolina State University, Department of Physics, Raleigh, NC 27695, USA
Technische Universität Darmstadt, Department of Physics, 64289 Darmstadt,
Germany  ExtreMe Matter Institute EMMI, GSI Helmholtzzentrum für
Schwerionenforschung GmbH, 64291 Darmstadt, Germany  A. Schwenk
<EMAIL_ADDRESS>Technische Universität Darmstadt, Department of
Physics, 64289 Darmstadt, Germany  ExtreMe Matter Institute EMMI, GSI
Helmholtzzentrum für Schwerionenforschung GmbH, 64291 Darmstadt, Germany  Max-
Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117 Heidelberg, Germany

###### Abstract

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 current-current 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 low-energy 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/\Lambda_{\text{b}})^{\nu}$, with $Q$ the typical momentum scale governing
processes in the nucleus, which is of the order of the pion mass $m_{\pi}$,
and $\Lambda_{\text{b}}$ the breakdown scale $\Lambda_{\rm b}=500-600$ MeV.
Leading order (LO) terms, i.e., $\nu=0$ for nuclear forces and $\nu=-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 (N2LO) 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 (N3LO) 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_{\mu}=(\rho,\mathbf{j})$, with
charge operator $\rho$ and three-vector current operator $\mathbf{j}$, between
initial and final states of the nuclear system $|i\rangle$ and $|f\rangle$.
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
$\mathbf{j}$ can be expressed in terms of magnetic multipole operators, the
long-wavelength limits of which correspond to the magnetic moment operators,
which are dominated by the magnetic dipole contribution. With the magnetic
dipole operator, one can calculate ground-state properties like the nuclear
magnetic moment, defined by

$\mu\equiv\matrixelement{\xi JM=J\,}{\,\mu_{z}\,}{\,\xi JM=J},$ (1)

where $J$ and $M$ are the nuclear spin and its projection, respectively, and
$\xi$ represents all other quantum numbers relevant to describe the state. In
addition, one can calculate magnetic 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]:

$\Gamma_{\gamma,i\rightarrow
f}=\frac{2\,(\lambda+1)}{\lambda\big{[}(2\lambda+1)!!\big{]}^{2}}\,E_{\gamma}^{2\lambda+1}B(M\lambda;\xi_{i}J_{i}\rightarrow\xi_{f}J_{f}),$
(2)

where $\lambda$ is the angular momentum of the photon with energy
$E_{\gamma}$. We use units with $\hbar=c=\varepsilon_{0}=1$, and the last term
above is the transition strength

$B(M\lambda;\xi J_{i}\rightarrow\xi
J_{f})\equiv\frac{1}{2J_{i}+1}\big{|}\langle\xi_{f}J_{f}||\,O^{\text{mag}}_{\lambda}\,||\xi_{i}J_{i}\rangle\big{|}^{2}.$
(3)

Here, $O^{\text{mag}}_{\lambda}$ represents the magnetic multipole operator in
the long-wavelength limit (momentum transfer $Q\to 0$) [13, 14], with the
reduced matrix element
$\langle\xi_{f}J_{f}||\,O^{\text{mag}}_{\lambda}\,||\xi_{i}J_{i}\rangle$. In
the case $\lambda=1$ the multipole operator is the magnetic dipole operator
$O^{\text{mag}}_{1z}=\sqrt{\frac{3}{4\pi}}\mu_{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
phenomenlogical 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 N3LO in Refs. [15, 16, 17, 18, 19, 20, 21,
22, 21, 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\leq 9$ have been calculated with a hybrid
method, combining phenomenological wave functions with chiral magnetic dipole
operators up to N3LO, 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 6Li 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 momentum-space 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.

Figure 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.

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.

###  II.1 Magnetic form factor normalization

The one-body current operator at LO ($\sim eQ^{-2}$) is given in momentum
space by [21]

$\displaystyle\mathbf{j}^{(-2)}(\mathbf{Q})=\frac{e}{2m_{N}}\bigg{(}2\,e_{N}\big{(}Q^{2}\big{)}\,\mathbf{K}+\mathrm{i}\,\mu_{N}\big{(}Q^{2}\big{)}\,\bm{\sigma}\times\mathbf{Q}\bigg{)},$
(4)

where $e$ is the elementary charge, $m_{N}$ the nucleon mass, $\mathbf{k}$ and
$\mathbf{k}^{\prime}$ are the initial and final nucleon momenta,
$\mathbf{K}=(\mathbf{k}+\mathbf{k}^{\prime})/2$, $\bm{\sigma}$ is the vector
of Pauli spin matrices, and $\mathbf{Q}$ is the spatial part of the momentum
transfer associated with the photon $Q^{2}$ is $Q^{\mu}Q_{\mu}$. Momentum
conservation requires that the relation
$\mathbf{k}^{\prime}=\mathbf{k}+\mathbf{Q}$ holds. The functions
$e_{N}\big{(}Q^{2}\big{)}$ and $\mu_{N}\big{(}Q^{2}\big{)}$ are given by

$\displaystyle e_{N}\big{(}Q^{2}\big{)}$
$\displaystyle=\frac{G_{\textnormal{E}}^{S}\big{(}Q^{2}\big{)}+G_{\textnormal{E}}^{V}\big{(}Q^{2}\big{)}\tau_{z}}{2},$
(5) $\displaystyle\mu_{N}\big{(}Q^{2}\big{)}$
$\displaystyle=\frac{G_{\textnormal{M}}^{S}\big{(}Q^{2}\big{)}+G_{\textnormal{M}}^{V}\big{(}Q^{2}\big{)}\tau_{z}}{2},$
(6)

with $G_{\textnormal{E}}^{S/V}$ ($G_{\textnormal{M}}^{S/V}$) the isoscalar $S$
and isovector $V$ nucleon electric (magnetic) form factors, respectively. At
zero momentum transfer, the form factors are known to be
$G_{\textnormal{E}}^{S}(0)=G_{\textnormal{E}}^{V}(0)=1$,
$G_{\textnormal{M}}^{S}(0)=0.880$ $\mu_{N}$, and
$G_{\textnormal{M}}^{V}(0)=4.706$ $\mu_{N}$, where $\mu_{N}=e\hbar/2m_{\rm
proton}$ is the nuclear magneton. For all form-factor calculations in this
work we employ the nucleon parametrization derived by Ye et al. [35]. This
parametrization includes two-photon exchange corrections as well as
information from new high-precision electron-nucleon scattering data,
including uncertainties.

At NLO ($\sim 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]

$\displaystyle\mathbf{j}^{(-1)}(\mathbf{Q})=$
$\displaystyle-\mathrm{i}\,e\frac{g_{A}^{2}}{4F_{\pi}^{2}}G^{V}_{\textnormal{E}}(Q^{2})(\bm{\tau}_{1}\times\bm{\tau}_{2})_{z}\bigg{[}\bm{\sigma}_{1}-\mathbf{q}_{1}\frac{\bm{\sigma}_{1}\cdot\mathbf{q}_{1}}{q_{1}^{2}+m_{\pi}^{2}}\bigg{]}$
$\displaystyle\times\frac{\bm{\sigma}_{2}\cdot\mathbf{q}_{2}}{q_{2}^{2}+m_{\pi}^{2}}+1\leftrightharpoons
2,$ (7)

with the axial coupling $g_{A}=1.27$, the pion decay constant $F_{\pi}=92.3$
MeV, the averaged pion mass $m_{\pi}=138.039$ MeV, the momentum transfers
$\mathbf{q}_{i}=\mathbf{k}_{i}^{\prime}-\mathbf{k}_{i}$, the two-body center
of mass momenta
$\mathbf{P}^{(\prime)}=(\mathbf{k}_{1}^{(\prime)}+\mathbf{k}_{2}^{(\prime)})/2$,
and the Pauli isospin matrices $\bm{\tau}_{i}$, operating on nucleon $i$. Here
momentum conservation implies $\mathbf{P}^{\prime}=\mathbf{P}+\mathbf{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 operator have been derived see,
e.g., Refs. [7, 21, 9, 10, 11]. 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].

Figure 2: Diagrams for the LO (top row) and NLO (bottom row) contributions to
the electromagnetic current operator, indicated by their scaling according to
$e\mathcal{Q}^{\nu}$. 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
$e\mathcal{Q}^{-3}$.

For the triton (t) and helion (h), the magnetic form factor is given by [27]

$\displaystyle F_{\textnormal{M}}(Q)$
$\displaystyle=-\frac{2m_{N}}{Q}\langle\Psi^{\mathcal{J},\mathcal{T}}_{\mathcal{M}^{\prime}_{\mathcal{J}}=+\mathcal{J},\mathcal{M}_{\mathcal{T}}}|\,\mathrm{j}_{+}(\mathbf{Q})\,|\Psi^{\mathcal{J},\mathcal{T}}_{\mathcal{M}_{\mathcal{J}}=-\mathcal{J},\mathcal{M}_{\mathcal{T}}}\rangle,$
(8)

where $Q=|\mathbf{Q}|$, $\mathcal{J}=1/2$ and $\mathcal{M}_{\mathcal{J}}$ are
the total three-body angular momentum and its projection, $\mathcal{T}=1/2$
and $\mathcal{M}_{\mathcal{T}}$ are the three-body isospin and its projection,
$|\Psi\rangle$ 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:

$\mu_{\textnormal{t}/\textnormal{h}}=F_{\textnormal{M}}(0).$ (9)

###  II.2 Magnetic moment operator

The magnetic moment operator is determined from the nuclear current operator
in momentum space by [7]

$\bm{\mu}=-\frac{\mathrm{i}}{2}\lim_{\mathbf{Q}\rightarrow
0}\,\bm{\nabla}_{\mathbf{Q}}\times\mathbf{j}(\mathbf{Q}).$ (10)

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

$\bm{\mu}=\sum_{i}^{A}\bm{\mu}_{\textnormal{1b},i}+\sum_{i<j}^{A}\bm{\mu}_{\textnormal{2b},ij}+\ldots,$
(11)

where $\bm{\mu}_{\textnormal{1b},i}$ is the single-nucleon contribution and
$\bm{\mu}_{\textnormal{2b},ij}$ 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]

$\displaystyle\bm{\mu}_{\textnormal{1b},i}=$
$\displaystyle\sum_{\kappa=0}^{1}(\bm{\mu}^{\kappa}_{\textnormal{spin},i}+\bm{\mu}^{\kappa}_{\textnormal{orb},i}),$
(12) $\displaystyle\bm{\mu}^{0}_{\textnormal{spin},i}$
$\displaystyle=G^{S}_{M}(0)\bm{\sigma}_{i},$ (13)
$\displaystyle\bm{\mu}^{1}_{\textnormal{spin},i}$
$\displaystyle=G^{V}_{M}(0)\tau_{i,z}\bm{\sigma}_{i},$ (14)
$\displaystyle\bm{\mu}^{0}_{\textnormal{orb},i}$
$\displaystyle=\frac{\mu_{N}}{2}\bm{\ell}_{i},$ (15)
$\displaystyle\bm{\mu}^{1}_{\textnormal{orb},i}$
$\displaystyle=\frac{\mu_{N}}{2}\tau_{i,z}\bm{\ell}_{i},$ (16)

with the orbital angular momentum $\bm{\ell}_{i}$.

Because 2BC operators are translationally invariant with respect to the two-
body center of mass $\mathbf{R}_{ij}=(\mathbf{r}_{i}+\mathbf{r}_{j})/2$, the
center-of-mass motion can be factored out as

$\mathbf{j}_{{\rm
2b},ij}(\mathbf{Q},\mathbf{R}_{ij})=e^{\mathrm{i}\mathbf{Q}\cdot\mathbf{R}_{ij}}\mathbf{j}_{{\rm
2b},ij}(\mathbf{Q}).$ (17)

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]

$\bm{\mu}_{{\rm
2b},ij}^{\textnormal{int}}=-\frac{\mathrm{i}}{2}\lim_{\mathbf{Q}\rightarrow
0}\,\bm{\nabla}_{\mathbf{Q}}\times\mathbf{j}_{{\rm 2b},ij}(\mathbf{Q}),$ (18)

whereas the center-of-mass dependence is contained in the so-called “Sachs”
term [38]

$\bm{\mu}_{{\rm
2b},ij}^{\textnormal{Sachs}}=\frac{1}{2}\,\mathbf{R}_{ij}\times\mathbf{j}_{{\rm
2b},ij}(\mathbf{Q}).$ (19)

This division into two parts can be made for 2BC operators at any order.

Summing the contributions of the seagull and the pion-in-flight terms yields
for the total NLO intrinsic operator

$\displaystyle\bm{\mu}^{\textnormal{NLO, int}}_{{\rm
2b},ij}(\mathbf{r}_{ij})=$ $\displaystyle-\frac{g_{A}^{2}m_{\pi}}{32\pi
F_{\pi}^{2}}(\bm{\tau}_{i}\times\bm{\tau}_{j})_{z}\,\bigg{\\{}f(r_{ij})\,\big{[}(\bm{\sigma}_{i}\times\bm{\sigma}_{j})\cdot\hat{\mathbf{r}}_{ij}\big{]}\,\hat{\mathbf{r}}_{ij}$
$\displaystyle\mathopen{\hphantom{-\frac{g_{A}^{2}m_{\pi}}{8\pi
F_{\pi}^{2}}(\bm{\tau}_{1}\times\bm{\tau}_{2})_{z}\,\bigg{\\{}f}}-(\bm{\sigma}_{i}\times\bm{\sigma}_{j})\bigg{\\}}e^{-m_{\pi}r_{ij}},$
(20)

with $f(r_{ij})=1+1/(m_{\pi}r_{ij})$ and the unit vector
$\hat{\mathbf{r}}_{ij}$ of $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$.
The result for the NLO Sachs term is given by

$\bm{\mu}^{\textnormal{NLO, Sachs}}_{{\rm
2b},ij}(\mathbf{r}_{ij},\mathbf{R}_{ij})=-\frac{1}{2}(\bm{\tau}_{i}\times\bm{\tau}_{j})_{z}\,V_{1\pi}(\mathbf{r}_{ij})\,\mathbf{R}_{ij}\times\mathbf{r}_{ij},$
(21)

where $V_{1\pi}(\mathbf{r}_{ij})$ is the coordinate-space one-pion-exchange
potential without isospin dependence:

$\displaystyle V_{1\pi}(\mathbf{r}_{ij})=$
$\displaystyle\frac{m_{\pi}^{2}}{12\pi}\frac{g_{A}^{2}}{4F_{\pi}^{2}}\bigg{\\{}\bigg{[}S_{ij}(\hat{\mathbf{r}}_{ij})\,h(r_{ij})+\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}\bigg{]}\frac{e^{-m_{\pi}r_{ij}}}{r_{ij}}$
$\displaystyle\mathopen{\hphantom{\frac{m_{\pi}^{2}}{12\pi}\frac{g_{A}^{2}}{F_{\pi}^{2}}(\bm{\tau}_{1}\cdot\bm{\tau}_{2})\bigg{\\{}\bigg{[}}}-\frac{4\pi}{3}\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}\,\delta(\mathbf{r}_{ij})\bigg{\\}}.$
(22)

Here,
$S_{ij}(\hat{\mathbf{r}}_{ij})=3(\hat{\mathbf{r}}_{ij}\cdot\bm{\sigma}_{i})(\hat{\mathbf{r}}_{ij}\cdot\bm{\sigma}_{j})-\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}$
and $h(r_{ij})=1+3/(m_{\pi}r_{ij})+3/(m_{\pi}r_{ij})^{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 center-of-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

$\displaystyle\ket{k_{i}\,(\ell_{i}s_{i})j_{i}m_{j_{i}}\,t_{i}m_{t_{i}}},$
(23)

with the absolute value of the single-particle momentum
$k_{i}=|\mathbf{k}_{i}|$, orbital angular momentum $l_{i}$, spin
$s_{i}=\tfrac{1}{2}$, total angular momentum $j_{i}$ and its projection
$m_{j_{i}}$, and isospin $t_{i}=\tfrac{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

$\displaystyle\ket{p\,\alpha_{\text{2b}}}\equiv\ket{p\,(LS)JM_{J}\,TM_{T}},$
(24)

with the relative momentum
$\mathbf{p}=\frac{1}{2}(\mathbf{k}_{1}-\mathbf{k}_{2})$ of two nucleons,
$p=|\mathbf{p}|$, the relative orbital angular momentum $L$, two-body spin
$S$, total angular momentum $J$ and its projection $M_{J}$, and total isospin
$T$ and its projection $M_{T}$. The collective index $\alpha_{\text{2b}}$
defines the set of two-body quantum numbers
$\alpha_{\text{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 $\mathbf{p}$:

$\displaystyle\ket{p\,q\,\alpha}\equiv\ket{p\,q\,[(LS)J\,(\ell
s)j]\mathcal{J}\mathcal{M}_{\mathcal{J}}\,(Tt)\mathcal{T}\mathcal{M}_{\mathcal{T}}}.$
(25)

Here
$\mathbf{q}=\frac{2}{3}\left(\mathbf{k}_{3}-\frac{1}{2}(\mathbf{k}_{1}+\mathbf{k}_{2})\right)$
is the second Jacobi momentum, whereas $\ell$, $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 $\mathcal{J}$ and
$\mathcal{T}$, respectively. Here again, the collective index
$\alpha=\\{L,S,J,T,l,s,j,\mathcal{J},\mathcal{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

$\displaystyle\ket{n_{i}(\ell_{i}s_{i})j_{i}m_{j_{i}}\,t_{i}m_{t_{i}}},$ (26)

while the two-body basis states become

$\displaystyle\ket{N\,\alpha_{\text{2b}}}\equiv\ket{N(LS)JM_{J}\,TM_{T}},$
(27)

and the three-body HO basis states are specified by

$\displaystyle\ket{N\,n\,\alpha}\equiv\ket{N\,n\,[(LS)J\,(\ell
s)j]\mathcal{J}\mathcal{M}_{\mathcal{J}}\,(Tt)\mathcal{T}\mathcal{M}_{\mathcal{T}}}.$
(28)

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.

###  III.1 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 partial-wave basis defined in Eq. (25) [40, 39]:

$\displaystyle\left\langle
p^{\prime}q^{\prime}\alpha^{\prime}\left|\mathbf{j}(\mathbf{Q})\right|pq\alpha\right\rangle=$
$\displaystyle\sum_{M_{J}M_{J}^{\prime}m_{j}m_{j}^{\prime}}\mathcal{C}_{J^{\prime}M^{\prime}_{J}jm^{\prime}_{j}}^{\mathcal{J}^{\prime}\mathcal{M}_{\mathcal{J}}^{\prime}}\mathcal{C}_{JM_{J}jm_{j}}^{\mathcal{J}\mathcal{M}_{\mathcal{J}}}\sum_{M_{T}M^{\prime}_{T}m_{t}m^{\prime}_{t}}\mathcal{C}_{T^{\prime}M^{\prime}_{T}t^{\prime}m^{\prime}_{t}}^{\mathcal{T}^{\prime}\mathcal{M}^{\prime}_{\mathcal{T}}}\mathcal{C}_{TM_{T}tm_{t}}^{\mathcal{T}\mathcal{M}_{\mathcal{T}}}\mathcal{P}_{\left(L^{\prime}S^{\prime}\right)J^{\prime}T^{\prime}(LS)JT}^{M_{J}^{\prime}M_{T}^{\prime}M_{J}M_{T}}\left(\mathbf{Q},p,p^{\prime}\right)\mathcal{Q}_{\left(\ell^{\prime}s^{\prime}\right)j^{\prime}t^{\prime}(\ell
s)jt}^{m_{j}^{\prime}m_{t}^{\prime}m_{j}m_{t}}\left(\mathbf{Q},q,q^{\prime}\right),$
(29)

with

$\displaystyle\mathcal{P}_{\left(L^{\prime}S^{\prime}\right)J^{\prime}T^{\prime}(LS)JT}^{M_{J}^{\prime}M_{T}^{\prime}M_{J}M_{T}}\left(\mathbf{Q},p,p^{\prime}\right)=$
$\displaystyle\frac{1}{(2\pi)^{3}}\int\mathrm{d}\mathbf{p}_{1}\mathrm{d}\mathbf{p}_{1}^{\prime}\frac{\delta\left(p_{1}^{\prime}-p^{\prime}\right)}{p_{1}^{\prime}p^{\prime}}\mathcal{Y}_{L^{\prime}S^{\prime}}^{*J^{\prime}M_{J}^{\prime}}\left(\hat{\mathbf{p}}_{1}^{\prime}\right)\left\langle\mathbf{p}_{1}^{\prime}T^{\prime}M_{T}^{\prime}\left|\mathbf{j}(\mathbf{Q})\right|\mathbf{p}_{1}TM_{T}\right\rangle\frac{\delta\left(p-p_{1}\right)}{pp_{1}}\mathcal{Y}_{LS}^{JM_{J}}\left(\hat{\mathbf{p}}_{1}\right),$
(30)
$\displaystyle\mathcal{Q}_{\left(l^{\prime}s^{\prime}\right)j^{\prime}t^{\prime}(\ell
s)jt}^{m_{j}^{\prime}m_{t}^{\prime}m_{j}m_{t}}\left(\mathbf{Q},q,q^{\prime}\right)=$
$\displaystyle\frac{1}{(2\pi)^{3}}\int\mathrm{d}\mathbf{q}_{1}\mathrm{d}\mathbf{q}_{1}^{\prime}\frac{\delta\left(q_{1}^{\prime}-q^{\prime}\right)}{q_{1}^{\prime}q^{\prime}}\mathcal{Y}_{\ell^{\prime}s^{\prime}}^{*j_{j}^{\prime}m^{\prime}_{j}}\left(\hat{\mathbf{q}}_{1}^{\prime}\right)\left\langle\mathbf{q}_{1}^{\prime}t^{\prime}m_{t}^{\prime}\left|\mathbf{j}(\mathbf{Q})\right|\mathbf{q}_{1}tm_{t}\right\rangle\frac{\delta\left(q-q_{1}\right)}{qq_{1}}\mathcal{Y}_{\ell
s}^{jm_{j}}\left(\hat{\mathbf{q}}_{1}\right),$ (31)

and the spinor spherical harmonics

$\displaystyle\mathcal{Y}_{LS}^{JM_{J}}\left(\hat{\mathbf{a}}\right)=\sum_{M_{L},M_{S}}\mathcal{C}_{LM_{L}SM_{S}}^{JM_{J}}Y_{LM_{L}}\left(\hat{\mathbf{a}}\right)\left|SM_{S}\right>.$
(32)

This factorized representation is very useful in practice as one-body operator
can be represented in a natural way in terms of the quantity
$\mathcal{Q}\left(\mathbf{Q},q,q^{\prime}\right)$, while for two-body
operators the dynamics of the interaction with the external probe can be
parametrized naturally via $\mathcal{P}\left(\mathbf{Q},p,p^{\prime}\right)$.

####  III.1.1 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 $\mathbf{q}$ resp.
$\mathbf{q}^{\prime}$ describing its motion relative to the subsystem of the
other two nucleons. Specifically, we have:

$\displaystyle\mathbf{k}^{\prime}_{1}=\mathbf{k}_{1},\quad\mathbf{k}^{\prime}_{2}=\mathbf{k}_{2},\quad\mathbf{k}^{\prime}_{3}=\mathbf{k}_{3}+\mathbf{Q},$
(33)

which implies for the Jacobi momenta:

$\mathbf{p}^{\prime}=\mathbf{p},\quad\mathbf{q}^{\prime}=\mathbf{q}+\frac{2}{3}\mathbf{Q}\,.$
(34)

Consequently, the matrix elements of $\mathcal{P}(\mathbf{Q},p,p^{\prime})$
are trivial since both particles with the relative momentum
$\mathbf{p}=\mathbf{p}^{\prime}$ are just spectator particles:

$\displaystyle\mathcal{P}_{\left(L^{\prime}S^{\prime}\right)J^{\prime}T^{\prime}(LS)JT}^{M_{J}^{\prime}M_{T}^{\prime}M_{J}M_{T}}\left(\mathbf{Q},p,p^{\prime}\right)=\delta_{M_{T}M_{T}^{\prime}}\delta_{M_{J}M^{\prime}_{J}}\delta_{JJ^{\prime}}\delta_{LL^{\prime}}\delta_{SS^{\prime}}\delta_{TT^{\prime}}\frac{\delta\left(p-p^{\prime}\right)}{p^{\prime}p}\,.$
(35)

The quantity $\mathcal{Q}(\mathbf{Q},q,q^{\prime})$ on the other hand contains
the dynamics of the probe described by the one-body current and is given by:

$\displaystyle\mathcal{Q}_{\left(l^{\prime}s^{\prime}\right)j^{\prime}t^{\prime}(\ell
s)jt}^{m_{j}^{\prime}m_{t}^{\prime}m_{j}m_{t}}\left(\mathbf{Q},q,q^{\prime}\right)=\int\mathrm{d}\hat{\mathbf{q}}\frac{\delta\left(q^{\prime}-|\mathbf{q}+\frac{2}{3}\mathbf{Q}|\right)}{q^{\prime
2}}\mathcal{Y}_{\ell^{\prime}s^{\prime}}^{*j_{j}^{\prime}m^{\prime}_{j}}\left(\widehat{\mathbf{q}+\frac{2}{3}\mathbf{Q}}\right)\left\langle(\mathbf{q}+2/3\mathbf{Q})t^{\prime}m_{t}^{\prime}\left|\mathbf{j}(\mathbf{Q})\right|\mathbf{q}tm_{t}\right\rangle\mathcal{Y}_{\ell
s}^{jm_{j}}\left(\hat{\mathbf{q}}\right)\,.$ (36)

The evaluation is straightforward and can be performed numerically or
partially analytically, depending on the specific form of the current
operator.

####  III.1.2 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 $\mathbf{p}$ resp. $\mathbf{p}^{\prime}$
to interact with the external probe, i.e.,
$\mathbf{k}^{\prime}_{3}=\mathbf{k}_{3}$. We first express the single-particle
momenta in terms of the Jacobi and center-of-mass momenta [39]:

$\displaystyle\mathbf{k}_{1}$
$\displaystyle=\mathbf{p}-\frac{\mathbf{q}}{2}+\frac{1}{3}\mathbf{P}_{\mathrm{3N}}$
(37) $\displaystyle\mathbf{k}_{2}$
$\displaystyle=-\mathbf{p}-\frac{\mathbf{q}}{2}+\frac{1}{3}\mathbf{P}_{\mathrm{3N}}$
(38) $\displaystyle\mathbf{k}_{3}$
$\displaystyle=\mathbf{q}+\frac{1}{3}\mathbf{P}_{\mathrm{3N}}.$ (39)

As a result of the interaction process, the center-of-mass momentum changes,
$\mathbf{P}^{\prime}_{\mathrm{3N}}=\mathbf{P}_{\mathrm{3N}}+\mathbf{Q}$, and
hence:

$\displaystyle\mathbf{q}^{\prime}=\mathbf{q}+\frac{1}{3}\mathbf{Q},$ (40)

while the Jacobi momentum $\mathbf{p}=(\mathbf{k}_{2}-\mathbf{k}_{1})/2$ and
$\mathbf{p}^{\prime}=(\mathbf{k}^{\prime}_{2}-\mathbf{k}^{\prime}_{1})/2$ are
again unaffected by the interaction with the external probe. Overall, we
obtain:

$\displaystyle\mathcal{P}_{\left(L^{\prime}S^{\prime}\right)J^{\prime}T^{\prime}(LS)JT}^{M_{J}^{\prime}M_{T}^{\prime}M_{J}M_{T}}\left(\mathbf{Q},p,p^{\prime}\right)$
$\displaystyle=\frac{1}{(2\pi)^{3}}\int\mathrm{d}\hat{\mathbf{p}}^{\prime}d\hat{\mathbf{p}}\,\mathcal{Y}_{L^{\prime}S^{\prime}}^{*J^{\prime}M_{J}^{\prime}}\left(\hat{\mathbf{p}}^{\prime}\right)\left\langle\mathbf{p}^{\prime}T^{\prime}M_{T}^{\prime}\left|\mathbf{j}(\mathbf{Q})\right|\mathbf{p}TM_{T}\right\rangle\mathcal{Y}_{LS}^{JM_{J}}\left(\hat{\mathbf{p}}\right)$
(41)
$\displaystyle\mathcal{Q}_{\left(l^{\prime}s^{\prime}\right)j^{\prime}t^{\prime}(\ell
s)jt}^{m_{j}^{\prime}m_{t}^{\prime}m_{t}m_{t}}\left(\mathbf{Q},q,q^{\prime}\right)$
$\displaystyle=\delta_{m_{t}m^{\prime}_{t}}\delta_{ss^{\prime}}\delta_{m_{s}m^{\prime}_{s}}\int\mathrm{d}\hat{\mathbf{q}}\frac{\delta\left(q^{\prime}-|\mathbf{q}+\frac{1}{3}\mathbf{Q}|\right)}{q^{\prime
2}}\mathcal{Y}_{\ell^{\prime}s^{\prime}}^{*j_{j}^{\prime}m_{j}}\left(\widehat{\mathbf{q}+\frac{1}{3}\mathbf{Q}}\right)\mathcal{Y}_{\ell
s}^{jm_{j}}\left(\hat{\mathbf{q}}\right)\,.$ (42)

The quantity $\mathcal{Q}$ is a current-independent function that in this case
depends only on the kinematics specified by the momentum $\mathbf{Q}$, while
the two-body quantity $\mathcal{P}$ contains all the information about the
current operator.

###  III.2 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 momentum $\mathbf{Q}$ to be along the $z$ direction. In
the Jacobi NCSM, a wave function $|\Psi_{\rm 3b}\rangle$ is obtained through
diagonalization of the Hamiltonian and expressed as a superposition of
antisymmetrized HO basis states:

$|\Psi_{\rm 3b}\rangle=\sum_{i}c_{i}|i\rangle.$ (43)

Note that the antisymmetrized HO basis $|i\rangle$ is not the same as the
state defined in Eq. (28) and computed as

$|i\rangle=\sum_{Nn\alpha}|Nn\alpha\rangle\langle Nn\alpha|i\rangle,$ (44)

with the coefficient of fractional parentage $\langle Nn\alpha|i\rangle$ [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

$\langle\Psi_{\rm 3b}^{\prime}|\sum_{i}\mu_{{\rm 1b},i}|\Psi_{\rm
3b}\rangle=3\sum_{kl}\sum_{N^{\prime}n^{\prime}\alpha^{\prime}}\sum_{Nn\alpha}c^{\prime*}_{k}c_{l}\\\
\times\langle k|N^{\prime}n^{\prime}\alpha^{\prime}\rangle\langle
Nn\alpha|l\rangle\langle N^{\prime}n^{\prime}\alpha^{\prime}|\mu_{{\rm
1b},3}|Nn\alpha\rangle,$ (45)

for one-body operators. A similar expression can be found for two-body
operators:

$\langle\Psi_{\rm 3b}^{\prime}|\sum_{i<j}\mu_{{\rm 2b},ij}|\Psi_{\rm
3b}\rangle=3\sum_{kl}\sum_{N^{\prime}n^{\prime}\alpha^{\prime}}\sum_{Nn\alpha}c^{\prime*}_{k}c_{l}\\\
\times\langle k|N^{\prime}n^{\prime}\alpha^{\prime}\rangle\langle
Nn\alpha|l\rangle\langle N^{\prime}n^{\prime}\alpha^{\prime}|\mu_{{\rm
2b},12}|Nn\alpha\rangle.$ (46)

The main tasks are to find expressions for

$\langle N^{\prime}n^{\prime}\alpha^{\prime}|||\mu^{\kappa}_{{\rm
1b},3}|||Nn\alpha\rangle,$ (47)

and

$\langle N^{\prime}n^{\prime}\alpha^{\prime}|||\mu^{\kappa}_{{\rm
2b},12}|||Nn\alpha\rangle.$ (48)

Here, we introduced doubly reduced matrix element with respect to spin and
isospin, where $\kappa$ 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.

####  III.2.1 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
$\tilde{R}_{n\ell}(k)=\innerproduct{k\,\ell}{n\,\ell}$, defined by the overlap
between momentum and HO eigenstates, given by the slightly modified definition
from Ref. [43]

$\tilde{R}_{n\ell}(k)=\sqrt{\frac{2n!b^{3}}{\Gamma(n+\ell+3/2)}}\,\big{(}kb\big{)}^{\ell}\,e^{-\frac{1}{2}k^{2}b^{2}}\,L_{n}^{\ell+1/2}(k^{2}b^{2}),$
(49)

with $b\equiv 1/\sqrt{m_{N}\omega}$ the oscillator length in terms of the
oscillator frequency $\omega$ and the nucleon mass $m_{N}$, and
$L_{n}^{\ell}(x)$ are generalized Laguerre polynomials. Similarly, coordinate-
space radial wave functions $R_{n\ell}(r)=\innerproduct{r\,\ell}{n\,\ell}$ are
given by

$\displaystyle R_{n\ell}(r)=$
$\displaystyle(-1)^{n}\sqrt{\frac{2n!}{\Gamma(n+\ell+3/2)b^{3}}}\,\Bigg{(}\frac{r}{b}\Bigg{)}^{\ell}\,e^{-\frac{1}{2}(\frac{r}{b})^{2}}L_{n}^{\ell+1/2}\Bigg{(}\frac{r^{2}}{b^{2}}\Bigg{)},$
(50)

which are connected to the momentum-space functions through a Fourier-Bessel
transform

$\displaystyle
R_{n\ell}(r)=\int\\!\mathrm{d}k\,k^{2}\innerproduct{r\,\ell}{k\,\ell}\tilde{R}_{n\ell}(k),$
(51)

where the overlap is described by spherical Bessel functions
$\innerproduct{r\,\ell}{k\,\ell}=\sqrt{2/\pi}j_{\ell}(kr)$.

####  III.2.2 One-body operator

Matrix elements of the one-body magnetic dipole operator defined in Eqs.
(13)-(16) are given by

$\displaystyle\langle n^{\prime}\ell^{\prime}j^{\prime}|||\mu^{0}_{\rm
spin}|||n\ell j\rangle$
$\displaystyle=2\sqrt{3}G^{S}_{M}(0)\hat{j}^{\prime}\hat{j}(-1)^{\ell^{\prime}+j^{\prime}+3/2}$
$\displaystyle\quad\times\left\\{\begin{array}[]{ccc}1/2&1/2&1\\\
j^{\prime}&j&\ell\end{array}\right\\}\delta_{nn^{\prime}}\delta_{\ell\ell^{\prime}},$
(54) $\displaystyle\langle n^{\prime}\ell^{\prime}j^{\prime}|||\mu^{1}_{\rm
spin}|||n\ell j\rangle$
$\displaystyle=6G^{V}_{M}(0)\hat{j}^{\prime}\hat{j}(-1)^{\ell^{\prime}+j^{\prime}+3/2}$
$\displaystyle\quad\times\left\\{\begin{array}[]{ccc}1/2&1/2&1\\\
j^{\prime}&j&\ell\end{array}\right\\}\delta_{nn^{\prime}}\delta_{\ell\ell^{\prime}},$
(57) $\displaystyle\langle n^{\prime}\ell^{\prime}j^{\prime}|||\mu^{0}_{\rm
orb}|||n\ell j\rangle$
$\displaystyle=\frac{\mu_{N}}{\sqrt{2}}\hat{j}^{\prime}\hat{j}\hat{\ell}\sqrt{\ell(\ell+1)}(-1)^{\ell+j+3/2}$
$\displaystyle\quad\times\left\\{\begin{array}[]{ccc}\ell&\ell&1\\\
j^{\prime}&j&1/2\end{array}\right\\}\delta_{nn^{\prime}}\delta_{\ell\ell^{\prime}},$
(60) $\displaystyle\langle n^{\prime}\ell^{\prime}j^{\prime}|||\mu^{1}_{\rm
orb}|||n\ell j\rangle$
$\displaystyle=\sqrt{\frac{3}{2}}\mu_{N}\hat{j}^{\prime}\hat{j}\hat{\ell}\sqrt{\ell(\ell+1)}(-1)^{\ell+j+3/2}$
$\displaystyle\quad\times\left\\{\begin{array}[]{ccc}\ell&\ell&1\\\
j^{\prime}&j&1/2\end{array}\right\\}\delta_{nn^{\prime}}\delta_{\ell\ell^{\prime}}.$
(63)

For the spin term $\mu^{\kappa}_{\rm spin}$, the required reduced matrix
element for Eq. (47) is

$\displaystyle\langle N^{\prime}n^{\prime}\alpha|||$
$\displaystyle\mu^{\kappa}_{{\rm
spin},3}|||Nn\alpha\rangle=(-1)^{\mathcal{J}^{\prime}+J+j+1}\hat{\mathcal{J}}^{\prime}\hat{\mathcal{J}}\left\\{\begin{array}[]{ccc}j^{\prime}&\mathcal{J}^{\prime}&J\\\
\mathcal{J}&j&1\end{array}\right\\}$ (64)
$\displaystyle\times(-1)^{\mathcal{T}^{\prime}+T+3/2}\hat{\mathcal{T}}^{\prime}\hat{\mathcal{T}}\left\\{\begin{array}[]{ccc}1/2&\mathcal{T}^{\prime}&T\\\
\mathcal{T}&1/2&\kappa\end{array}\right\\}$ $\displaystyle\times\langle
n^{\prime}\ell^{\prime}j^{\prime}|||\mu^{\kappa}_{\rm spin}|||n\ell
j\rangle\delta_{N^{\prime}N}\delta_{L^{\prime}L}\delta_{S^{\prime}S}\delta_{J^{\prime}J}\delta_{T^{\prime}T}.$

A similar expression can be found for the orbital contribution
$\mu^{\kappa}_{\rm orb}$:

$\displaystyle\langle N^{\prime}n^{\prime}\alpha|||$
$\displaystyle\mu^{\kappa}_{{\rm
orb},3}|||Nn\alpha\rangle=\frac{2}{3}(-1)^{\mathcal{J}^{\prime}+J+j+1}\hat{\mathcal{J}}^{\prime}\hat{\mathcal{J}}\left\\{\begin{array}[]{ccc}j^{\prime}&\mathcal{J}^{\prime}&J\\\
\mathcal{J}&j&1\end{array}\right\\}$ (65)
$\displaystyle\times(-1)^{\mathcal{T}^{\prime}+T+3/2}\hat{\mathcal{T}}^{\prime}\hat{\mathcal{T}}\left\\{\begin{array}[]{ccc}1/2&\mathcal{T}^{\prime}&T\\\
\mathcal{T}&1/2&\kappa\end{array}\right\\}$ $\displaystyle\times\langle
n^{\prime}\ell^{\prime}j^{\prime}|||\mu^{\kappa}_{\rm orb}|||n\ell
j\rangle\delta_{N^{\prime}N}\delta_{L^{\prime}L}\delta_{S^{\prime}S}\delta_{J^{\prime}J}\delta_{T^{\prime}T}.$

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.

####  III.2.3 Intrinsic magnetic dipole operator

After angular momentum recoupling, the reduced matrix element of the intrinsic
magnetic dipole operator (II.2) with respect to a relative two-body HO basis
(27) can be computed through

$\displaystyle\langle$ $\displaystyle N^{\prime}\alpha^{\prime}_{\rm
2b}|||\mu^{\rm NLO,int}_{\rm 2b}|||N\alpha_{\rm
2b}\rangle=-\frac{g^{2}_{A}m_{\pi}}{32\pi F^{2}_{\pi}}\sqrt{\frac{8\pi}{3}}$
(66) $\displaystyle\times(i\langle
T^{\prime}||(\bm{\tau}_{1}\times\bm{\tau}_{2})||T\rangle)$
$\displaystyle\times\int
drr^{2}R_{N^{\prime}L^{\prime}}(r)R_{NL}(r)e^{-\sqrt{2}m_{\pi}r}$
$\displaystyle\times\sum_{w=0,2}A^{L^{\prime}S^{\prime}J^{\prime}}_{LSJ}(w,1)\left[f(\sqrt{2}r)\mathcal{C}_{1010}^{w0}\,+\sqrt{3}\delta_{w0}\right],$

with

$\displaystyle A^{L^{\prime}S^{\prime}J^{\prime}}_{LSJ}(w,x)$
$\displaystyle\equiv\langle(L^{\prime}S^{\prime})J^{\prime}||[Y_{w}(\hat{\mathbf{r}})[\sigma_{1}\sigma_{2}]_{x}]_{1}||(LS)J\rangle$
(67)
$\displaystyle=6\sqrt{\frac{3}{4\pi}}\hat{J^{\prime}}\hat{J}\hat{L}\hat{S^{\prime}}\hat{S}\hat{w}\hat{x}$
$\displaystyle\quad\times\left\\{\begin{array}[]{ccc}L^{\prime}&L&w\\\
S^{\prime}&S&x\\\
J^{\prime}&J&1\end{array}\right\\}\left\\{\begin{array}[]{ccc}1/2&1/2&1\\\
1/2&1/2&1\\\
S^{\prime}&S&x\end{array}\right\\}\mathcal{C}_{L0w0}^{L^{\prime}0}\,.$

Also, the matrix element of the intrinsic two-body magnetic moment operator
for Eq. (48) can be written as

$\displaystyle\langle N^{\prime}n^{\prime}\alpha^{\prime}|||$
$\displaystyle\mu^{\rm NLO,int}_{{\rm
2b},12}|||Nn\alpha\rangle=(-1)^{J^{\prime}+j^{\prime}+\mathcal{J}+1}(-1)^{T^{\prime}+\mathcal{T}+3/2}$
(68)
$\displaystyle\times\hat{\mathcal{J}}^{\prime}\hat{\mathcal{J}}\left\\{\begin{array}[]{ccc}\mathcal{J}^{\prime}&J^{\prime}&j\\\
J&\mathcal{J}&1\end{array}\right\\}$
$\displaystyle\times\hat{\mathcal{T}}^{\prime}\hat{\mathcal{T}}\left\\{\begin{array}[]{ccc}\mathcal{T}^{\prime}&T^{\prime}&1/2\\\
T&\mathcal{T}&1\end{array}\right\\}$ $\displaystyle\times\langle
N^{\prime}\alpha^{\prime}_{2b}|||\mu^{\rm NLO,int}_{\rm
2b}|||N\alpha_{2b}\rangle\delta_{n^{\prime}n}\delta_{\ell^{\prime}\ell}\delta_{j^{\prime}j}.$

####  III.2.4 Sachs operator

In the previous section, we described how to embed a two-body 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

$\ket{N_{\text{NN}}N_{12}\,\alpha_{12}}\\\
\equiv\ket{N_{\text{NN}}N_{12}[L_{\text{NN}}(L_{12}S_{12})J_{12}]J_{\text{rc}}M_{J_{\text{rc}}}T_{12}M_{T_{12}}},$
(69)

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-of-mass angular momenta generates the total
angular momentum of the two-body system $J_{\text{rc}}$ with projection
$M_{J_{\text{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$.

Figure 3: Schematic of two different coordinate systems representing a three-
nucleon system. The left part represents the Jacobi coordinate system
characterized by ($\mathbf{P}_{\text{3N}}$, $\mathbf{p}$, $\mathbf{q}$), while
the right part shows the coordinate system described by
($\mathbf{P}_{\text{NN}}$, $\mathbf{p}$, $\mathbf{k}_{3}$).

The matrix element of the Sachs operator with respect to the states defined in
Eq. 69 is given by

$\displaystyle\langle$ $\displaystyle N^{\prime}_{\rm
NN}N^{\prime}_{12}\alpha^{\prime}_{12}|||\mu^{\rm NLO,Sachs}_{\rm 2b}|||N_{\rm
NN}N_{12}\alpha_{12}\rangle=\frac{g^{2}_{A}m^{2}_{\pi}}{24F^{2}_{\pi}}$ (70)
$\displaystyle\times i\langle
T^{\prime}||(\bm{\tau}_{1}\times\bm{\tau}_{2})||T\rangle$
$\displaystyle\times\int
drr^{2}R_{N^{\prime}_{12}L^{\prime}_{12}}(r)R_{N_{12}L_{12}}(r)e^{-\sqrt{2}m_{\pi}r}$
$\displaystyle\times\sum_{x=0,2}\sum_{w}\left[\left(\mathcal{C}_{1010}^{x0}\,\mathcal{C}_{10x0}^{w0}\,+\sqrt{\frac{1}{3}}\delta_{w1}\delta_{x0}\right)h(\sqrt{2}r)-\sqrt{\frac{1}{3}}\delta_{w1}\delta_{x0}\right]$
$\displaystyle\times\int dRR^{3}R_{N_{\rm NN}^{\prime}L_{\rm
NN}^{\prime}}(R)R_{N_{\rm NN}L_{\rm NN}}(R)B^{L_{\rm
NN}^{\prime}L_{12}^{\prime}S_{12}^{\prime}J_{12}^{\prime}J_{\rm
rc}^{\prime}}_{L_{\rm NN}L_{12}S_{12}J_{12}J_{\rm rc}}(w,x),$

with

$\displaystyle B^{L_{\rm
NN}^{\prime}L_{12}^{\prime}S_{12}^{\prime}J_{12}^{\prime}J_{\rm
rc}^{\prime}}_{L_{\rm NN}L_{12}S_{12}J_{12}J_{\rm rc}}(w,x)$
$\displaystyle\equiv\langle[L_{\rm
NN}^{\prime}(L_{12}^{\prime}S_{12}^{\prime})J_{12}^{\prime}]J_{\rm
rc}^{\prime}||[Y_{1}(\hat{\mathbf{R}})[Y_{w}(\hat{\mathbf{r}})[\sigma_{1}\sigma_{2}]_{x}]_{1}]_{1}||[L_{\rm
NN}(L_{12}S_{12})J_{12}]J_{\rm rc}^{\prime}\rangle$ (71)
$\displaystyle=3\sqrt{\frac{1}{4\pi}}\hat{J}_{\rm rc}^{\prime}\hat{J}_{\rm
rc}\hat{L}_{\rm NN}\left\\{\begin{array}[]{ccc}L_{\rm NN}^{\prime}&L_{\rm
NN}&1\\\ J_{12}^{\prime}&J_{12}&1\\\ J_{\rm rc}^{\prime}&J_{\rm
rc}&1\end{array}\right\\}\mathcal{C}_{L_{\rm NN}010}^{L_{\rm
NN}^{\prime}0}\,A^{L_{12}^{\prime}S_{12}^{\prime}J_{12}^{\prime}}_{L_{12}S_{12}J_{12}}(w,x).$

In addition to the integration over $r$, the matrix elements of the Sachs
operator are integrated over the center-or-mass coordinate $R$. Same as the
intrinsic contribution, the isospin rank of the operator $\kappa$ 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. 69 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. 70 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
$\ket{n_{3}(\ell_{3}\tfrac{1}{2})j_{3}}$, which is defined with respect to the
origin, and couple its angular momentum $j_{3}$ with the total two-body
angular momentum $J_{\text{rc}}$ to the total three-body angular momentum
$\mathcal{J}_{\text{tot}}$:

$\ket{N_{\text{NN}}N_{12}n_{3}\,\alpha_{12}\,\alpha_{3}\,\mathcal{J}_{\text{tot}}\mathcal{M}_{\mathcal{J}_{\text{tot}}}\mathcal{T}\mathcal{M}_{\mathcal{T}}}\\\
\begin{split}\equiv&\ket{N_{\text{NN}}N_{12}n_{3}\\{[L_{\text{NN}}(L_{12}S_{12})J_{12}]J_{\text{rc}}(\ell_{3}\tfrac{1}{2})j_{3}\\}\mathcal{J}_{\text{tot}}\mathcal{M}_{\mathcal{J}_{\text{tot}}}}\\\
&\times\ket{(T_{12}\tfrac{1}{2})\mathcal{T}\mathcal{M}_{\mathcal{T}}}.\end{split}$
(72)

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. (72)
can be represented in terms of reduced two-body matrix elements from Eq. 70,
angular momentum and isospin coupling factors, and a third particle which is
diagonal in all its quantum numbers:

$\displaystyle\langle$ $\displaystyle N^{\prime}_{\rm
NN}N^{\prime}_{12}n^{\prime}_{3}\alpha^{\prime}_{12}\alpha^{\prime}_{3}\mathcal{J}^{\prime}_{\rm
tot}\mathcal{T}^{\prime}|||\mu^{\rm NLO,Sachs}_{{\rm 2b},12}|||N_{\rm
NN}N_{12}n_{3}\alpha_{12}\alpha_{3}\mathcal{J}_{\rm tot}\mathcal{T}\rangle$
(73) $\displaystyle=(-1)^{J^{\prime}_{\rm rc}+j^{\prime}_{3}+\mathcal{J}_{\rm
tot}+1}\hat{\mathcal{J}}^{\prime}_{\rm tot}\hat{\mathcal{J}}_{\rm
tot}\left\\{\begin{array}[]{ccc}\mathcal{J}^{\prime}_{\rm tot}&J^{\prime}_{\rm
rc}&j_{3}\\\ J_{\rm rc}&\mathcal{J}_{\rm tot}&1\end{array}\right\\}$
$\displaystyle\times(-1)^{T^{\prime}+\mathcal{T}+3/2}\hat{\mathcal{T}}^{\prime}\hat{\mathcal{T}}\left\\{\begin{array}[]{ccc}\mathcal{T}^{\prime}&T^{\prime}&1/2\\\
T&\mathcal{T}&1\end{array}\right\\}$ $\displaystyle\times\langle
N^{\prime}_{\rm NN}N^{\prime}_{12}\alpha^{\prime}_{12}|||\mu^{\rm
NLO,Sachs}_{\rm 2b}|||N_{\rm
NN}N_{12}\alpha_{12}\rangle\delta_{n^{\prime}_{3}n_{3}}\delta_{\ell^{\prime}_{3}\ell_{3}}\delta_{j^{\prime}_{3}j_{3}}$

This result shows that the majority of the work to calculate three-body matrix
elements consists in determining the two-body matrix elements.

To calculate the overlap between the three-body basis states defined in Eq. 72
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 $\mathcal{J}$ with the center-of-mass orbital angular momentum
$\mathcal{L}_{\text{3N}}$ to the total angular momentum
$\mathcal{J}_{\text{tot}}$, so that the basis from Eq. 28 is extended to

$\ket{\mathcal{N}_{\text{3N}}Nn\,\alpha_{\text{3N}}\alpha\,\mathcal{J}_{\text{tot}}\mathcal{M}_{\mathcal{J}_{\text{tot}}}}\equiv\\\
\ket{\mathcal{N}_{\text{3N}}Nn\\{\mathcal{L}_{\text{3N}}[(LS)J(\ell\tfrac{1}{2})j]\mathcal{J}\\}\mathcal{J}_{\text{tot}}\mathcal{M}_{\mathcal{J}_{\text{tot}}}(T\tfrac{1}{2})\mathcal{T}\mathcal{M}_{\mathcal{T}}},$
(74)

where the subscript 3N denotes quantities related to the three-body 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:

$\langle\mathcal{N}_{\rm 3N}^{\prime}N^{\prime}n^{\prime}\alpha_{\rm
3N}^{\prime}\alpha^{\prime}\mathcal{J}_{\rm
tot}^{\prime}\mathcal{T}^{\prime}|||\mu^{\rm NLO,Sachs}_{{\rm
2b},12}|||\mathcal{N}_{\rm 3N}Nn\alpha_{\rm 3N}\alpha\mathcal{J}_{\rm
tot}\mathcal{T}\rangle\\\
\begin{split}=&\sum_{N_{12}^{\prime}\,N_{\text{NN}}^{\prime}\,n_{3}^{\prime}}\sum_{\alpha^{\prime}_{12}\,\alpha^{\prime}_{3}}\sum_{N_{12}\,N_{\text{NN}}\,n_{3}}\sum_{\alpha_{12}\,\alpha_{3}}\innerproduct{\mathcal{N}_{\text{3N}}^{\prime}N^{\prime}n^{\prime}\,\alpha^{\prime}_{\text{3N}}\,\alpha^{\prime}\,\mathcal{J}_{\text{tot}}}{N_{\text{NN}}^{\prime}N_{12}^{\prime}n_{3}^{\prime}\,\alpha^{\prime}_{12}\,\alpha^{\prime}_{3}\,\mathcal{J}_{\text{tot}}}\\\
&\times\langle N^{\prime}_{\rm
NN}N^{\prime}_{12}n^{\prime}_{3}\alpha^{\prime}_{12}\alpha^{\prime}_{3}\mathcal{J}^{\prime}_{\rm
tot}\mathcal{T}^{\prime}|||\mu^{\rm NLO,Sachs}_{{\rm 2b},12}|||N_{\rm
NN}N_{12}n_{3}\alpha_{12}\alpha_{3}\mathcal{J}_{\rm tot}\mathcal{T}\rangle\\\
&\times\innerproduct{N_{\text{NN}}N_{12}n_{3}\,\alpha_{12}\,\alpha_{3}\,\mathcal{J}_{\text{tot}}}{N_{\text{3N}}Nn\,\alpha_{\text{3N}}\,\alpha\,\mathcal{J}_{\text{tot}}}.\end{split}$
(75)

The remaining task is to determine the basis transformation brakets
$\innerproduct{N_{\text{NN}}N_{12}n_{3}\,\alpha_{12}\,\alpha_{3}\,\mathcal{J}_{\text{tot}}}{\mathcal{N}_{\text{3N}}Nn\,\alpha_{\text{3N}}\,\alpha\,\mathcal{J}_{\text{tot}}}$,
as the matrix element with respect to the three-body basis is already given in
Eq. 73. Again, since the intrinsic and center-of-mass motions are exactly
decoupled, we are free to choose any $\mathcal{N}_{\rm 3N}$ and
$\mathcal{L}_{\rm 3N}$. The most convenient choice is $\mathcal{N}_{\rm
3N}=\mathcal{L}_{\rm 3N}=0$. Then, the overlap is found to be

$\displaystyle\langle N_{\rm NN}N_{12}$ $\displaystyle
n_{3}\alpha_{12}\alpha_{3}\mathcal{J}_{\rm tot}|0Nn\alpha\mathcal{J}_{\rm
tot}\rangle=(-1)^{L_{\rm NN}-J_{\rm rc}+\ell_{3}+1/2+j_{3}+j+\mathcal{J}}$
(76) $\displaystyle\quad\times\hat{J}_{\rm
rc}\hat{j}\hat{j}_{3}\hat{l}\left\\{\begin{array}[]{ccc}J&L_{\rm NN}&J_{\rm
rc}\\\
j_{3}&\mathcal{J}&j\end{array}\right\\}\left\\{\begin{array}[]{ccc}L_{\rm
NN}&\ell_{3}&\ell\\\ 1/2&j&j_{3}\end{array}\right\\}$
$\displaystyle\quad\times\langle 00n\ell|N_{\rm NN}L_{\rm
NN}n_{3}\ell_{3}\rangle_{d=2}$
$\displaystyle\quad\times\delta_{NN_{12}}\delta_{LL_{12}}\delta_{SS_{12}}\delta_{JJ_{12}}\delta_{TT_{12}}\delta_{\mathcal{J}\mathcal{J}_{\rm
tot}}.$

Here, the Talmi-Moshinsky brackets are necessary to transform between the
three- and two-body center-of-mass systems. 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. (75) takes the required form
in Eq. (48) as we set $\mathcal{N}_{\rm 3N}=\mathcal{L}_{\rm 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.1
with corresponding partial-wave expanded wave functions, which are obtained by
solving the three-body Faddeev equations. Magnetic moments are calculated as
the zero-momentum-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.2, 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 N3LO with cutoffs $\Lambda=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 $\Lambda_{\text{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 N3LO, with a cutoff
$\Lambda=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
$\Lambda_{\text{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

$\displaystyle
F_{\text{M}}(Q)=-\frac{2m_{N}}{Q}\matrixelement{\Psi^{\mathcal{J},\mathcal{T}}_{\mathcal{M}_{\mathcal{J}}+1,\mathcal{M}_{\mathcal{T}}}({\rm
F})}{\,j_{+}(Q)\,}{\Psi^{\mathcal{J},\mathcal{T}}_{\mathcal{M}_{\mathcal{J}},\mathcal{M}_{\mathcal{T}}}({\rm
F})},$ (77)

and

$\displaystyle\mu_{\text{t/h}}=\matrixelement{\Psi^{\mathcal{J},\mathcal{T}}_{\mathcal{M}_{\mathcal{J}},\mathcal{M}_{\mathcal{T}}}({\rm
NCSM})}{\,\mu_{z}\,}{\Psi^{\mathcal{J},\mathcal{T}}_{\mathcal{M}_{\mathcal{J}},\mathcal{M}_{\mathcal{T}}}({\rm
NCSM})},$ (78)

respectively. Here,
$|\Psi^{\mathcal{J},\mathcal{T}}_{\mathcal{M}_{\mathcal{J}}},({\rm F})\rangle$
and $|\Psi^{\mathcal{J},\mathcal{T}}_{\mathcal{M}_{\mathcal{J}}},({\rm
NCSM})\rangle$ are the Faddeev and NCSM wave functions, respectively, with
total three-body spin and isospin specified by $\mathcal{J}=\mathcal{T}=1/2$
and maximally projected total angular momentum states, i.e.,
$\mathcal{M}_{\mathcal{J}}=\pm\mathcal{J}$. The isospin projection
$\mathcal{M}_{\mathcal{T}}={-}1/2$ or $\mathcal{M}_{\mathcal{T}}=1/2$
determines whether the wave function represents a triton or a helion,
respectively.

Figure 4: Triton (left column) and helion (right column) magnetic form
factors, in units of $\mu_{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 N3LO
with cutoff $\Lambda=500$, $450$, and $420$ MeV, respectively. Light and dark
shaded bands represent the $95\%$ and $68\%$ DoBs.

In order to perform Faddeev calculations, we truncate the basis by choosing a
maximal value for the relative total two-body angular momentum $J$. Our
calculations include partial waves up to $J\leq 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.

###  IV.1 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
$\mu_{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 $\Lambda=420$, $450$, and $500$ MeV at N3LO 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\sim 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.

Figure 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 $\Lambda=500$ MeV to $\Lambda=600$ MeV.

The magnetic form-factor results for 3He 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 3H, 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\sim 3$ fm-1; nevertheless they
remain within the $68\%$ DoB intervals. At $Q=0$, the results for
$\Lambda=500$ MeV, $F_{\text{M, t}}(0)=2.583$ $\mu_{N}$ and $F_{\text{M,
h}}(0)=-1.767$ $\mu_{N}$, deviate from the experimental values for the triton
($2.9789624659(59)$ $\mu_{N}$) and helion ($-2.127625307(25)$ $\mu_{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
sizeable at $Q=0$ [53, 21], 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 consequence 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 one-body 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.

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 one-body 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 $\Lambda=500$ MeV to $\Lambda=600$
MeV of the employed chiral interaction and therefore have a different
interpretation compared to our bands.

Our one-body current results for 3He with cutoffs $\Lambda=450$ and $500$ MeV
fall within the yellow band over the entire momentum range, while our results
for 3H 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 N3LO, 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_{\textnormal{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 circles 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 N2LO.

Figure 6: Triton and helion magnetic moments in units of $\mu_{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, N2LO + NLO 2BC, and N3LO + 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 N2LO current operator
corrections, respectively, from Ref. [21]. The black diamonds represent the
experimental values for both nuclei [52].

The magnetic moments can be understood as one-nucleon hole with respect to
4He, and the single-particle limit for triton (helion) of 2.793 (-1.913)
$\mu_{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 ${\sim}10\%$ and improves agreement with experiment. Our
final result, labeled “N3LO + 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 N3LO 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 $\mu_{S}$ and the isovector combination $\mu_{V}$ of
the trinucleon magnetic moments to constrain the new LECs [7, 21].111The
isoscalar ($\mu_{S}$) and isovector ($\mu_{V}$) combinations of the trinucleon
magnetic moments are defined by
$\mu_{S}\equiv\mu_{\textnormal{t}}+\mu_{\textnormal{h}}$ and
$\mu_{V}\equiv\mu_{\textnormal{t}}-\mu_{\textnormal{h}}$. As a result, the
experimental 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$.

###  IV.2 Test of magnetic dipole operator

Figure 7: Convergence of the triton (left column) and helion (right column)
magnetic moments as a function of $N_{\text{max}}$ for four different HO
frequencies $\hbar\omega=10,20,30,$ and $40$ MeV obtained with the EM NN
interaction with cutoff $\Lambda=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 $\mu_{N}$.

In this section we present the trinucleon magnetic moments obtained from the
magnetic dipole operator with the Jacobi NCSM, as discussed in Section III.2.
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 N3LO
interaction with cutoff $\Lambda=500$ MeV, without 3N interactions.

Figure 7 shows the convergence of the triton (left) and helion (right)
magnetic moments as a function of $N_{\text{max}}$ for four different HO
frequencies $\hbar\omega=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 “$\mu_{\text{HO,
1b}}+\mu_{\text{HO, 2b}}$”, with the corresponding form-factor normalization
in absolute terms. The two-body contribution consists of the intrinsic and
Sachs terms, i.e., $\mu_{\text{HO, 2b}}=\mu^{\text{NLO,
intrinsic}}_{2b}+\mu^{\text{NLO, Sachs}}_{2b}$, where the majority of the
correction is contributed by the intrinsic component.

Table 1: Triton and helion magnetic moments and their (cumulative) contributions, in $\mu_{N}$, from the form factor normalization, the magnetic dipole operator, as well as experimental values [52]. The NCSM wave function is computed at $N_{\rm max}=40$ and $\hbar\omega=20$ MeV. |  $\mu_{\textnormal{t}}$ $[\mu_{N}]$ |  |  $\mu_{\textnormal{h}}$ $[\mu_{N}]$
---|---|---|---
| $F_{\textnormal{M}}(0)$ | $\hat{\bm{\mu}}$ |  | $F_{\textnormal{M}}(0)$ | $\hat{\bm{\mu}}$
LO | $2.622$ | $2.622$ |  | $-1.783$ | $-1.783$
NLO | $2.838$ | $2.837$ |  | $-1.995$ | $-1.994$
intrinsic | - | $0.195$ |  | - | $-0.191$
Sachs | - | $0.021$ |  | - | $-0.021$
Exp. | $2.979$ |  | $-2.128$

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_{\text{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 1 displays the magnetic moments obtained from both methods. The left
column contains the results for the triton magnetic moment $\mu_{\text{t}}$
from the form-factor normalization and the magnetic dipole operator and the
right column gives the same results for the helion $\mu_{\text{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.2,
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
$\hat{\bm{\mu}}^{\textnormal{NLO, intrinsic}}_{\textnormal{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 $\ll 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 6Li 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. (66) and (70). 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 partial-wave 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\leq 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 partial-wave 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.

###### Acknowledgements.

We thank K. Wendt for useful discussions and checks. The three-body NCSM
calculations are done with NuHamil code [54]. This work was supported in part
by the BMBF Contract No. 05P21RDFNB, by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – Project-ID 279384907 – SFB 1245, the
Cluster of Excellence “Precision Physics, Fundamental Interactions, and
Structure of Matter” (PRISMA+ EXC 2118/1, Project ID 39083149), by the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant Agreement No. 101020842), by the Max
Planck Society, and by the National Science Foundation under Grant No.
PHY–2044632. This material is based upon work supported by the U.S. Department
of Energy, Office of Science, Office of Nuclear Physics, under the FRIB Theory
Alliance, award DE-SC0013617.

##  Appendix A Three-body Jacobi coordinate

For our Faddeev calculations, the Jacobi momenta are defined as

$\left(\begin{array}[]{c}\mathbf{K}_{\rm cm}\\\ \mathbf{p}\\\
\mathbf{q}\end{array}\right)=\left(\begin{array}[]{ccc}1&1&1\\\
\frac{1}{2}&-\frac{1}{2}&0\\\
-\frac{1}{3}&-\frac{1}{3}&\frac{2}{3}\end{array}\right)\left(\begin{array}[]{c}\mathbf{k}_{1}\\\
\mathbf{k}_{2}\\\ \mathbf{k}_{3}\end{array}\right),$ (79)

with the single-particle momenta $\mathbf{k}_{i}$, $i=1,2,3$. The
corresponding conjugate coordinates are given by

$\left(\begin{array}[]{c}\mathbf{R}_{\rm cm}\\\ \mathbf{r}_{p}\\\
\mathbf{r}_{q}\end{array}\right)=\left(\begin{array}[]{ccc}\frac{1}{3}&\frac{1}{3}&\frac{1}{3}\\\
1&-1&0\\\
-\frac{1}{2}&-\frac{1}{2}&1\end{array}\right)\left(\begin{array}[]{c}\mathbf{r}_{1}\\\
\mathbf{r}_{2}\\\ \mathbf{r}_{3}\end{array}\right),$ (80)

with the single-particle coordinates $\mathbf{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

$\left(\begin{array}[]{c}\bm{\pi}_{0}\\\ \bm{\pi}_{1}\\\
\bm{\pi}_{2}\end{array}\right)=\left(\begin{array}[]{ccc}\sqrt{\frac{1}{3}}&\sqrt{\frac{1}{3}}&\sqrt{\frac{1}{3}}\\\
\sqrt{\frac{1}{2}}&-\sqrt{\frac{1}{2}}&0\\\
\sqrt{\frac{1}{6}}&\sqrt{\frac{1}{6}}&-\sqrt{\frac{2}{3}}\end{array}\right)\left(\begin{array}[]{c}\mathbf{k}_{1}\\\
\mathbf{k}_{2}\\\ \mathbf{k}_{3}\end{array}\right),$ (81)

and

$\left(\begin{array}[]{c}\bm{\xi}_{0}\\\ \bm{\xi}_{1}\\\
\bm{\xi}_{2}\end{array}\right)=\left(\begin{array}[]{ccc}\sqrt{\frac{1}{3}}&\sqrt{\frac{1}{3}}&\sqrt{\frac{1}{3}}\\\
\sqrt{\frac{1}{2}}&-\sqrt{\frac{1}{2}}&0\\\
\sqrt{\frac{1}{6}}&\sqrt{\frac{1}{6}}&-\sqrt{\frac{2}{3}}\end{array}\right)\left(\begin{array}[]{c}\mathbf{r}_{1}\\\
\mathbf{r}_{2}\\\ \mathbf{r}_{3}\end{array}\right).$ (82)

##  Appendix B Three-body overlap for Sachs term evaluation

Here we give a derivation of Eq. (76). To this end, we consider the following
recoupling so that we can factor out the orbital part:

$\displaystyle|N_{\rm NN}$ $\displaystyle
N_{12}n_{3}\alpha_{12}\alpha_{3}\mathcal{J}_{\rm tot}\rangle$ (83)
$\displaystyle=(-1)^{j_{3}+\mathcal{J}_{\rm tot}-J_{\rm rc}+L_{\rm
NN}+\ell_{3}+1/2}$
$\displaystyle\quad\times\sum_{\Lambda}(-1)^{\Lambda}\hat{J_{\rm
rc}}\hat{\Lambda}\left\\{\begin{array}[]{ccc}J_{12}&L_{\rm NN}&J_{\rm rc}\\\
j_{3}&\mathcal{J}_{\rm tot}&\Lambda\end{array}\right\\}$
$\displaystyle\quad\times\sum_{\lambda}\hat{\lambda}\hat{j_{3}}\left\\{\begin{array}[]{ccc}L_{\rm
NN}&\ell_{3}&\lambda\\\ 1/2&\Lambda&j_{3}\end{array}\right\\}$
$\displaystyle\quad\times|N_{\rm NN}N_{12}n_{3}\\{J_{12}\,[(L_{\rm
NN}\ell_{3})\lambda\,\tfrac{1}{2}]\Lambda\\}\mathcal{J}_{\text{tot}}\rangle.$

To factorize the three-body center-of-mass part, one can use the coordinate
transformation

$\left(\begin{array}[]{c}\bm{\xi}_{0}\\\
\bm{\xi}_{2}\end{array}\right)=\left(\begin{array}[]{cc}\sqrt{\frac{2}{3}}&\sqrt{\frac{1}{3}}\\\
\sqrt{\frac{1}{3}}&-\sqrt{\frac{2}{3}}\end{array}\right)\left(\begin{array}[]{c}\sqrt{\frac{1}{2}}(\mathbf{r}_{1}+\mathbf{r}_{2})\\\
\mathbf{r}_{3}\end{array}\right).$ (84)

The above coordinate transformation ensures the following transformation using
the Talmi-Moshinsky bracket [55, 56], with the notation given in Ref. [57]:

$\displaystyle|N_{\rm NN}$ $\displaystyle N_{12}n_{3}\\{J_{12}[(L_{\rm
NN}\ell_{3})\lambda\frac{1}{2}]\Lambda\\}\mathcal{J}_{\rm tot}\rangle$ (85)
$\displaystyle=\sum_{\mathcal{N}_{\rm 3N}^{\prime}\mathcal{L}_{\rm
3N}^{\prime}n^{\prime}\ell^{\prime}}|\mathcal{N}_{\rm
3N}^{\prime}N_{12}n^{\prime}\\{J_{12}[(\mathcal{L}_{\rm
3N}\ell^{\prime})\lambda\frac{1}{2}]\Lambda\\}\mathcal{J}_{\rm tot}\rangle$
$\displaystyle\quad\times\langle\mathcal{N}_{\rm 3N}^{\prime}\mathcal{L}_{\rm
3N}^{\prime},n^{\prime}\ell^{\prime}:\lambda|N_{\rm NN}L_{\rm
NN},n_{3}\ell_{3}:\lambda\rangle_{d=2}.$

For the Jacobi three-body basis, the recoupling is given by

$\displaystyle|\mathcal{N}_{\rm 3N}^{\prime}$ $\displaystyle
N_{12}n^{\prime}\\{J_{12}[(\mathcal{L}_{\rm
3N}\ell^{\prime})\lambda\frac{1}{2}]\Lambda\\}\mathcal{J}_{\rm tot}\rangle$
(86)
$\displaystyle=\sum_{j^{\prime}}(-1)^{\ell^{\prime}+1/2+j^{\prime}}\hat{\lambda}\hat{j^{\prime}}\left\\{\begin{array}[]{ccc}\mathcal{L}_{\rm
3N}^{\prime}&\ell^{\prime}&\lambda\\\
1/2&\Lambda&j^{\prime}\end{array}\right\\}$
$\displaystyle\quad\times\sum_{\mathcal{J}^{\prime}}(-1)^{J_{12}+j^{\prime}+\mathcal{J}^{\prime}}\hat{\mathcal{J}}^{\prime}\hat{\Lambda}\left\\{\begin{array}[]{ccc}J_{12}&j^{\prime}&\mathcal{J}^{\prime}\\\
\mathcal{L}_{\rm 3N}^{\prime}&\mathcal{J}_{\rm
tot}&\Lambda\end{array}\right\\}$ $\displaystyle\quad\times|\mathcal{N}_{\rm
3N}^{\prime}N_{12}n^{\prime}[\mathcal{L}_{\rm
3N}^{\prime}(J_{12}j^{\prime})\mathcal{J}^{\prime}]\mathcal{J}_{\rm
tot}\rangle.$

Combining Eqs. (83), (85), and (86), and using the diagonality of the states,
we obtain

$\displaystyle\langle\mathcal{N}_{\rm 3N}$ $\displaystyle Nn\alpha_{\rm
3N}\alpha\mathcal{J}_{\rm tot}|\mathcal{N}_{\rm
3N}^{\prime}N_{12}n^{\prime}[\mathcal{L}_{\rm
3N}^{\prime}(J_{12}j^{\prime})\mathcal{J}^{\prime}]\mathcal{J}_{\rm
tot}\rangle$ (87) $\displaystyle=\delta_{\mathcal{N}_{\rm 3N}\mathcal{N}_{\rm
3N}^{\prime}}\delta_{\mathcal{L}_{\rm 3N}\mathcal{L}_{\rm
3N}^{\prime}}\delta_{NN_{12}}\delta_{LL_{12}}\delta_{SS_{12}}\delta_{JJ_{12}}\delta_{nn^{\prime}}\delta_{\ell\ell^{\prime}}\delta_{jj^{\prime}}\delta_{\mathcal{J}\mathcal{J}^{\prime}},$

and one can find the overlap:

$\displaystyle\langle N_{\rm NN}$ $\displaystyle
N_{12}n_{3}\alpha_{12}\alpha_{3}\mathcal{J}_{\rm tot}|\mathcal{N}_{\rm
3N}Nn\alpha_{\rm 3N}\alpha\mathcal{J}_{\rm tot}\rangle$ (88)
$\displaystyle=(-1)^{J_{12}+L_{12}-J_{\rm
rc}+\ell_{3}+\ell+\mathcal{J}+\mathcal{J}_{\rm tot}}$
$\displaystyle\quad\times\sum_{\Lambda\lambda}(-1)^{j_{3}+\Lambda}\hat{J}_{\rm
rc}\hat{\mathcal{J}}\hat{j}_{3}\hat{j}\hat{\Lambda}^{2}\hat{\lambda}^{2}$
$\displaystyle\quad\times\left\\{\begin{array}[]{ccc}J_{12}&L_{\rm NN}&J_{\rm
rc}\\\ j_{3}&\mathcal{J}_{\rm
tot}&\Lambda\end{array}\right\\}\left\\{\begin{array}[]{ccc}J_{12}&j&\mathcal{J}\\\
\mathcal{L}_{\rm 3N}&\mathcal{J}_{\rm tot}&\Lambda\end{array}\right\\}$
$\displaystyle\quad\times\left\\{\begin{array}[]{ccc}L_{\rm
NN}&\ell_{3}&\lambda\\\
1/2&\Lambda&j_{3}\end{array}\right\\}\left\\{\begin{array}[]{ccc}\mathcal{L}_{\rm
3N}&\ell&\lambda\\\ 1/2&\Lambda&j\end{array}\right\\}$
$\displaystyle\quad\times\langle\mathcal{N}_{\rm 3N}^{\prime}\mathcal{L}_{\rm
3N}^{\prime},n^{\prime}\ell^{\prime}:\lambda|N_{\rm NN}L_{\rm
NN},n_{3}\ell_{3}:\lambda\rangle_{d=2}$
$\displaystyle\quad\times\delta_{NN_{12}}\delta_{LL_{12}}\delta_{SS_{12}}\delta_{JJ_{12}}.$

Inserting $\mathcal{N}_{\rm 3N}=\mathcal{L}_{\rm 3N}=0$ into Eq. (88) leads to
Eq. (76).

##  Appendix C Transformation to single-particle basis

For applications to medium-mass and heavier systems, current developments will
need to be combined with basis-expansion 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 center-of-mass coordinates. They are
related by Talmi-Moshinsky transformations, already mentioned in Appendex B.

Our goal here is to show the transformation for

$\langle
n^{\prime}_{1}\ell^{\prime}_{1}j^{\prime}_{1},n^{\prime}_{2}\ell^{\prime}_{2}j^{\prime}_{2}:J_{\rm
tot}^{\prime}||\mathcal{O}^{\lambda}||n_{1}\ell_{1}j_{1},n_{2}\ell_{2}j_{2}:J_{\rm
tot}\rangle.$ (89)

Here the state $|n_{1}\ell_{1}j_{1},n_{2}\ell_{2}j_{2}:J_{\rm tot}\rangle$ is
the antisymmetrized product of $\\{n_{1},\ell_{1},j_{1}\\}$ and
$\\{n_{2},\ell_{2},j_{2}\\}$ states coupling to the total angular momentum
$J_{\rm tot}$ in the proton-neutron ($p$-$n$) basis. One can find the
following transformation

$\displaystyle|n_{1}\ell_{1}j_{1},$ $\displaystyle n_{2}\ell_{2}j_{2}:J_{\rm
tot}\rangle=\sum_{N_{12}N_{\rm NN}\alpha_{12}}|N_{\rm
NN}N_{12}\alpha_{12}\rangle$ (90) $\displaystyle\times\langle N_{\rm
NN}N_{12}\alpha_{12}|n_{1}\ell_{1}j_{1},n_{2}\ell_{2}j_{2}:J_{\rm
tot}\rangle,$

where the overlap is

$\displaystyle\langle N_{\rm NN}N_{12}$
$\displaystyle\alpha_{12}|n_{1}\ell_{1}j_{1},n_{2}\ell_{2}j_{2}:J_{\rm
tot}\rangle=(-1)^{L+L_{\rm NN}+S_{12}+J_{\rm tot}}f_{12}$ (91)
$\displaystyle\times\hat{j}_{1}\hat{j}_{2}\hat{S}_{12}\hat{J}_{12}\sum_{\Lambda}\hat{\Lambda}^{2}\left\\{\begin{array}[]{ccc}\ell_{1}&s_{1}&j_{1}\\\
\ell_{2}&s_{2}&j_{2}\\\ \Lambda&S_{12}&J_{\rm tot}\end{array}\right\\}$
$\displaystyle\times\left\\{\begin{array}[]{ccc}S_{12}&L_{12}&J_{12}\\\ L_{\rm
NN}&J_{\rm tot}&\Lambda\end{array}\right\\}\delta_{J_{\rm rc}J_{\rm tot}}$
$\displaystyle\times\langle N_{\rm NN}L_{\rm
NN},N_{12}L_{12}:\Lambda|n_{1}\ell_{1},n_{2}\ell_{2}:\Lambda\rangle_{d=1},$

with antisymmetrization factor

$f_{12}=\left\\{\\!\\!\\!\begin{array}[]{cc}\sqrt{\frac{1}{2(1+\delta_{n_{1}n_{2}}\delta_{\ell_{1}\ell_{2}}\delta_{j_{1}j_{2}})}}[1+(-1)^{L_{12}+S_{12}}]&pp\textrm{
or }nn\\\ 1&pn\end{array}\right..$ (92)

Note that we have removed the isospin part from $\alpha_{12}$ as we work in
the proton-neutron formalism. Then, the matrix elements in the laboratory
frame can be computed as

$\displaystyle\langle$ $\displaystyle
n_{1}^{\prime}\ell_{1}^{\prime}j_{1}^{\prime},n_{2}^{\prime}\ell_{2}^{\prime}j_{2}^{\prime}:J_{\rm
tot}^{\prime}||O^{\lambda}||n_{1}\ell_{1}j_{1},n_{2}\ell_{2}j_{2}:J_{\rm
tot}\rangle$ (93) $\displaystyle=\sum_{N_{12}^{\prime}N_{\rm
NN}^{\prime}\alpha_{12}^{\prime}}\sum_{N_{12}N_{\rm NN}\alpha_{12}}\langle
n_{1}^{\prime}\ell_{1}^{\prime}j_{1}^{\prime},n_{2}^{\prime}\ell_{2}^{\prime}j_{2}^{\prime}:J_{\rm
tot}^{\prime}|N_{\rm NN}^{\prime}N_{12}^{\prime}\alpha_{12}^{\prime}\rangle$
$\displaystyle\quad\times\langle N_{\rm
NN}^{\prime}N_{12}^{\prime}\alpha_{12}^{\prime}||O^{\lambda}||N_{\rm
NN}N_{12}\alpha_{12}\rangle$ $\displaystyle\quad\times\langle N_{\rm
NN}N_{12}\alpha_{12}|n_{1}\ell_{1}j_{1},n_{2}\ell_{2}j_{2}:J_{\rm
tot}\rangle.$

If the operator does not depend on the center-of-mass coordinate, the matrix
element can be evaluated as

$\displaystyle\langle N_{\rm NN}^{\prime}N_{12}^{\prime}$
$\displaystyle\alpha_{12}^{\prime}||O^{\lambda}||N_{\rm
NN}N_{12}\alpha_{12}\rangle=(-1)^{L_{\rm NN}+J_{\rm
rc}+J^{\prime}_{12}+\lambda}$ (94) $\displaystyle\times\hat{J}_{\rm
rc}^{\prime}\hat{J}_{\rm rc}\left\\{\begin{array}[]{ccc}J_{\rm
rc}^{\prime}&J_{12}^{\prime}&L_{\rm NN}\\\ J_{12}&J_{\rm
rc}&\lambda\end{array}\right\\}\langle N^{\prime}_{12}\alpha_{\rm
2b}^{\prime}||O^{\lambda}||N_{12}\alpha_{\rm 2b}\rangle$
$\displaystyle\times\delta_{N_{\rm NN}N_{\rm NN}^{\prime}}\delta_{L_{\rm
NN}L_{\rm NN}^{\prime}}.$

## References

  * Epelbaum  _et al._ [2009] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, Rev. Mod. Phys. 81, 1773 (2009).
  * Machleidt and Entem [2011] R. Machleidt and D. R. Entem, Phys. Rep. 503, 1 (2011).
  * Hammer  _et al._ [2020] H.-W. Hammer, S. König, and U. van Kolck, Rev. Mod. Phys. 92, 025004 (2020).
  * Epelbaum  _et al._ [2015] E. Epelbaum, H. Krebs, and U.-G. Meißner, Eur. Phys. J. A 51, 53 (2015).
  * Furnstahl  _et al._ [2015] R. J. Furnstahl, N. Klco, D. R. Phillips, and S. Wesolowski, Phys. Rev. C 92, 024005 (2015).
  * Pastore  _et al._ [2008] S. Pastore, R. Schiavilla, and J. L. Goity, Phys. Rev. C 78, 064002 (2008).
  * Pastore  _et al._ [2009] S. Pastore, L. Girlanda, R. Schiavilla, M. Viviani, and R. B. Wiringa, Phys. Rev. C 80, 034004 (2009).
  * Pastore  _et al._ [2011] S. Pastore, L. Girlanda, R. Schiavilla, and M. Viviani, Phys. Rev. C 84, 024001 (2011).
  * Kölling  _et al._ [2009] S. Kölling, E. Epelbaum, H. Krebs, and U.-G. Meißner, Phys. Rev. C 80, 045502 (2009).
  * Kölling  _et al._ [2011] S. Kölling, E. Epelbaum, H. Krebs, and U.-G. Meißner, Phys. Rev. C 84, 054008 (2011).
  * Krebs [2020] H. Krebs, Eur. Phys. J. A 56, 234 (2020).
  * Suhonen [2007] J. Suhonen,  _Anatol. Stud._, Theoretical and Mathematical Physics, Vol. 23 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007) pp. 191–196.
  * de Forest Jr. and Walecka [1966] T. de Forest Jr. and J. D. Walecka, Adv. Phys. 15, 1 (1966).
  * Donnelly and Walecka [1975] T. W. Donnelly and J. D. Walecka, Annu. Rev. Nucl. Sci. 25, 329 (1975).
  * Phillips and Cohen [2000] D. R. Phillips and T. D. Cohen, Nucl. Phys. A 668, 45 (2000).
  * Walzl and Meißner [2001] M. Walzl and U.-G. Meißner, Phys. Lett. B 513, 37 (2001).
  * Phillips [2003] D. R. Phillips, Phys. Lett. B 567, 12 (2003).
  * Phillips [2007] D. R. Phillips, J. Phys. G 34, 365 (2007).
  * Pavón Valderrama  _et al._ [2008] M. Pavón Valderrama, A. Nogga, E. Ruiz Arriola, and D. R. Phillips, Eur. Phys. J. A 36, 315 (2008).
  * Kölling  _et al._ [2012] S. Kölling, E. Epelbaum, and D. R. Phillips, Phys. Rev. C 86, 047001 (2012).
  * Piarulli  _et al._ [2013] M. Piarulli, L. Girlanda, L. E. Marcucci, S. Pastore, R. Schiavilla, and M. Viviani, Phys. Rev. C 87, 014006 (2013).
  * Epelbaum  _et al._ [2014] E. Epelbaum, A. M. Gasparyan, J. Gegelia, and M. R. Schindler, Eur. Phys. J. A 50, 51 (2014).
  * Nevo Dinur  _et al._ [2019] N. Nevo Dinur, O. J. Hernandez, S. Bacca, N. Barnea, C. Ji, S. Pastore, M. Piarulli, and R. B. Wiringa, Physical Review C 99, 034004 (2019).
  * Schiavilla  _et al._ [2019] R. Schiavilla, A. Baroni, S. Pastore, M. Piarulli, L. Girlanda, A. Kievsky, A. Lovato, L. E. Marcucci, S. C. Pieper, M. Viviani, and R. B. Wiringa, Phys. Rev. C 99, 034005 (2019).
  * Filin  _et al._ [2020] A. A. Filin, V. Baru, E. Epelbaum, H. Krebs, D. Möller, and P. Reinert, Phys. Rev. Lett. 124, 082501 (2020).
  * Filin  _et al._ [2021] A. A. Filin, D. Möller, V. Baru, E. Epelbaum, H. Krebs, and P. Reinert, Phys. Rev. C 103, 024313 (2021).
  * Pastore  _et al._ [2013] S. Pastore, S. C. Pieper, R. Schiavilla, and R. B. Wiringa, Phys. Rev. C 87, 035503 (2013).
  * Friman-Gayer  _et al._ [2021] U. Friman-Gayer, C. Romig, T. Hüther, K. Albe, S. Bacca, T. Beck, M. Berger, J. Birkhan, K. Hebeler, O. J. Hernandez, I. Isaak, S. König, N. Pietralla, P. C. Ries, J. Rohrer, R. Roth, D. Savran, M. Scheck, A. Schwenk, R. Seutin, and V. Werner, Phys. Rev. Lett. 126, 102501 (2021).
  * Gysbers  _et al._ [2019] P. Gysbers, G. Hagen, J. D. Holt, G. R. Jansen, T. D. Morris, P. Navrátil, T. Papenbrock, S. Quaglioni, A. Schwenk, S. R. Stroberg, and K. A. Wendt, Nature Phys. 15, 428 (2019).
  * Hagen  _et al._ [2014] G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean, Reports Prog. Phys. 77, 096302 (2014), arXiv:1312.7872 .
  * Hergert  _et al._ [2016] H. Hergert, S. Bogner, T. Morris, A. Schwenk, and K. Tsukiyama, Phys. Rep. 621, 165 (2016), arXiv:1512.06956 .
  * Barbieri and Carbone [2017] C. Barbieri and A. Carbone, Self-Consistent Green’s Function Approaches, in  _An Advanced Course in Computational Nuclear Physics_, edited by M. Hjorth-Jensen, M. P. Lombardo, and U. van Kolck (Springer International Publishing, Cham, 2017) p. 571.
  * Stroberg  _et al._ [2019] S. R. Stroberg, H. Hergert, S. K. Bogner, and J. D. Holt, Annu. Rev. Nucl. Part. Sci. 69, 307 (2019), arXiv:1902.06154 .
  * Tichai  _et al._ [2020] A. Tichai, R. Roth, and T. Duguet, Front. Phys. 8, 164 (2020).
  * Ye  _et al._ [2018] Z. Ye, J. Arrington, R. J. Hill, and G. Lee, Physics Letters B 777, 8 (2018).
  * [36] S. Pal, S. Sarker, P. J. Fasano, P. Maris, J. P. Vary, and M. A. Caprio,  arXiv:2304.01389 .
  * Blatt and Weisskopf [1979] J. M. Blatt and V. F. Weisskopf,  _Theoretical Nuclear Physics_ (Springer, New York, 1979).
  * Sachs [1948] R. G. Sachs, Phys. Rev. 74, 433 (1948).
  * Hebeler [2021] K. Hebeler, Phys. Rept. 890, 1 (2021), arXiv:2002.09548 [nucl-th] .
  * Golak  _et al._ [2005] J. Golak, R. Skibiński, H. Witała, W. GlÃ¶ckle, A. Nogga, and H. Kamada, Physics Reports 415, 89 (2005).
  * Navrátil and Barrett [1998] P. Navrátil and B. R. Barrett, Phys. Rev. C 57, 562 (1998), arXiv:nucl-th/9711027v1 [arXiv:nucl-th] .
  * Barrett  _et al._ [2013] B. R. Barrett, P. Navrátil, and J. P. Vary, Prog. Part. Nucl. Phys. 69, 131 (2013).
  * Caprio  _et al._ [2012] M. A. Caprio, P. Maris, and J. P. Vary, Phys. Rev. C 86, 034312 (2012).
  * Entem  _et al._ [2020] D. R. Entem, R. Machleidt, and Y. Nosyk, Front. Phys. 8, 57 (2020).
  * Drischler  _et al._ [2019] C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 122, 042501 (2019).
  * Seutin [2021] R. Seutin,  _Electromagnetic Properties of Light Nuclei with Chiral Effective Field Theory Currents_ , Ph.D. thesis, TU Darmstadt (2021).
  * Entem and Machleidt [2003] D. R. Entem and R. Machleidt, Phys. Rev. C 68, 041001 (2003).
  * Melendez  _et al._ [2017] J. A. Melendez, S. Wesolowski, and R. J. Furnstahl, Phys. Rev. C 96, 024003 (2017).
  * Melendez  _et al._ [2019] J. A. Melendez, R. J. Furnstahl, D. R. Phillips, M. T. Pratola, and S. Wesolowski, Phys. Rev. C 100, 044001 (2019).
  * Phillips [2016] D. R. Phillips, Annu. Rev. Nucl. Part. Sci. 66, 421 (2016).
  * Amroun  _et al._ [1994] A. Amroun, V. Breton, J.-M. Cavedon, B. Frois, D. Goutte, F. P. Juster, P. Leconte, J. Martino, Y. Mizuno, X.-H. Phan, S. K. Platchkov, I. Sick, and S. Williamson, Nucl. Phys. A 579, 596 (1994).
  * Mohr  _et al._ [2021] P. J. Mohr, B. N. Taylor, and D. B. Newell, Rev. Mod. Phys. 93, 025010 (2021).
  * Carlson and Schiavilla [1998] J. Carlson and R. Schiavilla, Rev. Mod. Phys. 70, 743 (1998).
  * Miyagi [2023] T. Miyagi, Eur. Phys. J. A 59, 150 (2023).
  * Talmi [1952] I. Talmi, Helv. Phys. Acta 25, 185 (1952).
  * Moshinsky [1959] M. Moshinsky, Nucl. Phys. 13, 104 (1959).
  * Kamuntavicius  _et al._ [2001] G. Kamuntavicius, R. Kalinauskas, B. Barrett, S. Mickevicius, and D. Germanas, Nucl. Phys. A 695, 191 (2001).