[a]A. Hannaford-Gunn, [a] E. Sankey,

# Investigating the Compton amplitude subtraction function in lattice QCD

K. U. Can  R. Horsley  H. Perlt  P. E. L. Rakow  G. Schierholz  K. Somfleth
H. Stüben  R. D. Young  and J. M. Zanotti  for the CSSM-QCDSF-UKQCD
Collaboration

###### Abstract

Theoretical predictions of the proton–neutron mass difference and measurements
of the proton’s charge radius require inputs from the Compton amplitude
subtraction function. Model-dependent and non-relativistic calculations of
this subtraction function vary significantly, and hence it contributes
sizeable uncertainties to the aforementioned physical quantities. We report on
the use of Feynman-Hellmann methods in lattice QCD to calculate the
subtraction function from first principles. In particular, our initial results
show anomalous high-energy behaviour that is at odds with the prediction from
the operator product expansion (OPE). Therefore, we investigate the
possibility that this unexpected behaviour is due to lattice artifacts, by
varying the lattice spacing and volume, and comparing different
discretisations of the vector current. Finally, we explore a Feynman-Hellmann
implementation that is less sensitive to short-distance contributions and show
that the subtraction function’s anomalous behaviour can be attributed to these
short-distance contributions. As such, this work represents the first steps in
achieving a complete understanding of the Compton amplitude subtraction
function.

##  1 Introduction

$p$$q$$q$$p$

Figure 1: Forward Compton scattering for a proton with momentum $p$ and
virtual photon with momentum $q$.

The forward Compton amplitude is a necessary input into two important physical
quantities: predictions for the proton–neutron mass difference, and the
hadronic background for measurements of the proton charge radius. All the
components of this Compton amplitude can be determined experimentally, except
for its subtraction function, $S_{1}(Q^{2})$. Therefore, a first-principles
lattice QCD calculation of this subtraction function is of great interest.

The mass difference of the proton and neutron has two sources: the different
masses of the up and down quarks, and the different charges of these
quarks—see Refs. [1, 2] for reviews. The leading electromagnetic contribution
to the mass difference (Fig. 2(a)) can be evaluated from the Cottingham sum
rule [3],

$\delta m^{\text{EM}}=-\frac{i}{2m_{p}}\frac{\alpha}{(2\pi)^{3}}\int
d^{4}q\frac{T^{\mu}_{\mu}(p,q)}{Q^{2}-i\epsilon},$ (1.1)

where $T_{\mu\nu}$ is the forward, spin-averaged Compton amplitude for a
proton:

$T_{\mu\nu}(p,q)\equiv\frac{i}{2}\sum_{s}\int d^{4}z\ e^{iq\cdot z}\langle
p,s|{T}\\{j_{\mu}(z)j_{\nu}(0)\\}|p,s\rangle.$ (1.2)

This amplitude describes the process of photon-proton scattering,
$\gamma^{(*)}(q)P(p)\to\gamma^{(*)}(q)P(p)$, with no momentum transfer between
initial and final states (Fig. 1).

Similarly, the Compton amplitude is required to constrain measurements of the
proton charge radius from the muonic-hydrogen Lamb shift. Crucially, recent
determinations of the charge radius from this Lamb shift conflict with
previous results obtained via electron–proton scattering [4] by seven standard
deviations [5]—the so-called ‘proton radius puzzle’ [6, 7]. The hadronic
corrections to the Lamb shift are dependent on the two-photon-exchange diagram
(Fig. 2(b)) [8, 9]:

$\mathcal{M}_{\text{TPE}}=-ie^{4}\int\frac{d^{4}q}{(2\pi)^{4}}\frac{T_{\mu\nu}L^{\mu\nu}}{Q^{4}-i\epsilon},$
(1.3)

where $L_{\mu\nu}$ is the leptonic contribution, which can be calculated from
QED, and $T_{\mu\nu}$ is the proton Compton amplitude. Since the Compton
subtraction function is poorly constrained, it contributes the dominant
uncertainty to the hadronic background [10, 11, 12]. Hence more precise
determinations of the subtraction function could help clarify the proton
radius puzzle.

It has been conjectured that the subtraction function receives contributions
from a $J=0$ fixed pole, which is independent of $Q^{2}$ [13]. Such a
singularity may arise from the exchange of a particle of spin-zero or from a
contact interaction term. If correct, this would have far-reaching
phenomenological consequences [14], and potentially thwart the Cottingham sum
rule. However, this conjecture has remained rather controversial [15].

The starting point for our study is the spin-averaged Compton amplitude, Eq.
(1.2), which can be decomposed into two structures:

$T_{\mu\nu}(p,q)=\left(-g_{\mu\nu}+\frac{q_{\mu}q_{\nu}}{q^{2}}\right)\mathcal{F}_{1}(\omega,Q^{2})+\left(p_{\mu}-\frac{p\cdot
q}{q^{2}}q_{\mu}\right)\left(p_{\nu}-\frac{p\cdot
q}{q^{2}}q_{\nu}\right)\frac{\mathcal{F}_{2}(\omega,Q^{2})}{p\cdot q},$ (1.4)

where $\mathcal{F}_{1,2}$ are the _Compton structure functions_ , defined in
terms of the photon virtuality, $Q^{2}={-q^{2}}$, and the inverse Bjorken
scaling variable, $\omega=2(p\cdot q)/Q^{2}$. These structure functions
satisfy the following fixed-$Q^{2}$ dispersion integrals [16]:

$\displaystyle\overline{\mathcal{F}}_{1}(\omega,Q^{2})=\mathcal{F}_{1}(\omega,Q^{2})-\mathcal{F}_{1}(0,Q^{2})$
$\displaystyle=2\omega^{2}\int_{0}^{1}dx\frac{2xF_{1}(x,Q^{2})}{1-x^{2}\omega^{2}-i\epsilon},$
(1.5) $\displaystyle\mathcal{F}_{2}(\omega,Q^{2})$
$\displaystyle=4\omega\int_{0}^{1}dx\frac{F_{2}(x,Q^{2})}{1-x^{2}\omega^{2}-i\epsilon},$
(1.6)

where $F_{1,2}(x=\omega^{-1},Q^{2})$ are measurable from deep inelastic
scattering (DIS) cross sections.

However, equation (1.5) features a contribution from the _Compton amplitude
subtraction function_ ,

$S_{1}(Q^{2})\equiv\mathcal{F}_{1}(\omega=0,Q^{2}),$ (1.7)

which is not experimentally accessible. Instead, this subtraction function has
been determined from model-dependent, dispersive, non-relativistic, and
effective theory calculations [17, 18, 19, 10, 20, 21, 12, 11, 22, 23, 24, 25,
26]. These calculations have sizeable errors and are not always consistent
with one another [27]. At large $Q^{2}$, the subtraction function can be
evaluated model-independently using the operator product expansion (OPE) [28,
29], and the following asymptotic behaviour is predicted:

$S_{1}(Q^{2})\sim\frac{m_{N}^{2}}{Q^{2}},\quad\text{for}\quad Q^{2}\gg
m_{N}^{2}.$ (1.8)

However, as can be seen in Eqs. (1.1) and (1.3), a determination for the whole
domain of $S_{1}(Q^{2})$ is necessary for inputs into the aforementioned
physical quantities. Therefore, a determination of this subtraction function,
particularly for low and intermediate $Q^{2}$, is of great interest.

$p$$p$

(a) Electromagnetic self-energy

$k$$k$$p$$p$$q$$-q$

(b) Two-photon exchange

Figure 2: Two processes which have dependence on hadronic structure given by
the Compton amplitude.

Recently, the Feynman-Hellmann method has been used to calculate the forward
Compton amplitude in lattice QCD [30, 31], and has been extended to off-
forward kinematics [32]. In this report, we apply this technique to calculate
the Compton subtraction function. In particular, our initial results display
large deviations from the asymptotic behaviour predicted by the OPE, Eq.
(1.8). Therefore, we conduct an investigation of lattice artifacts of the
subtraction function. We vary the lattice volume and spacing, and compare
different discretisations of the vector current. Finally, we present a novel
method to explore the effects of current-current contact terms in Feynman-
Hellmann methods. Hence this work lays the foundation for future lattice QCD
calculations of the Compton subtraction function, which would allow for better
theoretical predictions of the proton–neutron mass difference, and more
precise experimental determinations of the proton charge radius.

##  2 Feynman-Hellmann: Local vector current implementation

Feynman-Hellmann methods provide a feasible alternative to the direct
calculation of four-point functions. A background field is introduced to the
quark propagator, with a small unphysical coupling, $\lambda$. At order
$\lambda^{2}$, the perturbation to the propagator is a four-point function
with a sum over time-slices on which the background field is inserted. As
such, one set of Feynman-Hellmann inversions yields the sum over time-slices
necessary to calculate a discretisation of the Compton amplitude, Eq. (1.2).
By contrast, a direct four-point function evaluation of this amplitude would
require $\mathcal{O}(T^{2})$ inversions to get all possible insertion times
for a lattice of temporal extent $T$.

It was shown in Ref. [31] that for a perturbed nucleon propagator,
$\mathcal{G}_{\lambda}(t)\simeq A_{\lambda}e^{-E_{\lambda}t}$, the perturbed
energy, $E_{\lambda}$, can be related to the Compton amplitude subtraction
function:

$-m_{N}\frac{\partial^{2}E_{\lambda}}{\partial\lambda^{2}}\bigg{|}_{\lambda=0}=S_{1}(Q^{2}),$
(2.1)

in the nucleon’s rest frame: $\vec{p}=0$. The derivation of this Feynman-
Hellmann relation, Eq. (2.1), was carried out in terms of hadronic states with
continuous spacetime coordinates [31]. Since our aim here will be on
understanding and controlling several lattice artifacts, we instead focus on
the Feynman-Hellmann implementation at the level of lattice quark propagators.

A common discretisation of the electromagnetic current is the _local current_
:

$j_{\mu}^{\text{loc}}(n)=\bar{\psi}(n)\gamma_{\mu}\psi(n).$ (2.2)

To calculate the Compton amplitude, Eq. (1.2), with this discretisation, we
compute the following perturbed propagator:

$\begin{split}S_{{\lambda}}=\big{[}M-\lambda\mathcal{O}_{\vec{q}}\big{]}^{-1},\end{split}$
(2.3)

where the perturbing matrix is

$\big{[}\mathcal{O}_{\vec{q}}\big{]}_{n,m}=\delta_{n,m}\phi_{\vec{q}}(\vec{n})i\gamma_{3},\quad\text{with}\quad\phi_{\vec{q}}(\vec{n})=e^{i\vec{q}\cdot\vec{n}}+e^{-i\vec{q}\cdot\vec{n}},$
(2.4)

for $\lambda$ a small coupling. $M$ is the fermion matrix, which can be
written, up to the clover term [33], as

$M(n,m)\equiv\left(m_{0}+4\right)\delta_{n,m}-\frac{1}{2a}\sum_{\mu}\bigg{[}\vec{\mathcal{D}}_{\mu}(n,m)+\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{\mu}(n,m)\bigg{]},$
(2.5)

with

$\vec{\mathcal{D}}_{\mu}(n,m)\equiv(1-\gamma_{\mu})U_{\mu}(n)\delta_{n+\hat{\mu},m},\quad\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{\mu}(n,m)\equiv(1+\gamma_{\mu})U_{\mu}(n-\hat{\mu})^{\dagger}\delta_{n-\hat{\mu},m}.$

Therefore, the second derivative of the perturbed propagator is

$\frac{\partial^{2}}{\partial\lambda^{2}}{S_{{\lambda}}}\bigg{|}_{{\lambda}=0}=2M^{-1}\mathcal{O}_{{\vec{q}}}M^{-1}\mathcal{O}_{{\vec{q}}}M^{-1},$
(2.6)

which is a four-point function with the local current, Eq. (2.2).

In terms of the quark propagators, the nucleon propagators with the
perturbation to the $u$ or $d$ quark111Since the perturbation is only applied
to the fermion propagators and not to the sea quarks, these and all other
results in this report are connected only. are written, up to spin/colour
structure, as

$\mathcal{G}^{uu}_{\lambda}\sim\big{\langle}S^{u}_{\lambda}S^{u}_{\lambda}S^{d}\big{\rangle},\quad\mathcal{G}^{dd}_{\lambda}\sim\big{\langle}S^{u}S^{u}S^{d}_{\lambda}\big{\rangle},\quad\mathcal{G}^{ud}_{\lambda_{1},\lambda_{2}}\sim\big{\langle}S^{u}_{\lambda_{1}}S^{u}_{\lambda_{1}}S^{d}_{\lambda_{2}}\big{\rangle}.$
(2.7)

Recalling the spectral representation of the nucleon perturbed propagator,
$\mathcal{G}^{qq^{\prime}}_{\lambda}(t)\simeq
A^{qq^{\prime}}_{\lambda}e^{-E^{qq^{\prime}}_{\lambda}t}$, we have the
following flavour-dependent Feynman-Hellmann relations:

$-m_{N}\frac{\partial^{2}E^{qq}_{\lambda}}{\partial\lambda^{2}}\bigg{|}_{\lambda=0}=S^{qq}_{1}(Q^{2}),\quad-2m_{N}\frac{\partial^{2}E^{ud}_{\lambda}}{\partial\lambda_{1}\partial\lambda_{2}}\bigg{|}_{\lambda_{1}=\lambda_{2}=0}=S^{ud}_{1}(Q^{2})+S_{1}^{du}(Q^{2}).$
(2.8)

Therefore, to calculate the proton contribution to the subtraction function,
we take

$S_{1}^{p}(Q^{2})=\sum_{q,q^{\prime}=u,d}e_{q}e_{q^{\prime}}S^{qq^{\prime}}_{1}(Q^{2})=\frac{4}{9}S_{1}^{uu}(Q^{2})+\frac{1}{9}S_{1}^{dd}(Q^{2})-\frac{2}{9}\big{(}S_{1}^{ud}(Q^{2})+S_{1}^{du}(Q^{2})\big{)}.$
(2.9)

Table 1: Gauge ensemble details $\begin{array}[]{c c c c c c c c
c}N_{f}&L^{3}\times
T&L[\text{fm}]&a[\text{fm}]&\beta&\kappa&m_{\pi}[\text{GeV}]&Z_{V}\\\
\hline\cr\hline\cr 2+1&32^{3}\times 64&2.4&0.074&5.50&0.120900&0.47&0.86\\\
\hline\cr 2+1&48^{3}\times 96&3.3&0.068&5.65&0.122005&0.41&0.87\\\
&48^{3}\times 96&2.8&0.058&5.80&0.122810&0.43&0.88\\\ \end{array}$

### Results: Local Current Implementation

The simulations with the local current implementation were carried out on
three different gauge ensembles generated by the QCDSF/UKQCD Collaborations
[34] (Tab. 1), with varying volume, lattice spacing and quark masses. All
three ensembles are at the SU(3) flavour symmetric point:
$\kappa_{l}=\kappa_{s}=\kappa$. The inserted momentum is always chosen to be
of the form $q_{\mu}=(0,\vec{q})$. Therefore, to calculate the subtraction
function, for which $p\cdot q=0$, we simply choose our nucleon sink momentum
$\vec{p}=0$. For the $32^{3}\times 64$ lattice, we calculate a large range of
$Q^{2}$ values, while for the two larger volumes, there are fewer. For the
majority of points, the statistics are
$N_{\text{meas}}\sim\mathcal{O}(100)-\mathcal{O}(1000)$, with some higher
statistics points, $N_{\text{meas}}\sim\mathcal{O}(10000)$, for the smallest
volume.

The results are presented in Fig. 3, and show an asymptotic behaviour that
clearly deviates from the OPE prediction of $S_{1}(Q^{2})\sim Q^{-2}$, given
in Eq. (1.8). Instead of trending to zero for $Q^{2}\gg m_{N}^{2}$, the
subtraction function approaches a non-zero value.

Figure 3: Local current proton subtraction function results for a range of
lattices.

In Regge analysis, it has been pointed out that the Compton amplitude may
contain an OPE-breaking ‘fixed pole’ [35, 14], which could possibly account
for the results in Fig. 3.

However, since the OPE is such a successful tool, we must investigate whether
or not this anomalous asymptotic behaviour is due to a lattice artifact. As
such artifacts will vanish in the continuum limit, we might expect them to be
sensitive to variations in the lattice volume and spacing. However, the
results in Fig. 3 indicate only very minimal volume and spacing dependence in
the subtraction function. This is in agreement with recent calculations using
baryon chiral perturbation theory, where it was found that finite volume
corrections to the Compton amplitude subtraction function would indeed be
small [36].

Since the anomalous asymptotic behaviour of the subtraction function does not
vary greatly with changes in volume and spacing, we next investigate how it
depends on the discretisation of the vector current.

##  3 Feynman-Hellmann: Conserved vector current implementation

The results shown in the previous section were based on the _local_
discretisation of the vector current, Eq. (2.2). In this section, we will
repeat the same calculation with the _conserved_ vector current,

$\begin{split}j^{\text{con}}_{\mu}(n)=\frac{1}{2}\bar{\psi}(n)\left(\vec{\mathcal{D}}_{\mu}(n,m)-\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{\mu}(n,m)\right)\psi(m).\end{split}$
(3.1)

For the Wilson fermion action, this operator is a Noether current, with a
renormalisation factor of $Z_{V}=1$, in contrast to the local operator.
However, our implementation of this current followed here introduces an
unphysical contamination to the energy shift, which we refer to as the
_seagull term_.

We implement the conserved current by introducing a perturbation on the gauge
links:

$U_{\mu}(n)\to[1+\delta_{\mu 3}(e^{i\lambda\phi(n)}-1)]U_{\mu}(n).$ (3.2)

Note that different implementations are possible; however, by modifying the
gauge links, we can use existing optimised algorithms to more efficiently
invert the fermion matrix.

Figure 4: Plot of the ratio
$R(\Upsilon_{u},t)=\mathcal{G}_{3}(\Upsilon_{u},t)/\mathcal{G}_{2}(t)$ where
the baryon three-point function was calculated with fixed insertion at $t=10$.
The purple shaded horizontal region is the fitting window ($t\in[16,20]$).

Therefore, the perturbed fermion matrix is

$M_{\lambda}(n,m)=\left(m_{0}+4\right)\delta_{n,m}-\frac{1}{2a}\sum_{\mu}\bigg{[}\vec{\mathcal{D}}_{\mu}(n,m,\lambda)+\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{\mu}(n,m,\lambda)\bigg{]},$
(3.3)

with

$\displaystyle\vec{\mathcal{D}}_{\mu}(n,m,\lambda)$
$\displaystyle\to\vec{\mathcal{D}}_{\mu}(n,m)\left[1+\delta_{\mu
3}\left(e^{i\lambda\phi(n)}-1\right)\right],$
$\displaystyle\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{\mu}(n,m,\lambda)$
$\displaystyle\to\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{\mu}(n,m)\left[1+\delta_{\mu
3}\left(e^{-i\lambda\phi(n+\hat{\mu})}-1\right)\right].$

Figure 5: Proton subtraction function results using the conserved current
implementation, before (green) and after (yellow) removal of seagull term,
$\Upsilon$. The results in blue are from the local implementation seen in Fig.
3.

The perturbed propagator is then $S_{\lambda}=[M_{\lambda}]^{-1}$, which has
the second derivative

$\frac{\partial^{2}}{\partial\lambda^{2}}{S_{{\lambda}}}\bigg{|}_{{\lambda}=0}=\underbrace{M^{-1}\big{(}\phi_{\vec{q}}\vec{\mathcal{D}}_{3}+\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{3}\phi_{\vec{q}}\big{)}M^{-1}}_{\text{seagull
contribution}}+\underbrace{2M^{-1}\big{(}\phi_{\vec{q}}\vec{\mathcal{D}}_{3}-\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{3}\phi_{\vec{q}}\big{)}M^{-1}\big{(}\phi_{\vec{q}}\vec{\mathcal{D}}_{3}-\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{3}\phi_{\vec{q}}\big{)}M^{-1}}_{\text{four-
point function}}.$ (3.4)

Therefore, the Feynman-Hellmann relation, Eq. (2.1), has an additional
_seagull term_ for this conserved current implementation:

$-m_{N}\frac{\partial^{2}E_{\lambda}}{\partial\lambda^{2}}\bigg{|}_{\lambda=0}=S_{1}(Q^{2})-\Upsilon,$
(3.5)

where the seagull term is given by

$\Upsilon_{q}=2\langle
N|\bar{\psi}_{q}(n)\big{(}\vec{\mathcal{D}}_{3}+\reflectbox{$\vec{\reflectbox{$\mathcal{D}$}}$}_{3}\big{)}\psi_{q}(n)|N\rangle.$
(3.6)

Note that
$(\phi_{\vec{q}})^{2}=2+e^{2i\vec{q}\cdot\vec{n}}+e^{-2i\vec{q}\cdot\vec{n}}$,
and the $e^{\pm 2i\vec{q}\cdot\vec{n}}$ terms do not produce a contribution to
the energy shift. Hence $\Upsilon$ is independent of the transferred momentum
$Q^{2}$.

###  3.1 Seagull results

As discussed in the preceding section, the seagull term arises from the fact
that the perturbation to the fermion matrix is not a simple linear shift, but
was instead applied to the gauge links, Eq. (3.2). Since this seagull term is
purely an artifact of the implementation process, we must remove it from our
calculation of the energy shift to isolate the subtraction function.

The three-point function corresponding to the seagull term is computed using a
conventional sequential source through the operator. The required matrix
element is then extracted from the ratio of the three-point and two-point
functions:

$R(\Upsilon_{q},t)=\frac{\mathcal{G}_{3}(\Upsilon_{q},t)}{\mathcal{G}_{2}(t)}\overset{t\gg
a}{\simeq}\Upsilon_{q},$ (3.7)

for $\Upsilon_{q}$ as defined in Eq. (3.6). The results for the up-quark are
presented in Fig. 4. Once this matrix element is calculated we remove it from
the conserved current calculation, as per Eq. (3.5), to isolate the
subtraction function.

###  3.2 Results: Conserved current implementation

The proton subtraction function is calculated for the conserved current using
the same gauge ensemble as the $32^{3}\times 64$ volume local current seen in
Tab. 1 with $N_{\text{meas}}=1000$.

In comparison to the local current, we can see in Fig. 5 that the conserved
current subtraction function exhibits different behaviour in its intermediate
$Q^{2}$ values. Once we remove the seagull term to isolate the subtraction
function, the large $Q^{2}$ behaviour closely matches OPE-breaking behaviour
of the local current. Therefore, the anomalous asymptotic behaviour of the
subtraction function can not be attributed directly to the discretisation of
the current.

##  4 Temporal Interlacing

As outlined in the previous sections, the observed anomalous asymptotic
behaviour of $S_{1}(Q^{2})$ changes very little with different discretisations
of the current or with the variations in lattice spacing and volume. Here, we
explore a Feynman-Hellmann implementation that allows us to sample the
integration region of the Compton amplitude, Eq. (1.2), more coarsely. The
motivation for considering such an implementation is that we will remove any
contamination from a potential contact term, however it comes with the caveat
that this method will also remove any essential short-distance contributions,
such as the $Z$-graph which is proposed to be responsible for the elusive
fixed pole [13]. Keeping this in mind, however, by exploring such an
implementation we hope to gain insights into whether or not the observed
anomalous behaviour is due to a short-distance effect.

As previously discussed, our calculation has a sum over all time slices on
which the currents are inserted, including contributions for which the
Euclidean separation of currents is $|z|\sim a$, with $z_{\mu}$ as in Eq.
(1.2). As a first step for investigating the effects of these contributions,
we implement a coarser sampling of the temporal integration region by
inserting two currents on different sets of time slices thereby introducing a
minimum temporal separation, $\tau_{\text{min}}$. For instance, in the
simplest case where $\tau_{\text{min}}=1$, we insert one current on the even
time slices and the other on the odd:

$[\mathcal{O}_{1}]_{n,m}=\delta^{\text{even}}_{t_{n},t_{m}}\delta_{\vec{n},\vec{m}}\phi_{\vec{q}}(\vec{n})i\gamma_{3},\quad[\mathcal{O}_{2}]_{n,m}=\delta^{\text{odd}}_{t_{n},t_{m}}\delta_{\vec{n},\vec{m}}\phi_{\vec{q}}(\vec{n})i\gamma_{3},$
(4.1)

where we have defined $\delta^{\text{even}(\text{odd})}$ to be non-vanishing
only on even (odd) timeslices. The perturbed propagator is now

$\begin{split}S_{\vec{\lambda}}(z_{n};z_{m})&=\bigg{[}\big{[}M-\lambda_{1}\mathcal{O}_{1}-\lambda_{2}\mathcal{O}_{2}\big{]}^{-1}\bigg{]}_{n,m}.\end{split}$
(4.2)

Therefore, the Feynman-Hellmann relation for $\tau_{\text{min}}=1$ interlacing
is

$\begin{split}-E_{N}\frac{\partial^{2}E_{\vec{\lambda}}}{\partial\lambda_{1}\lambda_{2}}\bigg{|}_{\vec{\lambda}=0}=\sum_{t_{1}=0,2,4,6...}\sum_{t_{2}=1,3,5,...}\sum_{\vec{z}}e^{-i\vec{q}\cdot\vec{z}}\langle
N(\vec{p})|T\\{j_{3}(\vec{z},t_{1})j_{3}(0,t_{2})\\}|N(\vec{p})\rangle.\end{split}$
(4.3)

With a judicious choice of kinematics, the RHS of Eq. (4.3) is proportional to
a discretisation of $S_{1}(Q^{2})$. The interlacing in Eq. (4.3) changes the
measure of the two sums over time-slices from $a\to 2a$, which must be
accounted for by a factor of four. However, once this normalisation is
accounted for, in the continuum limit Eq. (4.3) approaches the same object as
our previous discretisations of the Compton amplitude. Similar results can be
derived for $\tau_{\text{min}}=2$.

Figure 6: Results for multiple interlacings with the local current.
Interlacing results have low statistics, $N_{\text{meas}}\approx 200$, while
the uninterlaced ($\tau_{\text{min}}=0$) use the implementation from section
2, and have $N_{\text{meas}}\approx 1000$. Left: The proton subtraction
function, with OPE prediction from Ref. [29]. Right: The up quark contribution
to the subtracted Compton structure function, $\overline{\mathcal{F}}_{1}$,
for $Q^{2}\approx 4.7\;\text{GeV}^{2}$.

###  4.1 Interlacing results

Our preliminary results for $S_{1}(Q^{2})$ with the interlacing method are
calculated on the $32^{3}\times 64$ gauge configurations (see Tab. 1), using
the local current, Eq. (2.2). We implement interlacings for two different
values of minimum time separation: $\tau_{\text{min}}=0$, which are simply the
results from section 2, and $\tau_{\text{min}}=1$, for which the new method
has been applied. Since our interlaced results are exploratory, they have
relatively low statistics: $N_{\text{meas}}\approx 200$.

In left plot of Fig. 6, we show the subtraction function with and without
interlacing. We observe that the anomalous asymptotic behaviour of the
subtraction function is reduced in the interlaced results. We also observe
that the $Q^{2}$ dependence of the two results are remarkably similar. Hence
it would appear that by removing the $t=0$ contribution to the integral, we
have essentially removed a contribution that is constant in $Q^{2}$. Whether
this is due to an unphysical contact term or a fixed pole remains a question
to be addressed in future work. By contrast, in the right plot of Fig. 6, we
observe that the subtracted Compton structure function,
$\overline{\mathcal{F}}_{1}$ defined in Eq. (1.5), which is independent of
$S_{1}(Q^{2})$, is largely unaffected by the interlacing.

This demonstrates that the anomalous asymptotic behaviour of the subtraction
function can be attributed to very short-distance contributions (i.e. $|z|\sim
a$), which are removed by the interlacing procedure. However, the $|z|\sim a$
contributions apparently do not affect the Compton structure function
$\overline{\mathcal{F}}_{1}$ significantly. These very short-distance
contributions could be lattice artifacts as suggested by Martinelli et al.
[37, 38], or they could be of a physical origin, such as the proposed
interactions giving rise to an OPE-breaking ‘fixed pole’. However, further
investigation, numerical and analytic, is needed before we draw any strong
conclusions.

##  5 Conclusion

In this report, we present several calculations of the Compton amplitude
subtraction function, $S_{1}(Q^{2})$, in lattice QCD. In contrast to the OPE
prediction of $S_{1}(Q^{2})\sim Q^{-2}$, our initial results trend to a large
non-zero value at high-energies. This anomalous behaviour was found to persist
even after varying the lattice spacing and changing the current
discretisation.

In the final section, we present a novel method, temporal interlacing, that
allows us to more coarsely sample the integration region of the Compton
amplitude. Using this method, we demonstrate that the anomalous behaviour of
the subtraction function can be attributed to very short-distance
contributions, which may be lattice artifacts or physical contributions—future
work will aim to clarify this. We are currently performing an investigation of
the Compton amplitude and its subtraction function on configurations with
gradient flow as an extension of [39]. This investigation is ongoing and
results will be reported elsewhere.

This work provides a foundation for future lattice QCD calculations of the
Compton amplitude subtraction function. This will allow us to reduce the
theoretical uncertainties in predictions for the proton–neutron mass
difference and improve determinations of the proton charge radius.

##  6 Acknowledgements

The numerical configuration generation (using the BQCD lattice QCD program
[40])) and data analysis (using the Chroma software library [41]) was carried
out on the DiRAC Blue Gene Q and Extreme Scaling (EPCC, Edinburgh, UK) and
Data Intensive (Cambridge, UK) services, the GCS supercomputers JUQUEEN and
JUWELS (NIC, Jülich, Germany) and resources provided by HLRN (The North-German
Supercomputer Alliance), the NCI National Facility in Canberra, Australia
(supported by the Australian Commonwealth Government) and the Phoenix HPC
service (University of Adelaide). AHG is supported by an Australian Government
Research Training Program (RTP) Scholarship. RH is supported by STFC through
grant ST/P000630/1. PELR is supported in part by the STFC under contract
ST/G00062X/1. GS is supported by DFG Grant No. SCHI 179/8-1. KUC, RDY and JMZ
are supported by the Australian Research Council grants DP190100297 and
DP220103098.

## References

  * [1] J. Gasser and H. Leutwyler, _Quark Masses_ , _Phys. Rept._ 87 (1982) 77.
  * [2] G.A. Miller, B.M.K. Nefkens and I. Slaus, _Charge symmetry, quarks and mesons_ , _Phys. Rept._ 194 (1990) 1.
  * [3] W.N. Cottingham, _The neutron proton mass difference and electron scattering experiments_ , _Annals Phys._ 25 (1963) 424.
  * [4] J.C. Bernauer et al., _High-precision determination of the electric and magnetic form factors of the proton_ , _Physical Review Letters_ 105 (2010)  [1007.5076].
  * [5] R. Pohl et al., _The size of the proton_ , _Nature_ 466 (2010) 213.
  * [6] R. Pohl, R. Gilman, G.A. Miller and K. Pachucki, _Muonic hydrogen and the proton radius puzzle_ , _Ann. Rev. Nucl. Part. Sci._ 63 (2013) 175 [1301.0905].
  * [7] A. Antognini, F. Hagelstein and V. Pascalutsa, _The proton structure in and out of muonic hydrogen_ , 2205.10076.
  * [8] O. Tomalak and M. Vanderhaeghen, _Two-photon exchange correction in elastic unpolarized electron-proton scattering at small momentum transfer_ , _Phys. Rev. D_ 93 (2016) 013023 [1508.03759].
  * [9] A. Afanasev, P.G. Blunden, D. Hasell and B.A. Raue, _Two-photon exchange in elastic electron–proton scattering_ , _Prog. Part. Nucl. Phys._ 95 (2017) 245 [1703.03874].
  * [10] C.E. Carlson and M. Vanderhaeghen, _Higher order proton structure corrections to the Lamb shift in muonic hydrogen_ , _Phys. Rev. A_ 84 (2011) 020102 [1101.5965].
  * [11] R.J. Hill and G. Paz, _Model independent analysis of proton structure for hydrogenic bound states_ , _Phys. Rev. Lett._ 107 (2011) 160402 [1103.4617].
  * [12] G.A. Miller, _Proton Polarizability Contribution: Muonic Hydrogen Lamb Shift and Elastic Scattering_ , _Phys. Lett. B_ 718 (2013) 1078 [1209.4667].
  * [13] S.J. Brodsky, F.E. Close and J.F. Gunion, _A gauge-invariant scaling model of current interactions with Regge behavior and finite fixed pole sum rules_ , _Phys. Rev. D_ 8 (1973) 3678.
  * [14] S.J. Brodsky, F.J. Llanes-Estrada and A.P. Szczepaniak, _Local Two-Photon Couplings and the J=0 Fixed Pole in Real and Virtual Compton Scattering_ , _Phys. Rev. D_ 79 (2009) 033012 [0812.0395].
  * [15] D. Müller and K.M. Semenov-Tian-Shansky, _$J=0$ fixed pole and $D$-term form factor in deeply virtual Compton scattering_, _Phys. Rev. D_ 92 (2015) 074025 [1507.02164].
  * [16] D. Drechsel, B. Pasquini and M. Vanderhaeghen, _Dispersion relations in real and virtual Compton scattering_ , _Phys. Rept._ 378 (2003) 99 [hep-ph/0212124].
  * [17] J.M. Alarcon, V. Lensky and V. Pascalutsa, _Chiral perturbation theory of muonic hydrogen Lamb shift: polarizability contribution_ , _Eur. Phys. J. C_ 74 (2014) 2852 [1312.1219].
  * [18] I. Caprini, _Constraints on the virtual Compton scattering on the nucleon in a new dispersive formalism_ , _Phys. Rev. D_ 93 (2016) 076002 [1601.02787].
  * [19] I. Caprini, _Analyticity and Regge asymptotics in virtual Compton scattering on the nucleon_ , _Eur. Phys. J. C_ 81 (2021) 309 [2102.08601].
  * [20] F.B. Erben, P.E. Shanahan, A.W. Thomas and R.D. Young, _Dispersive estimate of the electromagnetic charge symmetry violation in the octet baryon masses_ , _Phys. Rev. C_ 90 (2014) 065205 [1408.6628].
  * [21] J. Gasser and H. Leutwyler, _Implications of Scaling for the Proton - Neutron Mass - Difference_ , _Nucl. Phys. B_ 94 (1975) 269.
  * [22] A.W. Thomas, X.G. Wang and R.D. Young, _Electromagnetic Contribution to the Proton-Neutron Mass Splitting_ , _Phys. Rev. C_ 91 (2015) 015209 [1406.4579].
  * [23] O. Tomalak and M. Vanderhaeghen, _Two-photon exchange correction in elastic unpolarized electron-proton scattering at small momentum transfer_ , _Phys. Rev. D_ 93 (2016) 013023 [1508.03759].
  * [24] O. Tomalak, _Electromagnetic proton–neutron mass difference_ , _Eur. Phys. J. Plus_ 135 (2020) 411 [1810.02502].
  * [25] A. Walker-Loud, C.E. Carlson and G.A. Miller, _The Electromagnetic Self-Energy Contribution to $M_{p}-M_{n}$ and the Isovector Nucleon MagneticPolarizability_, _Phys. Rev. Lett._ 108 (2012) 232301 [1203.0254].
  * [26] A. Walker-Loud, _On the Cottingham formula and the electromagnetic contribution to the proton-neutron mass splitting_ , _PoS_ CD2018 (2019) 045 [1907.05459].
  * [27] J. Gasser, H. Leutwyler and A. Rusetsky, _Sum rule for the Compton amplitude and implications for the proton–neutron mass difference_ , _Eur. Phys. J. C_ 80 (2020) 1121 [2008.05806].
  * [28] J.C. Collins, _Renormalization of the Cottingham Formula_ , _Nucl. Phys. B_ 149 (1979) 90.
  * [29] R.J. Hill and G. Paz, _Nucleon spin-averaged forward virtual compton tensor at large $q^{2}$_, _Physical Review D_ 95 (2017)  [1611.09917].
  * [30] A. Chambers et al., _Nucleon structure functions from operator product expansion on the lattice_ , _Physical Review Letters_ 118 (2017)  [1703.01153].
  * [31] K.U. Can et al., _Lattice QCD evaluation of the Compton amplitude employing the Feynman-Hellmann theorem_ , _Phys. Rev. D_ 102 (2020) 114505 [2007.01523].
  * [32] A. Hannaford-Gunn, K.U. Can, R. Horsley, Y. Nakamura, H. Perlt, P.E.L. Rakow et al., _Generalized parton distributions from the off-forward Compton amplitude in lattice QCD_ , _Phys. Rev. D_ 105 (2022) 014502 [2110.11532].
  * [33] B. Sheikholeslami and R. Wohlert, _Improved Continuum Limit Lattice Action for QCD with Wilson Fermions_ , _Nucl. Phys._ B259 (1985) 572.
  * [34] W. Bietenholz et al., _Flavour blindness and patterns of flavour symmetry breaking in lattice simulations of up, down and strange quarks_ , _Phys. Rev. D_ 84 (2011) 054509 [1102.5300].
  * [35] M.J. Creutz, S.D. Drell and E.A. Paschos, _High-Energy Limit for the Real Part of Forward Compton Scattering_ , _Phys. Rev._ 178 (1969) 2300.
  * [36] J. Lozano, A. Agadjanov, J. Gegelia, U.G. Meißner and A. Rusetsky, _Finite volume corrections to forward Compton scattering off the nucleon_ , _Phys. Rev. D_ 103 (2021) 034507 [2010.10917].
  * [37] C. Dawson, G. Martinelli, G.C. Rossi, C.T. Sachrajda, S.R. Sharpe, M. Talevi et al., _New lattice approaches to the delta I = 1/2 rule_ , _Nucl. Phys. B_ 514 (1998) 313 [hep-lat/9707009].
  * [38] G. Martinelli, _Hadronic weak interactions of light quarks_ , _Nucl. Phys. B Proc. Suppl._ 73 (1999) 58 [hep-lat/9810013].
  * [39] K.U. Can et al., _Advances in lattice hadron physics calculations using the gradient flow_ , in _38th International Symposium on Lattice Field Theory_ , 12, 2021 [2112.06355].
  * [40] T.R. Haar, Y. Nakamura and H. Stüben, _An update on the BQCD Hybrid Monte Carlo program_ , _EPJ Web Conf._ 175 (2018) 14011 [1711.03836].
  * [41] R.G. Edwards and B. Joo, _The Chroma software system for lattice QCD_ , _Nucl.Phys.Proc.Suppl._ 140 (2005) 832 [hep-lat/0409003].