[a]Lorenzo Barca

# Progress on nucleon transition matrix elements with a lattice QCD
variational analysis

Gunnar S. Bali  Sara Collins

###### Abstract

Nucleon weak matrix elements can be extracted from nucleon correlation
functions with lattice QCD simulations. The signal-to-noise ratio prohibits
the analysis at large source-sink separations and as a consequence, excited
state contamination affects the extraction of the nucleon matrix elements.
Chiral perturbation theory (ChPT) suggests that the dominant contamination in
some of these channels is due to $N\pi$ states where the pion carries the same
momentum of the current. In this talk, we report updates on the variational
analysis with $qqq$-operators (nucleon-like) and $(qqq)(\bar{q}q)$-operators
(nucleon-pion-like) where we report for the first time some preliminary
results of $\langle N\pi|\mathcal{J}|N\rangle$, modulo some kinematic and
volume factors, and we compare the results against ChPT. This pilot study is
performed on a CLS ensemble with $N_{f}=3$, $m_{\pi}\approx
420~{}\mathrm{MeV}$, $a\approx 0.1~{}\mathrm{fm}$ and $T=2L\approx
4.8~{}\mathrm{fm}$.

##  1 Introduction

The next-generation neutrino oscillation experiments DUNE [1, 2] and Hyper-
Kamiokande [3] will address fundamental questions in particle physics and
cosmology, including whether neutrinos oscillate differently from
antineutrinos. CP violation in the lepton sector [4] is necessary for a
possible Standard Model explanation of the matter-antimatter asymmetry in our
universe [5]. Other prospects for physics searches at DUNE include more
precise studies of the neutrino mass hierarchy, supernova neutrino bursts [6],
and BSM physics [7]. Of fundamental importance are the parameters that
describe the cross sections of neutrino-nucleus scattering data, as very high
precision is required in the analysis of such rare weak processes. Nuclear
models relate the neutrino-nucleus scattering cross-section, whose events
occur at the detector level, to the neutrino-nucleon scattering cross-section,
which is parametrised in terms of nucleon form factors [8, 9].

MiniBooNE, a previous terrestrial neutrino oscillation experiment, observed a
significant production of nucleon-pion final states (CC-1$\pi$) in addition to
quasielastic charged current processes (CCQE) [10]. These additional states
are produced indirectly via resonance states, such as
$\nu_{\mu}+n\to\mu^{-}+\Delta^{+}(\to\pi^{+}+n)$ in the intermediate energy
range $E_{\nu}\approx 0.1-20~{}\mathrm{GeV}$ [11]. Achieving a few-percent
accuracy in determining (anti)neutrino oscillation parameters requires precise
knowledge of the electroweak transitions $N\to\Delta$, $N\to N^{*}$ and $N\to
N\pi$ in the non-perturbative regime [12].

First principles calculations based on lattice QCD simulations can provide the
form factors that enter in the cross-sections [13, 14]. Nucleon form factors
are extracted from ratios of nucleon three- and two-point functions computed
using Monte Carlo techniques. A major challenge of lattice nucleon structure
calculations is the problem of excited state contamination. In principle, this
can be overcome by evaluating nucleon correlation functions at large source-
sink separations, as the excited state contamination decays exponentially with
the Euclidean time separation. Unfortunately, the signal-to-noise ratio
deteriorates with the distance between nucleon operators, making it
impractical to analyse the data at long distances. A complete variational
analysis which involves constructing the nucleon spectrum using the Lüscher
method [15, 16] and then using the GEVP solutions to diagonalise the
interpolating operators [17, 18], can potentially solve this problem. This
more robust approach requires an analysis with both single- ($N$) and multi-
hadron ($N\pi$, $N\pi\pi$, …) operators, making lattice simulations demanding
due to the increased number of Wick contractions and Dirac propagators.

Chiral perturbation theory (ChPT) can guide us in understanding which states
contribute the most to the contamination of the nucleon three-point functions
and, consequently, which interpolating operators are most relevant for the
variational analysis. A reliable estimation of the excited state contamination
and other systematic effects is essential for the accurate extraction of the
$N\to N$ electroweak form factors. At leading order in ChPT, it was found that
$N\pi$ states in the kinematic configuration where the pion carries the same
momentum as the current represent the dominant contribution [19, 20]. Other
$N\pi$ states as well as $N\pi\pi$ and nucleon excitations will also be
relevant to reach few percent accuracy. Multistate fits informed by ChPT have
proven to be successful in determining form factors that satisfy the Pion Pole
Dominance (PPD) and the generalised Goldberger-Treiman (also known as PCAC)
relations [21, 22].

In this study, we carry out a variational analysis utilising just a single
$qqq$ (three-quark or $N$-like) interpolating operator and a single
$(qqq)(\bar{q}q)$ (five-quark or $N\pi$-like) interpolating operator.
Specifically, we consider the $N\pi$-like operator where the pion carries the
same momentum as the current, aiming to verify the ChPT prediction regarding
its large contribution. We compute nucleon matrix elements $\langle
N|\mathcal{J}|N\rangle$, previously presented in [23] and, for the first time,
we present preliminary results for positive-parity transition matrix elements
$\langle N\pi|\mathcal{J}|N\rangle$ in a finite volume, up to some kinematic
factors. Here, $|N\rangle$ denotes the nucleon ground state and $|N\pi\rangle$
a physical nucleon-pion finite volume state. This work provides an update of
the results of [24, 23]. In sec. 2, we provide an overview of the extraction
of nucleon form factors from nucleon matrix elements using the standard ratio
method, along with a discussion of key results from ChPT in these channels. In
sec. 3, we discuss how we use the variational analysis to determine the matrix
elements by projecting the operators to the desired eigenstate. In sec. 4, we
present preliminary results for $\langle N\pi|\mathcal{J}|N\rangle$ at zero
spatial momentum transfer. Finally, in sec. 5, we draw conclusions and discuss
key observations of this analysis.

##  2 Nucleon matrix elements

###  2.1 Ratio method

The nucleon matrix elements $\langle
N(\mathbf{p}^{\prime})|\mathcal{J}(\mathbf{q})|N(\mathbf{p})\rangle$ appear in
the spectral decomposition of the nucleon three-point correlation functions

$C^{k,\mathcal{J}}_{3pt}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)~{}=~{}\mathbb{P}^{+}_{k}~{}\langle\mathrm{O}_{3q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle~{},$
(1)

as it can be seen by inserting two complete sets of states. In this
expression, $\mathrm{O}_{3q}$ is a $3$-quark operator that has the nucleon
quantum numbers $J^{PC}=\frac{1}{2}^{++}$ and $\mathcal{J}$ is the local
current projected to momentum $\mathbf{q}=\mathbf{p}^{\prime}-\mathbf{p}$. The
matrix $\mathbb{P}^{+}_{k}$ projects the operators to positive parity and it
aligns the spin along the direction $k$. In eq. (1), $t$ is the source-sink
separation between the two nucleon operators and $\tau$ is the intermediate
time where the current is located. Additionally, one computes the nucleon two-
point functions

$C_{2pt}(\mathbf{p},t)~{}=~{}\mathbb{P}^{+}~{}\langle\mathrm{O}_{3q}(\mathbf{p},t)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle~{},$
(2)

where $\mathbb{P}^{+}$ is a positive-parity projector and the spectral
decomposition reads explicitly

$C_{2pt}(\mathbf{p},t)~{}=~{}\frac{|\langle\Omega|\mathrm{O}^{+}_{3q}|N(\mathbf{p})\rangle|^{2}}{2E_{N}V}e^{-E_{N}t}~{}+~{}\sum_{\mathbf{k}}\frac{|\langle\Omega|\mathrm{O}^{+}_{3q}|N(\mathbf{p}+\mathbf{k})\pi(-\mathbf{k})\rangle|^{2}}{(2E_{N}V)(2E_{\pi}V)}e^{-E_{N\pi}t}~{}+~{}\dots$
(3)

where $\mathrm{O}^{+}_{3q}=\mathbb{P}^{+}\mathrm{O}_{3q}$ and $|\Omega\rangle$
is the vacuum state. In this expression, to simplify the notation we include
only the ground state $N$ and the tower of $N\pi$ states with total momentum
$\mathbf{p}$, and in the ellipses we include all the other states in the
spectrum. Using nucleon two- and three-point functions, we construct the
following ratio

$R^{k,\mathcal{J}}(\mathbf{p}^{\prime},t;\mathbf{q},\tau):=~{}\frac{C_{3pt}^{k,\mathcal{J}}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)}{C_{2pt}(\mathbf{p}^{\prime},t)}\
\sqrt{\frac{C_{2pt}(\mathbf{p}^{\prime},\tau)~{}C_{2pt}(\mathbf{p}^{\prime},t)~{}C_{2pt}(\mathbf{p},t-\tau)}{C_{2pt}(\mathbf{p},\tau)~{}C_{2pt}(\mathbf{p},t)~{}C_{2pt}(\mathbf{p}^{\prime},t-\tau)}}~{}.$
(4)

By employing the spectral decomposition of the nucleon three- and two-point
functions, it is straightfordward to show that the nucleon ground state
contribution to eq. (4) is independent of $t$ and $\tau$, and proportional to
$\langle N(\mathbf{p}^{\prime})|\mathcal{J}(\mathbf{q})|N(\mathbf{p})\rangle$.
The latter can be decomposed in terms of nucleon form factors through Lorentz
decomposition, which depends on the intermediate current. For pseudoscalar and
axial-vector currents, the Lorentz decompositions in Euclidean space-time are

$\displaystyle\langle
N(\mathbf{p}^{\prime})|{\mathcal{P}(\mathbf{q})}|N(\mathbf{p})\rangle=\bar{u}_{N}(\mathbf{p}^{\prime})~{}\left[\gamma_{5}G_{P}(Q^{2})\right]~{}u_{N}(\mathbf{p})~{},$
(5) $\displaystyle\langle
N(\mathbf{p}^{\prime})|{\mathcal{A}_{\mu}(\mathbf{q})}|N(\mathbf{p})\rangle=\bar{u}_{N}(\mathbf{p}^{\prime})~{}\left[\gamma_{\mu}\gamma_{5}G_{A}(Q^{2})+\frac{q_{\mu}}{2m_{N}}\gamma_{5}\widetilde{G}_{P}(Q^{2})\right]~{}u_{N}(\mathbf{p})~{},$
(6)

where $q_{\mu}=p^{\prime}_{\mu}-p_{\mu}$ and $Q^{2}=-q^{\mu}q_{\mu}$ is the
four-momentum transfer. In eq. (5), $G_{P}$ is the pseudoscalar form factor,
while in eq. (6), $G_{A}$ and $\widetilde{G}_{P}$ are the axial and induced
pseudoscalar form factors, respectively. If we write the Lorentz decomposition
like $\langle
N(\mathbf{p}^{\prime})|{\mathcal{J}(\mathbf{q})}|N(\mathbf{p})\rangle=\bar{u}_{N}(\mathbf{p}^{\prime})\left[FF_{N\mathcal{J}N}\right]u_{N}(\mathbf{p})$,
the spectral decomposition of the ratio in eq. (4) yields

$R^{k,\mathcal{J}}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)=\sqrt{\frac{(E^{\prime}_{N}V)(E_{N}V)}{(E^{\prime}_{N}+m_{N})(E_{N}+m_{N})}}\frac{\mathrm{Tr}\left\\{\mathbb{P}^{+}_{k}(\not{p}^{\prime}+m_{N})\left[FF_{N\mathcal{J}N}\right](\not{p}+m_{N})\right\\}}{(2E^{\prime}_{N}V)(2E_{N}V)}~{}+~{}\dots,$
(7)

where we are considering only the nucleon ground state, while the
contributions from excited states, which decay exponentially with time, are
included in the ellipses. The nucleon form factors in
$\left[FF_{N\mathcal{J}N}\right]$ can be extracted from the ratio in eq. (4)
at large $t$ and $\tau$, as the contribution from excited states is reduced.
However, the correlation functions in eqs. (1)–(2) are estimated with Monte
Carlo simulations and the error associated with the measurement decreases more
slowly than the signal. Thus, the degradation of the signal-to-noise ratio
prohibits the analysis at large source-sink separations. In principle, the
error can be systematically reduced as it scales with
$\mathcal{O}(1/\sqrt{N})$, where $N$ is the number of gauge configurations or
statistical samples. However, it is impractical with the current computers and
algorithms to compute nucleon correlation functions with
$N\gg\mathcal{O}(10^{4})$, and thus the statistical sample size is usually
limited to $10^{3}\lesssim N\lesssim 10^{4}$. Due to this constraint, the
extraction of nucleon form factors from the ratio in eq. (4) is affected by
excited state contamination (ESC). To resolve the nucleon ground state, one
needs a robust method that takes into account all the relevant excited states.

###  2.2 $N\pi$ states contamination in nucleon correlation functions from
ChPT

ChPT can provide guidance on the contribution of specific states to the
correlation functions. Since ChPT is an expansion for small pion masses and
momenta, its predictions are most reliable towards physical pion masses and
small momenta. In [25], it was estimated with leading order in ChPT (LO-ChPT)
that the contamination of the tower of $N(\mathbf{p})\pi(-\mathbf{p})$ states
to the zero-momentum nucleon two-point function is at the level of a few
percent for close to physical pion masses ($\approx 200~{}\mathrm{MeV}$) and
source-sink separations of $t\approx 0.5~{}\mathrm{fm}$. The contribution of
$N\pi\pi$ states to the nucleon two-point functions is estimated to be even
smaller at the per mille level [26]. Multiparticle states are usually
suppressed by a volume factor for each additional state. While LO-ChPT
predicts that $N\pi$ states do not contribute much to the nucleon two-point
functions at sufficently large distance, their contribution to the nucleon
three-point functions can be very large. In [19, 20], it was predicted with
LO-ChPT that the $N\pi$ states contribution depends on the momentum transfer
$Q^{2}$ and it can rise up to $40\%$ for physical pion masses at
$t=2~{}\mathrm{fm}$. Most of this contribution is due to $N\pi$ states which
are either created at the source or at the sink and where the pion has the
same momentum as the current. These current-enhanced diagrams are not volume
suppressed. More explicitly, in lattice simulations the nucleon operators at
the sink are usually projected to $\mathbf{p}^{\prime}=\mathbf{0}$, and thus
the momentum of the source operator is $\mathbf{p}=-\mathbf{q}$ by momentum
conservation. With these kinematics, the dominant $N\pi$ states at the sink
have total momentum $\mathbf{p}^{\prime}=\mathbf{0}$ and are in the momentum
configuration $N(-\mathbf{q})\pi(\mathbf{q})$, while at the source they are in
the momentum configuration $N(\mathbf{0})\pi(-\mathbf{q})$. In the forward
limit ($\mathbf{q}=\mathbf{0}$), ChPT predicts still a large $N\pi$ states
contamination in specific channels when the sink and source nucleon operators
have non-zero momentum $\mathbf{p}^{\prime}=\mathbf{p}\neq\mathbf{0}$, see
[21]. This was confirmed numerically in [23]. The remaining tower of
interacting $N\pi$ states with other relative momenta and
$N(\mathbf{0})\pi(\mathbf{0})\pi(\mathbf{0})$ states can further contribute to
contaminate the nucleon three-point correlation functions up to few percent at
the physical point, and they must be taken into account for high-precision
determination of the nucleon form factors. Resonance contributions in a finite
volume like the Roper in the positive-parity nucleon channel must also be
considered and effective field theories provide guidance on the interacting
energy levels [27].

##  3 Variational analysis

The variational analysis consists of constructing a basis of $n$ operators
$\mathrm{O}_{i}$ that have the quantum numbers of the nucleon and which are
used to compute a matrix of two-point correlation functions

$\mathbb{C}_{2pt}(\mathbf{p},t)_{ij}=~{}\mathbb{P}^{+}~{}\langle\mathrm{O}_{i}(\mathbf{p},t)~{}\bar{\mathrm{O}}_{j}(\mathbf{p},0)\rangle~{}.$
(8)

This matrix is then employed to solve the Generalised EigenValue Problem
(GEVP)

$\mathbb{C}_{2pt}(\mathbf{p},t)V(\mathbf{p},t;t_{0})=\mathbb{C}_{2pt}(\mathbf{p},t_{0})\Lambda(\mathbf{p},t;t_{0})V(\mathbf{p},t;t_{0})~{}$
(9)

to find the matrices of eigenvalues
$\Lambda(\mathbf{p},t;t_{0})=\mathrm{diag}\left(\lambda^{1}(\mathbf{p},t;t_{0}),...,\lambda^{n}(\mathbf{p},t;t_{0})\right)$
and eigenvectors
$V(\mathbf{p},t;t_{0})=\left(\mathbf{v}^{1}(\mathbf{p},t;t_{0}),\dots,\mathbf{v}^{n}(\mathbf{p},t;t_{0})\right)$.
Each eigenvalue decays exponentially with the energy of the specific state in
the spectrum $\lambda^{\alpha}\propto e^{-E_{\alpha}(t-t_{0})}$. The
$n$-dimensional eigenvectors are orthogonal with respect to the scalar product

$\left(\tilde{v}(\mathbf{p},t;t_{0}),\mathbb{C}_{2pt}(\mathbf{p},t_{0})\tilde{v}(\mathbf{p},t;t_{0})\right)=\sum_{i,j}^{n}\tilde{v}_{i}^{*\alpha}(\mathbf{p},t;t_{0})~{}\mathbb{C}_{2pt}(\mathbf{p},t_{0})_{ij}~{}\tilde{v}_{j}^{\beta}(\mathbf{p},t;t_{0})=\delta^{\alpha\beta}~{},$
(10)

where we have normalized
$\tilde{v}^{\alpha}=\left(v^{\alpha},\mathbb{C}_{2pt}v^{\alpha}\right)^{-1/2}v^{\alpha}$.
Since our basis is finite and cannot span the infinite dimensional Hilbert
space, the $\alpha$-th eigensolutions are correct up to
$\mathcal{O}(e^{-(E_{n+1}-E_{\alpha})t})$, where $(n+1)$ labels the first
state that is missing in our finite basis. The GEVP eigenvectors can be used
to diagonalise the correlation matrix and to construct an operator that has
overlap predominantly with a desired state $\alpha$ in the spectrum. To
construct the GEVP-improved operator the following linear combination of
eigenvectors and operators must be considered:

$\mathrm{O}^{\alpha}(\mathbf{p},t)=\sum_{i}^{n}v_{i}^{\alpha}(\mathbf{p},t_{\alpha};t_{0})\mathrm{O}_{i}(\mathbf{p},t)~{}.$
(11)

With a sufficiently large (and orthogonal) basis the system will be completely
diagonalised and the eigenvectors will be time independent already from
$t\approx 0.0~{}\mathrm{fm}$. However, due to the finite size of the operator
basis, the eigenvectors are usually not time independent as remaining
excitations are not completely removed at small distances. The linear
combination is then constructed using eigenvectors at a reference time
$t_{\alpha}$ such that a plateau in the effective masses of the eigenvalues
sets in with the energy of the state $\alpha$. After the eigenstate projection
in eq. (11), we expect that

$\bar{\mathrm{O}}^{\alpha}|\Omega\rangle\approx Z^{\alpha}|\alpha\rangle~{},$
(12)

where the approximation is due to the finite size of our operator basis. In
particular, the GEVP-improved operators are then used to compute GEVP-improved
two-point functions, which read

$C_{2pt}(\mathbf{p},t)^{\alpha}~{}=~{}\langle\mathrm{O}^{\alpha}(\mathbf{p},t)~{}\bar{\mathrm{O}}^{\alpha}(\mathbf{p},0)\rangle=\sum_{i,j}~{}v_{i}^{\alpha}(\mathbf{p},t_{\alpha};t_{0})\langle\mathrm{O}_{i}(\mathbf{p},t)~{}\bar{\mathrm{O}}_{j}(\mathbf{p},0)\rangle~{}v_{j}^{\alpha}(\mathbf{p},t_{\alpha};t_{0})~{}.$
(13)

Similarly, one can adopt the eigenstate-projected operators to compute the
GEVP-improved three-point functions

$C^{k,\mathcal{J}}_{3pt}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)^{\alpha\beta}~{}=\mathbb{P}^{+}_{k}~{}\langle\mathrm{O}^{\alpha}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}^{\beta}(\mathbf{p},0)\rangle~{}.$
(14)

In order for the GEVP to resolve a finite number of $n$ states in the nucleon
spectrum, the operator basis must be as orthogonal as possible, i.e., the
operators must be constructed to have major overlap with different states in
the nucleon spectrum. To extract nucleon matrix elements $\langle
N^{\prime}|\mathcal{J}|N\rangle$ with the variational method, we project the
correlation functions to the nucleon eigenstate, i.e. $\alpha=\beta=N$ in eqs.
(13)-(14), and we construct the ratio

$R^{k,\mathcal{J}}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)^{{NN}}=~{}\frac{C_{3pt}^{j,\mathcal{J}}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)^{{NN}}}{C_{2pt}(\mathbf{p}^{\prime},t)^{N}}\
\sqrt{\frac{C_{2pt}(\mathbf{p}^{\prime},\tau)^{N}~{}C_{2pt}(\mathbf{p}^{\prime},t)^{N}~{}C_{2pt}(\mathbf{p},t-\tau)^{N}}{C_{2pt}(\mathbf{p},\tau)^{N}~{}C_{2pt}(\mathbf{p},t)^{N}~{}C_{2pt}(\mathbf{p}^{\prime},t-\tau)^{N}}}~{}.$
(15)

This expression is the same as in eq. (4), except for the fact that the
correlation functions are now projected to the nucleon state ($N$). Therefore,
the exponential corrections from higher energy states included in the ellipses
of eq. (7) are removed, provided that the basis is good enough to resolve the
$n$ states in the spectrum. Similarly, to extract transition matrix elements
$\langle N\pi|\mathcal{J}|N\rangle$, we project the correlation functions both
to $N$ and $N\pi$ and construct the ratio

$R^{k,\mathcal{J}}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)^{{N\pi
N}}=~{}\frac{C_{3pt}^{j,\mathcal{J}}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)^{{N\pi
N}}}{C_{2pt}(\mathbf{p}^{\prime},t)^{N\pi}}\
\sqrt{\frac{C_{2pt}(\mathbf{p}^{\prime},\tau)^{N\pi}~{}C_{2pt}(\mathbf{p}^{\prime},t)^{N\pi}~{}C_{2pt}(\mathbf{p},t-\tau)^{N}}{C_{2pt}(\mathbf{p},\tau)^{N}~{}C_{2pt}(\mathbf{p},t)^{N}~{}C_{2pt}(\mathbf{p}^{\prime},t-\tau)^{N\pi}}}~{}.$
(16)

A similar ratio without GEVP-improved correlators and with
$\mathbf{p}^{\prime}=\mathbf{p}=\mathbf{q}=\mathbf{0}$ was employed in [28] to
determine the $N\pi$ contribution in electric polarizabilities
($N\gamma^{*}\to N\pi$). The three-point function in this expression is
projected at the source to the $N$ eigenstate and at the sink to the $N\pi$
eigenstate. By applying the spectral decomposition to the correlation
functions in eq. (16) with the normalisation convention defined in eq. (3), we
obtain the expression

$\displaystyle R^{k,\mathcal{J}}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)^{{N\pi
N}}=$
$\displaystyle~{}\sqrt{\frac{(2E^{\prime}_{N}V)(2E^{\prime}_{\pi}V)(2E_{N}V)}{2(E_{N}+m_{N})\sum_{\sigma}u^{\sigma}_{N\pi}(\mathbf{p}^{\prime})\bar{u}^{\sigma}_{N\pi}(\mathbf{p}^{\prime})}}~{}\times~{}$
$\displaystyle\times$
$\displaystyle~{}\frac{\mathrm{Tr}\left[\mathbb{P}^{+}_{k}\sum_{\sigma}u^{\sigma}_{N\pi}(\mathbf{p}^{\prime})\bar{u}^{\sigma}_{N\pi}(\mathbf{p}^{\prime})\left[FF_{N\pi\mathcal{J}N}\right](\not{p}+m_{N})\right]}{(2E^{\prime}_{N}V)(2E_{N}V)}~{}$
(17)

where $\left[FF_{N\pi\mathcal{J}N}\right]$ is the decomposition of $\langle
N\pi|\mathcal{J}|N\rangle$ in terms of form factors as in eqs. (5)-(6),

$\langle
N\pi|\mathcal{J}|N\rangle=u_{N\pi}\left[FF_{N\pi\mathcal{J}N}\right]u_{N}$
(18)

In principle, one should include all the relevant operators to reconstruct the
nucleon spectrum ($N$, $N^{*}$, $N\pi$, $N\pi\pi$, $\Delta\pi$, …) and use the
Lüscher method [29] to compute the phase shift of $N\pi$ scattering in the
Roper channel from the energy levels in a finite volume. In particular, multi-
hadron operators must be constructed with relative momenta up to some value,
above which their effect is expected to be negligible. This would give more
confidence that no relevant operator is missing in the variational basis and
that the eigenvectors and eigenvalues are correct up to negligible exponential
corrections. Only then, the GEVP eigenvectors can be used to construct GEVP-
projected two-point and three-point functions. An example of systematic
effects in the reconstruction of the spectra due to missing relevant operators
in the GEVP basis can be seen for example in Fig. 6 of [30] in the
$\rho$-resonance channel.

While this approach of reconstructing the spectrum is very sound, current
algorithms make it impractical to realise. The difficulty lies in the
computation of many all-to-all Dirac propagators that are needed to evaluate
two-point and three-point functions for two- and three-particle operators. In
addition, the computational expense of computing a symmetric matrix of two-
point or three-point functions increases with $n\times(n+1)/2$, with $n$ being
the number of operators in the basis. Therefore, efficient algorithms to
compute multi-hadron correlation functions are required. In this work, we use
the sequential method [31] along with stochastic methods [32] to compute the
Dirac propagators. Unfortunately, the former method is not the most efficient
way to compute these multi-hadron correlation functions as many sequential
propagators are required to take into account for instance the different
momentum combinations and polarizations. There are promising and versatile
alternatives to the methods used in this work like distillation [33],
stochastic distillation [34], and sparsening methods. The latter is employed
in [35] to include $N\pi$ and $N\pi\pi$ operators in the variational analysis.

Although the final goal must be to include all the relevant operators in the
variational basis, investigations of the variational analysis with the most
relevant operators can still provide quite accurate results on the final
quantities of interest. As we discussed in the previous section, ChPT can be
used to predict which are the most relevant operators for a given specific
channel. In particular, LO-ChPT predicts that $N\pi$ states play an important
role in the determination of pseudoscalar and axial-vector form factors from
nucleon three-point functions at source-sink separations that are can
currently be reached ($\approx 1.5~{}\mathrm{fm}$). Motivated by these
different studies, we construct a basis with just two operators that have very
large overlap with $N$ and $N\pi$ states, respectively. We then compute a
$2\times 2$ matrix of two-point functions to find the GEVP eigenvalues and
eigenvectors. The latter are then used to project the interpolating operators
to the GEVP eigenstates $N$ and $N\pi$, and to extract $\langle
N|\mathcal{J}|N\rangle$, $\langle N\pi|\mathcal{J}|N\rangle$, by means of eq.
(15) or (16), respectively.

##  4 Lattice simulations

###  4.1 Details of the simulation and construction of operators

We simulate QCD with $N_{f}=3$ degenerate quarks using an
$\mathcal{O}(a)$-improved Wilson-Clover action. The measurements are performed
on 800 gauge configurations of the CLS111See https://wiki-
zeuthen.desy.de/CLS/CLS for more information. ensemble A653, with
$m_{\pi}=420$ MeV, $N_{t}=2N_{s}=48a$ and lattice spacing $a=0.098$ fm, see
[36]. We investigate processes $\langle
I^{\prime}_{z};J^{PC}|\mathcal{J}^{-}|I_{z};J^{PC})\rangle$ with
$I_{z}=-\frac{1}{2}=-I^{\prime}_{z}$, and $J^{PC}=\frac{1}{2}^{++}$, i.e.
lattice irrep $G_{1}^{(g)}$. To this end, we construct two operators
$\mathrm{O}_{3q}$ and $\mathrm{O}_{5q}$ to have nucleon quantum numbers
($J^{PC}=\frac{1}{2}^{++}$). The three-quark operator is the standard nucleon
operator
$\mathrm{O}_{3q}^{\gamma}=\epsilon_{abc}\left[u^{\top}C\gamma_{5}d\right]q_{\gamma}$,
with $q=u,d$ for proton or neutron operators, while the five-quark operator is
a nucleon-pion-like operator
$\mathrm{O}_{5q}\sim\mathrm{O}_{3q}\mathrm{O}_{\bar{q}q}$ with
$\mathrm{O}_{\bar{q}q}\sim\bar{q}\gamma_{5}q$ . For the analysis, we construct
operators in the rest frame ($\mathbf{p}=\mathbf{0}$) and in six equivalent
moving frames with discretised unit lattice momentum
$\mathbf{p}=\mathbf{e_{i}}:=\pm\frac{2\pi}{L}\hat{n}_{i}$, where $\hat{n}_{i}$
is a unit vector along the direction $i=x,y,z$. The two-particle operator is
projected to zero momentum by constructing the composite operator according to
the lattice irrep $G_{1}^{+}$ of ${}^{2}O_{h}$, the double cover of the cubic
group $O_{h}$, see [37, 38]. The projection to unit lattice momentum is done
by constructing the two-particle operator according to the lattice irrep
$G_{1}$ of the little group $C_{4v}$, see [39]. Here is a summary of these
constructions:

##### Rest frame

$\displaystyle\mathrm{O}_{5q}^{G_{1}^{+},\uparrow}(\mathbf{p}=\mathbf{0})$
$\displaystyle=\mathrm{O}^{\downarrow}_{3q}(-\mathbf{e_{x}})\mathrm{O}_{\bar{q}q}(\mathbf{e_{x}})-\mathrm{O}^{\downarrow}_{3q}(\mathbf{e_{x}})\mathrm{O}_{\bar{q}q}(-\mathbf{e_{x}})-i\mathrm{O}^{\downarrow}_{3q}(-\mathbf{e_{y}})\mathrm{O}_{\bar{q}q}(\mathbf{e_{y}})~{}+$
$\displaystyle+i\mathrm{O}^{\downarrow}_{3q}(\mathbf{e_{y}})\mathrm{O}_{\bar{q}q}(-\mathbf{e_{y}})+\mathrm{O}^{\uparrow}_{3q}(-\mathbf{e_{z}})\mathrm{O}_{\bar{q}q}(\mathbf{e_{z}})-\mathrm{O}^{\uparrow}_{3q}(\mathbf{e_{z}})\mathrm{O}_{\bar{q}q}(-\mathbf{e_{z}})~{},$
(19)
$\displaystyle\mathrm{O}_{5q}^{G_{1}^{+},\downarrow}(\mathbf{p}=\mathbf{0})$
$\displaystyle=\mathrm{O}^{\uparrow}_{3q}(-\mathbf{e_{x}})\mathrm{O}_{\bar{q}q}(\mathbf{e_{x}})-\mathrm{O}^{\uparrow}_{3q}(\mathbf{e_{x}})\mathrm{O}_{\bar{q}q}(-\mathbf{e_{x}})+i\mathrm{O}^{\uparrow}_{3q}(-\mathbf{e_{y}})\mathrm{O}_{\bar{q}q}(\mathbf{e_{y}})~{}+$
$\displaystyle-i\mathrm{O}^{\uparrow}_{3q}(\mathbf{e_{y}})\mathrm{O}_{\bar{q}q}(-\mathbf{e_{y}})-\mathrm{O}^{\downarrow}_{3q}(-\mathbf{e_{z}})\mathrm{O}_{\bar{q}q}(\mathbf{e_{z}})+\mathrm{O}^{\downarrow}_{3q}(\mathbf{e_{z}})\mathrm{O}_{\bar{q}q}(-\mathbf{e_{z}})~{},$
(20)

##### Moving frame

$\displaystyle\mathrm{O}_{5q,1}^{G_{1},\uparrow/\downarrow}(\mathbf{p}=\mathbf{e_{i}})$
$\displaystyle=\pm\mathrm{O}^{\uparrow/\downarrow}_{3q}(\mathbf{e_{i}})\mathrm{O}_{\bar{q}q}(\mathbf{0})~{},$
(21)
$\displaystyle\mathrm{O}_{5q,2}^{G_{1},\uparrow/\downarrow}(\mathbf{p}=\mathbf{e_{i}})$
$\displaystyle=\pm\mathrm{O}^{\uparrow/\downarrow}_{3q}(\mathbf{0})\mathrm{O}_{\bar{q}q}(\mathbf{e_{i}})~{}.$
(22)

The superscripts $\uparrow/\downarrow$ refer to the polarization of the spin.
The momentum projection of $3q$\- and $\bar{q}q$-operators is performed via
Fourier transform

$\mathrm{O}_{3q}(\mathbf{p}_{1})\mathrm{O}_{\bar{q}q}(\mathbf{p}_{2})=\sum_{\mathbf{x}}e^{-i\mathbf{p}_{1}\cdot\mathbf{x}}\mathrm{O}_{3q}(\mathbf{x})\sum_{\mathbf{y}}e^{-i\mathbf{p}_{2}\cdot\mathbf{y}}\mathrm{O}_{\bar{q}q}(\mathbf{y})~{},$
(23)

which requires the two operators to be located at different spatial position
on the same timeslice. While nucleons have isospin $I_{N}=\frac{1}{2}$, pions
have $I_{\pi}=1$. In general, the five-quark interpolator has overlap with
isospin states $I=\frac{1}{2}\oplus 1=\frac{1}{2},\frac{3}{2}$. We use the
Clebsch-Gordan coefficients to project the two-hadron operator to isospin
$I=\frac{1}{2}$ and $I_{z}=\pm\frac{1}{2}$, to obtain

$\displaystyle\left(I=\frac{1}{2},I_{z}=+\frac{1}{2}\right)\hskip 56.9055pt$
$\displaystyle\mathrm{O}_{5q}^{p}=+\sqrt{\frac{2}{3}}\mathrm{O}_{3q}^{p}\mathrm{O}_{\bar{q}q}^{\pi^{0}}-\frac{1}{\sqrt{3}}\mathrm{O}_{3q}^{n}\mathrm{O}_{\bar{q}q}^{\pi^{+}}~{},\hskip
56.9055pt$ (24)
$\displaystyle\left(I=\frac{1}{2},I_{z}=-\frac{1}{2}\right)\hskip 56.9055pt$
$\displaystyle\mathrm{O}_{5q}^{n}=-\frac{1}{\sqrt{3}}\mathrm{O}_{3q}^{n}\mathrm{O}_{\bar{q}q}^{\pi^{0}}+\sqrt{\frac{2}{3}}\mathrm{O}_{3q}^{p}\mathrm{O}_{\bar{q}q}^{\pi^{-}}~{},\hskip
56.9055pt$ (25)

where the superscripts $p,n,\pi$ are short-hand notations for proton, neutron,
and pion isospin quantum numbers. The expressions for these standard operators
read

$\displaystyle\mathrm{O}_{3q}^{p}(x)^{\gamma}=\epsilon_{abc}\left[u^{\top}(x)^{a}_{\alpha}(C\gamma_{5})^{\alpha\beta}d(x)^{b}_{\beta}\right]q(x)^{c}_{\gamma}~{},$
(26)
$\displaystyle\mathrm{O}_{3q}^{n}(x)^{\gamma}=\epsilon_{abc}\left[u^{\top}(x)^{a}_{\alpha}(C\gamma_{5})^{\alpha\beta}d(x)^{b}_{\beta}\right]d(x)^{c}_{\gamma}~{},$
(27)
$\displaystyle\mathrm{O}_{\bar{q}q}^{\pi^{-}}(y)=\bar{u}(y)\gamma_{5}d(y)~{},$
(28)
$\displaystyle\mathrm{O}_{\bar{q}q}^{\pi^{-}}(0)=\frac{1}{\sqrt{2}}\left[\bar{u}(y)\gamma_{5}u(y)-\bar{d}(y)\gamma_{5}d(y)\right]~{},$
(29)
$\displaystyle\mathrm{O}_{\bar{q}q}^{\pi^{+}}(y)=-\bar{d}(y)\gamma_{5}u(y)~{}.$
(30)

The minus sign in eq. (30) is derived by applying the ladder isospin operator
$I^{+}$ to the operator in eq. (28). We combine the isospin projections with
the lattice irrep projections in eqs. (19)-(22) to construct the final two-
hadron operator with the quantum numbers of protons and neutrons. To construct
extended $N$\- and $N\pi$-like operators we apply Gaussian smearing to the
quark fields [40] and APE smoothing to the gauge fields [41]. We use the
smearing parameters in [36].

###  4.2 Computation of correlation functions

Using the single- and two-hadron operators in eqs. (24)-(27) we compute the
matrix of three-point correlation functions

$\mathbb{C}_{3pt}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)=\begin{pmatrix}\langle\mathrm{O}_{3q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle&\langle\mathrm{O}_{3q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{5q}(\mathbf{p},0)\rangle\\\
\langle\mathrm{O}_{5q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle&\langle\mathrm{O}_{5q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{5q}(\mathbf{p},0)\rangle\end{pmatrix}~{}.$
(31)

$u(x)$$u(x)$$d(x)$$p$$\bar{d}(z)$$u(z)$$\bar{u}(y)$$d(y)$$\pi^{-}$$\bar{u}(0)$$\bar{u}(0)$$\bar{d}(0)$$p$A-like

$d(x)$$u(x)$$u(x)$$p$$\bar{d}(z)$$u(z)$$\bar{u}(y)$$d(y)$$\pi^{-}$$\bar{u}(0)$$\bar{u}(0)$$\bar{d}(0)$$p$B-like

$d(x)$$u(x)$$u(x)$$p$$\bar{d}(z)$$u(z)$$\bar{u}(y)$$d(y)$$\pi^{-}$$\bar{d}(0)$$\bar{u}(0)$$\bar{u}(0)$$p$C-like

$u(x)$$u(x)$$d(x)$$p$$\bar{d}(z)$$u(z)$$\bar{u}(y)$$d(y)$$\pi^{-}$$\bar{u}(0)$$\bar{u}(0)$$\bar{d}(0)$$p$D-like

Figure 1:
This schematic plot represents the topologies $A$, $B$, $C$ and $D$ of the
Wick contractions in the correlation functions
$\langle\mathrm{O}_{5q}^{n}(\mathbf{p}^{\prime},t)~{}\mathcal{J}^{-}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}^{p}_{3q}(\mathbf{p},0)\rangle$.
The red lines are all-to-all propagators, two for each diagram, and the black
lines are point-to-all propagators. The red cross is the current insertion.

We neglect the lower diagonal entry
$\langle\mathrm{O}_{5q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{5q}(\mathbf{p},0)\rangle$
because it corresponds to $N\pi\to N\pi$ terms, which are expected to be
volume suppressed in the standard three-point functions
$\langle\mathrm{O}_{3q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle$,
and therefore their contribution will be much smaller. Also, it is numerically
more demanding to compute these correlation functions with our methods as it
involves several Wick contractions and all-to-all propagators. The off-
diagonal entries in eq. (31) correspond to the current-enhanced terms (when
allowed) in the standard three-point functions, i.e., $N\to N\pi$ and $N\pi\to
N$ states. The Wick contractions of
$\langle\mathrm{O}_{5q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle$
were first discussed in [24], and they can be grouped in 4 diagrams: A, B, C,
and D, see Fig. 1. The diagrams contain 2 all-to-all propagators and we adopt
the sequential method [31] with 2 sequential propagators to compute the
diagrams A, B, and C. The diagram D is much cheaper to compute because we use
the point-to-all method for the nucleon part and the one-end-trick [32] with
12 stochastic vectors for the pion-to-current term. In particular, our lattice
simulations show that for the diagram D

$\langle\mathrm{O}_{5q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle_{D}\approx\langle\mathrm{O}_{3q}(\mathbf{p}^{\prime}_{N},t)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle~{}\langle\mathrm{O}_{\bar{q}q}(\mathbf{p}^{\prime}_{\pi},t)~{}\mathcal{J}(\mathbf{q},\tau)\rangle~{}.$
(32)

This factorisation, together with the observation that the diagram D has the
largest signal, explains which channels have the large $N\pi$ contamination in
standard nucleon three-point functions. In the case of the axial-vector
current $\mathcal{J}=\mathcal{A}_{\mu}$, the diagram D reads

$\displaystyle\langle\mathrm{O}_{5q}(\mathbf{p}^{\prime},t)~{}\mathcal{A}_{\mu}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle_{D}$
$\displaystyle\propto\delta_{\mathbf{p}^{\prime},\mathbf{p}}~{}\delta_{\mathbf{p}^{\prime}_{\pi},\mathbf{q}}~{}|\langle\Omega|\mathrm{O}_{3q}|N\rangle|^{2}if_{\pi}q_{\mu}~{}e^{-E_{N}t}e^{-E_{\pi}(t-\tau)}$
(33)

which is non-zero only when $\mathbf{p}^{\prime}_{N}=\mathbf{p}$,
$\mathbf{p}^{\prime}_{\pi}=\mathbf{q}$ and $q_{\mu}\neq 0$. This explains the
large ($N\pi$) excited state contamination in the $\mathcal{A}_{4}$ channel,
and the milder excited state contamination in the channel with
$\mathcal{A}_{i}$ and $q_{i}=0$, where the momentum of the current is zero
along the direction where the spin is aligned and the current component, see
Fig. 4 of [21]. Similar arguments hold with a pseudoscalar current.

Figure 2:
Two-point correlation functions on the left and respective effective masses on
the right. The diagram B and C, which consists of $3$-quark propagators from
time $0$ to $t$ are decaying immediately with the nucleon energy, while
diagrams A and D, which consists of $4$-quark and $1$-antiquark propagators
from $0$ to $t$ are decaying with an interacting energy close to the lowest
non-interacting $N\pi$ (and $N\pi\pi$).

In order to use the variational method, we compute the $2\times 2$ GEVP matrix

$\mathbb{C}_{2pt}(\mathbf{p},t)=\begin{pmatrix}\langle\mathrm{O}_{3q}(\mathbf{p},t)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle&\langle\mathrm{O}_{3q}(\mathbf{p},t)~{}\bar{\mathrm{O}}_{5q}(\mathbf{p},0)\rangle\\\
\langle\mathrm{O}_{5q}(\mathbf{p},t)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle&\langle\mathrm{O}_{5q}(\mathbf{p},t)~{}\bar{\mathrm{O}}_{5q}(\mathbf{p},0)\rangle\\\
\end{pmatrix}~{},$ (34)

see discussion in sec. (3). The correlation functions
$\langle\mathrm{O}_{3q}(\mathbf{p},t)~{}\bar{\mathrm{O}}_{5q}(\mathbf{p},0)\rangle$
are constructed from
$\langle\mathrm{O}_{3q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle$
with $\mathcal{J}=\mathcal{P}$, $\tau=0$, and at different source-sink
separations $t$. Similarly, the ($N\pi$-like) two-point correlation functions
$\langle\mathrm{O}_{5q}(\mathbf{p},t)~{}\bar{\mathrm{O}}_{5q}(\mathbf{p},0)\rangle$
are constructed from the three-point functions
$\langle\mathrm{O}_{5q}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{3q}(\mathbf{p},0)\rangle$;
to visualise this construction see Fig. 1 with the current
$\mathcal{J}=\mathcal{P}$ at the source ($\tau=0$). The sequential method,
which is used in this pilot study, is not the most efficient approach and in
the future other methods like distillation [33] will be considered, especially
when $N\pi\pi$ operators will be included in the basis. In Fig. 2, there is a
comparison plot between the diagrams A, B, C, and D that make up the
$N\pi$-like two-point functions and the standard nucleon-like two-point
functions. On the right plot, the effective energies of the diagrams and their
sum are compared against $N$, $N\pi$ and $N\pi\pi$ energies in the non-
interacting limit. The operators are projected to total momentum zero in these
plots. While the diagrams A and D decay with some interacting $N\pi$ energy
which is very close to the non-interacting $N\pi$ energy (or $N\pi\pi$ energy
as they are degenerate in this finite volume), the diagrams B and C decay
quite fast with the nucleon energy. This behaviour is expected, see e.g. Fig.
7 of [42], as the diagrams A and D are constructed from 5 propagators
travelling from time $0$ to time $t$, while the diagrams B and C contain just
3 propagators from time $0$ to time $t$ and 2 propagators from time $t$ to
time $t$, see Fig. 1. Eventually, the sum of all the diagrams that compose the
$N\pi$-like two-point function will decay at large time with the nucleon
energy as the $N\pi$ operators are constructed to have the quantum numbers of
the nucleon.

###  4.3 GEVP results

We solve the GEVP in eq. (9) using the matrix in eq. (34) to determine the
matrix of eigenvectors $V(\mathbf{p},t;t_{0})$ and eigenvalues
$\Lambda(\mathbf{p},t;t_{0})$. In Fig. 3, the GEVP results at
$t_{0}=0.2~{}\mathrm{fm}$ are plotted for $\mathbf{p}=\mathbf{0}$ and
$\mathbf{p}=\mathbf{e_{i}}$ ($|\mathbf{p}|\approx 520~{}\mathrm{MeV}$). For
the kinematic frame with momentum $\mathbf{p}$, we use both the $5q$-operators
in eqs. (21)-(22). On the left plots, the effective energies of the 2 GEVP
eigenvalues are compared to the $N$, $N\pi$ and $N\pi\pi$ energies in the non-
interacting limit. The nucleon and pion energies are extracted at large time
from nucleon and pion two-point functions. In the rest frame with zero total
momentum, the energy of the first eigenvalue matches the nucleon energy. The
second eigenvalue has an energy very close to that of the non-interacting
$N\pi$ P-wave state and the $N\pi\pi$ S-wave state, which are
(quasi-)degenerate in this ensemble. This comparison between $N\pi$ and
$N\pi\pi$ states in the non-interacting limit varies with the pion mass and
the physical volume, see e.g. Fig. 2 in [43]. In the moving frames with total
momentum $\mathbf{p}$, the first eigenvalue is again consistent with the
nucleon state and the second eigenvalue is more clearly a $N\pi$ state, due to
the non-degeneracy of the $N\pi$ and $N\pi\pi$ states in the non-interacting
limit. We thus refer to the two GEVP eigenvalues as $\lambda^{N}$ and
$\lambda^{N\pi}$. The GEVP results for the moving frames in Fig. 3 are
averaged over equivalent momenta
($\mathbf{p}=\hat{e}_{x},\hat{e}_{y},\hat{e}_{z},-\hat{e}_{x},-\hat{e}_{y},-\hat{e}_{z}$).
We have to stress that the ensemble is at the $SU(3)$ symmetric point and this
higher symmetries imply more degeneracy with other physical states, e.g.
$\Lambda K$ are degenerate to $N\pi$. This degeneracy could in principle
affect the GEVP eigensolutions and must be taken into account to reduce
systematic effects.

On the right plots, the components of the normalised GEVP eigenvectors are
presented. The eigenvectors are normalised according to eq. (10) and we refer
to them as $\tilde{v}^{N}$ and $\tilde{v}^{N\pi}$. In particular, the
components are $\tilde{v}^{N}=(\tilde{v}^{N}_{3q},\tilde{v}^{N}_{5q})$ and
$\tilde{v}^{N\pi}=(\tilde{v}^{N\pi}_{3q},\tilde{v}^{N\pi}_{5q})$.

Figure 3: Effective energies of GEVP eigenvalues (green and red data points)
on the left, compared to $N\pi$ and $N\pi\pi$ non-interacting energies, as
well as nucleon energies extracted from nucleon two-point functions. The pion
energies are extracted from pion two-point functions. (top) Energies in the
rest frame with total momentum zero. (middle) GEVP energies in the moving
frame using the $5q$ operator in eq. (22), where the nucleon operator is
projected to zero momentum. (bottom) GEVP energies in the moving frame using
the $5q$ operator in eq. (21), where the pion operator is projected to zero
momentum. The relative momentum is $\mathbf{p}=\frac{2\pi}{L}\hat{n}$, with
$\hat{n}$ a unit vector along the 3 spatial directions. On the right, GEVP
eigenvectors normalised w.r.t. $\mathbb{C}_{2pt}(t_{0})$ as in eq. (10).

###  4.4 GEVP ratios $\langle N|\mathcal{J}|N\rangle$

Using the eigenvectors $v^{N}(\mathbf{p})$ in the plateau region, i.e.
$0.5~{}\mathrm{fm}\lesssim t\lesssim 0.9~{}\mathrm{fm}$ in Fig. 3, we
construct the GEVP-projected interpolating operator

$\mathrm{O}^{N}(\mathbf{p})=\sum_{i}v_{i}^{N}(\mathbf{p})\mathrm{O}_{i}(\mathbf{p})\hskip
60.00009pt\bar{\mathrm{O}}^{N}(\mathbf{p})=\sum_{i}v_{i}^{N}(\mathbf{p})\bar{\mathrm{O}}_{i}(\mathbf{p})$
(35)

where we emphasize that the eigenvectors are fairly constant with $t$, and use
eqs. (13)-(14) to construct the GEVP-improved two- and three-point correlation
functions

$\displaystyle C_{2pt}(\mathbf{p},t)^{N}$
$\displaystyle=\langle\mathrm{O}^{N}(\mathbf{p},t)~{}\bar{\mathrm{O}}^{N}(\mathbf{p},0)\rangle~{},$
(36) $\displaystyle
C^{k,\mathcal{J}}_{3pt}(\mathbf{p}^{\prime},t;\mathbf{q},\tau)^{NN}$
$\displaystyle=\langle\mathrm{O}^{N}_{k}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{k}^{N}(\mathbf{p})\rangle~{}.$
(37)

The subscript $k$ refers to the polarization of the interpolating operators.
After the GEVP-projection, one can check that the two-point functions decay
with an energy close to the nucleon energy and that the three-point functions
decay with $\tau$ with the source-sink energy difference
$E_{N}-E^{\prime}_{N}$. It is important to emphasize that the inclusion of
$N\pi$ operators in the variational basis did not improve much the plateau of
the nucleon two-point functions, but it makes a substantial change in the
nucleon three-point functions.

#### Forward limit

The GEVP-improved correlation functions are used to construct the GEVP ratios
in eq. (15). In the forward limit, we set $\mathbf{p}^{\prime}=\mathbf{p}$ and
$\mathbf{q}=\mathbf{0}$, and investigate both $\mathbf{p}=\mathbf{0}$ and
$\mathbf{p}=\hat{e}_{i}$. We expect from LO-ChPT that only the channels
$\mathcal{J}=\mathcal{A}_{4},\mathcal{P}$ with
$\mathbf{p}=\mathbf{p}^{\prime}=\hat{e}_{i}$ have current-enhanced $N\pi$
contamination, while in the other channels, namely $\mathcal{A}_{4}$ and
$\mathcal{P}$ with $\mathbf{p}=\mathbf{0}$ and $\mathcal{A}_{i}$ with
$\mathbf{p}=\mathbf{0},\hat{e}_{i}$ the excited state contamination is overall
much smaller. This expectation is also confirmed by our numerical study, where
in [23] we report the channels affected by large ESC. The inclusion in our
basis of $N\pi$-like interpolating operators with unit relative momenta has
little effect in the determination of $g_{A}$ with
$\mathcal{J}=\mathcal{A}_{i}$ and with $\mathbf{p}=\mathbf{0},\hat{e}_{i}$. As
discussed in the previous section, this is mainly due to the diagram D which
is not contributing to these channel when $q_{i}=0$, see eq. (33), but the
other diagrams have a non-zero but rather small amplitude.

#### Off-forward limit

In the off-forward limit, we expect from LO-ChPT a large $N\pi$ contamination
in the channels with $\mathcal{J}=\mathcal{A}_{\mu},\mathcal{P}$ and with
$q_{i}\neq 0$ where $i$ is the direction of the polarization. Our variational
analysis with $N$\- and $N\pi$-like operators confirm again these predictions
and it shows that the $N\pi$ operators included in this work have little
effect in the off-forward limit channel $\mathcal{A}_{i}$ with $q_{i}=0$, e.g.
$\mathcal{A}_{i\neq 3}$, $\mathbb{P}^{+}_{i\neq 3}$ and
$\mathbf{q}=\frac{2\pi}{L}(0,0,1)$ like in Fig. 4 of [21], where the excited
state contamination is much smaller. This can be seen again from the behaviour
of the diagram D in eq. (33). A comparison of standard and GEVP-improved
ratios in the off-forward limit is presented in [23]. Full results in the
forward and off-forward limit will appear in a separate publication.

###  4.5 GEVP ratios for the extraction of $\langle
N\pi(\mathbf{p}^{\prime})|\mathcal{J}(\mathbf{q})|N(\mathbf{p})\rangle$

We compute for the first time the transition matrix elements $\langle
N\pi(\mathbf{p})|\mathcal{J}(\mathbf{0})|N(\mathbf{p})\rangle$ by constructing
the GEVP-improved operators

$\displaystyle\mathrm{O}^{N}(\mathbf{p})=\sum_{i=3q,5q}v_{i}^{N}(\mathbf{p})\mathrm{O}_{i}(\mathbf{p})\hskip
60.00009pt\bar{\mathrm{O}}^{N\pi}(\mathbf{p})=\sum_{i=3q,5q}v_{i}^{N\pi}(\mathbf{p})\bar{\mathrm{O}}_{i}(\mathbf{p})$
(38)

and by computing the GEVP-improved correlator in eq. (36) together with

$\displaystyle C_{2pt}(\mathbf{p},t)^{N\pi}$
$\displaystyle=\langle\mathrm{O}^{N\pi}(\mathbf{p},t)~{}\bar{\mathrm{O}}^{N\pi}(\mathbf{p},0)\rangle~{},$
(39) $\displaystyle
C^{k,\mathcal{J}}_{3pt}(\mathbf{p},t;\mathbf{0},\tau)^{N\pi N}$
$\displaystyle=\langle\mathrm{O}_{k}^{N\pi}(\mathbf{p}^{\prime},t)~{}\mathcal{J}(\mathbf{q},\tau)~{}\bar{\mathrm{O}}_{k}^{N}(\mathbf{p})\rangle~{}.$
(40)

Then we construct the GEVP-ratio in eq. (16) and in Fig. (4) we show the
renormalised GEVP-ratios with $\mathcal{J}=\mathcal{A}_{4},\mathcal{P}$ for
different source-sink separations and with
$\mathbf{p}^{\prime}=\mathbf{p}=\hat{e}_{k}$. The resulting GEVP ratios are
much flatter than the ones for the extraction of $\langle
N|\mathcal{J}|N\rangle$, which seems very promising for future calculations.
These lattice results can be compared against ChPT, which predicts that

$\frac{\langle
N(\mathbf{p})\pi(\mathbf{0})|\mathcal{P}(\mathbf{0})|N(\mathbf{p})\rangle}{\langle
N(\mathbf{p})\pi(\mathbf{0})|\mathcal{A}_{4}(\mathbf{0})|N(\mathbf{p})\rangle}=\frac{m_{\pi}}{2m_{\ell}}~{}.$
(41)

This relation should hold up to $\mathcal{O}(a)$-effects as the intermediate
currents are not $\mathcal{O}(a)$-improved. The ratio of matrix elements is
equivalent to the double ratio
$\frac{R^{k,\mathcal{P}}(\mathbf{p},t;\mathbf{0},\tau)^{N\pi
N}}{R^{k,\mathcal{A}_{4}}(\mathbf{p},t;\mathbf{0},\tau)^{N\pi
N}}=\frac{C^{k,\mathcal{P}}_{3pt}(\mathbf{p},t;\mathbf{0},\tau)^{N\pi
N}}{C^{k,\mathcal{A}_{4}}_{3pt}(\mathbf{p},t;\mathbf{0},\tau)^{N\pi N}}$, and
the two ratios are plotted on the right panel in Fig. 4. It is reassuring to
see that they are not completely in disagreement with ChPT as
$\mathcal{O}(a)$-effects and $SU(3)$ symmetry effects must be still taken into
account. These transition matrix elements are at a momentum transfer
$Q^{2}=m_{\pi}^{2}\approx 0.175~{}\mathrm{GeV}^{2}$. Results for the matrix
elements with $\mathbf{q}\neq\mathbf{0}$, which are of paramount importance
for the analysis of the neutrino oscillation data and more, will be presented
in a separate publication.

Figure 4:
Renormalised GEVP-ratios in eq. (16) at $\mathbf{q}=\mathbf{0}$,
$\mathbf{p}=\hat{e}_{i}$ using the pseudoscalar (left) and temporal axial-
vector (middle) current. The GEVP ratios for $\langle
N\pi|\mathcal{J}|N\rangle$ in these plots result much flatter than the GEVP
ratios for the extraction of $\langle N|\mathcal{J}|N\rangle$, see [23]. A
constant fit using only the largest source-sink separations give for these
ratios: $R^{i,\mathcal{P}}=14.22(39)$ and $R^{i,\mathcal{A}_{4}}=2.43(8)$. On
the right plot, there is a comparison between the double ratio, which is
equivalent to the ratio of transition matrix elements, and its ChPT
prediction.

##  5 Conclusions

Elastic nucleon form factors $N\to N$, as well as transition form factors for
$N\to\Delta$, $N\to N^{*}$ and $N\to N\pi$ are of paramount importance for
analysis of neutrino-nucleus scattering data in the resonance region
production. In this study, we used a variational analysis to explicitly remove
the $N\to N\pi$ matrix elements from the standard nucleon-like three-point
functions. As a by-product, we present some preliminary results for the
transition matrix elements $N\to N\pi$ at vanishing spatial momentum transfer
for the isospin transition $1/2\to 1/2$. Results for the elastic nucleon form
factors were presented in [23]. We find that with the GEVP analysis the ratios
tend to be surprisingly flat at source-sink separations $\gtrsim
0.8~{}\mathrm{fm}$, which allowed a simple constant fit to the data at the
largest available source-sink separation. The results suffer from
$\mathcal{O}(a)$ effects as the currents are not $\mathcal{O}(a)$-improved.
The ratio of $\mathcal{P}$ and $\mathcal{A_{4}}$ matrix elements is found to
be within $\sim 10\%$ of the ChPT prediction. Other results at non-vanishing
momentum transfer and with other currents will be presented in a separate
publication. Our lattice simulations are performed on a single $N_{f}=3$
ensemble with lattice spacing $a=0.1~{}\mathrm{fm}$,
$m_{\pi}=420~{}\mathrm{MeV}$ and $T=2L=4.8~{}\mathrm{fm}$. The $SU(3)$
symmetry introduces another small complication as the $N\pi$ states are
degenerate to $\Lambda K$ states. Although we expect this to not have any
large effect in the three-point functions, they might have some few percent
effect in the nucleon two-point functions and therefore they can introduce
some small systematic effects in the GEVP eigensolutions. This work represents
only the first step towards the calculation of transition form factors with a
more complete variational analysis, where single-hadron and multi-hadron
operators must be included in a Lüscher formalism to reconstruct the spectra.
The results confirm the ChPT prediction of very large $N\pi$ states
contamination in some nucleon-like three-point functions and we emphasize that
one particular lattice diagram (D) is the most important one.

#### Acknowledgements

L. B. would like to thank the staff members of the lattice group at the
Lawrence Berkeley National Laboratory, namely K.-F. Liu, W. Bigeng, R.
Briceño, A. Walker-Loud, and D. Pefkou, for their warm hospitality and for the
interesting discussions. L. B. has benefited from discussions with O. Bär, R.
Briceño, J. Green, R. Gupta, L. Leskovec, and K.-F. Liu. The work is supported
by the German Research Foundation (DFG) research unit FOR5269 "Future methods
for studying confined gluons in QCD". Simulations were performed on the QPACE
3 computer of SFB/TR-55, using an adapted version of the Chroma [44] software
package.

## References

  * [1] R. Acciarri et al.  Long-Baseline Neutrino Facility (LBNF) and Deep Underground Neutrino Experiment (DUNE): Conceptual Design Report, Volume 2: The Physics Program for DUNE at LBNF.  12 2015.
  * [2] Babak Abi et al.  Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume I Introduction to DUNE.  JINST, 15(08):T08008, 2020.
  * [3] J. Bian et al.  Hyper-Kamiokande experiment: a Snowmass white paper.  In 2022 Snowmass Summer Study, 3 2022.
  * [4] B. Pontecorvo.  Neutrino Experiments and the Problem of Conservation of Leptonic Charge.  Zh. Eksp. Teor. Fiz., 53:1717–1725, 1967.
  * [5] A. D. Sakharov.  Violation of CP Invariance, C asymmetry, and baryon asymmetry of the universe.  Pisma Zh. Eksp. Teor. Fiz., 5:32–35, 1967.
  * [6] B. Abi et al.  Supernova neutrino burst detection with the Deep Underground Neutrino Experiment.  Eur. Phys. J. C, 81(5):423, 2021.
  * [7] B. Abi et al.  Prospects for beyond the Standard Model physics searches at the Deep Underground Neutrino Experiment.  Eur. Phys. J. C, 81(4):322, 2021.
  * [8] Joachim Kopp, Noemi Rocco, and Zahra Tabrizi.  Unleashing the Power of EFT in Neutrino-Nucleus Scattering.  1 2024.
  * [9] Alessandro Lovato, Noemi Rocco, and Noah Steinberg.  One and Two-Body Current Contributions to Lepton-Nucleus Scattering.  12 2023.
  * [10] A. A. Aguilar-Arevalo, C. E. Anderson, Bazarko, et al.  First measurement of the muon neutrino charged current quasielastic double differential cross section.  Phys. Rev. D, 81:092005, May 2010.
  * [11] J. A. Formaggio and G. P. Zeller.  From ev to eev: Neutrino cross sections across energy scales.  Rev. Mod. Phys., 84:1307–1341, Sep 2012.
  * [12] Daniel Simons, Noah Steinberg, Alessandro Lovato, Yannick Meurice, Noemi Rocco, and Michael Wagman.  Form factor and model dependence in neutrino-nucleus cross section predictions.  10 2022.
  * [13] L. Alvarez Ruso et al.  Theoretical tools for neutrino scattering: interplay between lattice QCD, EFTs, nuclear physics, phenomenology, and neutrino event generators.  3 2022.
  * [14] Andreas S. Kronfeld et al.  Lattice QCD and Particle Physics.  7 2022.
  * [15] Martin Luscher.  Two particle states on a torus and their relation to the scattering matrix.  Nucl. Phys. B, 354:531–578, 1991.
  * [16] Martin Luscher.  Signatures of unstable particles in finite volume.  Nucl. Phys. B, 364:237–251, 1991.
  * [17] Benoit Blossier, Michele Della Morte, Georg von Hippel, Tereza Mendes, and Rainer Sommer.  On the generalized eigenvalue method for energies and matrix elements in lattice field theory.  JHEP, 04:094, 2009.
  * [18] John Bulava, Michael Donnellan, and Rainer Sommer.  On the computation of hadron-to-hadron transition matrix elements in lattice QCD.  JHEP, 01:140, 2012.
  * [19] Oliver Bar.  $N\pi$-state contamination in lattice calculations of the nucleon axial form factors.  Phys. Rev. D, 99(5):054506, 2019.
  * [20] Oliver Bar.  $N\pi$-state contamination in lattice calculations of the nucleon pseudoscalar form factor.  Phys. Rev. D, 100(5):054507, 2019.
  * [21] Gunnar S. Bali, Lorenzo Barca, Sara Collins, Michael Gruber, Marius Löffler, Andreas Schäfer, Wolfgang Söldner, Philipp Wein, Simon Weishäupl, and Thomas Wurm.  Nucleon axial structure from lattice QCD.  JHEP, 05:126, 2020.
  * [22] Yong-Chull Jang, Rajan Gupta, Boram Yoon, and Tanmoy Bhattacharya.  Axial Vector Form Factors from Lattice QCD that Satisfy the PCAC Relation.  Phys. Rev. Lett., 124(7):072002, 2020.
  * [23] Lorenzo Barca, Gunnar Bali, and Sara Collins.  Toward N to N$\pi$ matrix elements from lattice QCD.  Phys. Rev. D, 107(5):L051505, 2023.
  * [24] Lorenzo Barca, Gunnar S. Bali, and Sara Collins.  Investigating $N\to N\pi$ axial matrix elements.  PoS, LATTICE2021:359, 2022.
  * [25] Oliver Bar.  Nucleon-pion-state contribution to nucleon two-point correlation functions.  Phys. Rev. D, 92(7):074504, 2015.
  * [26] Oliver Bär.  Three-particle $N\pi\pi$ state contribution to the nucleon two-point function in lattice QCD.  Phys. Rev. D, 97(9):094507, 2018.
  * [27] Daniel Severt, Maxim Mai, and Ulf-G. Meißner.  Particle-dimer approach for the Roper resonance in a finite volume.  JHEP, 04:100, 2023.
  * [28] Xuan-He Wang, Cong-Ling Fan, Xu Feng, Lu-Chang Jin, and Zhao-Long Zhang.  Nucleon electric polarizabilities and nucleon-pion scattering at physical pion mass.  10 2023.
  * [29] M. Luscher.  Volume Dependence of the Energy Spectrum in Massive Quantum Field Theories. 2. Scattering States.  Commun. Math. Phys., 105:153–188, 1986.
  * [30] David J. Wilson, Raul A. Briceno, Jozef J. Dudek, Robert G. Edwards, and Christopher E. Thomas.  Coupled $\pi\pi,K\bar{K}$ scattering in $P$-wave and the $\rho$ resonance from lattice QCD.  Phys. Rev. D, 92(9):094502, 2015.
  * [31] G. Martinelli and C. T. Sachrajda.  A Lattice Study of Nucleon Structure.  Nucl. Phys. B, 316:355–372, 1989.
  * [32] Philippe Boucaud et al.  Dynamical Twisted Mass Fermions with Light Quarks: Simulation and Analysis Details.  Comput. Phys. Commun., 179:695–715, 2008.
  * [33] Michael Peardon, John Bulava, Justin Foley, Colin Morningstar, Jozef Dudek, Robert G. Edwards, Balint Joo, Huey-Wen Lin, David G. Richards, and Keisuke Jimmy Juge.  A Novel quark-field creation operator construction for hadronic physics in lattice QCD.  Phys. Rev. D, 80:054506, 2009.
  * [34] Colin Morningstar, John Bulava, Justin Foley, Keisuke J. Juge, David Lenkner, Mike Peardon, and Chik Him Wong.  Improved stochastic estimation of quark propagation with Laplacian Heaviside smearing in lattice QCD.  Phys. Rev. D, 83:114505, 2011.
  * [35] Anthony V. Grebe and Michael Wagman.  Nucleon-Pion Spectroscopy from Sparsened Correlators.  11 2023.
  * [36] Gunnar S. Bali, Sara Collins, Peter Georg, Daniel Jenkins, Piotr Korcyl, Andreas Schäfer, Enno E. Scholz, Jakob Simeth, Wolfgang Söldner, and Simon Weishäupl.  Scale setting and the light baryon spectrum in Nf = 2 + 1 QCD with Wilson fermions.  JHEP, 05:035, 2023.
  * [37] S. Prelovsek, U. Skerbis, and C. B. Lang.  Lattice operators for scattering of particles with spin.  J. High Energy Phys., 01:129, 2017.
  * [38] C. B. Lang, L. Leskovec, M. Padmanath, and S. Prelovsek.  Pion-nucleon scattering in the Roper channel from lattice QCD.  Phys. Rev. D, 95(1):014510, 2017.
  * [39] M. Göckeler, R. Horsley, M. Lage, U.-G. Meißner, P. E. L. Rakow, A. Rusetsky, G. Schierholz, and J. M. Zanotti.  Scattering phases for meson and baryon resonances on general moving-frame lattices.  Phys. Rev. D, 86:094513, 2012.
  * [40] S. Güsken.  A study of smearing techniques for hadron correlation functions.  Nucl. Phys. B, 17:361–364, 1990.
  * [41] (APE Collaboration) M. Albanese et al.  Glueball Masses and String Tension in Lattice QCD.  Phys. Lett. B, 192:163–169, 1987.
  * [42] Keh-Fei Liu.  PDF in PDFs from Hadronic Tensor and LaMET.  Phys. Rev. D, 102(7):074502, 2020.
  * [43] Jeremy Green.  Systematics in nucleon matrix element calculations.  PoS, LATTICE2018:016, 2018.
  * [44] Robert G. Edwards and Balint Joó.  The Chroma software system for Lattice QCD.  Nucl. Phys. B Proc. Suppl., 140:832, 2005.


11institutetext:  ConsenSys, New York, USA
11email<EMAIL_ADDRESS>
# Formal and Executable Semantics of the Ethereum Virtual Machine in Dafny

Franck Cassez  Joanne Fuller  Milad K. Ghale  David J. Pearce Horacio M. A.
Quiles

###### Abstract

The Ethereum protocol implements a replicated state machine. The network
participants keep track of the system state by: 1) agreeing on the sequence of
transactions to be processed and 2) computing the state transitions that
correspond to the sequence of transactions. Ethereum transactions are
programs, called _smart contracts_ , and computing a state transition requires
executing some code. The Ethereum Virtual Machine (EVM) provides this
capability and can execute programs written in EVM _bytecode_. We present a
formal and executable semantics of the EVM written in the verification-
friendly language Dafny: it provides $(i)$ a readable, formal and verified
specification of the semantics of the EVM; $(ii)$ a framework to formally
reason about bytecode.

##  1 Introduction

A distinctive feature of Ethereum is that transactions are programs, _smart
contracts_ , and computing a state transition requires to run the contract
code to compute the next state. This capability is provided by the Ethereum
Virtual Machine (EVM) that can execute programs written in EVM _bytecode_. The
original and informal specification of the EVM is in the Yellow Paper [28].

As a decentralised platform, Ethereum encourages _client diversity_ : network
participants are free to choose which implementation of the EVM they want to
run, and there are several implementations to choose from written in different
languages e.g., Go, Java. All the EVM implementations must agree on the state
transitions, otherwise the network would split and the blockchain would
_fork_. However, the original specification in the Yellow Paper [28] has some
known shortcomings: ($i$) it is hard to read and does not provide a formal
semantics of the EVM and the bytecode; ($ii$) the lack of a formal semantics
makes it hard for Ethereum client developers to guarantee that they interpret
the Yellow Paper in a consistent way; ($iii$) designing compilers from high-
level languages (e.g., Solidity111The most popular language to write smart
contracts.) to EVM bytecode without a formal semantics is error-prone and,
without a precise semantics of the EVM, it is hard to design _certified_
compilers (preserving of semantics from a high to a low-level language.).

One can argue that existing implementations of the EVM (e.g., in Go, Java)
provide a de facto semantics for it. Whilst this is true to some extent, such
implementations do not enable _formal reasoning_ about bytecode. Furthermore,
whilst smart contracts can be written in high-level languages like Solidity,
they must be compiled into EVM bytecode before being executed on the EVM.
Tools for checking safety properties (e.g., absence of overflow, division by
zero, etc) at the Solidity level are problematic if they cannot guarantee such
properties hold at the bytecode level. One solution is to design a provably
correct compiler, but this is a complex and long-term endeavour [17].
Alternatively we can provide techniques, supported by tools, to reason about
properties of the bytecode. This is what we propose to do in this work.

### Our Contribution.

We present a complete and formal specification of the EVM in Dafny, available
at https://github.com/ConsenSys/evm-dafny. We provide a formal semantics where
the meaning of an instruction is given as a partial function that maps states
to states. Our semantics is language-agnostic, readable and can be used as a
reference for developers of EVMs or to aid compiler writers. Moreover, it is a
complete and usable framework for formally reasoning about correctness of EVM
bytecode using Dafny.

##  2 Background & Motivation

In this section we give an overview of the EVM and show how our formal
specification in Dafny can be used to verify properties of bytecode programs.

The Ethereum blockchain stores the bytecode of the _contracts_ into a database
and each contract has its own permanent storage. In what follows, we assume a
given contract and refer to storage as that allocated to this contract.

### Instructions and States.

The EVM [28] is a _stack-based_ machine [28] which supports 142 instructions:
arithmetic operations (e.g., ADD, MUL), comparisons and bitwise operations
(e.g., ISZERO, NOT), cryptographic primitives (e.g., SHA3), environment
information (e.g., BALANCE, CALLVALUE), block information (e.g., NUMBER,
GASLIMIT), stack/memory/control flow (e.g., PUSH, POP, MSTORE, SLOAD, JUMP),
logging (e.g., LOG1), and system operations (e.g., CREATE, CALL,
DELEGATECALL). An _executing state_ of the EVM is a tuple containing several
components. We restrict our attention to the following subset of these
components:

code:

a sequence of $n$ bytes indexed from $0$ to $n-1$; The byte at index $0\leq
k<n$ is either an instruction opcode or an immediate operand. For instance the
sequence of opcodes given by $s=[{\tt 0x60,0x01,0x60,0x02,0x01,0x50,0x00}]$
corresponds to the more readable _program_ “PUSH1 0x01; PUSH1 0x02;ADD; POP;
STOP”. Here, the byte at $s[1]$ (${\tt 0x01}$) is the operand of the
instruction at $s[0]$ (PUSH1).

pc:

the _program counter_ (initially $0$) identifies the next instruction to
execute. For example, if pc is $4$, executing the instruction at $s[4]$ (ADD)
increments it by $1$ so $s[5]$ is the next instruction to execute. When
executing instructions with operands (e.g., “PUSH1 0x01” at $s[0]s[1]$) the pc
is incremented by $1+v$ where $v$ is the number of operands.

stack:

a _stack_ of 256-bit words (initially empty); instructions can push or pop the
stack. For example, starting from an empty stack $[\,]$, executing the
instructions “PUSH1 0x1; PUSH1 0x2” gives $[{\tt 0x02,0x01}]$. Executing the
ADD instruction from the stack $[{\tt 0x02,0x01}]$ pops 2 operands, adds them
and pushes the result yielding a new stack $[\tt 0x03=0x01+0x02]$.

memory:

a 256-bit addressable, contiguous array of bytes (initially empty). Memory is
volatile and only available during the current program execution. Memory
expands on-demand when a value is read or written to a given location (which
incurs some cost in gas). Values can be read from/written to memory using the
instructions MLOAD, MSTORE, MLOAD8 or MSTORE8.

storage:

a map from 256-bit addresses to 256-bit words which constitutes the contract’s
permanent storage. Storage can be read/written using the instructions SLOAD or
SSTORE.

gas:

the fuel left for future computations. Executing an instruction consumes _gas_
in the EVM, and this ensures that no infinite computation can occur.

In the EVM, program execution may _abort_ under exceptional cases including:

Out-of-gas:

the gas left in the current state does not cover the cost of executing the
next instruction (including cost of memory expansion if any);

Stack exceptions:

the stack size cannot exceed $1024$. Moreover, some instructions (e.g., POP)
can only be executed if the stack has enough elements and otherwise the
execution should abort.

The EVM has _failure states_ to capture aborted computations. As a result, a
_state_ of the EVM is either a failure state or a non-failure state.

### Bytecode Verification.

Using our formal semantics, we can guarantee security properties of bytecode
programs using the Dafny verifier. Dafny is a verification-friendly language
and as such the code can be instrumented with predicates and pre- and
postconditions that are checked by the verifier at _compile time_. We use this
feature to prove properties on the bytecode. The following simple Dafny
program illustrates a proof:

⬇

1method AddBytes(x: u8, y: u8) {

2 // Initialise an EVM with some gas and the bytecode to execute.

3 var st := InitEmpty(gas:=1000, code:=[PUSH1,x,PUSH1,y,ADD]);

4 // Execute 3 compute steps

5 st := ExecuteN(st,3);

6 // Check that the top of the stack is the sum of x and y

7 assert st.Peek(0) == (x as u256) + (y as u256);

8}

This simple code snippet illustrates several aspects of the verification
process. First we can verify _family of programs_ as the parameters x,y are
arbitrary unsigned integers over 8 bits. This is done by creating an EVM and
stepping through the code, e.g., using the ExecuteN function. Second, we
specify the expected property of the code using the assert statement (line 7)
which is a _verification_ statement: it is not executed at runtime as in
conventional programming languages but _checked_ at _compile-time_ , and must
hold for all inputs. For this program Dafny can prove automatically that the
assert statement is never violated. The proof uses the semantics of opcodes
that are invoked in the computation of ExecuteN. Note that if we change u8 to
u256 the property does not hold as an overflow can occur in the execution of
ADD: this is flagged by the Dafny verifier with “Cannot prove assertion at
line 7”. Another set of checks that are performed automatically are related to
pre- and postconditions. For instance the ADD instructions requires at least
two elements on the stack. This is specified by a precondition in the function
that defines the semantics of ADD. If the code above had only one PUSH1
instruction Dafny would flag that the ADD cannot be performed as a
precondition is violated. Overall, this short code snippet demonstrates that
we can specify and verify functional correctness properties of bytecode, and
thanks to the pre- and postconditions used to specify the semantics of the
instructions, we can detect/fix possible exceptions (e.g., stack overflow)
before runtime.

The example in listing A.1 shows how we can reason about storage updates and
exceptions (aborted computations).

This contract code maintains a counter at storage location $0$ which is
incremented by one on every contract call. Initially, the contract storage is
unconstrained in the input state st and, hence, any location can contain any
value. The code of the contract aims to capture overflows and to revert if an
overflow occurs. The intent is that either the contract reverts (overflow
detected) or the counter is incremented by $1$. Listing A.1 gives a Dafny
proof of this.222The code in the paper may not compile or verify as we have
simplified it for clarity. The code in https://github.com/ConsenSys/evm-dafny
compiles and verifies. The preconditions (lines 10–11) ensure that st is an
execution (non-failure) state with pc == 0, empty stack, enough gas, and has
the contract code to execute.

The postconditions (lines 12-14) specify that the computation either
increments the counter (at storage location $0$) or the computation reverts.
The proof divides up into two essential parts:

  1. 1.

Execute the first 7 bytecodes and store the intermediate state in nst.

  2. 2.

An overflow occurs when the result of the addition is $0$. So depending on the
result at the top of the stack, nst.Peek(0), we decide whether the rest of the
computation will either succeed or revert.

Dafny successfully verifies this code and guarantees the postconditions on
lines 12-14 for all input states st satisfying the preconditions (lines
10-11). This provides strong guarantees about the bytecode: ($i$) it either
reverts or computes the increment but never runs out of gas, nor ends up in an
invalid state (e.g., stack overflow or underflow), ($ii$) the program
terminates normally _if and only if_ the initial value stored at location $0$
is strictly less than MAX_U256 (line 13), ($iii$) on normal termination, the
value in storage location $0$ is incremented by one (line 14).

Listing A.1: Verifying bytecode with Reverts.

⬇

1const INC_CONTRACT := Code.Create([

2 // Put STORAGE[0] on stack and increment by one

3 PUSH1, 0x0, SLOAD, PUSH1, 1, ADD,

4 // If result non-zero branch to JUMPDEST, else REVERT

5 DUP1, PUSH1, 0xf, JUMPI, PUSH1, 0x0, PUSH1, 0x0, REVERT,

6 // Write result back to STORAGE[0] and return

7 JUMPDEST, PUSH1, 0x0, SSTORE, STOP]);

8

9method IncProof(st: State) returns (st’: State)

10 requires st.OK? && st.PC() == 0 && st.Gas() >= 40000 ...

11 requires st.evm.code == INC_CONTRACT

12 ensures st’.REVERTS? || st’.RETURNS?

13 ensures st’.RETURNS? <==> (st.Load(0) as nat) < MAX_U256

14 ensures st’.RETURNS? ==> st’.Load(0) == (st.Load(0) + 1) {

15 // Execute upto (and including) JUMPI.

16 var nst := ExecuteN(st,7);

17 // Consider branches separately

18 if nst.Peek(0) == 0 { // test top of the stack

19 assert nst.PC() == 0xa;

20 nst := ExecuteN(nst,3);

21 assert nst.REVERTS?;

22 } else {

23 assert nst.PC() == 0xf;

24 nst := ExecuteN(nst,4);

25 assert nst.RETURNS?;

26 }

27 return nst;

28}

##  3 The Dafny-EVM

Our EVM is written in Dafny and provides a definition of the semantics as a
function mapping states to states. A key design decision made early on was to
develop a functionally pure formalisation of the EVM. In this section we
describe the main components of the Dafny-EVM and conclude with some
observations.

### Machine State.

Line numbers hereafter refer to Listing A.3. A state of the EVM is a record
containing various fields such as gas, pc, stack, code, memory.

Each module (state, stack, memory, …) provides a datatype, possibly
incorporating some contraints (e.g., EvmState.T). For brevity, we omit some
fields which contain information about the enclosing transaction and the so-
called substate. The State datatype (line 7) models normal execution (OK),
failure (INVALID), returning (RETURNS), reverting (REVERTS), etc.

Figure 1: Source files of the Dafny-EVM. Top group contains bytecode semantics
and top-level types. Middle group contains abstractions of the main
components. Bottom group are fundamental primitives (e.g. for manipulating
bytes and ints). “Loc” (lines of codes) at the time of writing.

### Stack, Memory and Storage.

We have implemented several submodules to provide operations on
stack/memory/storage. This is summarised in Figure 1. We lift the operations
on stack/memory/storage into the State datatype. In Dafny this is done by
adding the functions right after the definition of a datatype (line 11). This
allows us to compose them easily and improves readability. For instance the
Add function that implements the semantics of opcode ADD is defined using a
sequence of operations st.Pop().Pop().Push(...).Next() where st is an
executing state (e.g. OK). We employ preconditions (requires) to ensure lifted
operations are limited to applicable states only (typically executing states,
such as OK), and also that preconditions of the functions on
stack/memory/storage are satisfied (e.g., for Pop() the stack size must be
large enough); for Push() (line 20) the stack cannot be full (stack size is
limited to 1024).

In Dafny, preconditions are checked by the verifier and must provably hold at
each call site. Notice that Dafny enforces the constraints on integer types so
every time we compute (e.g., ADD) and store the result in a 256bit word, we
must prove that the value is less than $2^{256}$ (the EVM dictates modulo
arithmetic for this). The pre-/post-conditions and type checks enforced by the
Dafny verifier help ensure that our EVM specification is consistent and that
functions are well-defined.

Listing A.2: Semantics of MLOAD, Bytecode module

⬇

1function method MLoad(st: State) : State

2 requires st.IsExecuting() {

3if st.Operands() >= 1 then

4 var loc := st.Peek(0) as nat;

5 var nst := st.Expand(loc,32); // Break out expanded state

6 nst.Pop().Push(nst.Read(loc)).Next() // Read from expanded state

7else

8 State.INVALID(STACK_UNDERFLOW)

9}

Memory operations are provided by the Memory module, with various functions
being attached to State, e.g., Read, Write lines 26–28. A key observation is
that, in both cases, address addr + 31 must be within allocated memory. This
is because memory in the EVM is byte addressable and we are reading/writing
u256 values (i.e., which are 32 bytes long). The semantics of MLOAD (Listing
A.2) highlights the complexity of memory operations. Since Read(loc) (line 6)
has the precondition loc + 31 < Memory.Size (line 26 of Listing A.3), this
must hold for state nst. In fact, this follows because the call to Expand()
(line 5) ensures sufficient memory. If the call to Expand() within MLoad was
not enforcing this constraint, then Dafny would raise a precondition violation
on nst.Read(loc).

### Gas.

In our design, we chose to split out the gas calculation from the instruction
semantics. Whilst this does introduce some repetition, we argue it reduces
cognitive load. In particular, since this avoids interweaving the gas
calculation throughout the instruction semantics which (for performance
reasons) is commonly done in actual implementations (including the _execution
specs_ 333https://github.com/ethereum/execution-specs).

### Contract Calls.

Various instructions (e.g. CALL, DELEGATECALL) enable one contract to call
another. These differ from others as they can involve executing arbitrarily
many instructions in the called contract. We implement this using a mechanism
akin to continuations but, for brevity, omit the details here.

Listing A.3: The EvmState module (partial)

⬇

1module EvmState {

2 datatype Raw = EVM(gas:nat, pc:nat, stack:Stack.T, code:Code.T,

3 mem:Memory.T, world:WorldState.T, ...)

4

5 type T = c:Raw | c.context.address in c.world.accounts

6

7 datatype State = OK(evm:T) | REVERTS(gas:nat,data:seq<u8>)

8 | RETURNS(gas:nat,data:seq<u8>,...) | INVALID(Error) | ...

9 {

10 // Predicates

11 predicate method IsExecuting(): bool { ... }

12

13 // Stack functions

14 function method Capacity(): nat

15 requires IsExecuting() { Stack.Capacity(evm.stack) }

16 function method Peek(k: nat): u256

17 requires IsExecuting() && k < Stack.Size(evm.stack) { ... }

18 function method Pop(): State

19 requires IsExecuting() && 0 < Stack.Size(evm.stack) { ... }

20 function method Push(v: u256) : State

21 requires IsExecuting()

22 requires Capacity() > 0 {

23 OK(evm.(stack:=Stack.Push(evm.stack,v)))

24 }

25 // Memory functions

26 function method Read(address: nat): u256

27 requires IsExecuting() && (addr+31) < Memory.Size(evm.mem) {...}

28 function method Write(address: nat, val: u256): State

29 requires IsExecuting() && (addr+31) < Memory.Size(evm.mem) {...}

30 ...

31 function method Expand(addr: nat, n: nat): (s’: State)

32 requires IsExecuting()

33 ensures s’.IsExecuting() && MemSize() <= s’.MemSize()

34 ensures (addr + n) < MemSize() ==> (evm.mem == s’.evm.mem) {...}

35 }

36 ...

37}

### Observations.

The Dafny-EVM Code provides a readable and executable specification of the
EVM. There are several benefits of using a verification-friendly language:
using pre- and postconditions to write the semantics provides a high level of
assurance; furthermore, the code is executable and can be compiled into
several target languages including Java, C#, Go. We now highlight some
observations based on our experiences from this project.

  * •

Specification. Dafny treats function calls within expressions as interpreted,
but treats method calls as uninterpreted [5, 15]. Roughly speaking this means
that, when verifying a function call, the verifier has free access to the
function’s body. In contrast, for method calls, the verifier can only access
what is given in the specification (i.e. its pre- and postconditions). As
such, we consider methods ill-suited for formalising specifications (such as
for the EVM). This is because we cannot abstract a specification any further
than already done (i.e. we cannot specify a specification).

  * •

Verification. Functions can have preconditions that restrict the domain of
their inputs. In Dafny preconditions are enforced at each call site. We argue
that this results in better code by enforcing consistency across function
calls. Dafny enforces that every function must have a proof of termination
which guarantees the absence of infinite loops in our state transition
function. We believe that this degree of assurance is hard to attain with non
verification-friendly languages.

  * •

Performance. Code generated from the functionally pure subset of Dafny can
perform poorly because of the need to clone compound structures (e.g. maps and
arrays) to preserve purity (i.e., referential transparency). Dafny does not,
for example, employ clone elimination [16, 26, 19] or mutable value semantics
[22, 21]. Performance was not a critical concern given our aim of developing a
formal specification rather than an efficient implementation and in practice,
we did not encounter any significant issue here.

During the project, a number of issues and challenges arose. For example, the
lack of an exponentiation operator in Dafny meant that, for the EXP bytecode,
we had to implement this as a recursive function. Some low level operations
involving bits & bytes (e.g., shifting) present significant challenges as the
native int type does not support bitwise operators. One can use a conversion
from (e.g. u256) into the bitvector types (e.g. bv256) provided by Dafny which
do support bitwise operations — however, this can lead to problems verifying
code.

##  4 Practical Experiences

From the outset of this project, we were unsure whether Dafny would be
practical for this sizeable formalisation task. Overall, however, we are
pleased to report that Dafny has, for the most part, proven itself more than
capable. Of course, it was not all plain sailing and we encountered several
challenges which required developing techniques and/or workarounds.

### Code Generation.

Dafny can generate code for a variety of targets, including: C#, Go, Java,
C++, Python and JavaScript. Furthermore, whilst Dafny does not support I/O
operations per se, these can be implemented on the target side. We took
advantage of this to embed the Dafny-generated code into a thin Java wrapper
that performs I/O and allows us to test our EVM against existing
implementations. Note that the generated code is not proved to be equivalent
to the original Dafny code. For various reasons (e.g., knowledge within the
team) we chose Java as the target language with gradle managing the build.
This worked well enough, though there are some points to make:

  * •

Foreign Function Interface. Code generated from Dafny does not conform to the
stylistic norms of Java, but is otherwise relatively easy to interface with. A
runtime library is provided by Dafny against which generated code must be
compiled. This provides (amongst other things) alternative collection
implementations (e.g. DafnySequence, DafnyMap, etc).

  * •

External Code. For the semantics of KECCAK256 and some precompiled contracts,
we preferred to call out to native Java code (i.e. rather than implement e.g.
sha256 in Dafny itself). However, whilst Dafny does support extern
declarations, these are not (at the time of writing) well supported by the
Java code generator. Instead, we had to give default implementations (e.g.
returning 0) and employ build trickery to make it work.

  * •

Target language idiosyncrasies. Translation to a target language introduces
risks. E.g., Dafny employs Euclidean Division for its integer division
operator (i.e. always rounds down rather than towards zero), which is a trap
for the unwary and by chance we identified a bug in the Java code generator
where sometimes standard division was being
applied.444https://github.com/dafny-lang/dafny/issues/2367 We also encountered
unsoundness in the translation of Dafny collections (e.g. seq<u8>) to
Java555https://github.com/dafny-lang/dafny/issues/2859, and buggy
implementation of datatype in C#. 666https://github.com/dafny-
lang/dafny/issues/1412

### Verification and Testing.

For completeness, we developed many unit tests for various components of our
formalism. The Ethereum Common Tests also provide tens of thousands of tests
for ensuring EVM compatibility.777https://github.com/ethereum/tests As such,
we have been using these to check our formalisation against existing
implementations. This required generating executable code from our
specification which presented several challenges (discussion of which is
unfortunately omitted for brevity). At the time of writing, we have selected
around 7500 representative tests out of the 13K Common tests (Berlin hardfork)
and 6900 are passing ($92\%$). Of the 143 failing tests, the majority (100)
are failing because: some precompiled contracts are not yet fully implemented
(44); we do not currently check for branches into instruction operands (56).
The remaining (approx. 450) tests are skipped for various reasons e.g.,
timeout or breaking the testing system. Finally, we note that all of our tests
are run as part of Continuous Integration before a pull request can be merged.

##  5 Related Work

Initial attempt at a formal specification of the EVM may be attributed to
Hirai [13] with a formalisation of the EVM in the programming development
environment Lem [18]. The formalisation in [13] is restricted to a single
contract execution and proving bytecode is limited in terms of automation.
Later, Amani  _et al._ [3] built upon Hirai’s formalisation and proposed an
Isabelle/HOL formalisation. Their contribution introduces a program logic to
reason about bytecode (restricted to a subset of 36/142 EVM instructions) but
they rely on the construction of a control flow graph to define the semantics
of a program. Reasoning about bytecode is limited to linear sequences of
instructions (blocks) and not fully automated. Another Isabelle/HOL
specification was also developed in [9] specialised for gas consumption
analysis and for proving termination of bytecode.

More recently, Grishchenko  _et al._ [10] have proposed a partial (not all
opcodes are supported and the gas cost semamtics is incomplete) formalisation
of the EVM in $F^{\star}$ targeting verification of security properties.

The most advanced formalisation is probably the KEVM [12] using the
$\mathbb{K}$ Framework [23]. It provides a formal and executable specification
of the syntax and semantics of EVM bytecode. Using the built-in automated
tools of the $\mathbb{K}$ Framework, it is possible to generate an
interpreter, compiler, debugger and to some extent a verifier that can be used
to check the bytecode of some contracts [20]. The default input format (used
for KEVM) of the $\mathbb{K}$ Framework is XML-based which may not be the most
developer-friendly format. Similarly, IELE [24] attempts to design a more
readable language than EVM bytecode and to be the target of high-level
languages including Solidity, Vyper, Plutus. IELE is defined using the
$\mathbb{K}$ Framework and uses LLVM tools (compiler) as a backend.

There are several implementations of the EVM in different languages and
clients e.g., Geth888https://geth.ethereum.org,
Besu999https://github.com/hyperledger/besu, and more recently the _execution-
specs_ in Python3. The implementations in Geth and Besu are respectively in Go
and Java and cannot be used to reason about bytecode. The Python
implementation relies on specific imperative language features of Python
(mutability, exceptions) and does not provide a functional definition of the
instructions semantics nor an explicit specification of exceptional cases: for
instance the Python code does not provide preconditions or explicit handling
of exceptions, and exceptions can happen deep in the call stack which may
hinder readability.

There are several tools Oyente [4], EtherIR [1], eThor [25], Rattle [27], and
Certora [14] to perform static analysis of EVM bytecode. There are also
extensions to specifically analyse the gas consumption like GASTAP [2],
GasReducer [8]. Those tools build an abstract representation of the bytecode
and it is unclear whether the abstraction is semantics preserving.

In contrast to the formalisations, implementations and tools referenced above,
our formal semantics is language-agnostic (defines the state transition
function as a function), easy to read and developer-friendly, provides
mathematical and verified pre- and postconditions for the semantics of
instructions. Moroever, our semantics can be used to perform _deductive
reasoning_ about bytecode including gas consumption using standard invariants.

##  6 Conclusion

We have proposed a formal semantics of the EVM in a pure functional subset of
Dafny. Our semantics is human readable, machine checked and executable, and
provides a sound framework to formally reason about bytecode.

This opens up the door for several direct applications:101010Examples are
available in https://github.com/ConsenSys/evm-dafny.

  * •

complete smart contract verification: in practice, this can be a costly
process and may require specific verification skills or familiarity with
Dafny.

  * •

correctness of compiler optimisations: several gas optimisation patterns e.g.,
a sequence of opcodes SWAP1 POP POP optimised in POP POP can now be verified.

  * •

correctness of under/overflow detection: to detect an overflow in arithmetic
modulo ADD(x, y) it is common to first compute the result r = ADD(x, y) and
then check that r >= x. We can formally prove that this is sound.

  * •

synthesise verified bytecode: we have designed a methodology [6] to specify
and verify smart contracts directly in Dafny. We are exploring _refinement
proof techniques_ to synthesise bytecode from the verified Dafny code of a
contract. Ultimately we may develop a Dafny-to-EVM _certified compiler_.

Although the benefits of our approach are evident in the formal methods’
community, adoption of these techniques in the Ethereum ecosystem is still
challenging. Whilst established techniques, e.g., using Solidity to write
contracts, or using Python to write specifications, can be questionable [11],
they are still prevalent in the Ethereum community. The main hurdles for
mainstream adoption of our approach are probably two-fold: ($i$) provide
developer-friendly tools to write contracts; Dafny and the tool support around
it (e.g., verification performance improvement, counter example generation
[7], VSCode integration) already partially solves this issue; and ($ii$)
educate the Ethereum community to understand the long-term benefits of formal
verification for the Ethereum ecosystem.

## References

  * [1] Albert, E., Gordillo, P., Livshits, B., Rubio, A., Sergey, I.: Ethir: A framework for high-level analysis of ethereum bytecode. In: Lahiri, S.K., Wang, C. (eds.) Automated Technology for Verification and Analysis. pp. 513–520. Springer International Publishing, Cham (2018)
  * [2] Albert, E., Gordillo, P., Rubio, A., Sergey, I.: Running on fumes - preventing out-of-gas vulnerabilities in ethereum smart contracts using static resource analysis. In: Ganty, P., Kaâniche, M. (eds.) Verification and Evaluation of Computer and Communication Systems - 13th International Conference, VECoS 2019, Porto, Portugal, October 9, 2019, Proceedings. Lecture Notes in Computer Science, vol. 11847, pp. 63–78. Springer (2019). https://doi.org/10.1007/978-3-030-35092-5_5
  * [3] Amani, S., Bégel, M., Bortin, M., Staples, M.: Towards verifying ethereum smart contract bytecode in Isabelle/HOL. In: Andronick, J., Felty, A.P. (eds.) Proceedings of the 7th ACM SIGPLAN International Conference on Certified Programs and Proofs, CPP 2018, Los Angeles, CA, USA, January 8-9, 2018\. pp. 66–77. ACM (2018). https://doi.org/10.1145/3167084
  * [4] Badruddoja, S., Dantu, R., He, Y., Upadhayay, K., Thompson, M.: Making smart contracts smarter. In: 2021 IEEE International Conference on Blockchain and Cryptocurrency (ICBC). pp. 1–3 (2021). https://doi.org/10.1109/ICBC51069.2021.9461148
  * [5] Bradley, A.R., Manna, Z.: The calculus of computation - decision procedures with applications to verification. Springer (2007)
  * [6] Cassez, F., Fuller, J., Anton Quiles, H.M.: Deductive verification of smart contracts with Dafny. In: Groote, J.F., Huisman, M. (eds.) Formal Methods for Industrial Critical Systems - 27th International Conference, FMICS 2022, Warsaw, Poland, September 14-15, 2022, Proceedings. Lecture Notes in Computer Science, vol. 13487, pp. 50–66. Springer (2022). https://doi.org/10.1007/978-3-031-15008-1_5
  * [7] Chakarov, A., Fedchin, A., Rakamaric, Z., Rungta, N.: Better counterexamples for dafny. In: Fisman, D., Rosu, G. (eds.) Tools and Algorithms for the Construction and Analysis of Systems - 28th International Conference, TACAS 2022, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2022, Munich, Germany, April 2-7, 2022, Proceedings, Part I. Lecture Notes in Computer Science, vol. 13243, pp. 404–411. Springer (2022). https://doi.org/10.1007/978-3-030-99524-9_23
  * [8] Chen, T., Li, Z., Zhou, H., Chen, J., Luo, X., Li, X., Zhang, X.: Towards saving money in using smart contracts. In: Zisman, A., Apel, S. (eds.) Proceedings of the 40th International Conference on Software Engineering: New Ideas and Emerging Results, ICSE (NIER) 2018, Gothenburg, Sweden, May 27 \- June 03, 2018. pp. 81–84. ACM (2018). https://doi.org/10.1145/3183399.3183420
  * [9] Genet, T., Jensen, T.P., Sauvage, J.: Termination of ethereum’s smart contracts. In: Samarati, P., di Vimercati, S.D.C., Obaidat, M.S., Ben-Othman, J. (eds.) Proceedings of the 17th International Joint Conference on e-Business and Telecommunications, ICETE 2020 - Volume 2: SECRYPT, Lieusaint, Paris, France, July 8-10, 2020. pp. 39–51. ScitePress (2020). https://doi.org/10.5220/0009564100390051
  * [10] Grishchenko, I., Maffei, M., Schneidewind, C.: A semantic framework for the security analysis of ethereum smart contracts. In: Bauer, L., Küsters, R. (eds.) Principles of Security and Trust - 7th International Conference, POST 2018, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2018, Thessaloniki, Greece, April 14-20, 2018, Proceedings. Lecture Notes in Computer Science, vol. 10804, pp. 243–269. Springer (2018). https://doi.org/10.1007/978-3-319-89722-6_10
  * [11] Guido, D.: Episode 6: What the hell are the blockchain people doing, and why isn’t it a dumpster fire? (2021), https://galois.com/blog/2020/11/introducing-the-building-better-systems-podcast/, In Building Better Systems (podcast), Joey Dodds, Shpat Morina, Galois
  * [12] Hildenbrandt, E., Saxena, M., Rodrigues, N., Zhu, X., Daian, P., Guth, D., Moore, B.M., Park, D., Zhang, Y., Stefanescu, A., Rosu, G.: KEVM: A complete formal semantics of the ethereum virtual machine. In: 31st IEEE Computer Security Foundations Symposium, CSF 2018, Oxford, United Kingdom, July 9-12, 2018. pp. 204–217. IEEE Computer Society (2018). https://doi.org/10.1109/CSF.2018.00022
  * [13] Hirai, Y.: Defining the ethereum virtual machine for interactive theorem provers. In: Brenner, M., Rohloff, K., Bonneau, J., Miller, A., Ryan, P.Y.A., Teague, V., Bracciali, A., Sala, M., Pintore, F., Jakobsson, M. (eds.) Financial Cryptography and Data Security - FC 2017 International Workshops, WAHC, BITCOIN, VOTING, WTSC, and TA, Sliema, Malta, April 7, 2017, Revised Selected Papers. Lecture Notes in Computer Science, vol. 10323, pp. 520–535. Springer (2017). https://doi.org/10.1007/978-3-319-70278-0_33
  * [14] Jackson, D., Nandi, C., Sagiv, M.: Certora technology white paper. Medium Post (2022), https://medium.com/certora/certora-technology-white-paper-cae5ab0bdf1
  * [15] Kroening, D., Strichman, O.: Decision Procedures - An Algorithmic Point of View, Second Edition. Springer-Verlag (2016)
  * [16] Lameed, N., Hendren, L.J.: Staged static techniques to efficiently implement array copy semantics in a MATLAB JIT compiler. In: Proceedings of the confererence on Compiler Construction (CC). pp. 22–41 (2011)
  * [17] Leroy, X.: A formally verified compiler back-end. J. Autom. Reason. 43(4), 363–446 (2009). https://doi.org/10.1007/s10817-009-9155-4
  * [18] Mulligan, D.P., Owens, S., Gray, K.E., Ridge, T., Sewell, P.: Lem: reusable engineering of real-world semantics. In: Jeuring, J., Chakravarty, M.M.T. (eds.) Proceedings of the 19th ACM SIGPLAN international conference on Functional programming, Gothenburg, Sweden, September 1-3, 2014. pp. 175–188. ACM (2014). https://doi.org/10.1145/2628136.2628143
  * [19] Odersky, M.: How to make destructive updates less destructive. In: Proceedings of the ACM symposium on the Principles Of Programming Languages (POPL). pp. 25–36 (1991)
  * [20] Park, D., Zhang, Y., Rosu, G.: End-to-end formal verification of ethereum 2.0 deposit smart contract. In: Lahiri, S.K., Wang, C. (eds.) Computer Aided Verification - 32nd International Conference, CAV 2020, Los Angeles, CA, USA, July 21-24, 2020, Proceedings, Part I. Lecture Notes in Computer Science, vol. 12224, pp. 151–164. Springer (2020). https://doi.org/10.1007/978-3-030-53288-8_8
  * [21] Pearce, D.J., Groves, L.: Designing a verifying compiler: Lessons learned from developing Whiley. Science of Computer Programming pp. 191–220 (2015)
  * [22] Racordon, D., Shabalin, D., Zheng, D., Abrahams, D., Saeta, B.: Implementation strategies for mutable value semantics. Journal of Object Technology 21(2) (2022)
  * [23] Rosu, G.: $\mathbb{K}$: A semantic framework for programming languages and formal analysis tools. In: Pretschner, A., Peled, D., Hutzelmann, T. (eds.) Dependable Software Systems Engineering, NATO Science for Peace and Security Series - D: Information and Communication Security, vol. 50, pp. 186–206. IOS Press (2017). https://doi.org/10.3233/978-1-61499-810-5-186
  * [24] Runtime Verification: The iele virtual machine. Blog post (2022), https://runtimeverification.com/the-iele-virtual-machine/
  * [25] Schneidewind, C., Grishchenko, I., Scherer, M., Maffei, M.: ethor: Practical and provably sound static analysis of ethereum smart contracts. In: Ligatti, J., Ou, X., Katz, J., Vigna, G. (eds.) CCS ’20: 2020 ACM SIGSAC Conference on Computer and Communications Security, Virtual Event, USA, November 9-13, 2020. pp. 621–640. ACM (2020). https://doi.org/10.1145/3372297.3417250
  * [26] Shankar, N.: Static analysis for safe destructive updates in a functional language. In: Proceedings of the Symposium Logic-Based Program Synthesis and Transformation (LOPSTR). pp. 1–24 (2001)
  * [27] Trail of Bits: Rattle – an Ethereum EVM binary analysis framework. Medium Post (2018), https://blog.trailofbits.com/2018/09/06/rattle-an-ethereum-evm-binary-analysis-framework/
  * [28] Wood, G.: Ethereum: a secure decentralised generalised transaction ledger. ethereum project yellow paper (2022), https://ethereum.github.io/yellowpaper/paper.pdf, berlin version d77a387 - 2022-04-26