The ArgoNeuT Collaboration

# First Constraints on Heavy QCD Axions with a Liquid Argon Time Projection
Chamber using the ArgoNeuT Experiment

R. Acciarri  Fermi National Accelerator Laboratory, Batavia, Illinois 60510,
USA  C. Adams  Argonne National Laboratory, Lemont, Illinois 60439, USA  B.
Baller  Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA
V. Basque  Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA
F. Cavanna  Fermi National Accelerator Laboratory, Batavia, Illinois 60510,
USA  R. T. Co  School of Physics and Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  William I. Fine Theoretical Physics Institute,
University of Minnesota, Minneapolis, MN 55455, USA  R. S. Fitzpatrick
University of Michigan, Ann Arbor, Michigan 48109, USA  B. Fleming  Yale
University, New Haven, Connecticut 06520, USA  P. Green
<EMAIL_ADDRESS>University of Manchester, Manchester M13 9PL,
United Kingdom  University of Oxford, Oxford OX1 3RH, United Kingdom  R.
Harnik  Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA
K. J. Kelly  CERN, Esplande des Particules, 1211 Geneva 23, Switzerland  S.
Kumar  University of California, Berkeley, California 94720, USA  Lawrence
Berkeley National Laboratory, Berkeley, California 94720, USA  K. Lang
University of Texas at Austin, Austin, Texas 78712, USA  I. Lepetic  Rutgers
University, Piscataway, New Jersey 08854, USA  Z. Liu<EMAIL_ADDRESS>School
of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
X. Luo  University of California, Santa Barbara, California, 93106, USA  K. F.
Lyu  School of Physics and Astronomy, University of Minnesota, Minneapolis, MN
55455, USA  O. Palamara  Fermi National Accelerator Laboratory, Batavia,
Illinois 60510, USA  G. Scanavini  Yale University, New Haven, Connecticut
06520, USA  M. Soderberg  Syracuse University, Syracuse, New York 13244, USA
J. Spitz  University of Michigan, Ann Arbor, Michigan 48109, USA  A. M. Szelc
University of Edinburgh, Edinburgh EH9 3FD, United Kingdom  W. Wu  Fermi
National Accelerator Laboratory, Batavia, Illinois 60510, USA  T. Yang  Fermi
National Accelerator Laboratory, Batavia, Illinois 60510, USA

###### Abstract

We present the results of a search for heavy QCD axions performed by the
ArgoNeuT experiment at Fermilab. We search for heavy axions produced in the
NuMI neutrino beam target and absorber decaying into dimuon pairs, which can
be identified using the unique capabilities of ArgoNeuT and the MINOS near
detector. This decay channel is motivated by a broad class of heavy QCD axion
models that address the strong CP and axion quality problems with axion masses
above the dimuon threshold. We obtain new constraints at a 95% confidence
level for heavy axions in the previously unexplored mass range between 0.2-0.9
GeV, for axion decay constants around tens of TeV.

††preprint: FERMILAB-PUB-22-527-ND-T, UMN-TH-4128/22, FTPI-MINN-22/19

Introduction.—The QCD axion was proposed [1, 2] to address the strong CP
problem [3, 4, 5, 6, 7]. However, in the simplest implementations, this
mechanism suffers from the axion quality problem [8, 9, 10]. Heavy QCD axions,
defined as those with a coupling to gluons but with a much larger mass than
the QCD axion, are motivated by their potential role resolving the axion
quality problem while preserving the solution to the strong CP problem [11,
12, 13, 14, 15, 16, 17, 18]. Furthermore, they can explain various phenomena
in astrophysics [19] and cosmology [20, 21]. In these models, the larger axion
mass, $m_{a}$, and smaller axion decay constant, $f_{a}$, also open up various
decay channels involving Standard Model (SM) particles and enhance the axion
interaction strengths. This enables searches for these heavy QCD axions in
beam-dump and collider experiments.

In this _Letter_ , we perform a search for heavy QCD axions with
$200~{}\text{MeV}\lesssim m_{a}\lesssim 1~{}\text{GeV}$ using the ArgoNeuT
experiment [22] – a 0.24 ton Liquid Argon Time Projection Chamber (LArTPC)
neutrino detector that collected five months of data in 2009-2010 in the
Neutrinos at the Main Injector (NuMI) beamline [23] at Fermilab. The data used
corresponds to $1.25\text{\times}{10}^{20}$ protons-on-target (POT) acquired
while the NuMI beam was in anti-neutrino mode [22]. The axions can be produced
via couplings with SM mesons and protons when the 120 GeV proton beam strikes
the graphite target or the hadron absorber located 1033 m and 318 m upstream
of ArgoNeuT, respectively. The produced axions can then propagate to ArgoNeuT
where the decay signature $a\to\mu^{+}\mu^{-}$ is searched for. The muon pair
is reconstructed in ArgoNeuT as minimally ionizing tracks that can then be
matched to a pair of tracks with opposite charges in the magnetised MINOS near
detector (MINOS-ND) [24] located immediately downstream of ArgoNeuT.

Heavy QCD axions.— Heavy QCD axions must couple to gluons to solve the strong
CP problem. Furthermore, consistent with Grand Unified Theories, these axions
can also couple to the other gauge bosons of the SM. These considerations
motivate the following couplings [25],

$\displaystyle\mathcal{L}_{\rm gauge}=\frac{c_{3}\alpha_{3}}{8\pi
f_{a}}aG\tilde{G}+\frac{c_{2}\alpha_{2}}{8\pi
f_{a}}aW\tilde{W}+\frac{c_{1}\alpha_{1}}{8\pi f_{a}}aB\tilde{B}.$ (1)

Here $G$ is SM gluon field strength and $\tilde{G}$ is its dual. Couplings to
$SU(2)$ and $U(1)$ gauge fields, $W\tilde{W}$ and $B\tilde{B}$, are defined
analogously. The coefficients $\alpha_{i}=g_{i}^{2}/(4\pi)$ are given in terms
of the three gauge couplings $g_{i}$ defined at the scale $m_{a}$. We will set
$c_{3}=c_{2}=c_{1}=1$ hereafter. We note, however, that the results presented
in this _Letter_ primarily depend on the $aG\tilde{G}$ coupling and would
still apply if $c_{1},c_{2}\ll c_{3}$.

Along with couplings to gauge bosons, axions can also couple to SM fermions,
as in the DFSZ models [26, 27]. While coupling to both quarks and leptons can
appear, we consider axion coupling to only SM leptons in a flavor diagonal way
in order to focus on a parameter space that is complementary to the multitude
of flavor searches (see, e.g., [28] for a recent summary and references) and
is theoretically well motivated. Therefore, we consider [25]

$\displaystyle\mathcal{L}_{\rm
lepton}=\sum_{\ell=e,\mu,\tau}\frac{\partial_{\mu}a}{2f_{a}}\left(c_{V\ell}\bar{\ell}\gamma^{\mu}\ell+c_{A\ell}\bar{\ell}\gamma^{\mu}\gamma_{5}\ell\right).$
(2)

Here $c_{V\ell},c_{A\ell}$ control the flavor universal vector and axial
coupling of the axion to SM charged leptons.

Given the axion couplings in Eqs. (1) and (2), the relevant decay modes of the
axions for masses $0.2~{}\text{GeV}<m_{a}<1~{}\text{GeV}$ are into photons,
muons and some exclusive hadronic modes. Considering these decay modes, we
show the behaviour of the decay length of the axion in its rest frame for
$f_{a}=10$ TeV in Fig. 1 (top). Discussion of the individual contributions of
the various decay channels involved can be found in the Supplemental Material
[29]. Hereafter, we focus on two theory benchmarks with $c_{\ell}\approx
c_{A\ell}=1/36$ and $c_{\ell}\approx c_{A\ell}=1/100$ to illustrate the
constraining power in the $f_{a}$-$m_{a}$ plane. Choosing $c_{\ell}$ smaller
than $c_{i}=1$ is motivated because, from a theoretical perspective, a
suppressed leptonic coupling can naturally emerge in models where axions
directly couple to some new heavy leptons, which in turn mix with SM leptons
giving $c_{\ell}\propto\theta_{\rm mix}^{2}\ll 1$ [30, 31, 32]. Our choice of
$c_{\ell}=1/36$ and $c_{\ell}=1/100$ then corresponds to small mixing angle
benchmarks $\theta_{\rm mix}\approx 1/6$ and $\theta_{\rm mix}\approx 1/10$,
respectively, that can arise from loop-induced or vector-like fermion-induced
models [33]. From the perspective of axion searches in ArgoNeuT, a smaller
$c_{\ell}<c_{i}$ also makes the axion sufficiently long-lived to reach the
detector while not suppressing its production via the gluon coupling.

Figure 1: Lifetime (top) and dimuon branching ratio (bottom) of the axion for
the two benchmark scenarios as a function of the axion mass.

Since our search is based on muon final states, we also show in Fig. 1
(bottom) the branching ratio of the axion into two muons for the two
benchmarks. The dimuon mode is a dominant decay channel for most of the mass
range, enabling ArgoNeuT to be uniquely sensitive to these scenarios. In
regions where dimuon decays are subdominant, future searches in other
channels, e.g. $\gamma\gamma$ and multi-hadron states, could provide
complimentary coverage.

Generation and simulation.— Due to the presence of the gluon coupling, axions
mix with SM pseudoscalar mesons $\pi,\eta,\eta^{\prime}$ 111Other modes of
axion production, such as from meson decays, are subdominant unless for
example additional flavor-violating couplings are introduced. Such models are
not considered in this analysis. A recent work [61] focused on such a model,
but without a gluon coupling. with the corresponding mixing angles
$\theta_{aM}$, $M=\pi,\eta,\eta^{\prime}$, given by [25, 35, 36]

$\displaystyle\theta_{aM}=\dfrac{f_{\pi}}{f_{a}}\frac{C_{aM}}{m_{a}^{2}-m_{M}^{2}}.$
(3)

Here $f_{\pi}=93$ MeV is the pion decay constant and $C_{a\pi}=m_{a}^{2}/6$,
$C_{a\eta}=(m_{a}^{2}-4m_{\pi}^{2}/9)/\sqrt{6}$ and
$C_{a\eta^{\prime}}=(m_{a}^{2}-16m_{\pi}^{2}/9)/(2\sqrt{3})$. Our subsequent
analysis assumes that $\theta_{aM}\ll 1$, and therefore we mask the parameter
space for which $m_{a}\approx m_{M}$, as shown by vertical gray bands in Fig.
1. To compute the spectrum of axions produced in the NuMI beam we first
simulate the spectrum of SM mesons using Pythia8 [37, 38]. We find that on
average 2.89 $\pi^{0}$, 0.33 $\eta$ and 0.03 $\eta^{\prime}$ are produced per
proton collision. Subsequently, applying the mixing angles in Eq. (3), we can
compute the total number of axions produced considering both the 87% of
protons interacting in the NuMI target and the 10% that reach the downstream
hadron absorber with energies $\sim$120 GeV [23]. Using the geometrical
acceptance of ArgoNeuT and the axion branching ratio into muons (Fig. 1), we
then compute the axion decays to a dimuon final state that would be seen in
the detector. The LArSoft software framework [39] is used to simulate these
muons in ArgoNeuT. LArSoft propagates the particles using GEANT4 [40] and then
performs detector response modelling and reconstruction [22, 41]. The standard
MINOS simulation and reconstruction framework is then used to model the
particles that exit ArgoNeuT and reach the MINOS-ND [22, 24].

Signature and selection.— In ArgoNeuT the axion decay
$a\rightarrow\mathit{{\mu}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{+}}}\mathit{{\mu}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{-}}}$
can be detected as a pair of minimally ionizing particles (MIPs). The parent
axion energies and the kinematics of the resulting muons for two different
axion masses are shown in Fig. 2. The muons are highly energetic, with
$\langle E_{\mu^{\pm}}\rangle\approx 20$ GeV, resulting in them typically
exiting ArgoNeuT and propagating to the downstream MINOS-ND. They are also
highly forward-going with an average angle with respect to the beam direction
of $\langle\theta_{\rm beam}\rangle\approx 0.75\degree$ to $\langle\theta_{\rm
beam}\rangle\approx 2.5\degree$ and an average opening angle between them of
$\langle\theta_{\rm opening}\rangle\approx 1.5\degree$ to $\langle\theta_{\rm
opening}\rangle\approx 5.0\degree$, in each case depending on $m_{a}$. As a
result of ArgoNeuT’s angular reconstruction resolution of between $1\degree$
and $3\degree$ (depending on track orientation) [42], the muons may overlap
and appear as one track in ArgoNeuT. Once reaching the MINOS-ND, the muons
then separate due to their opposite charges in the MINOS-ND magnetic field of
about 1 T [24].

Figure 2: Parent axion energy (top), energy of the resulting individual muons
(middle) and opening angle between them (bottom) for axions with $m_{a}=360$
MeV (black) and $m_{a}=720$ MeV (blue).

The axion decay signature is similar to the signature in ArgoNeuT’s previous
heavy neutral lepton search [43] and an analogous selection strategy can be
followed. However, it is a resonant two-body decay as opposed to the non-
resonant three-body decay previously considered. Therefore, given the
different nature of this decay and the different production mechanism and
resulting kinematics, the selection has been modified and re-optimised. We
consider two different scenarios, depending on how forward-going the muon pair
are. In the first, the pair of muons are sufficiently separated to be
reconstructed as two distinct minimally ionizing tracks. In the second, the
separation between them is too small and they are instead reconstructed as a
single track with twice the minimally ionizing particle $dE/dx$. These will be
referred to subsequently as two-track and double-MIP type events,
respectively. Due to the highly forward-going nature of the muon pairs
resulting from axion decays, the double-MIP signature dominates. In both
cases, the tracks in ArgoNeuT can be matched with two tracks in the downstream
MINOS-ND that are reconstructed with opposite charges. In the second scenario,
axion decays occurring in the upstream cavern are also considered. The
resulting muons then pass through the ArgoNeuT detector and can be matched to
the MINOS-ND as before. The ArgoNeuT physics run coincided with the
construction of the upstream MINERvA detector [44]. Therefore, only the 63 cm
region between MINERvA and ArgoNeuT is taken into account.

A series of pre-selection cuts are first applied. Events with incomplete
reconstruction are removed by requiring that at least 80% of energy
depositions are associated with reconstructed objects. Next, events are
identified that are compatible with axion decays: requiring a maximum of two
tracks in ArgoNeuT, where short tracks ($L\leq 5$ cm) are ignored to ensure
events are not removed due to the presence of $\delta$-rays, and at least two
tracks in the MINOS-ND. Reconstructed tracks are also classified based on
whether they originate from within or outside of the ArgoNeuT fiducial volume:
$1\leq x\leq 46$ cm (drift), $-19\leq y\leq 19$ cm (vertical) and $z\geq 3$ cm
(beam direction).

Two selection paths are then followed in ArgoNeuT depending on the event
topology. In the two-track case events are selected that have two minimally-
ionising tracks originating from a common vertex within the ArgoNeuT fiducial
volume. Tracks are considered to be minimally-ionising if they have an average
$dE/dx<3.1$ MeV/cm and are considered to have originated from a common vertex
if they have a separation of less than 4 cm between the track starts. The pair
of tracks are also required to have an opening angle between them of
$\theta_{opening}\leq 15\degree$ and exit towards the MINOS-ND. In the double-
MIP case events are selected that have a single track, originating from either
within ArgoNeuT or from upstream of the detector, that has an average $dE/dx$
consistent with two overlapping minimally ionising particles. The region near
the track vertex prior to the potential separation of the two muons is
assessed: considering the average $dE/dx$ over the first 5 cm, where any high
energy depositions ($dE/dx>10$ MeV/cm) from $\delta$-rays are removed. The
track is required to have an average $dE/dx>3.1$ MeV/cm. Additionally, the
track is required to be forward-going with an angle with respect to the beam
direction $\theta_{beam}\leq 10\degree$ and exit towards the MINOS-ND.

Next, matching to the MINOS-ND is performed. In the two-track scenario each of
the two tracks are individually matched to tracks reconstructed in the MINOS-
ND, with matching tolerances of $r_{\rm diff}\leq 12.0$ cm (radial) and
$\theta_{\rm diff}\leq 0.17$ rad (angular) accounting for the expected maximum
deflection of the muons between the two detectors [41]. The matched tracks are
also both required to start within the calorimeter region of the MINOS-ND,
which is directly downstream of the ArgoNeuT detector, and start within $20$
cm of the first instrumented plane. In the double-MIP scenario, since a single
track in ArgoNeuT now being matched to two tracks in the MINOS-ND, the
matching tolerances are loosened to double the two-track case. Several further
selection cuts are applied to the matched tracks in the MINOS-ND. They are
required to both be muon-like based both on the track length and $dE/dx$:
$L_{MINOS-ND}\geq 1$ m and $6\leq dE/dx_{MINOS-ND}\leq 12$ MeV/cm. They are
also required to have timestamps consistent with having originated from a
single decay in ArgoNeuT, $|\Delta t_{0}|<20$ ns, and be reconstructed with
opposite charges.

The selection efficiency is shown in Fig. 3 for simulated decays of axions
with mass $m_{a}=500$ MeV as a function of the axion energy, $E_{a}$. Decays
occurring inside ArgoNeuT (black) and at 25 cm (blue) and 50 cm (red) upstream
of the detector are shown. For decays inside the detector the efficiency is
around 50% above $E_{a}\sim 15$ GeV. However, at lower $E_{a}$ one or both of
the muons may have insufficient energy to propagate to the MINOS-ND before
stopping. This causes the matching to fail resulting in a sharp decline in the
selection efficiency. For decays upstream of the detector, where only the
double-MIP topology is considered, the selection efficiency is lower. As the
distance the muons have to propagate before reaching the detector increases,
the less likely they are to remain overlapping resulting in the events being
rejected. At lower $E_{a}$ the muons are less boosted, further increasing the
likelihood they separate.

Figure 3: Selection efficiency for $m_{a}=500$ MeV axion decays as a function
of axion energy, $E_{a}$. Decays occurring inside ArgoNeuT (black) and at 25
cm (blue) and 50 cm (red) upstream of the detector are shown.

For selected axion candidate events the invariant mass of the parent axion can
be reconstructed. This is achieved using the trajectories of the tracks
reconstructed in ArgoNeuT (or matching between ArgoNeuT and the MINOS-ND in
the double-MIP scenario), combined with the momentum reconstruction in the
MINOS-ND. The axion invariant mass can be reconstructed with a resolution of
$\sim$100 MeV when both muons are contained within the MINOS-ND, and $\sim$200
MeV if exiting. A constraint on the invariant mass is not applied in the
selection since the search is performed across a significant axion mass range.
However, if a signal were to be observed in the data the invariant mass could
be a powerful tool to further characterise it.

Backgrounds and systematic uncertainties.—

The primary backgrounds to this search originate from beam neutrino
interactions within the detector or the surrounding materials. These are
modelled using the GENIE neutrino generator [45] along with a data-driven
model of the rate and kinematics of beam-induced through-going muons [42, 46,
47]. Modelling of the NuMI beam flux developed by the MINERvA collaboration is
used [48]. The most significant backgrounds arise from charged-current muon
neutrino interactions: either single charged pion production, resulting in two
approximately minimally ionising tracks that could be mistaken for the two-
track signature; or quasi-elastic scattering producing a single muon with low
energy protons near the interaction vertex causing it to mimic the double-MIP
signature. In each case, these could then be incorrectly matched with other
nearby muons in the MINOS-ND. However, the majority of these types of
interactions are removed during the selection due to either the event topology
in ArgoNeuT or the MINOS-ND or as a result of the precise timing resolution of
the MINOS-ND [24]. The total background expectation for the ArgoNeuT data-set
is $0.1\pm 0.1$ events.

The systematic uncertainties are dominated by theoretical uncertainties in the
axion flux prediction. One way to estimate these is from the strong coupling,
$\alpha_{S}$, uncertainty of around $20\%$ [49]. Another way is from the
uncertainty on the measured pion flux of the NuMI beam that is at the 5-10%
level [50]. Furthermore, the production is modeled by chiral perturbation
theory where the matched effective coefficients also have uncertainties. The
current best fit for the pion coupling has an uncertainty of 10-20% [51]. To
conservatively illustrate the impact of these uncertainties, the theory
prediction of the axion flux is varied by 30% in the final sensitivity
evaluation. Future refinements in the axion flux modeling would improve the
accuracy of the derived limit. The impact of experimental systematic
uncertainties is also evaluated. These arise from uncertainties in the muon
reconstruction and resulting selection efficiency in ArgoNeuT (3.3%) [41, 42,
22, 52] and the MINOS-ND (0.4%) [53], the POT evaluation (1%) [47], and the
determination of the electron drift velocity and hence the total volume of
instrumented argon (2.2%) [42]. The experimental systematic uncertainties have
a combined impact of 4.1%.

Results.— Zero events pass the selection in ArgoNeuT’s full
$1.25\text{\times}{10}^{20}$ POT anti-neutrino-mode data-set, consistent with
the background prediction of $0.1\pm 0.1$ events. Our exclusion of parameter
space at a 95% confidence level is shown in Fig. 4 for both the
$c_{\ell}=1/36$ and $c_{\ell}=1/100$ scenarios, evaluated using a Bayesian
approach with a uniform prior [54]. The impact of the uncertainties on the
expected constraint, dominated by the theoretical uncertainty, is also shown
for the $c_{\ell}=1/36$ case 222The impact of the uncertainties on the
$c_{\ell}=1/100$ case is similar, but is not shown to aid clarity.. In the
presence of an axion there are new decay modes for SM mesons, such as
$B^{+}\rightarrow K^{+}a$, where $a$ can subsequently decay into $\mu\mu$.
Searches for such rare decay modes place important constraints on our
parameter space. We find the following searches give significant constraints:
$K^{+}\rightarrow\pi^{+}\nu\bar{\nu}$ by the NA62 collaboration [56];
$K^{\pm}\rightarrow\pi\mu\mu$ by the NA48/2 collaboration [57];
$B^{0}\rightarrow K^{*0}\mu\mu$ by the LHCb collaboration [58]; and
$B^{+}\rightarrow K^{+}\eta\pi\pi$ by the BaBar collaboration [59]. For all
these cases, we recast the presented bounds as appropriate for the axion
lifetime in our scenario. The strongest resulting constraints are shown in
Fig. 4 for each benchmark model.

This measurement leads to new constraints on previously unexplored parameter
space for heavy QCD axions with masses above the dimuon threshold and below 1
GeV (where hadronic decays would dominate). The coverage of the axion decay
constant for the benchmark model with $c_{\ell}=1/36$ is around $f_{a}\sim 50$
TeV for masses up to 0.65 GeV. The coverage of the axion decay constant for
the benchmark model with $c_{\ell}=1/100$ is around $f_{a}\sim 20$ TeV for
masses up to 0.84 GeV. For these benchmark couplings ArgoNeuT provides
significant improvement on existing constraints. ArgoNeuT also has
constraining power for couplings ranging between approximately
$c_{\ell}=\mathcal{O}(10^{-1}\mathchar 45\relax 10^{-3})$. For much larger
$c_{\ell}$ the axions would predominantly decay before reaching the detector
reducing the constraining power, and for much smaller $c_{\ell}$ the axions
would no longer dominantly decay to muons but to mesons and photons that are
not searched for in this analysis. To explore these regions of phase space
would require future searches with detectors at shorter baselines than
ArgoNeuT, or probing different decay modes.

Figure 4: Constraints on the axion model parameter space at 95% CL from the
exposure of $1.25\times 10^{20}$ POT at ArgoNeuT (blue shaded region and black
contours). The derived limits for $c_{\ell}=1/36$ and $c_{\ell}=1/100$ are
shown by the solid and dashed contours. The uncertainty on the expected
constraint, predominantly arising from the theoretical uncertainties, is shown
by the dark blue band. The red and orange contours show the strongest existing
constraints evaluated for the two benchmark scenarios. The gray-shaded band
indicates a region with increased theoretical uncertainty around the $\eta$
mass.

Conclusions.— We have presented the first search for heavy QCD axions in a
LArTPC using the ArgoNeuT experiment. This type of axion is particularly
motivated by the strong CP puzzle and the axion quality problem. We search for
such axions produced in the NuMI neutrino beam and then decaying with a dimuon
signature within, or close to, the ArgoNeuT detector. In ArgoNeuT’s
$1.25\text{\times}{10}^{20}$ POT anti-neutrino mode data-set zero events pass
the selection, consistent with the background prediction. This measurement
leads to a significant new exclusion region for heavy axions in the mass range
between 0.2-0.9 GeV for an axion decay constant around tens of TeV, over a
broad range of axion-lepton couplings around
$c_{\ell}=\mathcal{O}(10^{-1}\mathchar 45\relax 10^{-3})$. The search can be
extended to various new heavy QCD axion models and paves the way for heavy QCD
axion searches at future neutrino facilities. The techniques developed could
also be used to constrain other dark sector particle models with long-lived
resonance decays into dimuons.

Acknowledgements.— This manuscript has been authored by Fermi Research
Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of
Energy, Office of Science, Office of High Energy Physics. We gratefully
acknowledge the cooperation of the MINOS Collaboration in providing their data
for use in this analysis. We wish to acknowledge the support of Fermilab, the
Department of Energy, and the National Science Foundation in ArgoNeuT’s
construction, operation, and data analysis. This project has received funding
from the Science and Technology Facilities Council (STFC), part of the United
Kingdom Research and Innovation; and from the Royal Society UK grants:
RGF\EA\180209 and URF\R\201022. R.C. was supported in part by U.S. Department
of Energy (DOE) grant DE-SC0011842. S.K. was supported in part by the National
Science Foundation (NSF) grant PHY-1915314 and the DOE contract DE-
AC02-05CH11231. Z.L. and K.F.L. were supported in part by the DOE grant DE-
SC0022345. For the purpose of Open Access, the author has applied a CC BY
public copyright licence to any Author Accepted Manuscript version arising
from this submission.

## References

  * Peccei and Quinn [1977a] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn [1977b] R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977b).
  * ’t Hooft [1976] G. ’t Hooft, Phys. Rev. Lett. 37, 8 (1976).
  * Baker  _et al._ [2006] C. A. Baker _et al._ , Phys. Rev. Lett. 97, 131801 (2006), arXiv:hep-ex/0602020 [hep-ex] .
  * Pendlebury  _et al._ [2015] J. M. Pendlebury _et al._ , Phys. Rev. D 92, 092003 (2015), arXiv:1509.04411 [hep-ex] .
  * Graner  _et al._ [2016] B. Graner, Y. Chen, E. G. Lindahl, and B. R. Heckel, Phys. Rev. Lett. 116, 161601 (2016), [Erratum: Phys.Rev.Lett. 119, 119901 (2017)], arXiv:1601.04339 [physics.atom-ph] .
  * Abel  _et al._ [2020] C. Abel _et al._ (nEDM), Phys. Rev. Lett. 124, 081803 (2020), arXiv:2001.11966 [hep-ex] .
  * Giddings and Strominger [1988] S. B. Giddings and A. Strominger, Nucl. Phys. B307, 854 (1988).
  * Coleman [1988] S. R. Coleman, Nucl. Phys. B310, 643 (1988).
  * Gilbert [1989] G. Gilbert, Nucl. Phys. B328, 159 (1989).
  * Rubakov [1997] V. A. Rubakov, JETP Lett. 65, 621 (1997), arXiv:hep-ph/9703409 .
  * Berezhiani  _et al._ [2001] Z. Berezhiani, L. Gianfagna, and M. Giannotti, Phys. Lett. B 500, 286 (2001), arXiv:hep-ph/0009290 .
  * Hook [2015] A. Hook, Phys. Rev. Lett. 114, 141801 (2015), arXiv:1411.3325 [hep-ph] .
  * Fukuda  _et al._ [2015] H. Fukuda, K. Harigaya, M. Ibe, and T. T. Yanagida, Phys. Rev. D 92, 015021 (2015), arXiv:1504.06084 [hep-ph] .
  * Gherghetta  _et al._ [2016] T. Gherghetta, N. Nagata, and M. Shifman, Phys. Rev. D 93, 115010 (2016), arXiv:1604.01127 [hep-ph] .
  * Agrawal and Howe [2018] P. Agrawal and K. Howe, JHEP 12, 029 (2018), arXiv:1710.04213 [hep-ph] .
  * Hook  _et al._ [2020] A. Hook, S. Kumar, Z. Liu, and R. Sundrum, Phys. Rev. Lett. 124, 221801 (2020), arXiv:1911.12364 [hep-ph] .
  * Gherghetta  _et al._ [2020] T. Gherghetta, V. V. Khoze, A. Pomarol, and Y. Shirman, JHEP 03, 063 (2020), arXiv:2001.05610 [hep-ph] .
  * Nomura and Thaler [2009] Y. Nomura and J. Thaler, Phys. Rev. D 79, 075008 (2009), arXiv:0810.5397 [hep-ph] .
  * Takahashi and Yin [2021] F. Takahashi and W. Yin, JCAP 10, 057 (2021), arXiv:2105.10493 [hep-ph] .
  * Co  _et al._ [2022] R. T. Co, T. Gherghetta, and K. Harigaya, (2022), arXiv:2206.00678 [hep-ph] .
  * Anderson  _et al._ [2012a] C. Anderson _et al._ , JINST 7, P10019 (2012a), arXiv:1205.6747 [physics.ins-det] .
  * Adamson  _et al._ [2016] P. Adamson _et al._ , Nucl. Instrum. Meth. A 806, 279 (2016), arXiv:1507.06690 [physics.acc-ph] .
  * Michael  _et al._ [2008] D. G. Michael _et al._ (MINOS), Nucl. Instrum. Meth. A 596, 190 (2008), arXiv:0805.3170 [physics.ins-det] .
  * Bauer  _et al._ [2017] M. Bauer, M. Neubert, and A. Thamm, JHEP 12, 044 (2017), arXiv:1708.00443 [hep-ph] .
  * Dine  _et al._ [1981] M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. 104B, 199 (1981).
  * Zhitnitsky [1980] A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980), [Yad. Fiz.31,497(1980)].
  * Goudzovski  _et al._ [2022] E. Goudzovski _et al._ , (2022), arXiv:2201.07805 [hep-ph] .
  * [29] See Supplemental Material at [URL will be inserted by publisher] for discussion of the various axion decay modes, which includes [60].
  * Kim  _et al._ [2005] J. E. Kim, H. P. Nilles, and M. Peloso, JCAP 01, 005 (2005), arXiv:hep-ph/0409138 .
  * Kaplan and Rattazzi [2016] D. E. Kaplan and R. Rattazzi, Phys. Rev. D 93, 085007 (2016), arXiv:1511.01827 [hep-ph] .
  * Buen-Abad  _et al._ [2021] M. A. Buen-Abad, J. Fan, M. Reece, and C. Sun, JHEP 09, 101 (2021), arXiv:2104.03267 [hep-ph] .
  * Co  _et al._ [2023] R. T. Co, S. Kumar, and Z. Liu, JHEP 02, 111 (2023), arXiv:2210.02462 [hep-ph] .
  * Note [1] Other modes of axion production, such as from meson decays, are subdominant unless for example additional flavor-violating couplings are introduced. Such models are not considered in this analysis. A recent work [61] focused on such a model, but without a gluon coupling.
  * Ertas and Kahlhoefer [2020] F. Ertas and F. Kahlhoefer, JHEP 07, 050 (2020), arXiv:2004.01193 [hep-ph] .
  * Kelly  _et al._ [2021] K. J. Kelly, S. Kumar, and Z. Liu, Phys. Rev. D 103, 095002 (2021), arXiv:2011.05995 [hep-ph] .
  * Bierlich  _et al._ [2022] C. Bierlich _et al._ , (2022), arXiv:2203.11601 [hep-ph] .
  * Sjöstrand  _et al._ [2015] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, Comput. Phys. Commun. 191, 159 (2015), arXiv:1410.3012 [hep-ph] .
  * Snider and Petrillo [2017] E. Snider and G. Petrillo, J. Phys. Conf. Ser. 898, 042057 (2017).
  * Agostinelli  _et al._ [2003] S. Agostinelli _et al._ , Nucl. Instrum. Meth. A 506, 250 (2003).
  * Acciarri  _et al._ [2018] R. Acciarri _et al._ (ArgoNeuT), Phys. Rev. D 98, 052002 (2018), arXiv:1804.10294 [hep-ex] .
  * J. Spitz [2011] J. Spitz, “Ph.D. Thesis,” (2011), FERMILAB-THESIS-2011-36.
  * Acciarri  _et al._ [2021] R. Acciarri _et al._ (ArgoNeuT), Phys. Rev. Lett. 127, 121801 (2021), arXiv:2106.13684 [hep-ex] .
  * Aliaga  _et al._ [2014] L. Aliaga _et al._ (MINERvA), Nucl. Instrum. Meth. A 743, 130 (2014), arXiv:1305.5199 [physics.ins-det] .
  * Andreopoulos  _et al._ [2010] C. Andreopoulos _et al._ , Nucl. Instrum. Meth. A 614, 87 (2010), arXiv:0905.2517 [hep-ph] .
  * Anderson  _et al._ [2012b] C. Anderson _et al._ (ArgoNeuT), Phys. Rev. Lett. 108, 161802 (2012b), arXiv:1111.0103 [hep-ex] .
  * Acciarri  _et al._ [2014] R. Acciarri _et al._ (ArgoNeuT), Phys. Rev. D 89, 112003 (2014), arXiv:1404.4809 [hep-ex] .
  * Aliaga  _et al._ [2016] L. Aliaga _et al._ (MINERvA), Phys. Rev. D 94, 092005 (2016), arXiv:1607.00704 [hep-ex] .
  * Deur  _et al._ [2016] A. Deur, S. J. Brodsky, and G. F. de Teramond, Nucl. Phys. 90, 1 (2016), arXiv:1604.08082 [hep-ph] .
  * Paley  _et al._ [2014] J. M. Paley _et al._ (MIPP), Phys. Rev. D 90, 032001 (2014), arXiv:1404.5882 [hep-ex] .
  * Borasoy [1999] B. Borasoy, Phys. Rev. D 59, 054021 (1999), arXiv:hep-ph/9811411 .
  * Anderson  _et al._ [2012c] C. Anderson _et al._ (ArgoNeuT), JINST 7, P10020 (2012c), arXiv:1205.6702 [physics.ins-det] .
  * Adamson  _et al._ [2010] P. Adamson _et al._ (MINOS), Phys. Rev. D 81, 072002 (2010), arXiv:0910.2201 [hep-ex] .
  * Workman  _et al._ [2022] R. Workman _et al._ (Particle Data Group), Prog. Theor. Exp. Phys. 2022 (2022), 083C01.
  * Note [2] The impact of the uncertainties on the $c_{\ell}=1/100$ case is similar, but is not shown to aid clarity.
  * Cortina Gil  _et al._ [2021] E. Cortina Gil _et al._ (NA62), JHEP 06, 093 (2021), arXiv:2103.15389 [hep-ex] .
  * Batley  _et al._ [2017] J. R. Batley _et al._ (NA48/2), Phys. Lett. B 769, 67 (2017), arXiv:1612.04723 [hep-ex] .
  * Aaij  _et al._ [2015] R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 115, 161802 (2015), arXiv:1508.04094 [hep-ex] .
  * Aubert  _et al._ [2008] B. Aubert _et al._ (BaBar), Phys. Rev. Lett. 101, 091801 (2008), arXiv:0804.0411 [hep-ex] .
  * Aloni  _et al._ [2019] D. Aloni, Y. Soreq, and M. Williams, Phys. Rev. Lett. 123, 031803 (2019), arXiv:1811.03474 [hep-ph] .
  * Bertuzzo  _et al._ [2022] E. Bertuzzo, A. L. Foguel, G. M. Salla, and R. Z. Funchal, (2022), arXiv:2202.12317 [hep-ph] .


# Revisiting general dark matter-bound-electron interactions

Jin-Han Liang<EMAIL_ADDRESS>Yi Liao<EMAIL_ADDRESS>Xiao-
Dong Ma<EMAIL_ADDRESS>Hao-Lin Wang<EMAIL_ADDRESS>Key Laboratory
of Atomic and Subatomic Structure and Quantum Control (MOE), Guangdong Basic
Research Center of Excellence for Structure and Fundamental Interactions of
Matter, Institute of Quantum Matter, South China Normal University, Guangzhou
510006, China  Guangdong-Hong Kong Joint Laboratory of Quantum Matter,
Guangdong Provincial Key Laboratory of Nuclear Science, Southern Nuclear
Science Computing Center, South China Normal University, Guangzhou 510006,
China

###### Abstract

In this letter we revisit general dark matter (DM)-bound-electron interactions
studied previously in the influential work of [Catena et al., Phys. Rev. Res.
2, 033195 (2020)]. We derive the DM-electron response functions and find a
crucial minus sign was missed for the second atomic response function $W_{2}$
defined in that work. The minus sign has significant phenomenological
consequences when explaining experimental bounds on specific DM scenarios.
Furthermore, for the most general DM-electron nonrelativistic or relativistic
interactions for DM with spin up to one, we find there are three DM response
functions ($a_{0,1,2}$) whose corresponding atomic response functions
($\widetilde{W}_{0,1,2}$) are linear combinations of the four response
functions ($W_{1,2,3,4}$) given in that work,

$\widetilde{W}_{0}=W_{1},\,\widetilde{W}_{2}=|\bm{v}_{0}^{\perp}|^{2}W_{1}+W_{3}-2{m_{e}\,\bm{q}\cdot\bm{v}_{0}^{\perp}\over\bm{q}^{2}}W_{2},\,\widetilde{W}_{3}={(\bm{q}\cdot\bm{v}_{0}^{\perp})^{2}\over\bm{q}^{2}}W_{1}+{m_{e}^{2}\over\bm{q}^{2}}W_{4}-2{m_{e}\,\bm{q}\cdot\bm{v}_{0}^{\perp}\over\bm{q}^{2}}W_{2}.$

Due to the minus sign correction for $W_{2}$, there can be significant
cancellations between the $W_{2}$ and $W_{3,4}$ terms, so that
$\widetilde{W}_{2,3}$ are dominated by the usual response function $W_{1}$ in
some cases. Ignoring the sign could thus result in misinterpretation of the
experimental data in some DM scenarios. As an example, we show that the recent
XENON1T constraint on the fermionic DM anapole moment is weakened by a factor
of 2 or so. Many DM scenarios involving DM or electron axial-vector current
can yield $W_{2}$ and thus are potentially affected by the sign.

Introduction.—The nature of dark matter (DM) is one of the biggest puzzles in
particle astrophysics. Particle DM is a well-motivated scenario that fits very
well within the standard cosmological model ($\Lambda$CDM), and in the
meanwhile has a detectable possibility in terrestrial experiments. The mass
scale for potential DM candidates spans a very wide range of energy scales,
from the superlight axion or dark photon candidate Duffy:2009ig ;
Fabbrichesi:2020wbt ; Agrawal:2021dbo , to the keV sterile neutrino
Drewes:2016upu ; Boyarsky:2018tvu , to relatively more massive particles that
are generally termed weakly interacting massive particles (WIMPs)
Arcadi:2017kky ; Roszkowski:2017nbc . Among them, the WIMPs are the most
prominent ones in consideration of their origin in particle physics models,
the power of fitting the observational data, and the promising discovery
potential in direct detection experiments.

The search for WIMP DM begins with the nuclear recoil effect. The null results
from large xenon and argon detectors (such as PandaX PandaX-II:2017hlx , XENON
XENON:2023cxc , LUX-ZEPLIN (LZ) LZ:2022lsv , and DarkSide DarkSide:2018bpj )
in the past two decades have put stringent bounds on the DM-nucleus cross
section with DM mass above ${\cal O}(10\,\rm GeV)$ Roszkowski:2017nbc ;
Schumann:2019eaa . Due to the kinematic restriction of DM-nucleus (or DM-
nucleon) elastic scattering as well as detection threshold limitation, the DM-
nucleon interaction remains less constrained for DM mass below ${\cal
O}(10\,\rm GeV)$ from nuclear recoil signals. This limitation could be
circumvented by exploiting inelastic processes like the bremsstrahlung process
Kouvaris:2016afs  and the Migdal effect Migdal:1941 ; Ibe:2017yqa , or novel
low-threshold detectors in condensed matter systems (see Kahn:2021ttr ;
Mitridate:2022tnv  and references cited therein). In contrast to nuclear
recoil experiments searching for the DM-nucleon interaction, the DM-electron
interaction offers a more powerful alternative probe of low mass DM particles
below ${\cal O}(10\,\rm GeV)$ through the electron recoil signal LZ:2023poo ;
XENON:2019gfn ; XENON:2021qze ; XENON:2022ltv ; PandaX:2022xqx ;
DarkSide:2018ppu ; PandaX-II:2021nsg . This is due to the significantly
smaller mass of the target electron, allowing it to gain more easily the
recoil energy from a light DM particle. For instance, the single-electron
search conducted by XENON1T has the capability of exploring DM as light as
approximately 5 MeV XENON:2021qze .

Since the search for light DM through DM-electron scattering has become more
and more important than ever before, a thorough understanding of the response
functions of the electron in a target material is urgent. The dual-phase
detectors, with liquid xenon or argon as the target material, are the most
promising DM direct detection experiments to explore the parameter space of
WIMPs in a wide mass region, and can effectively distinguish between the
signals induced by the DM-electron and DM-nucleus scattering. Theoretically,
the DM-electron scattering in a liquid xenon target was first calculated in
Kopp:2009et ; Essig:2011nj ; Essig:2012yx , in which some simplified DM models
were studied, with the atomic effect of bound electrons being encoded in the
ionization form factor Essig:2011nj  or the $K$-factor Roberts:2016xfw . The
relativistic correction and many-body effect of the ionization form factor
were further considered in Roberts:2015lga ; Roberts:2016xfw ; Pandey:2018esq
, which could be sizable for a large transfer momentum. Recently, a general
study on the atomic response functions in the liquid xenon and argon targets
was carried out in Catena:2019gfa  in the framework of nonrelativistic (NR)
effective field theory (EFT), in which three new atomic response functions
were introduced in addition to the usual ionization form factor. These new
atomic response functions appear for general DM-electron interactions and are
particularly relevant for some special DM scenarios. They have been used in
subsequent theoretical Catena:2022fnk  and experimental XENON:2021qze  papers.
In this work, we rederive these response functions, and find that a crucial
minus sign was missed for the second response function $W_{2}$ defined in
Catena:2019gfa . The minus sign is physically important when explaining
experimental bounds in terms of specific DM scenarios. For the purpose of
illustration and the interest in its own right, we study the constraints on
the DM anapole operator which were considered previously in Catena:2022fnk ;
XENON:2021qze . We find that the constraints are weakened by a factor of 2 or
so upon correcting the sign of $W_{2}$. Along the line we formulate the
squared matrix element of the DM-atom scattering from the most general NR EFT
operators for DM up to spin one, and find that it can be compactly organized
into three terms (see Eq. 7), each of which is a product of a DM response
function and an atomic response function.

Formalism for DM-atom scattering.—For the DM-electron scattering in atomic
target, we label the initial state of the bound electron by the usual atomic
quantum numbers $(n,\ell,m)\equiv|1\rangle$ and the final state of the ionized
electron by $(k^{\prime},\ell^{\prime},m^{\prime})\equiv|2\rangle$ with
$k^{\prime}$ being its momentum. Denoting the initial (final) momentum of DM
by $p\,(p^{\prime})$, the transition amplitude for
$|p,1\rangle\to|p^{\prime},2\rangle$ can be written as Catena:2019gfa

$\displaystyle{\cal M}_{1\to 2}=\int\frac{{\rm
d}^{3}\bm{k}}{(2\pi)^{3}}\tilde{\psi}_{2}^{*}(\bm{k}+\bm{q}){\cal
M}(\bm{q},{\boldsymbol{v}_{\rm el}^{\perp}})\tilde{\psi}_{1}(\bm{k}),$ (1)

where $\tilde{\psi}_{1}\equiv\tilde{\psi}_{n\ell m}$ and
$\tilde{\psi}_{2}\equiv\tilde{\psi}_{k^{\prime}\ell^{\prime}m^{\prime}}$ are
the initial and final electron wavefunctions in momentum space respectively,
and are related to the wavefunctions in configuration space by the Fourier
transformation, $\tilde{\psi}(\bm{k})=\int{\rm
d}^{3}\bm{r}\,e^{-i\bm{k}\cdot\bm{r}}\psi(\bm{r})$. ${\cal
M}(\bm{q},{\boldsymbol{v}_{\rm el}^{\perp}})$ is the matrix element for free
electron-DM scattering that is only a function of the momentum transfer
$\bm{q}\equiv\bm{p}-\bm{p}^{\prime}$ and the DM-electron transverse velocity
${\boldsymbol{v}_{\rm el}^{\perp}}=\bm{v}-{\bm{q}/(2\mu_{xe})}-{\bm{k}/m_{e}}$
in the nonrelativistic framework. Here $\bm{v}$ is the incoming DM velocity in
the lab frame, $\bm{k}$ is the momentum of the initial electron of mass
$m_{e}$, and $\mu_{xe}$ the reduced DM-electron mass, respectively.

For the most general DM-electron interactions at leading order in the NR EFT
framework given in Table 1 in Appendix, which are at most linear order in
${\boldsymbol{v}_{\rm el}^{\perp}}$ or $\bm{k}$, the free amplitude is
decomposed as Catena:2019gfa , ${\cal M}={\cal
M}_{\boldsymbol{k}=0}+\boldsymbol{k}\cdot(\nabla_{\bm{k}}{\cal
M})_{\boldsymbol{k}=0}$. This results in two atomic form factors $f_{1\to
2}(\bm{q})$ and $\boldsymbol{f}_{1\to 2}(\bm{q})$ in Eq. 1:

$\displaystyle f_{1\to 2}(\bm{q})$ $\displaystyle=\int\frac{{\rm
d}^{3}\bm{k}}{(2\pi)^{3}}\tilde{\psi}_{2}^{*}(\bm{k}+\bm{q})\tilde{\psi}_{1}(\bm{k}),$
(2a) $\displaystyle\boldsymbol{f}_{1\to 2}(\bm{q})$
$\displaystyle=\int\frac{{\rm
d}^{3}\bm{k}}{(2\pi)^{3}}\tilde{\psi}_{2}^{*}(\bm{k}+\bm{q}){\bm{k}\over
m_{e}}\tilde{\psi}_{1}(\bm{k}).$ (2b)

The first one is the usual ionization form factor, while the second was first
recognized in Catena:2019gfa . We find that a minus sign was missed in the
vector form factor $\bm{f}_{1\to 2}(\bm{q})$ of Catena:2019gfa  when
transforming from momentum to configuration space,

$\displaystyle f_{1\to 2}(\bm{q})$ $\displaystyle=\int
d^{3}\bm{r}\,\psi^{*}_{2}(\bm{r})e^{i\bm{q}\cdot\bm{r}}\psi_{1}(\bm{r}),$ (3a)
$\displaystyle\boldsymbol{f}_{1\to 2}(\bm{q})$ $\displaystyle=\int
d^{3}\bm{r}\,\psi^{*}_{2}(\bm{r})e^{i\bm{q}\cdot\bm{r}}{-i\boldsymbol{\nabla}\over
m_{e}}\psi_{1}(\bm{r}).$ (3b)

This will cause a considerable impact on the matrix element squared as spelled
out below.

The decomposition of ${\cal M}$ in terms of $\bm{k}$ leads the authors in
Catena:2019gfa  to identify three additional atomic response functions
$W_{2,3,4}$ beyond the usual $W_{1}$. From their definitions reproduced in
Appendix, one sees that the minus sign in $\boldsymbol{f}_{1\to 2}(\bm{q})$ is
brought into the function $W_{2}$. To appreciate the consequence of this sign,
we decompose the free matrix element into the ${\boldsymbol{v}_{\rm
el}^{\perp}}$ independent and linearly-dependent parts as suggested by the
defining form of NR operators,

$\displaystyle{\cal M}(\bm{q},{\boldsymbol{v}_{\rm el}^{\perp}})$
$\displaystyle={\cal M}(\bm{q},0)+\bm{v}_{\rm
el}^{\perp}\cdot\nabla_{\bm{v}_{\rm el}^{\perp}}{\cal
M}(\bm{q},{\boldsymbol{v}_{\rm el}^{\perp}})$ $\displaystyle={\cal M}_{\tt
S}+\left(\bm{v}_{0}^{\perp}-{{\bm{k}}\over m_{e}}\right)\cdot{\boldsymbol{\cal
M}}_{\tt V},$ (4)

where ${\cal M}_{\tt S}\equiv{\cal M}(\bm{q},0)$, ${\boldsymbol{\cal M}}_{\tt
V}\equiv\nabla_{\bm{v}_{\rm el}^{\perp}}{\cal M}(\bm{q},{\boldsymbol{v}_{\rm
el}^{\perp}})$, and the $\bm{k}$-independent part of $\bm{v}_{\rm el}^{\perp}$
is denoted by $\bm{v}_{0}^{\perp}\equiv\bm{v}-{\bm{q}/(2\mu_{xe})}$. We find
the spin-averaged atomic matrix element squared takes the form,

$\displaystyle\overline{|{\cal M}_{1\to 2}|^{2}}$ $\displaystyle=a_{0}|f_{\tt
S}|^{2}+a_{1}|\bm{f}_{\tt V}|^{2}+{a_{2}\over x_{e}}\left|{\bm{q}\over
m_{e}}\cdot\bm{f}_{\tt V}\right|^{2}$ (5) $\displaystyle+i\,a_{3}{\bm{q}\over
m_{e}}\cdot(\bm{f}_{\tt V}\times\bm{f}_{\tt V}^{*})+2\,{\rm
Im}\left[a_{4}f_{\tt S}\bm{f}_{\tt V}^{*}\cdot{\bm{q}\over m_{e}}\right],$

where $x_{e}\equiv\bm{q}^{2}/m_{e}^{2}$ and
$y_{e}\equiv\bm{q}\cdot\bm{v}_{0}^{\perp}/m_{e}$, $f_{\tt S}(\bm{q})\equiv
f_{1\to 2}(\bm{q})$, and $\bm{f}_{\tt
V}(\bm{q})\equiv\bm{v}_{0}^{\perp}f_{1\to 2}(\bm{q})-\bm{f}_{1\to 2}(\bm{q})$.
$a_{0,1,2,3,4}$ are DM response functions depending only on the Wilson
coefficients (WCs) of NR operators and $x_{e}$. The explicit forms of
$a_{0,1,2}$ are provided in Eq. 12, Eq. 13, and Eq. 14 for the scalar,
fermion, and vector DM cases in Appendix, respectively. The other two terms
$a_{3,4}$ have no contribution to the usual DM-atom scattering after summing
over the atomic magnetic quantum numbers $(m,m^{\prime})$, and will be
presented in the accompanying long paper liao:2024lp .

With $\overline{|{\cal M}_{1\to 2}|^{2}}$ in Eq. 5, the total ionization rate
for the initial atomic orbital ($n,\ell$), ${\cal R}_{\rm ion}^{n\ell}$, can
be obtained by summing over the initial magnetic quantum number $m$ and
integrate/sum over all the allowed final electron states
($k^{\prime},\ell^{\prime},m^{\prime}$). Neglecting the effect arising from
the non-spherically symmetric nature of the DM velocity distribution in the
lab frame, the differential rate can be generally written as Essig:2015cda ;
Catena:2019gfa

$\displaystyle\frac{{\rm d}{\cal R}_{{\rm ion}}^{n\ell}}{{\rm d}\ln
E_{e}}=\frac{n_{\rm dm}}{128\pi m_{\rm dm}^{2}m_{e}^{2}}\int{\rm
d}q~{}q\int\frac{{\rm d}^{3}\bm{v}}{v}f_{\rm
dm}(\boldsymbol{v})\overline{\left|{\cal M}_{{\rm ion}}^{n\ell}\right|^{2}},$
(6)

where $n_{\rm dm}$ is the local DM number density with velocity distribution
$f_{\rm dm}(\bm{v})$ and $m_{\rm dm}$ is the DM mass. The integration over $v$
is bounded below by $v_{\rm min}$ due to the kinematics restriction. Using the
decomposition in Revisiting general dark matter-bound-electron interactions
and the squared matrix element in Eq. 5, our main result in this letter is the
recognition of the averaged ionization matrix element squared,

$\displaystyle\overline{\left|{\cal
M}_{\rm{ion}}^{n\ell}\right|^{2}}=a_{0}\widetilde{W}_{0}+a_{1}\widetilde{W}_{1}+a_{2}\widetilde{W}_{2},$
(7)

where the atomic response functions $\widetilde{W}_{0,1,2}$ are related to
$W_{1,2,3,4}$ in Catena:2019gfa  as,111Notice that both $\widetilde{W}_{i}$
and $W_{i}$ depend on the initial atomic quantum numbers $(n,\ell)$, the final
electron momentum $k^{\prime}$, and the momentum transfer $|\bm{q}|$. We do
not show them explicitly for brevity.

$\displaystyle\widetilde{W}_{0}$ $\displaystyle=W_{1},$ (8a)
$\displaystyle\widetilde{W}_{1}$
$\displaystyle=|\bm{v}_{0}^{\perp}|^{2}W_{1}-2{y_{e}\over x_{e}}W_{2}+W_{3},$
(8b) $\displaystyle\widetilde{W}_{2}$ $\displaystyle={y_{e}^{2}\over
x_{e}}W_{1}-2{y_{e}\over x_{e}}W_{2}+{1\over x_{e}}W_{4}.$ (8c)

We find the decomposition in Revisiting general dark matter-bound-electron
interactions is more amenable to obtain the square of the ionization matrix
element with $a_{0}$ and $a_{1,2}$ capturing different NR contributions, as
can be seen explicitly from Eq. 12, Eq. 13, and Eq. 14 in Appendix.
Furthermore, the recognition of the three new combinations of atomic response
functions is useful to understand which NR operators contribute dominantly to
the scattering. In Eq. 7, $\widetilde{W}_{0}$ (or $W_{1}$) is the usual atomic
$K$-factor in the literature Roberts:2016xfw ; Roberts:2019chv . Its
associated DM response function $a_{0}$ incorporates the most common DM
scenarios through the WCs $c_{1}$ and $c_{4}$ of the NR operators ${\cal
O}_{1}$ and ${\cal O}_{4}$, which are the well-known spin-independent and
spin-dependent NR interactions extensively studied in the literature.
$\widetilde{W}_{1,2}$ are mainly related to the NR operators involving
transverse velocity ${\boldsymbol{v}_{\rm el}^{\perp}}$, which can be
generated by DM-electron interactions involving electron or DM axial-vector
current such as the anapole moment of fermionic DM. Thus, a sign mistake in
$W_{2}$ will result in a wrong interpretation of experimental data for such
types of interactions.

Figure 1: Individual terms in $\widetilde{W}_{1,2}$ are shown as a function of
$|\boldsymbol{q}|$. For better layout we show $-2(y_{e}/x_{e})W_{2}$ (black
dotted curve) for $|\bm{q}|\lesssim 11.5\,\rm keV$ and $+2(y_{e}/x_{e})W_{2}$
(black solid) for $|\bm{q}|\gtrsim 11.5\,\rm keV$. The lower limit of
$|\bm{q}|$ is determined by kinematic consideration.

To assess the numerical importance of the sign, we show in Fig. 1 the
individual terms of $\widetilde{W}_{1,2}$ as a function of the momentum
transfer $|\bm{q}|$ for a xenon target in the orbital 5p. We choose a typical
light DM mass $m_{\rm dm}=500\,\rm MeV$ and fix the ionized electron momentum
to be $k^{\prime}=10\,\rm keV$. From Fig. 1, it is clear that the $W_{2}$ term
is positive for $|\bm{q}|\gtrsim 11.5\,\rm keV$, and approaches $W_{3}$ and
$W_{4}/x_{e}$ as $|\bm{q}|$ increases. With the minus sign accompanying
$W_{2}$ in $\widetilde{W}_{1,2}$, it implies there is a significant
cancellation between the $W_{2}$ and $W_{3,4}$ terms in $\widetilde{W}_{1,2}$,
leaving the $W_{1}$ term as the potentially dominant contribution across a
wide region of $|\bm{q}|$. This behavior is common to other atomic orbitals,
DM mass, and the ionized electron momentum.

Figure 2: The ionization spectrum from the NR interaction ${\cal O}_{7}$ with
a short distance origin at the xenon target from different response functions.
The DM can be either a scalar, fermion, or vector particle. Note that the
solid black and dotted green curves almost overlap.

Next we examine the influence on the differential rate due to the NR operator
${\cal O}_{7}$, which is a typical operator employed in Catena:2019gfa  to
show the importance of the three new response functions ($W_{2},W_{3},W_{4}$).
To facilitate comparison with Catena:2019gfa , we consider the xenon and argon
targets and focus on the contact (or short-distance) interaction with
$c_{10}^{\tt SD}=10^{-3}$ and $m_{\rm dm}=10\,\rm MeV$. In Fig. 2, we show for
the xenon target the individual contribution to the differential rate from
$W_{1}$, $W_{2}$, and $W_{3}$ (including all coefficients but not the minus
sign in front of $W_{2}$ in Eq. 8 for better layout) by the dotted green,
solid blue, and dotted magenta curves. From Fig. 2, it can be seen that the
contributions from $W_{2}$ and $W_{3}$ are of a similar magnitude and much
larger than that from $W_{1}$, so that they almost exactly cancel out, leading
to a total contribution (solid black curve) that is completely dominated by
$W_{1}$. For comparison, we also show the total rate reported in
Catena:2019gfa  by the dashed black curve.222There is about a global factor of
3/4 difference than the one given in Catena:2019gfa , which has been included
in our plot. It can be seen that the sign difference can cause up to an order
of magnitude difference in the differential rate. A similar plot for the argon
target is provided in Fig. 3. In that case, the cancellation between $W_{2}$
and $W_{3}$ is less complete, leading to a slightly higher total contribution
than that from $W_{1}$ alone. Note that the spectrum is the same for ${\cal
O}_{8}$ with fermion DM as indicated in Eq. 13, which will be used in the
following analysis.

Figure 3: Same as Fig. 2 except for an argon target.

Constraints on DM anapole.—In this part, we use the DM anapole operator as an
example to illustrate the importance of the sign for $W_{2}$. The anapole
interaction Lagrangian for a fermionic DM $\chi$ is

$\displaystyle\mathcal{L}_{\rm
anapole}=\frac{g}{2\Lambda^{2}}\bar{\chi}\gamma^{\mu}\gamma^{5}\chi\partial^{\nu}F_{\mu\nu},$
(9)

where $F_{\mu\nu}$ is the electromagnetic field strength tensor, $g$ is a
dimensionless coupling, and $\Lambda$ is a heavy particle scale. We consider
$\chi$ as a Majorana particle for direct comparison with the results in the
literature Catena:2019gfa ; XENON:2021qze . The anapole moment is the unique
electromagnetic property of a Majorana fermion since all others (such as
charge radius and electromagnetic dipole moments) do not exist. In the NR
limit, the DM anapole operator matches to two NR operators, $\mathcal{O}_{8}$
and $\mathcal{O}_{9}$, with the WCs,

$\displaystyle c_{8}=-c_{9}=8\,e\,m_{e}m_{\chi}{g\over\Lambda^{2}}.$ (10)

Since the sign of $W_{2}$ has a significant impact on the contribution from
${\cal O}_{8}$, as illustrated in Fig. 2, it is expected that it also plays a
crucial role in determining the constraints on the DM anapole operator. The
previous constraints were established using data from various experiments,
including XENON10, DarkSide-50, the XENON1T S2-only search (calculated in
Catena:2019gfa ), and the XENON1T single electron (SE) search conducted by the
XENON1T collaboration XENON:2021qze . All of these constraints were obtained
with a wrong sign associated with $W_{2}$. In the following, we analyze these
constraints from xenon target experiments. In addition, we calculate the
constraint using the latest S2-only data from the PandaX-4T experiment
PandaX:2022xqx .

To compare with the S2-only data observed in the experiments, we utilize the
constant-W model Essig:2017kqs ; Essig:2012yx  to convert the ${\rm d}{\cal
R}_{{\rm ion}}^{n\ell}/{\rm d}E_{e}$ spectrum into that of the number of
ionized electrons, ${\rm d}{\cal R}_{{\rm ion}}^{n\ell}/{\rm d}n_{e}$. For
each ionized electron, we assume that the number of photoelectrons (PEs)
induced by it follows a Gaussian distribution with a mean value $g_{2}$ and a
width value $\sigma_{\rm S2}$ as in Essig:2017kqs ; Essig:2012yx .
Specifically, for XENON10 (XENON1T), we take $g_{2}=27\,(33)$ and $\sigma_{\rm
S2}=6.7\,(7)$ Essig:2017kqs ; Essig:2019xkx . By convoluting with the Gaussian
distribution Essig:2012yx ; Essig:2017kqs ; XENON:2019gfn , we obtain the
spectrum of S2 signals (the number of PEs), ${\rm d}{\cal R}_{{\rm
ion}}^{n\ell}/{\rm dPE}$. For the constraints from the XENON10 and XENON1T
S2-only searches, we employ the same efficiency factors and statistical
methods as in Catena:2019gfa . With a wrong sign of $W_{2}$ we were able to
reproduce the constraints in Catena:2019gfa  to check our numerical analysis.

XENON1T-SE: Since the XENON1T collaboration did not provide all the necessary
information in their paper XENON:2021qze , it is challenging to reproduce its
constraint with high consistency. Therefore, we adopt a rescaling method to be
described below to update the constraint with the corrected sign of $W_{2}$.
The signal region in the XENON1T SE search is selected as $\rm
S2\in[14,150]\,\rm PE$, corresponding to signals with 1–5 ionized electrons
XENON:2021qze . For a specific DM mass, we find that the spectrum ${\rm
d}{\cal R}_{{\rm ion}}^{n\ell}/{\rm d}n_{e}$ with the correct-sign $W_{2}$
differs from that with the wrong-sign $W_{2}$ by nearly a global factor. For
instance, for $m_{\chi}=1\,\rm GeV$, this global factor is approximately 3 for
the anapole operator. Hence, for each DM mass, we calculate
$\sum_{n_{e}=1}^{5}{\rm d}{\cal R}_{{\rm ion}}^{n\ell}/{\rm d}n_{e}$ with both
the wrong-sign and correct-sign $W_{2}$. Then, we deliver our new constraint
by rescaling the constraint of the XENON1T collaboration.

PandaX-4T: With the exposure of 0.55 ton$\cdot$yr, 103 events (including 95.8
background events estimated) were observed in the signal region of S2
$\in[60,200]\,\rm PE$ PandaX:2022xqx . According to the Poisson distribution,
a 90 % C.L. requires the number of signal events to be less than 21.5. We
employ the total efficiency (red solid line) given in Fig. (1) of
PandaX:2022xqx . Since the parameters $g_{2}$ and $\sigma_{\rm S2}$ of the
PandaX-4T S2-only search were not given explicitly in the paper, we read from
its Fig. (3) $g_{2}=17.7$ and assume $\sigma_{\rm S2}=\sqrt{g_{2}}$.

To ensure consistent comparison with previous calculations based on a wrong-
sign $W_{2}$, we adopt the same DM parameters as Catena:2019gfa :
$\rho_{\chi}=0.4\,\rm GeV/cm^{3}$ for the local DM density, $v_{0}=220\,\rm
km/s$ for the Sun’s circular velocity, $v_{\rm esc}=544\,\rm km/s$ for the
galactic escape velocity, and $v_{\oplus}=244\,\rm km/s$ for the speed of
Earth in the galactic rest frame. Our results on XENON10, XENON1T S2, and
XENON1T SE, and the latest PandaX-4T data are shown in Fig. 4 with solid
curves. Correcting the $W_{2}$ sign, the constraints are weakened by a factor
of $\sim 2$, depending on the DM mass. With a larger DM mass, the discrepancy
tends to be larger. The constraint derived from the latest PandaX-4T data
imposes the most stringent limit for DM mass $m_{\chi}\gtrsim 20\,\rm MeV$,
reaching $g/\Lambda^{2}\lesssim 0.01\,\rm GeV^{-2}$ for $m_{\chi}\gtrsim
0.1\,\rm GeV$. In the low-mass region with $m_{\chi}\lesssim 20\,\rm MeV$, the
stronger limit is set by the previous XENON10 and XENON1T-SE experiments,
which have a low energy threshold down to one single ionized electron.

Figure 4: Constraints on DM anapole moment from DM direct detection
experiments, including XENON10, XENON1T S2-only search, XENON1T-SE search, and
PandaX-4T. The solid lines show our results with a correct-sign $W_{2}$, while
the dashed lines show the previous constraints Catena:2019gfa ; XENON:2021qze
with a wrong-sign $W_{2}$.

Conclusion.—In this letter, we studied the general DM-electron interactions
and found, for the most general leading-order nonrelativistic dark matter-
electron interactions, that those interactions naturally organize themselves
into the differential rate in the form of three DM response functions
$a_{0,1,2}$. Each of $a_{0,1,2}$ is multiplied by an atomic response function
$\widetilde{W}_{0,1,2}$ which is a linear combination of the four atomic
response functions $W_{1,2,3,4}$ defined in Catena:2019gfa . Most importantly,
we found a crucial sign mistake associated with $W_{2}$. Our results show that
this leads to important phenomenological implications. As far as general NR
operators are concerned, it will influence the interpretation of experimental
data on the NR operators entering $a_{1,2}$, i.e., ${\cal O}_{7,8,12,13,14}$
and ${\cal O}_{17,21,22,23,25,26}$ for DM up to spin one. Furthermore, at the
more fundamental level of relativistic field theory, several interesting DM
scenarios are particularly affected by the sign of $W_{2}$. Especially, the
interactions with electron or DM axial-vector or tensor currents can induce
contributions related to $W_{2}$, and thus would be affected by the sign
problem. The calculation details in this letter, together with the systematic
constraints on other nonrelativistic and relativistic interactions from the
latest experiments, will be provided in the accompanying long paper
liao:2024lp .

## Acknowledgements

This work was supported in part by the Guangdong Major Project of Basic and
Applied Basic Research No. 2020B0301030008, and by the Grants No.
NSFC-12035008, No. NSFC-12247151, No. NSFC-12305110, and No. NSFC-12347121.

## References

  * (1) L. D. Duffy and K. van Bibber, “Axions as Dark Matter Particles,” New J. Phys. 11 (2009) 105008, [arXiv:0904.3346 [hep-ph]].
  * (2) M. Fabbrichesi, E. Gabrielli, and G. Lanfranchi, “The Dark Photon,” [arXiv:2005.01515 [hep-ph]].
  * (3) P. Agrawal et al., “Feebly-interacting particles: FIPs 2020 workshop report,” Eur. Phys. J. C 81 no. 11, (2021) 1015, [arXiv:2102.12143 [hep-ph]].
  * (4) M. Drewes et al., “A White Paper on keV Sterile Neutrino Dark Matter,” JCAP 01 (2017) 025, [arXiv:1602.04816 [hep-ph]].
  * (5) A. Boyarsky, M. Drewes, T. Lasserre, S. Mertens, and O. Ruchayskiy, “Sterile neutrino Dark Matter,” Prog. Part. Nucl. Phys. 104 (2019) 1–45, [arXiv:1807.07938 [hep-ph]].
  * (6) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, “The waning of the WIMP? A review of models, searches, and constraints,” Eur. Phys. J. C 78 no. 3, (2018) 203, [arXiv:1703.07364 [hep-ph]].
  * (7) L. Roszkowski, E. M. Sessolo, and S. Trojanowski, “WIMP dark matter candidates and searches—current status and future prospects,” Rept. Prog. Phys. 81 no. 6, (2018) 066201, [arXiv:1707.06277 [hep-ph]].
  * (8) PandaX-II Collaboration, X. Cui et al., “Dark Matter Results From 54-Ton-Day Exposure of PandaX-II Experiment,” Phys. Rev. Lett. 119 no. 18, (2017) 181302, [arXiv:1708.06917 [astro-ph.CO]].
  * (9) XENON Collaboration, E. Aprile et al., “First Dark Matter Search with Nuclear Recoils from the XENONnT Experiment,” Phys. Rev. Lett. 131 no. 4, (2023) 041003, [arXiv:2303.14729 [hep-ex]].
  * (10) LZ Collaboration, J. Aalbers et al., “First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment,” Phys. Rev. Lett. 131 no. 4, (2023) 041002, [arXiv:2207.03764 [hep-ex]].
  * (11) DarkSide Collaboration, P. Agnes et al., “Low-Mass Dark Matter Search with the DarkSide-50 Experiment,” Phys. Rev. Lett. 121 no. 8, (2018) 081307, [arXiv:1802.06994 [astro-ph.HE]].
  * (12) M. Schumann, “Direct Detection of WIMP Dark Matter: Concepts and Status,” J. Phys. G 46 no. 10, (2019) 103003, [arXiv:1903.03026 [astro-ph.CO]].
  * (13) C. Kouvaris and J. Pradler, “Probing sub-GeV Dark Matter with conventional detectors,” Phys. Rev. Lett. 118 no. 3, (2017) 031803, [arXiv:1607.01789 [hep-ph]].
  * (14) A. B. Migdal, “Ionization of atoms accompanying $\alpha$-and $\beta$-decay,” J. Phys. 4, 449 (1941) (1941) .
  * (15) M. Ibe, W. Nakano, Y. Shoji, and K. Suzuki, “Migdal Effect in Dark Matter Direct Detection Experiments,” JHEP 03 (2018) 194, [arXiv:1707.07258 [hep-ph]].
  * (16) Y. Kahn and T. Lin, “Searches for light dark matter using condensed matter systems,” Rept. Prog. Phys. 85 no. 6, (2022) 066901, [arXiv:2108.03239 [hep-ph]].
  * (17) A. Mitridate, T. Trickle, Z. Zhang, and K. M. Zurek, “Snowmass white paper: Light dark matter direct detection at the interface with condensed matter physics,” Phys. Dark Univ. 40 (2023) 101221, [arXiv:2203.07492 [hep-ph]].
  * (18) LZ Collaboration, J. Aalbers et al., “Search for new physics in low-energy electron recoils from the first LZ exposure,” Phys. Rev. D 108 no. 7, (2023) 072006, [arXiv:2307.15753 [hep-ex]].
  * (19) XENON Collaboration, E. Aprile et al., “Light Dark Matter Search with Ionization Signals in XENON1T,” Phys. Rev. Lett. 123 no. 25, (2019) 251801, [arXiv:1907.11485 [hep-ex]].
  * (20) XENON Collaboration, E. Aprile et al., “Emission of single and few electrons in XENON1T and limits on light dark matter,” Phys. Rev. D 106 no. 2, (2022) 022001, [arXiv:2112.12116 [hep-ex]].
  * (21) XENON Collaboration, E. Aprile et al., “Search for New Physics in Electronic Recoil Data from XENONnT,” Phys. Rev. Lett. 129 no. 16, (2022) 161805, [arXiv:2207.11330 [hep-ex]].
  * (22) PandaX Collaboration, S. Li et al., “Search for Light Dark Matter with Ionization Signals in the PandaX-4T Experiment,” Phys. Rev. Lett. 130 no. 26, (2023) 261001, [arXiv:2212.10067 [hep-ex]].
  * (23) DarkSide Collaboration, P. Agnes et al., “Constraints on Sub-GeV Dark-Matter–Electron Scattering from the DarkSide-50 Experiment,” Phys. Rev. Lett. 121 no. 11, (2018) 111303, [arXiv:1802.06998 [astro-ph.CO]].
  * (24) PandaX-II Collaboration, C. Cheng et al., “Search for Light Dark Matter-Electron Scatterings in the PandaX-II Experiment,” Phys. Rev. Lett. 126 no. 21, (2021) 211803, [arXiv:2101.07479 [hep-ex]].
  * (25) J. Kopp, V. Niro, T. Schwetz, and J. Zupan, “DAMA/LIBRA and leptonically interacting Dark Matter,” Phys. Rev. D 80 (2009) 083502, [arXiv:0907.3159 [hep-ph]].
  * (26) R. Essig, J. Mardon, and T. Volansky, “Direct Detection of Sub-GeV Dark Matter,” Phys. Rev. D 85 (2012) 076007, [arXiv:1108.5383 [hep-ph]].
  * (27) R. Essig, A. Manalaysay, J. Mardon, P. Sorensen, and T. Volansky, “First Direct Detection Limits on sub-GeV Dark Matter from XENON10,” Phys. Rev. Lett. 109 (2012) 021301, [arXiv:1206.2644 [astro-ph.CO]].
  * (28) B. M. Roberts, V. A. Dzuba, V. V. Flambaum, M. Pospelov, and Y. V. Stadnik, “Dark matter scattering on electrons: Accurate calculations of atomic excitations and implications for the DAMA signal,” Phys. Rev. D 93 no. 11, (2016) 115037, [arXiv:1604.04559 [hep-ph]].
  * (29) B. M. Roberts, V. V. Flambaum, and G. F. Gribakin, “Ionization of atoms by slow heavy particles, including dark matter,” Phys. Rev. Lett. 116 no. 2, (2016) 023201, [arXiv:1509.09044 [physics.atom-ph]].
  * (30) M. K. Pandey, L. Singh, C.-P. Wu, J.-W. Chen, H.-C. Chi, C.-C. Hsieh, C. P. Liu, and H. T. Wong, “Constraints from a many-body method on spin-independent dark matter scattering off electrons using data from germanium and xenon detectors,” Phys. Rev. D 102 no. 12, (2020) 123025, [arXiv:1812.11759 [hep-ph]].
  * (31) R. Catena, T. Emken, N. A. Spaldin, and W. Tarantino, “Atomic responses to general dark matter-electron interactions,” Phys. Rev. Res. 2 no. 3, (2020) 033195, [arXiv:1912.08204 [hep-ph]].
  * (32) R. Catena, D. Cole, T. Emken, M. Matas, N. Spaldin, W. Tarantino, and E. Urdshals, “Dark matter-electron interactions in materials beyond the dark photon model,” JCAP 03 (2023) 052, [arXiv:2210.07305 [hep-ph]].
  * (33) J.-H. Liang, Y. Liao, X.-D. Ma, and H.-L. Wang, “A systematic investigation on dark matter-electron scattering in effective field theories,”. (In preparation).
  * (34) R. Essig, M. Fernandez-Serra, J. Mardon, A. Soto, T. Volansky, and T.-T. Yu, “Direct Detection of sub-GeV Dark Matter with Semiconductor Targets,” JHEP 05 (2016) 046, [arXiv:1509.01598 [hep-ph]].
  * (35) B. M. Roberts and V. V. Flambaum, “Electron-interacting dark matter: Implications from DAMA/LIBRA-phase2 and prospects for liquid xenon detectors and NaI detectors,” Phys. Rev. D 100 no. 6, (2019) 063017, [arXiv:1904.07127 [hep-ph]].
  * (36) R. Essig, T. Volansky, and T.-T. Yu, “New Constraints and Prospects for sub-GeV Dark Matter Scattering off Electrons in Xenon,” Phys. Rev. D 96 no. 4, (2017) 043017, [arXiv:1703.00910 [hep-ph]].
  * (37) R. Essig, J. Pradler, M. Sholapurkar, and T.-T. Yu, “Relation between the Migdal Effect and Dark Matter-Electron Scattering in Isolated Atoms and Semiconductors,” Phys. Rev. Lett. 124 no. 2, (2020) 021801, [arXiv:1908.10881 [hep-ph]].
  * (38) X.-G. He, X.-D. Ma, and G. Valencia, “FCNC B and K meson decays with light bosonic Dark Matter,” JHEP 03 (2023) 037, [arXiv:2209.05223 [hep-ph]].

##  Appendix A Appendix

Nonrelativistic operators for DM-electron scattering—The NR DM-electron
interactions can be built by using the identity and spin operators. For the
electron (and fermion or vector DM), we denote the identity and spin operators
by $\mathds{1}_{e}$ and $\bm{S}_{e}$ ($\mathds{1}_{x}$ and $\bm{S}_{x}$),
$x=\phi,\chi,X$ for the scalar, fermion, vector DM cases respectively. For the
vector DM, besides $\mathds{1}_{x}$ and $\bm{S}_{x}$, a rank-two traceless
spin tensor operator $\boldsymbol{\tilde{\cal S}}_{x}$ is needed:

$\displaystyle\tilde{{\cal
S}}_{x}^{ij}=\frac{1}{2}\left(S_{x}^{i}S_{x}^{j}+i\leftrightarrow
j\right)-\frac{2}{3}\delta^{ij}.$ (11)

Together with the momentum transfer $\boldsymbol{q}$ and DM-electron
transverse velocity ${\boldsymbol{v}_{\rm el}^{\perp}}$, one can construct the
relevant NR operators. In Table 1, we summarize the most general NR operators
for DM-electron scattering at leading order. The NR operator basis is given
for operators that are at most quadratic in $\bm{q}$ and linear in
${\boldsymbol{v}_{\rm el}^{\perp}}$ for a DM particle up to spin one. The
operators from ${\cal O}_{1}$ to ${\cal O}_{15}$ were given previously in
Refs. Catena:2019gfa ; Catena:2022fnk , and the operators from ${\cal O}_{17}$
to ${\cal O}_{26}$ are specific to vector DM-electron scattering which can be
generated by the relativistic interactions in He:2022ljo ; liao:2024lp . We
note in passing that the basis of operators for DM-quark interactions can be
obtained by changing the electron label to the flavor of quarks.

Operator | Operator | Operator
---|---|---
${\cal O}_{1}=\mathds{1}_{x}\mathds{1}_{e}$ | ${\cal O}_{10}=\mathds{1}_{x}\,\frac{i\bm{q}}{m_{e}}\cdot\bm{S}_{e}$ | ${\cal O}_{19}=\frac{\bm{q}}{m_{e}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot\frac{\bm{q}}{m_{e}}\mathds{1}_{e}$
${\cal O}_{3}=\mathds{1}_{x}\Big{(}\frac{i\bm{q}}{m_{e}}\times{\boldsymbol{v}_{\rm el}^{\perp}}\Big{)}\cdot\bm{S}_{e}$ | ${\cal O}_{11}=\bm{S}_{x}\cdot\frac{i\bm{q}}{m_{e}}\mathds{1}_{e}$ | ${\cal O}_{20}=-\frac{\bm{q}}{m_{e}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot\Big{(}\frac{\bm{q}}{m_{e}}\times\bm{S}_{e}\Big{)}$
${\cal O}_{4}=\bm{S}_{x}\cdot\bm{S}_{e}$ | ${\cal O}_{12}=-\bm{S}_{x}\cdot({\boldsymbol{v}_{\rm el}^{\perp}}\times\bm{S}_{e})$ | ${\cal O}_{21}={\boldsymbol{v}_{\rm el}^{\perp}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot\bm{S}_{e}$
${\cal O}_{5}=\bm{S}_{x}\cdot\Big{(}\frac{i\bm{q}}{m_{e}}\times{\boldsymbol{v}_{\rm el}^{\perp}}\Big{)}\mathds{1}_{e}$ | ${\cal O}_{13}=(\bm{S}_{x}\cdot{\boldsymbol{v}_{\rm el}^{\perp}})\Big{(}\frac{i\bm{q}}{m_{e}}\cdot\bm{S}_{e}\Big{)}$ | ${\cal O}_{22}=\Big{(}\frac{i\bm{q}}{m_{e}}\times{\boldsymbol{v}_{\rm el}^{\perp}}\Big{)}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot\bm{S}_{e}+{\boldsymbol{v}_{\rm el}^{\perp}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot\Big{(}\frac{i\bm{q}}{m_{e}}\times\bm{S}_{e}\Big{)}$
${\cal O}_{6}=\Big{(}\bm{S}\cdot\frac{\bm{q}}{m_{e}}\Big{)}\Big{(}\frac{\bm{q}}{m_{e}}\cdot\bm{S}_{e}\Big{)}$ | ${\cal O}_{14}=(\bm{S}_{x}\cdot\frac{i\bm{q}}{m_{e}})({\boldsymbol{v}_{\rm el}^{\perp}}\cdot\bm{S}_{e})$ | ${\cal O}_{23}=-\frac{i\bm{q}}{m_{e}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot({\boldsymbol{v}_{\rm el}^{\perp}}\times\bm{S}_{e})$
${\cal O}_{7}=\mathds{1}_{x}\,{\boldsymbol{v}_{\rm el}^{\perp}}\cdot\bm{S}_{e}$ | ${\cal O}_{15}=\bm{S}_{x}\cdot\frac{\bm{q}}{m_{e}}\Big{[}\frac{\bm{q}}{m_{e}}\cdot({\boldsymbol{v}_{\rm el}^{\perp}}\times\bm{S}_{e})\Big{]}$ | ${\cal O}_{24}={\bm{q}\over m_{e}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot\Big{(}{\bm{q}\over m_{e}}\times{\boldsymbol{v}_{\rm el}^{\perp}}\Big{)}$
${\cal O}_{8}=\bm{S}_{x}\cdot{\boldsymbol{v}_{\rm el}^{\perp}}\,\mathds{1}_{e}$ | ${\cal O}_{17}=\frac{i\bm{q}}{m_{e}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot{\boldsymbol{v}_{\rm el}^{\perp}}\,\mathds{1}_{e}$ | ${\cal O}_{25}=\Big{(}{\bm{q}\over m_{e}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot{\boldsymbol{v}_{\rm el}^{\perp}}\Big{)}\Big{(}{\bm{q}\over m_{e}}\cdot\bm{S}_{e}\Big{)}$
${\cal O}_{9}=-\bm{S}_{x}\cdot\Big{(}\frac{i\bm{q}}{m_{e}}\times\bm{S}_{e}\Big{)}$ | ${\cal O}_{18}=\frac{i\bm{q}}{m_{e}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot\bm{S}_{e}$ | ${\cal O}_{26}=\Big{(}{\bm{q}\over m_{e}}\cdot\boldsymbol{\tilde{\cal S}}_{x}\cdot{\bm{q}\over m_{e}}\Big{)}({\boldsymbol{v}_{\rm el}^{\perp}}\cdot\bm{S}_{e})$
Table 1: The leading-order nonrelativistic DM-electron interactions. Here
$\boldsymbol{q}=\boldsymbol{p}-\boldsymbol{p}^{\prime}$ with $\boldsymbol{p}$
($\boldsymbol{p}^{\prime}$) being the momentum of initial (final) DM state.
${\cal O}_{2}$ and ${\cal O}_{16}$ being quadratic in ${\boldsymbol{v}_{\rm
el}^{\perp}}$ are not shown.

Dark matter response functions—The calculation details for DM response
functions are given in the accompanying long paper liao:2024lp . For the DM up
to spin one, the DM response functions are

  * •

Scalar DM:

$\displaystyle a_{0}$ $\displaystyle=c_{1}^{2}+{1\over
4}c_{10}^{2}{\color[rgb]{0,1,1}x_{e}},$ $\displaystyle a_{1}$
$\displaystyle={1\over 4}c_{7}^{2}+{1\over
4}c_{3}^{2}{\color[rgb]{0,1,1}x_{e}},$ $\displaystyle a_{2}$
$\displaystyle=-{1\over 4}c_{3}^{2}{\color[rgb]{0,1,1}x_{e}}.$ (12)

  * •

Fermion DM:

$\displaystyle a_{0}$
$\displaystyle=c_{1}^{2}+\frac{3}{16}c_{4}^{2}+\left({1\over
8}c_{9}^{2}+{1\over 4}c_{10}^{2}+{1\over 4}c_{11}^{2}+{1\over
8}c_{4}c_{6}\right){\color[rgb]{0,1,1}x_{e}}+{1\over
16}c_{6}^{2}{\color[rgb]{0,1,1}x_{e}^{2}},$ (13a) $\displaystyle a_{1}$
$\displaystyle={1\over 4}c_{7}^{2}+{1\over 4}c_{8}^{2}+{1\over
8}c_{12}^{2}+\left({1\over 4}c_{3}^{2}+{1\over 4}c_{5}^{2}+{1\over
16}c_{13}^{2}+{1\over 16}c_{14}^{2}-{1\over
8}c_{12}c_{15}\right){\color[rgb]{0,1,1}x_{e}}+\frac{1}{16}c_{15}^{2}{\color[rgb]{0,1,1}x_{e}^{2}},$
(13b) $\displaystyle a_{2}$ $\displaystyle=-\left({1\over 4}c_{3}^{2}+{1\over
4}c_{5}^{2}-{1\over 8}c_{13}c_{14}-{1\over
8}c_{12}c_{15}\right){\color[rgb]{0,1,1}x_{e}}-\frac{1}{16}c_{15}^{2}{\color[rgb]{0,1,1}x_{e}^{2}}.$
(13c)

  * •

Vector DM:

$\displaystyle a_{0}$
$\displaystyle=c_{1}^{2}+\frac{1}{2}c_{4}^{2}+\left({1\over
3}c_{9}^{2}+{1\over 4}c_{10}^{2}+{2\over 3}c_{11}^{2}+{5\over
36}c_{18}^{2}+{1\over
3}c_{4}c_{6}\right){\color[rgb]{0,1,1}x_{e}}+\left({1\over 6}c_{6}^{2}+{2\over
9}c_{19}^{2}+{1\over 12}c_{20}^{2}\right){\color[rgb]{0,1,1}x_{e}^{2}},$ (14a)
$\displaystyle a_{1}$ $\displaystyle={1\over 4}c_{7}^{2}+{2\over
3}c_{8}^{2}+{1\over 3}c_{12}^{2}+{5\over 36}c_{21}^{2}+\left({1\over
4}c_{3}^{2}+{2\over 3}c_{5}^{2}+{1\over 6}c_{13}^{2}+{1\over
6}c_{14}^{2}+{1\over 6}c_{17}^{2}+{3\over 8}c_{22}^{2}+{7\over
72}c_{23}^{2}-{1\over 3}c_{12}c_{15}\right.$ $\displaystyle\left.+{1\over
12}c_{22}c_{23}+{1\over 12}c_{21}c_{25}-{1\over
18}c_{21}c_{26}\right){\color[rgb]{0,1,1}x_{e}}+\left({1\over
6}c_{15}^{2}+{1\over 6}c_{24}^{2}+{1\over 24}c_{25}^{2}+{1\over
18}c_{26}^{2}\right){\color[rgb]{0,1,1}x_{e}^{2}},$ (14b) $\displaystyle
a_{2}$ $\displaystyle=-\left({1\over 4}c_{3}^{2}+{2\over 3}c_{5}^{2}-{1\over
18}c_{17}^{2}+{7\over 24}c_{22}^{2}+{1\over 72}c_{23}^{2}-{1\over
3}c_{13}c_{14}-{1\over 3}c_{12}c_{15}-{1\over 36}c_{21}c_{25}-{1\over
6}c_{21}c_{26}+{1\over 4}c_{22}c_{23}\right){\color[rgb]{0,1,1}x_{e}}$
$\displaystyle-\left({1\over 6}c_{15}^{2}+{1\over 6}c_{24}^{2}-{1\over
72}c_{25}^{2}-{1\over 9}c_{25}c_{26}\right){\color[rgb]{0,1,1}x_{e}^{2}}.$
(14c)

In the above, we highlight the dependence on $x_{e}=\bm{q}^{2}/m_{e}^{2}$ in
cyan that helps to identify the NR operators with leading order contributions.
The above results are for the case of real WCs $c_{i}$ that correspond to
Hermitian operators for the elastic DM scenario. For the complex WCs, the
results are obtained by changing $c_{i}c_{j}\to{\rm Re}(c_{i}^{*}c_{j})$ in
each term. The case with complex WCs could appear in the inelastic DM
scenario, all of which will be provided in liao:2024lp .

Atomic response functions—The definition for the atomic response functions
$W_{1,2,3,4}$ is as follows Catena:2019gfa ,

$\displaystyle W_{1}$ $\displaystyle=V{4k^{\prime
3}\over{(2\pi)^{3}}}\sum_{m=-\ell}^{\ell}\sum_{\ell^{\prime}=0}^{\infty}\sum_{m^{\prime}=-\ell^{\prime}}^{\ell^{\prime}}|f_{1\to
2}(\bm{q})|^{2},$ (15a) $\displaystyle W_{2}$ $\displaystyle=V{4k^{\prime
3}\over{(2\pi)^{3}}}\sum_{m=-\ell}^{\ell}\sum_{\ell^{\prime}=0}^{\infty}\sum_{m^{\prime}=-\ell^{\prime}}^{\ell^{\prime}}{\bm{q}\over
m_{e}}\cdot f_{1\to 2}(\bm{q})\bm{f}_{1\to 2}^{*}(\bm{q}),$ (15b)
$\displaystyle W_{3}$ $\displaystyle=V{4k^{\prime
3}\over{(2\pi)^{3}}}\sum_{m=-\ell}^{\ell}\sum_{\ell^{\prime}=0}^{\infty}\sum_{m^{\prime}=-\ell^{\prime}}^{\ell^{\prime}}|\bm{f}_{1\to
2}(\bm{q})|^{2},$ (15c) $\displaystyle W_{4}$ $\displaystyle=V{4k^{\prime
3}\over{(2\pi)^{3}}}\sum_{m=-\ell}^{\ell}\sum_{\ell^{\prime}=0}^{\infty}\sum_{m^{\prime}=-\ell^{\prime}}^{\ell^{\prime}}|{\bm{q}\over{m_{e}}}\cdot\bm{f}_{1\to
2}(\bm{q})|^{2}.$ (15d)

Here the summations have been performed over the magnetic quantum number ($m$)
for the initial electron state and the angular and magnetic quantum numbers
($\ell^{\prime}$ and $m^{\prime}$) for the final electron state, and $V$ is
the normalization volume. Notice the definition of $W_{1,2,3,4}$ are
independent of specific DM properties. Although the three new combinations of
atomic response functions ($\widetilde{W}_{0,1,2}$) are related to the DM mass
through $y_{e}$, we emphasize that they are useful for understanding which NR
operators give dominant contributions, as discussed in the paper.