Dark matter scattering off 4He in
chiral effective field theory

J. de Vriesa,b, C. Körberc, A. Noggad, S. Shaina,b
a Institute for Theoretical Physics Amsterdam and Delta Institute for
Theoretical Physics, University of Amsterdam, Science Park 904,
1098 XH Amsterdam, The Netherlands
b Nikhef, Theory Group, Science Park 105, 1098 XG, Amsterdam, The Netherlands

c Institut für Theoretische Physik II, Ruhr-Universität Bochum,
D-44780 Bochum, Germany
d Institut für Kernphysik (IKP-3), Institute for Advanced Simulation (IAS-4),
Center for Advanced Simulation and Analytics (CASA),
and Jülich Center for Hadron Physics (JCHP),
Forschungszentrum Jülich, D-52425 Jülich, Germany

We study dark matter scattering off 4He and other light nuclei using chiral
effective field theory. We consider scalar DM interactions and include both
one- and two-nucleon scattering processes. The DM interactions and nuclear
wave functions are obtained from chiral effective field theory and we work up
to fourth order in the chiral expansion for the latter to investigate the
chiral convergence. The results for the scattering rates can be used to
determine the sensitivity of planned experiments to detect relatively light
dark matter particles using 4He. We find that next-to-leading-order scalar
currents are smaller than expected from power counting for scattering off 4He
confirming earlier work. However, the results for two-nucleon corrections
exhibit a linear regulator dependence indicating potential problems in the
applied power counting. We observe a linear correlation between the, in
principle not observable, D-wave probability of various light nuclei and the
scalar two-nucleon matrix elements, again pointing towards potentially missing
contributions.

##  1 Introduction

While no direct evidence for dark matter (DM) has been found in a laboratory
experiment, there exists strong observational evidence for its existence [1].
Next-generation direct detection experiments aim to probe uncharted parameter
space for a wide range of DM masses, going beyond the ‘traditional’ WIMP
(Weakly Interacting Massive Particles) regime of GeV-to-TeV DM masses [2, 3,
4]. To interpret these direct detection searches in terms of the underlying DM
models and to connect to cosmological aspects of DM, such as the relic
density, it is important to have well-controlled theoretical predictions for
WIMP cross-sections off atomic nuclei. These predictions are complicated due
to the presence of widely separated energy scales associated with particle,
hadronic, and nuclear processes. The last decade has seen the development of
effective field theory (EFT) approaches to overcome this difficulty [5, 6, 7,
8, 9, 10, 11]. For example, by assuming that DM fields are singlets under the
Standard Model (SM) gauge symmetries, the DM interactions with SM fields can
be captured by a series of effective operators which are dominated by the
operators of the lowest dimension. After renormalization-group evolution to
lower energies [12, 13], these interactions can be matched to effective
interactions between DM and hadrons and nuclei, which, in turn, can be used to
compute DM-nucleus scattering rates.

In this work, we apply the framework of chiral EFT, the low-energy EFT of QCD,
to perform the last steps of this program. We consider scalar-mediated DM-SM
interactions, which we match onto chiral EFT hadronic interactions. The chiral
EFT power counting, often called Weinberg counting, predicts that DM-nucleon
interactions (or currents) provide the leading-order (LO) contribution,
whereas two-nucleon currents appear at next-to-leading-order (NLO) [14, 7, 8,
9, 15]. This is potentially interesting as these NLO two-body currents are
predicted to be relatively large and could potentially be used to unravel the
underlying DM interactions from DM scattering off various nuclear isotopes.

Unfortunately, explicit computations of the two-body corrections are found to
be inconclusive. While shell-model computations for 132Xe found
$\mathcal{O}(10\%$-$20\%)$ corrections [16], of the expected size,
calculations on lighter nuclei found smaller effects [17, 15, 18]. More
problematic is that two-body corrections in DM scattering off the deuteron,
3H, and 3He, strongly depend on the details of the applied wave function and
thus on the applied nucleon-nucleon potential used to generate the wave
functions [15]. Similar conclusions were drawn in Ref. [18], which applied
phenomenological wave functions of various light nuclei and observed a large
dependence on the applied regulator used in intermediate steps of the quantum
Monte Carlo calculations. These results are worrisome as they indicate a
potential problem in the chiral EFT power counting for two-nucleon currents
and potentially jeopardize the interpretation of WIMP-nucleus scattering.
Similar problems were recently identified in other nuclear probes of beyond-
the-Standard Model physics, such as neutrinoless double beta decay [19] and
indicate the presence of additional short-range two-nucleon currents.

Our main purpose here is to extend the calculations of Ref. [15] to scattering
off 4He nuclei for two reasons. First of all, experiments have been proposed
using a liquid 4He detector [3, 20]. The main motivation to use relatively
light target nuclei is to get sensitivity to lighter WIMPs to which more
conventional experiments, for instance those involving Xe nuclei, have less
sensitivity. Second, compared to $A=2$ and $A=3$ nuclei, the binding energy
per nucleon is much larger and comparable to heavier isotopes. We are
particularly interested in determining whether the conclusions of [15]
regarding two-nucleon currents are related to the diluteness of the deuteron
and 3H and 3He nuclei.

In Section 2, our computational framework is described. Specifically, we
introduce the relation of the DM scattering cross sections to a set of
response functions which are calculated in following based on chiral nuclear
wave functions. The calculation of the wave functions is summarized in Section
3. In this section, we also discuss our Bayesian approach to quantify
uncertainties due to the limiting the wave function calculations to finite
order. Finally, we briefly introduce nuclear transition densities that connect
the wave function to the application to DM currents. Our results are given in
Section 4. Special attention is given here to the dependence on the regulator
applied in the calculations. We conclude in Section 5. Details of the
definition of the densities and the partial wave decomposition of the current
are deferred to the appendix.

##  2 Computational framework

(a)

(b)

(c)

Figure 1: Diagrams contribution to DM-nucleus scattering. Solid lines
correspond to nucleons, single dashed lines to pions and double-dashed lines
to DM.

In this work, we follow up the computations of Ref. [15] and perform first-
principle computations of DM scattering off 4He isotopes. We apply chiral EFT
to both generate the 4He wave functions as well as the DM currents and
systematically compute the resulting DM-4He scattering rate. We select 4He
because of two reasons as mentioned above: it is more sensitive for light WIMP
searches, see e.g. [21, 20] and the binding energy per nucleon is much larger
than for $A<4$ nuclei and comparable to heavier isotopes. We are particularly
interested to determine whether the conclusions of [15] regarding two-nucleon
currents are related to specific spin-isospin properties and/or the diluteness
of the deuteron and 3H and 3He nuclei. Much of the calculation framework
follows that of Ref. [15] (see also Ref. [18] which uses a very similar setup)
and here we only give the main ingredients. Our starting point are scalar DM-
SM interactions of the form

$\mathcal{L}_{\chi}=\bar{\chi}\chi\left(c_{u}\,m_{u}\,\bar{u}u+c_{d}\,m_{d}\,\bar{d}d+c_{s}\,m_{s}\,\bar{s}s+c_{G}\,\alpha_{s}G_{\mu\nu}^{a}G^{\mu\nu\,a}\right)\
,$ (1)

where $\chi$ denotes a spin-$\tfrac{1}{2}$ DM fermion (for other DM spins, the
computations are almost identical). $u$, $d$, and $s$ denote, respectively,
quark fields with quark masses $m_{u,d,s}$, and $G_{\mu\nu}^{a}$ is the gluon
field strength. We factored out one power of $\alpha_{s}=g_{s}^{2}/(4\pi)$.
$c_{u,d,s,G}$ describe three unknown coupling constants of mass dimension
$(-3)$ that parametrize the couplings strengths of DM with quarks and gluons.
They can be computed in specific DM models, for example in Higgs portal models
[22]. We consider the Lagrangian in Eq. (1) to be valid at relatively low
energies ($\mu=1$ GeV) where we match to hadronic DM interactions. Other
interactions beyond those in Eq. (1) are certainly possible, but tend to lead
to suppressed two-body currents. See for instance Ref. [23] for recent first-
principle computations for spin-dependent cross sections.

By application of chiral perturbation theory, the interactions in Eq. (1) can
be matched to interactions between DM and nucleons, pions, and heavier
hadrons. The leading-order (LO) currents involve a single nucleon and are
conveniently written as

$\displaystyle J^{(\mathrm{one-body})}(\boldsymbol{q})$ $\displaystyle=$
$\displaystyle\left[\sigma_{\pi N}-\frac{9g_{A}^{2}\pi m_{\pi}^{3}}{4(4\pi
f_{\pi})^{2}}F\left(\frac{|\boldsymbol{q}|}{2m_{\pi}}\right)\right]\,\bar{c}_{q^{\mathrm{(is)}}}-\frac{\delta
m_{N}}{4}\bar{c}_{q^{\mathrm{(iv)}}}\,\tau^{3}_{i}$ (2)
$\displaystyle+c_{s}\sigma_{s}-c_{G}\frac{8\pi}{9}m_{N}^{G}\,.$

Here we have introduced the isoscalar and isovector combinations of couplings

$\bar{c}_{q^{\mathrm{(is)}}}=\frac{1}{2}\left[c_{u}(1-\varepsilon)+c_{d}(1+\varepsilon)\right]\,,\qquad\bar{c}_{q^{\mathrm{(iv)}}}=\left[c_{u}(1-\frac{1}{\varepsilon})+c_{d}(1+\frac{1}{\varepsilon})\right]\,,$
(3)

in terms of $\varepsilon=(m_{d}-m_{u})/(m_{d}+m_{u})=0.36\pm 0.03$[24], and
the loop function associated to Diagram 1(b)

$F(x)=\frac{-x+(1+2x^{2})\arctan x}{3x}\,,$ (4)

and $F(x\ll 1)\simeq 5x^{2}/9+\dots$. The various low-energy constants are
given by [25]

$\sigma_{\pi N}=(59.1\pm 3.5)\,\mathrm{MeV}\ ,\qquad\delta m_{N}=(2.32\pm
0.17)\,\mathrm{MeV}\,,$ (5)

where we used a Roy-Steiner extraction of the pion nucleon sigma term [26].
Lattice QCD tends to predict somewhat smaller values but might be plagued by
excited-state contamination [27] and must be corrected for isospin breaking
[28]. For the strange matrix elements we use [29]

$\sigma_{s}=m_{s}\frac{dm_{N}}{dm_{s}}=(37\pm 3)\mbox{ MeV}\,.$ (6)

We defined $m_{N}^{G}\equiv m_{N}-\sigma_{\pi N}-\sigma_{s}$ as the gluonic
contribution to the nucleon mass.

At next-to-leading-order (NLO) Diagram 1(c) induces a two-body current given
by

$J^{\mathrm{two-
body}}(\boldsymbol{q})=-m_{\pi}^{2}\left(\frac{g_{A}}{2f_{\pi}}\right)^{2}\frac{(\boldsymbol{\sigma}_{1}\cdot\boldsymbol{q}_{1})(\boldsymbol{\sigma}_{2}\cdot\boldsymbol{q}_{2})}{(\boldsymbol{q}_{1}^{\,2}+m_{\pi}^{2})(\boldsymbol{q}_{2}^{\,2}+m_{\pi}^{2})}\tau_{1}\cdot\tau_{2}\,\bar{c}_{q^{\mathrm{(is)}}}\
,$ (7)

where $\boldsymbol{q}_{i}=\boldsymbol{p}^{\,\prime}_{i}-\boldsymbol{p}_{i}$ is
the difference between the outgoing and incoming momentum of nucleon $i$, and
$\sigma_{i}$ ($\tau_{i}$) the spin (isospin) of nucleon $i$.

At next order in the expansion (N2LO), several new diagrams appears, but they
can all be absorbed into the LO one-body currents and are effectively absorbed
in the values of the low-energy constants used in Eq. (2). At N3LO there
appear new two-nucleon currents both from tree- and loop-level diagrams [30]
but also two-nucleon contact interactions with presently undetermined low-
energy constants. One outstanding issue is how to consistently regularize the
loop diagrams entering the current and the loop diagrams arising from
iterating the nucleon-nucleon potential that lead to the nuclear wave
functions [31]. For these reasons we stick to the scalar currents up to N2LO
which only involve tree-level diagrams and no unknown low-energy constants.

###  2.1 Scattering cross section

We investigate scattering processes
$\chi(\boldsymbol{p}_{\chi})+T(\boldsymbol{p}_{T})\rightarrow\chi(\boldsymbol{p}_{\chi}^{\prime})+T(\boldsymbol{p}_{T}^{\prime})$,
where $T$ denotes the target nucleus with mass $m_{T}$ consisting of $A$
nucleons. The elastic unpolarized differential cross section is given by

$\frac{\text{d}{\sigma}\,}{\text{d}{\boldsymbol{q}^{2}}\,}=\frac{1}{4\pi
v_{\chi}^{2}}\frac{1}{2j+1}\sum_{m_{j},m_{j}^{\prime}=-j}^{j}\bigg{|}\left\langle\Psi_{T},jm_{j}^{\prime}|\
\hat{J}(\boldsymbol{q}^{2})\ |\Psi_{T},jm_{j}\right\rangle\bigg{|}^{2}\ ,$ (8)

where $\boldsymbol{q}$ is the momentum transfer from DM to the target nucleus,
and $v_{\chi}$ the DM velocity. The wave function of the target nucleus
$\ket{\Psi_{T},jm_{j}}$ corresponds to a nucleus with total angular momentum
$j$ and polarization $m_{j}$.

To discuss the various one- and two-body nuclear response functions it is
convenient to factor out the isoscalar piece. We follow Ref. [15] and classify
the cross section in terms of response functions,
$\mathcal{F}^{(\nu)}_{i,\,a}(\boldsymbol{q}^{\,2})$, that depend on

  * •

$i=\\{\mathrm{is},\,\mathrm{iv},\,s,\,G\\}$ which label the dependence on the
DM-quark and DM-gluon couplings.

  * •

The chiral order $\nu=0,\ 1,\dots$ of the current where the dominant current
for every DM interaction
($\\{\bar{c}_{q^{\mathrm{(is)}}},\,\bar{c}_{q^{\mathrm{(iv)}}},\,c_{s},\,c_{G}\\}$)
begins at order $0$.

  * •

We divide the NLO contributions into two-nucleon, $a=2b$, and radius
corrections, $a=r$.

We then decompose

$\displaystyle\frac{\text{d}{\sigma}\,}{\text{d}{\boldsymbol{q}^{2}}\,}=\bar{c}_{q^{\mathrm{(is)}}}^{2}\frac{\sigma_{\pi
N}^{2}A^{2}}{4\pi v_{\chi}^{2}}\bigg{|}$
$\displaystyle\left(\mathcal{F}^{\left(0\right)}_{\mathrm{is}}{\left(\boldsymbol{q}^{2}\right)}+\mathcal{F}^{\left(1\right)}_{\mathrm{is},\,2b}\left(\boldsymbol{q}^{2}\right)+\mathcal{F}^{\left(1\right)}_{\mathrm{is},\,r}\left(\boldsymbol{q}^{2}\right)+\dots\right)$
$\displaystyle+\alpha_{\mathrm{iv}}\left(\mathcal{F}^{\left(0\right)}_{\mathrm{iv}}{\left(\boldsymbol{q}^{2}\right)}+\dots\right)+\alpha_{s}\left(\mathcal{F}^{\left(0\right)}_{s}{\left(\boldsymbol{q}^{2}\right)}+\dots\right)$
$\displaystyle+\alpha_{G}\left(\mathcal{F}^{\left(0\right)}_{G}{\left(\boldsymbol{q}^{2}\right)}+\dots\right)\bigg{|}^{2}\,,$
(9)

where we kept terms up to NLO in the chiral power counting, and

$\alpha_{\mathrm{iv}}=-\left(\frac{\delta m_{N}}{4\sigma_{\pi
N}}\right)\frac{\bar{c}_{q^{\mathrm{(iv)}}}}{\bar{c}_{q^{\mathrm{(is)}}}},\qquad\alpha_{s}=\left(\frac{\sigma_{s}}{\sigma_{\pi
N}}\right)\frac{c_{s}}{\bar{c}_{q^{\mathrm{(is)}}}},\qquad\alpha_{G}=\left(-\frac{8\pi}{9}\frac{m^{G}_{N}}{\sigma_{\pi
N}}\right)\frac{c_{G}}{\bar{c}_{q^{\mathrm{(is)}}}}\ .$ (10)

Many response functions are actually equal up to this order in the power
counting and we have

$\displaystyle\mathcal{F}^{\left(0\right)}_{\mathrm{is}}(\boldsymbol{q}^{2})$
$\displaystyle=$
$\displaystyle\mathcal{F}^{\left(0\right)}_{s}(\boldsymbol{q}^{2})\,\,=\,\,\mathcal{F}^{\left(0\right)}_{G}(\boldsymbol{q}^{2})\,,$
(11)

which are normalized to $\mathcal{F}^{\left(0\right)}_{\mathrm{is}}(0)=1$. The
radius correction $\mathcal{F}^{\left(1\right)}_{\mathrm{is},\
r}\left(\boldsymbol{q}^{2}\right)$ also involves the same nuclear information
as $\mathcal{F}^{\left(0\right)}_{\mathrm{is}}$ except for an additional
overall dependence on $\boldsymbol{q}^{2}$. We therefore need to perform
nuclear calculations of
$\mathcal{F}^{\left(0\right)}_{\mathrm{is}}(\boldsymbol{q}^{2})$,
$\mathcal{F}^{\left(0\right)}_{\mathrm{iv}}(\boldsymbol{q}^{2})$, and
$\mathcal{F}^{\left(1\right)}_{\mathrm{is},\
2b}\left(\boldsymbol{q}^{2}\right)$, which we will present below for various
target nuclei.

##  3 Nuclear wave function and scattering matrix elements

###  3.1 Nuclear wave functions

We use the momentum-space basis to calculate the wave function and matrix
elements involving the DM interactions discussed above. 2H, 3He, and 4He wave
functions are evaluated by solving the non-relativistic Schrödinger equation
in momentum-space. Ref. [15] and Ref. [32] shows detailed calculation of 2H
and 3He, and 4He wave functions, respectively. In this section, we summarize
the relevant parts of the wave function calculations.

The wave function of deuteron ($d$) is easily evaluated by solving

$\displaystyle\ket{\psi_{d}}=\frac{1}{E_{d}-T}V_{12}\ket{\psi_{d}}\,,$ (12)

where $E_{d}$ is the deuteron binding energy, $T$ is the two-nucleon ($N\\!N$)
kinetic energy, and $V_{12}$ is the $N\\!N$ potential. It is convenient to use
the partial wave basis ($\ket{p\alpha}$) for the wave function, where $p$ is
the momentum of the deuteron and $\alpha$ represents the partial wave
contributing to the deuteron bound state. The deuteron bound state has a total
orbital angular momentum $l_{12}=0,2$, a spin $s_{12}=1$, and a total angular
momentum $j_{12}=1$.

We use modern phenomenological $N\\!N$ interactions and $\chi$PT to calculate
the nuclear wave functions. We calculate the former by using the standard AV18
[33] and CD-Bonn [34] phenomenological interactions. Since these two
interactions have different non-observable kinetic energies, potential
energies, and D-state probabilities, they provide an added benefit of checking
model dependence [15].

We mainly use semilocal coordinate-space regularized (SCS) chiral interactions
from Ref. [35] to calculate the chiral wave functions up to N4LO chiral order.
For 4He, we also show some results based on the semilocal momentum-space (SMS)
regularized interactions from Ref. [36] which include N5LO contact
interactions in the F-waves, these are labeled N4LO+. The approach
successfully describes various $N\\!N$ observables to a high degree of
accuracy provided that the interactions are regulated using finite value cut-
offs. The SCS interactions use a regularization defined in the configuration
space and parameterized by a short-distance scale $R$. In Ref. [37], an
optimal range of cut-offs $R$ between $0.8\mbox{ fm}$ and $1.2\mbox{ fm}$ has
been identified. For this range of cutoffs, the LECs are of natural size and
no spurious bound states exists. It was also found that the best description
of $N\\!N$ data can be achieved for $R=0.9\mbox{ fm}$. Since an
electromagnetic (EM) interaction accompanies AV18, we included the same EM
interactions for the CD-Bonn and the chiral potentials except for the neutron-
proton system since the relevant parts are here already included in the strong
part.

The calculation of the 3He and 3H wave functions was summarized in Ref. [15].
The wave functions for 4He ($\ket{\Psi}$) are obtained by rewriting the
Schrödinger equation into a set of two Yakobovsky equations [32]

$\displaystyle\ket{\psi_{1A}}$
$\displaystyle\equiv\ket{\psi_{(12)3,4}}=G_{0}t_{12}P[(1-P_{34})\ket{\psi_{1A}}+\ket{\psi_{2A}}]+(1+G_{0}t_{12})V^{(3)}_{123}\ket{\Psi}\,,$
(13) $\displaystyle\ket{\psi_{2A}}$
$\displaystyle\equiv\ket{\psi_{(12)34}}=G_{0}t_{12}\tilde{P}[(1-P_{34})\ket{\psi_{1A}}+\ket{\psi_{2A}}],$
(14)

where we have introduced two Yakubovsky components $\ket{\psi_{1A}}$ and
$\ket{\psi_{2A}}$. The $N\\!N$ interactions enter into this equation by the
$N\\!N$ T-matrix $t_{12}$ and $G_{0}$ denotes the free 4N propagator. From the
three-nucleon force (3NF) the part $V_{123}^{(3)}$ symmetric in subsystem (12)
is applied. The total 3NF is obtained by cyclic permutations as

$\displaystyle V_{123}=V_{123}^{(1)}+V_{123}^{(2)}+V_{123}^{(3)}.$ (15)

In our previous study, we found that the contribution of the 3NFs to the
matrix elements relevant for dark matter scattering is minor [15]. Therefore,
we will neglect this part in the calculations below. Using the antisymmetry of
the four-nucleon wave functions, it is possible to relate different Yakubovsky
components using the permutation operators
$P=P_{12}P_{23}+P_{13}P_{23},P_{34}$ and $\tilde{P}=P_{14}P_{23}$. The total
4N wave function can then be obtained applying these permutation operators to
the Yakubovsky components

$\displaystyle\ket{\Psi}=[1-(1+P)P_{34}](1+P)\ket{\psi_{1A}}+(1+P)(1+\tilde{P})\ket{\psi_{2A}}$
(16)

The advantage of using Yakubovsky components for these bound state
calculations is that the partial wave decomposition converges faster for the
Yakubovsky components than for the wave function and also because $\psi_{1}$
and $\psi_{2}$ are expanded independently in 3+1 and 2+2 Jacobi momenta,
respectively.

Fig. 8 of the appendix defines these Jacobi coordinates. Explicitly, the basis
states for 3+1 coordinates are given by

$\displaystyle\ket{p_{12}p_{3}p_{4}\alpha}=\ket{p_{12}p_{3}p_{4}\left[\left((l_{12}s_{12})j_{12}(l_{3}\,1/2)I_{3}\right)j_{3}(l_{4}\,1/2)I_{4}\right]J\
\left((t_{12}\,1/2)\tau_{3}1/2\right)T}$ (17)

where $p_{12},p_{3}$ and $q_{4}$ are magnitudes of the Jacobi relative momenta
between particle (12), the third particle and the (12) subsystem and the
fourth particle and the (123) subsystem, respectively. The orbital angular
momenta related to these momenta are denoted by $l_{12},l_{3}$, and $l_{4}$
and are coupled to the spin of the (12) system $s_{12}$, and the
spin-$\tfrac{1}{2}$ of particles 3 and 4 to intermediate angular momenta
$j_{12},I_{3}$, and $I_{4}$. We indicate the coupling scheme using the
brackets in Eqs. (17) and (18). $j_{12}$ and $I_{3}$ are coupled to the total
angular momentum of the (123) system, $j_{3}$, which is when coupled with
$I_{4}$ to the total angular momentum of the 4He state $J=0$. The isospin of
the (12) subsystem $t_{12}$ is coupled with the isospin of the third particle
to $\tau_{3}$ and with the fourth particle to the total isospin of 4He $T=0$.

For 2+2 coordinates the basis states are

$\displaystyle\ket{p_{12}p_{34}q\beta}=\ket{p_{12}p_{34}q\left(\left[(l_{12}s_{12})j_{12}\lambda\right]I\,(l_{34}s_{34})j_{34}\right)J\
(t_{12}\,t_{34})T}$ (18)

where $p_{12},p_{34}$, and $q$ are magnitudes of the Jacobi relative momenta
between particle (12) and (34) and relative momentum of (12) and (34)
subsystems, respectively. The additional orbital angular momenta related to
these momenta are denoted by $l_{34}$ and $\lambda$. $l_{12}$ and $l_{34}$ are
coupled to the spin of the (12) and (34) systems, $s_{12}$ and $s_{34}$, to
subsystem angular momenta $j_{12}$ and $j_{34}$. $j_{12}$ and $\lambda$ are
coupled to the intermediate angular momentum $I$ which is then coupled with
$j_{34}$ to the total angular momentum $J=0$. The isospin of the (12) and (34)
subsystems are coupled to the total isospin $T=0$.

For our calculation, we restrict the orbital angular momentum to 6, $j_{12}$
to 5 and the sum of the orbital angular moment $l_{12}+l_{3}+l_{4}$ and
$l_{12}+l_{34}+\lambda$ to 10. With this restriction, we achieve binding
energy converged up to 20 keV. More details can be found in Ref. [32].

###  3.2 Chiral expansion and uncertainty estimation

On top of the chiral expansion scheme associated with the dark matter current
operators, also the wave functions are obtained by chiral power counting
scheme. The expansion of the wave functions follows chiral effective field
theory[36] for describing the target nucleus in the absence of dark matter.
The parameter for this expansion is associated with a momentum scale at the
order of binding momenta and the pion mass $\sim Q_{\psi}\sim
m_{\pi}/\Lambda_{\chi}\sim\gamma/\Lambda_{\chi}$, where $\Lambda_{\chi}\sim 1$
GeV, is the chiral EFT breakdown scale and $\gamma$ denotes the binding
momentum. The current expansion considers an external scalar source, dark
matter [38], which also has an an expansion in $Q_{\psi}$ but also obtains
corrections scaling as $Q_{J}\sim q/\Lambda_{\chi}$, where $q$ denotes the
momentum transfer. For the currents we will collectively denote the scaling
with $Q_{J}$ to also account for corrections scaling as
$m_{\pi}/\Lambda_{\chi}\sim\gamma/\Lambda_{\chi}$.

The truncation error of matrix elements are dominated by the largest
corrections to either source. That is, performing the computation at a fixed
order of the wave functions $\nu$ and the current operator $\mu$, the matrix
element is expected to scale as

$\displaystyle\ket{\psi^{(\nu)}}$
$\displaystyle\equiv\ket{\psi}[1+\mathcal{O}(Q_{\psi}^{\nu+1})]\,,$ (19)
$\displaystyle\hat{J}^{(\mu)}$
$\displaystyle\equiv\hat{J}[1+\mathcal{O}(Q_{J}^{\mu+1})]\,,$ (20)
$\displaystyle\Rightarrow\braket{\psi^{(\nu)}}{\hat{J}^{(\mu)}}{\psi^{(\nu)}}$
$\displaystyle=\braket{\psi}{\hat{J}}{\psi}[1+\max(O(Q_{J}^{\mu+1}),O(Q_{\psi}^{\nu+1}))]\,.$
(21)

For external currents at a momentum scale close to the nuclear scale, it makes
sense to combine power countings [39] so that the maximum of both corrections
can be related to the minimum of respective orders. As such, one has to
simultaneously increase both orders to make more accurate predictions.

To test the accuracy in the description of the scalar-DM matrix elements, we
test the convergence behavior for DM-currents in dependence of the chiral
order of the wave functions

$|\braket{\psi^{(\nu)}}{\hat{J}^{(\mu_{\rm fixed})}}{\psi^{(\nu)}}|^{2}\equiv
J^{(\nu)}(\mu_{\rm fixed})\,,$ (22)

This allows us to employ standard chiral uncertainty estimation schemes[40]
for the matrix element, where the converged matrix element $O_{J}^{(\infty)}$
is equal to the matrix element at a finite order plus corrections

$J^{(\infty)}(\mu_{\rm fixed})=J^{(\nu)}(\mu_{\rm fixed})+\delta
J^{(\nu)}(\mu_{\rm fixed})\,.$ (23)

This correction, the truncation error, follows the expansion

$J^{(\nu)}(\mu_{\rm fixed})\leftrightarrow J^{\rm
ref}\sum_{n=0}^{\nu}c_{n}Q_{\psi}^{n}\,,\qquad\delta J^{(\nu)}(\mu_{\rm
fixed})\leftrightarrow J^{\rm ref}\sum_{n=\nu+1}^{\infty}c_{n}Q_{\psi}^{n}\,.$
(24)

For a fixed reference scale $J^{\rm ref}$ and expansion parameter $Q_{\psi}$,
knowing the values of all expansion coefficients $c_{n}$ allows to infer the
size of the truncation error. In particular, we apply a Bayesian inference
scheme where we relate the distribution of expansion coefficients to the final
observable. As such the posterior predictive distribution of the truncation
error is obtained by integrating over (marginalizing) the posterior
distribution of the expansion coefficients $P(\\{c_{n}\\}|M,D)$ times the
distribution of the observable given the coefficients (i.e., a delta
distribution)

$\displaystyle P(\delta|\nu;M,D)$ $\displaystyle\sim\left(\prod_{n}\int
dc_{n}\right)\delta(\delta-\delta
J^{(\nu)}(\\{c_{n}\\}))P(\\{c_{n}\\}|M,D)\,.$ (25)

The data points $D$ correspond to observables computed at different chiral
orders $D=\\{J^{(\nu)}(\mu_{\rm fixed})\\}$, times a prior distribution
$P(\\{c_{n}\\}|M)$ for the expansion coefficients

$\displaystyle P(\\{c_{n}\\}|M,D)$
$\displaystyle\sim\exp\left\\{-\chi^{2}(\\{c_{n}\\},D)/2\right\\}P(\\{c_{n}\\}|M)\,.$
(26)

The letter $M$ indicates that your probabilities also depend on the model,
i.e. the ansatz in Eq. (24). The posterior distribution is given by a data
likelihood distribution ($\sim\exp\\{-\chi^{2}/2\\}$) with

$\chi^{2}=\sum_{\nu}\left(J^{\nu}-J^{\rm
ref}\sum_{n=0}^{\nu}c_{n}Q_{\psi}^{n}\right)^{2}\,.$ (27)

The $\chi^{2}$ distribution employed here does not accommodate any additional
uncertainty as the only uncertainty for the data points $D$, in this case, are
numerical errors which we assume to be uncorrelated and of equal size for each
computation. Thus they correspond to an overall normalization factor not
important for this analysis. The employed prior, following the prescription in
the literature, aims a quantifying the naturalness assumption. This estimation
is enabled by the software gsum[41]111With df=0.6, scale=0.8 and the expansion
parameter $Q_{\psi}=(m_{\pi}^{k}+q^{k})/(\Lambda_{b}(m_{\pi}^{k-1}+q^{k-1}))$,
where $q$ is the scattering momentum, $m_{\pi}$ the pion mass and
$\Lambda_{b}\sim 600$ MeV the breakdown scale. The choice of $k$ is not fixed
and we follow [41] and set $k=8$..

###  3.3 Density matrix formalism

In this section, we describe the theoretical framework for calculating matrix
elements in Eq. (8). We can separate this calculation into two parts: the
first part is the interaction kernel, which captures the one- and few-body
currents, and the second part is the structure, which contains the accurate
eigenstates of the Hamiltonian of the nucleus.

In the traditional method, we calculate the cross-section by integrating the
interaction kernel directly with the wave function of the nucleus.

$\displaystyle\braket{\Psi^{\prime}}{O}{\Psi}_{\rm traditional}$
$\displaystyle=\int\\{\mathrm{d}p_{i}\\}\Psi^{{}^{\prime}\dagger}(\\{p_{i}\\})\Psi(\\{p_{i}\\})O(\\{p_{i}\\})\,,$
(28)

where $\\{p_{i}\\}$ are the momentum of the internal nucleons $i$, the quantum
numbers of the nucleons are suppressed for simplicity.

Ref. [42] introduced an alternative method using transition density amplitudes
$\rho$ (referred to as densities)

$\displaystyle\braket{\Psi^{\prime}}{O}{\Psi}_{\rm density\ formalism}$
$\displaystyle=\sum_{\\{n_{i}\\}}\rho(\\{n_{i}\\})\otimes O(\\{n_{i}\\})\,,$
(29)

where $\\{n_{i}\\}$ are the quantum number of the internal nucleons $i$. This
approach significantly improves the computational time; for Compton scattering
the computation time gets reduced by a factor of ten [42]. It also separates
the calculations of structure part and operator part so that it is more easy
to apply the same operators to different nuclei. Finally, the densities and
operators can be independently tested mitigating the danger of computational
mistakes. The densities are freely available using a python package for
accessing the files [43]. We have used this method to calculate the matrix
elements and we found an order-of-magnitude speed-up. This works because the
densities do not depend on the interaction probes (in our case, scalar dark
matter interaction) and the nuclear-structure part of the calculation gets
factorized from the interaction kernel part.

For any external-probe matrix element, we need to read these nuclear densities
and do convolution with the appropriate interaction kernels that encode the
one- and two-body current operators on the momentum-spin basis. The
computational effort associated with this nuclear structure piece of the
calculations increases significantly with $A$, but highly parallelized and
optimized codes exist that solve the Schrödinger equation and obtain the wave
functions of light nuclei. Constructing densities from those wave functions is
straightforward. In short, once we generate the nuclear densities using a
supercomputer, all future matrix element calculations can be performed on
personal computers by reading these densities. A detailed calculation of one-
and two-body densities is presented in Ref. [42] for $A=3$. In Appendix A, we
briefly introduce the generalization to 4He.

##  4 Results and Discussion

We calculated the response functions
$\mathcal{F}_{\mathrm{is}}^{(0,1)}(\mbox{\boldmath$q$}^{2})$ for transfer
momenta in the range $q=(0-200)$ MeV using wave functions of different chiral
orders and a range of cut-offs
$\\{R_{2},R_{3},R_{4},R_{5}\\}=\\{0.9,1.0,1.1,1.2\\}$ fm (we have omitted
$R_{1}=0.8$ fm as it leads to a poorer description of the $N\\!N$ data than
$R_{2}$). Before we discuss the wave function dependence and theoretical
uncertainties of our computations, we first analyze the general features using
a fixed wave function obtained from the cut-off $R_{2}$.

In Fig. 2, we present the response functions for 2H, 3He, and 4He nuclei as
function of the transferred momentum. The maximum momentum transfer is related
to the DM and target mass and the DM velocity distribution. Using the standard
distribution and the DM escape velocity $v_{\chi}^{\mathrm{esc}}\simeq 550$
km/s [44], and assuming $m_{\chi}\sim m_{{}^{4}\mathrm{He}}$ the maximum
momentum transfer is roughly $|\mbox{\boldmath$q$}|\leq 3\times A$ MeV, so
about $12$ MeV for 4He [15]. We show results for larger momentum transfer as
well to analyze the accuracy and cut-off dependence of our results, and for
non-standard DM velocity distributions.

The LO response function
$\left|\mathcal{F}_{\mathrm{is}}^{(0)}(\mbox{\boldmath$q$}^{2})\right|^{2}$,
including just the one-nucleon DM interactions, is shown for the three light
nuclei in Fig. 2 in blue. The result use N4LO wave functions and a fixed cut-
off. At zero momentum transfer, the LO response functions equal unity because
of the normalization condition. The NLO response function includes the radius
and two-body currents and is defined as

$\displaystyle\left|\mathcal{F}_{\mathrm{is}}^{(0+1)}(\mbox{\boldmath$q$}^{2})\right|^{2}\equiv\left|\mathcal{F}_{\mathrm{is}}^{(0)}(\mbox{\boldmath$q$}^{2})+\mathcal{F}_{\mathrm{is},r}^{(1)}(\mbox{\boldmath$q$}^{2})+\mathcal{F}_{\mathrm{is},2b}^{(1)}(\mbox{\boldmath$q$}^{2})\right|^{2}\,.$
(30)

and is shown in orange. We conclude that
$\mathcal{F}_{\mathrm{is}}^{(0)}(\mbox{\boldmath$q$}^{2})$ is a good
approximation for the considered range of transfer momentum. The NLO
contributions only change the LO results by a few percent for 2H, 3He, and
4He. Despite the larger binding energies, the NLO corrections does not grow
for 4He and remain significantly smaller than power-counting estimates. It is
intriguing that much larger corrections, up to $20$-$30$%, are found for
heavier systems [16], and it would be interesting to find for which size of
nucleus this behavour kicks in. For 2H and 3He, the NLO corrections decrease
for larger transfer momentum due to a cancellation between radius and two-body
NLO corrections. This is not true for 4He where both terms have the same sign
although, as discussed below, this does depend on the applied regulator.

(a) 2H

(b) 3He

(c) 4He

Figure 2: Isoscalar structure functions of 2H, 3He, and 4He as a function of
transfer momentum $q$. We apply N4LO wave functions with a fixed cut-off
$R_{2}=0.9$ fm. The chiral order $\nu=0\ (0+1)$ is shown by the blue (orange)
line.  Figure 3: The percentage relative contribution of radius corrections of
4He as a function of transfer momentum $q$ for cut-off $R_{2}$ and the N4LO
wave function.  Figure 4: The percentage relative contribution of two-body
corrections of 2H, 3He, and 4He as a function of transfer momentum $q$ for
various cut-offs for N2LO (left) and N4LO (right) chiral wave functions.
Figure 5: Isoscalar one-body structure functions as a function of transfer
momentum $q$ for LO to N4LO chiral wave functions. The left (right) panel
matrix elements are calculated using cut-off $R_{5}=1.2\mbox{ fm}$
($R_{2}=0.9\mbox{ fm}$).

To capture the details of the NLO contribution, we introduce the relative
radius and two-body corrections

$\displaystyle\Delta^{(r)}$
$\displaystyle=\frac{\left|\mathcal{F}_{\mathrm{is}}^{(0+1)}(\mbox{\boldmath$q$}^{2})\right|^{2}-\left|\mathcal{F}_{\mathrm{is}}^{(0)}(\mbox{\boldmath$q$}^{2})+\mathcal{F}_{\mathrm{is},2b}^{(1)}(\mbox{\boldmath$q$}^{2})\right|^{2}}{\left|\mathcal{F}_{\mathrm{is}}^{(0+1)}(\mbox{\boldmath$q$}^{2})\right|^{2}}\,,$
$\displaystyle\Delta^{(2b)}$
$\displaystyle=\frac{\left|\mathcal{F}_{\mathrm{is}}^{(0+1)}(\mbox{\boldmath$q$}^{2})\right|^{2}-\left|\mathcal{F}_{\mathrm{is}}^{(0)}(\mbox{\boldmath$q$}^{2})+\mathcal{F}_{\mathrm{is},r}^{(1)}(\mbox{\boldmath$q$}^{2})\right|^{2}}{\left|\mathcal{F}_{\mathrm{is}}^{(0+1)}(\mbox{\boldmath$q$}^{2})\right|^{2}}\,.$
(31)

The radius contribution is given by

$\displaystyle\Delta^{(r)}\simeq\frac{2\mathcal{F}_{\mathrm{is},r}^{(1)}(\mbox{\boldmath$q$}^{2})}{\mathcal{F}_{\mathrm{is}}^{(0+1)}(\mbox{\boldmath$q$}^{2})}\simeq-\frac{2}{\sigma_{\pi
N}}\frac{9g_{A}^{2}\pi m_{\pi}^{3}}{4(4\pi
F_{\pi})^{2}}F\left(\frac{\left|\mbox{\boldmath$q$}\right|}{2m_{\pi}}\right)\xrightarrow{\left|q\right|\ll
2m_{\pi}}-\frac{5g_{A}^{2}\pi
m_{\pi}\left|\mbox{\boldmath$q$}\right|^{2}}{8(4\pi F_{\pi})^{2}\sigma_{\pi
N}}\,.$ (32)

The radius contribution vanishes at lower momentum but becomes dominant at
higher momentum and makes the net NLO correction negative. Furthermore, all
the nuclear effects drops out and $\Delta^{(r)}$ shows identical momentum
dependence for all nuclei, this has been verified by our results. In Fig. 3,
we show the momentum dependence of $\Delta^{(r)}$ for the 4He nucleus. We
omitted other nuclei and cutoffs because of the identical behavior.

In Fig. 4, we present our results on $\Delta^{(2b)}$ for various nuclei. In
the left panels, we apply N2LO chiral wave functions and show results for the
four applied regulators. The two-body corrections are small and lie well below
power-counting estimates as was found in Refs. [15, 18]. What is worrisome is
that the regulator dependence is large, ranging from
$\Delta^{(2b)}({}^{3}\mathrm{He})=(-0.01\pm 0.13)\%$ to
$\Delta^{(2b)}({}^{4}\mathrm{He})=(0.4\pm 1.9)\%$ at $q^{2}=0$. Even the signs
are uncertain. One might hope this would improve once higher-order wave
functions are applied and the corresponding results are shown in the right
panels. The 2H and 3He results are indeed essentially regulator independent
for N4LO wave functions, confirming results found in Ref. [15]. However, no
such convergence is shown for 4He as clearly seen in bottom-right panel of
Fig. 4.

This is puzzling. The power counting indicates that at this order the only NLO
currents are coming from one-body radius corrections and the pion-range two-
body corrections. There are no free counter terms so the calculations should
show regulator-independent matrix elements. We stress that similar cut-off
dependence is seen in quantum Monte Carlo calculations of DM response
functions for light nuclei using phenomenological $N\\!N$ potentials [18]. In
Fig. 5, we demonstrate that the power-counting issues only affect the two-body
currents. The one-body corrections show good convergence for all nuclei even
for N2LO wave functions.

Figure 6: Isoscalar structure functions of 4He up to chiral order $\nu=0+1$ as
a function of cut-offs for different transfer momentum $q$ . The first and
second panel indicate the leading order ($\nu=0$) one-nucleon isoscalar
structure functions and the full order ($\nu=0+1$) isoscalar structure
functions, respectively. The third (fourth) panel indicate the percent radius
(two-body) corrections. We presents the results for chiral wave functions
(N2LO to N4LO+). We omitted the results for LO and NLO wave functions due to
the large uncertainty bands.

In Fig. 6, we present our findings in a different way. We plot LO and NLO
isoscalar structure functions of 4He for three values of transfer momenta $q$
and four cut-offs. We have included different chiral order wave functions to
check the convergence of the matrix elements, LO and NLO wave functions are
not included. Apart from
$\mathcal{F}_{\mathrm{is},2b}^{(1)}(\mbox{\boldmath$q$}^{2})$, the isoscalar
structure functions display order-by-order convergence for wave functions and
exhibit regulator independence. The
$\mathcal{F}_{\mathrm{is}}^{(0)}(\mbox{\boldmath$q$}^{2})$ shows the expected
decrease in values for higher values of transfer momentum, which is in
agreement with Ref. [15]. The increase in uncertainty of
$\mathcal{F}_{\mathrm{is}}^{(0)}(\mbox{\boldmath$q$}^{2})$ for higher transfer
momentum is the reflection of the complicated nature of four-nucleon wave
functions.

NLO contributions occur only at few percent level. The relatively minute
uncertainty of $\Delta^{(r)}(\mbox{\boldmath$q$}^{2})$ can be understood as
the correction is proportional to the
$\mathcal{F}_{\mathrm{is}}^{(0)}(\mbox{\boldmath$q$}^{2})$, see Eq. (32). The
two-body corrections accounts for almost all the uncertainties associated with
the NLO corrections. The bottom panel of Fig. 6 shows that
$\mathcal{F}_{\mathrm{is},2b}^{(1)}(\mbox{\boldmath$q$}^{2})$ fails to
converge to a fixed value for different regulators and for highest-order wave
functions shows a linear dependence on the applied cut-off for all values of
transfer momenta.

The uncertainty bands in Fig. 6 have been obtained from the prescription in
Sect. 3.2. Each uncertainty band is centered at the average of the amplitude
values for different cut-offs. For properly renormalized results, we would
expect the amplitudes of different regulators to lie within the uncertainty
margins. As we can see in Fig. 6, this is observed for the LO and one-body NLO
corrections. For $\Delta^{(2b)}(\mbox{\boldmath$q$}^{2})$, on the other hand,
the amplitude values span beyond the uncertainty bands. In fact, for the N4LO+
wave function, none of the $\Delta^{(r)}(\mbox{\boldmath$q$}^{2})$ values fall
within the band.

What does this regulator dependence of two-nucleon matrix elements imply? It
could indicate a problem with the applied Weinberg power counting for scalar
currents. In principle, regulator dependence implies sensitivity to short-
distance physics which has to be absorbed by short-range operators, in this
case a DM-$N\\!N$ interaction. Such terms appear in Weinberg power counting
only at N3LO but might have to be promoted to lower order to ensure order-by-
order renormalization. There is by now a lot of evidence that this is the case
for neutrinoless double beta decay where a formally subleading contact term
must be included at leading order [19, 45, 46, 47]. If this is the case here,
it implies that including NLO scalar currents for DM-nucleus scattering
depends at least on one unknown LEC, not included in existing many-body
nuclear structure calculations [16, 48]. It would also impact our
understanding of the quark mass dependence on nuclear binding energies through
the nuclear sigma terms. Interestingly, given sufficiently accurate lattice
QCD data on nuclear binding energies as a function of the quark mass, it
should be possible to fit the scalar contact counter term to lattice data and
determining whether indeed these terms must be promoted to leading order [17,
49]. While it is tempting to argue that the missing counter term is not too
relevant considering that the two-nucleon matrix elements are found to be
small, this is dangerous as there is no guarantee that the finite part of the
counter term is of the same size and it could instead be a genuine
$\mathcal{O}(30\%)$ correction.

On the other hand, it might be that due to some accidental or so far not-
understood cancellation, that the NLO two-nucleon matrix elements are small
and essentially demoted to N2LO contributions. The observed regulator
dependence is indeed of that size (however, we are only varying the regulators
in a small range and the uncertainty might grow with a larger range). If true,
it would be interesting to investigate the N2LO response functions by
including higher-order scalar currents derived in Ref. [30]. Unfortunately,
these currents come with unknown low-energy constants.

Figure 7: The plots shows the correlation between the isoscalar structure
functions and the D-wave probability for 2H, 3He, and 4He represented by
circle, square, and diamond data points. The data points are assigned
different colors corresponding to the chiral order of the wave function. The
cut-off $R_{i}$ used for the calculation is indicated by label $i$ on the data
point. For 2H and 3He, we have included the cut-off $R_{1}=0.8\mbox{ fm}$ and
results from phenomenological wave functions from Ref. [15].

Finally, we note that the observed cut-off dependence of our results can be
cast in a different light. In Fig. 7, we present an observed linear
correlation between the two-body matrix element and the D-wave probability for
2H, 3He, and 4He. The D-wave probability describes how much of the nuclear
wave function involves $l=2$ $N\\!N$ states and is a consequence of $S-D$
mixing from the nuclear tensor force. It was pointed out by Friar [50] already
that the D-wave probability is not an observable: it can be changed by unitary
transformations of the $N\\!N$ potential. Therefore, the correlation with the
unobservable quantity also points to a potential flaw in the applied power
counting.

##  5 Conclusions

We have investigated DM scattering of light nuclei using the framework of
chiral effective field theory. Our main focus has been on scattering off the
4He nucleus, a potential target for direct-detection experiments, but we have
also investigated 2H and 3He nuclei. The nuclear wave functions were
calculated from chiral nuclear forces up to fourth order in the chiral
expansion. The light nuclei allowed us to calculate bound-state and scattering
equations directly with a high degree of accuracy. This allowed us to
investigate the uncertainties and convergence of observables for different
chiral orders of the wave functions with a high degree of precision.

We focused on the scalar interaction of DM and quarks. Using $\chi$PT, we
calculated the ensuing scalar hadronic current up to NLO (actually N2LO since
no new corrections appear at this order). The NLO current consists of radius
and two-body corrections. For 2H, 3He, and 4He, the NLO currents modify LO
results by a few percent only despite larger power-counting expectations.
Radius and two-nucleon corrections are found to be of similar size in all
three nuclei under investigation. The one-body and radius corrections are in
good agreement with earlier work [15, 18].

Taken at face value, our results imply that NLO corrections are small and the
LO response functions can be directly used in the analysis of experiments
using light nuclei as targets. Unfortunately, our results for NLO two-nucleon
corrections exhibit a clear regulator dependence. While, for 2H and 3He we
somewhat recover regulator independence when using high-order chiral wave
functions, this is not the case for 4He. Furthermore, in all cases we find a
linear correlation between the two-body matrix element and the D-wave
probability of the nucleus. Because there are no contact terms at NLO to
absorb the regulator dependence, the regulator dependence could indicate that
the scalar currents are not properly renormalized. The LO current and NLO
radius contributions do not show any such regulator dependence, indicating
that all regulator dependence arises from the two-body effects. A potential
remedy could be to promote formally subleading DM-nucleon-nucleon interactions
to lower order, as is done, for instance, in the description of neutrinoless
double beta decay. This would also impact existing predictions for dark matter
scattering of heavier nuclei. It could also be that two-nucleon effects should
be demoted to N2LO. More work is needed to clarify this issue for instance by
investigating the impact of higher-order scalar currents derived in Ref. [30]
in the light nuclei under consideration and beyond.

## Acknowledgments

We thank Martin Hoferichter and Achim Schwenk for useful discussions and
encouragement. JdV acknowledges support from the Dutch Research Council (NWO)
in the form of a VIDI grant. This work is further supported in part by the DFG
and the NSFC through funds provided to the Sino-German CRC 110 “Symmetries and
the Emergence of Structure in QCD” (NSFC Grant No. 12070131001, Project-ID
196253076 - TRR 110), and by the MKW NRW under the funding code NW21-024-A.
The numerical calculations were performed on JURECA of the Jülich
Supercomputing Centre, Jülich, Germany.

##  Appendix A Density matrix formalism

###  A.1 Kinematics and Partial-Wave Decomposition

We assume that the nucleus has $A$ nucleons, which in the initial state have
total angular momentum $J$, spin-projection $M$ onto the $z$-axis, and
isospin-projection $M_{T}$, several isospins $T$ may contribute in general.
Here, we will be mostly concerned with 4He for which $A=4,J=0,M_{T}=0$, the
dominant contribution comes from $T=0$, but isospin breaking could in
principle introduce small $T=1$ and $2$ contributions to the 4He wave
function. It turns out that these contributions are not important for the dark
matter interactions studied here.

Additionally, we also assume that the scattering happens in the center of mass
frame of the nucleus-DM system. We denote the momentum of the incoming
(outgoing) nucleon as $-\mbox{\boldmath$k$}$
($-\mbox{\boldmath$k$}^{\prime}$), and the incoming (outgoing) DM will have
the opposite momentum. The momentum of incoming DM is chosen to lie along the
$z$-axis $\mbox{\boldmath$k$}=k\boldsymbol{e}_{z}$, which also is the
quantization axis for the spin and angular momentum projections. We
furthermore assume that the scattering is happening in $x$-$z$ plane so that
$\mbox{\boldmath$k$}^{\prime}=k^{\prime}\left(\sin\vartheta\boldsymbol{e}_{x}+\cos\vartheta\boldsymbol{e}_{z}\right)$.
The elastic scattering ensures
$\left|\mbox{\boldmath$k$}\right|=\left|\mbox{\boldmath$k$}^{\prime}\right|$,
and the transfer momentum into nucleus is
$\mbox{\boldmath$Q$}=\mbox{\boldmath$k$}-\mbox{\boldmath$k$}^{\prime}$.

$1$$2$$3$$4$$\boldsymbol{p}_{12}$$\boldsymbol{q}_{4}$$\boldsymbol{p}_{3}$

$\displaystyle
1$$2$$3$$4$$\boldsymbol{p}_{12}$$\boldsymbol{q}$$\boldsymbol{p}_{34}$

Figure 8:  3+1 and 2+2 Jacobi coordinates of the four-nucleon system.

.

Relative Jacobi coordinates can be defined in various ways be defining
different subsystems. For the 4N system, except for permutations of particles,
one finds only two distinct ways for arranging these subclusters, the so-
called ‘3+1’ and ‘2+2’ coordinates, as depicted in Fig. 8. The ‘3+1’
coordinates singles out the last nucleon and provides a relative momentum of
this last nucleon with respect to the remaining $A-1$ nucleons. This
coordinate will be most useful for the definition of one-nucleon densities
below. The ‘2+2’ coordinates single out the pair of nucleons $1$ and $2$
moving relative to the remaining $A-2$ nucleons. We will later use this
coordinate for defining the two-nucleon densities.

Denoting the nucleon momenta by $\mbox{\boldmath$k$}_{i}$, we can define the
Jacobi momentum of the spectator particle in the $A-1$+$1$ coordinates as

$\mbox{\boldmath$q$}_{A}=\frac{A-1}{A}\mbox{\boldmath$k$}_{A}-\frac{1}{A}\sum_{i=1}^{A-1}\mbox{\boldmath$k$}_{i}=\mbox{\boldmath$k$}_{A}+\frac{1}{A}\mbox{\boldmath$k$}$
(33)

which, for 4He ($A$=4) nucleus, is

$\mbox{\boldmath$q$}_{4}=\frac{3}{4}\mbox{\boldmath$k$}_{4}-\frac{1}{4}(\mbox{\boldmath$k$}_{1}+\mbox{\boldmath$k$}_{2}+\mbox{\boldmath$k$}_{3})=\mbox{\boldmath$k$}_{4}+\frac{1}{4}\mbox{\boldmath$k$}.$
(34)

This relation will be important for the implementation of one-body densities.

For the two-body densities, the relative and total pair momenta of the (12)
nucleon system

$\displaystyle\mbox{\boldmath$p$}_{12}$
$\displaystyle=\frac{1}{2}(\mbox{\boldmath$k$}_{1}-\mbox{\boldmath$k$}_{2}),$
$\displaystyle\mbox{\boldmath$k$}_{12}$
$\displaystyle=\mbox{\boldmath$k$}_{1}+\mbox{\boldmath$k$}_{2}$ (35)

are the momenta of interest. These can be easiest defined for $A-2$+$2$
coordinates. Conveniently, $\mbox{\boldmath$p$}_{12}$ is already part of this
set of Jacobi coordinates. The total pair momentum is related to the relative
momentum of the (12) subsystem and the remaining nucleons

$q$
$\displaystyle=\frac{A-2}{A}\mbox{\boldmath$k$}_{12}-\frac{2}{A}\sum_{i=3}^{A}\mbox{\boldmath$k$}_{i}=\mbox{\boldmath$k$}_{12}+\frac{2}{A}\mbox{\boldmath$k$}=\mbox{\boldmath$k$}_{12}+\frac{1}{2}\mbox{\boldmath$k$}$
(36)

where the last relation holds for $A=4$.

We will be focusing on 4He for the remaining discussions. We denote the spin-
projection state of 4He as $\ket{M}$, where we suppress the labels of $JM_{T}$
and the bound-state energy. This state is the eigenstate for the nucleus at
rest, total angular momentum, and its z-component. For 4He, the total angular
momentum of $J=0$ implies that $M=0$.

The wave function can be expressed in both kind of coordinates, ‘3+1’ and
‘2+2’. In this work, we use Yakubovsky equations to get the wave functions. In
this case, both representations of the wave functions are naturally obtained.
For calculations that are based on, e.g., the Jacobi no-core shell model [51,
52], wave functions are usually only obtained for ‘$A$-$1$–$1$’ coordinates.
However, also in this case, transformation to ‘2+2’ coordinates can be easily
performed. In following, we will assume that the wave functions are given in
‘3+1’ coordinates by

$\psi_{\alpha}(p_{12}p_{3}q_{4})=\braket{p_{12}p_{3}q_{4}\alpha}{M}$ (37)

where $\alpha$ represents the set of angular momentum and isospin quantum
numbers defined in Eq. (17). Similarly, the wave function in ‘2+2’ coordinates
reads

$\psi_{\beta}(p_{12}p_{34}q)=\braket{p_{12}p_{34}q\beta}{M}.$ (38)

The corresponding set of quantum numbers is defined in Eq. (18).

The labels of outgoing particles are primed. The states $\ket{M}$ must be
multiplied by an eigenstate of the nuclear cm momentum operator to conserve
the momentum. This will result in the
$\delta^{(3)}(\mbox{\boldmath$k$}-\mbox{\boldmath$k$}^{\prime}-\mbox{\boldmath$Q$})$.
These momentum delta functions are not included in the following calculations
to shorten the notation.

###  A.2 One-Body Density

Assuming that the active nucleon is the fourth one, the relevant matrix
element for one-body operator reads

$\displaystyle\bra{\mbox{\boldmath$k$}_{4}^{\prime}}{\bra{l_{4}^{\prime}m_{4}^{l\prime}}}\bra{s_{4}^{\prime}m_{4}^{s\prime}}\bra{t_{4}^{\prime}m_{4}^{t\prime}}\hat{O}_{4}(\mbox{\boldmath$k$},\mbox{\boldmath$Q$})\ket{t_{4}m^{t}_{4}}\ket{s_{4}m^{s}_{4}}{\ket{l_{4}m_{4}^{l}}}\ket{\mbox{\boldmath$k$}_{4}}$
$\displaystyle=\delta_{m^{t^{\prime}}_{4},m^{t}_{4}}\delta^{(3)}(\mbox{\boldmath$k$}_{4}^{\prime}-\mbox{\boldmath$k$}_{4}-\mbox{\boldmath$Q$})O_{4}(m^{s^{\prime}}_{4}m^{s}_{4}m^{t}_{4};\mbox{\boldmath$k$}_{4};\mbox{\boldmath$k$},\mbox{\boldmath$Q$}).$
(39)

Here $m^{t}_{4}$ is conserved because the DM interaction does not change the
charge of the struck nucleon. For the one-body density, we allow for a
dependence of the operator on the momentum $\mbox{\boldmath$k$}_{4}$ of the
active nucleon which is necessary to take boost corrections into account as
they appear in subleading operators. This explicit dependence on the momentum
of the active nucleon is included in terms of a multipole expansion which we
truncate at order $K_{\rm max}$.

$O_{4}(m_{4}^{s\prime}m^{s}_{4}m^{t}_{4};\mbox{\boldmath$k$}_{4};\mbox{\boldmath$k$},\mbox{\boldmath$Q$})\equiv\sum_{K=0}^{K_{\rm
max}}\sum_{\kappa=-K}^{K}\sqrt{\frac{4\pi}{2K+1}}(k_{4})^{K}Y_{K\kappa}(\hat{k}_{4})\tilde{O}_{4}(m_{4}^{s\prime}m^{s}_{4}m^{t}_{4};K\kappa;\mbox{\boldmath$k$},\mbox{\boldmath$Q$}),$
(40)

where $\hat{k}_{4}$ is the unit vector of $\mbox{\boldmath$k$}_{4}$. In
practice, we only use $K_{\rm max}=1$.

Now let us look at the matrix element of $\hat{O}_{4}$ that we want to
calculate

$\displaystyle\bra{M^{\prime}}$
$\displaystyle\hat{O}_{4}(\mbox{\boldmath$k$},\mbox{\boldmath$Q$})\ket{M}=$
$\displaystyle\sum_{\alpha\alpha^{\prime}}\int\mathrm{d}p_{12}\ p_{12}^{2}\
\mathrm{d}p_{3}\ p_{3}^{2}\ \mathrm{d}q_{4}\ q_{4}^{2}\
\mathrm{d}p_{12}^{\prime}\ p_{12}^{{}^{\prime}2}\ \mathrm{d}p_{3}^{\prime}\
p_{3}^{{}^{\prime}2}\ \mathrm{d}q_{4}^{\prime}\ q_{4}^{{}^{\prime}2}\
\psi^{\dagger}_{\alpha^{\prime}}(p_{12}^{\prime}p_{3}^{\prime}q_{4}^{\prime})\psi_{\alpha}(p_{12}p_{3}q_{4})$
$\displaystyle\times\bra{p_{12}^{\prime}p_{3}^{\prime}q_{4}^{\prime},\left[[l_{12}^{\prime}s_{12}^{\prime}]j_{12}^{\prime}(l_{3}^{\prime}s_{3}^{\prime})I_{3}^{\prime}]j_{3}^{\prime}(l_{4}^{\prime}s_{4}^{\prime})I_{4}^{\prime}\right]J^{\prime}M^{\prime},\left[(t_{12}^{\prime}t_{3}^{\prime})\tau_{3}^{\prime}t_{4}^{\prime}\right]T^{\prime}M_{T}^{\prime}}\hat{O}_{4}(\mbox{\boldmath$k$},\mbox{\boldmath$Q$})$
$\displaystyle\hskip
28.45274pt\ket{p_{12}p_{3},\left[[l_{12}s_{12}]j_{12}(l_{3}s_{3})I_{3}]j_{3}(l_{4}s_{4})I_{4}\right]JM,\left[(t_{12}t_{3})\tau_{3}t_{4}\right]TM_{T}}.$
(41)

Using Clebsch-Gordan coefficients $\left\langle
j_{1}m_{1},j_{2}m_{2}\middle|jm\right\rangle$ we can explicitly decompose
$\alpha$ so as to separate the spin-isospin quantum numbers of the active
nucleon 4 from those of the other three nucleons:

$\displaystyle\bra{M^{\prime}}\hat{O}_{4}(\mbox{\boldmath$k$},\mbox{\boldmath$Q$})\ket{M}=$
$\displaystyle\sum_{\alpha\alpha^{\prime}}\int\mathrm{d}p_{12}\ p_{12}^{2}\
\mathrm{d}p_{3}\ p_{3}^{2}\ \mathrm{d}q_{4}\ q_{4}^{2}\
\mathrm{d}p_{12}^{\prime}\ p_{12}^{{}^{\prime}2}\ \mathrm{d}p_{3}^{\prime}\
p_{3}^{{}^{\prime}2}\ \mathrm{d}q_{4}^{\prime}\ q_{4}^{{}^{\prime}2}\
\psi^{\dagger}_{\alpha^{\prime}}(p_{12}^{\prime}p_{3}^{\prime}q_{4}^{\prime})\psi_{\alpha}(p_{12}p_{3}q_{4})$
$\displaystyle\times\sum_{m_{4}^{\prime}}\left\langle
j_{3}^{\prime}(M^{\prime}-m_{4}^{\prime}),I_{4}^{\prime}m_{4}^{\prime}\middle|J^{\prime}M^{\prime}\right\rangle\sum_{m_{4}^{t\prime}}\left\langle\tau_{3}^{\prime}(M_{T}^{\prime}-m_{4}^{t\prime}),t_{4}^{\prime}m_{4}^{t\prime}\middle|T^{\prime}M_{T}^{\prime}\right\rangle$
$\displaystyle\times\sum_{m_{4}}\left\langle
j_{3}(M-m_{4}),I_{4}m_{4}\middle|JM\right\rangle\sum_{m^{t}_{4}}\left\langle\tau_{3}(M_{T}-m^{t}_{4}),t_{4}m^{t}_{4}\middle|TM_{T}\right\rangle$
$\displaystyle\times\bra{p_{12}^{\prime}p_{3}^{\prime},[l_{12}^{\prime}s_{12}^{\prime}]j_{12}^{\prime}(l_{3}^{\prime}s_{3}^{\prime})I_{3}^{\prime}]j_{3}^{\prime}(M^{\prime}-m_{4}^{\prime}),(t_{12}^{\prime}t_{3}^{\prime})\tau_{3}^{\prime}(M_{T}^{\prime}-m_{4}^{t\prime})}$
$\displaystyle\qquad\ket{p_{12}p_{3},[l_{12}s_{12}]j_{12}(l_{3}s_{3})I_{3}]j_{3}(M-m_{4}),(t_{12}t_{3})\tau_{3}(M_{T}-m^{t}_{4})}$
$\displaystyle\times\braket{q_{4}^{\prime},(l_{4}^{\prime}s_{4}^{\prime})I_{4}^{\prime}m_{4}^{\prime},t_{4}^{\prime}m_{4}^{t\prime}}{\hat{O}_{4}(\mbox{\boldmath$k$},\mbox{\boldmath$Q$})}{q_{4},(l_{4}s_{4})I_{4}m_{4},t_{4}m^{t}_{4}}$
(42)

Here we used identities from the spin projections:
$M=m_{4}+m_{123},M_{T}=m^{t}_{4}+m^{t}_{123},\ldots$. Now we expand the last
term in the above equation using Eq. (39), the momentum conserving delta
function can be rewritten using Eq. (34) as

$\delta^{(3)}(\mbox{\boldmath$k$}_{4}^{\prime}-\mbox{\boldmath$k$}_{4}-\mbox{\boldmath$Q$})=\delta^{(3)}(\mbox{\boldmath$q$}_{4}^{\prime}-\mbox{\boldmath$q$}_{4}-\frac{3}{4}\mbox{\boldmath$Q$})$
(43)

We further expand the last line of Eq. (42) by inserting the solid angles
$\hat{q}_{4},\hat{q}_{4}^{\prime}$ of the third particles’ momentum. The state
of the active nucleon is further decoupled into an orbital and spin part which
introduces summations of spin projections $m_{4}^{s\prime}$ and $m^{s}_{4}$ .
Then we use the result in Eq. (40), and perform the integrations and
summations corresponding to the appropriate delta functions to arrive at the
final form

$\displaystyle\bra{M^{\prime}}\hat{O}_{4}(\mbox{\boldmath$k$},\mbox{\boldmath$Q$})\ket{M}=\sum_{K=0}^{K_{\rm
max}}\sum_{\kappa=-K}^{K}\sum_{\begin{subarray}{c}m_{4}^{s^{\prime}}m_{4}^{s}\\\
m^{t}_{4}\end{subarray}}\tilde{O}_{4}(m_{4}^{s\prime}m^{s}_{4}m^{t}_{4};K\kappa;\mbox{\boldmath$k$},\mbox{\boldmath$Q$})\rho^{K\kappa;m^{t}_{4}M_{T},M^{\prime}M}_{m_{4}^{s\prime}m^{s}_{4}}(\mbox{\boldmath$k$},\mbox{\boldmath$Q$})$

where we define the one-body (transition) density by summing over those
quantum numbers that are not involved in the interaction:

$\displaystyle\rho^{K\kappa;m^{t}_{4}M_{T},M^{\prime}M}_{m_{4}^{s\prime}m^{s}_{4}}(\mbox{\boldmath$k$},\mbox{\boldmath$Q$})$
$\displaystyle=\sum_{\alpha\alpha^{\prime}}\int\mathrm{d}p_{12}\ p_{12}^{2}\
\mathrm{d}p_{3}\ p_{3}^{2}\ \mathrm{d}q_{4}\ q_{4}^{2}\ $
$\displaystyle\times\delta_{l_{12}^{\prime}l_{12}}\,\delta_{s_{12}^{\prime}s_{12}}\,\delta_{j_{12}^{\prime}j_{12}}\,\delta_{t_{12}^{\prime}t_{12}}\,\delta_{l_{3}^{\prime}l_{3}}\,\delta_{I_{3}^{\prime}I_{3}}\,\delta_{j_{3}^{\prime}j_{3}}\delta_{\tau_{3}^{\prime}\tau_{3}}\,\delta_{M_{T}^{\prime}M_{T}}\,$
$\displaystyle\times\sum_{m_{4}}\left\langle
j_{3}^{\prime}(M-m_{4}),I_{4}^{\prime}M^{\prime}-M+m_{4}\middle|J^{\prime}M^{\prime}\right\rangle\left\langle\tau_{3}^{\prime}(M_{T}-m^{t}_{4}),t_{4}^{\prime}m^{t}_{4}\middle|T^{\prime}M_{T}\right\rangle$
$\displaystyle\times\left\langle
j_{3}(M-m_{4}),I_{4}m_{4}\middle|JM\right\rangle\left\langle\tau_{3}(M_{T}-m^{t}_{4}),t_{4}m^{t}_{4}\middle|TM_{T}\right\rangle$
$\displaystyle\times\left\langle
l_{4}^{\prime}(M^{\prime}-M+m_{4}-m_{4}^{s\prime}),s_{4}^{\prime}m_{4}^{s\prime}\middle|I_{4}^{\prime}(M^{\prime}-M+m_{4})\right\rangle\left\langle
l_{4}(m_{4}-m^{s}_{4}),s_{4}m^{s}_{4}\middle|I_{4}m_{4}\right\rangle$
$\displaystyle\times\int\mathrm{d}\hat{q}_{4}Y^{\dagger}_{l_{4}^{\prime}(M^{\prime}-M+m_{4}-m_{4}^{s\prime})}(\widehat{\mbox{\boldmath$q$}_{4}+\tfrac{3}{4}\mbox{\boldmath$Q$}})Y_{l_{4}(m_{4}-m^{s}_{4})}(\hat{}\mbox{\boldmath$q$}_{4})$
$\displaystyle\times\sum_{K=0}^{K_{\rm
max}}\sum_{\kappa=-K}^{K}\sqrt{\frac{4\pi}{2K+1}}|\mbox{\boldmath$q$}_{4}-\tfrac{\mbox{\boldmath$k$}}{4}|^{K}Y_{K\kappa}(\widehat{\mbox{\boldmath$q$}_{4}-\tfrac{\mbox{\boldmath$k$}}{4}})\psi^{\dagger}_{\alpha^{\prime}}(p_{12}\left|\mbox{\boldmath$q$}_{4}+\tfrac{3}{4}\mbox{\boldmath$Q$}\right|)\psi_{\alpha}(p_{12}p_{3}q_{4})\
.$ (44)

This definition separates the calculation of the operator from the calculation
of the wave functions. Based on this definition, we can define transition
densities for arbitrary nuclear wave functions that can then be combined with
separately calculated matrix elements of operators.

###  A.3 Two-Body Density

In this work, we take into account one- and two-nucleon operators contributing
to elastic DM scattering. For the application of two-nucleon operators, the
relevant matrix element is given by

$\displaystyle\bra{\mbox{\boldmath$p$}_{12}^{\prime}\mbox{\boldmath$k$}_{12}^{\prime}}\bra{s_{12}^{\prime}m_{12}^{s\prime}}\bra{t_{12}^{\prime}m_{12}^{t\prime}}\hat{O}_{12}(\mbox{\boldmath$k$},\mbox{\boldmath$Q$})\ket{t_{12}m^{t}_{12}}\ket{s_{12}m^{s}_{12}}\ket{{\mbox{\boldmath$p$}_{12}\mbox{\boldmath$k$}_{12}}}$
$\displaystyle=\delta^{(3)}(\mbox{\boldmath$k$}_{12}^{\prime}-\mbox{\boldmath$k$}_{12}-\mbox{\boldmath$Q$})\delta_{m_{12}^{t\prime}m^{t}_{12}}O_{12}(s_{12}^{\prime}t_{12}^{\prime}m_{12}^{s\prime}s_{12}t_{12}m^{t}_{12};\mbox{\boldmath$p$}_{12}^{\prime},\mbox{\boldmath$p$}_{12};\mbox{\boldmath$Q$}).$
(45)

For these operators, that generally contribute only in higher orders, it will
not be necessary to take boost corrections into account. Therefore, we assume
that the operator does not have an explicit dependence on
$\boldsymbol{k}_{12}$. The two-nucleon operator is usually represented in
terms of partial-wave states, which can be written as,

$\displaystyle\braket{\alpha_{12}^{\prime}p_{12}^{\prime}\mbox{\boldmath$k$}_{12}^{\prime}}{\hat{O}_{12}}{\alpha_{12}p_{12}\mbox{\boldmath$k$}_{12}}=\delta_{m_{12}^{t\prime}m^{t}_{12}}\delta^{(3)}(\mbox{\boldmath$k$}_{12}^{\prime}-\mbox{\boldmath$k$}_{12}-\mbox{\boldmath$Q$})$
$\displaystyle\times\sum_{m_{12}^{s\prime}m^{s}_{12}}\left\langle
l_{12}(m_{12}-m^{s}_{12}),s_{12}m^{s}_{12}\middle|j_{12}m_{12}\right\rangle\left\langle
l_{12}^{\prime}(m_{12}^{\prime}-m_{12}^{s\prime}),s_{12}^{\prime}m_{12}^{s\prime}\middle|j_{12}^{\prime}m_{12}^{\prime}\right\rangle$
$\displaystyle\times\int\mathrm{d}\hat{p}_{12}^{\prime}\mathrm{d}\hat{p}_{12}Y^{\dagger}_{l_{12}^{\prime}(m_{12}^{\prime}-m_{12}^{s\prime})}(\hat{p}_{12}^{\prime})Y_{l_{12}(m_{12}-m^{s}_{12})}(\hat{p}_{12})$
$\displaystyle\times
O_{12}(s_{12}^{\prime}t_{12}^{\prime}m_{12}^{s\prime}s_{12}t_{12}m^{s}_{12}m^{t}_{12};\mbox{\boldmath$p$}_{12}^{\prime},\mbox{\boldmath$p$}_{12};\mbox{\boldmath$Q$})$
$\displaystyle\equiv\delta_{m_{12}^{t\prime}m^{t}_{12}}\delta^{(3)}(\mbox{\boldmath$k$}_{12}^{\prime}-\mbox{\boldmath$k$}_{12}-\mbox{\boldmath$Q$})O_{12}^{\alpha_{12}^{\prime}\alpha_{12}}(p_{12}^{\prime},p_{12};\mbox{\boldmath$Q$}).$
(46)

For the calculation of matrix elements of such an operator, it is advisable to
use coordinates that separate of a two-nucleon cluster from a cluster of the
remaining nucleons. For an $A=4$ system, such coordinates are usually labeled
as ‘2+2’ coordinates and are defined in Eq. (18). The scattering process
changes the momenta within the subcluster (12) and the relative momentum $q$
of this subcluster and the remaining nucleons. Using Eq. (36) and Eq. (35)
along with $\sum_{i}\mbox{\boldmath$k$}_{i}=-\mbox{\boldmath$k$}$, we can
rewrite the delta function as,

$\delta^{(3)}(\mbox{\boldmath$k$}_{12}^{\prime}-\mbox{\boldmath$k$}_{12}-\mbox{\boldmath$Q$})=\delta^{(3)}(\mbox{\boldmath$q$}-\mbox{\boldmath$q$}^{\prime}+\frac{1}{2}\mbox{\boldmath$Q$}),$
(47)

we can use this delta function to integrate $q$ while calculating the matrix
element. Now we can write the two-body matrix element following the same
procedure as one-body matrix element and we get the factorization,

$\bra{M^{\prime}}\hat{O}_{12}(\mbox{\boldmath$q$})\ket{M}\equiv\sum_{\alpha_{12}^{\prime}\alpha_{12}}\int\mathrm{d}p_{12}p_{12}^{2}\mathrm{d}p_{12}^{\prime}\mathrm{d}p_{12}^{{}^{\prime}2}O_{12}^{\alpha_{12}^{\prime}\alpha_{12}}(p_{12}^{\prime},p_{12};\mbox{\boldmath$Q$})\rho_{\alpha_{12}^{\prime}\alpha_{12}}^{M_{T},M^{\prime},M}(p_{12}^{\prime},p_{12};\mbox{\boldmath$Q$})$
(48)

where we define the two-body (transition) density as:

$\displaystyle\rho_{\alpha_{12}^{\prime}\alpha_{12}}^{M_{T},M^{\prime},M}(p_{12}^{\prime}p_{12};\mbox{\boldmath$Q$})$
$\displaystyle=\sum_{\alpha^{\prime}(\alpha_{12}^{\prime})\alpha(\alpha_{12})}\sum_{m_{I}}\left\langle
Im_{I},j_{34}M-m_{I}\middle|JM\right\rangle\left\langle
I^{\prime}m_{I}+M^{\prime}-M,j_{34}M-m_{I}\middle|J^{\prime}M^{\prime}\right\rangle\delta_{j_{34}j_{34}^{{}^{\prime}}}$
$\displaystyle\times\left\langle
j_{12}m_{12},\lambda(m_{I}-m_{12})\middle|Im_{I}\right\rangle\left\langle
j_{12}^{\prime}m_{12}^{\prime},\lambda^{\prime}(M^{\prime}-M+m_{I}-m_{12}^{\prime})\middle|I^{\prime}(m_{I}+M^{\prime}-M)\right\rangle$
$\displaystyle\times\left\langle
t_{12}m^{t}_{12},t_{34}(M_{T}-m^{t}_{12})\middle|TM_{T}\right\rangle\left\langle
t_{12}^{\prime}m^{t}_{12},t_{34}^{\prime}(M_{T}-m^{t}_{12})\middle|T^{\prime}M_{T}\right\rangle$
$\displaystyle\times\int\mathrm{d}q^{\prime}q^{{}^{\prime}2}\int\mathrm{d}\hat{q}^{\prime}\int\mathrm{d}p_{34}^{\prime}p_{34}^{{}^{\prime}2}\
\psi^{\dagger}_{\alpha^{\prime}}(p_{12}^{\prime}p_{34}^{\prime}q^{\prime})\psi_{\alpha}(p_{12}p_{34}^{\prime}\left|\mbox{\boldmath$q$}^{\prime}+\tfrac{1}{2}\mbox{\boldmath$Q$}\right|)$
$\displaystyle\times
Y^{\dagger}_{\lambda^{\prime}(M^{\prime}-M+m_{I}^{\prime}-m_{12}^{\prime})}(\hat{q}^{\prime})Y_{\lambda(m_{I}-m_{12})}(\widehat{\mbox{\boldmath$q$}^{\prime}+\tfrac{1}{2}\mbox{\boldmath$Q$}}).$
(49)

##  Appendix B Partial wave decomposition

In this section, we give the partial wave decomposition of two-body operator
$O_{12}^{\alpha_{12}^{\prime}\alpha_{12}}(p_{12}^{\prime},p_{12};\mbox{\boldmath$Q$})$
in (48). For notational convenience, we will drop the subscript (12) from now
on. We can decompose the operator into total-spin part and isospin part

$\displaystyle
O^{\alpha^{\prime}\alpha}_{m_{\chi}^{\prime}m_{\chi}}(p^{\prime},p;\mbox{\boldmath$Q$})$
$\displaystyle=\left(\sum_{\xi=0}^{\infty}\sum_{|m_{\xi}|\leq\xi}\left\langle
jm_{j},\xi
m_{\xi}\middle|j^{\prime}m_{j}^{\prime}\right\rangle\mathcal{O}_{(\alpha^{\prime}\alpha)\xi
m_{\xi},m_{\chi}^{\prime}m_{\chi}}^{(\text{spin})}\right)$
$\displaystyle\hskip
28.45274pt\times\left(\sum_{\tau=0}^{2}\sum_{|m_{\tau}|\leq\tau}\left\langle
tm_{t},\tau
m_{\tau}\middle|t^{\prime}m_{t}\right\rangle\mathcal{O}_{(t^{\prime}t)\tau
m_{\tau}}^{(\text{iso})}\right)\,,$ (50)

with

$\displaystyle\mathcal{O}_{(\alpha^{\prime}\alpha)\xi
m_{\xi},m_{\chi}^{\prime}m_{\chi}}^{(\text{spin})}=\sqrt{\hat{j}\hat{l}^{\prime}\hat{s}^{\prime}\hat{\xi}}\sum_{\lambda\sigma}\left(\begin{array}[]{ccc}j&\xi&j^{\prime}\\\
l&\lambda&l^{\prime}\\\
s&\sigma&s^{\prime}\end{array}\right)\mathcal{O}_{((l^{\prime}l)\lambda,(s^{\prime}s)\sigma)\xi
m_{\xi},m_{\chi}^{\prime}m_{\chi}}(p^{\prime},p,\mbox{\boldmath$Q$})\,,$ (54)

and

$\displaystyle\mathcal{O}_{(t^{\prime}t)\tau
m_{\tau}}^{(\text{iso})}=\frac{\hat{\tau}}{\hat{t}^{\prime}}\sum_{m_{t}}\left\langle
tm_{t},\tau m_{\tau}\middle|t^{\prime}m_{t}\right\rangle\langle
t^{\prime}m_{t}|O^{(\text{iso})}|tm_{t}\rangle\,,$ (55)

where $m_{\chi}(m_{\chi}^{\prime})$ is the spin quantum number of the incoming
(outgoing) dark matter, $\hat{x}\equiv 2x+1$, and $\left\langle
j_{1}m_{1},j_{2}m_{2}\middle|jm\right\rangle$ are the Clebsch-Gordan
coefficients. The partial wave decomposition of the isospin part is given in
Table 1.

$t^{\prime}$ | $t$ | $\tau$ | $m_{\tau}$ | $\tilde{\mathcal{O}}_{(t^{\prime}t)\tau m_{\tau}}^{(\text{iso})}$
---|---|---|---|---
0 | 0 | 0 | 0 | -3
1 | 1 | 0 | 0 | 1
Table 1: Partial wave decomposition of isospin operator $O^{(\text{iso})}$ as
given in Eq. (55) for non-zero channels.

The decomposed operator can be factorized in an angular momentum- and spin-
dependent parts

$\displaystyle\mathcal{O}_{((l^{\prime}l)\lambda,(s^{\prime}s)\sigma)\xi
m_{\xi},m_{\chi}^{\prime}m_{\chi}}$
$\displaystyle(p^{\prime},p;\mbox{\boldmath$Q$})$
$\displaystyle\equiv\sum_{m_{\lambda}m_{\sigma}}\left\langle\lambda
m_{\lambda},\sigma m_{\sigma}\middle|\xi m_{\xi}\right\rangle$
$\displaystyle\hskip
28.45274pt\times\int\mathrm{d}\theta^{\prime}\int\mathrm{d}\phi^{\prime}\int\mathrm{d}\theta\int\mathrm{d}\phi\mathcal{O}_{(l^{\prime}l)\lambda
m_{\lambda}}^{(L)}(\theta^{\prime},\phi^{\prime},\theta,\phi)\mathcal{O}_{(s^{\prime}s)\sigma
m_{\sigma},m_{\chi}^{\prime}m_{\chi}}^{(S)}(\mbox{\boldmath$p$}^{\prime},\mbox{\boldmath$p$};\mbox{\boldmath$Q$})\,,$
(56)

the angular momentum-dependent part is a function of spherical harmonics while
the spin-dependent part evaluates the remaining matrix-element

$\displaystyle\mathcal{O}_{(l^{\prime}l)\lambda
m_{\lambda}}^{(L)}(\theta^{\prime},\phi^{\prime},\theta,\phi)\equiv\frac{\hat{\lambda}}{\hat{l}^{\prime}}\sum_{m_{l}^{\prime}m_{l}}\left\langle
lm_{l},\lambda m_{\lambda}\middle|l^{\prime}m_{l}^{\prime}\right\rangle
Y_{l^{\prime}m_{l}^{\prime}}^{*}(\theta^{\prime},\phi^{\prime})Y_{lm_{l}}(\theta,\phi)$
(57)

and

$\displaystyle\mathcal{O}_{(s^{\prime}s)\sigma
m_{\sigma},m_{\chi}^{\prime}m_{\chi}}^{(S)}(\mbox{\boldmath$p$}^{\prime},\mbox{\boldmath$p$};\mbox{\boldmath$Q$})\equiv\frac{\hat{\sigma}}{\hat{s}^{\prime}}\sum_{m_{s}^{\prime}m_{s}}\left\langle
sm_{s},\sigma m_{\sigma}\middle|s^{\prime}m_{s}^{\prime}\right\rangle\langle
s^{\prime}m_{s}^{\prime}m_{\chi}^{\prime}|O(\mbox{\boldmath$p$}^{\prime},\mbox{\boldmath$p$};\mbox{\boldmath$Q$})|sm_{s}m_{\chi}\rangle$
(58)

Conservation laws which simplify the angular-integration can be applied
depending on the spin matrix-elements of the operator
$\mathcal{O}_{(s^{\prime}s)\sigma
m_{\sigma},m_{\chi}^{\prime}m_{\chi}}^{(S)}$.

###  B.1 Isoscalar two-body correction

For two body operators, we can write the Jacobi momenta in terms of nucleon
momentum $\mbox{\boldmath$k$}_{i},\mbox{\boldmath$k$}_{i}^{\prime}$ as

$\displaystyle\mbox{\boldmath$p$}_{12}$
$\displaystyle=\frac{\mbox{\boldmath$k$}_{1}-\mbox{\boldmath$k$}_{2}}{2}\,,$
$\displaystyle\mbox{\boldmath$p$}_{12}^{\prime}$
$\displaystyle=\frac{\mbox{\boldmath$k$}_{1}^{\prime}-\mbox{\boldmath$k$}_{2}^{\prime}}{2}\,,$
(59) $\displaystyle\mbox{\boldmath$q$}_{1}$
$\displaystyle=\frac{\mbox{\boldmath$Q$}}{2}+\mbox{\boldmath$p$}_{12}-\mbox{\boldmath$p$}_{12}^{\prime}\,,$
$\displaystyle\mbox{\boldmath$q$}_{2}$
$\displaystyle=\frac{\mbox{\boldmath$Q$}}{2}-\mbox{\boldmath$p$}_{12}+\mbox{\boldmath$p$}_{12}^{\prime}\,.$
(60)

For scalar dark matter interactions, the spin matrix element in (58) is
proportional to
$\left|\mbox{\boldmath$k$}\right|\equiv\left|\mbox{\boldmath$p$}_{12}+\tilde{\alpha}\mbox{\boldmath$p$}_{12}^{\prime}+\tilde{\beta}\mbox{\boldmath$Q$}\right|$
and $m_{\xi}=0$. In the structure $\left|\mbox{\boldmath$k$}\right|$,
azimuthal angles appears as $\phi_{p_{12}}-\phi_{p_{12}^{\prime}}$ (using
$\mbox{\boldmath$Q$}=Q\boldsymbol{e}_{z}$), it is convenient to use the
following transformation

$\displaystyle\int\mathrm{d}\phi_{p_{12}}\int\mathrm{d}\phi_{p_{12}^{\prime}}\rightarrow\int\mathrm{d}\varphi\int\mathrm{d}\Phi\,,$
(61)

with
$\varphi=\phi_{p_{12}}-\phi_{p_{12}^{\prime}},\Phi=\tfrac{\phi_{p_{12}}+\phi_{p_{12}^{\prime}}}{2}$.
This transformation modifies the spherical harmonics in Eq. (57) as

$\displaystyle
Y_{l^{\prime}m_{l}^{\prime}}^{*}(\theta^{\prime},\phi^{\prime})Y_{lm_{l}}(\theta,\phi)\rightarrow
Y_{l^{\prime}m_{l}^{\prime}}^{*}(\theta^{\prime},0)Y_{lm_{l}}(\theta,0)e^{im_{l}\varphi}e^{im_{\lambda}\left(\frac{\varphi}{2}-\Phi\right)}$
(62)

where we used $m_{l}+m_{\lambda}=m_{l}^{\prime}$ and $m_{\xi}=0$, thereby
$m_{\lambda}=-m_{\sigma}$. This allows us to absorb one integration in the
spin matrix-element

$\displaystyle\mathcal{O}_{((l^{\prime}l)\lambda,(s^{\prime}s)\sigma)\xi
m_{\xi},m_{\chi}^{\prime}m_{\chi}}(p_{12}^{\prime},p_{12},Q)$
$\displaystyle\hskip
28.45274pt\equiv\sum_{m_{\lambda}m_{\sigma}}\left\langle\lambda
m_{\lambda},\sigma m_{\sigma}\middle|\xi m_{\xi}\right\rangle$
$\displaystyle\hskip
56.9055pt\times\int\mathrm{d}\theta^{\prime}\int\mathrm{d}\theta\int\mathrm{d}\varphi\tilde{\mathcal{O}}_{(l^{\prime}l)\lambda
m_{\lambda}}^{(L)}(\theta_{p_{12}^{\prime}},\theta_{p_{12}},\varphi)\tilde{\mathcal{O}}_{(s^{\prime}s)\sigma
m_{\sigma},m_{\chi}^{\prime}m_{\chi}}^{(S)}(p_{12}^{\prime},p_{12},Q,\varphi)\,,$
(63)

with

$\displaystyle\tilde{\mathcal{O}}_{(l^{\prime}l)\lambda
m_{\lambda}}^{(L)}(\theta_{p_{12}^{\prime}},\theta_{p_{12}},\varphi)\equiv\frac{\hat{\lambda}}{\hat{l}^{\prime}}\sum_{m_{l}^{\prime}m_{l}}\left\langle
lm_{l},\lambda m_{\lambda}\middle|l^{\prime}m_{l}^{\prime}\right\rangle
Y_{l^{\prime}m_{l}^{\prime}}^{*}(\theta_{p_{12}^{\prime}},0)Y_{lm_{l}}(\theta_{p_{12}},0)e^{im_{l}\varphi}\,,$
(64) $\displaystyle\tilde{\mathcal{O}}_{(s^{\prime}s)\sigma
m_{\sigma},m_{\chi}^{\prime}m_{\chi}}^{(S)}(p_{12}^{\prime},p_{12},Q,\varphi)$
$\displaystyle\qquad\equiv\frac{\hat{\sigma}}{\hat{s}^{\prime}}\sum_{m_{s}^{\prime}m_{s}}\left\langle
sm_{s},\sigma
m_{\sigma}\middle|s^{\prime}m_{s}^{\prime}\right\rangle\int\mathrm{d}\Phi\langle
s^{\prime}m_{s}^{\prime}m_{\chi}^{\prime}|O(\mbox{\boldmath$p$}_{12}^{\prime},\mbox{\boldmath$p$}_{12},\mbox{\boldmath$Q$})|sm_{s}m_{\chi}\rangle
e^{im_{\sigma}(\Phi-\tfrac{\varphi}{2})}\,.$ (65)

For scalar dark matter scattering, the last expression is either $\Phi$
independent or easy enough to perform the integration analytically, this
significantly reduces the numerical scaling and is eventually used for
computation.

$s^{\prime}$ | $s$ | $\sigma$ | $m_{\sigma}$ | $\tilde{\mathcal{O}}_{(s^{\prime}s)\sigma m_{\sigma},m_{\chi}^{\prime}m_{\chi}}^{(S)}(p^{\prime},p^{\prime}q,\varphi)$
---|---|---|---|---
0 | 0 | 0 | 0 | -2$\pi(\alpha-2\gamma)$
0 | 1 | 1 | -1 | $\sqrt{6}\pi\beta Q$
0 | 1 | 1 | 1 | $\sqrt{6}\pi\beta^{*}Q$
1 | 0 | 1 | -1 | $\sqrt{2}\pi\beta Q$
1 | 0 | 1 | 1 | $\sqrt{2}\pi\beta^{*}Q$
1 | 1 | 0 | 0 | $\frac{2}{3}\pi(\alpha-2\gamma)$
1 | 1 | 2 | -2 | $2\sqrt{\frac{5}{3}}\pi\beta^{2}$
1 | 1 | 2 | -1 | $4\sqrt{\frac{5}{3}}\pi\beta\delta$
1 | 1 | 2 | 0 | $\frac{1}{3}\sqrt{10}\pi(\alpha+3\delta^{2}-3\left|\beta\right|^{2}-2\gamma-\frac{3}{4}Q^{2})$
1 | 1 | 2 | 1 | $-4\sqrt{\frac{5}{3}}\pi\beta^{*}\delta$
1 | 1 | 2 | 2 | $2\sqrt{\frac{5}{3}}\pi(\beta^{*})^{2}$
Table 2: Partial wave decomposition of spin operator in Eq. (66) for non-zero
channels. The coefficients $\alpha,\beta,\gamma,\delta$ are defined in Eq.
(67).

Here we perform the above decomposition for Diagram 1(c) given by Eq. (7). For
simplicity, we are omitting the overall
$m_{\pi}^{2}\left(\frac{g_{A}}{2f_{\pi}}\right)^{2}$ factor and the $\Phi$
independent denominator
$(\mbox{\boldmath$q$}_{1}^{2}+m_{\pi}^{2})(\mbox{\boldmath$q$}_{2}^{2}+m_{\pi}^{2})$.
The spin-momentum decomposition for this operator

$\displaystyle\tilde{\mathcal{O}}_{(s^{\prime}s)\sigma
m_{\sigma},m_{\chi}^{\prime}m_{\chi}}^{(S)}(p_{12}^{\prime},p_{12},Q,\varphi)$
$\displaystyle\qquad=\frac{\hat{\sigma}}{\hat{s}^{\prime}}\sum_{m_{s}^{\prime}m_{s}}\left\langle
sm_{s},\sigma
m_{\sigma}\middle|s^{\prime}m_{s}^{\prime}\right\rangle\int\mathrm{d}\Phi\langle
s^{\prime}m_{s}^{\prime}|-(\mbox{\boldmath$\sigma$}_{1}\cdot\mbox{\boldmath$q$}_{1})(\mbox{\boldmath$\sigma$}_{2}\cdot\mbox{\boldmath$q$}_{2})|sm_{s}\rangle
e^{im_{\sigma}(\Phi-\tfrac{\varphi}{2})}\,$ (66)

is given in Table 2. The coefficient used in the table are

$\displaystyle\alpha$ $\displaystyle=p_{12}^{2}+p_{12}^{\prime
2}-\left(\frac{Q}{2}\right)^{2}\,,$ $\displaystyle\beta$
$\displaystyle=e^{i\varphi}p_{12}\sin(\theta_{p_{12}})-p_{12}^{\prime}\sin\theta_{p_{12}^{\prime}}\,,$
$\displaystyle\gamma$
$\displaystyle=p_{12}p_{12}^{\prime}(\cos\theta_{p_{12}}\cos\theta_{p_{12}^{\prime}}+\cos\varphi\sin\theta_{p_{12}}\sin\theta_{p_{12}^{\prime}})\,,$
$\displaystyle\delta$
$\displaystyle=p_{12}\cos\theta_{p_{12}}-p_{12}^{\prime}\cos\theta_{p_{12}^{\prime}}\,.$
(67)

The momenta $\mbox{\boldmath$q$}_{1,2}$ and the denominator are given by

$\displaystyle\mbox{\boldmath$q$}_{1,2}^{2}=\alpha-2\gamma\pm\delta
Q+\frac{1}{2}Q^{2}\Rightarrow(\mbox{\boldmath$q$}_{1}^{2}+m_{\pi}^{2})(\mbox{\boldmath$q$}_{2}^{2}+m_{\pi}^{2})=\left(\alpha-2\gamma+\frac{1}{2}Q^{2}+m_{\pi}^{2}\right)^{2}-\delta^{2}Q^{2}\,.$
(68)

This expression is spin independent. All these structures are combined and
convoluted with the density matrix which is used to numerically calculate the
scattering matrix elements.

## References

  * [1] G. Bertone and D. Hooper,  Rev. Mod. Phys. 90, 045002 (2018), 1605.04909.
  * [2] M. Battaglieri et al.,  US Cosmic Visions: New Ideas in Dark Matter 2017: Community Report,  in U.S. Cosmic Visions: New Ideas in Dark Matter, 2017, 1707.04591.
  * [3] S. A. Hertel, A. Biekert, J. Lin, V. Velan, and D. N. McKinsey,  Phys. Rev. D 100, 092007 (2019), 1810.06283.
  * [4] M. Schumann,  J. Phys. G 46, 103003 (2019), 1903.03026.
  * [5] J. Fan, M. Reece, and L.-T. Wang,  JCAP 11, 042 (2010), 1008.1591.
  * [6] A. L. Fitzpatrick, W. Haxton, E. Katz, N. Lubbers, and Y. Xu,  JCAP 02, 004 (2013), 1203.3542.
  * [7] V. Cirigliano, M. L. Graesser, and G. Ovanesyan,  JHEP 10, 025 (2012), 1205.2695.
  * [8] P. Klos, J. Menéndez, D. Gazit, and A. Schwenk,  Phys. Rev. D 88, 083516 (2013), 1304.7684,  [Erratum: Phys.Rev.D 89, 029901 (2014)].
  * [9] M. Hoferichter, P. Klos, and A. Schwenk,  Phys. Lett. B 746, 410 (2015), 1503.04811.
  * [10] J. Brod, A. Gootjes-Dreesbach, M. Tammaro, and J. Zupan,  JHEP 10, 065 (2018), 1710.10218.
  * [11] T. R. Richardson, X. Lin, and S. T. Nguyen,  Phys. Rev. C 106, 044003 (2022), 2110.15959.
  * [12] F. D’Eramo, B. J. Kavanagh, and P. Panci,  JHEP 08, 111 (2016), 1605.04917.
  * [13] F. Bishara, J. Brod, B. Grinstein, and J. Zupan,  JHEP 03, 089 (2020), 1809.03506.
  * [14] G. Prezeau, A. Kurylov, M. Kamionkowski, and P. Vogel,  Phys. Rev. Lett. 91, 231301 (2003), astro-ph/0309115.
  * [15] C. Körber, A. Nogga, and J. de Vries,  Phys. Rev. C 96, 035805 (2017), 1704.01150.
  * [16] M. Hoferichter, P. Klos, J. Menéndez, and A. Schwenk,  Phys. Rev. D 94, 063505 (2016), 1605.08043.
  * [17] S. R. Beane, S. D. Cohen, W. Detmold, H. W. Lin, and M. J. Savage,  Phys. Rev. D 89, 074505 (2014), 1306.6939.
  * [18] L. Andreoli, V. Cirigliano, S. Gandolfi, and F. Pederiva,  Phys. Rev. C 99, 025501 (2019), 1811.01843.
  * [19] V. Cirigliano et al.,  Phys. Rev. Lett. 120, 202001 (2018), 1802.10097.
  * [20] SPICE/HeRALD, A. Biekert et al.,  Phys. Rev. D 105, 092005 (2022), 2108.02176.
  * [21] J. Liao et al.,  (2021), 2103.02161.
  * [22] G. Arcadi, A. Djouadi, and M. Raidal,  Phys. Rept. 842, 1 (2020), 1903.03616.
  * [23] B. S. Hu et al.,  Phys. Rev. Lett. 128, 072502 (2022), 2109.00193.
  * [24] Particle Data Group, R. L. Workman et al.,  PTEP 2022, 083C01 (2022).
  * [25] D. A. Brantley et al.,  (2016), 1612.07733.
  * [26] M. Hoferichter, J. Ruiz de Elvira, B. Kubis, and U.-G. Meißner,  Phys. Rev. Lett. 115, 092301 (2015), 1506.04142.
  * [27] R. Gupta et al.,  Phys. Rev. Lett. 127, 242002 (2021), 2105.12095.
  * [28] M. Hoferichter, J. R. de Elvira, B. Kubis, and U.-G. Meißner,  Phys. Lett. B 843, 138001 (2023), 2305.07045.
  * [29] S. Borsanyi et al.,  (2020), 2007.03319.
  * [30] H. Krebs, E. Epelbaum, and U.-G. Meißner,  Eur. Phys. J. A 56, 240 (2020), 2005.07433.
  * [31] H. Krebs, E. Epelbaum, and U. G. Meißner,  Few Body Syst. 60, 31 (2019), 1902.06839.
  * [32] A. Nogga, H. Kamada, W. Gloeckle, and B. R. Barrett,  Phys. Rev. C 65, 054003 (2002), nucl-th/0112026.
  * [33] R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla,  Phys. Rev. C 51, 38 (1995), nucl-th/9408016.
  * [34] R. Machleidt,  Phys. Rev. C 63, 024001 (2001), nucl-th/0006014.
  * [35] E. Epelbaum, H. Krebs, and U. G. Meißner,  Phys. Rev. Lett. 115, 122301 (2015), 1412.4623.
  * [36] P. Reinert, H. Krebs, and E. Epelbaum,  Eur. Phys. J. A 54, 86 (2018), 1711.08821.
  * [37] E. Epelbaum, H. Krebs, and U. G. Meißner,  Eur. Phys. J. A 51, 53 (2015), 1412.0142.
  * [38] H. Krebs,  Eur. Phys. J. A 56, 234 (2020), 2008.00974.
  * [39] A. A. Filin et al.,  Phys. Rev. C 103, 024313 (2021), 2009.08911.
  * [40] R. J. Furnstahl, N. Klco, D. R. Phillips, and S. Wesolowski,  Phys. Rev. C 92, 024005 (2015), 1506.01343.
  * [41] J. A. Melendez, R. J. Furnstahl, D. R. Phillips, M. T. Pratola, and S. Wesolowski,  Phys. Rev. C 100, 044001 (2019), 1904.10581.
  * [42] H. W. Grießhammer, J. A. McGovern, A. Nogga, and D. R. Phillips,  Few Body Syst. 61, 48 (2020), 2005.12207.
  * [43] A. Nogga,  NUCDENS: A package to access momentum space densities for light nuclei.; v1.1.0, 2023.
  * [44] M. C. Smith et al.,  Mon. Not. Roy. Astron. Soc. 379, 755 (2007), astro-ph/0611671.
  * [45] V. Cirigliano, W. Dekens, J. de Vries, M. Hoferichter, and E. Mereghetti,  Phys. Rev. Lett. 126, 172002 (2021), 2012.11602.
  * [46] R. Wirth, J. M. Yao, and H. Hergert,  Phys. Rev. Lett. 127, 242502 (2021), 2105.05415.
  * [47] L. Jokiniemi, P. Soriano, and J. Menéndez,  Phys. Lett. B 823, 136720 (2021), 2107.13354.
  * [48] M. Hoferichter, P. Klos, J. Menéndez, and A. Schwenk,  Phys. Rev. D 99, 055031 (2019), 1812.05617.
  * [49] NPLQCD, E. Chang et al.,  Phys. Rev. Lett. 120, 152002 (2018), 1712.03221.
  * [50] J. L. Friar,  Phys. Rev. C 20, 325 (1979).
  * [51] S. Liebig, U. G. Meißner, and A. Nogga,  Eur. Phys. J. A 52, 103 (2016), 1510.06070.
  * [52] H. Le, J. Haidenbauer, U.-G. Meißner, and A. Nogga,  Eur. Phys. J. A 56, 301 (2020), 2008.11565.