Measuring the Migdal effect in semiconductors for dark matter detection

The Migdal effect has received much attention from the dark matter direct detection community, in particular due to its power in setting limits on sub-GeV particle dark matter. Currently, there is no experimental confirmation of the Migdal effect through nuclear scattering using Standard Model probes. In this work, we extend existing calculations of the Migdal effect to the case of neutron-nucleus scattering, with a particular focus on neutron scattering angle distributions in silicon. We identify kinematic regimes wherein the assumptions present in current calculations of the Migdal effect hold for neutron scattering, and demonstrate that these include viable neutron calibration schemes. We then apply this framework to propose an experimental strategy to measure the Migdal effect in cryogenic silicon detectors using an upgrade to the NEXUS facility at Fermilab.

A proliferation of direct detection experiments searching for sub-GeV dark matter (DM) has been matched by a suite of theoretical work to better understand the kinematics of low-energy scattering in the regime where particle physics and condensed matter intersect [1]. This kinematic regime primarily differs from traditional WIMP scattering in that the energy and momentum transfers involved are comparable to the fundamental scales of the target (set by the gap energy and inverse atomic size, respectively), meaning that standard elastic scattering approximations [2] no longer hold. Indeed, the primary scattering channel of interest for sub-GeV DM searches has long been DM-electron scattering [3], which must account for both the inherent binding energy of the scattered electron and the band structure of the target. More recently, several theoretical advancements have uncovered yet another inelastic scattering channel of interest for sub-GeV DM, nuclear recoils that directly ionize the scattered atom, a process denoted the "Migdal effect" (ME).
The theoretical underpinnings of the ME go back to the early work of Arkady Migdal [4,5], who calculated the probability that a radioactive decay would directly ionize the daughter nucleus. Such ionization has been measured in radioactive decay, and is more commonly referred to as "electron shake-off" [6][7][8]. Though a handful of papers [9][10][11] pointed out the likely relevance of this effect for DM-nucleus scattering, progress on the ME accelerated after Ref. [12] derived the necessary electronic excitation probabilities relevant for DM experiments. A flurry of theoretical activity [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27] followed, expanding the theory of the ME for galactic DM scattering in iso-lated atom [18], molecular [27], and solid-state [20,21,26] targets, as well as for solar coherent elastic neutrinonucleus scattering [14]. The ME was also shown to dominate over another important inelastic channel, namely the bremsstrahlung process [14,28]. Several experimental collaborations have since used these theoretical results to set what are currently the strongest limits on DM-nuclear scattering below ∼1 GeV [29][30][31][32][33][34][35]. Out of all of this work from the DM community, only Refs. [36][37][38][39][40] have so far explicitly considered the ME for neutron scattering. A parallel effort in chemistry focused on neutron scattering in isolated atoms and molecules [41,42], in part to explain anomalously large neutron cross sections on hydrides [43][44][45]. Our work differs from these proposals by focusing on the angular distribution of neutrons scattered from solid-state targets. 1 While the existence of the ME is well founded, the magnitude of the effect must be measured to understand the expected DM signal in a direct detection experiment.
In this Letter, we highlight many of the subtle differences in the ME between the cases of sub-GeV DM and traditional neutron probes. These differences arise because sub-GeV DM is lighter than the neutron and thus carries less momentum than a neutron of the same kinetic energy. We carefully delineate the theoretical approximations made when calculating ME rates to define the regime where they continue to hold for neutron scattering, which differs considerably from the regime of validity for DM scattering depending on the neutron energy. We expand the framework of Ref. [38] to include the angular dependence of neutron scattering, as is used in standard neutron calibration experiments involving neutron detec-tor backing arrays. A key finding of this study is that, for an isotropic target in the limit of small momentum transferred to the electronic system (the "soft limit"), the angular distribution of the scattered neutron factorizes from the electronic matrix element, allowing for a direct calibration of this matrix element for DM scattering. We demonstrate that the electronic matrix element for the ME can be measured with a greatly reduced or even completely absent elastic scattering background by a judicious choice of the neutron beam energy and the neutron scattering angle. However, the expected rates we find for such a measurement are at the boundary of what is currently feasible with existing setups and techniques. We conclude that, in order to measure the ME with neutrons in semiconductors, a dedicated low-energy neutron calibration setup is required, and propose one such experiment using modifications to the existing NEXUS facility [47] at Fermi National Accelerator Laboratory (Fermilab).
The ME is defined as the ionization or excitation of an atomic electron accompanying the recoil of the atom's nucleus [5]. For sub-GeV DM, the ME greatly enhances the sensitivity of direct detection experiments to the DMnucleon cross section [15,16,48] because the electronic excitations are observable even when the nuclear recoil is below threshold. For both isolated atoms and semiconductors, under various sets of assumptions (which will be discussed further below and delineated in detail in Appendices B and C), the ME rate spectrum R M factorizes into a quasi-elastic nuclear recoil rate R el and an electronic excitation probability d P e /dω, such that Here, E r is the nuclear recoil energy, ω is the total energy deposited in the electronic system (excitation or ionization), and q = qq is the momentum transfer from the neutron probe to the target. For both classes of targets, the electronic spectrum scales as q 2 , and we have explicitly factored out this scaling. The goal of this Letter is to devise a scheme to measure d P e /dω in semiconductors. For isolated atoms, the electronic ionization spectrum is [12] d P e dω where m e and m N are the electron and nucleus mass, respectively, r e is the electron position operator, the sum runs over initial and final single-electron orbital quantum numbers, and the final state is a spherical wave with wavenumber k = 2m e (ω − |E b |) where E b is the binding energy of the initial state. Eq. (2) was derived within the context of the Born-Oppenheimer approximation [12], but in fact does not require this assumption and is correct to O(m e /m N ) 2 [1,41]. On the other hand, Eq. (2) does assume q m N /(m e a 0 ) 200 MeV m N 26 GeV where a 0 is the Bohr radius, since it was derived from the dipole approximation to the exponential exp i me m N q · r e . Note that the dependence on q in Eq. (2) drops out when summing over sphericallysymmetric filled electron shells, such that d P e /dω is isotropic.
In a solid-state system, the ME derivation must be modified because the constituent atoms are no longer free, which gives a characteristic energy scale ω ph for optical phonon excitation (≈ 10 − −100 meV in typical solid-state systems [1]) and also removes the constraint of exact momentum conservation for the nucleus because the atoms are no longer in momentum eigenstates. However, the form of Eq. (1) can be recovered under the following three assumptions: • Impulse approximation: q 2m N ω ph • Free-ion approximation: initial nucleus state is a zero-momentum plane wave • Soft limit: k q and q · k m N ω, where k is the momentum transferred to the electronic system.
The result for a solid-state system is [20] d P e dω sol.

=
4α where Z ion (k) is an effective momentum-dependent charge of the nucleus plus inner-shell electrons, α 1/137 is the fine structure constant, and W( k, ω) is the energy loss function (ELF) of the target which measures its response to charge perturbations. If the ELF is isotropic and only depends on the magnitude of k, the dependence onq drops out, as in the atomic case. We have verified that for all of the kinematic configurations we will consider, the impulse approximation is valid.
To demonstrate the main kinematic features of the ME, we model the valence shell of silicon with Eq. (3), using the isotropic GPAW ELF from DarkELF [49], which has a regime of validity ω 75 eV, k 22 keV. For larger ω, we model the inner-shell electrons with Eq. (2) [12]. This division is somewhat artificial, and we discuss its limitations in Appendix C. As we show in Appendices A and B, in the soft limit, we can convert the nuclear recoil spectrum into an angular spectrum: where E n is the kinetic energy of the incident neutron, θ n is the lab-frame angle of the scattered neutron,P is a kinematic prefactor containing all angular dependence, and we have expressed the spectrum as a differential probability P M of Migdal scattering per incident neutron (rather than a flux-dependent rate). Fig. 1 illustrates these kinematics and the experimental setup.
FIG. 1. A diagram of an ideal neutron scattering experiment with a backing array, which consists of a series of active (red) and passive (grayscale) elements. Neutrons are generated isotropically from a source placed inside of a shield with a small opening to collimate the beam. These neutrons then enter the vacuum chamber (often a dilution refrigerator) with energy En, and scatter with lab-frame angle θn into a circular backing array element after transferring Er of energy to the nuclear recoil and energy ω to electrons, which together are detected as ionization energy Eion. Unscattered neutrons, meanwhile, pass through a capture detector (e.g. 3 He counter) to help normalize the simulated beam flux before arresting in a beamstop.
In Eq. (4), we have explicitly noted the separation of the ionization (ME) probability and the kinematic prefactor containing the angular dependence, and absorbed the q 2 scaling from Eq. (1) into P such that it now has units of (events/neutron)×[eV] 2 . This unconventional choice of normalization allows us to group all the terms depending on the experimentally-controllable variables E n and cos θ n together in the explicit expression where m n is the neutron mass, µ ≡ m n m N /(m n +m N ) is the reduced mass, σ el is the elastic neutron cross section on a target material with density ρ T and thickness L in the beam direction, A N is the target's atomic mass number, N 0 is Avogadro's number, and In Appendix C, we demonstrate that the factorization of Eq. (4) that allows this separation does not hold outside of the soft limit for a semiconductor. We therefore note that calibrating the semiconductor ME outside of the soft limit, for example with high-energy (MeV-scale) neutrons, is fundamentally no longer probing the same regime as sub-GeV DM-nucleus scattering, where the soft limit approximations always hold. For the proposed calibrations discussed in the rest of this Letter, we will fall safely within the soft limit (see Appendix B), and thus the results of any such calibration are effectively a measurement of d P e /dω that can be directly translated to DM-nucleus scattering.
Since we do not measure ω directly, we change variables again to the observable E ion , the total amount of energy available as ionization, defined as where f n (E r ) denotes the ionization efficiency for elastic nuclear recoils as a function of the nuclear recoil energy E r (ω, θ n ). For the purposes of this work we consider the Sarkis model [50] as a best theoretical approximation for the ionization efficiency in the mostly uncalibrated regime of small E r [51]. In general, calibrations of the ME will be dependent on this ionization efficiency (quenching) model; however, we propose two specific calibration schemes in this Letter designed to minimize this dependence. A more complete treatment would also consider the systematic or theoretical fluctuations in f n (E r ), but this is outside the scope of this work.
To predict the number of electron-hole pairs n e as a function of E ion , we use the charge production model presented in Ref. [52]. This is a data-driven model of impact ionization in silicon that more accurately models the response for low n e than a model of Fano statistics alone (i.e. Ref. [53]). Ref. [52] provides a set of functions p ne (E ion ) for the probability of producing n e pairs for energy deposit E ion . Thus, to compute measured ionization rates as a function of angle, we integrate Eq. (4) against p ne to find the differential angular probability of Migdal events binned in n e , The inherent widths of the p ne leads to a smearing effect that can affect our signal (even before considering experimental factors, such as non-monochromaticity of the beam). We show two examples of observable spectra in Fig. 2 for different choices of E n and θ n , where we decompose the spectral contributions to the rate from elastic, valence band ME, and inner shell ME scatters.

FIG. 2.
Differential probability spectra dPn e /dθn (in units of events/neutron/degree of angular coverage) are shown per detectable charge quanta ne in the left (right) plot for an ideal 1 cm thick silicon detector in a En = 24 (2) keV monoenergetic neutron beam at a fixed scattering angle of θn = 72 (10) degrees, assuming the Sarkis ionization efficiency (quenching) model [50] and Ramanathan charge production model [52]. In both cases, we assume perfect backing detector with full azimuthal coverage. Left: for higher neutron energies and wide angles, the contribution from the inner shell [12] is distinct above the elastic peak. Right: for low neutron energies and shallow angles, the contribution from the valence band [49] separates from the elastic peak. Fig. 2 illustrates two possible strategies for calibrating the ME in the correct kinematic regime. The q 2 scaling of the Migdal probabilities translates to an enhancement approximately proportional to E n (1 − cos θ n ). Larger momentum transfers (which lead to larger nuclear recoil energies), achieved either by raising the neutron beam energy or by looking at a larger scattering angle, will therefore give a higher rate of Migdal events. However, to avoid the aforementioned difficulties of the inherent smearing due to the Fano statistics, it is important to keep the nuclear recoil energy scale small enough that the elastic spectrum does not smear too much into the Migdal tail. Thus, the first experimental strategy is to target setups that balance the q 2 rate enhancement with low recoil energy, in order to clearly isolate the highside Migdal rate tail. As can be seen in the left plot of Fig. 2, this strategy is particularly useful for calibrating Eq. (2), the ME contribution for inner shell electrons, but care must be taken not to increase the neutron energy and scattering angle outside of the kinematic regime of interest (see Appendix B).
The second strategy is to employ low-energy neutrons scattering at low angles such that the quenched nuclear recoils are too small to produce any secondary ionization, effectively eliminating the observable elastic contribution (the second term of Eq. (7)). This strategy is challenging in that it involves novel neutron source development, but is able to calibrate Eq. (3), the ME contribution from valence electrons independent of other contributions, as shown in the right plot of Fig. 2. Since sub-GeV DM will typically only produce single-to few-electron events, this setup more closely mimics what a sub-GeV DM sig-nal would look like in a single-electron threshold charge detector.
In a real experiment, no neutron beam will be perfectly monochromatic, such that contamination from higherenergy neutrons and gamma backgrounds must be taken carefully into account through validated simulations. There are a number of common neutron sources and methods that are implemented in the lab, each of which can be turned into a fairly monochromatic beam with careful application. These include deuterium-deuterium (D-D) and deuterium-tritium (D-T) generators [54], proton accelerators incident on a 7 Li [55] or 51 V [56] target, and photoneutron sources that exploit the 9 Be disintegration threshold of 1.67 MeV [57,58]. Each of these options comes with its own advantages and disadvantages, so we will emphasize that the fairly rare probability of a Migdal scatter, even in an ideal setup, necessitates (a) a low-background environment, as is typically achieved with significant overburden, thus complicating the use of proton accelerators, and (b) a high flux of low-energy neutrons, which can be achieved with either photoneutron sources or moderated D-D (or D-T) generators. Collimation in both of these cases is achieved through robust 4π shielding minus a small beam hole, typically around 1 cm in diameter (dependent on detector size).
Because of logistical challenges associated with having a sufficiently high-activity gamma source to produce a high flux of neutrons from a photoneutron source, we will focus the rest of the Letter on using a D-D neutron generator, as is employed by NEXUS. A D-D generator leverages fusion reactions to generate isotropic 2.5 MeV neutrons without any primary gamma backgrounds [54] (although secondary gammas will be produced by neutrons interacting in surrounding shielding materials). Using clever application of "filters," it is possible to prune, or even adjust, the beam energy spectrum by exploiting anti-resonances in the neutron scattering cross section [59]. Filters have the added advantage of removing unwanted secondary gamma backgrounds from neutron interactions in the shield materials and any primary xrays produced by the generator, which can be shielded by even a small amount of material. Of note, prominent anti-resonances in iron and scandium can be used to select 24 keV [60] and 2 keV [61] neutrons, respectively. A downside of using filters to select an optimal beam energy is a substantial reduction in neutron flux, thus requiring longer exposures, hotter sources, and lower ambient backgrounds. Another option for reducing the neutron energy is to employ neutron reflectors [62], but this would require more substantial modification to the NEXUS setup, and so is not the focus of this study. Lower-energy neutron beams also mandate progress in low-energy neutron backing detectors, which is an active area of study [63] but outside the scope of this Letter.
As a schematic setup, the NEXUS facility at Fermilab is designed to provide a D-D generator neutron beam incident on a 10 mK, single-electron resolution detector (e.g. SuperCDMS HVeV [64][65][66][67]) in a ∼100 cts/kg/day/keV radiation environment. Crucially, the chosen detector should be thin compared to the mean free path in silicon of ∼10 cm for a <50 keV neutron (for which the cross section is constant [68]), to ensure that neutrons scatter only once on average and thus that angular smearing from multiple scattering is not a concern. The NEXUS D-D generator (Adelphi model DD108) produces up to ∼10 9 neutrons/s isotropically, with a collimated ∼10 3 neutrons/s rate incident on a ∼1 cm 2 detector area. Filters should be able to modulate a higherenergy neutron source (such as the D-D generator) down to the respective anti-resonance with roughly 10% energy width, minimal higher-energy contamination, and a ∼10 3 reduction in overall flux. This means that, with minimal modification, NEXUS could produce a filtered keV-scale neutron beam with ∼1 neutron/s incident on a singleelectron threshold semiconductor detector with a custom backing array in a low-background environment; to the best of our knowledge, no other such facility currently exists.
In such a setup with a series of ∆θ n ≈ 10 • wide backing arrays at different angles (including at θ n ±5 • around the central angles shown in Fig. 2), one would expect to see only a handful of neutron-induced Migdal events with roughly one month of exposure in the case of the 24 keV beam. This should be sufficient to calibrate the normalization of the electronic matrix element in Eq. (2), but upgrades to increase the rate would be required to fully reconstruct it. Meanwhile, the 2 keV beam setup requires more exposure than is practical without a more substantial upgrade to NEXUS in order to calibrate even the normalization of the matrix element in Eq. (3). To ac-complish this, the rate could be increased with a hotter D-D (or D-T) neutron source or by deploying multiple silicon detectors in the beam, but each of these improvements comes with trade-offs and complications. One of the biggest challenges in any of these setups would be to sufficiently eliminate higher-energy neutron contamination in the beam from the energy region of interest, which can hopefully be achieved through careful angular tagging.
In this work we have extended previous calculations of the ME to study the angular distributions of the neutron-induced ME in silicon.
These results can also be applied to an atomic target calibration by using Eq. (2) for the valence shell as well as the inner shells. We have demonstrated that Migdal scatters leave a distinct pattern in ionization measurements at fixed angles, providing a clear experimental target for calibration studies. We further emphasize that inherent spreading in the energy resolution in silicon strongly motivates the use of lower-energy neutrons and angular selection for a clean measurement. Lower-energy neutrons are also better kinematically tuned to mimic sub-GeV DM scattering, allowing direct calibration of ME probabilities in the kinematic regime of interest. Practical applications of this work will need to account for detector-specific backgrounds and non-ideal beam effects in their design, as well as the backgrounds from inelastic nuclear scattering. This work lays out necessary steps toward the calibration of the ME with neutrons in silicon (and germanium), which will be crucial to validate both existing (e.g. EDELWEISS [30] CDEX [31], SuperCDMS [34], DAMIC at SNOLAB [69], and SENSEI [70]) and next-generation (e.g. SENSEI at SNOLAB, DAMIC-M [71,72], Oscura [73], and Super-CDMS SNOLAB [74]) limits on sub-GeV DM-nuclear scattering. In this Appendix, we derive the kinematics for inelastic 2-body scattering in the soft and free-ion limits, where the initial-state nucleus is at rest and the electron system takes energy ω but no momentum. In particular, the derivation in the lab frame is necessary to keep track of the scattered neutron angle; previous derivations from e.g. Ref. [12] are in the center-of-mass frame, with all angular dependence integrated out, whereas we want to preserve the angular dependence in the lab frame.
In the lab frame, and under the assumptions of the soft and free-ion limits, energy conservation gives where p i and p f are the initial and final momentum of the neutron, respectively, while momentum conservation gives since the momentum q transferred to the target goes entirely to the recoiling nucleus, which gets momentum q N = q. We can thus rewrite the energy conservation equation as Rearranging terms and plugging in for | p f |, we find Finally, we can simplify to Inverting this to solve for E r yields a quadratic equation with solution In principle, there are two roots, but for m n < m N (as is always the case in our setup) the second root is spurious. The choice of root is consistent with the limit ω → 0, where our expression reduces to standard textbook results for 2 → 2 elastic scattering of unequal masses (e.g. [75]): We note that many of these formulas simplify somewhat in the center-of-mass frame, where the neutron scattering angle θ n = θ n + O(m n /m N ) is almost identical to the lab frame for m n m N . However, for small angles the corrections are significant, so we work in the lab frame for consistency. We note that in this limit, these results are analogous to the more simplified center-of-mass Eq. (94) from Ref. [12].
Finally, in order to translate the standard result into angular coordinates, we must take the derivative of Eq. (A6) with respect to cos θ n , which yields In this Appendix, we adapt the formalism of Ref. [20] to derive the angular spectrum of the scattered neutron in the soft limit, restoring the dependence on the momentum k transferred to the electronic system, and keeping careful track of any assumptions or approximations made along the way. We begin with the general expression for the electronic energy spectrum in the soft limit, Eq. (A33) in Ref. [20]: where K is a reciprocal lattice vector, F is a form factor parametrizing the zero-point momentum spread of the initialstate nucleus, and is the dielectric function of the target. As in Appendix A above, the momentum transferred to the target q is equal to the momentum of the recoiling nucleus q N in the soft limit, so we use q instead of q N to maintain the distinction with the full calculation outside the soft limit in Appendix C below. In the case of DM scattering, we typically integrate over the unobserved DM momentum p f , but for neutron scattering we want to keep p f and integrate over the unobserved nuclear recoil momentum q. The prefactor also changes, with where v n is the initial neutron velocity and b n is the neutron scattering length in silicon. Note that the convention in the neutron scattering literature is typically to define b n such that the neutron mass rather than the reduced mass appears in Eq. (B2); for 28 Si the corresponding value of b n is 4.1 fm [68]. For simplicity of notation, we will often abbreviate k + K ≡ k and d 3 k K → d 3 k , since the lattice structure will not be essential to our arguments. We will also write W( k , ω) ≡ Im(− −1 ( k , ω)) for the ELF. Following Ref. [20], we make the free-ion approximation where F can be replaced by a momentum-conserving delta function, The form factor F 2 is a Gaussian with width q 0 = 2m N ω ph 56 keV in Si, where ω ph 60 meV is a typical optical phonon energy. Note that the impulse approximation already requires q q 0 , so the spread in F 2 will typically not induce a large deviation from exact momentum conservation when the impulse approximation is satisfied. The impulse approximation will be valid for the kinematic regime we consider (p i q 0 , θ n not too small), but may fail for small neutron energies and/or very forward scattering. However, see Refs. [26,76], which demonstrate that the impulse approximation may be extended below its nominal regime of validity and coincides almost exactly with a full treatment using the phonon density of states.
Performing the q integral in Eq. (B1) using the delta function amounts to the replacement q = p i − p f . This leaves The radial part of the p f integral can be performed using the energy-conserving delta function, for which the algebra is equivalent to the derivation in Appendix A. The azimuthal integral is trivial and gives a factor of 2π, leaving where The square root in the definition of p f can, in principle, restrict the range of scattering angles: The kinematic threshold where scattering is forbidden occurs when the right-hand side is greater than 1.
For ω E n and m n < m N , which will always be the case for the kinematics we consider, the right-hand side is negative and there is no angular restriction, but the p − f solution is negative and therefore spurious. We will thus relabel p + f → p f . There is a very narrow range of energies close to threshold, ω ∈ where both roots are allowed. This is an extremely fine-tuned kinematical region, with the difference between the lower and upper boundaries being ∆ω/E n = m 2 n m 2 n +m 2 N 10 −3 for Si, and thus it is outside the regime of relevance for these studies (both because our proposed neutron source does not have this precision on the initial energy, and because we are never considering order-1 fractions of the initial energy taken by the electrons). However, it may be relevant for neutron scattering on very light targets such as helium.
As a check on these results, consider the elastic limit ω → 0. The angular restriction from the square root is which is always satisfied for any θ as long as m N ≥ m n . The solution for p f becomes which recovers the classical elastic scattering results.
To simplify the dot products in Eq. (B5), note from the original form of the energy delta function that so where θ k is the angle between the momentum transfer p i − p f and the momentum in the electron system k . Combining everything, we now have At this point, we are able to perform the remaining angular integrals if we assume isotropy of the target, such that Z 2 ion and W depend only on | k |. This assumption does not hold exactly for any lattice structure, but is likely a reasonable approximation for the highly-symmetric diamond cubic crystal structure of silicon (or germanium). Assuming isotropy, the azimuthal integral trivially gives a factor of 2π, and treating cos θ k as the polar angle of the k integral, we pick up a factor of d cos θ k cos 2 θ k = 2 3 . Thus, we are left with Eq. (B13) shows that, under the assumptions of isotropy and the soft limit, the only kinematic dependence of the integrand is carried by ω, and thus the Migdal rate factorizes as claimed in the main text. Indeed, the integral in Eq. (B13) is proportional to the electronic spectrum, Eq.(3) (repeated here for convenience), where in the second equality we have used the assumed isotropy of the ELF to integrate over the angles. Restoring the prefactor C n , this gives the desired factorization, where d P e /dω depends only on ω and may be calculated using the DarkELF code package [49] independent of the neutron scattering experimental parameters. To convert the cross section to a probability per neutron P M , we use the relation where ρ T is the mass density of the target, L is the thickness of the target, N 0 is Avogadro's number, and A N ≈ m N /m n is the atomic mass number of the target. This expression is valid when L is much less than the neutron mean free path, which (as discussed in the main text) is necessary to prevent angular smearing from multiple scattering. Using the definition of the elastic cross section in terms of the scattering length, σ el ≡ 4πb 2 n , and substituting p 2 i = 2m n E n and Eq. (B6) for p f , gives Eq. (5) in the main text which is an explicit expression for the angular spectrum in terms of the experimental variables E n , ω, and cos θ n .
As a final check on our results, we recover the original form of the Migdal rate as an energy spectrum, Eq. (1) as follows. First, use Eq. (B10) to replace p 2 i − p 2 f − 2m n ω in the numerator with q 2 (m n /m N ). Converting the scattering probability to a rate R M by multiplying by ΦA, where Φ is neutron flux in neutrons/cm 2 /s and A is target area, gives where β is defined in Eq. (6). We can rewrite this using the Jacobian computed in Eq. (A8), where in the first equality we have replaced AL with V , the volume of the target. The second equality follows from the standard formulas for elastic nuclear recoil because we have defined σ el using the neutron scattering convention where the target is treated as infinitely heavy, µ → m n . Fig. 3 shows the spectrum resulting from Eq. (4) for a broader range of kinematics than shown in the main text. For nuclear quenching, we consider the Sarkis model [77] and the Lindhard model [78] as lower and upper bounds on the quenching factor, respectively. At the highest neutron energies and wide scattering angles, the soft-limit approximation clearly fails because q is too large, such that the scaling with q 2 is unphysical and the Migdal rate appears enhanced compared to the elastic scattering rate. We will discuss other failures of the soft limit in Appendix C below.
FIG. 3. Each panel shows the differential probability spectra dPn e /dθn (in units of events/neutron/degree of angular coverage) per detectable charge quanta ne for elastic (solid) and Migdal (dashed) scattering expected for the Lindhard [78] (blue), Sarkis [77] (red), and Sarkis 2022 [50] (orange) ionization efficiency models and Ramanathan charge production model [52] in a setup with an ideal neutron beam of energy En incident on a 1 cm thick silicon detector with perfect detection efficiency. The left (right) column shows measured spectra for wide-(low-) angle neutron scattering at θn = 72 (10) degrees and the rows show, from top to bottom, the spectra for incident monoenergetic neutrons of 2.5 MeV (as from an unmoderated D-D generator), 24 keV (as from an iron filter), and 2 keV (as from a Sc filter). In all cases, we assume a perfect backing detector with full azimuthal coverage. Note that the high-energy, wide-angle case (top left) is outside of the soft-limit regime where our derivations are valid, giving the nonphysical amplification shown. On the other hand, the low-energy, low-angle case (bottom right) clearly demonstrates the scenario wherein the ME can be probed and measured independent of charge yield model using a single-electron sensitive detector.
In this Appendix, we investigate how the kinematics of the Migdal effect change outside the soft limit. We start now with (A32) from Ref. [20], suitably modified for neutron scattering: As in Eq. (B2), we change the prefactor to C n appropriate for neutron scattering, abbreviate k ≡ k + K, and assume the free-ion approximation. Here, we write q N instead of q for the momentum transferred to the nucleus, because the electron system momentum k appears in the momentum delta function. To facilitate comparison to the soft-limit derivation, we can first rewrite the expression in parentheses as (C2)

Soft and low-momentum limits
The two components of the soft limit correspond to dropping k from two different parts of the integrand in Eq. (C1). Specifically, Assuming only condition (SA) and isotropy of the target, we can perform identical manipulations to those in Appendix B, and Eq. (C1) reduces to dσ dωd cos θ n = C n m n m N ω 4 where | p i − p f | = p 2 i + p 2 f − 2p i p f cos θ n and p f is given by p + f in Eq. (B6). We refer to the result obtained using only (SA) as the low-momentum limit. At this point it is clear that unlike in the case of the soft limit, the angular dependence of the scattered neutron does not factorize from the electronic spectrum, even in the limit of an isotropic material, because of the presence of the term coupling cos θ k and ω, which cannot be ignored without assumption (SB).
Since k is not observable and is, in fact, integrated over in the rate, a strict application of the soft limit effectively restricts the range of integration of k for fixed momentum transfer q and electronic energy ω. Since | q · k | = O(qk ) unless q and k are very nearly orthogonal (which, for an isotropic material, would suffer a 1/(4π) suppression in the angular part of the k integral), the soft limit is roughly equivalent to where we have replaced q with √ 2m N E r . To make contact with the kinematics discussed in the main text, we can write E r as a function of E n , θ n , and ω using Eq. (A6). In Fig. 4, we explicitly plot the soft limit condition (C4) for the kinematics considered in Fig. 3. The two branches of k max correspond to condition (SA) at small angles and (SB) at large angles, and the red and green horizontal lines correspond to the regimes of validity of two models for the ELF in silicon, an isotropic free-electron gas (see e.g. Ref [79]) and the GPAW model implemented in DarkELF, respectively. As noted in Ref. [80], the free-electron gas model with Fermi velocity v F = 8.6 × 10 −3 c and plasmon frequency ω p = 18.5 eV is a reasonable approximation for the measured ELF in silicon for k 10 keV and 5 eV ω 30 eV, capturing, in particular, the Fermi-broadened free-electron peak at k √ 2m e ω. Deviations from the soft limit are expected when the blue curve drops below the red and green lines, which occurs either for very small scattering angles or for wider angles with sufficiently large E n and small ω (top row and middle left).
We now argue that the low-momentum limit matches the full result without (SA) very closely for our entire parameter space, which is convenient since Eq. (C3) is easily amenable to numerical integration given a loss function W(k , ω) . The relevance of (SB) but not (SA) for the soft limit can already be seen from Fig. 4, where (SA) is only violated for θ n 1 • for all choices of E n and ω. To be somewhat more quantitative, we adopt the simple free-electron gas model of the ELF described above, which has a closed-form analytic expression [79]. However, it features an unphysical vanishing of the ELF at large k where core electron wave functions should have nonvanishing support. A proper treatment of the ELF in silicon would include the effects of core electrons through, for example, "all-electron reconstruction" [81], which would eliminate the need to use the isolated atom formalism to compute the Migdal spectrum at large ω. We leave this for future work. For the k -dependent ion charge, we use an atomic form factor model, where λ TF = v F /( √ 3ω p ) is the Thomas-Fermi screening length and Z 0 = 4 is the charge of the silicon ion excluding the valence shell. Fig. 5 shows the effect on the angular spectrum dσ/d cos θ n of the extra term that would vanish in the full soft limit (SB). We have deliberately chosen kinematics that maximize the effect of (SB). The soft limit is expected to fail when the largest value of k allowed by (SB) drops below the region where the ELF has large support. In the case of both the free-electron gas and GPAW ELFs, the ELF vanishes identically outside the regime of validity shown in Fig. 4 above, and thus deviations from the soft limit are largest when the bound from (SB) is close to the limit of validity of the model. A more detailed calculation that accounts for inner electron shells in the ELF would not feature a hard cutoff in k but rather something closer to a power-law falloff from momentum-space atomic orbitals, but the general phenomenon we illustrate here will still hold.
The blue curve in Fig. 5 shows the soft limit from Appendix B, the green curve shows the low-momentum limit from Eq. (C3), and the orange curve shows the full calculation of Eq. (C1), the details of which are given in Sec. C 2 below. The kinematics are chosen to match Fig. 4, middle left and top right. There are order-1 deviations from the soft limit at wide scattering angles, but the low-momentum limit only differs from the full calculation at the percent level for any scattering angle. Furthermore, the soft limit underestimates the full result, because when (SB) does not hold, the nucleus propagator is closer to on shell.

Full calculation outside the low-momentum limit
For completeness, we now evaluate the full expression for the Migdal angular spectrum without assumption (SA). Starting from Eq. (C1), we can immediately perform the q N integral with the momentum delta function by making the replacement q N = p i − p f − k . This leaves To evaluate the remaining delta function, we use the nonrelativistic dispersion relation for the neutron, so that the delta function enforces or equivalently