# Atomic binding corrections for high energy fixed target experiments

Ryan Plestid  Mark B. Wise  Walter Burke Institute for Theoretical Physics,
California Institute of Technology, Pasadena, CA 91125, USA

###### Abstract

High energy beams incident on a fixed target may scatter against atomic
electrons. To a first approximation, one can treat these electrons as at rest.
For precision experiments, however, it is important to be able to estimate the
size of, and when necessary calculate, sub-leading corrections. We discuss
atomic binding corrections to relativistic lepton-electron scattering. We
analyze hydrogen in detail, before generalizing our analysis to multi-electron
atoms. Using the virial theorem, and many-body sum rules, we find that the
corrections can be reduced to measured binding energies, and the expectation
value of a single one-body operator. We comment on the phenomenological impact
for neutrino flux normalization and an extraction of hadronic vacuum
polarization from elastic muon electron scattering at MUonE.

††preprint: CALT-TH/2024-010

##  I Introduction

Leptonic interactions provide a clean laboratory for precision physics. In
this work we consider high energy leptons (specifically neutrinos and muons)
scattering on atomic electrons in a fixed target. When considering fixed
target experiments the asymptotic in-state includes atomic electrons bound to
nuclei. For certain fixed target experiments, high levels of precision are
necessary. Precision goals at these experiments can range from modest to
extremely demanding. As an example of a modest goal, the future DUNE
experiment will aim for a percent-level determination of their neutrino flux
normalization using $\nu e\rightarrow\nu e$ scattering Hewes  _et al._ (2021).
A more stringent example is provided by the MUonE experiment, which aims to
extract hadronic vacuum polarization from $\mu e\rightarrow\mu e$ scattering
with an error budget at the level of 10 ppm Abbiendi (2022); Banerjee  _et
al._ (2020).

It is intuitively obvious that these systems are well approximated by treating
the atomic electrons as at rest, and working in the lab frame. For experiments
with demanding precision goals, approximating the electrons as free and at
rest may not be sufficient at which point binding corrections should be
included. At bare minimum a parametric estimate for the size of these
corrections should be established, and if necessary they should be computed. A
number of basic questions immediately arise: what controls the binding
expansion? how large are these corrections? how do they scale with atomic
number? how does one include relativistic corrections?

To answer these questions concretely, we focus on neutrino and muon scattering
on atomic electrons. For example, consider neutrinos scattering from a
hydrogen atom at rest (i.e., with zero momentum),

$\nu(\mathbf{k})+~\\!\\!{}^{1}{\rm H}({\bf\it
0})\rightarrow\nu(\mathbf{k}^{\prime})+e^{-}(\mathbf{p}^{\prime})+p^{+}(\mathbf{h}^{\prime})~{}.$
(1)

For a generic atom $A$ we consider $\nu A\rightarrow\nu e^{-}B^{+}$ and
perform a sum over all final states $B^{+}$. We focus on kinematics such that
$|\mathbf{p}^{\prime}|\sim|\mathbf{k}^{\prime}|\sim|\mathbf{k}|$ are large,
and $|\mathbf{h}^{\prime}|\sim\sqrt{2m_{e}\epsilon_{A}}$ is small, where
$\epsilon_{A}$ is the binding energy of the atom; for hydrogen
$\epsilon_{H}\simeq\alpha^{2}m_{e}/2$. Our goal is to provide a simple
universal formula for the binding corrections to the cross section.

A natural question is what controls the size of binding corrections
parametrically. For example when matrix elements depend on the Mandelstam
variable $s=(k+p)^{2}$ with $k_{\mu}=(\omega,0,0,\omega)$, considering an
electron at rest vs. in motion we find $s=m_{e}^{2}+2m_{e}\omega$ vs.
$s=m_{e}^{2}+2(E_{e}\omega-p_{e}\omega\cos\theta)$. This would suggest
corrections of the size $p_{e}/m_{e}\sim\sqrt{\epsilon_{A}/m_{e}}$.
Alternatively, one may guess that since energy conservation is modified by the
binding energy $\epsilon_{A}$, and the available energy is of order $\omega$
that corrections will be parametrically given by $\epsilon_{A}/\omega$ a
dramatically smaller correction. As we will see in what follows, the leading
order corrections are of order $\epsilon_{A}/m_{e}$.

We study these binding corrections in detail for $\nu e\rightarrow\nu e$
scattering from an argon atom, and for $\mu e\rightarrow\mu e$ scattering from
a carbon atom. We find that the parametric scaling mentioned above holds in
general, but that a accidental cancellation makes the binding correction to
$\nu e\rightarrow\nu e$ scattering very small. By way of contrast, for $\mu
e\rightarrow\mu e$ scattering the corrections are of their natural size. For
MUonE, we find that these corrections are roughly $5\times$ larger than their
stated error budget goal, and must be included in future analyses.

The rest of the paper is organized along the following lines. First in Section
II we define, without approximation, the kinematics and phase space for
scattering from hydrogen. Next in Section III we describe how, at lowest order
in the fine structure constant $\alpha$ (i.e. neglecting final state
interactions), the relevant matrix element in the case of hydrogen can be
related to well known free electron matrix elements. Next, in Section IV we
outline a consistent expansion scheme to compute binding corrections. In
Section V we discuss realistic multi-electron atoms and how to use sum rules
and the virial theorem to reduce the problem to the atomic matrix element of a
single one-body operator. Next in Section VI we discuss applications at high
energy fixed target experiments including DUNE and MUonE. Finally in Section
VII we summarize our findings and outline potentially interesting future
directions.

##  II Kinematics and phase space

Consider a particle scattering on hydrogen producing a proton and electron in
the final state. For definiteness, we will take this to be a neutrino. The
relevant matrix element is written in terms of the weak Hamiltonian density
and is given by,

$\begin{split}T&=\int\differential^{4}x\matrixelement{\nu
e^{-}p^{+}}{\mathcal{H}_{W}(x)}{~\\!\\!{}^{1}{\rm H}\nu}~{}\\\
&=(2\pi)^{4}\delta^{(4)}(\Sigma\tilde{p})\matrixelement{\nu
e^{-}p^{+}}{\mathcal{H}_{W}}{~\\!\\!{}^{1}{\rm H}\nu}~{}.\end{split}$ (2)

The tilded energy and momentum conservation reads

$\displaystyle\Sigma\tilde{E}$
$\displaystyle=m_{e}-\epsilon_{H}+\omega-\omega^{\prime}-E^{\prime}-T_{p}^{\prime}~{}$
(3) $\displaystyle\Sigma\tilde{\mathbf{p}}$
$\displaystyle=\mathbf{k}-\mathbf{k}^{\prime}-\mathbf{h}^{\prime}-\mathbf{p}^{\prime}~{},$
(4)

where $T_{p}^{\prime}$ is the proton’s kinetic energy, $\omega$ is the
incident neutrino’s energy, $\omega^{\prime}$ the outgoing neutrino’s energy,
$E^{\prime}$ the outgoing electron’s energy, and $\epsilon_{H}$ the binding
energy of hydrogen. We will focus on scattering from an unpolarized target. We
average over the orientation of all atoms, which is equivalent to averaging
over the quantum number $m_{J}$ (i.e., the magnetic quantum number for total
angular momentum). This gives the same spin-averaged matrix elements as would
be obtained by summing over free-particle spin and averaging over $m_{\ell}$
(i.e., the orbital quantum number) Bechler (1993). Then, suppressing all spin
indices for simplicity, the cross section is given by

$\sigma=\frac{1}{2\omega}\frac{1}{2M_{{\rm
H}}}\int\differential\Pi_{\nu}\differential\Pi_{e}\differential\Pi_{p}~{}(2\pi)^{4}\delta^{(4)}(\Sigma\tilde{p})~{}|\mathcal{M}|^{2}~{},$
(5)

where $T=(2\pi)^{4}\delta^{(4)}(\Sigma\tilde{p})\mathcal{M}$ and Lorentz
invariant phase space is written for a particle of species-$a$ as
$\differential\Pi_{a}=\differential^{3}p_{a}/[(2E_{a})(2\pi)^{3}]$; this
formula assumes a relativistic normalization of states. It is convenient to
eliminate the neutrino’s phase space using the 3-momentum delta function,

$\sigma=\frac{1}{2\omega}\frac{1}{2M_{{\rm
H}}}\int\differential\Pi_{p}~{}\int\differential\Pi_{e}\frac{1}{2\omega^{\prime}}~{}(2\pi)\delta(\Sigma\tilde{E})~{}|\mathcal{M}|^{2}~{}.$
(6)

The energy conserving delta function can be used to eliminate the angle of the
outgoing electron, $\cos\theta^{\prime}$, relative to the beam axis. The
explicit solution is given by

$\cos\theta^{\prime}=\frac{2\omega(E^{\prime}-m_{e}+|\mathbf{h}^{\prime}|\cos\theta+\epsilon_{H}-T_{p})-(E^{\prime}-m_{e}+\epsilon_{H}-T_{p})^{2}+\mathbf{h}^{\prime
2}-2|\mathbf{h}^{\prime}||\mathbf{p}^{\prime}|\sin\phi\sin\theta^{\prime}\sin\theta+\mathbf{p}^{\prime
2}}{2|\mathbf{p}^{\prime}|(\omega+|\mathbf{h}^{\prime}|\cos\theta)}~{},$ (7)

where $\theta$ and $\phi$ label the directions of the outgoing proton.
Including the relevant Jacobian and integrating over the azimuthal angle
(trivially by symmetry) we then arrive at

$\begin{split}\sigma=\frac{1}{2\omega}\frac{1}{2M_{{\rm
H}}}\int\differential\Pi_{p}&\quantity(1+\frac{|\mathbf{h}^{\prime}|\cos\theta}{\omega})^{-1}\\\
&\hskip 43.36464pt\int\frac{\differential
E^{\prime}}{8\pi\omega}\quantity(1+\frac{T_{p}}{\omega^{\prime}})~{}|\mathcal{M}|^{2}~{}.\end{split}$
(8)

The limits of integration run from $E_{e}=m_{e}$ (the threshold for
ionization) to $E_{e}^{\rm max}$ which is given by solving

$E^{\prime}+\sqrt{\mathbf{q}^{2}+2\mathbf{p}^{\prime}\cdot\mathbf{q}+\mathbf{p}^{\prime
2}}+T_{p}^{\prime}=m_{e}-\epsilon_{H}+\omega~{}.$ (9)

where $\mathbf{q}=\mathbf{k}-\mathbf{h}^{\prime}$. At this stage all of our
expressions are exact, and we now begin to discuss sensible approximations.

We drop all terms suppressed by
$\lambda\sim\epsilon_{H}/m_{p},\epsilon_{H}/\omega$ (for example the proton’s
recoil energy).111Terms that go like
$\mathbf{h}^{\prime}/\omega\sim\lambda^{1/2}$ vanish upon integration over the
final-state proton’s phase space for a rotationally invariant state. Then,

$\begin{split}\sigma=\frac{1}{2\omega}\frac{1}{2M_{{\rm
H}}}\int\differential\Pi_{p}\int\frac{\differential
E^{\prime}}{8\pi\omega}~{}|\mathcal{M}|^{2}+O\quantity(\lambda)~{},\end{split}$
(10)

with

$\begin{split}\cos\theta^{\prime}&=\frac{2\omega(E^{\prime}-m_{e})-(E^{\prime}-m_{e})^{2}+\mathbf{p}^{\prime
2}}{2|\mathbf{p}^{\prime}|\omega}+O\quantity(\lambda)~{},\end{split}$ (11)

and

$E_{\rm
max}=\frac{1}{2}\left(\frac{m_{e}^{2}}{2\omega+m_{e}}+2\omega+m_{e}\right)~{}+O\quantity(\lambda)~{}.$
(12)

##  III Reduction to plane waves

The matrix element of the weak Hamiltonian in Eq. 2 can be written in terms of
plane-wave (i.e., free-electron) states. For the outgoing scattering state,
this is trivial at lowest order in $\alpha$ since $\ket{e^{-}p^{+}}_{\rm
out}\simeq\ket{e^{-}p^{+}}_{\rm PW}$ with the subscript denoting a plane-wave
state. Final state interactions can be added perturbatively order-by-order in
$\alpha$ but are not included in what we present below. The bound states
demand the inclusion of the Coulomb interaction at zeroth order, but may still
be expanded in plane wave states. At leading order in $\alpha$ a hydrogen atom
at rest can be written as222The state $\ket{e^{-}(\mathbf{p})}$ has a
polarization spinor $u(\mathbf{p})$. This accounts for the leading
relativistic correction to $\psi(\mathbf{p})$ (see Appendix C for a
discussion).

$\ket{~\\!\\!{}^{1}{\rm
H}}\simeq\sqrt{2M_{H}}\int\frac{\differential^{3}p}{(2\pi)^{3}}\frac{\tilde{\psi}(\mathbf{p})}{\sqrt{2E_{e}}\sqrt{2E_{p}}}~{}\ket{p^{+}e^{-}}_{{\rm
PW}}~{},$ (13)

with the proton carrying momentum $-\mathbf{p}$ and the electron carrying
momentum $\mathbf{p}$. We suppress spin indices for simplicity anticipating
spin averaging in the initial state. The function $\tilde{\psi}(\mathbf{p})$
is the Fourier transform of the coordinate space wavefunction
$\psi(\mathbf{x})$. The factors in square roots are necessary when working
with relativistically normalized states. The energies that appear are those of
free particles, $E_{e}=\sqrt{\mathbf{p}^{2}+m_{e}^{2}}$ and
$E_{p}=\sqrt{\mathbf{p}^{2}+m_{p}^{2}}$. One can easily check that provided
$\int\differential^{3}p/(2\pi)^{3}|\tilde{\psi}(\mathbf{p})|^{2}=1$ and the
electron and proton states are normalized relativistically, then it follows
the hydrogen atom has the correct normalization also.

Using
$\innerproduct{p^{+}(\mathbf{h}^{\prime})}{p^{+}(-\mathbf{p})}=2E_{p}(\mathbf{p})(2\pi)^{3}\delta^{(3)}(\mathbf{p}+\mathbf{h}^{\prime})$,
we may therefore write, to leading order in $\alpha$,

$\begin{split}\mathcal{M}&=\matrixelement{\nu
e^{-}p^{+}}{\mathcal{H}_{W}}{~\\!\\!{}^{1}{\rm H}\nu}\\\
&\simeq\sqrt{\frac{2M_{H}2E_{p}}{2E_{e}}}\tilde{\psi}(\mathbf{p})~{}{\sf
M}(\mathbf{p}^{\prime},\mathbf{k}^{\prime},\mathbf{p},\mathbf{k})\end{split}$
(14)

where we have introduced the matrix element for $2\rightarrow 2$ scattering of
an incident neutrino from a free electron with momentum $\mathbf{p}$,

${\sf
M}(\mathbf{p}^{\prime},\mathbf{k}^{\prime},\mathbf{p},\mathbf{k})=\bra{\nu(\mathbf{k}^{\prime})e^{-}(\mathbf{p}^{\prime})}{\mathcal{H}_{W}}\ket{\nu(\mathbf{k})e^{-}(\mathbf{p})}_{\rm
PW}~{}.$ (15)

Note that this matrix element is a function of the three momenta only and
carries no ambiguity.

Writing the cross section in terms of the free-electron matrix element we then
find,

$\begin{split}\sigma=\int\frac{\differential^{3}p}{(2\pi)^{3}}|\tilde{\psi}(\mathbf{p})|^{2}~{}\frac{1}{2E}\frac{1}{2\omega}\int\frac{\differential
E^{\prime}}{8\pi\omega}~{}\frac{1}{2}\sum_{\rm spins}|{\sf
M}|^{2}~{},\end{split}$ (16)

where spin sums and averages which were previously left implicit have been
made explicit. The factor of $1/2$ in this formula results from averaging over
the initial electron spin and assumes the incident neutrino is polarized (as
it must be).

Let us now evaluate the free-electron matrix element for the example at hand
of $\nu e\rightarrow\nu e$ scattering using

$\begin{split}{\sf
M}=-\sqrt{8}G_{F}&\quantity[\bar{u}(\mathbf{k}^{\prime})\gamma^{\mu}P_{L}u(\mathbf{k})]\quantity[\bar{u}(\mathbf{p}^{\prime})\gamma^{\mu}(c_{L}P_{L}+c_{R}P_{R})u(\mathbf{p})]~{},\end{split}$
(17)

$P_{L,R}=(1\mp\gamma_{5})/2$ are chiral-projectors. Summing over spins and
performing traces over Dirac matrices we obtain

$\sum_{\rm spins}|{\sf M}|^{2}=128G_{F}^{2}\times\big{[}c_{L}^{2}(p\cdot
k)(p^{\prime}\cdot k^{\prime})+c_{R}^{2}(p\cdot k^{\prime})(k\cdot
p^{\prime})-4c_{L}c_{R}m_{e}^{2}(k\cdot k^{\prime})\big{]}~{}.$ (18)

We have used the notation $k\cdot
p=k^{\mu}p_{\mu}=\omega\sqrt{\mathbf{p}^{2}+m_{e}^{2}}-\mathbf{k}\cdot\mathbf{p}$;
these invariants are only functions of three-momenta.

At this point it is important to emphasize that the kinematic relations for
our process differ from those of $2\rightarrow 2$ electron scattering. In
particular when using energy and momentum conservation to relate
$(p^{\prime}\cdot k^{\prime})$ to $(p\cdot k)$ we must account for the binding
energy. Let us define $\langle V\rangle\equiv
m_{e}-\epsilon_{H}-E(\mathbf{p})$, and introduce the four-velocity of the
hydrogen atom $u_{\mu}=(1,0,0,0)$ in the lab frame. We then have

$p_{\mu}+k_{\mu}+\langle V\rangle
u_{\mu}=p^{\prime}_{\mu}+k^{\prime}_{\mu}~{}.$ (19)

This allows us to trade different Lorentz invariants for one another,

$\displaystyle\begin{split}p^{\prime}\cdot k^{\prime}&=p\cdot k+\langle
V\rangle(p+k)\cdot u+\tfrac{1}{2}\langle V\rangle^{2}\\\ &\simeq p\cdot
k+\langle V\rangle\omega~{},\end{split}$ (20)
$\displaystyle\begin{split}p^{\prime}\cdot k&=p\cdot k^{\prime}+\langle
V\rangle(k^{\prime}-p)\cdot u-\tfrac{1}{2}\langle V\rangle^{2}~{}\\\ &\simeq
p\cdot k^{\prime}+\langle V\rangle\omega^{\prime}~{}.\end{split}$ (21)
$\displaystyle\begin{split}k\cdot k^{\prime}&=p\cdot
p^{\prime}+m_{\ell}^{2}-m_{e}^{2}+\langle V\rangle(p^{\prime}-p)\cdot
u-\tfrac{1}{2}\langle V\rangle^{2}~{}\\\ &\simeq p\cdot
p^{\prime}+m_{\ell}^{2}+\langle V\rangle E^{\prime}~{},\end{split}$ (22)

where $m_{\ell}=0$ for $\nu e\rightarrow\nu e$ scattering, but will be equal
to the muon mass when we consider $\mu e\rightarrow\mu e$ scattering.

##  IV Weak binding expansion

In the limit of weak binding we may estimate
$|\mathbf{p}|\sim\sqrt{2\epsilon_{H}m_{e}}$ and expand in
$\epsilon_{H}/m_{e}$. We need only the first terms in $\mathbf{p}^{2}$ and
$\langle V\rangle$ separately. Setting $m_{\ell}^{2}=0$,

$\displaystyle(p\cdot k)(p^{\prime}\cdot k^{\prime})$
$\displaystyle\simeq(p\cdot k)(p\cdot k+\omega\langle V\rangle)$ (23)
$\displaystyle\simeq(p\cdot
k)_{0}^{2}\quantity[1+\frac{\mathbf{p}^{2}}{m_{e}^{2}}(1+\cos^{2}\theta_{e\nu})+\frac{\langle
V\rangle}{m_{e}}]~{},$ $\displaystyle(p\cdot k^{\prime})(p^{\prime}\cdot k)$
$\displaystyle\simeq(p\cdot k^{\prime})(p\cdot
k^{\prime}+\omega^{\prime}\langle V\rangle)$ (24) $\displaystyle\simeq(p\cdot
k^{\prime})_{0}^{2}\quantity[1+\frac{\mathbf{p}^{2}}{m_{e}^{2}}(1+\cos^{2}\theta_{e\nu^{\prime}})+\frac{\langle
V\rangle}{m_{e}}]~{},$ $\displaystyle k\cdot k^{\prime}$ $\displaystyle\simeq
p\cdot p^{\prime}+\omega^{\prime}\langle V\rangle$ (25)
$\displaystyle\simeq(p\cdot
p^{\prime})_{0}\quantity[1+\frac{\mathbf{p}^{2}}{2m_{e}^{2}}+\frac{\langle
V\rangle}{m_{e}}]~{},$

where we have dropped terms linear in $\mathbf{p}$ that vanish upon averaging
over the electron’s wavefunction. The outgoing electron and neutrino are
nearly aligned with the beam axis, with angles set by
$1/\gamma\sim\sqrt{m_{e}/\omega}$ where $\gamma$ is the boost to the center of
mass frame. We may therefore set
$\cos\theta_{e\nu^{\prime}}\simeq\cos\theta_{ee^{\prime}}\simeq\cos\theta_{e\nu}\equiv\cos\theta$.333If
one worked with e.g., $k\cdot k^{\prime}$ rather than $p\cdot p^{\prime}$,
then large cancellations occur
$\omega\omega^{\prime}-\mathbf{k}\cdot\mathbf{k}^{\prime}\sim O(m_{e}\omega)$.
One cannot set
$\cos\theta_{e\nu^{\prime}}\simeq\cos\theta_{ee^{\prime}}\simeq\cos\theta_{e\nu}$
and must retain higher order terms.

The subscript zero denotes the Lorentz invariant evaluated for an elecron at
rest e.g., $(p\cdot k)^{2}_{0}=(m_{e}^{2}+2m_{e}\omega)^{2}$. We observe that
the terms multiplying $c_{L}^{2}$ and $c_{R}^{2}$ in Eq. 18 will be multiplied
by a common factor, whereas the term proportional to $c_{L}c_{R}$ is not.
Nevertheless this cross term is lepton-mass suppressed ($m_{e}/\omega\sim
10^{-3}$ for $\omega\sim 1~{}{\rm GeV}$), and we can neglect it when computing
the corrections.

The cross section can now be written in the simple form

$\begin{split}\sigma\simeq&\int\frac{\differential^{3}p}{(2\pi)^{3}}|\tilde{\psi}(\mathbf{p})|^{2}\quantity(1+\frac{1}{3}\frac{\mathbf{p}^{2}}{m_{e}^{2}}-\frac{\epsilon_{H}}{m_{e}})\\\
&\hskip 108.405pt\times\frac{1}{2m}\frac{1}{2\omega}\int\frac{\differential
E^{\prime}}{8\pi\omega}~{}\frac{1}{2}\sum_{\rm spins}|{\sf
M}|_{0}^{2}~{}.\end{split}$ (26)

where we have expanded $1/E=1/m(1-\mathbf{p}^{2}/2m_{e})$. We recognize the
second line as the cross section for scattering with an electron at rest. We
have exchanged $\langle V\rangle$ for $\epsilon_{H}$ using
$\mathbf{p}^{2}/2m_{e}+\langle V\rangle=-\epsilon_{H}$. We define the binding
correction $\differential\sigma_{\mathcal{B}}$ via
$\differential\sigma\simeq\differential\sigma_{0}+\differential\sigma_{\mathcal{B}}$.
We therefore find at the level of the differential cross section, the binding
correction can be written as

$\begin{split}\differential\sigma_{\mathcal{B}}&=\frac{1}{m_{e}}\frac{\matrixelement{~\\!\\!{}^{1}{\rm
H}}{\tfrac{2}{3}\hat{T}+\hat{H}}{~\\!\\!{}^{1}{\rm
H}}}{\innerproduct{~\\!\\!{}^{1}{\rm H}}{~\\!\\!{}^{1}{\rm
H}}}\times\differential\sigma_{0}\\\
&=\frac{1}{m_{e}}\langle\tfrac{2}{3}\hat{T}+\hat{H}\rangle\differential\sigma_{0}~{},\end{split}$
(27)

where $\langle\ldots\rangle$ means expectation value divided by the states
norm. This simple factorized form is a consequence of the symmetric nature of
the leading terms proportional to $c_{L}^{2}$ and $c_{R}^{2}$ in the tree-
level matrix element. Note that $\langle\hat{T}\rangle$ is related to
$\langle\hat{H}\rangle$ by the virial theorem, and so the entire correction is
fixed in terms of hydrogen’s binding energy.

##  V Beyond hydrogen

When considering more complicated atoms the analysis above must be
generalized. We are in fact interested in the sum over all final atomic
states. Let us consider a reaction

$\nu(\mathbf{k})+A({\bf\it
0})\rightarrow\nu(\mathbf{k}^{\prime})+e^{-}(\mathbf{p}^{\prime})+B^{+}(\mathbf{h}^{\prime})~{},$
(28)

where $B$ is an atomic system with one fewer electron than $A$, with center of
mass momentum $\mathbf{h}^{\prime}$, but otherwise unconstrained. We can
separate out the final state phase space of $B^{+}$ into its center of mass
motion and “the rest” $\sum_{B}$. We therefore have

$\begin{split}\sigma=\frac{1}{2M_{A}}\frac{1}{2\omega}\sum_{B}\int\differential\Pi_{B}&\int\differential\Pi_{\nu}\differential\Pi_{e}~{}\\\
&~{}~{}(2\pi)^{4}\delta^{(4)}(\Sigma\tilde{p})~{}|\mathcal{M}_{B}|^{2}~{}.\end{split}$
(29)

The energy conserving delta function now carries a dependence on the final
state $B$

$\begin{split}\Sigma\tilde{E}&=E_{A}+\omega-E^{\prime}-\omega^{\prime}-E_{B}~{}.\end{split}$
(30)

As above, we neglect the recoil energy of the system $B$

$E_{A}-E_{B}\simeq(m_{e}-\epsilon_{A})+\epsilon_{B}~{},$ (31)

where $\epsilon_{A}$ and $\epsilon_{B}$ are the atomic binding energies of $A$
and $B^{+}$ respectively. Introducing

${\cal E}=m_{e}+\omega-(E^{\prime}+\omega^{\prime})~{},$ (32)

(the energy conservation condition for an electron at rest) we have

$\Sigma\tilde{E}={\cal E}-(\epsilon_{A}-\epsilon_{B})~{}.$ (33)

Using the mode expansion of the field $\psi_{e}$, the matrix element can be
written in the form

$\begin{split}\matrixelement{e^{-}B^{+}}{\bar{\psi}_{e}\Gamma\psi_{e}}{A}=\int\frac{\differential^{3}q}{(2\pi)^{3}\sqrt{2E_{q}}}\bar{u}&(\mathbf{p}^{\prime})\Gamma
u(\mathbf{q})\\\
&\times\matrixelement{B^{+}}{\hat{a}_{\mathbf{q}}}{A}~{}.\end{split}$ (34)

with $\Gamma$ an appropriate Dirac bilinear. Notice that
$\matrixelement{B^{+}(\mathbf{h}^{\prime})}{\hat{a}_{\mathbf{q}}}{A({\bf\it
0})}\propto(2\pi)^{3}\delta^{(3)}(\mathbf{q}-\mathbf{h}^{\prime})$. We can
therefore replace $u(\mathbf{q)}\leftrightarrow u(\mathbf{h}^{\prime})$. It is
convenient to insert
$\int\differential\epsilon~{}\delta(\epsilon-\epsilon_{A}+\epsilon_{B})$ and
such that $\sum_{B}$ can be pulled through the other integrals,

$\begin{split}\sigma=\int\differential\epsilon&\int\frac{1}{2\omega}\int\frac{\differential\Pi_{e}}{2\omega^{\prime}}(2\pi)\delta({\cal
E}-\epsilon)\\\
&\sum_{B}\int\differential\Pi_{B}\int\frac{\differential^{3}q}{(2\pi)^{3}\sqrt{2E_{q}}}\frac{\differential^{3}p}{(2\pi)^{3}\sqrt{2E_{p}}}\delta(\epsilon_{A}-\epsilon_{B}-\epsilon)~{}\matrixelement{A}{\hat{a}_{\mathbf{q}}^{\dagger}}{B^{+}}\matrixelement{B^{+}}{\hat{a}_{\mathbf{p}}}{A}~{}\frac{1}{2}\sum_{\rm
spins}|{\sf
M}(\mathbf{p},\mathbf{k},\mathbf{p}^{\prime},\mathbf{k}^{\prime})|^{2}~{},\end{split}$
(35)

where the free-electron matrix element, ${\sf M}$ is defined above and we have
anticipated that the atomic matrix elements are spin diagonal after averaging
over initial spins of $A$ and summing over the finals spins of $B$. This may
be conveniently re-written in terms of the spectral function using
$\sumint_{B}=\sum_{B}\int\differential\Pi_{B}$,

$\displaystyle S_{A}(\epsilon,\mathbf{p})=$
$\displaystyle\int\frac{\differential^{3}q}{(2\pi)^{3}}\sumint_{B}\matrixelement{A}{\hat{a}^{\dagger}_{\mathbf{p}}}{B}\delta(\epsilon_{A}-\epsilon_{B}-\epsilon)\matrixelement{B}{\hat{a}_{\mathbf{q}}}{A}$
$\displaystyle=$
$\displaystyle\int\frac{\differential^{3}q}{(2\pi)^{3}}\sumint_{B}\matrixelement{A}{\hat{a}^{\dagger}_{\mathbf{p}}\delta(\epsilon_{A}-\hat{H}-\epsilon)}{B}\matrixelement{B}{\hat{a}_{\mathbf{q}}}{A}$
$\displaystyle=$
$\displaystyle\int\frac{\differential^{3}q}{(2\pi)^{3}}\matrixelement{A}{\hat{a}^{\dagger}_{\mathbf{p}}\delta(\epsilon_{A}-\hat{H}-\epsilon)\hat{a}_{\mathbf{q}}}{A}~{}.$
(36)

The Hamiltonian $\hat{H}$ is defined with all rest masses of elementary
particles subtracted. It may be written as
$\hat{H}=\hat{T}+\hat{V}_{1}+\hat{V}_{2}$ i.e., as a kinetic energy, electron-
nucleus potential, and inter-electron potential contribution. We then have

$\displaystyle\sigma$
$\displaystyle=\int\differential\epsilon\int\frac{\differential^{3}p}{(2\pi)^{3}}S_{A}(\epsilon,\mathbf{p})$
(37) $\displaystyle\hskip
65.04034pt\frac{1}{2\omega}\frac{1}{2E}\int\differential\Pi_{e}\frac{1}{2\omega^{\prime}}\delta({\cal
E}-\epsilon)\frac{1}{2}\sum_{\rm spins}|{\sf M}|^{2}.$

The energy and momentum conservation conditions now read

$E+\omega+\tilde{V}=E^{\prime}+\omega^{\prime}~{},$ (38)

with $\tilde{V}\simeq-\epsilon-\mathbf{p}^{2}/2m_{e}$. We then find for the
cross section

$\begin{split}\sigma=&\int\differential\epsilon\int\frac{\differential^{3}p}{(2\pi)^{3}}S_{A}(\epsilon,\mathbf{p})~{}\quantity(1+\frac{1}{3}\frac{\mathbf{p}^{2}}{m_{e}^{2}}-\frac{\epsilon}{m_{e}})\\\
&\hskip 108.405pt\times\frac{1}{2m}\frac{1}{2\omega}\int\frac{\differential
E^{\prime}}{8\pi\omega}~{}\frac{1}{2}\sum_{\rm spins}|{\sf
M}|_{0}^{2}~{},\end{split}$ (39)

When computing the corrections from $\mathbf{p}^{2}$ and $\epsilon$, it is
sufficient to use a non-relativistic spectral function which satisfies energy
weighted sum rules Polls  _et al._ (1994) (the latter of which is due to
Koltun Koltun (1972))

$\displaystyle\int\differential\epsilon\int\\!\frac{\differential^{3}p}{(2\pi)^{3}}\frac{\mathbf{p}^{2}}{2m_{e}}S_{A}(\epsilon,\mathbf{p})=\langle\hat{T}\rangle_{A}=\epsilon_{A}~{},$
(40)
$\displaystyle\int\differential\epsilon\\!\int\\!\frac{\differential^{3}p}{(2\pi)^{3}}(-\epsilon)~{}S_{A}(\epsilon,\mathbf{p})=\langle\hat{T}\rangle_{A}+\langle\hat{V}_{1}\rangle_{A}+2\langle\hat{V}_{2}\rangle_{A}$
$\displaystyle\hskip
216.81pt=\quantity(-3\epsilon_{A}-\langle\hat{V}_{1}\rangle_{A})~{}.$ (41)

To relate the matrix elements to binding energies we have used the virial
theorem: $\langle T\rangle_{A}=\epsilon_{A}$ and
$\langle\hat{V}_{1}+{V}_{2}\rangle_{A}=-2\epsilon_{A}$. For a sketch
derivation of the sum rules see Appendix A. We then find for the binding
correction

$\sigma_{\mathcal{B}}=\frac{1}{m_{e}}\frac{\matrixelement{A}{\tfrac{5}{3}\hat{T}+\hat{V}_{1}+2\hat{V}_{2}}{A}}{Z_{A}\innerproduct{A}{A}}\sigma_{0}~{}.$
(42)

Therefore, at the level of differential cross section for $\nu e\rightarrow\nu
e$ scattering we have

$\begin{split}\differential\sigma_{\mathcal{B}}&=\frac{1}{Z_{A}m_{e}}\quantity(-\frac{7}{3}\epsilon_{A}-\langle\hat{V}_{1}\rangle_{A})\differential\sigma_{0}~{}.\end{split}$
(43)

The relevant atomic input for the leading order binding corrections is
therefore the measurable binding energy, and a single one-body matrix element,
$\matrixelement{A}{\hat{V}_{1}}{A}/\innerproduct{A}{A}$. To evaluate this
matrix element for an atom with $Z_{A}$ electrons the radial density field of
electrons $n_{A}(\mathbf{r})=\langle\hat{n}(\mathbf{r})\rangle_{A}$ is
sufficient

$\langle\hat{V}_{1}\rangle_{A}=Z_{A}\times\int\differential^{3}r~{}~{}n_{A}(\mathbf{r})\quantity(\frac{-Z_{A}\alpha}{|\mathbf{r}|})~{},$
(44)

where $\int\differential^{3}r~{}n_{A}(\mathbf{r})=1$. This is particularly
interesting because $n_{A}(\mathbf{r})$ is amenable to calculation, and is in
general easier to predict than e.g., the full many-body wavefunction.

To estimate parametric scaling with $Z_{A}$ we can make use of the Thomas-
Fermi model (see Appendix B), which should be accurate at the 10% level or
better.444More accurate many-body calculations and the Thomas-Fermi model
differ by $\sim 10\%$ at the level of $n_{A}(\mathbf{r})$, but these
differences tend to largely cancel for integrated quantities such as
$\langle\hat{V}_{1}\rangle_{A}$. One finds,

$\langle\hat{V}_{1}\rangle_{A}\approx-Z_{A}^{7/3}\times(30~{}{\rm eV})~{}.$
(45)

In this expression one factor of $Z_{A}$ comes from the total number of
electrons in the atom. Another factor of $Z_{A}$ comes from the nuclear
Coulomb potential. The final factor of $Z_{A}^{1/3}$ comes from the Thomas-
Fermi length scale $b_{0}\approx 0.88Z_{A}^{-1/3}a_{0}$ with $a_{0}=1/(\alpha
m_{e})$ the Bohr radius.

##  VI Applications

In this section we discuss some simple applications of our methods to problems
of physical interest. We focus DUNE and MUonE specifically.

###  VI.1 Neutrino flux measurements

An immediate application of our results is the measurement of neutrino flux
normalization at future high intensity neutrino experiments. For concreteness
we will focus on scattering event rates at the DUNE near detector with a
liquid argon time projection chamber (LArTPC). This is interesting both
because of DUNE’s anticipated high statistical sample of $\nu e\rightarrow\nu
e$ scattering, but also because argon is a noble element. Therefore, the
treatment used above in terms of a spherically symmetric isolated atom (i.e.,
without complications due to molecular physics) applies immediately.

Recent estimates suggest that a measurement of the $\nu e$ elastic scattering
rate in a 30 ton LArTPC in the near detector hall can supply a measurement of
the neutrino flux normalization with a $2\%$ uncertainty after five years of
operation Marshall  _et al._ (2020). The authors of Ref. Marshall  _et al._
(2020) assume a $\nu e$ elastic scattering rate of $\sim 120\times
30=3.6\times 10^{3}$ events per year for a 30 ton liquid argon detector and a
five year run time. The DUNE conceptual design report suggests roughly twice
as large a statistical sample will be available per year Hewes  _et al._
(2021). When combined with a ten vs. five year run time this could reduce
statistical errors by a factor of two. Higher statistics, and improved control
over detector systematics may then allow for determination of the neutrino
flux with $\lesssim 1\%$ precision. It is important to have theoretical
control over all possible corrections at a sub-percent level. Estimating the
effect of binding corrections a priori one would expect a correction of size
$\sim Z_{\rm Ar}^{4/3}\alpha^{2}\approx 2.5\times 10^{-3}$.

The total binding energy of argon is given by $\epsilon_{\rm Ar}=14.40~{}{\rm
keV}$ Kramida  _et al._ (2023). Estimating the average inverse radius using
tabulated Hartree-Fock calculations from Ref. Mann (1968), and averaging over
atomic orbitals, we find $\langle V_{1}\rangle_{\rm Ar}=-34.2~{}{\rm keV}$.

We then find for the binding correction to the cross section for scattering on
atomic electrons in argon,

$\begin{split}\differential\sigma_{\cal
B}&=\quantity(-\quantity[\frac{7}{3}\quantity(1.57)-3.71]\times
10^{-3})\differential\sigma_{0}\\\ &\approx(-6\times
10^{-5})\times\differential\sigma_{0}~{}.\end{split}$ (46)

This is much smaller than the previously estimated error on the cross section
($\pm 0.37\%$) Tomalak and Hill (2020). We find due to the delicate
cancellation between binding energy and one-body potential terms, that this
estimate is very sensitive to the precise value of each term. Nevertheless,
the size is always much smaller than $10^{-3}$ and can therefore be safely
neglected for neutrino flux normalization measurements with an argon target.

###  VI.2 Muon electron scattering

MUonE will measure muon electron scattering with tagged incoming muons, and
excellent angular resolution for outgoing electrons and muons. Bound state
effects can be, and indeed are, important for MUonE. The purpose of the
experiment is to measure the running of $\alpha(\mu)$ where $\alpha$ is the
electromagnetic fine structure constant. The strategy may be summarized
schematically as Banerjee  _et al._ (2020)

$\derivative{\sigma}{t}=\derivative{\sigma^{0}}{t}\quantity|\frac{\alpha(t)}{\alpha(0)}|^{2}\times(1+\delta_{R}(t))~{},$
(47)

where $\delta_{R}$ contains radiative corrections. Since hadronic vacuum
polarization is a small effect, the demands for $\delta_{R}(t)$ are
extraordinarily stringent (10 ppm). As we have discussed above in the case of
neutrino scattering, bound-state corrections should also be included i.e., one
should make the replacement
$(1+\delta_{R})\rightarrow(1+\delta_{R}+\delta_{\cal B})$. As we will see, the
bound state corrections to $\mu e\rightarrow\mu e$ have important
phenomenological implications for an extraction of hadronic vacuum
polarization at MUonE.

The angular resolution at MUonE is superb; $0.02~{}{\rm mrad}$ Abbiendi
(2022). It is therefore interesting to ask if the effects from atomic binding
can significantly shift the kinematics of $\mu e\rightarrow\mu e$ scattering
events. We have checked the modified kinematic constraints, and find that for
the $\sim 150~{}{\rm GeV}$ muon beam at MUonE, the shifts in $\theta_{\mu}$
and $\theta_{e}$ are suppressed by $\epsilon_{A}/E_{\mu}\ll 0.02~{}{\rm mrad}$
and are therefore irrelevant for practical considerations.

At tree-level $\mu e\rightarrow\mu e$ is mediated by an off-shell photon,
whose virtuality is given by $q^{2}=(k-k^{\prime})^{2}$ with $k$ the four-
momentum of the incoming muon, and $k^{\prime}$ the momentum of the outgoing
muon. The relevant matrix element in our problem is then
$\matrixelement{e^{-}B^{+}}{\hat{J}_{\mu}}{A}$ with $\hat{J}_{\mu}$ the
(leptonic) electromagnetic current. Following identically to above, for the
reaction

$\mu^{\pm}(\mathbf{k})+A\rightarrow\mu^{\pm}(\mathbf{k}^{\prime})+e^{-}(\mathbf{p}^{\prime})+B(-\mathbf{p})~{},$
(48)

we arrive at

$\begin{split}\sigma=&\int\differential\epsilon\int\frac{\differential^{3}p}{(2\pi)^{3}}S_{A}(\epsilon,\mathbf{p})\\\
&\hskip
108.405pt\times\frac{1}{2E}\frac{1}{2|\mathbf{k}|}\int\frac{\differential
E^{\prime}}{8\pi\omega}~{}\frac{1}{4}\sum_{\rm spins}|{\sf
M}|^{2}~{},\end{split}$ (49)

where the factor of $1/4$ comes from averaging over both electron and muons
spins in the initial state. Neglecting $m_{e}/\omega$ and
$m_{e}^{2}/m_{\mu}^{2}$, the spin-summed matrix element is given by

$\sum_{\rm spins}|{\sf M}|^{2}\simeq 32e^{4}\frac{(k\cdot
p^{\prime})(k^{\prime}\cdot p)+(k\cdot p)(k^{\prime}\cdot
p^{\prime})-m_{\mu}^{2}(p\cdot
p^{\prime})}{\quantity[(k-k^{\prime})^{2}]^{2}}~{}.$

When expanding Lorentz invariants in terms of $\mathbf{p}$ we must be more
careful than in the case of neutrino scattering. Terms proportional to
$\cos\theta$ can appear from both the numerator and denominator. One should
therefore keep terms linear in $\mathbf{p}$, and the appropriate identities
are

$\displaystyle(p\cdot k)$ $\displaystyle\simeq(p\cdot
k\phantom{{}^{\prime}})_{0}\quantity[1-\frac{|\mathbf{p}|\cos\theta}{m_{e}}+\frac{\mathbf{p}^{2}}{2m_{e}^{2}}]~{},$
(50) $\displaystyle(p\cdot k^{\prime})$ $\displaystyle\simeq(p\cdot
k^{\prime})_{0}\quantity[1-\frac{|\mathbf{p}|\cos\theta}{m_{e}}+\frac{\mathbf{p}^{2}}{2m_{e}^{2}}]~{},$
(51) $\displaystyle(p\cdot p^{\prime})$ $\displaystyle\simeq(p\cdot
p^{\prime})_{0}\quantity[1-\frac{|\mathbf{p}|\cos\theta}{m_{e}}+\frac{\mathbf{p}^{2}}{2m_{e}^{2}}]~{},$
(52)

where we (again) use
$\cos\theta_{ee^{\prime}}\eqsim\cos\theta_{e\nu^{\prime}}\eqsim\cos\theta_{e\nu}\equiv\cos\theta$.
Next, we expand in $\mathbf{p}^{2}$ and $\epsilon$. After averaging over
angles in the bound state, we may make the replacement,

Figure 1: The function $f(t)$, as defined in Eq. 55, for $\sqrt{s}=405~{}{\rm
MeV}$, plotted over the interval $t_{\rm min}\leq t\leq 0$ where the largest
momentum transfer is given by $t_{\rm min}=-(s-m_{\mu}^{2})^{2}/s$.

$\begin{split}\frac{1}{2E}\sum_{\rm spins}|{\sf
M}|^{2}\rightarrow&\frac{1}{2m}\bigg{\\{}~{}\sum_{\rm spins}|{\sf
M}|^{2}_{0}\quantity[1+\frac{\epsilon}{m_{e}}]\\\ &\hskip
65.04034pt-\frac{8e^{4}m_{\mu}^{2}}{(p\cdot
p^{\prime})_{0}}\quantity[\frac{\epsilon}{m_{e}}+\frac{1}{3}\frac{\mathbf{p}^{2}}{m_{e}^{2}}]~{}\bigg{\\}}.\end{split}$
(53)

Upon integration against the spectral function and using the relevant sum
rules, we find,

$\begin{split}\derivative{\sigma}{t}&\simeq\derivative{\sigma^{(0)}}{t}\bigg{[}1+\frac{1}{Z_{A}m_{e}}\quantity[3\epsilon_{A}+\langle\hat{V}_{1}\rangle_{A}]\\\
&\hskip
130.08731pt-\frac{f(t)}{Z_{A}m_{e}}\quantity(\frac{11}{3}\epsilon_{A}+\langle\hat{V}_{1}\rangle_{A})\bigg{]}~{},\end{split}$
(54)

where the Mandelstam variable $t$ for an electron at rest is defined as
$t\simeq-2(p\cdot p^{\prime})_{0}=-2m_{e}E^{\prime}$ in the limit of small
electron mass. The function $f(t)$ (plotted in Fig. 1 for $\sqrt{s}=405~{}{\rm
MeV}$) is obtained from Eq. 53 after dividing through by the leading order
answer,

$\displaystyle\begin{split}f(t)&=\frac{m_{\mu}^{2}(p\cdot
p^{\prime})_{0}}{(p\cdot k)_{0}^{2}+(p\cdot
k^{\prime})_{0}^{2}-m_{\mu}^{2}(p\cdot p^{\prime})_{0}}\\\
&=\frac{-2m_{\mu}^{2}t}{(s+t-m_{\mu}^{2})^{2}+(s-m_{\mu})^{2}+2m_{\mu}^{2}t}~{}.\end{split}$
(55)

For a shape-only measurement, the relevant quantity is

$\begin{split}\frac{1}{\sigma}\derivative{\sigma}{t}&\simeq\frac{1}{\sigma^{(0)}}\derivative{\sigma^{(0)}}{t}\bigg{(}1-\frac{f(t)}{Z_{A}m_{e}}\quantity[\frac{11}{3}\epsilon_{A}+\langle\hat{V}_{1}\rangle_{A}]\bigg{)}~{},\end{split}$
(56)

To get a sense of the size of this correction we may assume a carbon target,
taking $\epsilon_{\rm C}=1.03~{}{\rm keV}$ Kramida  _et al._ (2023) and
$\langle V_{1}\rangle_{\rm C}=-2.40~{}{\rm keV}$ Mann (1968). We find,

$\frac{1}{Z_{A}m_{e}}\quantity[\frac{11}{3}\epsilon_{A}+\langle\hat{V}_{1}\rangle_{A}]\approx
45\times 10^{-5}~{}.$ (57)

This is almost two orders of magnitude larger than the theory-error target of
MUonE. Using $\sqrt{s}=405~{}{\rm MeV}$ and $\sqrt{-t}=\sqrt{-t_{\rm
peak}}=330~{}{\rm MeV}$ Abbiendi (2022) we find that this correction amounts
to a $5\times 10^{-5}$ shift which is a sizeable effect, when compared to the
$10^{-5}$ error budget demanded by MUonE Banerjee  _et al._ (2020). Notice
that there is no delicate cancellation in Eq. 57 in contrast to the accidental
cancellation in Eq. 46, and the binding correction is of natural size.

##  VII Discussion and conclusions

Using a simple perturbative analysis, the virial theorem, and the Koltun sum
rule Koltun (1972), we have identified a new model independent relationship
between the expectation value of the single body potential operator, and
binding corrections to high energy scattering. We have studied these
corrections in detail for $\nu e\rightarrow\nu e$ and $\mu e\rightarrow\mu e$
scattering. We find that due to rotational symmetry, corrections begin at
$O(\epsilon_{A}/m_{e})$ where $\epsilon_{A}$ is the binding energy of the
target atom. Our main results are Eq. 43 for $\nu e\rightarrow\nu e$
scattering, and Eqs. 54 and 56 for $\mu e\rightarrow\mu e$ scattering.

In our analysis we have treated the final state electron as a free-particle
solution (which amounts to using the impulse approximation). In general there
will be perturbative corrections from Coulomb exchange with the final state
system $\ket{B}$ that are not included in the radiative corrections to $\nu
e\rightarrow\nu e$ and $\mu e\rightarrow\mu e$ scattering involving free
electrons. The effect of final state interactions can be estimated in the case
of hydrogen, but it is presently unclear if any universal form exists for
many-body atoms. These corrections should be computed in the future and added
to the binding corrections discussed in this paper.

The results presented here are important for ultra-precise measurements of
muon electron scattering, as is planned for the MUonE experiment. Unlike in
neutrino flux measurements, where the total cross section is most important,
for MUonE the most relevant observable is the differential distribution with
respect to the angles of outgoing muon and electron. The binding corrections
discussed here will impact extractions of hadronic vacuum polarization from
MUonE. Importantly, the error budget there is $\sim 10^{-5}$ and even for
light nuclei (e.g., carbon with $Z^{4/3}\approx 11$) atomic binding
corrections must be incorporated.

Finally, we note that while we have studied atoms, the lessons learned here
may also be useful for nuclei. For instance, binding corrections must be
applied to neutrino deuterium scattering data to extract the axial form factor
Meyer  _et al._ (2016). Deuterium is a weakly bound system whose dynamics are
reasonably well understood, and an expansion in $\epsilon_{D}/m_{N}$, with
$\epsilon_{D}$ the deuterium binding energy and $m_{N}$ the nucleon mass,
should be achievable; the inclusion of final state interactions may be more
challenging. It would be interesting to understand if similar energy weighted
sum rules can help supply a model independent set of deuterium binding
corrections.

## Acknowledgements

RP thanks Alexis Nikolakopoulos for helpful discussions. RP is supported by
the Neutrino Theory Network under Award Number DEAC02-07CH11359. RP and MW are
supported by the U.S. Department of Energy, Office of Science, Office of High
Energy Physics under Award Number DE-SC0011632, and by the Walter Burke
Institute for Theoretical Physics.

##  Appendix A Sum rules

Since the Koltun sum rule Koltun (1972) is very simple, but may be unfamiliar
to some readers, we provide a self contained discussion (see also Ref. Polls
_et al._ (1994)). Let us consider the hole spectral function for a state $A$

$\displaystyle S_{A}(\epsilon,\mathbf{p})=$
$\displaystyle\int\frac{\differential^{3}q}{(2\pi)^{3}}\sumint_{B}\matrixelement{A}{\hat{a}^{\dagger}_{\mathbf{p}}}{B}\delta(\epsilon_{A}-\epsilon_{B}-\epsilon)\matrixelement{B}{\hat{a}_{\mathbf{q}}}{A}~{}.$
(58)

The first (trivial) sum rule stems from the identity
$\int\differential\epsilon\delta(x-\epsilon)f(x)=f(x)$ for any test function.
Therefore,

$\begin{split}\int\differential\epsilon
S_{A}(\epsilon,\mathbf{p})&=\int\frac{\differential^{3}q}{(2\pi)^{3}}\bra{A}{\hat{a}^{\dagger}_{\mathbf{p}}}\sumint_{B}\outerproduct{B}{B}{\hat{a}_{\mathbf{q}}}\ket{A}~{}\\\
&=\int\frac{\differential^{3}q}{(2\pi)^{3}}\matrixelement{A}{\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{q}}}{A}~{}.\end{split}$
(59)

When weighted against $\mathbf{p}^{2}/2m_{e}$ and integrated over
$\differential^{3}p$ this gives $\langle\hat{T}\rangle$.

Next, for the Koltun sum rule, we re-write $\epsilon$ in terms of the
Hamiltonian via (suppressing $\int\differential^{3}q/(2\pi)^{3}$ for brevity’s
sake),

$\displaystyle~{}~{}~{}~{}\epsilon\sumint_{B}\matrixelement{A}{\hat{a}^{\dagger}_{\mathbf{q}}}{B}\delta(\epsilon_{A}-\epsilon_{B}-\epsilon)\matrixelement{B}{\hat{a}_{\mathbf{p}}}{A}~{}$
(60)
$\displaystyle=(\epsilon_{A}-\epsilon_{B})\sumint_{B}\matrixelement{A}{\hat{a}^{\dagger}_{\mathbf{q}}}{B}\delta(\epsilon_{A}-\epsilon_{B}-\epsilon)\matrixelement{B}{\hat{a}_{\mathbf{q}}}{A}~{},$
$\displaystyle=\sumint_{B}\matrixelement{A}{\hat{a}^{\dagger}_{\mathbf{p}}}{B}\delta(\epsilon_{A}-\epsilon_{B}-\epsilon)\matrixelement{B}{\hat{a}_{\mathbf{q}}\epsilon_{A}-\epsilon_{B}\hat{a}_{\mathbf{q}}}{A}~{},$
$\displaystyle=\sumint_{B}\matrixelement{A}{\hat{a}^{\dagger}_{\mathbf{p}}}{B}\delta(\epsilon_{A}-\epsilon_{B}-\epsilon)\matrixelement{B}{\hat{a}^{\dagger}_{\mathbf{p}}[\hat{H},\hat{a}_{\mathbf{q}}]}{A}~{}.$

where we have used $\hat{H}\ket{A,B}=-\epsilon_{A,B}\ket{A,B}$. Then, using
$(-\epsilon)$ instead of $\epsilon$ and integrating we obtain

$\int\differential\epsilon~{}(-\epsilon)S_{A}(\epsilon,\mathbf{p})=\int\frac{\differential^{3}q}{(2\pi)^{3}}\matrixelement{A}{\hat{a}^{\dagger}_{\mathbf{p}}[\hat{a}_{\mathbf{q}},\hat{H}]}{A}~{}.$
(61)

Taking a Hamiltonian with one- and two-body operators,

$\begin{split}\hat{H}=&\sum_{\mathbf{p}}\frac{\mathbf{p}^{2}}{2m_{e}}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}+\sum_{\mathbf{p},\mathbf{p}^{\prime}}V_{1}(\mathbf{p},\mathbf{p}^{\prime})\hat{a}_{\mathbf{p}^{\prime}}^{\dagger}\hat{a}_{\mathbf{p}}\\\
&~{}~{}~{}~{}+\sum_{\mathbf{p},\mathbf{p}^{\prime},\mathbf{q},\mathbf{q}^{\prime}}V_{2}(\mathbf{p},\mathbf{p}^{\prime},\mathbf{q},\mathbf{q}^{\prime})\hat{a}_{\mathbf{p}^{\prime}}^{\dagger}\hat{a}_{\mathbf{q}^{\prime}}^{\dagger}\hat{a}_{\mathbf{p}}\hat{a}_{\mathbf{q}}~{},\end{split}$
(62)

evaluating the commutator, and integrating over $\differential\epsilon$ to
remove the delta function, we arrive at

$\int\frac{\differential^{3}p}{(2\pi)^{3}}\int\\!\differential\epsilon(-\epsilon)S_{A}(\epsilon,\mathbf{p})=\frac{\matrixelement{A}{\hat{T}+\hat{V}_{1}+2\hat{V}_{2}}{A}}{\innerproduct{A}{A}}.$
(63)

For a more general Hamiltonian written as an expansion $n$-body operators as
$\hat{H}=\sum_{n}\hat{\mathcal{O}}_{n}$ one finds Koltun (1972)

$\int\frac{\differential^{3}p}{(2\pi)^{3}}\int\differential\epsilon(-\epsilon)S_{A}(\epsilon,\mathbf{p})=\sum_{n}n\frac{\matrixelement{A}{\hat{\cal
O}_{n}}{A}}{\innerproduct{A}{A}}~{}.$ (64)

##  Appendix B Thomas-Fermi estimate

The Thomas-Fermi model offers a crude description of bulk atomic properties
March and Haar (2016). Since we have reduced the problem of binding
corrections to evaluating $\langle\hat{V}_{1}\rangle_{A}$ the quality of the
estimate depends only on the accuracy with which we can model
$n_{A}(\mathbf{r})$.

The Thomas-Fermi model can be summarized as $r^{2}n(\mathbf{r})\propto
x^{1/2}[f(x)]^{3/2}$ where $x=r/b_{0}$ with $b_{0}=0.88Z^{-1/3}a_{0}$. The
function $f(x)$ is determined by solving the differential equation,

$xf^{\prime\prime}(x)=x^{1/2}[f(x)]^{3/2}~{},$ (65)

with $f(0)=1$ and $f^{\prime}(0)\simeq 1.588$. Using this solution, and
integrating numerically, one finds

$\frac{\int\differential
x~{}\quantity(\frac{1}{x})x^{1/2}f(x)}{\int\differential
x~{}x^{1/2}f(x)}\approx 1.79~{}.$ (66)

Converting to physical units gives Eq. 45.

##  Appendix C Relativistic corrections

In Eq. 13 we have written a hydrogen (or hydrogen-like) atom in terms of an
entangled proton (nucleus) and a free-electron state. In the limit of
$m_{p}/m_{e}\rightarrow\infty$ the proton (nucleus) acts as a static Coulomb
field and the exact solution (to all orders in $Z\alpha$) is given by the
solution of the Dirac equation in a central potential. These wavefunctions
generically cannot be written in the form of Eq. 13. Nevertheless, as we
dicuss below, these effects begin first at $O(p^{3}/m_{e}^{3})$ and are
therefore sub-dominant in our expansion.

Consider the Dirac equation in a central potential. Using the “Bjorken and
Drell” basis this can be written as

$\begin{pmatrix}\phantom{2m}-\epsilon_{H}+V(\mathbf{x})&-\mathrm{i}{\bf\it\sigma}\cdot{\bf\it\nabla}\\\
\mathrm{i}{\bf\it\sigma}\cdot{\bf\it\nabla}&2m-\epsilon_{H}+V(\mathbf{x})\end{pmatrix}\begin{pmatrix}U(\mathbf{x})\\\
L(\mathbf{x})\end{pmatrix}=0~{}.$ (67)

Non-relativistic bound states satisfy the counting ${\bf\it\nabla}\sim
O(\lambda)$, $V(\mathbf{x})\sim O(\lambda^{2})$, and $\epsilon_{H}\sim
O(\lambda^{2})$, as required by the virial theorem. We may then expand order-
by-order in $\lambda$,

$\displaystyle U(\mathbf{x})$
$\displaystyle=U^{(0)}(\mathbf{x})+U^{(2)}(\mathbf{x})+U^{(4)}(\mathbf{x})+\ldots~{},$
(68) $\displaystyle L(\mathbf{x})$
$\displaystyle=L^{(1)}(\mathbf{x})+L^{(3)}(\mathbf{x})+L^{(5)}(\mathbf{x})+\ldots~{},$
(69) $\displaystyle\epsilon_{H}$
$\displaystyle=\epsilon_{H}^{(2)}+\epsilon_{H}^{(4)}+\epsilon_{H}^{(6)}+\ldots~{},$
(70)

where we have anticipated the vanishing of alternating orders. The fact that
${\bf\it\nabla}$ is the dominant term in the expansion implies the constraint
at leading order,

$L^{(1)}(\mathbf{x})=\frac{-\mathrm{i}{\bf\it\sigma}\cdot{\bf\it\nabla}}{2m}U^{(0)}(\mathbf{x})~{}.$
(71)

The zeroth order upper component satisfies the Schrödinger equation,

$\quantity[-\frac{1}{2m}\nabla^{2}+V(\mathbf{x})]U^{(0)}(\mathbf{x})=-\epsilon_{H}^{(0)}U^{(0)}(\mathbf{x})~{}.$
(72)

The second order correction $U^{(2)}(\mathbf{x})$ incorporates fine-structure
effects including, for example, spin-orbit coupling. We may therefore, through
$O(\lambda^{2})=O(\epsilon_{H}/m_{e})$, write the hydrogen-like Dirac
wavefunction as,

$\Psi(\mathbf{x})=\begin{pmatrix}\psi_{H}(\mathbf{x})\\\
\frac{-\mathrm{i}{\bf\it\sigma}\cdot{\bf\it\nabla}}{2m_{e}}\psi_{H}(\mathbf{x})\end{pmatrix}+O(\lambda^{3})~{},$
(73)

where $\psi_{H}(\mathbf{x)}\equiv U^{(0)}(\mathbf{x})+U^{(2)}(\mathbf{x})$ is
a two-component spinor. Equivalently in momentum space we find

$\Psi(\mathbf{p})=\begin{pmatrix}1\\\
\frac{{\bf\it\sigma}\cdot\mathbf{p}}{2m_{e}}\end{pmatrix}\tilde{\psi}_{H}(\mathbf{p})+O(\lambda^{3})~{}.$
(74)

We see that up to a normalization of states (i.e., relativistic vs. non-
relativistic conventions) this is equivalent to using Eq. 13 up to
$O(\lambda^{3})$ corrections.

When computing corrections to the matrix element, the contributions from
$U^{(2)}$ enter only via
$\int\differential^{3}p/(2\pi)^{3}|\tilde{\psi}_{H}(\mathbf{p})|^{2}|{\sf
M}_{0}|^{2}$. Since ${\sf M}_{0}$ is independent of $\mathbf{p}$ (by
definition) there is no correction, because the norm of the wavefunction is
fixed to all orders
$\int\differential^{3}p/(2\pi)^{3}|\tilde{\psi}(\mathbf{p})|^{2}=1$. When
computing the terms proportional to $\mathbf{p}^{2}/m_{e}^{2}$ and
$\epsilon/m_{e}$ in e.g., Eq. 26 it is sufficient to use the zeroth-order
approximation to the wavefunction.

The above analysis makes it clear that this approximations (likely) holds to
the same accuracy in many-body atoms. One can write the full QED Hamiltonian,
and similarly expand the result order-by-order in $\lambda$. The constraint
Eq. 71 will appear independent of the details of the interactions between
electrons, all of which are counted as $O(\lambda^{2})$.

The fully relativistic spectral function satisfies the exact sum rule
$\int\differential^{3}p/(2\pi)^{3}\int\differential\epsilon~{}S(\epsilon,\mathbf{p})=Z$.
When computing corrections due to bound states, i.e., for terms proportional
to $\mathbf{p}^{2}/m_{e}^{2}$ or $\epsilon/m_{e}$, it suffices to use the non-
relativistic approximation of the spectral function.

## References

  * Hewes  _et al._ (2021) V. Hewes _et al._ (DUNE), Instruments 5, 31 (2021), arXiv:2103.13910 [physics.ins-det] .
  * Abbiendi (2022) G. Abbiendi, Phys. Scripta 97, 054007 (2022), arXiv:2201.13177 [physics.ins-det] .
  * Banerjee  _et al._ (2020) P. Banerjee _et al._ , Eur. Phys. J. C 80, 591 (2020), arXiv:2004.13663 [hep-ph] .
  * Bechler (1993) A. Bechler, Journal of Physics A: Mathematical and General 26, 6039 (1993).
  * Polls  _et al._ (1994) A. Polls, A. Ramos, J. Ventura, S. Amari, and W. H. Dickhoff, Phys. Rev. C 49, 3050 (1994).
  * Koltun (1972) D. S. Koltun, Phys. Rev. Lett. 28, 182 (1972).
  * Marshall  _et al._ (2020) C. M. Marshall, K. S. McFarland, and C. Wilkinson, Phys. Rev. D 101, 032002 (2020), arXiv:1910.10996 [hep-ex] .
  * Kramida  _et al._ (2023) A. Kramida, Yu. Ralchenko, J. Reader, and and NIST ASD Team, NIST Atomic Spectra Database (ver. 5.11), [Online]. https://physics.nist.gov/asd [2024, March 12]. National Institute of Standards and Technology, Gaithersburg, MD. (2023).
  * Mann (1968) J. B. Mann,  _Atomic Structure Calculations II. Hartree-Fock Wavefunction and Radial Expectation Values: Hydrogen to Lawrencium_, Tech. Rep. LA-3691 (Los Alamos National Laboratory, 1968).
  * Tomalak and Hill (2020) O. Tomalak and R. J. Hill, Phys. Rev. D 101, 033006 (2020), arXiv:1907.03379 [hep-ph] .
  * Meyer  _et al._ (2016) A. S. Meyer, M. Betancourt, R. Gran, and R. J. Hill, Phys. Rev. D 93, 113015 (2016), arXiv:1603.03048 [hep-ph] .
  * March and Haar (2016) N. March and D. Haar,  _Self-Consistent Fields in Atoms: Hartree and Thomas–Fermi Atoms_ (Elsevier Science, 2016).