graphs

# A novel event generator for the automated simulation of neutrino scattering

Joshua Isaacson  Stefan Höche  Theoretical Physics Department, Fermi National
Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA  Diego Lopez
Gutierrez  Physics Department, Harvard University, 17 Oxford Street,
Cambridge, MA 02138, USA  Noemi Rocco  Theoretical Physics Department, Fermi
National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA

###### Abstract

An event generation framework is presented that enables the automatic
simulation of events for next-generation neutrino experiments in the Standard
Model or extensions thereof. The new generator combines the calculation of the
leptonic current based on an automated matrix element generator, and the
computation of the hadronic current based on a state-of-the-art nuclear
physics model. The approach is validated in Standard-Model simulations for
electron scattering and neutrino scattering. Furthermore, the first fully-
differential neutrino trident production results are shown in the quasielastic
region.

††preprint: FERMILAB-PUB-21-537-T, MCNET-21-31

##  I Introduction

Neutrino physics is entering an era of precision measurements. The development
of the Deep Underground Neutrino Experiment (DUNE) Abi  _et al._ (2021, 2020)
and of the Hyperkamiokande (HyperK) Abe  _et al._ (2018) detector will allow
us to probe neutrino interactions to unprecedented precision, and will allow
for tests of beyond the Standard Model (BSM) scenarios such as proton decay,
sterile neutrinos and heavy neutral leptons, as well as rare Standard Model
(SM) processes such as neutrino tridents Abi  _et al._ (2021); Yano (2021).
The excess of electron-like events in measurements at LSND Aguilar  _et al._
(2001) and MiniBooNE Aguilar-Arevalo  _et al._ (2007, 2009, 2010, 2013, 2018,
2021) has demonstrated the need for a toolkit which allows the physics
community to quickly test different hypotheses against experimental data.
While calculations can be done on a case-by-case basis, such as in Refs.
Gninenko (2009); McKeen and Pospelov (2010); Gninenko (2011); Dib  _et al._
(2011); Gninenko (2012); Masip  _et al._ (2013); Ballett  _et al._ (2017);
Magill  _et al._ (2018); Fischer  _et al._ (2020); Bertuzzo  _et al._ (2018a,
b); Ballett  _et al._ (2019a); Argüelles  _et al._ (2019); Ballett  _et al._
(2019b, 2020); Abdullahi  _et al._ (2020); Abdallah  _et al._ (2020a, b); de
Gouvêa  _et al._ (2020); Dentler  _et al._ (2020); Chang  _et al._ (2021), a
dedicated analysis for each new physics model becomes impractical when there
is a wealth of ideas and a large number of measurements that present
simultaneous constraints, such as the recent MicroBooNE results Abratenko . In
such cases an end-to-end simulation should be attempted, starting with the
Lagrangian of the underlying new physics hypothesis, and leading to particle-
level events that are generated fully differentially in the many-particle
phase space by means of Monte-Carlo methods. In addition, it must be possible
to implement this simulation pipeline as part of experimental analysis
frameworks, such that detector effects can be fully included.

The accurate simulation of particle-level events based on an underlying
physics model has been a cornerstone of high-energy physics experiments for
many decades. In the context of Large Hadron Collider (LHC) physics, the
problem of an end-to-end simulation pipeline has been encountered long ago,
and it was successfully addressed through the development of FeynRules
Christensen and Duhr (2009); Alloul  _et al._ (2014), specific event generator
interfaces Christensen  _et al._ (2011, 2012) and eventually the UFO file
format Degrande  _et al._ (2012). A complete toolchain has been developed,
starting with the conversion of the Lagrangian to Feynman rules Christensen
and Duhr (2009); Alloul  _et al._ (2014), followed by the perturbative
calculation of cross sections and event simulation by Amegic Krauss  _et al._
(2002), Comix Gleisberg and Höche (2008), Herwig++ Bähr  _et al._ (2008);
Bellm  _et al._ (2016), MadGraph Alwall  _et al._ (2014) or Whizard Kilian
_et al._ (2007), and finally leading to particle-level event simulation with
Herwig Bähr  _et al._ (2008); Bellm  _et al._ (2016), Pythia Sjöstrand  _et
al._ (2006, 2015) or Sherpa Gleisberg  _et al._ (2009); Bothmann  _et al._
(2019). The flexibility of the aforementioned tools plays a major role in the
exploration of model and parameter space by the LHC experiments. In
anticipation of similar needs at next generation neutrino experiments, we
develop a simulation framework that allows for the generation of particle-
level events in arbitrary new physics models, while at the same time
appropriately including nuclear effects. We validate our new framework in
Standard Model electron- and neutrino-nucleus scattering. We also calculate
the first results for fully differential neutrino trident production in the
quasielastic region.

This manuscript is structured as follows: Section II presents an introduction
to the problem of neutrino-nucleus scattering and lays out the general
framework for the calculation. Section III introduces our technique to compute
the leptonic tensor, and Sec. IV reviews the techniques for phase-space
integration. Section V outlines the computation of the hadronic tensor. First
numerical results are presented in Sec. VI and Sec. VII contains an outlook.

##  II Theory Overview

In neutrino experiments, the majority of events of interest can be described
via an exchange of a single gauge boson. In this case, the differential cross
section factorizes as:

$\frac{{\rm d}\sigma}{{\rm d}\Omega}\propto L_{\mu\nu}W^{\mu\nu},$ (1)

where $L_{\mu\nu}\,(W^{\mu\nu})$ is the leptonic (hadronic) tensor.

The leptonic tensor can be calculated using methods developed for collider
event generation and will be discussed in detail in Sec. III. On the other
hand, when dropping spin dependent terms, the hadronic tensor can be written
using the most general Lorentz structure as:

$W^{\mu\nu}=\left(-g^{\mu\nu}+\frac{q^{\mu}q^{\nu}}{q^{2}}\right)W_{1}+\frac{\hat{p}_{a}^{\mu}\hat{p}_{a}^{\nu}}{p_{a}\cdot
q}W_{2}-i\epsilon^{\mu\nu\alpha\beta}\frac{q_{\alpha}p_{a\beta}}{2p_{a}\cdot
q}W_{3},$ (2)

where $p_{a}$ is the momentum of the initial nucleon, $q$ is the momentum of
the probe, $\hat{p}_{a}^{\mu}=p_{a}^{\mu}-\frac{p_{a}\cdot q}{q^{2}}q^{\mu}$,
$Q^{2}=-q^{2}$, and $W_{i}$ are the nuclear structure functions.

Factorizing the differential cross section in this manner allows us to
separate the BSM effects that would be dominant in the leptonic tensor from
the complicated nuclear physics calculations of the hadronic tensor. Moreover,
neutrino event generators have already implemented many nuclear models into
their codes Juszczak  _et al._ (2006); Andreopoulos (2009); Hayato (2009);
Buss  _et al._ (2012), making the calculation of $W^{\mu\nu}$ a
straightforward process. Therefore, we can leverage the work done on the
nuclear side to make accurate predictions for different BSM scenarios at DUNE,
HyperK, MicroBooNE, and MiniBooNE.

While Eq. (1) works when a single type of gauge boson exchange dominates the
total cross section (i.e. the photon for electron scattering, the $W$ boson
for charged-current neutrino scattering, and the $Z$ boson for neutral-current
neutrino scattering), this approximation breaks down when interference is
important. In this case, the differential cross section can be obtained
through an extended factorization formula that reads

$\frac{{\rm d}\sigma}{{\rm
d}\Omega}=\sum_{i,j}L^{(ij)}_{\mu\nu}W^{(ij){\mu\nu}},$ (3)

where the sum over $i,j$ is for each allowed boson in the process. For
example, adding in the $Z$ boson contributions to electron scattering would
result in

$\frac{{\rm d}\sigma}{{\rm
d}\Omega}=L^{(\gamma\gamma)}_{\mu\nu}W^{(\gamma\gamma){\mu\nu}}+L^{(\gamma
Z)}_{\mu\nu}W^{(\gamma
Z){\mu\nu}}+L^{(Z\gamma)}_{\mu\nu}W^{(Z\gamma){\mu\nu}}+L^{(ZZ)}_{\mu\nu}W^{(ZZ){\mu\nu}}.$
(4)

This formulation of the problem quickly becomes unwieldy as the number of
allowed bosons in the process increases. In these situations, it is better to
write the differential cross section as the square of a product of leptonic
and hadronic currents:

$\frac{{\rm d}\sigma}{{\rm
d}\Omega}=\bigg{|}\sum_{i}L^{(i)}_{\mu}W^{(i)\mu}\bigg{|}^{2},$ (5)

where the interference terms are automatically included by taking the square
of the sum of amplitudes. In this work, we will use Eq. (5) for our
calculations, but in the cases when the nuclear physics cannot be written as a
current, we provide the option to work with the leptonic tensor in Eq. (3).

The BSM calculations we intend to automate involve only the leptonic current
and are independent of the exact nuclear model. We will therefore focus on the
quasielastic region in this work for concreteness. Specifically, the hadronic
current will be computed in the impulse approximation (IA), using the spectral
function formalism as discussed in Refs. Benhar  _et al._ (2008); Rocco  _et
al._ (2019a). In this approximation, the nuclear current can be expressed as a
sum of one-body currents, and the nuclear final state can be written as a
struck nucleon and an $A-1$ spectator nuclear remnant. Given the initial
nuclear ground state $\lvert\psi_{0}^{A}\rangle$ and the nuclear final state
$\langle\psi_{f}^{A}\rvert$, the current operator can be rewritten through the
insertion of a complete set of states ($\sum_{k}\lvert k\rangle\langle
k\rvert$) in the IA as:

$W^{\mu}=\langle\psi_{f}^{A}\rvert\mathcal{J}^{\mu}\lvert\psi_{0}^{A}\rangle\rightarrow\sum_{k}\left[\langle\psi_{f}^{A-1}\rvert\otimes\langle
k\rvert\right]\rvert\psi_{0}^{A}\rangle\langle
p_{a}+q\rvert\sum_{i}\mathcal{J}_{i}^{\mu}\lvert k\rangle,$ (6)

where $\mathcal{J}^{\mu}$ is the current operator of the boson probing the
nucleus, and $\mathcal{J}_{i}^{\mu}$ is the current operator contribution from
the $i^{\rm th}$ nucleon in the IA. The incoherent contribution to the
hadronic tensor can then be calculated as:

$W^{\mu\nu}=\int\frac{{\rm d}^{3}p_{a}}{(2\pi)^{3}}{\rm
d}E\,S(\mathbf{p_{a}},E_{r})\sum_{i}\langle
p_{a}\rvert\mathcal{J}_{i}^{\mu}\lvert p_{a}+q\rangle\langle
p_{a}+q\rvert\mathcal{J}_{i}^{\nu\dagger}\lvert k\rangle\delta\left(\omega-
E_{r}+m-E^{\prime}\right),$ (7)

where $\omega$ is the lepton energy loss, $E_{r}$ is the excitation energy of
the $A-1$ nucleon system, $E^{\prime}$ is the outgoing nucleon energy and
$S(\mathbf{p}_{a},E_{r})$ is the spectral function which gives the probability
of removing a nucleon with three-momentum $\mathbf{p}_{a}$ and excitation
energy $E_{r}$. The spectral functions of finite nuclei have been obtained
using different nuclear many-body approaches Dickhoff and Charity (2019);
Dickhoff and Atkinson (2020); Barbieri and Hjorth-Jensen (2009); Benhar  _et
al._ (1994); Sick  _et al._ (1994). In particular, those computed within the
correlated-basis function (CBF) theory, Self Consistent Green’s Function
(SCGF), and variational Monte Carlo (VMC) methods have yielded consistent
inclusive lepton-scattering cross sections on different nuclear targets Rocco
_et al._ (2019a); Rocco and Barbieri (2018); Andreoli  _et al._ (2021);
Ankowski  _et al._ (2015). The CBF spectral function used in this work is
given as a sum of two terms. The first term, which describes the low
excitation energy and momentum region, uses as input spectroscopic factors
extracted from ($e,e^{\prime}p$) scattering measurements. The second term
accounts for unbound states of the $A-1$ spectator system in which at least
one of the spectator nucleons is in a continuum state. This contribution is
obtained by folding the correlation component of the nuclear matter spectral
function obtained within the CBF theory with the the density profile of the
nucleus Benhar  _et al._ (1989).

The nuclear one-body current operator can be expressed in general as the sum
of the vector (V) and axial-vector (A) terms given as:

$\displaystyle\mathcal{J}^{\mu}$
$\displaystyle=\left(\mathcal{J}_{V}^{\mu}+\mathcal{J}_{A}^{\mu}\right)$
$\displaystyle\mathcal{J}^{\mu}_{V}$
$\displaystyle=\gamma^{\mu}\mathcal{F}^{a}_{1}+i\sigma^{\mu\nu}q_{\nu}\frac{\mathcal{F}^{a}_{2}}{2M}$
(8) $\displaystyle\mathcal{J}^{\mu}_{A}$
$\displaystyle=-\gamma^{\mu}\gamma_{5}\mathcal{F}^{a}_{A}-q^{\mu}\gamma_{5}\frac{\mathcal{F}^{a}_{P}}{M},$

where $\mathcal{F}^{a}_{i}$ are the process dependent nuclear form factors for
the exchange of the boson $a$ and $M$ is the nucleon mass. Since the pseudo-
scalar form factor ($\mathcal{F}^{a}_{P}$) is multiplied by the mass of the
outgoing lepton in the cross section, we neglect the contribution from it in
this work. The values for the form factors considered here are detailed in
Sec. V.

##  III Calculation of the Leptonic Current

Employing a dedicated version of the general-purpose event generator Sherpa
Gleisberg  _et al._ (2004, 2009); Bothmann  _et al._ (2019), we construct an
interface to the Comix matrix element generator Gleisberg and Höche (2008) to
extract the leptonic current. In Comix, the matrix element is computed using
the color-dressed Berends-Giele recursive relations Duhr  _et al._ (2006),
which can be understood as a Dyson-Schwinger based technique Kanaki and
Papadopoulos (2000); Mangano  _et al._ (2003). In this technique, the full
tree-level scattering amplitude is determined from off-shell currents which
are composed of all sub-diagrams connecting a certain sub-set of external
particles. These currents depend on the momenta $p_{1},\ldots,p_{n}$ of
external particles $1,\ldots,n$, and on their helicities and colors.

The off-shell currents satisfy the recursive relations

$J_{\alpha}(\pi)=P_{\alpha}(\pi)\,\sum\limits_{\mathcal{V}_{\alpha}^{\;\alpha_{1},\,\alpha_{2}}}\sum\limits_{\mathcal{P}_{2}(\pi)}\mathcal{S}(\pi_{1},\pi_{2})\;\mathcal{V}_{\alpha}^{\,\alpha_{1},\,\alpha_{2}}(\pi_{1},\pi_{2})\,J_{\alpha_{1}}(\pi_{1})J_{\alpha_{2}}(\pi_{2})\;,$
(9)

where, $P_{\alpha}(\pi)$ denotes a propagator, which depends on the particle
type $\alpha$ and the set of particles, $\pi$. The term
$\mathcal{V}_{\alpha}^{\,\alpha_{1},\alpha_{2}}(\pi_{1},\pi_{2})$ represents a
vertex, which depends on the particle types $\alpha$, $\alpha_{1}$ and
$\alpha_{2}$ and the decomposition of the set of particle labels, $\pi$, into
disjoint subsets $\pi_{1}$ and $\pi_{2}$. The quantity
$\mathcal{S}(\pi_{1},\pi_{2})$ is the symmetry factor associated with the
decomposition of $\pi$ into $\pi_{1}$ and $\pi_{2}$ Gleisberg and Höche
(2008). Superscripts refer to incoming particles at the vertex, subscripts to
outgoing particles. The sums run over all vertices in the interaction model
and over all unordered partitions $\mathcal{P}_{2}$ of the set $\pi$ into two
disjoint subsets. Since we will use the above recursion only to compute the
leptonic current, we can suppress color indices. Helicity labels are implied.
The initial currents $J_{\alpha_{i}}(p_{i})$ can be determined based on the
spin of the $i^{\text{th}}$ external particle. The external currents are given
by:

$J_{\alpha}(p_{i})=\begin{cases}1&\text{spin}=0\\\ u(p_{i})\ \text{or}\
v(p_{i})&\text{spin}=1/2\\\
\epsilon_{\alpha}(p_{i},k)&\text{spin}=1\end{cases},$ (10)

where $u(p_{i})$ and $v(p_{i})$ are the solutions to the Dirac equation,
$\epsilon_{\alpha}$ is the polarization vector, and $k$ is an auxiliary vector
to define the polarization.

A complete list of interaction vertices for the Standard Model can be found in
Ref. Gleisberg and Höche (2008). For completeness, we also list the nontrivial
expressions for the natively implemented Lorentz structures in App. A. Comix
allows to include Feynman rules for nearly arbitrary new physics scenarios
into Eq. (9) by means of a generic interface to the UFO file format Degrande
_et al._ (2012), and an automated generator for the Lorentz (and color)
structures Höche  _et al._ (2015); Krauss  _et al._ (2017). This extension of
Comix to BSM calculations also discusses the extension of Eq. (9) from
containing only 3-point vertices to arbitrary $N$-point vertices. The
generation of UFO files can be accomplished through the use of the FeynRules
program Christensen and Duhr (2009); Alloul  _et al._ (2014), which can take
any Lagrangian and obtain the needed Feynman rules for tree-level
calculations.

In order to compute the leptonic current needed to evaluate Eq. (1), we
introduce point-like nucleons into the theory, which act as auxiliary
particles that are needed only for bookkeeping purposes and to construct a
formally complete scattering amplitude. The off-shell currents coupling to
these nucleons are then used to define the leptonic current or the leptonic
tensor as

$\begin{split}L_{\mu}^{(i)}(1,\ldots,m)=&\;J^{(i)}_{\mu}(1,\ldots,m)\;,\\\
L_{\mu\nu}^{(ij)}(1,\ldots,m)=&\;J^{(i)}_{\mu}(1,\ldots,m)J^{(j)\dagger}_{\nu}(1,\ldots,m)\;,\end{split}$
(11)

where the particles $1,\ldots,m$ are the non strongly interacting particles
that contribute to the leptonic current.

Through the use of the recursion relations, and the interface with UFO files
generated from FeynRules, this generator can calculate almost any leptonic
current that may be of interest. The only limitations on the leptonic currents
that can be calculated are: can not currently handle any spin $>$ 1 particles,
can not have any color charged particles, only spin-1 particles can probe the
nucleus, and only tree-level diagrams can be calculated. Of these limitations,
only the exclusion of the color charged particles can not be resolved with
future work. Allowing color charged particles to interact with the nucleus
breaks the assumption that the degrees of freedom in QCD are protons and
neutrons. Implementing other probes of the nucleus involves updating the
nuclear physics to include the appropriate form factors for the different spin
probes. Handling particles with spin $>$ 1 requires implementing the needed
external particle states and appropriate propagators. The automation of one-
loop diagrams has been discussed in detail in Refs. Hirschi  _et al._ (2011);
Denner  _et al._ (2018); Buccioni  _et al._ (2019).

##  IV Phase-space integration

We employ the recursive phase-space generation techniques of Ref. Byckling and
Kajantie (1969) in combination with the multi-channel method of Ref. Kleiss
and Pittau (1994) and the adaptive multi-dimensional integration algorithm
Vegas Lepage (1980) in order to perform the phase-space integrals. Considering
a $2\to n$ scattering process with incoming particles $a$ and $b$ and outgoing
particles $1\ldots n$, the $n$-particle differential phase space element reads

${\rm d}\Phi_{n}(a,b;1,\ldots,n)=\left[\,\prod\limits_{i=1}^{n}\frac{{\rm
d}^{4}p_{i}}{(2\pi)^{3}}\,\delta(p_{i}^{2}-m_{i}^{2})\Theta(p_{i0})\,\right]\,(2\pi)^{4}\delta^{(4)}\left(p_{a}+p_{b}-\sum_{i=1}^{n}p_{i}\right)\;,$
(12)

where $p_{i}$ and $m_{i}$ are the momentum and on-shell masses of outgoing
particles, respectively.

For concreteness, consider the $2\to 2$ process $l+{}^{12}C\to l^{\prime}+N+X$
in the quasielastic regime. In this case, we not only have to consider the
final state lepton and nucleon, but also the initial state of the system. In
the quasielastic regime, the initial momentum of the nucleon can be generated
by considering an isotropic three-momentum and a removal energy of the nucleon
inside the nucleus. Furthermore, if the lepton is not monochromatic, then the
momentum of the lepton can be generated according to the initial flux. Putting
all the components together, the full differential phase space element is
given as ${\rm d}\Phi_{2}(a,b;1,2){\rm d}^{3}\mathbf{p}_{a}{\rm d}E_{r}{\rm
d}^{3}\mathbf{p}_{b}$, where ${\rm d}\Phi_{2}(a,b;1,2)$ is given in Eq. (12),
$\mathbf{p}_{a}(E_{r})$ is the three-momentum (removal energy) of the initial
nucleon, and $\mathbf{p}_{b}$ is the three-momentum of the initial lepton.
This process will proceed through a $t$-channel exchange of a gauge boson.
Therefore, it is efficient to rewrite the two-body final state phase space
such that it appropriately samples a $t$-channel momentum distribution as:

${\rm
d}\Phi_{2}(a,b;1,2)=\frac{\lambda\left(s_{ab},s_{1},s_{2}\right)}{16\pi^{2}2s_{ab}}\,{\rm
d}\cos\theta_{1}{\rm d}\phi_{1},$ (13)

where $\theta_{1}(\phi_{1})$ is the polar (azimuthal) angle with respect to
the axis formed by $p_{a}+p_{b}$. Here we have introduced the Källen function:

$\lambda^{2}\left(s_{a},s_{b},s_{c}\right)={\left(s_{a}-s_{b}-s_{c}\right)^{2}-4s_{b}s_{c}},$
(14)

where $s_{i}$ denotes the invariant mass of the particles $i$. The
differential phase-space element ${\rm d}^{3}\mathbf{p}_{a}$ can be rewritten
as:

${\rm d}^{3}\mathbf{p}_{a}=|\mathbf{p}_{a}|^{2}{\rm d}\mathbf{p}_{a}{\rm
d}\cos\theta_{a}{\rm d}\phi_{a},$ (15)

where $E_{r}$ is the removal energy, $\theta_{a}(\phi_{a})$ is the polar
(azimuthal) angle with respect to the beam direction, and the energy of the
nucleon is given as $E_{a}=m_{N}-E_{r}$ to convert the quasielastic energy
conservation $\delta$-function given in Eq. (7) to
$\delta\left(E_{a}+E_{b}-E_{1}-E_{2}\right)$, which is consistent with the
energy conservation of Eq. (12). The initial distribution for the lepton
depends on the experiment being considered, and is typically given as a
probability distribution given some initial momentum (and position in the case
of fluxes from NuMI beam simulations Adamson  _et al._ (2016)). Events can be
generated according to this probability distribution with various MC
techniques. In this work, we consider only monochromatic beams, so details of
these methods are left to a future work.

To generalize the phase-space integration to higher multiplicities, we follow
Ref. James , factorizing Eq. (12) as

${\rm d}\Phi_{n}(a,b;1,\ldots,n)={\rm
d}\Phi_{n-m+1}(a,b;\pi,m+1,\ldots,n)\,\frac{{\rm d}s_{\pi}}{2\pi}\,{\rm
d}\Phi_{m}(\pi;1,\ldots,m)\;.$ (16)

In this context, $\pi=\\{1,\ldots,m\\}$ corresponds to a subset of particle
indices, similar to the notation in Sec. III. We use overlined letters to
denote the missing subset, e.g.
$\overline{\alpha}=\\{a,b,1,\ldots,n\\}\setminus\alpha$. Equation (16) allows
the decomposition of the complete phase space into building blocks
corresponding to the propagator-like term ${\rm d}s_{\pi}/(2\pi)$ and the
$s$\- and $t$-channel decay processes

$\begin{split}S_{\pi}^{\,\rho,\pi\setminus\rho}&=\frac{\lambda(s_{\pi},s_{\rho},s_{\pi\setminus\rho})}{16\pi^{2}\,2\,s_{\pi}}\;{\rm
d}\cos\theta_{\rho}\,{\rm d}\phi_{\rho}\;,\\\
T_{\alpha,b}^{\,\pi,\overline{\alpha b\pi}}&=\frac{\lambda(s_{\alpha
b},s_{\pi},s_{\;\overline{\alpha b\pi}})}{16\pi^{2}\,2s_{\alpha b}}\;{\rm
d}\cos\theta_{\pi}\,{\rm d}\phi_{\pi}\;,\\\ \end{split}$ (17)

where $T_{\alpha,b}^{\,\pi,\overline{\alpha b\pi}}$ is the same function as in
Eq. (13). In addition, an overall factor of $(2\pi)^{4}\,{\rm
d}^{4}p_{\,\overline{\alpha
b}}\;\delta^{(4)}(p_{\alpha}+p_{b}-p_{\,\overline{\alpha b}})$ guarantees
four-momentum conservation. The basic idea of Ref. Byckling and Kajantie
(1969) is to match the indices in the virtuality integrals of Eq. (16) and the
lower left indices in the $T$ functions to indices of $t$-channel currents in
the hard matrix element. In this manner, one obtains an optimal integrator for
a particular Feynman diagram in a scalar theory that contains propagators with
all these indices. An optimal integrator for all scalar diagrams can then be
constructed using the multi-channel technique Kleiss and Pittau (1994).
Finally, the spin dependence of the hard matrix element can be mapped out more
carefully with the help of adaptive MC algorithms Lepage (1978); Bothmann  _et
al._ (2020); Gao  _et al._ (2020).

##  V Nuclear Matrix Element Interface

While the UFO file is extremely flexible, it is missing a way to easily define
the nuclear form factors involved within the interaction. We propose to use an
extension to the UFO file format that provides a way to consistently interface
with the form factors used in the neutrino event generators and other
neutrino-nucleus interaction codes. This extension will only be valid as long
as the Conserved Vector Current (CVC) hypothesis is valid. The CVC provides a
method to relate the vector form factor for an arbitrary model to the
electromagnetic form factors ($F_{1},F_{2}$). The electomagnetic form factors
$F_{1}$ and $F_{2}$ can be defined in terms of the electric and magnetic form
factors as:

$\displaystyle F_{1}^{p,n}$ $\displaystyle=\frac{G_{E}^{p,n}+\tau
G_{M}^{p,n}}{1+\tau}$ $\displaystyle F_{2}^{p,n}$
$\displaystyle=\frac{G_{E}^{p,n}-G_{M}^{p,n}}{1+\tau},$ (18)

with $\tau=-q^{2}/(2M)^{2}$. In this work, we consider the Kelly
parameterization for the electric and magnetic form factors Kelly (2004).

Additionally, we will consider a global axial-vector form factor as well.
While it is straightforward to implement any form factor, such as the
“z-expansion” Meyer  _et al._ (2016), we choose to use the dipole
parameterization to allow for a fair comparison to Ref. Rocco  _et al._
(2019a). The standard dipole parameterization is given as

$F_{A}=\frac{g_{A}}{(1-q^{2}/M_{A}^{2})^{2}},$ (19)

where the nucleon axial-vector coupling constant is taken to be $g_{A}=1.2694$
and the axial mass is taken to be $M_{A}$ = 1.0 GeV. Studying the impact of
different parameterizations of the axial-vector form factor are beyond the
scope of this work.

The form factors for an arbitrary process can be expressed in the form factors
given above through the use of their isospin dependent interactions. Here we
work out the form factors for the Standard Model photon, $Z$ boson, and $W$
boson, but the process can be applied similarly to any BSM scenario. Since the
form factors are defined in terms of the coupling to the photon, the photon
form factors ($\mathcal{F}^{\gamma}_{i}$; see Eq. (II)) are straightforward,
and given as:

$\mathcal{F}^{\gamma(p,n)}_{i}=F^{p,n}_{i},\quad\quad\mathcal{F}^{\gamma}_{A}=0.$
(20)

The $W^{\pm}$ boson couples via isospin through the $\tau^{\pm}$ operator,
which leads to the form factors:

$\mathcal{F}^{W(p,n)}_{i}=F^{p}_{i}-F^{n}_{i},\quad\quad\mathcal{F}^{W}_{A}=F_{A}.$
(21)

The $Z$ boson also couples via isospin, and the form factors are

$\displaystyle\mathcal{F}_{i}^{Z(p)}$
$\displaystyle=\left(\frac{1}{2}-2\sin^{2}\theta_{W}\right)F_{i}^{p}-\frac{1}{2}F_{i}^{n},\quad\quad\mathcal{F}_{A}^{Z(p)}=\frac{1}{2}F_{A},$
$\displaystyle\mathcal{F}_{i}^{Z(n)}$
$\displaystyle=\left(\frac{1}{2}-2\sin^{2}\theta_{W}\right)F_{i}^{n}-\frac{1}{2}F_{i}^{p},\quad\quad\mathcal{F}_{A}^{Z(n)}=-\frac{1}{2}F_{A},$
(22)

where $\sin\theta_{W}$ is the weak mixing angle. We leave the examples of
working out the form factors for specific BSM scenarios to a future work.

##  VI Results

We consider the processes of electron and neutrino scattering off Carbon, and
neutrino trident production off Carbon. The electron scattering and neutrino
scattering are used as a means to validate the code. The trident production is
the first fully-differential calculation of tridents in the quasielastic
region using nuclear spectral functions, to demonstrate the ability to extend
the predictions beyond $2\to 2$ processes. For simplicity, we only consider a
monochromatic beam for the incoming particle, but it is straightforward to
include appropriate fluxes. Additionally, we do not include any final state
interaction effects, which are known to shift the location of the peak of the
distribution, reduce the height of the peak and broaden the tails of the
distribution. For all the calculations, we set our electroweak parameters in
the $(\alpha,G_{F},M_{Z})$-scheme, and the values are chosen as
$\alpha=1/137$, $G_{F}=1.16637\times 10^{-5}$ GeV-2, and $M_{Z}=91.1876$ GeV.
Throughout the calculations, all leptons are considered massless particles.
Furthermore, we include the full propagator and do not take the infinite mass
limit for the contributions from $W$ or $Z$ bosons. This differs from previous
works in which the four-point Fermi interactions are used to model the
neutrino-nucleus interactions.

###  VI.1 Electron Scattering

There is a plethora of electron scattering data on a Carbon target that has
been collected. Here we compare to a select few experimental results, along
with the predictions based on Ref. Rocco  _et al._ (2019b). The data is the
double differential cross section in solid angle and the energy transfer
($\omega$) versus the energy transfer to the nucleus.

First, we consider the electron scattering off Carbon with an electron beam
energy of 961 MeV and 37.5 degrees scattering angle Sealock  _et al._ (1989).
The comparison between the data, the results based on Ref. Rocco  _et al._
(2019b), and our work can be seen in Fig. 1. We observe that the calculation
is consistent with the original theory calculation. Furthermore, the
difference in the quasielastic peak region between the theory calculation and
the data can be explained by including final state interaction effects, as
shown in Ref. Rocco  _et al._ (2019b). The high energy transfer tail can be
explained through the inclusion of two-body currents, resonance production,
and deep-inelastic scattering Rocco  _et al._ (2019b).

Figure 1: A comparison between our calculation, the data from Ref. Sealock
_et al._ (1989), and the theory calculation based on Ref. Rocco  _et al._
(2019b) for a 961 MeV electron scattering off Carbon at an angle of 37.5
degrees. The two theory calculations are consistent with each other, and the
difference with the data can be explained by including final state interaction
effects. Figure 2: A comparison between our calculation, the data from Ref.
Sealock  _et al._ (1989), and the theory calculation based on Ref. Rocco  _et
al._ (2019b) for a 1300 MeV electron scattering off Carbon at an angle of 37.5
degrees. The two theory calculations are consistent with each other, and the
difference with the data can be explained by including final state interaction
effects.

Figure 2 displays the comparison between the two theory approaches and data
for electron scattering off Carbon with an electron beam energy of 1300 MeV
and 37.5 degrees scattering angle Sealock  _et al._ (1989). Again, we observe
that the two theory calculations are consistent. There is an overall agreement
between the theory prediction based on Ref. Rocco  _et al._ (2019b) and the
data in the quasielastic region, the small difference at low $\omega$ has to
be ascribed to final state interaction effects as in Fig. 2. These can be
accounted for as demonstrated in Ref. Rocco  _et al._ (2019b) but their
inclusion is beyond the scope of this work.

###  VI.2 Neutrino Scattering

For neutrino scattering, we consider both neutrino and anti-neutrino beams.
Since there are no high-energy monochromatic beams of neutrinos, we only
compare our results to those from Ref. Rocco  _et al._ (2019a), but do not
compare to any experimental data. For both neutrino and anti-neutrino beams,
we consider both total cross section as a function of incoming neutrino
energy, along with the double differential cross section in outgoing angle and
energy of the outgoing lepton versus the energy transfer to the nucleus. We
consider only charged-current (CC) interactions here. In order to compare to
the results of Ref. Rocco  _et al._ (2019a), we include the effects from Pauli
blocking by ensuring that the outgoing nucleon has a momentum greater than the
Fermi momentum $k_{F}=225$ MeV.

As shown in Fig. 3, the calculation using the current approach agrees (up to
statistical uncertainties) with the calculation based on Ref. Rocco  _et al._
(2019a).

Figure 3:  A comparison between our calculation and the SF IA based on Ref.
Rocco  _et al._ (2019a) for the total cross section for charged current
scattering off a Carbon nucleus versus neutrino energy. On the left is the
neutrino channel, while on the right is the anti-neutrino channel.

Additionally, the differential cross section is compared in Fig. 4 for the CC
interaction at a fixed outgoing lepton angle of $30^{\circ}$ (left) and
$70^{\circ}$ (right) for an incoming 1 GeV neutrino and anti-neutrino. Here we
again see consistency between our method and that of Ref. Rocco  _et al._
(2019a).

Figure 4: A comparison between our calculation and the SF IA based on Ref.
Rocco  _et al._ (2019a) for the double differential cross section for a 1 GeV
neutrino on a Carbon target, with a fixed outgoing lepton angle of 30 degrees
on the left and 70 degrees on the right.

###  VI.3 Neutrino Tridents

$\nu_{\mu}$$\nu_{\mu}$$e^{+}$$e^{-}$$N$$N$$Z$$Z/\gamma$$\nu_{\mu}$$\nu_{\mu}$$e^{-}$$e^{+}$$N$$N$$Z$$Z/\gamma$

Figure 5: The two Feynman diagrams representing the neutrino trident cross
section.

To demonstrate that this generator is not limited to $2\to 2$ process only and
that it can handle interference terms, we consider the neutrino trident
process $\nu_{\mu}\,^{12}{\rm C}\to\nu_{\mu}e^{+}e^{-}X$ with a fixed neutrino
beam energy of 1 GeV, including both the $Z$ and photon interactions with the
nucleon. Additionally, this process is an important background to understand
for multiple lepton final state BSM explanations of the MiniBooNE excess. The
Feynman diagrams for this process can be found in Fig. 5. In order to regulate
the electron propagator pole, we set a minimum opening angle between the two
electrons of $5^{\circ}$. In addition to this cut, we also require that the
electrons have a minimum energy of 30 MeV, and have an angle from the neutrino
beam axis greater than $10^{\circ}$. The electroweak parameters and form
factors are identical to those used in the electron and neutrino scattering
results section. With this setup, we obtain a total cross section of
$3.973\times 10^{-14}\pm 2.764\times 10^{-17}$ nb, which is consistent with
the results of Ref. Ballett  _et al._ (2019c). In addition, we show the
results for the electron pair opening angle, the leading electron energy, the
sub-leading electron energy, and the invariant mass of the electron pair. The
electron pair opening angle can be seen in Fig. 6, and is defined as the angle
between the two electrons. This observable is important in understanding the
ability of the next-generation experiments to observe this process for the
first time based on their resolution for separating the two electrons from a
single electron. The leading and sub-leading electron energy can be seen in
Fig. 7. In comparing the leading and sub-leading energies, we see that one
electron tends to be significantly softer than the other electron. Finally,
the electron pair invariant mass can be seen in Fig. 8. This is an important
observable to distinguish neutrino trident processes from BSM scenarios that
have a pair of electrons in the final state, such as Ref. Bertuzzo  _et al._
(2018b). In the case of Ref. Bertuzzo  _et al._ (2018b), the $Z^{\prime}$ will
create a bump in the invariant mass spectrum, which would be easily
distinguishable from the standard neutrino trident prediction. We leave a
detailed analysis of separating the standard neutrino trident from
$Z^{\prime}$ models to a future work.

Figure 6: The opening angle between the two electrons in the neutrino trident
process with an incoming neutrino beam of 1 GeV. A cut is placed on the
minimum opening angle of $5^{\circ}$, and the electrons are required to have a
minimum energy of 30 MeV.

Figure 7: The leading (left) and sub-leading (right) energy for the outgoing
electrons in the neutrino trident process with an incoming neutrino beam of 1
GeV. A cut is placed on the minimum opening angle between the electrons of
$5^{\circ}$, and the electrons are required to have a minimum energy of 30
MeV. Figure 8: The invariant mass of the electron pair from the neutrino
trident process with an incoming neutrino beam of 1 GeV. A cut is placed on
the minimum opening angle between the electrons of $5^{\circ}$, and the
electrons are required to have a minimum energy of 30 MeV.

##  VII Conclusions

We developed a novel event generation framework for the automated simulation
of neutrino scattering at next-generation neutrino experiments. The framework
takes inspiration from similar tools developed for the automatic simulation of
events for the LHC community. The major complication that does not exist at
the LHC is the handling of the nuclear physics effects, which we address by
interfacing to a dedicated code for nuclear physics models in the quasielastic
scattering region. Adding in the two-body currents, resonance production,
shallow inelastic scattering, and deep inelastic scattering contributions is a
logical extension of this work. This would require updating the handling of
the initial states for the nuclear side, modifying the energy conserving delta
function such that the phase space techniques described here can be applied,
and implementing the additional nuclear effects. We leave this to a future
work.

We demonstrate that we reproduce the expected results for the commonly studied
processes of electron scattering and neutrino scattering off nuclei.
Furthermore, we demonstrate the ability of the framework to compute processes
beyond $2\to 2$ scattering by studying neutrino tridents. This process is
important for multi-lepton final state explanations of the MiniBooNE excess.
We show a variety of differential distributions demonstrating that this
framework is capable of simulating full-differential events for subsequent
analysis in the experimental simulation pipelines.

With the development of our new event generator, it becomes straightforward to
study possible beyond Standard Model physics scenarios in a rigorous manner.
Since the results obtained from the generator are fully differential in the
many-body phase space and include the complete nuclear effects, our framework
can assist the experiments in defining improved search strategies to separate
various BSM scenarios from the Standard Model and from each other. We leave
the details of this procedure to a future work.

The source code for our event generator will be provided upon request, and
made public upon publication of this work.

##  VIII Acknowledgments

We thank Pedro Machado, William Jay, Alessandro Lovato, and Gil Paz for useful
discussions and advice throughout the development of this tool. This
manuscript has been authored by Fermi Research Alliance, LLC under Contract
No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science,
Office of High Energy Physics.

##  Appendix A Natively supported Lorentz structures

To exemplify the generality of our code, we list in this appendix the
nontrivial expressions for the natively implemented Lorentz structures, not
including those that correspond to simple contractions of external
polarization vectors or spinors. The building blocks given here can be
extended to nearly arbitrary interactions (also including higher-point
functions) by means of the interface to FeynRules and UFO published in Ref.
Höche  _et al._ (2015).

$\displaystyle\parbox{56.9055pt}{\begin{center} \fmfset{thick}{1.25thin}
\fmfset{arrow_len}{2mm}\fmfset{curly_len}{1.5mm}
\fmfset{wiggly_len}{1.5mm}\fmfset{decor_size}{2mm}
\fmfset{dot_len}{1mm}\fmfset{dot_size}{2thick} \fmfset{dash_len}{1.5mm}
\fmfgraph*(25,20) \fmftop{i,j} \fmfbottom{ij} \fmf{fermion}{v1,i}
\fmf{photon}{j,v1} \fmf{fermion,tension=1.5}{v2,v1}
\fmf{phantom,tension=5}{ij,v2}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle\bar{u}$}{i}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle j$}{j}
\fmfv{d.shape=circle,d.size=3,l.a=0}{v2}
\end{center}}=\bar{u}j^{\mu}\gamma_{\mu}$
$\displaystyle=\left(0,0,\bar{u}_{0}j^{-}-\bar{u}_{1}j_{\perp},-\bar{u}_{0}j_{\perp}^{*}+\bar{u}_{1}j^{+}\right)\;,$
(23) $\displaystyle\parbox{56.9055pt}{\begin{center} \fmfset{thick}{1.25thin}
\fmfset{arrow_len}{2mm}\fmfset{curly_len}{1.5mm}
\fmfset{wiggly_len}{1.5mm}\fmfset{decor_size}{2mm}
\fmfset{dot_len}{1mm}\fmfset{dot_size}{2thick} \fmfset{dash_len}{1.5mm}
\fmfgraph*(25,20) \fmftop{i,j} \fmfbottom{ij} \fmf{photon}{v1,i}
\fmf{fermion}{j,v1} \fmf{fermion,tension=1.5}{v1,v2}
\fmf{phantom,tension=5}{ij,v2} \fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle
j$}{i} \fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle v$}{j}
\fmfv{d.shape=circle,d.size=3,l.a=0}{v2}
\end{center}}=j^{\mu}\gamma_{\mu}\frac{1-\gamma^{5}}{2}v$
$\displaystyle=\left(\begin{array}[]{c}0\\\ 0\\\
j^{+}v_{0}+j_{\perp}^{*}v_{1}\\\
j_{\perp}v_{0}+j^{-}v_{1}\end{array}\right)\;,$ (28)
$\displaystyle\parbox{56.9055pt}{\begin{center} \fmfset{thick}{1.25thin}
\fmfset{arrow_len}{2mm}\fmfset{curly_len}{1.5mm}
\fmfset{wiggly_len}{1.5mm}\fmfset{decor_size}{2mm}
\fmfset{dot_len}{1mm}\fmfset{dot_size}{2thick} \fmfset{dash_len}{1.5mm}
\fmfgraph*(25,20) \fmftop{i,j} \fmfbottom{ij} \fmf{fermion}{i,v1,j}
\fmf{photon,tension=1.5}{v1,v2} \fmf{phantom,tension=5}{ij,v2}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle\bar{u}$}{i}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle v$}{j}
\fmfv{d.shape=circle,d.size=3,l.a=0}{v2}
\end{center}}=\bar{u}\gamma^{\mu}\frac{1-\gamma^{5}}{2}v=\left(\begin{array}[]{c}\bar{u}_{0}v_{2}+\bar{u}_{1}v_{3}\\\
\bar{u}_{0}v_{3}+\bar{u}_{1}v_{2}\\\
i\left(\bar{u}_{1}v_{2}-\bar{u}_{0}v_{3}\right)\\\
\bar{u}_{0}v_{2}-\bar{u}_{1}v_{3}\end{array}\right)\;,$ (33)
$\displaystyle\parbox{56.9055pt}{\begin{center} \fmfset{thick}{1.25thin}
\fmfset{arrow_len}{2mm}\fmfset{curly_len}{1.5mm}
\fmfset{wiggly_len}{1.5mm}\fmfset{decor_size}{2mm}
\fmfset{dot_len}{1mm}\fmfset{dot_size}{2thick} \fmfset{dash_len}{1.5mm}
\fmfgraph*(25,20) \fmftop{i,j} \fmfbottom{ij} \fmf{fermion}{v1,i}
\fmf{photon}{j,v1} \fmf{fermion,tension=1.5}{v2,v1}
\fmf{phantom,tension=5}{ij,v2}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle\bar{u}$}{i}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle j$}{j}
\fmfv{d.shape=circle,d.size=3,l.a=0}{v2}
\end{center}}=\bar{u}j^{\mu}\gamma_{\mu}\frac{1+\gamma^{5}}{2}$
$\displaystyle=\left(\bar{u}_{2}j^{+}+\bar{u}_{3}j_{\perp},\bar{u}_{2}j_{\perp}^{*}+\bar{u}_{3}j^{-},0,0\right)\;,$
(34) $\displaystyle\parbox{56.9055pt}{\begin{center} \fmfset{thick}{1.25thin}
\fmfset{arrow_len}{2mm}\fmfset{curly_len}{1.5mm}
\fmfset{wiggly_len}{1.5mm}\fmfset{decor_size}{2mm}
\fmfset{dot_len}{1mm}\fmfset{dot_size}{2thick} \fmfset{dash_len}{1.5mm}
\fmfgraph*(25,20) \fmftop{i,j} \fmfbottom{ij} \fmf{photon}{v1,i}
\fmf{fermion}{j,v1} \fmf{fermion,tension=1.5}{v1,v2}
\fmf{phantom,tension=5}{ij,v2} \fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle
j$}{i} \fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle v$}{j}
\fmfv{d.shape=circle,d.size=3,l.a=0}{v2}
\end{center}}=j^{\mu}\gamma_{\mu}\frac{1+\gamma^{5}}{2}v$
$\displaystyle=\left(\begin{array}[]{c}j^{-}v_{2}-j_{\perp}^{*}v_{3}\\\
-j_{\perp}v_{2}+j^{+}v_{3}\\\ 0\\\ 0\end{array}\right)\;,$ (39)
$\displaystyle\parbox{56.9055pt}{\begin{center} \fmfset{thick}{1.25thin}
\fmfset{arrow_len}{2mm}\fmfset{curly_len}{1.5mm}
\fmfset{wiggly_len}{1.5mm}\fmfset{decor_size}{2mm}
\fmfset{dot_len}{1mm}\fmfset{dot_size}{2thick} \fmfset{dash_len}{1.5mm}
\fmfgraph*(25,20) \fmftop{i,j} \fmfbottom{ij} \fmf{fermion}{i,v1,j}
\fmf{photon,tension=1.5}{v1,v2} \fmf{phantom,tension=5}{ij,v2}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle\bar{u}$}{i}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle v$}{j}
\fmfv{d.shape=circle,d.size=3,l.a=0}{v2}
\end{center}}=\bar{u}\gamma^{\mu}\frac{1+\gamma^{5}}{2}v=\left(\begin{array}[]{c}\bar{u}_{2}v_{0}+\bar{u}_{3}v_{1}\\\
-\bar{u}_{2}v_{1}-\bar{u}_{3}v_{0}\\\
i\left(\bar{u}_{2}v_{1}-\bar{u}_{3}v_{0}\right)\\\
-\bar{u}_{2}v_{0}+\bar{u}_{3}v_{1}\end{array}\right)\;,$ (44)
$\displaystyle\parbox{56.9055pt}{\begin{center} \fmfset{thick}{1.25thin}
\fmfset{arrow_len}{2mm}\fmfset{curly_len}{1.5mm}
\fmfset{wiggly_len}{1.5mm}\fmfset{decor_size}{2mm}
\fmfset{dot_len}{1mm}\fmfset{dot_size}{2thick} \fmfset{dash_len}{1.5mm}
\fmfgraph*(25,20) \fmftop{i,j} \fmfbottom{ij} \fmf{photon}{i,v1,j}
\fmf{photon,tension=1.5}{v1,v2} \fmf{phantom,tension=5}{ij,v2}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle j(p)$}{i}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle j^{\prime}(q)$}{j}
\fmfv{d.shape=circle,d.size=3,l.a=0}{v2}
\end{center}}=\Gamma^{\mu\lambda\kappa}\left(p,q\right)j_{\lambda}j^{\prime}_{\kappa}$
$\displaystyle=jj^{\prime}\left(p-q\right)^{\mu}+j^{\prime}\left(2q+p\right)j^{\mu}-j\left(2p+q\right)j^{\prime\mu}\;,$
(45) $\displaystyle\parbox{56.9055pt}{\begin{center} \fmfset{thick}{1.25thin}
\fmfset{arrow_len}{2mm}\fmfset{curly_len}{1.5mm}
\fmfset{wiggly_len}{1.5mm}\fmfset{decor_size}{2mm}
\fmfset{dot_len}{1mm}\fmfset{dot_size}{2thick} \fmfset{dash_len}{1.5mm}
\fmfgraph*(25,20) \fmftop{i,j,k} \fmfbottom{ijk} \fmf{photon}{v1,i}
\fmf{photon}{v1,j} \fmf{photon}{v1,k} \fmf{photon,tension=3}{v1,v2}
\fmf{phantom,tension=5}{ijk,v2} \fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle
j$}{i} \fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle j^{\prime}$}{j}
\fmfv{d.size=0,l.d=3.0,l.a=90,l=$\displaystyle j^{\prime\prime}$}{k}
\fmfv{d.shape=circle,d.size=3,l.a=0}{v2}
\end{center}}=(2g^{\mu\nu}g^{\lambda\kappa}-g^{\mu\lambda}g^{\nu\kappa}-g^{\mu\kappa}g^{\nu\lambda})j_{\nu}j^{\prime}_{\kappa}j^{\prime\prime}_{\lambda}$
$\displaystyle=2(j^{\prime}j^{\prime\prime})j^{\mu}-(jj^{\prime})j^{\prime\prime\mu}-(jj^{\prime\prime})j^{\prime\mu}\;,$
(46)

## References

  * Abi  _et al._ (2021) B. Abi _et al._ (DUNE), Eur. Phys. J. C 81, 322 (2021), arXiv:2008.12769 [hep-ex] .
  * Abi  _et al._ (2020) B. Abi _et al._ (DUNE), (2020), arXiv:2002.03005 [hep-ex] .
  * Abe  _et al._ (2018) K. Abe _et al._ (Hyper-Kamiokande), (2018), arXiv:1805.04163 [physics.ins-det] .
  * Yano (2021) T. Yano (Hyper-Kamiokande), PoS ICRC2021, 1193 (2021).
  * Aguilar  _et al._ (2001) A. Aguilar _et al._ (LSND), Phys. Rev. D 64, 112007 (2001).
  * Aguilar-Arevalo  _et al._ (2007) A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. Lett. 98, 231801 (2007).
  * Aguilar-Arevalo  _et al._ (2009) A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. Lett. 102, 101802 (2009).
  * Aguilar-Arevalo  _et al._ (2010) A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. Lett. 105, 181801 (2010).
  * Aguilar-Arevalo  _et al._ (2013) A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. Lett. 110, 161801 (2013).
  * Aguilar-Arevalo  _et al._ (2018) A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. Lett. 121, 221801 (2018), arXiv:1805.12028 [hep-ex] .
  * Aguilar-Arevalo  _et al._ (2021) A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. D 103, 052002 (2021), arXiv:2006.16883 [hep-ex] .
  * Gninenko (2009) S. Gninenko, Phys. Rev. Lett. 103, 241802 (2009), arXiv:0902.3802 [hep-ph] .
  * McKeen and Pospelov (2010) D. McKeen and M. Pospelov, Phys. Rev. D 82, 113018 (2010), arXiv:1011.3046 [hep-ph] .
  * Gninenko (2011) S. N. Gninenko, Phys. Rev. D 83, 015015 (2011), arXiv:1009.5536 [hep-ph] .
  * Dib  _et al._ (2011) C. Dib, J. C. Helo, S. Kovalenko, and I. Schmidt, Phys. Rev. D 84, 071301 (2011), arXiv:1105.4664 [hep-ph] .
  * Gninenko (2012) S. Gninenko, Phys. Lett. B 710, 86 (2012), arXiv:1201.5194 [hep-ph] .
  * Masip  _et al._ (2013) M. Masip, P. Masjuan, and D. Meloni, JHEP 01, 106 (2013), arXiv:1210.1519 [hep-ph] .
  * Ballett  _et al._ (2017) P. Ballett, S. Pascoli, and M. Ross-Lonergan, JHEP 04, 102 (2017), arXiv:1610.08512 [hep-ph] .
  * Magill  _et al._ (2018) G. Magill, R. Plestid, M. Pospelov, and Y.-D. Tsai, Phys. Rev. D 98, 115015 (2018), arXiv:1803.03262 [hep-ph] .
  * Fischer  _et al._ (2020) O. Fischer, A. Hernández-Cabezudo, and T. Schwetz, Phys. Rev. D 101, 075045 (2020), arXiv:1909.09561 [hep-ph] .
  * Bertuzzo  _et al._ (2018a) E. Bertuzzo, S. Jana, P. A. N. Machado, and R. Zukanovich Funchal, (2018a), arXiv:1808.02500 [hep-ph] .
  * Bertuzzo  _et al._ (2018b) E. Bertuzzo, S. Jana, P. A. N. Machado, and R. Zukanovich Funchal, Phys. Rev. Lett. 121, 241801 (2018b), arXiv:1807.09877 [hep-ph] .
  * Ballett  _et al._ (2019a) P. Ballett, S. Pascoli, and M. Ross-Lonergan, Phys. Rev. D 99, 071701 (2019a), arXiv:1808.02915 [hep-ph] .
  * Argüelles  _et al._ (2019) C. A. Argüelles, M. Hostert, and Y.-D. Tsai, Phys. Rev. Lett. 123, 261801 (2019), arXiv:1812.08768 [hep-ph] .
  * Ballett  _et al._ (2019b) P. Ballett, M. Hostert, and S. Pascoli, Phys. Rev. D99, 091701 (2019b), arXiv:1903.07590 [hep-ph] .
  * Ballett  _et al._ (2020) P. Ballett, M. Hostert, and S. Pascoli, Phys. Rev. D 101, 115025 (2020), arXiv:1903.07589 [hep-ph] .
  * Abdullahi  _et al._ (2020) A. Abdullahi, M. Hostert, and S. Pascoli, (2020), arXiv:2007.11813 [hep-ph] .
  * Abdallah  _et al._ (2020a) W. Abdallah, R. Gandhi, and S. Roy, (2020a), arXiv:2010.06159 [hep-ph] .
  * Abdallah  _et al._ (2020b) W. Abdallah, R. Gandhi, and S. Roy, JHEP 12, 188 (2020b), arXiv:2006.01948 [hep-ph] .
  * de Gouvêa  _et al._ (2020) A. de Gouvêa, O. L. G. Peres, S. Prakash, and G. V. Stenico, JHEP 07, 141 (2020), arXiv:1911.01447 [hep-ph] .
  * Dentler  _et al._ (2020) M. Dentler, I. Esteban, J. Kopp, and P. Machado, Phys. Rev. D 101, 115013 (2020), arXiv:1911.01427 [hep-ph] .
  * Chang  _et al._ (2021) C.-H. V. Chang, C.-R. Chen, S.-Y. Ho, and S.-Y. Tseng, Phys. Rev. D 104, 015030 (2021), arXiv:2102.05012 [hep-ph] .
  * (33) P. Abratenko (MicroBooNE), .
  * Christensen and Duhr (2009) N. D. Christensen and C. Duhr, Comput. Phys. Commun. 180, 1614 (2009), arXiv:0806.4194 [hep-ph] .
  * Alloul  _et al._ (2014) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, Comput.Phys.Commun. 185, 2250 (2014), arXiv:1310.1921 [hep-ph] .
  * Christensen  _et al._ (2011) N. D. Christensen, P. de Aquino, C. Degrande, C. Duhr, B. Fuks, M. Herquet, F. Maltoni, and S. Schumann, Eur. Phys. J. C71, 1541 (2011), arXiv:0906.2474 [hep-ph] .
  * Christensen  _et al._ (2012) N. D. Christensen, C. Duhr, B. Fuks, J. Reuter, and C. Speckner, Eur. Phys. J. C 72, 1990 (2012), arXiv:1010.3251 [hep-ph] .
  * Degrande  _et al._ (2012) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer, and T. Reiter, Comput.Phys.Commun. 183, 1201 (2012), arXiv:1108.2040 [hep-ph] .
  * Krauss  _et al._ (2002) F. Krauss, R. Kuhn, and G. Soff, JHEP 02, 044 (2002), hep-ph/0109036 .
  * Gleisberg and Höche (2008) T. Gleisberg and S. Höche, JHEP 12, 039 (2008), arXiv:0808.3674 [hep-ph] .
  * Bähr  _et al._ (2008) M. Bähr, S. Gieseke, M. A. Gigg, D. Grellscheid, K. Hamilton, O. Latunde-Dada, S. Plätzer, P. Richardson, M. H. Seymour, A. Sherstnev, and B. R. Webber, The European Physical Journal C 58, 639 (2008).
  * Bellm  _et al._ (2016) J. Bellm _et al._ , Eur. Phys. J. C 76, 196 (2016), arXiv:1512.01178 [hep-ph] .
  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H.-S. Shao, T. Stelzer, P. Torrielli, M. Zaro, and et al., Journal of High Energy Physics 2014 (2014), 10.1007/jhep07(2014)079.
  * Kilian  _et al._ (2007) W. Kilian, T. Ohl, and J. Reuter, Eur. Phys. J. C71, 1742 (2007), arXiv:0708.4233 [hep-ph] .
  * Sjöstrand  _et al._ (2006) T. Sjöstrand, S. Mrenna, and P. Skands, JHEP 05, 026 (2006), hep-ph/0603175 .
  * Sjöstrand  _et al._ (2015) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, Comput. Phys. Commun. 191, 159 (2015), arXiv:1410.3012 [hep-ph] .
  * Gleisberg  _et al._ (2009) T. Gleisberg, S. Höche, F. Krauss, M. Schönherr, S. Schumann, F. Siegert, and J. Winter, JHEP 02, 007 (2009), arXiv:0811.4622 [hep-ph] .
  * Bothmann  _et al._ (2019) E. Bothmann _et al._ (Sherpa), SciPost Phys. 7, 034 (2019), arXiv:1905.09127 [hep-ph] .
  * Juszczak  _et al._ (2006) C. Juszczak, J. A. Nowak, and J. T. Sobczyk, Nucl. Phys. B Proc. Suppl. 159, 211 (2006), arXiv:hep-ph/0512365 .
  * Andreopoulos (2009) C. Andreopoulos (GENIE), Acta Phys. Polon. B 40, 2461 (2009).
  * Hayato (2009) Y. Hayato, Acta Phys. Polon. B 40, 2477 (2009).
  * Buss  _et al._ (2012) O. Buss, T. Gaitanos, K. Gallmeister, H. van Hees, M. Kaskulov, O. Lalakulich, A. B. Larionov, T. Leitner, J. Weil, and U. Mosel, Phys. Rept. 512, 1 (2012), arXiv:1106.1344 [hep-ph] .
  * Benhar  _et al._ (2008) O. Benhar, D. day, and I. Sick, Rev. Mod. Phys. 80, 189 (2008), arXiv:nucl-ex/0603029 .
  * Rocco  _et al._ (2019a) N. Rocco, C. Barbieri, O. Benhar, A. De Pace, and A. Lovato, Phys. Rev. C 99, 025502 (2019a), arXiv:1810.07647 [nucl-th] .
  * Dickhoff and Charity (2019) W. H. Dickhoff and R. J. Charity, Prog. Part. Nucl. Phys. 105, 252 (2019), arXiv:1811.03111 [nucl-th] .
  * Dickhoff and Atkinson (2020) W. H. Dickhoff and M. C. Atkinson, J. Phys. Conf. Ser. 1643, 012082 (2020), arXiv:1909.13750 [nucl-th] .
  * Barbieri and Hjorth-Jensen (2009) C. Barbieri and M. Hjorth-Jensen, Phys. Rev. C 79, 064313 (2009), arXiv:0902.3942 [nucl-th] .
  * Benhar  _et al._ (1994) O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, Nucl. Phys. A 579, 493 (1994).
  * Sick  _et al._ (1994) I. Sick, S. Fantoni, A. Fabrocini, and O. Benhar, Phys. Lett. B 323, 267 (1994).
  * Rocco and Barbieri (2018) N. Rocco and C. Barbieri, Phys. Rev. C 98, 025501 (2018), arXiv:1803.00825 [nucl-th] .
  * Andreoli  _et al._ (2021) L. Andreoli, J. Carlson, A. Lovato, S. Pastore, N. Rocco, and R. B. Wiringa, (2021), arXiv:2108.10824 [nucl-th] .
  * Ankowski  _et al._ (2015) A. M. Ankowski, O. Benhar, and M. Sakuda, Phys. Rev. D 91, 033005 (2015), arXiv:1404.5687 [nucl-th] .
  * Benhar  _et al._ (1989) O. Benhar, A. Fabrocini, and S. Fantoni, Nucl. Phys. A 497, 423C (1989).
  * Gleisberg  _et al._ (2004) T. Gleisberg, S. Höche, F. Krauss, A. Schälicke, S. Schumann, and J. Winter, JHEP 02, 056 (2004), hep-ph/0311263 .
  * Duhr  _et al._ (2006) C. Duhr, S. Höche, and F. Maltoni, JHEP 08, 062 (2006), hep-ph/0607057 .
  * Kanaki and Papadopoulos (2000) A. Kanaki and C. G. Papadopoulos, Comput. Phys. Commun. 132, 306 (2000), hep-ph/0002082 .
  * Mangano  _et al._ (2003) M. L. Mangano, M. Moretti, F. Piccinini, R. Pittau, and A. D. Polosa, JHEP 07, 001 (2003), hep-ph/0206293 .
  * Höche  _et al._ (2015) S. Höche, S. Kuttimalai, S. Schumann, and F. Siegert, Eur. Phys. J. C75, 135 (2015), arXiv:1412.6478 [hep-ph] .
  * Krauss  _et al._ (2017) F. Krauss, S. Kuttimalai, and T. Plehn, Phys. Rev. D95, 035024 (2017), arXiv:1611.00767 [hep-ph] .
  * Hirschi  _et al._ (2011) V. Hirschi, R. Frederix, S. Frixione, M. V. Garzelli, F. Maltoni, and R. Pittau, JHEP 05, 044 (2011), arXiv:1103.0621 [hep-ph] .
  * Denner  _et al._ (2018) A. Denner, J.-N. Lang, and S. Uccirati, Comput. Phys. Commun. 224, 346 (2018), arXiv:1711.07388 [hep-ph] .
  * Buccioni  _et al._ (2019) F. Buccioni, J.-N. Lang, J. M. Lindert, P. Maierhöfer, S. Pozzorini, H. Zhang, and M. F. Zoller, Eur. Phys. J. C 79, 866 (2019), arXiv:1907.13071 [hep-ph] .
  * Byckling and Kajantie (1969) E. Byckling and K. Kajantie, Nucl. Phys. B9, 568 (1969).
  * Kleiss and Pittau (1994) R. Kleiss and R. Pittau, Comput. Phys. Commun. 83, 141 (1994), arXiv:hep-ph/9405257 [hep-ph] .
  * Lepage (1980) G. P. Lepage, (1980), cLNS-80/447.
  * Adamson  _et al._ (2016) P. Adamson _et al._ , Nucl. Instrum. Meth. A 806, 279 (2016), arXiv:1507.06690 [physics.acc-ph] .
  * (77) F. James, CERN-68-15.
  * Lepage (1978) G. P. Lepage, J. Comput. Phys. 27, 192 (1978).
  * Bothmann  _et al._ (2020) E. Bothmann, T. Janßen, M. Knobbe, T. Schmale, and S. Schumann, SciPost Phys. 8, 069 (2020), arXiv:2001.05478 [hep-ph] .
  * Gao  _et al._ (2020) C. Gao, S. Höche, J. Isaacson, C. Krause, and H. Schulz, Phys. Rev. D 101, 076002 (2020), arXiv:2001.10028 [hep-ph] .
  * Kelly (2004) J. J. Kelly, Phys. Rev. C 70, 068202 (2004).
  * Meyer  _et al._ (2016) A. S. Meyer, M. Betancourt, R. Gran, and R. J. Hill, Phys. Rev. D 93, 113015 (2016), arXiv:1603.03048 [hep-ph] .
  * Rocco  _et al._ (2019b) N. Rocco, S. X. Nakamura, T. S. H. Lee, and A. Lovato, Phys. Rev. C 100, 045503 (2019b), arXiv:1907.01093 [nucl-th] .
  * Sealock  _et al._ (1989) R. M. Sealock _et al._ , Phys. Rev. Lett. 62, 1350 (1989).
  * Ballett  _et al._ (2019c) P. Ballett, M. Hostert, S. Pascoli, Y. F. Perez-Gonzalez, Z. Tabrizi, and R. Zukanovich Funchal, JHEP 01, 119 (2019c), arXiv:1807.10973 [hep-ph] .


# A Comment on the “Photon position operator with commuting components” by
Margaret Hawton

A. Jadczyk<EMAIL_ADDRESS>Laboratoire de Physique Théorique,
Université de Toulouse III & CNRS
and
Ronin Institute, Montclair, NJ 0704  A. M. Schlichtinger  Faculty of Physics
and Astronomy, University of Wroclaw
pl. M. Borna 9 50-204 Wroclaw, Poland

###### Abstract

We show that the position operator with commuting components proposed by M.
Hawton [M. Hawton, Phys. Rev. A 59, 954 (1999)] and developed in subsequent
papers, including the recent ones, does not have the properties required for a
photon position operator. Depending on the exact interpretation of the results
it either concerns a triplet of massless spin zero particles rather than the
photon, or does not have the required covariance properties under space
rotations.

##  I Introduction

The problem of localizability of photons is usually discussed in term of
existence of a self-adjoint position operator $\bf{Q}$ having the required
natural commutation relations with the generators of space translations, space
rotations as well as the space and time inversions within a unitary
representation of the Poincaré group acting on the Hilbert space of the photon
states. It is well known (cf. e.g. [1] and references therein) that there is a
unique solution to this problem: the so-called Pryce operator [2] with
noncommuting components. 111The nonexistence of such a position operator with
commuting components has been proven with a complete mathematical rigor (using
only bounded operators: unitary operators and orthogonal projection operators)
in [3, Theorem 2, p. 157].

Notwithstanding these rigorous mathematical results, in several papers (cf.
[4, 5, 6]), authored by M. Hawton et al., a statement has been made that there
exists a whole family of self-adjoint position operators satisfying the
requirements set in [1, 3], yet having, nevertheless, commuting components. We
have examined carefully this evident contradiction and have found that it has
its source in insufficient mathematical and logical precision of the
interpretation of the results obtained in Ref. [4]. It is our intent to fill
this precision gap, and to prevent further propagation of the error made in
[4] and employed recently in the recent studies of the photon Berry’s phase
problem [6]and of photon quantum mechanics in real Hilbert space [7].222The
error propagated also into PhD theses, seee.g. Ref. [8].

In Sec. II.1 we introduce the Lie algebra of the Poincaré group and its
particular realization in the Hilbert space ${\mathcal{H}}$ of vector valued
complex functions ${\mathbf{f}}({\mathbf{p}})$ with Lorentz non-invariant
measure $d^{3}p.$ We define the Pauli-Lubanski pseudovector $W^{\mu}$ and the
helicity operator $\Lambda$, defined by $W^{\mu}=\Lambda P^{\mu}$, and we
split ${\mathcal{H}}$ into photon subspace ${\mathcal{H}}_{ph}$ and its
orthogonal complement ${\mathcal{H}}_{0}$ of spin zero states. In Sec. II.2 we
state the standard conditions for photon position operator and define the
Pryce operator (with noncommuting components) - the unique solution of these
conditions. In Sec. II.3 we define the polarisation basis and introduce the
position operator as proposed in Ref. [4]. In Sec. III we identify the problem
with the the false staements made in Ref. [4]. Since [4] lacks clarity
concerning the selection of the particular representation of the Poincaré
group, we analyze all thee available representations and show that the
position operator defined in [4] does not have the required properties in all
three cases (though in each case for a different reason).

###  I.1 Notation

We use the Minkowski metric
$(\eta_{\mu\nu})=(\eta^{\mu\nu})=\mbox{diag}(+1,-1,-1,-1)$. We will use this
metric to rise and lower indices of vectors and tensors. Our $3$-vectors are,
as a rule, contravariant vectors with upper indices. Thus, for example, ${\bf
p}$ is a vector with components $(p^{1},p^{2},p^{3})$ etc. The three spin
matrices ${\mathbf{S}}$ are always written with lower indices
${\mathbf{S}}=(S_{1},S_{2},S_{3}).$ To simplify the notation we use a system
of units in which $c=1$ and $\hbar=1.$

##  II Photon wave function

Photons are defined through a particular unitary representation of the
Poincaré group characterized by mass zero and helicity $\pm 1.$ It is,
however, convenient to include states with helicity zero as well.

###  II.1 The Lie algebra

In any unitary representation the ten self-adjoint generators of the Poincaré
group satisfy the following commutation relations:

$\displaystyle[P_{\mu},P_{\nu}]=0,\,[P_{\mu},M_{\rho\sigma}]=i(\eta_{\mu\rho}P_{\sigma}-\eta_{\mu\sigma}P_{\rho}),$
$\displaystyle[M_{\mu\nu},M_{\rho\sigma}]=i(\eta_{\nu\rho}M_{\mu\sigma}-\eta_{\mu\rho}M_{\nu\sigma}+\eta_{\mu\sigma}M_{\nu\rho}-\eta_{\nu\sigma}M_{\mu\rho}).$

The rotation and the boost generators ${\bf M}$ and ${\bf N}$ are then defined
as $M^{i}=\frac{1}{2}\epsilon^{ijk}M_{jk},\,N^{i}=M^{0i}=M_{i0}.$ Thus, for
instance, $M^{1}=M^{23}=M_{23}$, etc.

Let ${\cal H}$ be the Hilbert space of square integrable $3$–component complex
functions on $\mathbb{R}^{3}$ with a scalar product333 In Ref. [4] a family of
scalar products with $d^{3}p/p_{0}$ replaced by $d^{3}p/p_{0}^{2\alpha}$ is
being considered. They all lead to unitarily equivalent representations of the
Poincaré group. The standard Lorentz invariant scalar product is obtained by
setting $\alpha=\frac{1}{2}.$ Here we choose $\alpha=0.$

$({\bf f},{\bf g})=\sum_{i=1}^{3}\int
f_{i}^{\star}({\mathbf{p}})g_{i}({\mathbf{p}})d^{3}p.$ (1)

With the notation ${\mathbf{p}}=(p^{1},p^{2},p^{3}),$
$p^{0}=p_{0}=\left(\sum_{i=1}^{3}p^{i}\right)^{\frac{1}{2}},$
$\boldsymbol{\pi}={\mathbf{p}}/p^{0},\ \ {\bf S}=(S_{1},S_{2},S_{3}),$

$S_{1}=\left(\begin{smallmatrix}0&0&0\\\ 0&0&-i\\\
0&i&0\end{smallmatrix}\right),\,S_{2}=\left(\begin{smallmatrix}0&0&i\\\
0&0&0\\\
-i&0&0\end{smallmatrix}\right),\,S_{3}=\left(\begin{smallmatrix}0&-i&0\\\
i&0&0\\\ 0&0&0\end{smallmatrix}\right)$

the generators $P^{\mu},$ $M^{i}=\frac{1}{2}\epsilon_{ijk}M^{jk},and$
$N^{i}=M^{0i}$ of the Poincaré group are given by

$P^{\mu}=p^{\mu},\,\bf{M}=\bf{L}+\bf{S},\,{\bf N}=\bf{K}+\bf{n},$ (2)

where

${\bf L}=-i({\mathbf{p}}\times\partial/\partial{\mathbf{p}}),\,{\bf
K}=i(p_{0}\,\partial/\partial{\mathbf{p}}+\frac{1}{2}\boldsymbol{\pi}),$ (3)
$\bf{n}=\boldsymbol{\pi}\times\bf{S}.$ (4)

The above operators satisfy the following commutation relations:

$[N^{i},p^{j}]=i\delta^{ij}p^{0},\,[N^{i},p^{0}]=ip^{i},\,[M^{i},p^{j}]=i\epsilon^{ijk}p^{k},$
(5) $[M^{i},p^{0}]=0,\,[N^{i},N^{j}]=-i\epsilon^{ijk}M^{k},$ (6)
$[M^{i},N^{j}]=i\epsilon^{ijk}N^{k},\,[M^{i},M^{j}]=i\epsilon^{ijk}M^{k}.$ (7)

Replacing $\bf{M}$ by $\bf{L}$ and $\bf{N}$ by $\bf{K}$ we obtain the same
commutation relations.

The unitary space inversion operator $\Pi$ and antiunitary time inversion
operator $\Theta$ for this representation are given by (cf. [1, Eq. (2.3)]
$(\Pi{\bf f})({\mathbf{p}})={\bf f}(-{\mathbf{p}}),$ $(\Theta{\bf
f})({\mathbf{p}})={\bf f}^{*}(-{\mathbf{p}}).$ For a general representation of
the Poincaré group one defines the four-dimensional Pauli-Lubanski
pseudovector $W^{\mu}$ as
$W_{\mu}=\frac{1}{2}\epsilon_{\nu\rho\sigma\mu}\,P^{\nu}M^{\rho\sigma}.$ We
have $W^{0}={\bf P}\cdot{\bf M},\,{\bf W}=P^{0}\,{\bf M}-{\bf P}\times{\bf
N}.$

It follows from the very definition that $\eta_{\mu\nu}P^{\mu}W^{\nu}=0.$ For
mass zero representations $P^{\mu}$ is lightlike, and therefore $W^{\mu}$ is
proportional to $P^{\mu}$: $W^{\mu}=\Lambda P^{\mu}.$ The proportionality
operator $\Lambda$ commutes with all the generators and is called the helicity
operator (see e.g. [9, p.64]). In our case one easily finds that for
generators $P^{\mu},\bf{L},\bf{K}$ we have $W^{\mu}=0$, therefore $\Lambda=0,$
while for the generators $P^{\mu},\bf{M},\bf{N}$ we have

$\Lambda=\boldsymbol{\pi}\cdot{\bf M}=\boldsymbol{\pi}\cdot\bf{S}.$ (8)

From the explicit form $(\Lambda{\bf f})=-i\boldsymbol{\pi}\times{\bf f}$ it
is easy to see that the spectrum of $\Lambda$ is discrete and consists of
three points $\lambda=\pm 1$ and $0.$ Therefore $\Lambda^{2}$ is a projection
onto the subspace ${\mathcal{H}}_{ph}$ of ${\mathcal{H}}.$
${\mathcal{H}}_{ph}$ is a direct sum of eigenspaces ${\mathcal{H}}_{\pm}$ of
$\Lambda$ corresponding to eigenvalues $\lambda=\pm 1.$ The photon states are
represented by vectors in ${\mathcal{H}}_{ph}.$ The orthogonal complement
${\mathcal{H}}_{0}$ of ${\mathcal{H}}_{ph}$ in ${\mathcal{H}}$ describes a
spinless particle. It easily follows from these definitions that

${\mathcal{H}}_{ph}=\\{{\bf f}\in{\mathcal{H}}:\,{\mathbf{p}}\cdot{\bf
f}({\mathbf{p}})=0\\},$ and that ${\mathcal{H}}_{0}=\\{{\bf
f}\in{\mathcal{H}}:\,{\bf f}({\mathbf{p}})=c({\mathbf{p}}){\mathbf{p}}\\}$ for
some scalar function $c({\mathbf{p}}).$
It is obvious from the above that $P_{0}$ and $M_{\mu\nu}$ leave the subspaces
${\mathcal{H}}_{ph}$ and ${\mathcal{H}}_{0}$ invariant.

###  II.2 The Pryce position operator

The standard requirements for the selfadjoint photon position operator
$X^{i}=X^{i*}$ involve covariance under space translations and rotations as
well as invariance under space and time reflections, and that $X^{i}$ leave
the photon subspace ${\mathcal{H}}_{ph}$ invariant. They are expressed by the
following formulas:

$\begin{split}[M^{i},X^{j}]&=i\epsilon^{ijk}X^{k},\,[P^{i},X^{j}]=-i\delta^{ij},\\\
\Pi X^{i}\Pi&=-X^{i},\,\Theta
X^{i}\Theta=X^{i},\,[X^{i},\Lambda^{2}]=0.\end{split}$ (9)

The last condition is equivalent to: if ${\bf
f}({\mathbf{p}})=c({\mathbf{p}})\,{\mathbf{p}},$ where $c({\mathbf{p}})$ is
any scalar valued complex function of ${\mathbf{p}}$, then $X^{i}{\bf
f}({\mathbf{p}})=c^{i}({\mathbf{p}})\,{\mathbf{p}},$ for some functions
$c^{i}({\mathbf{p}}),$ which is easier to verify in practice.

It is shown in Ref. [1] that the above problem has a unique solution. It is
the so-called Pryce operator ${\bf X}_{P}$ given by

$X_{P}^{i}=i\frac{\partial}{\partial
p^{i}}+\sum_{j,k}\epsilon_{ijk}\frac{P^{j}S_{k}}{P_{0}^{2}}.$ (10)

###  II.3 Polarization basis

Let ${\mathbf{e}}_{i}({\mathbf{p}})$ the triplet of vectors tangent to the
coordinate lines $p_{1},p_{2},p_{3}$

${\mathbf{e}}_{1}=\left(\begin{smallmatrix}1\\\ 0\\\
0\end{smallmatrix}\right),\,{\mathbf{e}}_{2}=\left(\begin{smallmatrix}0\\\
1\\\
0\end{smallmatrix}\right),\,{\mathbf{e}}_{3}=\left(\begin{smallmatrix}0\\\
0\\\ 1\end{smallmatrix}\right).$ (11)

The procedure of constructing the position operator with commuting components,
as proposed in Ref. [4] starts with a choice of another
${\mathbf{p}}$-dependent (orthonormal in $\mathbb{C}^{3}$) basis
$\tilde{{\mathbf{e}}}_{i}({\mathbf{p}})$, with
$\tilde{{\mathbf{e}}}_{1}({\mathbf{p}}),\tilde{{\mathbf{e}}}_{2}({\mathbf{p}})$
almost everywhere differentiable functions of ${\mathbf{p}}$ and orthogonal,
while ${\mathbf{p}}$ and
$\tilde{{\mathbf{e}}}_{3}({\mathbf{p}})={\mathbf{p}}/|{\mathbf{p}}|$. Apart of
these requirements the basis needs not be unspecified. 444In Refs. [4, 5, 6]
the basis is a real basis, and specified as related to polar coordinates in
the momentum space - see below. But the essence of the construction of the
‘position operator with commuting components’ does not depend on any
particular choice of the basis, except of the requirements mentioned above.
555Since the (complexified) tangent bundle of the sphere $S^{2}$ is nontrivial
(see e.g. [12]), the functions $\tilde{{\mathbf{e}}}_{i}({\mathbf{p}})$ are
defined only almost everywhere. Notice that they are ‘improper’
(unnormalizable) elements of ${\mathcal{H}}.$

Let $U({\mathbf{p}})$ be the unitary operator in $\mathbb{C}^{3}$ defined by

$U({\mathbf{p}}){\mathbf{e}}_{i}({\mathbf{p}})=\tilde{{\mathbf{e}}}_{i}({\mathbf{p}}),\,U^{*}({\mathbf{p}})\tilde{{\mathbf{e}}}_{i}({\mathbf{p}})={\mathbf{e}}_{i}({\mathbf{p}})\quad(i=1,2,3)$
(12)

for almost all ${\mathbf{p}}$. Then
$U^{*}({\mathbf{p}})\tilde{{\mathbf{e}}}_{i}({\mathbf{p}})={\mathbf{e}}_{i}({\mathbf{p}})$
and the matrix-valued function $U({\mathbf{p}})$ defines a unitary operator
$U$ in ${\mathcal{H}}.$ The operator $U$ maps the subspaces generated by
${\mathbf{e}}_{i}$ onto the subspaces generated by $\tilde{{\mathbf{e}}}_{i},$
while $U^{*}$ does the converse.

Let the hermitian operators $Q^{i}$ on ${\mathcal{H}}$ be defined as

$Q^{i}=i\frac{\partial}{\partial p^{i}}.$ (13)

The operators $Q^{i}$ evidently commute and leave the subspaces generated by
${\mathbf{e}}_{i}$ invariant. But they do not leave the subspace
${\mathcal{H}}_{ph}$ invariant, what constitutes the main problem with photon
localization.

We notice at this point that the operators $K^{i}$ defined by Eq. (3) can be
written as

$K^{i}=\frac{1}{2}(Q^{i}P^{0}+P^{0}Q^{i}).$ (14)

Following Ref. [4] (cf. also [6]) we introduce the operators $\hat{Q}^{i}$
defined as

$\hat{Q}^{i}=UQ^{i}U^{*}.$ (15)

Then $\hat{Q}^{i}$ are self-adjoint, commute and leave the subspaces generated
by $\tilde{{\mathbf{e}}}_{i}$ invariant. In particular $\hat{Q}_{i}$ leave the
subspace ${\mathcal{H}}_{ph}$ invariant.666They are denoted $\mathbf{r}_{i}$
in [4] and by $\mathbf{x}_{i}$ in [6].

Hawton [4, 5, 6] specifies the polarization basis to be the one induced by
spherical coordinate system in momentum space, with

$p^{1}=|{\mathbf{p}}|\sin\theta\cos\phi,\,p^{2}=|{\mathbf{p}}|\sin\theta\sin\phi,\,p^{3}=|{\mathbf{p}}|\cos\theta.$

Let
$\boldsymbol{\pi}={\mathbf{p}}/|{\mathbf{p}}|,\hat{\boldsymbol{\theta}},\hat{\boldsymbol{\phi}}$
be the unit tangent vectors to the coordinates of the spherical coordinate
system. Following Ref. [4] we call them
$\tilde{{\mathbf{e}}}_{3},\tilde{{\mathbf{e}}}_{1},\tilde{{\mathbf{e}}}_{2}.$
Explicite:

$\displaystyle\tilde{{\mathbf{e}}}_{1}({\mathbf{p}})=\left(\begin{smallmatrix}\cos\theta\cos\phi\\\
\cos\theta\sin\phi\\\
-\sin\theta\end{smallmatrix}\right),\,\tilde{{\mathbf{e}}}_{2}({\mathbf{p}})=\left(\begin{smallmatrix}-\sin\phi\\\
\cos\phi\\\ 0\end{smallmatrix}\right),$
$\displaystyle\tilde{{\mathbf{e}}}_{3}({\mathbf{p}})={\mathbf{p}}/|{\mathbf{p}}|=\left(\begin{smallmatrix}\sin\theta\cos\phi\\\
\sin\theta\sin\phi\\\ \cos\theta\end{smallmatrix}\right).$ (16)

We then define the matrix $U$ whose columns are these vectors:

$U=\begin{pmatrix}\cos\theta\cos\phi&-\sin\phi&\sin\theta\cos\phi\\\
\cos\theta\sin\phi&\cos\phi&\sin\theta\sin\phi\\\
-\sin\theta&0&\cos\theta\end{pmatrix}.$ (17)

Then $U$ satisfies Eqs. (12). The matrix $U$ is real, therefore $U^{*}=U^{t}.$
We notice that the third column describes a spinless particle. A
straightforward calculation, using Eq. (15), leads to the explicit form of
$\hat{\mathbf{Q}}$:

$\hat{\mathbf{Q}}=i\frac{\partial}{\partial{\mathbf{p}}}+\frac{1}{p_{0}^{2}}{\mathbf{p}}\times{\mathbf{S}}-\frac{\cot\theta}{p_{0}}\tilde{{\mathbf{e}}}_{2}S_{3},$
(18)

in agreement with Eq. (1) in [4].

##  III The problem and its solution

Following the definitions and the notation in Ref. [4] we introduce a new
representation of the Poincaré group defined by self-adjoint generators
$\hat{P}^{\mu},\,\hat{L}^{i},\,\hat{K}^{i}$ given by:

$\hat{P}^{\mu}=UP^{\mu}U^{*},\,\hat{L}^{i}=UL^{i}U^{*},\,\hat{K}^{i}=UK^{i}U^{*}.$
(19)

First of all we notice that, since $U$ is a multiplication by a
${\mathbf{p}}$-dependent matrix, it commutes with $P^{\mu}$. Therefore we have
that $\hat{P}^{\mu}=P^{\mu}.$ Now, since $\bf{L}=\bf{Q}\times\bf{P}$ (cf. Eq
(3)), we have

$\hat{\bf{L}}=\hat{\bf{Q}}\times\hat{\bf{P}},$ (20)

in agreement with Eq (17) of Ref. [4]. Similarly, from Eq. (14), we obtain

$\hat{\bf{K}}=\frac{1}{2}(\hat{Q}P^{0}+P^{0}\hat{Q})$ (21)

in agreement with Eq. (18) in [4].

The unitary representation of the Poincaré group defined by the generators
$\hat{P}^{\mu},\,\hat{L}^{i},\,\hat{K}^{i}$ describes a triplet of particles
of mass zero (since $\hat{P}^{2}=P^{2}=0$) and spin zero (since the Pauli-
Lubanski pseudovector $\hat{W}_{\mu},$ calculated for this representation, is
identically zero). I does nor represent photon’s boosts and rotations.

There is also another unitary representation determined by generators

$\tilde{P}^{\mu}=\hat{P}^{\mu}=P^{\mu},\,\tilde{{\bf M}}=U{\bf
M}U^{*},\,\tilde{{\bf N}}=U{\bf N}U^{*},$ (22)

unitarily equivalent to $P^{\mu},\bf{M},\bf{N},$ and describing photons of
both helicities and a scalar spin zero particle. From Eqs. (2)-(4) we
have777The operators $\tilde{S}_{i}$ are denoted $S_{{\mathbf{p}}i}$ in [4].
$\tilde{\bf{M}}=\hat{\bf{M}}+\tilde{\bf{S}},$ where
$\tilde{\bf{S}}=U{\bf{S}}U^{*}.$ The helicity operator $\tilde{\Lambda}$ for
this representation is given by $\tilde{\Lambda}=U\Lambda
U^{*}=\boldsymbol{\pi}\cdot\tilde{\bf{S}}.$ Its eigensubspaces are spanned by
the vectors $U\tilde{{\mathbf{e}}}_{i}=U^{2}\boldsymbol{\bf{e}}_{i}.$

It is stated explicitly in Ref. [4] that the operator $\hat{\bf Q}$ defined by
Eq. (13) satisfies the conditions placed on a position operator in [1, 2, 14,
15]. Since this statement is in evident contradiction with the results of the
cited references, we will now point out precisely where is the error in Ref.
[4] and in the subsequent papers.

Since the conditions on the photon position operator mentioned above and
written explicitly in Eqs. (9) are all group-theoretical, it is necessary to
list the relevant representations of the Poncaré group that are being under
consideration. We have three such representations. While the translation
generators $P^{\mu}$ are the same in all three cases, the generators of the
homogeneous Lorentz group are different. We have the original representation
with generators $M^{i},N^{i}$, the representation with “spinless” generators
$\hat{\bf L},\hat{\bf K}$ given by Eqs. (20),(21), and the complete
transformed representation with generators $\tilde{\bf M},\tilde{\bf N}$ given
by Eq. (22).

Ref. [4] is not explicit about which representation of the Poincaré group is
to be taken, therefore, in order to verify the conditions Eq. (9), we will
have to examine all three representations, and show which of these conditions
is not satisfied by Hawton’s position operator in each of the three cases.

#### The original representation ${\bf M},{\bf N}.$

In this case the operator $\hat{\bf Q}$ is not covariant under rotations,
$[M^{i},\hat{Q}^{j}]\neq i\epsilon^{ijk}\hat{Q}^{j}.$ Indeed, from Eq. (9) we
should have, in particular, $[M^{1},X^{1}]=0.$ But the explicit calculation
(using spherical coordinates in the momentum space) of this commutator on the
vector $f(p_{0},\theta,\phi)=(a(p_{0}),0,0),$ where $a$ is an arbitrary
function of $p_{0}$ is evidently $\neq 0.$ Thus Hawton position operator does
not satisfy the condition of covariance under rotations of the original
representation.

#### The representation $\hat{\bf M},\hat{\bf N}.$

Here all the conditions in Eq. 9) are satisfied, but the Pauli-Lubanski
pseudovector calculated for this representations is identically zero, thus the
helicity operator for this representation is identically zero. The
representation describes a triplet of spinless particles, not the photon.

#### The representation $\tilde{\bf M},\tilde{\bf N}.$

In this case $\hat{\bf Q}$ is covariant under translations and rotations, but
$\hat{\bf Q}$ does not commute with the helicity operator $\tilde{\Lambda}$ of
this representation  for the same reason for which ${\bf Q}$ does not commute
with $\Lambda$ defined by Eq. (8). In fact, we have $\hat{\bf Q}=U{\bf
Q}U^{*},$ $\tilde{\Lambda}=U\Lambda U^{*},$ therefore

$[\hat{\bf Q},\tilde{\Lambda}^{2}]=U[{\bf Q},\Lambda^{2}]U^{*}.$ (23)

Since $[{\bf Q},\Lambda^{2}]\neq 0,$ we also have $[\hat{\bf
Q},\tilde{\Lambda}^{2}]\neq 0.$

## References

  * [1] J. Mourad, Optimal photon localization, Phys. Lett. A 182, 319-322 (1993)
  * [2] M.H.L. Pryce, The Mass-Centre in the Restricted Theory of Relativity and its Connection with the Quantum Theory of Elementary Particles, Proc. Roy. Soc. Lond. A 195, 62-81 (1948)
  * [3] W.O. Amrein, Localizability for Particles of Mass Zero, Helv. Phys. Acta. 42, 149-190 (1969)
  * [4] M. Hawton, Photon position operator with commuting componetns, Phys. Rev. A 59, 954 (1999)
  * [5] M. Hawton, Photon wave mechanics and position eigenvectors, Phys. Rev. A 75, 062107 (2007)
  * [6] M. Hawton, V. Debierre, Photon position eigenvectors, Wigner’s little group, and Berry’s phase, J. Math. Phys. 60, 0521104 (2019)
  * [7] M. Hawton, Photon quantum mechanics in real Hilbert space, Phys. Rev. A 104, 052211 (2021)
  * [8] V. Debierre, The Photon Wave Function in Principle and Practice  PhD thesis, https://tel.archives-ouvertes.fr/tel-01406401 (2016)
  * [9] L.H. Ryder, Quantum Field Theory, (2nd ed., Cambridge University Press 1966)
  * [10] N.N. Bogolubov, A.A. Logunov, I.T. Todorov, Introduction to Axiomatic Quantum Field Theory (W.A. Benjamin, Inc., 1975)
  * [11] V.S. Varadarajan, Geometry of Quantum Theory (2nd ed., Springer 2007)
  * [12] Raymond Y. Ciao, Yong-Shi Wu, Manifestations of Berry’s Topological Phase for the Photon, Phys. Rev. Lett. 57, 933-936 (1986)
  * [13] I. Bialynicki-Birula, Z. Bialynicka-Birula, Berry’s phase in the relativistic theory of spinning particles, Phys. Rev. D 35, 2383-2387 (1987)
  * [14] T.D. Newton, E.P. Wigner, Localized States for Elementary Systems, Rev. Mod. Phys. 21, 400406 (1949)
  * [15] A.S. Wightman, On the Localizability of Quantum Mechanical Systems, Rev. Mod. Phys. 34, 845-872 (1962)
  * [16] H. Bacry, The Notions of Localizability and Space: From Eugene Wigner to Alain Connes, Nucl. Phys. B, (Proc. Suppl) 6, 222-230 (1989)