aainstitutetext: C.N. Yang Institute for Theoretical Physics, Stony Brook
University, NY 11794, USAbbinstitutetext: School of Natural Sciences,
Institute for Advanced Study, Princeton, NJ 08540, USAccinstitutetext:
Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2,
I-00185 Rome, Italyddinstitutetext: INFN Sezione di Roma, Piazzale Aldo Moro
2, I-00185 Rome, Italyeeinstitutetext: Department of Physics, University of
California, San Diego, CA 92093, USA

# The Migdal Effect in Semiconductors for Dark Matter with Masses below
$\sim\,$100 MeV

Kim V. Berghaus b,c,d Angelo Esposito a Rouven Essig e and Mukul Sholapurkar
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Dark matter scattering off a nucleus has a small probability of inducing an
observable ionization through the inelastic excitation of an electron, called
the Migdal effect. We use an effective field theory to extend the computation
of the Migdal effect in semiconductors to regions of small momentum transfer
to the nucleus, where the final state of the nucleus is no longer well
described by a plane wave. Our analytical result can be fully quantified by
the measurable dynamic structure factor of the semiconductor, which accounts
for the vibrational degrees of freedom (phonons) in a crystal. We show that,
due to the sum rules obeyed by the structure factor, the inclusive Migdal rate
and the shape of the electron recoil spectrum is well captured by
approximating the nuclei in the crystal as free ions; however, the exclusive
differential rate with respect to energy depositions to the crystal depends on
the phonon dynamics encoded in the dynamic structure function of the specific
material. Our results now allow the Migdal effect in semiconductors to be
evaluated even for the lightest dark matter candidates ($m_{\chi}\gtrsim 1$
MeV) that can kinematically excite electrons.

##  1 Introduction

Dark matter (DM) direct detection experiments seek evidence for non-
gravitational interactions between DM particles and ordinary matter. In
particular, the search for DM interactions with nuclei has seen tremendous
progress in the past few decades in probing DM masses above the proton, with
one- to multi-ton scale detectors now leading the search Aprile  _et al._
(2018); Aalbers  _et al._ (2022); Meng  _et al._ (2021). In this mass range,
typical searches look for DM particles scattering elastically off the nuclei
in the detector, as the resulting nuclear recoil energies are sufficiently
large to be above detector thresholds. However, for DM with masses below the
GeV-scale, the small recoil energies from elastic DM–nucleus scatterings is
challenging to detect, leading to much weaker constraints when compared to
those for heavier DM Alkhatib  _et al._ (2021); Abdelhameed  _et al._ (2019).
Instead, stronger constraints on DM–nucleus interactions can be obtained by
searching for inelastic DM processes, in which the DM can transfer an
$\mathcal{O}(1)$ fraction of its kinetic energy to the target Essig  _et al._
(2012); Kouvaris and Pradler (2017); Ibe  _et al._ (2018).

For DM–nucleus scattering, one can make use of the Migdal effect, in which the
small nuclear recoil is accompanied by a prompt electron, which can absorb
most of the DM kinetic energy Migdal (1939a, b); Migdal and Ballentine (1978);
Vergados and Ejiri (2005); Moustakidis  _et al._ (2005); Bernabei  _et al._
(2007); Ibe  _et al._ (2018). The first experimental results are promising
Akerib  _et al._ (2019); Liu  _et al._ (2019); Essig  _et al._ (2020);
Armengaud  _et al._ (2019); Aprile  _et al._ (2019a); Barak  _et al._ (2020);
Knapen  _et al._ (2022); Adhikari  _et al._ (2022); Al-Bakry  _et al._ (2022);
Armengaud  _et al._ (2022); Agnes  _et al._ (2022) and there is a clear path
for improvement Essig  _et al._ (2022). In addition, proposals exist to
discover and calibrate the Migdal effect in the laboratory with neutrons
Nakamura  _et al._ (2021); Liao  _et al._ (2021); Bell  _et al._ (2022);
Araújo  _et al._ (2022); Cox  _et al._ (2022); Adams  _et al._ (2022).
However, while the theoretical description of the Migdal effect has improved
significantly in the past few years Dolan  _et al._ (2018); Bell  _et al._
(2020); Baxter  _et al._ (2020); Essig  _et al._ (2020); Liang  _et al._
(2020); Liu  _et al._ (2020); Kahn  _et al._ (2021); Flambaum  _et al._
(2020); Bell  _et al._ (2021); Acevedo  _et al._ (2022); Wang  _et al._
(2022); Blanco  _et al._ (2022); Cox  _et al._ (2022); Liang  _et al._ (2021);
Knapen  _et al._ (2021); Liang  _et al._ (2022), a general description of the
Migdal effect for low-mass DM ($m_{\chi}\lesssim 50$ MeV) in solid-state
crystal (semiconductor) targets is still lacking. In this region, the DM
transfers very little energy and momentum to the nucleus and is unable to
dislodge it from its lattice site. A proper description of this region, which
includes the possibility of phonon production, is required. In this paper, we
will provide such a description based on an effective field theory (EFT).

The outline of the paper is as follows. In Sec. 2, we review the kinematics of
the Migdal effect in semiconductors, provide an overview of the various
theoretical calculations available in the literature, and contrast these with
our effective field theory approach. In Sec. 3, we derive the effective
Hamiltonian for the Migdal effect in semiconductors, and use it to derive an
expression for the DM scattering rate in terms of the dynamic structure
factor, $S(\bm{q},E)$, a measurable quantity that depends only on the
material’s properties. In Sec. 4, we discuss how to obtain the dynamic
structure factor from data and theoretical estimates. In Sec. 5, we show for
low DM masses the total DM–nucleus Migdal scattering rates, as well as the
differential rates, in electron recoil energy and in phonon energies for
silicon and germanium targets. In Sec. 6, we summarize our main findings and
conclude with future directions. In the Appendices, we present additional
derivations, as well as a discussion on the prospects of using neutron
scattering to obtain comprehensive data quantifying the dynamic structure
factor.

##  2 The Migdal effect in semiconductors in the low momentum regime

Semiconductors are excellent DM direct-detection targets. Their bandgaps,
denoted here by $\omega_{\rm g}$, are $\mathcal{O}(\text{eV})$, so that DM
with masses as low as $\mathcal{O}$(MeV) (which have kinetic energies
$\frac{1}{2}m_{\chi}v_{\chi}^{2}\sim\mathcal{O}(\text{eV})$, assuming
$v_{\chi}\sim 10^{-3}$) can deposit sufficient energy to excite an electron
across the bandgap and create an electron–hole pair. The excitation energies
are smaller than the ionization threshold of free atoms, which are
$\mathcal{O}(10\text{ eV})$. However, extending the description of the Migdal
effect for a free atom to semiconductor targets requires a careful treatment
of the fact that the nuclei are bound Knapen  _et al._ (2021); Liang  _et al._
(2022), the electrons are not localized Liang  _et al._ (2020), and the boost
operator inducing the ionization in the free-atom treatment cannot simply be
applied in a crystal with a preferred rest frame Essig  _et al._ (2020).
Recent work Knapen  _et al._ (2021); Liang  _et al._ (2021) has shown that the
Migdal effect in semiconductors can be described by a process of DM–nucleus
scattering in association with a nucleus–electron Coulomb interaction, the
latter being responsible for the ionization of one of the electrons in the
material. In particular, it has been shown that the well-studied energy loss
function (ELF) can be used to encode the screening from valence electrons.

The recent work Knapen  _et al._ (2021); Liang  _et al._ (2021) resolves two
of the three subtleties introduced by semiconductors, but does not fully
address the bound nature of the nuclei. In the treatment presented in Knapen
_et al._ (2021); Liang  _et al._ (2021), the final state of the nucleus is
treated as a plane wave, hence ignoring the effects coming from the rest of
the crystal lattice and any collective effects (phonons) that can be generated
along with the Migdal electron during the DM scattering event. The
approximation of treating the final state of the nucleus as a plane wave is
usually dubbed impulse approximation (ia), although within this approximation,
the treatment of the initial state of the nucleus in the literature varies. In
what Knapen  _et al._ (2021) call impulse approximation, the initial state of
the nucleus is treated as the ground state of an harmonic crystal. In Liang
_et al._ (2021), instead, impulse approximation refers to treating both the
initial and final states as plane waves. In this paper, our definition of the
impulse approximation is consistent with that in Knapen  _et al._ (2021),
whereas we label the treatment in Liang  _et al._ (2021) as the free-ion
approximation. To summarize, in this paper, we assume:

impulse approx.: | final (initial) state of nucleus is plane wave (harmonic crystal)
---|---
free-ion approx.: | final (initial) state of nucleus is plane wave (plane wave)

Both approximations are only good approximations when DM–nucleus collision
happens over time scales, $t_{\rm coll}$, much shorter than the typical time
scale characterizing the dynamics of the lattice, i.e., for

$\displaystyle t_{\rm coll}\sim 1/E_{r}\ll t_{\rm ph}\sim 1/\langle E_{\rm
ph}\rangle\,\qquad\textrm{(impulse approximation)\,,}$ (1)

with the free-ion approximation breaking down faster than the impulse
approximation (see Fig. 1). Here $\langle E_{\rm ph}\rangle$ is the average
phonon energy and $E_{r}$ is the energy of the recoiling nucleus,

$\displaystyle E_{r}={}$ $\displaystyle\frac{q^{2}}{2m_{\rm N}}$ (2a)
$\displaystyle\sim{}$ $\displaystyle\frac{m^{2}_{\chi}v^{2}_{\chi}}{2m_{\rm
N}}\simeq 0.04\leavevmode\nobreak\ \textrm{eV}\leavevmode\nobreak\
\left(\frac{m_{\chi}}{50\leavevmode\nobreak\
\textrm{MeV}}\right)^{2}\left(\frac{28\leavevmode\nobreak\
\textrm{GeV}}{m_{\rm N}}\right)\,,$ (2b)

where $q$ is the momentum lost by the DM, $m_{\rm N}$ is the mass of the
nucleus (taken to be the approximate mass of a silicon nucleus in the second
line), and we assume a DM velocity of $v_{\chi}\sim 10^{-3}$. Combining Eqs.
(1) and (2a), the impulse approximation is valid for
$q\gg\sqrt{2m_{N}\langle{E_{\textrm{ph}}}\rangle}$. For silicon and germanium,
the average phonon energies are $\langle
E^{\text{Si}}_{\text{ph}}\rangle\simeq 0.04\text{ eV}$ and $\langle
E^{\text{Ge}}_{\text{ph}}\rangle\simeq 0.025\text{ eV}$, respectively.
Combining Eqs. (1) and (2b), we find

$\displaystyle m_{\chi}\gg 50\leavevmode\nobreak\
\textrm{MeV}\sqrt{\frac{m_{\rm N}}{28\leavevmode\nobreak\
\textrm{GeV}}}\sqrt{\frac{\langle
E_{\text{ph}}\rangle}{0.04\,\text{eV}}}\qquad\textrm{(impulse
approximation)\,.}$ (3)

We see that the condition in Eq. (3) is no longer satisfied for DM candidates
with masses $m_{\chi}\lesssim 50$ MeV. Given that DM candidates with masses as
small as $\mathcal{O}$(MeV) are still kinematically able to excite an electron
in a semiconductor via the Migdal effect, it is imperative to have a robust
theoretical understanding of the Migdal effect also for $1\leavevmode\nobreak\
\textrm{MeV}\lesssim m_{\chi}\lesssim 50\leavevmode\nobreak\ \textrm{MeV}$,
where the impulse approximation fails.

In Fig. 1, we show a schematic representation of the break-down of the impulse
approximation as a function of the momentum transfer and of the dark matter
mass, and overlay it with the regime of validity of the EFT approach we employ
in this work to extend a description of the Migdal effect in semiconductors to
the low-mass regime. Fig. 1 also shows the regime of validity of the
incoherent as well as the harmonic approximations, introduced in Sec. 4, and
which are relevant for the numerical results presented in Sec. 5.

Figure 1: Schematic representation of the regimes of validity of the different
approximations discussed in this work. Here $m_{\rm N}$ is the nucleus mass,
$a$ the lattice spacing, $\langle E_{\rm ph}\rangle$ the typical average
phonon energy, and $\omega_{\rm g}$ the semiconductor gap. Both the impulse
and the free-ion approximation are valid when the momentum transfer is much
larger than the typical momentum scale set by the crystal, $q\gg\sqrt{2m_{\rm
N}\langle E_{\rm ph}\rangle}$. However, the free-ion approximation is
parametrically slower than the impulse approximation at converging to the
corresponding approximate value. Corrections to the free-ion approximations
are, in fact, proportional to $\sqrt{\langle E_{\rm ph}\rangle/E_{r}}$, while
corrections to the impulse approximation are proportional to $\langle E_{\rm
ph}\rangle/E_{r}$ (see Sec. 4.1.2 for details). The harmonic approximation is,
instead, valid in the opposite regime, $q\ll\sqrt{2m_{\rm N}\langle E_{\rm
ph}\rangle}$, where corrections to the atomic potential that are higher than
quadratic in the positions are neglected, and the displacement of the nucleus
by the DM scatter is small compared to the inter-atomic spacing. The
incoherent approximation applies to processes where the momentum exchanged is,
instead, much larger than the inverse lattice spacing, $q\gg 1/a$, when the
scattering particle is able to discriminate each single nucleus and hence does
not scatter coherently off multiple lattice sites. Finally, the EFT presented
in this work is valid when the separation between the lattice energy levels
excited by the DM, $\Delta E_{\lambda}$, is much smaller than the energy of
the Migdal electron, $\omega$—i.e., for momentum transfers $q\ll\sqrt{2m_{\rm
N}\omega_{\rm g}}$. Note that close to the boundaries of the indicated regions
of validity, the approximations will receive $\mathcal{O}(1)$ corrections. The
numerical values are estimated for the case of silicon, for which $m_{\rm
N}\simeq 28$ GeV, $\omega_{\rm g}\simeq 1.2$ eV, and $\langle E_{\rm
ph}\rangle\simeq 0.04$ eV. The typical DM velocity is taken to be
$v_{\chi}\simeq 10^{-3}$.

In this work, we reformulate the low-momentum Migdal effect in semiconductors
in terms of an effective field theory, which is insensitive to the complicated
short-distance details of the crystal lattice. We are able to do so by
exploiting the separation of scales between the Migdal ionization energy,
$\omega\gtrsim\mathcal{O}(1\text{ eV})$, and the typical energy of the lattice
excitations, $E_{\rm ph}\lesssim\mathcal{O}(0.1\text{ eV)}$. By integrating
out the short-distance (high-momentum) vibrational modes, we simplify the
description of the Migdal effect in semiconductors from a second-order effect
in old fashioned perturbation theory (a DM–nucleus interaction plus a
nucleus–electron Coulomb interaction) to a first order effect with an
effective Hamiltonian, the latter being independent of the complicated
interactions between the ions of the crystal.

The low-energy effective theory allows for a description of the final rate
that only depends on the dynamic structure factor of the target material,
$S(\bm{q},E)$—a measurable quantity that is fully characterized depends by the
momentum transferred to the lattice, $\bm{q}$, and the energy deposited to the
lattice, $E$. We make no assumptions about the initial, intermediate, or final
state of the nucleus in this EFT approach. We accurately quantify the double
differential Migdal emission rate $dR/(d\omega dE)$ in a semiconductor, where
$\omega$ is the energy of the ionized electron, an observable relevant for
next generation detectors with sensitivity to single- and multi-phonon
excitations. Interestingly, we use sum rules to show that the single
differential rate, $dR/d\omega$, is independent of the details of the lattice
dynamics, and equivalent to the same result obtained in the free-ion
approximation.

This extends the range of validity of the description of the Migdal effect in
semiconductors to all kinematically accessible DM masses, $m_{\chi}\gtrsim
1\text{ MeV}$. Additionally, the use of an effective field theory allows for
the systematic inclusion of higher order terms in the small frequency
expansion.

While our work was in progress, a proposal for a possible description on how
to take into effect the bound nucleus nature for small DM masses appeared in
Liang  _et al._ (2022). Inspired by standard field theory, the authors derive
Feynman rules for multi-phonon and electron–phonon interactions, and show
that, when the number of phonons in the final state is large, their results
are well reproduced by the impulse approximation. This approach, however, is
only applicable within the harmonic approximation for the lattice dynamics,
which treats the inter-atomic potential as purely quadratic in the positions,
and where the theory can be solved exactly. Our derivation, instead, holds
generally for any crystal. As we will show, when evaluated within the impulse
approximation, our result reduces to that of Knapen  _et al._ (2021), while
within the harmonic approximation, our analytical result agrees with Liang
_et al._ (2022).111The numerical results shown in v1 and v2 of Liang  _et al._
(2022) differ from our results in the low mass regime $(1\text{ MeV}-50\text{
MeV})$. The authors omitted a term in their numerical calculation, which they
included in v3 after our correspondence. We discuss the relevant term in more
detail in Section 5.

##  3 An effective theory for the Migdal effect

###  3.1 The effective Hamiltonian

As highlighted in Knapen  _et al._ (2021), the Migdal emission of an electron
in a solid state material formally involves two interactions, one between the
DM particle and the nuclei in the crystal, and another between the nuclei and
the electrons in the system, which induces their ionization. These two
interactions are mediated by a lattice vibrational degree of freedom. For a
realistic material, the intermediate lattice mode inducing the ionization is a
complicated eigenstate of the full lattice Hamiltonian, whose complete
description is prohibitive. For this reason, it is often convenient to
simplify the treatment by making some assumptions, most notably the impulse
Knapen  _et al._ (2021) and/or the harmonic Liang  _et al._ (2022)
approximations. While they do make the problem tractable, both these
approximations break down at given energy scales (Fig. 1) and have corrections
that are difficult to quantify.

Luckily, in the context of interest to us
($m_{\chi}\lesssim\mathcal{O}(100\text{ MeV})$), the problem is characterized
by a separation of scales: the energy of the Migdal electron is at least a few
eV,222The approximate thresholds to excite two electron–hole pairs in silicon
and germanium are $5\text{ eV}$ and $3.7\text{ eV}$, respectively. which is
much larger than the typical energy of the lattice degrees of freedom excited
in the process, which is of the order of tens to hundreds meV. As we will
show, in this regime, the lattice degree of freedom mediating the ionization
is highly off-shell. Following the standard EFT approach, one can integrate it
out, and the final effective interaction is independent of the complicated
microscopic dynamics of the crystal lattice. Operationally, this is achieved
by expanding the matrix element for large electron energy, $\omega$, reducing
the second order interaction to a first order interaction. As we show below,
the latter is such that it does not depend on the inter-atomic potential of
the crystal lattice. The schematic EFT procedure is represented in Fig. 2.

$\bm{p}_{e}$$\lambda_{i}$$\bm{p}_{i}$$\bm{p}_{e}+\bm{k}_{e}$$\lambda_{f}$$\bm{p}_{f}$$\lambda$$H_{\chi\rm
L}$$H_{e\rm
L}$$+$$\bm{p}_{e}$$\lambda_{i}$$\bm{p}_{i}$$\bm{p}_{e}+\bm{k}_{e}$$\lambda_{f}$$\bm{p}_{f}$$\lambda$$H_{e\rm
L}$$H_{\chi\rm
L}$$\simeq$$\bm{p}_{e}$$\lambda_{i}$$\bm{p}_{i}$$\bm{p}_{e}+\bm{k}_{e}$$\lambda_{f}$$\bm{p}_{f}$$H_{\rm
eff}$

Figure 2: Schematic representation of the EFT procedure applied to the result
of old-fashioned perturbation theory. Each line corresponds to a state in the
Hilbert spaces of the dark matter (dashed), electron (solid), and crystal
lattice (wavy). The intermediate lattice mode has high frequency and, when
integrated out, it leads to an effective Hamiltonian that is local in time and
independent on the complicated dynamics of the lattice.

Let us start by treating the problem in full generality. We will mostly adopt
the notation of Knapen  _et al._ (2021). In old-fashioned perturbation theory,
the Migdal emission is a second order process, whose matrix element is given
by,

$\displaystyle\begin{split}\mathcal{M}_{fi}={}&\sum_{\lambda}\bigg{[}\frac{\langle\lambda_{f},\bm{p}_{e}+\bm{k}_{e}|H_{e\rm
L}|\lambda,\bm{p}_{e}\rangle\langle\bm{p}_{f},\lambda|H_{\chi\rm
L}|\bm{p}_{i},\lambda_{i}\rangle}{\omega+E_{\lambda_{f}}-E_{\lambda}}\\\
&\quad\,+\frac{\langle\bm{p}_{f},\lambda_{f}|H_{\chi\rm
L}|\bm{p}_{i},\lambda\rangle\langle\lambda,\bm{p}_{e}+\bm{k}_{e}|H_{e\rm
L}|\lambda_{i},\bm{p}_{e}\rangle}{E_{\lambda_{i}}-E_{\lambda}-\omega}\bigg{]}\,.\end{split}$
(4)

Here, $\bm{p}_{i}$ and $\bm{p}_{f}$ are the initial and final DM momenta,
$\bm{p}_{e}$ and $\bm{p}_{e}+\bm{k}_{e}$ the initial and final electron
(crystal) momenta, and $\lambda_{i}$, $\lambda_{f}$, and $\lambda$ are the
initial, final, and intermediate lattice states, respectively. Moreover,
$\omega$ is the energy gained by the electron. The Hamiltonian $H_{e\rm{L}}$
represents the coupling between an electron and the lattice, while
$H_{\chi\rm{L}}$ analogously represents the coupling of the DM to the lattice.

For the range of DM masses of interest to us, the electron energy is always
substantially larger than the typical spacing between the lattice eigenvalues
entering in Eq. (4),
$\omega\gg|E_{\lambda_{f}}-E_{\lambda}|\,,|E_{\lambda_{i}}-E_{\lambda}|$. This
is because, for a given external kinematics, the matrix elements in the
numerator will suppress the contribution from intermediate states that are too
far away from $E_{\lambda_{i}}$ and $E_{\lambda_{f}}$.333This approximation
clearly breaks down for sufficiently heavy DM particles, which are able to
excite lattice modes whose energy is much higher than the electron’s recoil
energy—even kick a nucleus out of the crystal, in the extreme case. If we
estimate the energy given to the crystal as the recoil energy of a free
nucleus, $q^{2}/2m_{\rm N}$, the EFT breaks down when $q\sim\sqrt{2m_{\rm
N}\omega_{\rm g}}$, with $\omega_{\rm g}$ the semiconductor bandgap. In the
case of silicon and germanium, one has $\omega_{\rm g}\simeq 1.2$ eV and
$\omega_{\rm g}\simeq 0.7$ eV, respectively. This corresponds to
$m_{\chi}\gtrsim 250$ MeV and $m_{\chi}\gtrsim 320$ MeV. We note that, in
principle, our EFT is always valid for sufficiently energetic electrons,
$\omega\gg q^{2}/(2m_{\rm N})$. Here we take the most conservative value,
$\omega\sim\omega_{\rm g}$. We can then expand the matrix element, which reads

$\displaystyle\begin{split}\mathcal{M}_{fi}={}&\sum_{\lambda}\Bigg{[}\bigg{(}\frac{1}{\omega}-\frac{E_{\lambda_{f}}-E_{\lambda}}{\omega^{2}}\bigg{)}\langle\lambda_{f},\bm{p}_{e}+\bm{k}_{e}|H_{e\rm
L}|\lambda,\bm{p}_{e}\rangle\langle\bm{p}_{f},\lambda|H_{\chi\rm
L}|\bm{p}_{i},\lambda_{i}\rangle\\\
&-\bigg{(}\frac{1}{\omega}-\frac{E_{\lambda}-E_{\lambda_{i}}}{\omega^{2}}\bigg{)}\langle\bm{p}_{f},\lambda_{f}|H_{\chi\rm
L}|\bm{p}_{i},\lambda\rangle\langle\lambda,\bm{p}_{e}+\bm{k}_{e}|H_{e\rm
L}|\lambda_{i},\bm{p}_{e}\rangle\Bigg{]}+\mathcal{O}\big{(}1/\omega^{3}\big{)}\,,\\\
={}&\frac{1}{\omega^{2}}\bigg{[}\langle\bm{p}_{f},\lambda_{f},\bm{p}_{e}+\bm{k}_{e}|\big{(}H_{e\rm
L}H_{\rm L}-H_{\rm L}H_{e\rm L}\big{)}H_{\chi\rm L}-H_{\chi\rm
L}\big{(}H_{e\rm L}H_{\rm L}-H_{\rm L}H_{e\rm
L}\big{)}|\bm{p}_{i},\lambda_{i},\bm{p}_{e}\rangle\bigg{]}\\\
&+\mathcal{O}\big{(}1/\omega^{3}\big{)}\,,\end{split}$ (5)

where when going to the second equation we used the fact that the
$|\lambda\rangle$’s are eigenstates of the unperturbed lattice Hamiltonian,
$H_{\rm L}$, as well as the associated completeness relation,
$\sum_{\lambda}|\lambda\rangle\langle\lambda|=\mathbbm{1}$. One deduces the
effective interaction Hamiltonian,

$\displaystyle H_{\rm eff}=\frac{1}{\omega^{2}}\big{[}H_{\chi\rm L},[H_{\rm
L},H_{e\rm L}]\big{]}+\mathcal{O}\big{(}1/\omega^{3}\big{)}\,,$ (6)

such that
$\mathcal{M}_{fi}=\langle\lambda_{f},\bm{p}_{e}+\bm{k}_{e},\bm{p}_{f}|H_{\rm
eff}|\lambda_{i},\bm{p}_{e},\bm{p}_{i}\rangle+\mathcal{O}\big{(}1/\omega^{3}\big{)}$.
We can now specialize our result to the case at hand. The DM–lattice and
electron–lattice interactions are described by Knapen  _et al._ (2021),

$\displaystyle H_{\chi\rm L}={}$ $\displaystyle-\sum_{I}\frac{g_{\chi}g_{\rm
N}}{4\pi}\frac{e^{-m_{\phi}|\bm{x}_{\chi}-\bm{x}_{I}|}}{|\bm{x}_{\chi}-\bm{x}_{I}|}\equiv\sum_{I}H_{\chi\rm
L}^{(I)}\,,$ (7a) $\displaystyle H_{e\rm L}={}$
$\displaystyle-\frac{4\pi\alpha}{V}\sum_{I}\sum_{\bm{K},{\bm{K}}^{\prime}}\sum_{\bm{k}}\frac{\epsilon_{\bm{K}\bm{K}^{\prime}}^{-1}(\bm{k},\omega)Z(|\bm{k}+\bm{K}^{\prime}|)}{|\bm{k}+\bm{K}||\bm{k}+\bm{K}^{\prime}|}e^{i(\bm{k}+\bm{K})\cdot\bm{x}_{e}}e^{-i(\bm{k}+\bm{K}^{\prime})\cdot\bm{x}_{I}}\equiv\sum_{I}H_{e\rm
L}^{(I)}\,,$ (7b)

where $\alpha$ is the fine structure constant, $V$ the volume of the material,
and $Z(k)$ the effective atomic number accounting for the tightly bound core
electrons. Its momentum dependence encodes the fact that, at high momentum,
one probes deeper into the core electrons, effectively experiencing a larger
positive charge. Moreover,
$\epsilon_{\bm{K}\bm{K}^{\prime}}^{-1}(\bm{k},\omega)$ is the longitudinal
dielectric matrix (symmetrized over $\bm{K}$ and $\bm{K}^{\prime}$),
$\bm{K}^{(\prime)}$ is a reciprocal lattice vector, while $\bm{k}$ is confined
to the first Brillouin zone. Finally, $\bm{x}_{\chi}$, $\bm{x}_{e}$, and
$\bm{x}_{I}$ are respectively the positions of the DM, the Migdal electron,
and the $I$-th nucleus. For pedagogical reasons, we report the derivation of
the electron–lattice Hamiltonian in Appendix A. For the DM–lattice
interaction, we are assuming for concreteness a DM–nucleus interaction due to
a Yukawa coupling with a scalar (vector) mediator of mass $m_{\phi}$
($m_{A}$), hence allowing for both short and long range forces. Specifically,
we consider an interaction Lagrangian given by $\mathcal{L^{\phi}}_{\rm
int}=g_{\chi}\phi\bar{\chi}\chi+g_{\rm N}\phi\bar{N}N$ ($\mathcal{L^{A}}_{\rm
int}=g_{\chi}A_{\mu}\bar{\chi}\gamma^{\mu}\chi+g_{\rm
N}A_{\mu}\bar{N}\gamma^{\mu}N$,), where $N$ denotes the nucleus.

The lattice Hamiltonian can generically be written as

$\displaystyle H_{\rm L}=-\sum_{I}\frac{\nabla_{I}^{2}}{2m_{\rm
N}}+U\big{(}\\{\bm{x}_{I}\\}\big{)}\,,$ (8)

where $\bm{\nabla}_{I}$ is the gradient with respect to the $\bm{x}_{I}$
position and $m_{\rm N}$ the nucleus mass. From now on, we assume identical
atoms to keep the notation simple though our formalism can easily be
generalized to multi-atomic crystals such as GaAs Derenzo  _et al._ (2017,
2018). Importantly, for a realistic material the inter-atomic potential, $U$,
is in general a complicated anharmonic function. Nonetheless, the effective
Hamiltonian (6) only depends on commutators. Since $H_{e\rm L}$ and
$H_{\chi\rm L}$ are $c\,$-numbers, they commute with the potential, and the
effective interaction is therefore independent of the potential,

$\displaystyle H_{\rm eff}=\frac{1}{m_{\rm
N}\omega^{2}}\sum_{I}\bm{\nabla}_{I}H_{\chi\rm
L}^{(I)}\cdot\bm{\nabla}_{I}H_{e\rm
L}^{(I)}+\mathcal{O}\big{(}1/\omega^{3}\big{)}\,.$ (9)

###  3.2 The rate for Migdal emission

It is now straightforward to show that, given the above Hamiltonian, the
corresponding matrix element is,

$\displaystyle\begin{split}\mathcal{M}_{fi}={}&\\!-\frac{4\pi\alpha
g_{\chi}g_{\rm N}}{m_{\rm
N}\omega^{2}V^{2}}\sum_{I}\\!\sum_{\bm{K},\bm{K}^{\prime}}\\!\sum_{\bm{k}}\frac{\bm{q}\\!\cdot\\!(\bm{k}+\bm{K}^{\prime})}{q^{2}+m_{\phi}^{2}}\frac{\epsilon^{-1}_{\bm{K}\bm{K}^{\prime}}(\bm{k},\omega)}{|\bm{k}+\bm{K}||\bm{k}+\bm{K}^{\prime}|}Z(|\bm{k}+\bm{K}^{\prime}|)\\\
&\times\langle\lambda_{f}|e^{i(\bm{k}+\bm{K}^{\prime}-\bm{q})\cdot\bm{x}_{I}}|\lambda_{i}\rangle\langle\bm{p}_{e}+\bm{k}_{e}|e^{i(\bm{k}+\bm{K})\cdot\bm{x}_{e}}|\bm{p}_{e}\rangle\\\
={}&\\!-\frac{4\pi\alpha g_{\chi}g_{\rm N}}{m_{\rm
N}\omega^{2}V^{2}}\sum_{I}\\!\sum_{\bm{K},\bm{K}^{\prime}}\\!\\!\frac{\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{K}^{\prime})}{q^{2}+m_{\phi}^{2}}\frac{\epsilon^{-1}_{\bm{K}\bm{K}^{\prime}}(\bm{k}_{e},\omega)}{|\bm{k}_{e}+\bm{K}||\bm{k}_{e}+\bm{K}^{\prime}|}Z(|\bm{k}_{e}+\bm{K}^{\prime}|)\\\
&\times\langle\lambda_{f}|e^{i(\bm{k}_{e}+\bm{K}^{\prime}-\bm{q})\cdot\bm{x}_{I}}|\lambda_{i}\rangle[\bm{p}_{e}+\bm{k}_{e}|e^{i(\bm{k}_{e}+\bm{K})\cdot\bm{x}_{e}}|\bm{p}_{e}]_{\Omega}\,,\end{split}$

where $\bm{q}\equiv\bm{p}_{i}-\bm{p}_{f}$ is the momentum released by the DM,
and in the second line we have used Bloch’s theorem for the electron state
Knapen  _et al._ (2021). In particular, $|\bm{p}_{e}]_{\Omega}$ represents the
(periodic) Bloch wave function (e.g., Ashcroft and Mermin, 2020), and the
matrix element between two electronic states is to be computed over the volume
of a primitive cell, $\Omega$, as indicated by the subscript.

At zero temperature, $|\lambda_{i}\rangle$ is simply the ground state of the
lattice Hamiltonian, while the initial electron state,
$|\bm{p}_{e}]_{\Omega}$, follows a distribution $f_{b}(\bm{p}_{e})$, with $b$
labeling the electronic branch. Fermi’s golden rule then returns the following
rate,

$\displaystyle\frac{d\Gamma}{d\omega}={}$
$\displaystyle\frac{4\alpha}{V^{3}}\bigg{(}\frac{g_{\chi}g_{\rm N}}{m_{\rm
N}\omega^{2}}\bigg{)}^{\\!\\!2}\sum_{\bm{q}}\sum_{\bm{k}_{e}}\sum_{\bm{K},\bm{K}^{\prime}}\sum_{\bm{Q},\bm{Q}^{\prime}}\frac{\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{K}^{\prime})\,\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{Q}^{\prime})}{\big{(}q^{2}+m_{\phi}^{2}\big{)}^{2}}\frac{Z(|\bm{k}_{e}+\bm{K}^{\prime}|)Z(|\bm{k}_{e}+\bm{Q}^{\prime}|)}{|\bm{k}_{e}+\bm{K}^{\prime}||\bm{k}_{e}+\bm{Q}^{\prime}|}$
$\displaystyle\times\frac{4\pi^{2}\alpha}{V}\sum_{b,b^{\prime}}f_{b}(\bm{p}_{e})\frac{[\bm{p}_{e}|e^{-i\bm{Q}\cdot\bm{x}_{e}}|\bm{p}_{e}+\bm{k}_{e}]_{\Omega}\,[\bm{p}_{e}+\bm{k}_{e}|e^{i\bm{K}\cdot\bm{x}_{e}}|\bm{p}_{e}]_{\Omega}}{|\bm{k}_{e}+\bm{K}||\bm{k}_{e}+\bm{Q}|}$
$\displaystyle\times\delta(\omega-
E_{b^{\prime},\bm{p}_{e}+\bm{k}_{e}}-E_{b,\bm{p}_{e}})\epsilon_{\bm{K}\bm{K}^{\prime}}^{-1}(\bm{k}_{e},\omega)\epsilon_{\bm{Q}\bm{Q}^{\prime}}^{-1*}(\bm{k}_{e},\omega)$
(10a)
$\displaystyle\times\sum_{\lambda_{f}}\,\langle\lambda_{i}|\sum_{I}e^{-i(\bm{k}_{e}+\bm{Q}^{\prime}-\bm{q})\cdot\bm{x}_{I}}|\lambda_{f}\rangle\langle\lambda_{f}|\sum_{J}e^{i(\bm{k}_{e}+\bm{K}^{\prime}-\bm{q})\cdot\bm{x}_{J}}|\lambda_{i}\rangle$
$\displaystyle\times(2\pi)\delta\left(E_{\bm{p}_{i}}-E_{\bm{p}_{f}}+E_{\lambda_{i}}-E_{\lambda_{f}}-\omega\right)$
$\displaystyle\begin{split}={}&\frac{4\alpha}{V^{3}}\bigg{(}\frac{g_{\chi}g_{\rm
N}}{m_{\rm
N}\omega^{2}}\bigg{)}^{\\!\\!2}\sum_{\bm{q}}\sum_{\bm{k}_{e}}\sum_{\bm{K},\bm{Q}}\frac{\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{K})\,\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{Q})}{\big{(}q^{2}+m_{\phi}^{2}\big{)}^{2}}\frac{Z(|\bm{k}_{e}+\bm{K}|)Z(|\bm{k}_{e}+\bm{Q}|)}{|\bm{k}_{e}+\bm{K}||\bm{k}_{e}+\bm{Q}|}\\\
&\times\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{Q}}^{-1}(\bm{k}_{e},\omega)\big{)}\sum_{\lambda_{f}}\,\langle\lambda_{i}|\sum_{I}e^{-i(\bm{k}_{e}+\bm{Q}-\bm{q})\cdot\bm{x}_{I}}|\lambda_{f}\rangle\\\
&\times\langle\lambda_{f}|\sum_{J}e^{i(\bm{k}_{e}+\bm{K}-\bm{q})\cdot\bm{x}_{J}}|\lambda_{i}\rangle(2\pi)\delta\left(E_{\bm{p}_{i}}-E_{\bm{p}_{f}}+E_{\lambda_{i}}-E_{\lambda_{f}}-\omega\right)\,,\end{split}$
(10b)

where $E_{b,\bm{p}_{e}}$ and $E_{b^{\prime},\bm{p}_{e}+\bm{k}_{e}}$ are the
initial and final electron energies, and $E_{\bm{p}_{i}}$ and $E_{\bm{p}_{f}}$
the initial and final DM energies. In the second equality we used Lindhard’s
formula for the imaginary part of the dielectric matrix Lindhard (1954); Adler
(1962); Knapen  _et al._ (2021), corresponding precisely to the second and
third lines of Eq. (10a).444The generalization of Lindhard’s formula to the
reciprocal lattice actually includes an extra term, associated to the
possibility of absorbing an initially free electron (cf. Appendix A in Knapen
_et al._ (2021)). In the ground state the density of initially free electrons
is negligible, and this term does not contribute. The function
$\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{Q}}^{-1}(\bm{k}_{e},\omega)\big{)}$
is the ELF.

Now, in the space of reciprocal lattice vectors, the ELF has small off-
diagonal elements Knapen  _et al._ (2021). For this reason, we can set
$\bm{Q}\simeq\bm{K}$ in the last two lines of Eq. (10b).555For DM particle
with $m_{\chi}\gtrsim 10$ MeV, this approximation becomes even better, since
one has $q\gg|\bm{k}_{e}+\bm{K}|,|\bm{k}_{e}+\bm{Q}|$ and one can neglect the
electron momenta in the last lines of Eq. (10b). The electron momentum
supported by the ELF is, in fact, never larger than around $20$ keV Knapen
_et al._ (2022) for valence electrons which dominate the ionization
probability for $m_{\chi}\lesssim 100$ MeV.  After doing that, we obtain our
final expression for the rate of Migdal emission in a semiconductor:

$\displaystyle\frac{d\Gamma}{d\omega}={}$ $\displaystyle\frac{8\pi\alpha
N_{\rm T}}{V^{3}}\bigg{(}\frac{g_{\chi}g_{\rm N}}{m_{\rm
N}\omega^{2}}\bigg{)}^{\\!\\!2}\sum_{\bm{q}}\sum_{\bm{k}_{e}}\sum_{\bm{K},\bm{Q}}\frac{\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{K})\,\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{Q})}{\big{(}q^{2}+m_{\phi}^{2}\big{)}^{2}}\frac{Z(|\bm{k}_{e}+\bm{K}|)Z(|\bm{k}_{e}+\bm{Q}|)}{|\bm{k}_{e}+\bm{K}||\bm{k}_{e}+\bm{Q}|}$
$\displaystyle\times\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{Q}}^{-1}(\bm{k}_{e},\omega)\big{)}S(\bm{q}-\bm{k}_{e}-\bm{K},E_{\bm{p}_{i}}-E_{\bm{p}_{f}}-\omega)\,,$
(11)

where $N_{\rm T}$ denotes the number of ions in the target, and $S(\bm{q},E)$
is the partial dynamic structure factor of the lattice (e.g., Squires, 1996),
i.e.

$\displaystyle S(\bm{q},E)\equiv\frac{1}{N_{\rm
T}}\sum_{\lambda_{f}}\bigg{|}\langle\lambda_{f}|\sum_{I}e^{i\bm{q}\cdot\bm{x}_{I}}|\lambda_{i}\rangle\bigg{|}^{2}\delta\big{(}E_{\lambda_{i}}-E_{\lambda_{f}}-E\big{)}\,.$
(12)

The advantage of expressing the rate in terms of the dynamic structure factor
lies in the fact that, for many-body systems, it encodes the complicated
inter-particle correlations, and only depends on the active degrees of freedom
in the target system that are excited given a momentum transfer $\bm{q}$
Trickle  _et al._ (2020). In principle, one could define and attempt to
measure a dynamic structure factor that incorporates the Migdal effect
directly, by considering electronic degrees of freedom in its definition. This
is practically challenging, both theoretically and experimentally. For this
reason, in our discussion we keep the lattice and electronic degrees of
freedom separated. The former are encoded in the structure factor, while the
latter in the ELF.

##  4 The dynamic structure factor

In principle, the dynamic structure factor can be measured through neutron
scattering experiments, in which case no further assumptions would be
necessary—Eq. (3.2) can be used without detailed knowledge of the lattice
dynamics. In practice, however, comprehensive data over the relevant range in
$\bm{q}$ ($10\text{ keV}-100\text{ keV}$) does not yet exist.666For
$q<10\text{ keV}$, the dynamic structure factor is dominated by elastic
scattering and the phonon density of states, which has been measured through
neutron scattering experiments as well as calculated in ab initio calculations
for silicon and germanium.

In Appendix B, we discuss the experimental prospects of taking those
measurements. In this section, instead, we present an overview of the
different approximation schemes that one can employ to model the structure
factor.

###  4.1 The incoherent and harmonic approximations

In the absence of data, we get as much analytical understanding as possible by
evaluating $S(\bm{q},E)$ in a simple context. Specifically, we use the
incoherent and the harmonic approximations. The former neglects interference
effects coming from the scattering off different lattice points. This is a
fair assumption as long as the momentum transfer is substantially larger than
the (inverse) inter-particle separation, a condition that is satisfied for all
DM masses considered here (see Figure 1). In particular, the incoherent
approximation implies

$\displaystyle S(\bm{q},E)\simeq\frac{1}{N_{\rm
T}}\sum_{\lambda_{f}}\sum_{I}\left|\langle\lambda_{f}|e^{i\bm{q}\cdot\bm{x}_{I}}|\lambda_{i}\rangle\right|^{2}\delta\big{(}E_{\lambda_{i}}-E_{\lambda_{f}}-E\big{)}\,.$
(13)

The harmonic approximation, on the other hand, neglects corrections to the
atomic potential, $U(\\{\bm{x}_{I}\\})$, that are higher than quadratic in the
positions, and it holds when the atomic displacements are small compared to
the inter-atomic spacing. This is only true for the lightest DM particles
(say, $m_{\chi}\lesssim 30$ MeV Campbell-Deem  _et al._ (2022)), while for DM
particles with rates dominated by higher momentum transfers, anharmonicities
can give large corrections Campbell-Deem  _et al._ (2022). This latter
assumption is hence employed only to gain some analytical understanding of the
physics, and to estimate the event rates in a simplified setting.

Let us briefly review how to construct the dynamic structure factor within the
approximations outlined above. The position of each atom in the lattice can be
decomposed as $\bm{x}_{I}=\bm{x}_{I}^{0}+\bm{u}_{I}$, where $\bm{x}_{I}^{0}$
is the equilibrium location of the $I$-th atom, and $\bm{u}_{I}$ its relative
displacement. When evaluated for an harmonic crystal, Eq. (12) simplifies to
(e.g., Squires, 1996)

$\displaystyle S^{\rm(harm)}(\bm{q},E)=\frac{1}{N_{\rm
T}}\sum_{I}\int_{-\infty}^{\infty}\frac{dt}{2\pi}e^{-iEt}e^{-2W(\bm{q})}e^{\langle\bm{q}\cdot\hat{\bm{u}}_{I}(0)\,\bm{q}\cdot\hat{\bm{u}}_{I}(t)\rangle}\,,$
(14)

where $\hat{\bm{u}}_{I}(t)$ is the displacement operator in the Heisenberg
picture, responsible for the creation or annihilation of a single phonon
degree of freedom. Moreover,
$W(\bm{q})\equiv\frac{1}{2}\big{\langle}(\bm{q}\cdot\hat{\bm{u}}_{0}(0))^{2}\big{\rangle}$
is the Debye–Waller factor. The correlators appearing in Eq. (14) can both be
determined in terms of the phonon density of states, $g(E_{\text{ph}})$,
normalized as $\int_{0}^{\infty}dE_{\text{ph}}\,g(E_{\text{ph}})=1$. In
particular, within the assumptions of identical atoms and cubic symmetry, the
unequal time correlator is independent of the lattice point and of the
direction of the exchanged momentum, and it is given by Schober (2014);
Marshall  _et al._ (1971),

$\displaystyle\langle{\bm{q}\cdot\hat{\bm{u}}_{I}(0)\,{\bm{q}}\cdot\hat{\bm{u}}_{I}(t)\rangle}=\frac{q^{2}}{2m_{\rm
N}}\int_{-\infty}^{\infty}dE_{\text{ph}}^{\prime}\frac{g(E_{\text{ph}}^{\prime})}{E_{\text{ph}}^{\prime}}\big{(}n(E_{\text{ph}}^{\prime})+1\big{)}e^{iE_{\text{ph}}^{\prime}t}\equiv\frac{q^{2}}{2m_{\rm
N}}f(t)\,,$ (15)

where $n(E_{\text{ph}})$ is a Bose–Einstein distribution. At the low
temperature environments necessary for DM direct detection $n\simeq 0$.777Many
experiments operate at cryogenic temperatures, where $n\simeq 0$ is an
excellent approximation; even for SENSEI, which operates at a much higher
temperature of $T=130\text{ K}$, one gets (for a typical phonon energy of
$\langle E_{\rm ph}\rangle=40$ meV) $n\simeq 0.03$, which is also negligible.
From the equation above it follows that the Debye–Waller factor is

$\displaystyle 2W({\bm{q}})=\frac{q^{2}}{2m_{\rm
N}}\int_{0}^{\infty}dE_{\text{ph}}^{\prime}\frac{g(E_{\text{ph}}^{\prime})}{E_{\text{ph}}^{\prime}}\big{(}2n(E_{\text{ph}}^{\prime})+1\big{)}\equiv\frac{q^{2}}{2m_{\rm
N}}\big{\langle}E_{\rm ph}^{-1}\big{\rangle}\,,$ (16)

where, from now on, we indicate averages over phonon energies as $\langle
E_{\rm ph}^{n}\rangle$, for some $n$.

####  4.1.1 Exact recursive procedure

To compute the structure function, one can expand the factor
$e^{\langle{\bm{q}}\cdot\hat{\bm{u}}_{I}(0)\,{\bm{q}}\cdot\hat{\bm{u}}_{I}(t)\rangle}$
order by order, i.e. Marshall  _et al._ (1971); Schober (2014); Berghaus  _et
al._ (2022); Campbell-Deem  _et al._ (2022),

$\displaystyle S^{\rm(harm)}({\bm{q}},E)$
$\displaystyle=\int_{-\infty}^{\infty}\frac{dt}{2\pi}e^{-iEt}e^{-2W({\bm{q}})}\sum_{p}\frac{1}{p!}\left(\frac{q^{2}}{2m_{\rm
N}}f(t)\right)^{p}$ (17)
$\displaystyle=e^{-2W(\bm{q})}\left(\delta(E)+\frac{q^{2}}{2m_{\rm
N}}\left(\frac{g(E)}{E}\big{(}n(E)+1\big{)}\right)+\sum_{p\geq
2}\frac{1}{p!}\left(\frac{q^{2}}{2m_{\rm N}}\right)^{p}T_{p}(E)\right)\,.$

Each term in the sum above corresponds to the contribution due to a final
state with definite number of phonons—specifically, $p$ of them. The density
of states depends on the material properties of the semiconductor, and can be
obtained either experimentally (through neutron scattering, for example Nelin
and Nilsson (1972); Kim  _et al._ (2020)) or from ab initio calculations
Petretto  _et al._ (2018). The functions $T_{p}(E)$ are, instead, determined
by the recursion relation Marshall  _et al._ (1971); Schober (2014); Berghaus
_et al._ (2022)

$\displaystyle
T_{p}(E)=\int_{-\infty}^{+\infty}dE^{\prime}\,T_{1}(E-E^{\prime})T_{p-1}(E^{\prime})\,,$
(18)

where $T_{1}(E)=g(E)\big{(}n(E)+1\big{)}/E$. This is computationally
expensive, especially for large phonon numbers.

####  4.1.2 Free-ion approximation

The free-ion approximation assumes that the incoming and outgoing state of the
nucleus are well described by a plane wave, in analogy to the treatment of a
free atom. Strictly speaking this is never accurate for a nucleus in a
quadratic potential, yet for energy deposition of $E_{r}\gg\langle E_{\rm
ph}\rangle$ the approximation performs well. In this approximation the dynamic
structure function is simply given by:

$\displaystyle S^{\rm(\text{free ion})}({\bm{q}},E)\simeq\delta(E-E_{r})\,.$
(19)

Note that both the impulse and the free-ion approximations are valid for
$E_{r}\gg\langle E_{\rm ph}\rangle$. However, as we will show in Eqs. (20) and
Eq. (21), the corrections to the former scale by $\langle E_{\rm ph}\rangle
t\sim\langle E_{\rm ph}\rangle/E_{r}$, while the latter is rather controlled
by $\Delta/E_{r}=\sqrt{\langle E_{\rm ph}\rangle/E_{r}}$. Therefore, as a
function of recoil energy, the free-ion approximation becomes a good one
parametrically slower than the impulse approximation.

####  4.1.3 Impulse approximation

When the energy given to the lattice is sufficiently large ($E\gg\langle
E_{\text{ph}}\rangle$), it is possible to further simplify the harmonic
structure factor in Eq. 17 by using the so-called impulse approximation (ia)
Marshall  _et al._ (1971); Schober (2014); Knapen  _et al._ (2021). The
integrand in Eq. (15) has support on a range of energies close to the typical
phonon energy, $\langle E_{\rm ph}\rangle$. If the energy released to the
crystal is $E$, the scattering process happens over time scales $t\sim 1/E$.
When $E\gg\langle E_{\rm ph}\rangle$ one can expand Eq. (15) around $t=0$, to
obtain

$\displaystyle f(t)=\langle E^{-1}_{\text{ph}}\rangle+it-\frac{\langle
E_{\text{ph}}\rangle t^{2}}{2}+\mathcal{O}\left(\langle E_{\rm
ph}\rangle^{2}t^{3}\right)\,.$ (20)

One can justify this procedure more rigorously through a steepest descent
analysis Schober (2014); Marshall  _et al._ (1971); Campbell-Deem  _et al._
(2022). When used in Eq. (14), this returns the dynamic structure factor in
the impulse approximation, where it is well described by a Gaussian envelope
centered around the recoil energy expected in an elastic recoil of a free
nucleus, $E_{r}=q^{2}/2m_{\rm N}$. Specifically,

$\displaystyle
S^{\rm(ia)}({\bm{q}},E)\simeq\frac{1}{\sqrt{2\pi\Delta^{2}}}e^{-\frac{(E-E_{r}(q))^{2}}{2\Delta^{2}}}\,,$
(21)

where the width of the Gaussian is determined by $\Delta^{2}\equiv
E_{r}\langle E_{\text{ph}}\rangle$.

The impulse approximation can be intuitively understood due to the following
fact (see also our discussion in Section 2): when the collision happens over
times much shorter than the typical time scale characterizing the lattice
dynamics (the inverse phonon energy), the final nucleus does not have time to
probe the lattice potential. Its wave function can therefore be treated as a
free plane wave. While this considerably simplifies the problem, it breaks
down for small momentum transfers, i.e., when $q^{2}/2m_{\rm N}\lesssim\langle
E_{\rm ph}\rangle$. For lower momenta, the full lattice potential becomes
important, and phonon excitations must be properly accounted for in the
dynamic structure factor.

####  4.1.4 Improved impulse approximation

Within our simplifying assumption of an harmonic crystal, there exists an
additional method to calculate the dynamic structure factor that combines the
two strategies outlined in Subsections 4.1.1 and 4.1.3, and allows us to
extend the analytical calculation to small momentum transfers. This procedure
eases the computational burden necessary to recursively compute high orders of
$T_{p}(E)$ (Eq. (18)). The following result is known in the neutron scattering
literature Marshall  _et al._ (1971), but has not yet been applied in the DM
direct detection context. The idea is that, following the standard WKB
approach, the wave function of sufficiently excited final states closely
resembles a plane wave. For such states one then recovers the conditions of
applicability of the impulse approximation, as explained above. The operative
procedure is then to keep the exact elastic, single phonon, and up to
$p_{\text{ex}}$ multi-phonon response in Eq. (17), and to only expand $f(t)$
in small $t$ for higher order terms, i.e.

$\displaystyle\begin{split}\sum_{p\geq
p_{\text{ex}}}\frac{1}{p!}\left(\frac{q^{2}}{2m_{\rm
N}}f(t)\right)^{p}&\simeq\sum_{p\geq
p_{\text{ex}}}\frac{\big{(}2W(\bm{q})\big{)}^{p}}{p!}\left(1+i{\big{\langle}E^{-1}_{\text{ph}}\big{\rangle}}^{-1}t-{\big{\langle}E^{-1}_{\text{ph}}\big{\rangle}}^{-1}\frac{\langle
E_{\text{ph}}\rangle t^{2}}{2}\right)^{p}\\\ &\equiv\sum_{p\geq
p_{\text{ex}}}\frac{\big{(}2W(\bm{q})\big{)}^{p}}{p!}e^{xp}\,,\end{split}$
(22)

where
$x=i{\big{\langle}E^{-1}_{\text{ph}}\big{\rangle}}^{-1}t-\frac{1}{2}\tilde{\Delta}^{2}t^{2}+\mathcal{O}\big{(}{\langle
E_{\rm ph}\rangle}^{3}t^{3}\big{)}$, with $\tilde{\Delta}^{2}\equiv\langle
E_{\rm ph}\rangle\big{\langle}E_{\rm
ph}^{-1}\big{\rangle}^{-1}-\big{\langle}E_{\rm
ph}^{-1}\big{\rangle}^{-2}$.888Note that $\tilde{\Delta}^{2}\geq 0$ is always
guaranteed by the Cauchy–Schwarz inequality. Writing the total structure
factor as $S(\bm{q},E)=\sum_{p}S_{p}(\bm{q},E)$, this leads to an improved
impulse approximation (iia), such that

$\displaystyle\begin{split}S_{p\geq p_{\rm
ex}}^{\rm(iia)}(\bm{q},E)\simeq{}&e^{-2W(\bm{q})}\frac{{\big{(}2W(\bm{q})\big{)}}^{p}}{p!}\frac{1}{\sqrt{2\pi
p\tilde{\Delta}^{2}}}e^{-\frac{\left(E-p/\langle E_{\rm
ph}^{-1}\rangle\right)^{2}}{2p\tilde{\Delta}^{2}}}\\\
\equiv{}&e^{-2W(\bm{q})}\frac{{\big{(}2W(\bm{q})\big{)}}^{p}}{p!}T^{(n,{\rm
iia})}_{p}(E)\,,\end{split}$ (23)

where we defined the normalized $T^{(n)}$ function as $T^{(n)}_{p}(E)\equiv
T_{p}(E)/\big{\langle}E_{\rm ph}^{-1}\big{\rangle}^{p}$. The total structure
function can now be calculated easily up to arbitrary order by combining Eqs.
(17) and (23).

Figure 3: Comparison of the (normalized) multi-phonon terms,
${T}^{(n)}_{p}(E)$, obtained through the exact recursion relation in Eq. (17)
(solid) and through the improved impulse approximation (iaa) (dashed) for
silicon. The improved impulse approximation performs poorly for terms with
$p\leq 3$ (left panel); however, already around $p\geq 4$ the approximation
provides a fairly accurate description. Eventually, for $p\geq 10$ the exact
multi-phonon terms no longer show any features that differ from a Gaussian,
and are well described by the improved impulse approximation (right panel).

To compare this improved impulse approximation with the exact recursive
procedure we compute the ${T}^{(n)}_{p}(E)$ evaluated with these two methods
for the case of silicon, see Figure 3. The results highlight how the exact
higher order terms calculated by the recursive procedure in Eq. (17) lose
sensitivity to the features of the material-specific phonon density of states.
For silicon, multi-phonon terms larger than $p=10$ no longer exhibit features
distinguishing them from a simple Gaussian, and even smaller multi-phonon
terms, down to $p=4$, are decently well captured by the improved impulse
approximation.999The slight difference in peak position between the improved
impulse approximation and the exact result is also present in the regular
impulse approximation and can be analytically accounted for Gunn and Warner
(1984); Liang  _et al._ (2022).

The key difference between the improved impulse approximation and the standard
one lies in keeping more information about the phonon density of states. While
the standard approximation only depends on $\langle E_{\text{ph}}\rangle$, the
procedure presented here also depends on
$\big{\langle}E^{-1}_{\text{ph}}\big{\rangle}$, which enters through the width
of the individual Gaussian multi-phonon terms, $\tilde{\Delta}^{2}$.
Additionally, by keeping the exact terms up to a fixed order, the structure
factor is guaranteed to be accurate for all momentum transfers, $\bm{q}$. In
the large momentum limit the total structure function is the weighted sum of
Gaussians centered around each multi-phonon term, which matches onto the
regular impulse approximation.

###  4.2 The dynamic structure factor in silicon

Figure 4: Dynamic structure factor, $S(q,E)$, in silicon as a function of
energy deposition, $E$, and for different values of the momentum transfer,
$q$. We show the changes in $S(q,E)$ as the momentum transfer transitions from
the small ($W(q)\ll 1$) to the large momentum ($W(q)\gg 1$) regime. To
illustrate the relative contribution of the elastic term in Eq. (17) we
convolve it with a Gaussian envelope, $S_{0}(q,E)=e^{-2W(q)}\delta(E)\to
e^{-2W(q)}e^{-E^{2}/2\delta^{2}}/\sqrt{2\pi\delta^{2}}$, where we set
$\delta^{2}=10^{-3}\tilde{\Delta}^{2}$. Note that for the bottom right figure,
$q=90$ keV, the purely elastic scattering with $E=0$ is negligible, but it
becomes increasingly important for lower momentum transfers.

As input for our rate calculations, we evaluate $S(\bm{q},E)$ in the harmonic
approximation using Eq. (17), where we recursively calculate the multi-phonon
terms up to order $p=10$, and then use the improved impulse approximation for
all higher order terms in the sum. The results are reported in Figure 4. As
one can see, when transitioning towards small momentum transfers, the multi-
phonon contribution becomes subdominant with respect to the elastic one,
$S_{0}(q,E)=e^{-2W(q)}\delta(E)$, where the crystal recoils as a whole. Such
elastic response has received less attention in the existing DM literature on
the Migdal effect (although see e.g. Kahn  _et al._ (2021); Blanco  _et al._
(2022)). However, it is possible for a scattering event to lead only to
ionization, without exciting any phonons in the crystal. In fact, for DM for
masses $m_{\chi}\lesssim 50\,\text{MeV}$ those events dominate over scattering
events that are double inelastic processes, i.e., processes that produce both
a Migdal ionization as well as phonons. We emphasize the importance of this
contribution, noting that previous work targeting the low mass regime did not
include it in their numerical evaluation of the differential rates, leading to
a large discrepancy between the estimates obtained with the free-ion
approximation compared to the exact result accounting for lattice dynamics
Liang  _et al._ (2022). The latest version of Liang  _et al._ (2022) now
presents the corrected result.

##  5 Results

After deriving the expression for the Migdal rate within an effective field
theory in Sec. 3 in terms of the dynamic structure factor, which we discuss in
detail in Sec. 4, we now apply our analytical result to calculate the full and
differential rates we can expect in upcoming experiments. To obtain the total
rate per target mass, $R$, we divide by the detector mass $M_{\rm T}$, average
over the DM’s initial velocity, and multiply by the number of DM particles in
the detector such that

$\displaystyle R=n_{\chi}\frac{V}{M_{\rm T}}\int d\omega\int
d^{3}vf_{\chi}(\bm{v})\frac{d\Gamma(\bm{v})}{d\omega}\,,$ (24)

where $n_{\chi}$ is the local DM number density, and $f_{\chi}$ is the DM’s
velocity distribution. To evaluate Eq. (24), we drop the off-diagonal terms,
$\bm{K}\neq\bm{Q}$, in Eq. (3.2) to reduce computational burden.101010The
combined neglected contribution is at most comparable to the diagonal terms
Kozaczuk and Lin (2020). We also take the infinite volume limit, and replace
our sums with integrals. The rate per target mass can then be written as,

$\displaystyle R={}$ $\displaystyle\frac{n_{\chi}}{M_{\rm T}}\int
d^{3}vf_{\chi}(\bm{v})\int
d\omega\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\sum_{\bm{K}}8\pi\alpha
N_{\rm T}\bigg{(}\frac{g_{\chi}g_{\rm N}}{m_{\rm
N}\omega^{2}}\bigg{)}^{\\!\\!2}$ (25) $\displaystyle\times
Z^{2}(|\bm{k}_{e}+\bm{K}|)\frac{(\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{K}))^{2}}{\big{(}q^{2}+m_{\phi}^{2}\big{)}^{2}}\frac{\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}}{|\bm{k}_{e}+\bm{K}|^{2}}S\left(\bm{q}-\bm{k}_{e}-\bm{K},\bm{q}\\!\cdot\\!\bm{v}-\tfrac{q^{2}}{2m_{\chi}}-\omega\right)\,.$

Assuming for both the velocity distribution and structure factor to be
isotropic (see Appendix C for details), this equation simplifies to

$\displaystyle\begin{split}R={}&\int d\omega\frac{\alpha
A^{2}\bar{\sigma}_{n}n_{\chi}N_{\rm T}}{m_{\rm N}^{2}\mu^{2}_{\chi
n}\omega^{4}M_{\rm T}}\int d^{3}v\frac{f_{\chi}(v)}{v}\int
dq\,q^{3}F^{2}_{\text{DM}}(q)\int{d\cos}\theta_{qk}\cos^{2}\theta_{qk}\\\
&\times\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\,\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}\\\
&\times\int_{0}^{E^{\text{max}}}dE\,S\left(\sqrt{q^{2}+|\bm{k}_{e}+\bm{K}|^{2}-2q|\bm{k}_{e}+\bm{K}|\cos\theta_{qk}}\,,E\right).\end{split}$
(26)

In particular, we replaced the couplings with a reference cross section,
$\bar{\sigma}_{n}$, as follows: $\left(g_{\chi}g_{\rm N}\right)^{2}=\pi
A^{2}\bar{\sigma}_{n}\big{(}q^{2}_{0}+m^{2}_{\phi}\big{)}^{2}/\mu_{\chi
n}^{2}$, where $q_{0}=m_{\chi}v_{0}$ is the typical DM momentum, $A$ is the
mass number, and $\mu_{\chi n}$ the reduced DM–nucleon mass. We also
introduced the form factor $F_{\rm
DM}(q)\equiv\big{(}q_{0}^{2}+m_{\phi}^{2}\big{)}/\big{(}q^{2}+m_{\phi}^{2}\big{)}$.
Here we assume that the scalar mediator couples equally to protons and
neutrons. One obtains the result for a vector rather than scalar mediator by
replacing $m_{\phi}$ with $m_{A}$ in Eq. (25). For a dark photon we would make
the replacement $\left(g_{\chi}g_{\rm N}\right)^{2}=\pi
Z^{2}\bar{\sigma}_{p}\big{(}q^{2}_{0}+m^{2}_{A}\big{)}^{2}/\mu_{\chi p}^{2}$.

The electron momentum appearing in the structure factor satisfies
$|\bm{k}_{e}+\bm{K}|\lesssim 20\text{ keV}$ Knapen  _et al._ (2022). For
sufficiently heavy DM ($m_{\chi}\gtrsim 10\text{ MeV}$) one can then neglect
it with respect to the momentum transfer, $q\gg|\bm{k}_{e}+\bm{K}|$. In this
regime one can evaluate the angular integral analytically, which simplifies
the rate to a factorized expression,

$\displaystyle R\simeq{}$ $\displaystyle\int d\omega\frac{2\alpha
A^{2}\bar{\sigma}_{n}n_{\chi}N_{\rm T}}{3m_{\rm N}^{2}\mu^{2}_{\chi
n}\omega^{4}M_{\rm
T}}\int_{{v}_{\text{min}}}^{{v}_{\text{max}}}d^{3}{v}\frac{f_{\chi}(v)}{v}\int\frac{d^{3}\bm{k_{e}}}{(2\pi)^{3}}\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}$
$\displaystyle\times\int_{q_{\text{min}}}^{q_{\text{max}}}{d}q\,q^{3}F^{2}_{\text{DM}}(q)\int_{0}^{E^{\text{max}}}dE\leavevmode\nobreak\
S(q,E)\,.$ (27)

Near the kinematic endpoint, for low momenta when $q\sim|\bm{k}_{e}+\bm{K}|$,
this factorization breaks down and one should evaluate the full multi-
dimensional integral in Eq. (26).

To numerically compute Eq. (26), we quantify the partial dynamic structure
factor as described in Subsection 4.1.4 and illustrated in Fig. 4. For the
isotropic DM velocity distribution, we assume a standard truncated
Maxwell–Boltzmann distribution, boosted by the Earth velocity with respect to
the galactic rest frame which we average over the angular distribution

$\displaystyle f_{\chi}(v)=\frac{1}{4\pi}\int
d\Omega_{\bm{v}}\frac{1}{\pi^{\frac{3}{2}}v^{2}_{0}\left(v_{0}\erf(\frac{v_{\text{esc}}}{v_{0}})-\frac{2v_{\text{esc}}}{\sqrt{\pi}}e^{-\frac{v^{2}_{\text{esc}}}{v^{2}_{0}}}\right)}e^{-\frac{(\bm{v}+\bm{v}_{\rm
e})^{2}}{v^{2}_{0}}}\Theta(v_{\text{esc}}-|\bm{v}+\bm{v}_{\rm e}|)\,,$ (28)

where we set $v_{\rm e}=240\text{ km/s}$, $v_{0}=220\text{ km/s}$, and
$v_{\text{esc}}=500\text{ km/s}$. We also take the DM number density to be
$n_{\chi}=0.4\text{ cm}^{-3}\,(1\text{ GeV}/m_{\chi})$. For the electron
momentum integral we use the data tables provided in DarkELF Knapen  _et al._
(2022), which specify $Z(|\bm{k_{e}}+\bm{K}|)$ and
$\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}$,
calculated via time-dependent density functional theory methods with the GPAW
package Mortensen  _et al._ (2005); Enkovaara  _et al._ (2010) for silicon and
germanium.

Figure 5: The differential DM--nucleus Migdal scattering rate expected for 100
g of material and 1 year of exposure, computed for a heavy (left panel) and a
light (right panel) mediator, for silicon (cyan) and germanium (red), and for
two different DM masses, $m_{\chi}=5\text{ MeV}$ (dashed), and
$m_{\chi}=10\text{ MeV}$ (solid). The DM--nucleon reference cross-section
defined below Eq. (26) is set to $\bar{\sigma}_{n}=10^{-36}\text{ cm}^{2}$.
The vertical lines indicate the approximate ionization threshold required to
produce two electron-hole pairs in the detector for silicon (dotted,
$5.0\text{ eV}$) and germanium (dot-dashed, $3.7\text{ eV}$) Essig  _et al._
(2016).

Then we use quasi-adaptive Monte Carlo integration to evaluate the multi-
dimensional integral over $\bm{k}_{e}$, $\cos\theta_{qk}$, $E$, $q$, $v$ and
$\omega$, with the following limits of integration,

$\displaystyle q_{\rm
min}=m_{\chi}v\left(1-\sqrt{1-\frac{2\omega}{m_{\chi}v^{2}}}\right)\,,\quad
q_{\rm max}=m_{\chi}v\left(1+\sqrt{1-\frac{2\omega}{m_{\chi}v^{2}}}\right)\,,$
$\displaystyle v_{\rm min}=\sqrt{2\omega/m_{\chi}}\,,\quad v_{\rm
max}=v_{e}+v_{\rm esc}\,,\quad E^{\rm
max}=qv-\frac{q^{2}}{2m_{\chi}}-\omega\,.$

Our results for the differential ionization rate, ${dR}/{d\omega}$, for
silicon and germanium are shown in Fig. 5. The differential ionization rate
peaks at energies slightly above the band gap, and then drops with increasing
$\omega$, an expected feature due to the scaling with respect to ionization
energy $dR/d\omega\propto 1/\omega^{4}$. The decrease is more pronounced for
lighter DM masses as light DM has diminished phase-space available to impart
large ioniziation energies. Experiments measure an ionization signal $Q$, the
number of electron-hole pairs (denoted as $Qe^{-}$) produced in an event,
rather than $\omega$ directly. We map from one onto the other using Essig  _et
al._ (2016)

$Q(\omega)=1+\frac{(\omega-E_{\text{gap}})}{\epsilon}\,,$ (29)

where we take $\epsilon$, the mean energy per electron-hole pair, to be
$3.8\,\text{eV}$ and $2.9\,\text{eV}$ for silicon and germanium, respectively.
We assume the band gap energies $E^{\text{Si}}_{\text{gap}}=1.2\,\text{eV}$
and $E^{\text{Ge}}_{\text{gap}}=0.67\,\text{eV}$.

In Fig. 6 we, instead, show projections for the exclusion limits on the DM--
nucleon cross-section, $\bar{\sigma}_{n}$, in the low mass region for a
100-g-year detector made out of silicon and germanium for a heavy and a light
mediator. The exposure corresponds to the current plans of the SENSEI
detector, which consists of a silicon target. The germanium line depicts a
hypothetical detector. For our projections, we assume zero background events
for the bins with $\geq 2e^{-}$, but $N_{\text{bgk}}=10^{5}$ single-electron
background events; we assume the same number of background events for SENSEI
as for the hypothetical germanium detector. This background estimate
corresponds to the number of single-electron events anticipated at SENSEI for
the displayed exposure. The projections are dominated by the background free
$2e^{-}$-threshold shown in Fig. 5 for masses $m_{\chi}\gtrsim 2\text{ MeV}$,
where we integrate the differential rate shown from the displayed
$2e^{-}$-threshold to $\sim 20\,\text{eV}$ to obtain the projected number of
events. Our $90\%$ C.L. corresponds to an upper limit of 2.4 events. For
smaller masses, the $1e^{-}$ bin eventually delivers dominant constraints
(seen as the kinks in the curves). We also recast the constraints from SENSEI
taken at Fermilab near the MINOS hall (‘‘SENSEI@MINOS’’) Barak  _et al._
(2020). Since Barak  _et al._ (2020) presented results for the $1e^{-}$,
$2e^{-}$, $3e^{-}$, and $4e^{-}$ bins, we calculate the limit from each of
these bins and show the best one in Fig. 6. We find that the limit is
dominated by the $2e^{-}$-bin, for which we show the cross-section
$\overline{\sigma}_{n}$ that leads to more than $9.3$ events per $2.1$ g-day
exposure (we expect an analysis that combines multiple bins to yield a
slightly stronger constraint).

The shaded bands in Fig. 6 correspond to the treatment in Knapen  _et al._
(2021, 2022), where the dynamic structure factor is computed within the
impulse approximation, and the momentum integral IR cutoff is varied between
$q_{\rm min}=2\sqrt{2m_{\rm N}\langle E_{\text{ph}}\rangle}$ and $q_{\rm
min}=3\sqrt{2m_{\rm N}\langle E_{\text{ph}}\rangle}$, the threshold
approximating the break down of the impulse approximation. We find that this
prescription underestimates the rate, since it artificially excludes the low-
momentum region, which is still above the kinematic threshold. The effect is
more prominent for the light mediator projections, as the rate scales as $1/q$
leading to a large enhancement precisely in the small-momentum region, which
then dominates the rate.111111We stress that this underestimation is
inherently due to imposing a cutoff on the momentum integral, rather than to
the particular approximation used to describe the structure factor. We find,
in fact, a similar underestimation of the full rate when using the free-ion
approximation with a cutoff. This cutoff-dependent rate suppression is
illustrated in the right panel of Fig. 6, comparing the impulse approximation
for silicon and germanium. Although germanium has a larger rate than silicon
(solid lines), the uncertainty band obtained by varying the cutoff indicates a
smaller rate than silicon. This is due to the cutoff $q_{\text{min}}$ being
larger for the heavier germanium atom, artificially suppressing the rate.

Figure 6: Projections using the Migdal effect of the $90\%$ C.L. on the DM--
nucleon cross-section for a heavy (left panel) and light (right panel)
mediator for SENSEI (100-g-year silicon detector) (cyan) and a hypothetical
100-g-year germanium detector (red) with $1e^{-1}$ thresholds (see text for
details). We find no difference in the total rate between using the exact
dynamic structure factor (Eq. (14)) and the one obtained in the free-ion
approximation (Eq. (19))---see text for details. For comparison we show
projections for ‘‘LBECA,’’ a proposed 100 kg-year xenon detector, with a
$2e^{-}$ threshold (green, dashed) Essig  _et al._ (2020). The shaded cyan and
red regions correspond, respectively, to the result for silicon and germanium
reported in Knapen  _et al._ (2021), where the authors used the impulse
approximation and imposed an IR cutoff on the momentum integral when the
approximation breaks down. The uncertainty band corresponds to varying the
cutoff from $q_{\rm min}=2\sqrt{2m_{\rm N}\langle E_{\text{ph}}\rangle}$ to
$q_{\rm min}=3\sqrt{2m_{\rm N}\langle E_{\text{ph}}\rangle}$. The gray-shaded
regions correspond to current constraints from XENON10 Angle  _et al._ (2011),
XENON1T Aprile  _et al._ (2019b), and a recast of XENON1T data for cosmic-ray
up-scattered DM Bringmann and Pospelov (2019). The cyan shaded region
indicates our recast of the SENSEI@MINOS Barak  _et al._ (2020) constraints.

Interestingly, we find that both the total rate, $R$, and the differential
rate, ${dR}/{d\omega}$, are insensitive to whether the dynamic structure
factor is computed within the simple free-ion approximation of Eq. (19), or
the more accurate harmonic one, as in Eqs. (17) and (23). We also notice that
$R$ and $dR/d\omega$ are insensitive (up to sub-percent corrections) to
whether one used the complete expression in Eq. (26) or the approximate one in
Eq. (5). This is true even for small DM masses, such that
$q\sim|\bm{k}_{e}+\bm{K}|$. The reason for this lies in the fact that, for
most of the parameter space of interest, the maximum energy that can be
deposited to the crystal is much larger than the region where the structure
factor has non-zero support, i.e., $E^{\rm max}\gg E_{\rm ph},E_{r}$.
Consequently, one can use the exact sum rule obeyed by the structure factor,
$\int_{0}^{\infty}dE\,S(\bm{q},E)=1$.121212One can check this analytically
using, for example, the free-ion expression for the structure factor. In this
case, the energy integral simply changes the limits of integration over the
momentum transfer, $q_{\rm
min/max}=m_{\chi}v\left(1\mp\sqrt{1-\frac{2\omega}{m_{\chi}v^{2}}}\right)\to\mu_{\chi\rm
N}v\left(1\mp\sqrt{1-\frac{2\omega}{\mu_{\chi\rm N}v^{2}}}\right)$. For
$m_{\chi}\lesssim 100\text{ MeV}$ the two coincide with very good accuracy.
Moreover, for small momentum transfer, when $E^{\rm max}\sim E_{\rm ph}$, the
rate is anyway dominated by the elastic term (where the entire lattice
recoils), which still satisfies the previous sum rule. Therefore, for the
entire kinematic region of interest, the event rate can be written in a simple
factorized form,

$\frac{dR}{d\omega}\left(\omega\right)=\frac{n_{\chi}N_{T}}{m_{N}M_{T}}\int^{q_{\text{max}}}_{q_{\text{min}}}dq\int
d^{3}vvf_{\chi}(v)\frac{d\sigma}{dq}\left(q\right)\frac{dP}{d\omega}\left(q,\omega\right)\,,$
(30)

in terms of the differential cross section

$\frac{d\sigma}{dq}=\frac{A^{2}\bar{\sigma}_{n}}{\mu^{2}_{\chi
n}v^{2}}qF^{2}_{\text{DM}}(q)\,,$ (31)

and the ionization or shake-off probability (see also Knapen  _et al._ (2021,
2022)),

$\frac{dP}{d\omega}\left(q,\omega\right)=\frac{2\alpha
q^{2}}{3\omega^{4}m_{N}}\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\,\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}\,.$
(32)

As we just showed, the event rates inclusive in the energy released to the
lattice, $E$, are insensitive to the detailed properties of the structure
factor, as long as the structure factor obeys the sum rule. Extrapolating the
impulse approximation $S^{(\text{ia})}(q,\omega)$ to small momentum transfers
underestimates the rate by up to a factor of two since
$\int_{0}^{\infty}dES^{(\text{ia})}(q\to 0,\omega)\simeq 0.5$.

Nonetheless, our expression (26) can be used also to extract additional
information about the exclusive observable of the differential energy
distribution deposited to the lattice ${dR}/{dE}$, which instead strongly
depends on the structure factor. Specifically,

$\displaystyle\begin{split}\frac{dR}{dE}&=\int d\omega\frac{2\alpha
A^{2}\bar{\sigma}_{n}n_{\chi}N_{T}}{m_{N}^{2}\mu^{2}_{\chi
n}\omega^{4}M_{T}}\int
d^{3}{v}\frac{f_{\chi}(v)}{v}\int{d}q\,q^{3}F^{2}_{\text{DM}}(q)\int^{1}_{-1}{d\cos}\theta_{qk}\cos^{2}\theta_{qk}\\\
&\times\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\,\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}\\\
&\times
S\left(\sqrt{q^{2}+|\bm{k}_{e}+\bm{K}|^{2}-2q|\bm{k}_{e}+\bm{K}|\cos\theta_{qk}}\,,E\right)\,\Theta(E_{\text{max}}-E)\,.\end{split}$
(33)

Next-generation detectors, e.g., those employing transition edge sensor
technology D. N. McKinsey et al. [TESSERACT Collaboration] , aim at measuring
phonon and multi-phonon energy depositions. Such detectors are expected to be
sensitive to both the ionization as well as the lattice-energy deposition due
to a DM particle hitting a nucleus in a semiconductor. We expect this multi-
channel signal to be a powerful tool to discriminate between a DM signal and
backgrounds (see also Kahn  _et al._ (2021)). In Fig. 7, we show the
differential rate ${dR}/{dE}$ for silicon, for the case of a heavy mediator.
In particular, we contrast the spectrum obtained using the free-ion
approximation with the result obtained from the harmonic dynamic structure
factor. The two spectra are drastically different, both in shape and in
maximum energy deposits allowed, with the more accurate harmonic approximation
predicting much larger ones. We also show that neglecting or not
$|\bm{k}_{e}+\bm{K}|$ in the structure factor generally has no appreciable
effects on the differential rates. This is particularly advantageous, since
neglecting the electron momentum in Eq. (33) makes its evaluation much less
computationally demanding.

Figure 7: We show the differential rate ${dR}/{dE}$, measurable by next-
generation detectors, as function of energy $E$ deposited to the lattice,
contrasting the harmonic dynamic structure factor with the free-ion
approximation for DM with masses $m_{\chi}=10\text{ MeV}$ (left panel) and
$m_{\chi}=5\text{ MeV}$ (right panel). We set the reference cross section to
$\bar{\sigma}_{n}=10^{-36}\text{ cm}^{2}$ and the minimum electron energy to
$\omega_{\rm min}=5\text{ eV}$. The harmonic dynamic structure factor
accurately captures the non-zero probability of phonon excitations energies
above the allowed threshold in the free-ion approximation, illustrating the
break-down of the free-ion approximation. Neglecting the electron momentum in
the dynamic structure function results only in small differences even for the
lightest DM, $m_{\chi}=5\text{ MeV}$, where the electron momentum is of order
of the total momentum transfer. Consistent with Fig. 3 we convolve the elastic
scatterings with no energy deposition $S^{\text{(harm)}}_{0}=\delta(E)$ with a
Gaussian to illustrate the shift from phonon depositions to elastic
scatterings off the crystal.

##  6 Summary and discussion

When searching for sub-GeV dark matter with semiconductor targets, the
distinction between high- and low-energy physics is not sharp anymore,
requiring a framework that can quantify ionization
$\mathcal{O}(\gtrsim\text{eV})$ and phonon $\mathcal{O}(\sim 10\,\text{meV})$
signals at the same time. A robust theoretical understanding of these
detection channels is key for ensuring the continued success of experiments
such as SENSEI, DAMIC-M, SuperCDMS, CDEX, and others, as well as future
experiments such as Oscura Aguilar-Arevalo  _et al._ (2022) and SPICE D. N.
McKinsey et al. [TESSERACT Collaboration] . In this work, we extend the regime
of validity of the description of the Migdal effect in semiconductors to the
lowest mass regime ($m_{\chi}=1\,-50\text{ MeV}$) that is kinematically able
to give rise to ionization. We accurately incorporate the bound nature of the
nucleus in a crystal by taking advantage of the large separation of scales
between the final electron energy and the typical energies of the excitations
of the crystal lattice, formulating the problem in an effective field theory.
The effective field theory allows us to write the rate for Migdal emission in
a way that is universal and independent of the detailed knowledge of the
inter-atomic interactions. Specifically, the vibrational modes of the crystal
are all encoded in the dynamic structure factor, while the electronic dynamics
is encoded in the energy loss function. Both quantities are a priori directly
measurable. Comprehensive data quantifying the dynamic structure factor is not
yet available but is obtainable through neutron scattering experiments. In the
absence of data we utilize the harmonic approximation, which neglects higher
order corrections to the atomic potential, to illustrate how to model the
dynamic structure function in all regimes.

Whether the free-ion approximation accurately captures the Migdal rate in the
low mass regime had been an open question in the literature Knapen  _et al._
(2021); Liang  _et al._ (2021, 2022). We find that the exact inclusive rate
expression for the Migdal effect in a semiconductor agrees with the result one
obtains by extrapolating the free-ion approximation, which we prove to be a
consequence of an exact sum rule obeyed by the structure factor. When studying
the differential rate in the energy deposited to the crystal, $dR/dE$, the
above sum rule does not apply, and the harmonic structure factor predicts
energy depositions to the lattice substantially larger than what is found from
the naïve free-ion approximation. In the future, the small energy deposits to
the lattice may become an accessible observable, in which case it will be
imperative to use the exact result we have derived in this work.

###### Acknowledgements.

We are grateful to Gabriel Cuomo, Yoni Kahn, Simon Knapen, Tongyan Lin, Guido
Menichetti, Aditya Parikh, and Mauro Valli for useful discussions. We thank
Tien-Tien Yu and Joseph Pradler for early discussions on how to think about
the Migdal effect in semiconductors in the low-momentum regime. We also thank
Duncan Adam, Yoni Kahn, and Simon Knapen for comments on the manuscript. K.B.
thanks the Department of Energy’s Institute for Nuclear Theory at the
Universtiy of Washington for its hospitality and the Department of Energy for
partial support during the completion of this work, as well as the opportunity
to present and discuss with the participants of the Dark Matter in Compact
Objects, Stars, and in Low Energy Experiments workshop. A.E. and R.E. are
grateful to the organizers of the Pollica Summer Workshop supported by the
Regione Campania, Università degli Studi di Salerno, Università degli Studi di
Napoli ‘‘Federico II’’, the Physics departments ‘‘Ettore Pancini’’ and ‘‘E R
Caianiello’’, and the INFN. K.B. acknowledges the support of NSF Award
PHY1915093. For most of the development of this work A.E. has been a Roger
Dashen Member at the Institute for Advanced Study, whose work was also
supported by the U.S. Department of Energy, Office of Science, Office of High
Energy Physics under Award No. DE-SC0009988. R.E. acknowledges support from
DoE Grant DE-SC0009854, Simons Investigator in Physics Award 623940, the US-
Israel Binational Science Foundation Grant No. 2016153, and the Heising-Simons
Foundation Grant No. 79921. M.S. acknowledges support from Department of
Energy Grants DE-SC0009919 and DE-SC0022104.

##  Appendix A The electron--lattice Hamiltonian

In this Appendix, we show how to derive the general interaction between
valence electrons and the rest of the lattice, Eq. (7b). In our discussion, we
treat core and valence electrons separately, where the former are localized
around each nucleus, while the latter are not. The dielectric function is
defined through the response of the material to a given electric field.
Specifically, under the application of an external field, $\bm{E}_{\rm ext}$,
the material responds with a macroscopic field Jackson (2012),

$\displaystyle\bm{E}(\bm{x},\omega)=\int
d^{3}\bm{y}\,\epsilon^{-1}(\bm{x},\bm{y},\omega)\bm{E}_{\rm
ext}(\bm{y},\omega)\,,$ (34)

or, in Fourier space,131313We are defining the Fourier transform of the
dielectric function as $\epsilon^{-1}(\bm{p},\bm{q},\omega)=V^{-1}\int
d^{3}\bm{x}d^{3}\bm{y}\,\epsilon^{-1}(\bm{x},\bm{y},\omega)e^{i\bm{p}\cdot\bm{x}-i\bm{q}\cdot\bm{y}}$.

$\displaystyle\bm{E}(\bm{p},\omega)=\sum_{\bm{q}}\epsilon^{-1}(\bm{p},\bm{q},\omega)\bm{E}_{\rm
ext}(\bm{q},\omega)\,.$ (35)

We work for simplicity in the isotropic approximation, where the dielectric
function is just a scalar (rather than a tensor). The general case can be
treated in a similar way.

The dielectric function must be invariant under discrete translations by a
lattice vector, say $\bm{a}$:
$\epsilon^{-1}(\bm{x}+\bm{a},\bm{y}+\bm{a},\omega)=\epsilon^{-1}(\bm{x},\bm{y},\omega)$.
For its Fourier transform, $\epsilon^{-1}(\bm{p},\bm{q},\omega)$, this implies
that the difference $\bm{p}-\bm{q}$ is a reciprocal lattice vector.
Alternatively, one can write $\bm{p}=\bm{k}+\bm{K}$ and
$\bm{q}=\bm{k}+\bm{K}^{\prime}$, with $\bm{K}$ and $\bm{K}^{\prime}$ belonging
to the reciprocal lattice and $\bm{k}$ limited to the first Brillouin zone.
Therefore, Eq. (35) can be rewritten as

$\displaystyle\begin{split}\bm{E}(\bm{k}+\bm{K},\omega)={}&\sum_{\bm{K}^{\prime}}\epsilon^{-1}(\bm{k}+\bm{K},\bm{k}+\bm{K}^{\prime},\omega)\bm{E}_{\rm
ext}(\bm{k}+\bm{K}^{\prime},\omega)\,.\end{split}$ (36)

In order to find the Hamiltonian for a valence electron, we need to compute
the electric potential generating the field above. In momentum space, electric
field and potential, $\varphi$, are related by
$\bm{E}(\bm{p})=i\bm{p}\,\varphi(\bm{p})$, and from the equation above one
gets,

$\displaystyle\begin{split}\varphi(\bm{k}+\bm{K},\omega)={}&\sum_{\bm{K}^{\prime}}\frac{(\bm{k}+\bm{K})\cdot(\bm{k}+\bm{K}^{\prime})}{{|\bm{k}+\bm{K}|}^{2}}\epsilon^{-1}(\bm{k}+\bm{K},\bm{k}+\bm{K}^{\prime},\omega)\varphi_{\rm
ext}(\bm{k}+\bm{K}^{\prime},\omega)\\\ ={}&\sum_{\bm{K}^{\prime}}\epsilon_{\rm
LL}^{-1}(\bm{k}+\bm{K},\bm{k}+\bm{K}^{\prime},\omega)\varphi_{\rm
ext}(\bm{k}+\bm{K}^{\prime},\omega)\,,\end{split}$ (37)

where we have introduced the so-called longitudinal dielectric
function,141414We are grateful to Simon Knapen and Tongyan Lin for pointing
this quantity out to us. $\epsilon_{\rm
LL}^{-1}(\bm{p},\bm{q},\omega)\equiv(\bm{p}\cdot\bm{q}/p^{2})\epsilon^{-1}(\bm{p},\bm{q},\omega)$.
We now consider the nuclei and the core electrons as being a source of an
external potential which, following standard electromagnetism, is given by

$\displaystyle\varphi_{\rm ext}(\bm{x})=\frac{1}{4\pi}\int
d^{3}\bm{y}\frac{\sum_{I}\rho(|\bm{y}-\bm{x}_{I}|)}{|\bm{y}-\bm{x}|}\quad\Longrightarrow\quad\varphi_{\rm
ext}(\bm{p})=\sum_{I}\frac{\rho(p)}{p^{2}}e^{i\bm{p}\cdot\bm{x}_{I}}\,,$ (38)

with $\rho$ the charge density, which we assume is isotropic and centered
around each ion in the crystal. Introducing the symmetric dielectric matrix,
$\epsilon^{-1}_{\bm{K}\bm{K}^{\prime}}(\bm{k},\omega)\equiv\epsilon^{-1}_{\rm
LL}(\bm{k}+\bm{K},\bm{k}+\bm{K}^{\prime},\omega)|\bm{k}+\bm{K}|/|\bm{k}+\bm{K}^{\prime}|$,
the macroscopic electric potential reads

$\displaystyle\varphi(\bm{k}+\bm{K},\omega)=\sum_{I}\sum_{\bm{K}^{\prime}}\frac{\epsilon^{-1}_{\bm{K}\bm{K}^{\prime}}(\bm{k},\omega)}{|\bm{k}+\bm{K}||\bm{k}+\bm{K}^{\prime}|}\rho(|\bm{k}+\bm{K}^{\prime}|)e^{i(\bm{k}+\bm{K}^{\prime})\cdot\bm{x}_{I}}\,.$
(39)

To obtain the Hamiltonian of a valence electron we go back to real space, and
multiply the potential by the electron charge, $-e$. Defining also
$\rho(k)\equiv eZ(k)$, the result is

$\displaystyle H_{e\rm L}={}$
$\displaystyle-\frac{4\pi\alpha}{V}\sum_{I}\sum_{\bm{K},{\bm{K}}^{\prime}}\sum_{\bm{k}}\frac{\epsilon_{\bm{K}\bm{K}^{\prime}}^{-1}(\bm{k},\omega)Z(|\bm{k}+\bm{K}^{\prime}|)}{|\bm{k}+\bm{K}||\bm{k}+\bm{K}^{\prime}|}e^{i(\bm{k}+\bm{K})\cdot\bm{x}_{e}}e^{-i(\bm{k}+\bm{K}^{\prime})\cdot\bm{x}_{I}}\,.$
(40a)

##  Appendix B The dynamic structure factor from neutron spectroscopy

We expressed the Migdal scattering rate as Eq. (3.2), which depends on the
dynamic structure factor, $S(q,E)$, which incorporates the phonon degrees of
freedom in the material. While we discussed several simplifying assumptions
for $S(q,E)$ in Sec. 4, we discuss here how inelastic neutron scattering in
crystals can probe $S(q,E)$. Typically, the neutron scattering cross section
in a monatomic crystal can be written as Schober (2014),

$\displaystyle\frac{d^{2}\sigma}{d\Omega
dE}=\frac{\sigma_{a}}{4\pi}\frac{k_{f}}{k_{i}}S(q,E)\,,$ (41)

where $q$ and $E$ are the momentum and energy lost by the neutron, $k_{i}$
($k_{f}$) is the initial (final) momentum of the neutron, and $\sigma_{a}$ is
the total scattering cross section of a neutron with an individual atom in the
crystal. Thus, a measurement of the neutron scattering rate for a
monochromatic neutron beam in a crystal as a function of the final energy and
the direction of the neutron can directly provide a measurement of the
structure factor.

This measurement can be performed at a neutron spectroscopy experiment Hudson
(2006). Typically, a neutron beam is incident on the target material, and the
scattered neutrons are detected within a certain angular region. The energy
and the momentum lost by the initial neutron can be determined from the
measurement of the time-of-flight and the angular deflection caused by the
scattering.

Figure 8: Blue shaded region: The energy-momentum phase space that can be
probed kinematically with a 1.5 eV neutron beam at ARCS for a silicon sample.
Green shaded region: The relevant energy-momentum parameter space where the
dynamic structure factor in silicon is non-zero. Also shown are the momenta
that correspond to the Debye--Waller factor, $W(q)$, equal to 0.2, 1, and 5.

Here we consider the example of the spectrometer ARCS Abernathy  _et al._
(2012). ARCS can generate a monochromatic neutron beam with typical initial
energies in the range 15 meV -- 1500 meV. The angular range in the horizontal
plane is from $0\degree$ (forward scattering) to $135\degree$ (back-
scattering). Assuming perfect resolution on the final neutron energy, the
energy-momentum phase space that can be probed kinematically with a 1.5 eV
neutron beam at ARCS for a silicon sample is shown in Fig. 8. We also show the
region in energy-momentum space for which the dynamic structure factor in
silicon is expected to be appreciably different from zero. We estimate this
domain using the simple Gaussian description obtained within the impulse
approximation, as in Eq. (21). We also show the momenta that correspond to the
Debye--Waller factor $W$ equal to 0.2, 1, and 5, which represent the
transition from the harmonic regime to the impulse regime.

We see that with a 1.5 eV neutron, the energy-momentum parameter space where
the dynamic structure factor is non-zero in silicon can be probed up to
momenta of $\sim 100\leavevmode\nobreak\ \text{keV}$. Importantly, this covers
the region where the Debye--Waller factor transitions from values much smaller
than 1, where the dynamics are dominated by the elastic response of the entire
crystal, to values much higher than 1, where the dynamics are dominated by the
impulse approximation and the elastic nuclear recoil. The transition region
between these two extremes is susceptible to theory uncertainties and
anharmonic effects. Thus, a direct measurement of the dynamic structure factor
in this region would be valuable.

##  Appendix C The isotropic approximation

In this Appendix, we derive Eq. (26). We start with the rate per target mass
given in Eq. (25), and assuming an isotropic velocity distribution
$f_{\chi}(v)$, we first separate the integral over $\bm{v}$ into radial and
angular components. This gives,

$\displaystyle\begin{split}R=&\frac{n_{\chi}}{M_{\rm
T}}\int(2\pi)v^{2}dvf_{\chi}(v)\int
d\omega\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\,8\pi\alpha
N_{\rm T}\bigg{(}\frac{g_{\chi}g_{\rm N}}{m_{\rm N}\omega^{2}}\bigg{)}^{2}\\\
&\times\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\frac{(\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{K}))^{2}}{\big{(}q^{2}+m_{\phi}^{2}\big{)}^{2}}\frac{\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}}{|\bm{k}_{e}+\bm{K}|^{2}}\\\
&\times\int d\cos{\theta_{vq}}\leavevmode\nobreak\
S\left(\bm{q}-\bm{k}_{e}-\bm{K},qv\cos{\theta_{vq}}-\tfrac{q^{2}}{2m_{\chi}}-\omega\right)\,.\end{split}$
(42)

Changing variables in the last integral to the phonon energy,
$E=qv\cos{\theta_{vq}}-\frac{q^{2}}{2m_{\chi}}-\omega$, we get,

$\displaystyle\begin{split}R=&\frac{n_{\chi}}{M_{\rm
T}}\int(2\pi)v^{2}dvf_{\chi}(v)\int
d\omega\int\frac{d^{3}\bm{q}}{(2\pi)^{3}}\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\,8\pi\alpha
N_{\rm T}\bigg{(}\frac{g_{\chi}g_{\rm N}}{m_{\rm N}\omega^{2}}\bigg{)}^{2}\\\
&\times\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\frac{(\bm{q}\\!\cdot\\!(\bm{k}_{e}+\bm{K}))^{2}}{\big{(}q^{2}+m_{\phi}^{2}\big{)}^{2}}\frac{\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}}{|\bm{k}_{e}+\bm{K}|^{2}}\\\
&\times\frac{1}{qv}\int_{0}^{E^{\text{max}}}dE\leavevmode\nobreak\
S(\bm{q}-\bm{k}_{e}-\bm{K},E)\,,\end{split}$ (43)

where $E^{\text{max}}$ is the maximum possible phonon energy given by,

$\displaystyle E^{\text{max}}=qv-\frac{q^{2}}{2m_{\chi}}-\omega\,.$ (44)

We then also split the integral over $\bm{q}$ into radial and angular
components, and get,

$\displaystyle R={}$ $\displaystyle\frac{n_{\chi}}{M_{\rm
T}}\int(2\pi)v^{2}dvf_{\chi}(v)\int
d\omega\int\frac{(2\pi)q^{2}dq}{(2\pi)^{3}}\int_{-1}^{1}d\cos{\theta_{qk}}\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\,8\pi\alpha
N_{\rm T}\bigg{(}\frac{g_{\chi}g_{\rm N}}{m_{\rm N}\omega^{2}}\bigg{)}^{2}$
$\displaystyle\times\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\frac{q^{2}|\bm{k}_{e}+\bm{K}|^{2}\cos^{2}{\theta_{qk}}}{\big{(}q^{2}+m_{\phi}^{2}\big{)}^{2}}\frac{\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}}{|\bm{k}_{e}+\bm{K}|^{2}}$
(45a)
$\displaystyle\times\frac{1}{qv}\int_{0}^{E^{\text{max}}}dE\leavevmode\nobreak\
S(\bm{q}-\bm{k}_{e}-\bm{K},E)$ $\displaystyle={}$
$\displaystyle\frac{n_{\chi}}{M_{\rm T}}\int(2\pi)vdvf_{\chi}(v)\int
d\omega\int\frac{q^{3}dq}{(2\pi)^{2}}\int_{-1}^{1}d\cos{\theta_{qk}}\cos^{2}{\theta_{qk}}\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\,8\pi\alpha
N_{\rm T}\bigg{(}\frac{g_{\chi}g_{\rm N}}{m_{\rm N}\omega^{2}}\bigg{)}^{2}$
$\displaystyle\times\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\frac{\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}}{\big{(}q^{2}+m_{\phi}^{2}\big{)}^{2}}\int_{0}^{E^{\text{max}}}dE\leavevmode\nobreak\
S(\bm{q}-\bm{k}_{e}-\bm{K},E)\,.$ (45b)

Now, we replace $(g_{\chi}g_{\rm N})^{2}/(q_{0}^{2}+m_{\phi}^{2})^{2}={\pi
A^{2}\bar{\sigma}_{n}}/{\mu_{\chi n}^{2}}$, where $\bar{\sigma}_{n}$ is the
reference DM--nucleon cross section evaluated at $q_{0}=m_{\chi}v_{0}$, and
the momentum-dependence of the interaction is encoded in the DM form factor
$F_{\text{DM}}(q)=(q_{0}^{2}+m_{\phi}^{2})/(q^{2}+m_{\phi}^{2})$. This gives

$\displaystyle R=$ $\displaystyle\int d\omega\frac{2\alpha
A^{2}\bar{\sigma}_{n}n_{\chi}N_{\rm T}}{m_{\rm N}^{2}\mu^{2}_{\chi
n}\omega^{4}M_{\rm T}}\int(2\pi)vdvf_{\chi}(v)\int d\omega\int
q^{3}dq\leavevmode\nobreak\
F^{2}_{\text{DM}}(q)\int_{-1}^{1}d\cos{\theta_{qk}}\cos^{2}{\theta_{qk}}$
$\displaystyle\times\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}\int_{0}^{E^{\text{max}}}dE\leavevmode\nobreak\
S(\bm{q}-\bm{k}_{e}-\bm{K},E)\,.$ (46)

Finally, restoring the full integral over $\bm{v}$ assuming an isotropic
$f_{\chi}(v)$, and also assuming the dynamic structure factor only depends on
the magnitude of the momentum argument, we get

$\displaystyle\begin{split}R={}&\int d\omega\frac{\alpha
A^{2}\bar{\sigma}_{n}n_{\chi}N_{\rm T}}{m_{\rm N}^{2}\mu^{2}_{\chi
n}\omega^{4}M_{\rm T}}\int
d^{3}v\frac{f_{\chi}(v)}{v}\int{d}q\,q^{3}F^{2}_{\text{DM}}(q)\int^{1}_{-1}{d\cos}\theta_{qk}\cos^{2}\theta_{qk}\int\frac{d^{3}\bm{k}_{e}}{(2\pi)^{3}}\\\
&\times\sum_{\bm{K}}Z^{2}(|\bm{k}_{e}+\bm{K}|)\text{Im}\big{(}\\!-\epsilon_{\bm{K}\bm{K}}^{-1}(\bm{k}_{e},\omega)\big{)}\\\
&\times\int_{0}^{E^{\text{max}}}dE\,S\left(\sqrt{q^{2}+|\bm{k_{e}}+\bm{K}|^{2}-2q|\bm{k_{e}}+\bm{K}|\cos\theta_{qk}}\,,E\right)\,,\end{split}$
(47)

which is Eq. (26).

## References

  * Aprile  _et al._ (2018) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Aalbers  _et al._ (2022) J. Aalbers _et al._ (LZ), (2022), arXiv:2207.03764 [hep-ex] .
  * Meng  _et al._ (2021) Y. Meng _et al._ (PandaX-4T), Phys. Rev. Lett. 127, 261802 (2021), arXiv:2107.13438 [hep-ex] .
  * Alkhatib  _et al._ (2021) I. Alkhatib _et al._ (SuperCDMS), Phys. Rev. Lett. 127, 061801 (2021), arXiv:2007.14289 [hep-ex] .
  * Abdelhameed  _et al._ (2019) A. Abdelhameed _et al._ (CRESST), Phys. Rev. D 100, 102002 (2019), arXiv:1904.00498 [astro-ph.CO] .
  * Essig  _et al._ (2012) R. Essig, J. Mardon, and T. Volansky, Phys. Rev. D 85, 076007 (2012), arXiv:1108.5383 [hep-ph] .
  * Kouvaris and Pradler (2017) C. Kouvaris and J. Pradler, Phys. Rev. Lett. 118, 031803 (2017), arXiv:1607.01789 [hep-ph] .
  * Ibe  _et al._ (2018) M. Ibe, W. Nakano, Y. Shoji, and K. Suzuki, JHEP 03, 194 (2018), arXiv:1707.07258 [hep-ph] .
  * Migdal (1939a) A. Migdal, Sov. Phys. JETP 9, 1163 (1939a).
  * Migdal (1939b) A. Migdal, ZhETF 9, 1163 (1939b).
  * Migdal and Ballentine (1978) A. B. Migdal and L. E. Ballentine, Physics Today 31, 60 (1978).
  * Vergados and Ejiri (2005) J. D. Vergados and H. Ejiri, Phys. Lett. B606, 313 (2005), arXiv:hep-ph/0401151 [hep-ph] .
  * Moustakidis  _et al._ (2005) C. C. Moustakidis, J. D. Vergados, and H. Ejiri, Nucl. Phys. B 727, 406 (2005), arXiv:hep-ph/0507123 .
  * Bernabei  _et al._ (2007) R. Bernabei _et al._ , Int. J. Mod. Phys. A 22, 3155 (2007), arXiv:0706.1421 [astro-ph] .
  * Akerib  _et al._ (2019) D. S. Akerib _et al._ (LUX), Phys. Rev. Lett. 122, 131301 (2019), arXiv:1811.11241 [astro-ph.CO] .
  * Liu  _et al._ (2019) Z. Z. Liu _et al._ (CDEX), Phys. Rev. Lett. 123, 161301 (2019), arXiv:1905.00354 [hep-ex] .
  * Essig  _et al._ (2020) R. Essig, J. Pradler, M. Sholapurkar, and T.-T. Yu, Phys. Rev. Lett. 124, 021801 (2020), arXiv:1908.10881 [hep-ph] .
  * Armengaud  _et al._ (2019) E. Armengaud _et al._ (EDELWEISS), Phys. Rev. D 99, 082003 (2019), arXiv:1901.03588 [astro-ph.GA] .
  * Aprile  _et al._ (2019a) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 241803 (2019a), arXiv:1907.12771 [hep-ex] .
  * Barak  _et al._ (2020) L. Barak _et al._ (SENSEI), Phys. Rev. Lett. 125, 171802 (2020), arXiv:2004.11378 [astro-ph.CO] .
  * Knapen  _et al._ (2022) S. Knapen, J. Kozaczuk, and T. Lin, Phys. Rev. D 105, 015014 (2022), arXiv:2104.12786 [hep-ph] .
  * Adhikari  _et al._ (2022) G. Adhikari _et al._ (COSINE-100), Phys. Rev. D 105, 042006 (2022), arXiv:2110.05806 [hep-ex] .
  * Al-Bakry  _et al._ (2022) M. Al-Bakry _et al._ (SuperCDMS), (2022), arXiv:2203.02594 [hep-ex] .
  * Armengaud  _et al._ (2022) E. Armengaud _et al._ (EDELWEISS), (2022), arXiv:2203.03993 [astro-ph.GA] .
  * Agnes  _et al._ (2022) P. Agnes _et al._ (DarkSide), (2022), arXiv:2207.11967 [hep-ex] .
  * Essig  _et al._ (2022) R. Essig, G. K. Giovanetti, N. Kurinsky, D. McKinsey, K. Ramanathan, K. Stifter, and T.-T. Yu, in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.08297 [hep-ph] .
  * Nakamura  _et al._ (2021) K. D. Nakamura, K. Miuchi, S. Kazama, Y. Shoji, M. Ibe, and W. Nakano, PTEP 2021, 013C01 (2021), arXiv:2009.05939 [physics.ins-det] .
  * Liao  _et al._ (2021) J. Liao, H. Liu, and D. Marfatia, Phys. Rev. D 104, 015005 (2021), arXiv:2104.01811 [hep-ph] .
  * Bell  _et al._ (2022) N. F. Bell, J. B. Dent, R. F. Lang, J. L. Newstead, and A. C. Ritter, Phys. Rev. D 105, 096015 (2022), arXiv:2112.08514 [hep-ph] .
  * Araújo  _et al._ (2022) H. M. Araújo _et al._ , (2022), arXiv:2207.08284 [hep-ex] .
  * Cox  _et al._ (2022) P. Cox, M. J. Dolan, C. McCabe, and H. M. Quiney, (2022), arXiv:2208.12222 [hep-ph] .
  * Adams  _et al._ (2022) D. Adams, D. Baxter, H. Day, R. Essig, and Y. Kahn, (2022), arXiv:2210.04917 [hep-ph] .
  * Dolan  _et al._ (2018) M. J. Dolan, F. Kahlhoefer, and C. McCabe, Phys. Rev. Lett. 121, 101801 (2018), arXiv:1711.09906 [hep-ph] .
  * Bell  _et al._ (2020) N. F. Bell, J. B. Dent, J. L. Newstead, S. Sabharwal, and T. J. Weiler, Phys. Rev. D 101, 015012 (2020), arXiv:1905.00046 [hep-ph] .
  * Baxter  _et al._ (2020) D. Baxter, Y. Kahn, and G. Krnjaic, Phys. Rev. D 101, 076014 (2020), arXiv:1908.00012 [hep-ph] .
  * Liang  _et al._ (2020) Z.-L. Liang, L. Zhang, F. Zheng, and P. Zhang, Phys. Rev. D 102, 043007 (2020), arXiv:1912.13484 [cond-mat.mes-hall] .
  * Liu  _et al._ (2020) C. P. Liu, C.-P. Wu, H.-C. Chi, and J.-W. Chen, Phys. Rev. D 102, 121303 (2020), arXiv:2007.10965 [hep-ph] .
  * Kahn  _et al._ (2021) Y. Kahn, G. Krnjaic, and B. Mandava, Phys. Rev. Lett. 127, 081804 (2021), arXiv:2011.09477 [hep-ph] .
  * Flambaum  _et al._ (2020) V. V. Flambaum, L. Su, L. Wu, and B. Zhu, (2020), arXiv:2012.09751 [hep-ph] .
  * Bell  _et al._ (2021) N. F. Bell, J. B. Dent, B. Dutta, S. Ghosh, J. Kumar, and J. L. Newstead, Phys. Rev. D 104, 076013 (2021), arXiv:2103.05890 [hep-ph] .
  * Acevedo  _et al._ (2022) J. F. Acevedo, J. Bramante, and A. Goodman, Phys. Rev. D 105, 023012 (2022), arXiv:2108.10889 [hep-ph] .
  * Wang  _et al._ (2022) W. Wang, K.-Y. Wu, L. Wu, and B. Zhu, Nucl. Phys. B 983, 115907 (2022), arXiv:2112.06492 [hep-ph] .
  * Blanco  _et al._ (2022) C. Blanco, I. Harris, Y. Kahn, B. Lillard, and J. Pérez-Ríos, (2022), arXiv:2208.09002 [hep-ph] .
  * Liang  _et al._ (2021) Z.-L. Liang, C. Mo, F. Zheng, and P. Zhang, Phys. Rev. D 104, 056009 (2021), arXiv:2011.13352 [hep-ph] .
  * Knapen  _et al._ (2021) S. Knapen, J. Kozaczuk, and T. Lin, Phys. Rev. Lett. 127, 081805 (2021), arXiv:2011.09496 [hep-ph] .
  * Liang  _et al._ (2022) Z.-L. Liang, C. Mo, F. Zheng, and P. Zhang, (2022), arXiv:2205.03395 [hep-ph] .
  * Derenzo  _et al._ (2017) S. Derenzo, R. Essig, A. Massari, A. Soto, and T.-T. Yu, Phys. Rev. D96, 016026 (2017), arXiv:1607.01009 [hep-ph] .
  * Derenzo  _et al._ (2018) S. Derenzo, E. Bourret, S. Hanrahan, and G. Bizarri, J. Appl. Phys. 123, 114501 (2018), arXiv:1802.09171 [physics.ins-det] .
  * Ashcroft and Mermin (2020) N. Ashcroft and N. Mermin,  _Solid State Physics_ (Cengage., 2020).
  * Lindhard (1954) J. Lindhard, Dan. Vid. Selsk Mat.-Fys. Medd. 28, 8 (1954).
  * Adler (1962) S. L. Adler, Physical Review 126, 413 (1962).
  * Squires (1996) G. L. Squires,  _Introduction to the theory of thermal neutron scattering_ (Courier Corporation, 1996).
  * Trickle  _et al._ (2020) T. Trickle, Z. Zhang, K. M. Zurek, K. Inzani, and S. M. Griffin, JHEP 03, 036 (2020), arXiv:1910.08092 [hep-ph] .
  * Campbell-Deem  _et al._ (2022) B. Campbell-Deem, S. Knapen, T. Lin, and E. Villarama, Phys. Rev. D 106, 036019 (2022), arXiv:2205.02250 [hep-ph] .
  * Schober (2014) H. Schober, Journal of Neutron Research 17, 109 (2014).
  * Marshall  _et al._ (1971) W. Marshall, S. W. Lovesey, _et al._ ,  _Theory of thermal neutron scattering: the use of neutrons for the investigation of condensed matter_ (Clarendon Press, 1971).
  * Berghaus  _et al._ (2022) K. V. Berghaus, R. Essig, Y. Hochberg, Y. Shoji, and M. Sholapurkar, Phys. Rev. D 106, 023026 (2022), arXiv:2112.09702 [hep-ph] .
  * Nelin and Nilsson (1972) G. Nelin and G. Nilsson, Phys. Rev. B 5, 3151 (1972).
  * Kim  _et al._ (2020) D. S. Kim, O. Hellman, N. Shulumba, C. N. Saunders, J. Y. Y. Lin, H. L. Smith, J. E. Herriman, J. L. Niedziela, D. L. Abernathy, C. W. Li, and B. Fultz, Phys. Rev. B 102, 174311 (2020).
  * Petretto  _et al._ (2018) G. Petretto, S. Dwaraknath, H. Miranda, D. Winston, M. Giantomassi, M. van Setten, X. Gonze, K. Persson, G. Hautier, and G.-M. Rignanese, Scientific Data 5, 180065 (2018).
  * Gunn and Warner (1984) J. M. F. Gunn and M. Warner, Zeitschrift fur Physik B Condensed Matter 56, 13 (1984).
  * Kozaczuk and Lin (2020) J. Kozaczuk and T. Lin, Phys. Rev. D 101, 123012 (2020), arXiv:2003.12077 [hep-ph] .
  * Mortensen  _et al._ (2005) J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. B 71, 035109 (2005).
  * Enkovaara  _et al._ (2010) J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. Hansen, _et al._ , Journal of physics: Condensed matter 22, 253202 (2010).
  * Essig  _et al._ (2016) R. Essig, M. Fernandez-Serra, J. Mardon, A. Soto, T. Volansky, and T.-T. Yu, JHEP 05, 046 (2016), arXiv:1509.01598 [hep-ph] .
  * Angle  _et al._ (2011) J. Angle _et al._ (XENON10), Phys. Rev. Lett. 107, 051301 (2011), [Erratum: Phys.Rev.Lett. 110, 249901 (2013)], arXiv:1104.3088 [astro-ph.CO] .
  * Aprile  _et al._ (2019b) E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 251801 (2019b), arXiv:1907.11485 [hep-ex] .
  * Bringmann and Pospelov (2019) T. Bringmann and M. Pospelov, Phys. Rev. Lett. 122, 171801 (2019), arXiv:1810.10543 [hep-ph] .
  * (69) D. N. McKinsey et al. [TESSERACT Collaboration] (TESSERACT), ‘‘The TESSERACT Dark Matter Project, SNOWMASS LOI,’’ Availabe [Online].
  * Aguilar-Arevalo  _et al._ (2022) A. Aguilar-Arevalo _et al._ , (2022), arXiv:2202.10518 [astro-ph.IM] .
  * Jackson (2012) J. Jackson,  _Classical Electrodynamics_ (Wiley, 2012).
  * Hudson (2006) B. S. Hudson, Vibrational Spectroscopy 42, 25 (2006), a Collection of Papers Presented at the 3rd International Conference on Advanced Vibrational Spectroscopy (ICAVS-3), Delavan, WI, USA, 14-19 August 2005 - Part 1.
  * Abernathy  _et al._ (2012) D. L. Abernathy, M. B. Stone, M. Loguillo, M. Lucas, O. Delaire, X. Tang, J. Lin, and B. Fultz, Review of Scientific Instruments 83, 015114 (2012).