* (76) NNPDF collaboration, _Parton distributions from high-precision collider data_ , _Eur. Phys. J. C_ 77 (2017) 663 [1706.00428].
  * (77) T.-J. Hou et al., _New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC_ , 1912.10053.
  * (78) S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin and R. S. Thorne, _Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs_ , _Eur. Phys. J. C_ 81 (2021) 341 [2012.04684].
  * (79) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _Eur. Phys. J. C_ 75 (2015) 132 [1412.7420].
  * (80) D. Y. Bardin, A. Leike, T. Riemann and M. Sachwitz, _Energy Dependent Width Effects in e+ e- Annihilation Near the Z Boson Pole_ , _Phys. Lett. B_ 206 (1988) 539.
  * (81) M. Jezabek and J. H. Kuhn, _QCD Corrections to Semileptonic Decays of Heavy Quarks_ , _Nucl. Phys. B_ 314 (1989) 1.
  * (82) M. Cacciari, G. P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063 [0802.1189].
  * (83) A. Djouadi, J. Kalinowski and M. Spira, _HDECAY: A Program for Higgs boson decays in the standard model and its supersymmetric extension_ , _Comput. Phys. Commun._ 108 (1998) 56 [hep-ph/9704448].
  * (84) A. Djouadi, J. Kalinowski, M. Muehlleitner and M. Spira, _HDECAY: Twenty ++ years after_, _Comput. Phys. Commun._ 238 (2019) 214 [1801.09506].
  * (85) R. V. Harlander and W. B. Kilgore, _Higgs boson production in bottom quark fusion at next-to-next-to leading order_ , _Phys. Rev. D_ 68 (2003) 013001 [hep-ph/0304035].
  * (86) S. Frixione, _Isolated photons in perturbative QCD_ , _Phys. Lett. B_ 429 (1998) 369 [hep-ph/9801442].
  * (87) R. Gauld, A. Gehrmann-De Ridder, E. W. N. Glover, A. Huss and I. Majer, _Associated production of a Higgs boson decaying into bottom quarks and a weak vector boson decaying leptonically at NNLO in QCD_ , _JHEP_ 10 (2019) 002 [1907.05836].
  * (88) A. Rohatgi, _Webplotdigitizer: Version 4.5_ , 2021.
  * (89) NNPDF collaboration, _Parton distributions for the LHC Run II_ , _JHEP_ 04 (2015) 040 [1410.8849].
  * (90) S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky et al., _New parton distribution functions from a global analysis of quantum chromodynamics_ , _Phys. Rev. D_ 93 (2016) 033006 [1506.07443].
  * (91) L. Harland-Lang, A. Martin, P. Motylinski and R. Thorne, _Parton distributions in the LHC era: MMHT 2014 PDFs_ , _Eur. Phys. J. C_ 75 (2015) 204 [1412.3989].
  * (92) M. Czakon, A. Mitov and R. Poncelet, _NNLO QCD corrections to leptonic observables in top-quark pair production and decay_ , _JHEP_ 05 (2021) 212 [2008.11133].
  * (93) V. S. Fadin, V. A. Khoze and A. D. Martin, _How suppressed are the radiative interference effects in heavy instable particle production?_ , _Phys. Lett. B_ 320 (1994) 141 [hep-ph/9309234].
  * (94) G. Bevilacqua, H.-Y. Bi, H. B. Hartanto, M. Kraus, M. Lupattelli and M. Worek, _$t\bar{t}b\bar{b}$ at the LHC: on the size of corrections and $b$-jet definitions_, _JHEP_ 08 (2021) 008 [2105.08404].


# Early Performance Results on 4th Gen Intel(R) Xeon (R) Scalable Processors
with DDR and Intel(R) Xeon(R) processors, codenamed Sapphire Rapids with HBM
††thanks: This work was supported by the U.S. Department of Energy through the
Los Alamos National Laboratory (LANL). LANL is operated by Triad National
Security, LLC, for the National Nuclear Security Administration of U.S.
Department of Energy (Contract No. 89233218CNA000001). LA-UR-22-31875 Approved
for public release; distribution is unlimited.

‡Galen M. Shipman, Sriram Swaminarayan

Computer, Computational, and Statistical Sciences
Los Alamos National Laboratory
Los Alamos NM, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>Gary Grider, Jim Lujan

High Performance Computing
Los Alamos National Laboratory
Los Alamos NM, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>R. Joseph Zerr
X Computational Physics
Los Alamos National Laboratory
Los Alamos NM, USA
<EMAIL_ADDRESS>
###### Abstract

The Crossroads supercomputer was designed to simulate some of the most complex
physical devices in the world. These simulations routinely require 1/2
petabyte or more of system memory running on thousands of compute nodes for
months at a time on the most powerful supercomputers. Improvements in time to
solutions for these workloads can have major impact on our mission
capabilities. In this paper we present early results of representative
application workloads on 4th Gen Intel Xeon and Intel Xeon Processors
codenamed Sapphire Rapids with HBM. These results demonstrate an extremely
promising 8.57x improvement (node to node) over our prior generation Intel
Broadwell (BDW) based HPC systems. No code modifications were required to
achieve this speedup, providing a compelling path forward toward major
reductions in time to solution and the complexity of physical systems that can
be simulated in the future.

###### Index Terms:

hardware architecture, memory technology, performance analysis.

##  I Introduction

Los Alamos National Laboratory (LANL) has a long history in advancing the
state-of-the-art in High Performance Computing for complex multi-physics
simulations. The Crossroads supercomputer scheduled for deployment in 2023
will be the latest in a series of HPC systems developed in collaboration with
Intel and HPE/Cray. Crossroads will include 4th Gen Intel Xeon Scalable
processors, formerly codenamed Sapphire Rapids, based on the Intel 7 process
node and the novel Embedded Multi-Die Interconnect Bridge (EMIB) chiplet
integration technology. This CPU promises to deliver significant improvements
in simulation performance with little if any code modifications. As many of
our codes are significantly memory bound, LANL worked with Intel to develop a
Intel Xeon Scalable Processor with High Bandwidth Memory (HBM) 2e integrated
with EMIB technology.

Initial experience on pre-production 4th Gen Intel Xeon with Double Data Rate
(DDR) memory shows up to a 5.83x performance improvement in our multi-physics
simulations compared to our current Broadwell based production systems. Pre-
production Intel Xeon Processors codenamed Sapphire Rapids with HBM shows up
to a 8.57x performance improvement. These improvements were achieved with no
code changes, simply a recompile of the code to target the new
microarchitecture.

##  II 4th Gen Intel Xeon and Intel Xeon Processors codenamed Sapphire Rapids
with HBM Architectures

4th Gen Intel Xeon is the latest generation Xeon server processor based on
Golden Cove microarchitecture (using Intel 7 process technology). A number of
significant innovations from Intel are present in the processor including
EMIB, acceleration engines, CXL, and bandwidth improvements in Ultra Path
Interconnect. Microarchitecture and instructions per cycle (IPC) improvements
include a 2x increase in decode length, 50% increase in decode width, a 2.4x
increase in branch target size, 1.8x increase in the $\mu$op cache, a larger
number of execution ports and and a 1.5x increased reorder buffer. These
microarchitecture changes improve the IPC rate and enable increased
parallelism within the CPU.

4th Gen Intel Xeon is also the first CPU to incorporate Intel’s EMIB
technology providing a more scalable path to further performance improvements
and efficiency relative to traditional monolithic designs that are limited to
a single reticle. This ability to provide tight integration beyond the single
reticle limit is an important advancement, and it has other benefits including
the potential for higher yield and the ability to integrate disparate
processor and memory technologies in a single package. 4th Gen Intel Xeon is
composed of four 400$mm^{2}$ tiles each with golden cove cores with 48KB
instruction and 32KB data caches. L2 is private and is 2048KB. Our system
configurations are either dual socket 4th Gen Intel Xeon or dual socket Intel
Xeon Processors codenamed Sapphire Rapids with HBM each with 66 physical cores
for a total of 112 physical cores per node. 4th Gen Intel Xeon configuration
has 8 channels of 16 GB DDR-5 operating at 4800 MT/s per socket for a total of
256GB of memory per 2 CPU node. Intel Xeon Processors codenamed Sapphire
Rapids with HBM configuration has 4 banks of 8 high 16 Gbit HBM2e operating at
3200 MT/s per socket for a total of 128 GB of memory per 2 CPU node.

##  III Multi-Physics Simulation Codes

LANL develops and maintains a suite of multi-physics simulation codes used to
support our national security and science mission. These codes are used daily
by hundreds of scientists and engineers to answer questions of national
importance. Significant effort in validation and verification of these codes
underwrite their capabilities and provide confidence in their use to simulate
complex physical devices. As a result of this work, scientists are able to to
explore complex issues in Inertial Confinement Fusion [1], a regime of physics
that is only achievable at a single experimental device in the world, the
National Ignition Facility (NIF). The ability to simulate multi-material
reactive burn in high pressure regimes has enabled the simulation [2] of the
overdriven detonation states in the triple point overdrive experiment
conducted by LANL.

###  III-A xRAGE

xRAGE [3] is an Eulerian radiation/hydrodynamics code under active development
at LANL. xRAGE (Radiation Adaptive Grid Eulerian) supports simulations in 1D,
2D, and 3D, multiple materials, and coupling of radiation diffusion and
hydrodynamics. Adaptive Mesh Refinement (AMR) of unit aspect ratio cell
volumes supports large dynamic ranges of state space (centimeters to
kilometers, microns to meters) with refinement and derefinement of adjacent
$2^{d}$ cells (2 in 1D, 4 in 2D, and 8 in 3D) each timestep.

Most simulations using xRAGE are DRAM bandwidth bound with relatively low
arithmetic intensity. This may surprise some readers as conventional wisdom is
that HPC applications are largely FLOP bound, but prior analysis has shown
most subroutines in xRAGE simulations are completely or significantly memory
bound with an arithmetic intensity of 0.001 to 0.2 FLOP/byte using Intel
roofline [4] analysis. Analysis using BYFL [5] ”software based performance
counters” illustrated in Table I details the types of instructions and the
percentage issued from a representative xRAGE simulation. Sparse memory
operations result in load/store operations coupled with integer and array
indexing operations and represent 56% of all instructions issued in this
simulation as detailed in Table I in blue font. Fine- and coarse-grained
branching can be seen in the relatively high (25) percent of
conditional/branching instructions issued in this simulation, as detailed in
Table I in red font. Floating point operation represent only 6% of
instructions issued.

Instruction | Count | Percentage
---|---|---
Load | 6,775,030,849 | 18%
Branching | 6,063,697,707 | 16%
Integer Add | 5,334,155,682 | 15%
Array Indexing | 4,855,537,532 | 13%
Conditional | 3,299,248,274 | 9%
Store | 2,599,966,427 | 7%
Type cast | 1,959,938,043 | 5%
Sign extension | 1,541,094,404 | 4%
Stack frame allocation | 1,221,694,311 | 3%
FP multiplication | 1,171,615,897 | 3%
FP comparison | 1,141,415,386 | 3%
INT multiplication | 991,524,374 | 3%
TABLE I: Instruction breakout by type in a representative simulation [5]. Red
denotes operations attributable to branching (25%). Blue denotes operations
attributable to sparse data structures (56%).

This and deeper analysis [6] motivated a focused effort on addressing the
memory wall in the Crossroads architecture.

###  III-B FLAG

FLAG is a Lagrange-based radiation/hydrodynamics code that supports multiple
mesh optimization strategies in 1, 2, and 3D with support for fully
unstructured (polyhedral) grids. FLAG supports arbitrary Lagrangian-Eulerian
(ALE) and adaptive mesh refinement methods to handle mesh distortions /
tearing and multiple levels of resolution.

Similar to xRAGE simulations, most simulations using FLAG are also DRAM
bandwidth bound with relatively low arithmetic intensity. Memory access
patterns in FLAG are often sparse, which the Pennant [7] proxy application
attempts to capture. Analysis of Pennant, done as part of the Spatter [8]
work, shows a significant number of scatter/gather patterns with non-unit
strides. Some evidence gathered with gem5 simulations [6] further suggests
that loaded latency [8] can also impact the performance of these simulations
significantly.

The performance characteristics of these two codes, xRAGE and FLAG, are
representative of a large number of LANL codes that make use of similar data
structures. Improving the performance of these codes through architecture
innovation has a broad impact to wide range of codes and their applications at
LANL and therefore motivates our focus here. While both codes are capable of
simulating hydrodynamics coupled with radiation transport this early study
focuses exclusively on multi-material hydrodynamics test problems.

##  IV Experiments and Results

All experiments were conducted on Linux based dual-socket nodes. 4th Gen Intel
Xeon and Intel Xeon Processors codenamed Sapphire Rapids with HBM nodes were
configured as detailed above in section II. Intel Broadwell nodes were
configured with E5-2695 V4 with 4 channels per socket of DDR4 operating at
2400 MT/s. Linux processes were bound to individual cores and memory was
affinitized to the core to eliminate QPI traffic. All benchmarks were
conducted without hyperthreading and without over-subscription of cores. For
node-to-node comparisons, benchmark problems were limited to the 128 GB
available on our Intel Broadwell and Intel Xeon Processors codenamed Sapphire
Rapids with HBM.

The xRAGE code was used to benchmark an Asteroid impact problem using a
configuration similar to those described in the literature [9]. Three problem
sizes are given, one using 3 million mesh cells and a resident set size of up
to 46GB of system memory (Asteroid_3M), another using 6 million mesh cells and
a resident set size of up to 53GB of system memory (Asteroid_6M), and a third
using 16 million mesh cells and a resident set size of up to 104GB
(Asteroid_16M). Table II details the time to solution improvements of Asteroid
simulations using xRAGE on 4th Gen Intel Xeon relative to BDW DDR and Intel
Xeon Processors codenamed Sapphire Rapids with HBM relative to 4th Gen Intel
Xeon. Overall performance of both 4th Gen Intel Xeon and Intel Xeon Processors
codenamed Sapphire Rapids with HBM is significantly higher than that achieved
by BDW DDR. The smaller benchmark problem Asteroid_3M is 5.83x faster on 4th
Gen Intel Xeon relative to BDW DDR and 7.58x faster on Intel Xeon Processors
codenamed Sapphire Rapids with HBM relative to BDW DDR. Larger problems such
as Asteroid_6M and Asteroid_16M result in speedups of 5.14x and 3.39x on 4th
Gen Intel Xeon and 8.57x and 7.25x on Intel Xeon Processors codenamed Sapphire
Rapids with HBM demonstrating the significant performance improvements
resulting in the HBM memory bandwidth. Strong scaling plots from 1/4, 1/2 and
a full dual socket Broadwell, 4th Gen Intel Xeon and Intel Xeon Processors
codenamed Sapphire Rapids with HBM are illustrated in Figure 1.

Figure 1: Strong scaling of xRAGE Asteroid test problem at 1/4, 1/2, and a full dual socket Broadwell (BDW), 4th Gen Intel Xeon (SPR+DDR), and Intel Xeon Processors codenamed Sapphire Rapids with HBM (SPR+HBM) nodes | Relative time to solution
---|---
Benchmark | BDW+DDR / | BDW+DDR / | SPR+DDR /
| SPR+DDR | SPR+HBM / | SPR+HBM
Asteroid_3M | 5.83 | 7.58 | 1.29
Asteroid_6M | 5.14 | 8.57 | 1.66
Asteroid_16M | 3.39 | 7.25 | 2.14
TABLE II: Results of three simulations run using xRAGE utilizing 36 cores on
Intel Broadwell (BDW) and 112 cores on 4th Gen Intel Xeon (SPR+DDR) and Intel
Xeon Processors codenamed Sapphire Rapids with HBM (SPR+HBM) nodes. Figure 2:
Strong scaling of the FLAG Lagrange Hydro SOD test problem at 1/4, 1/2, and a
full dual socket Broadwell (BDW), 4th Gen Intel Xeon (SPR+DDR), and Intel Xeon
Processors codenamed Sapphire Rapids with HBM (SPR+HBM) nodes Figure 3: Strong
scaling of the FLAG ALE Hydro SOD test problem at 1/4, 1/2, and a full dual
socket Broadwell (BDW), 4th Gen Intel Xeon (SPR+DDR), and Intel Xeon
Processors codenamed Sapphire Rapids with HBM (SPR+HBM) nodes

The FLAG code was used to benchmark a Sod shock tube test problem described in
[10]. Two mesh management strategies are used, one using Lagrange
hydrodynamics where the mesh ”moves” during the simulation, labelled ”LH”. The
second uses ALE in which a Lagrange step is calculated and the material is
advected (Eulerian step) to the fixed mesh labelled ”ALE” for ALE with 2
materials. Table III details the time to solution improvements of Sod shock
tube simulations using FLAG on 4th Gen Intel Xeon relative to BDW DDR and
Intel Xeon Processors codenamed Sapphire Rapids with HBM relative to 4th Gen
Intel Xeon. Performance of Lagrange Hydro with a problem size of $130^{3}$ (3D
mesh) is 3.68x higher on 4th Gen Intel Xeon relative to BDW DDR and is 7.94x
higher on Intel Xeon Processors codenamed Sapphire Rapids with HBM. Smaller
problem sizes of $102^{3}$ and $80^{3}$ using LH show somewhat smaller but
still quite respectable performance improvements on 4th Gen Intel Xeon and
Intel Xeon Processors codenamed Sapphire Rapids with HBM. The smaller problem
sizes likely have a negligible cache effect here as they are still relatively
large. Performance of ALE Hydro with a problem size of $104^{3}$ (3D mesh) is
3.49x higher on 4th Gen Intel Xeon relative to BDW DDR and is 6.07x higher on
Intel Xeon Processors codenamed Sapphire Rapids with HBM. Smaller problem
sizes of $80^{3}$ and $64^{3}$ using ALE show a bit larger performance
improvements on 4th Gen Intel Xeon (3.58x and 3.6x) and a bit lower on Intel
Xeon Processors codenamed Sapphire Rapids with HBM (5.76x and 5.28x). ALE
requires significantly more instructions relative to LH for each hydrodynamics
cycle which likely accounts for some of these differences. Further analysis
using roofline modeling and frontend / backend analysis to determine potential
contributors to this result is warranted. Strong scaling plots from 1/4, 1/2
and a full dual socket Broadwell, 4th Gen Intel Xeon and Intel Xeon Processors
codenamed Sapphire Rapids with HBM for the FLAG Lagrange Hydro SOD test
problem are illustrated in Figure 2. Strong scaling results for the FLAG ALE
SOD Hydro test problem are illustrated in Figure 3.

| Relative time to solution
---|---
Benchmark | BDW+DDR / | BDW+DDR / | SPR+DDR /
| SPR+DDR | SPR+HBM | SPR+HBM
LH_NDIM_80 | 3.26 | 6.57 | 2.02
LH_NDIM_102 | 3.59 | 6.51 | 1.81
LH_NDIM_130 | 3.68 | 7.94 | 2.16
ALE_NDIM_64 | 3.60 | 5.28 | 1.47
ALE_NDIM_80 | 3.58 | 5.76 | 1.61
ALE_NDIM_104 | 3.49 | 6.07 | 1.74
TABLE III: Results of six simulations run using FLAG on two different test
problems utilizing 36 cores on Intel Broadwell (BDW) and 112 cores on 4th Gen
Intel Xeon (SPR+DDR) and Intel Xeon Processors codenamed Sapphire Rapids with
HBM (SPR+HBM) nodes.

##  V Conclusions

Multi-physics codes designed to simulate large state and phase spaces with
many materials are often memory bandwidth rather than FLOP bound. Addressing
this bottleneck has proven quite challenging as many architectures have
continued to grow in raw FLOP performance while memory bandwidth increases
have been more modest. To address this, LANL has worked closely with Intel to
integrate high-bandwidth memory in the Intel Xeon Processors codenamed
Sapphire Rapids with HBM as part of the Crossroads Supercomputer design. Early
experience with pre-production silicon shows extremely promising results with
performance improvements up to 8.57x over current HPC systems at LANL. While
other technologies have at times required significant refactoring or complete
rewrites of major portions of a code, a simple recompile was all that was
necessary when porting to Intel Xeon Processors codenamed Sapphire Rapids with
HBM. These results demonstrate a fruitful path towards efficient computing
technologies in the future for LANL’s most challenging applications.

## Acknowledgment

The authors thank members of the HPC division, the Darwin testbed team, the
Performance Engineering Team, the Applications Readiness Team, the Eulerian
Applications Project, and the Lagrangian Applications Project at Los Alamos
National Laboratory.

## References

  * [1] Fernando Grinstein, Vincent Chiravalle, and Brian Haines.  Transition and multiphysics in 3d icf capsule implosions.  Bulletin of the American Physical Society, 2022.
  * [2] Adam C Coleman, Carl E Johnson, and Matthew M Biss.  Modeling the lanl triple point overdrive experiment in the flag hydrocode.  In AIP Conference Proceedings, volume 2272, page 030007. AIP Publishing LLC, 2020.
  * [3] Michael Gittings, Robert Weaver, Michael Clover, Thomas Betlach, Nelson Byrne, Robert Coker, Edward Dendy, Robert Hueckstaedt, Kim New, W Rob Oakes, et al.  The rage radiation-hydrodynamic code.  Computational Science & Discovery, 1(1):015005, 2008.
  * [4] Diogo Marques, Helder Duarte, Aleksandar Ilic, Leonel Sousa, Roman Belenov, Philippe Thierry, and Zakhar A Matveev.  Performance analysis with cache-aware roofline model in intel advisor.  In 2017 International Conference on High Performance Computing & Simulation (HPCS), pages 898–907. IEEE, 2017.
  * [5] Scott Pakin and Patrick McCormick.  Hardware-independent application characterization.  In 2013 IEEE International Symposium on Workload Characterization (IISWC), pages 111–112. IEEE, 2013.
  * [6] Galen Shipman, Jered Dominguez-Trujillo, Kevin Sheridan, and Sriram Swaminarayan.  Assessing the memory wall in complex codes.  In 2022 IEEE/ACM Workshop on Memory Centric High Performance Computing (MCHPC), 2022.
  * [7] Charles R Ferenbaugh.  Pennant: an unstructured mesh mini-app for advanced architecture research.  Concurrency and Computation: Practice and Experience, 27(17):4555–4572, 2015.
  * [8] Patrick Lavin, Jeffrey Young, Richard Vuduc, Jason Riedy, Aaron Vose, and Daniel Ernst.  Evaluating gather and scatter performance on cpus and gpus.  In The International Symposium on Memory Systems, pages 209–222, 2020.
  * [9] Galen R Gisler, Robert P Weaver, Charles L Mader, and Michael L Gittings.  Two-and three-dimensional asteroid impact simulations.  Computing in Science & Engineering, 6(3):46–55, 2004.
  * [10] DE Burton, TC Carney, NR Morgan, and MA Kenamond.  Compendium of cch and xale test problems.  Los Alamos National Laboratory, 2013.


# Microscopic parametrization of the near threshold oscillations
of the nucleon time-like effective electromagnetic form factors

Francesco Rosini<EMAIL_ADDRESS>Simone Pacetti
Dipartimento di Fisica e Geologia, INFN Sezione di Perugia, 06123 Perugia,
Italy  Olga Shekhovtsova  INFN Sezione di Perugia, 06123 Perugia, Italy
National Science Centre, Kharkov Institute of Physics and Technology,
Akademicheskaya, Ukraine  Egle Tomasi-Gustafsson  DPhN, IRFU, CEA, Université
Paris-Saclay, 91191 Gif-sur-Yvette Cedex, France

###### Abstract

We present an analysis of the recent near threshold BESIII data for the
nucleon time-like effective form factors. The damped oscillation emerging from
the subtraction of the dipole formula is treated in non-perturbative-QCD,
making use of the light cone distribution amplitudes expansion. Non-
perturbative effects are accounted for by considering $Q^{2}$-dependent
coefficients in such expansions, whose free parameters are determined by
fitting to the proton and neutron data. Possible implications and future
analysis have been discussed.

††preprint: APS/123-QED

##  I  Introduction

The theoretical impossibility of describing the nucleon internal structure in
terms of strongly interacting quarks and gluons, which are the fundamental
fields of quantum chromodynamics, enhances the electromagnetic form factors
(EMFFs) to the role of unique and privileged tools to unravel the dynamics
underlying the electromagnetic interaction of nucleons. They provide the most
effective description of the mechanisms that determine and rule the dynamic
and static properties of nucleons. In specific reference frames, EMFFs
represent the Fourier transforms of spatial charge and magnetic momentum
densities.

Recently, the BESIII [1] experiment measured the time-like nucleon form
factors (FFs) at center-of-mass energies between 2.0 GeV and 3.5 GeV [2, 3, 4,
5, 6]. These data present an oscillating behavior [7, 8, 9, 10, 11, 12], which
manifests itself as a periodic, exponentially dumped component over the
typical dipolar carrier, usually identified as the only contribution. The
nature of such an oscillating component is still unknown. Possible
explanations rely either on the final state interaction between the baryon and
the antibaryon, or in a phenomenon intrinsic to the baryon structure. In the
latter case, the invoked phenomenon would be encoded by the EMFFs of nucleons.

In order to investigate this eventuality we propose a parametrization for the
EMFFs defined by considering the nucleons as triplets of collinear quarks
lying at light-like distances in the light-front framework [13].

The matrix element of the “+” component of the hadronic current $J^{\mu}$,
which depends directly on the EMFFs, evaluated between the baryon and
antibaryon particle states, can then be expanded using the Lorentz invariance
of the three quark Fock state’s matrix elements.

The resulting form depends on a set of functions of the four momentum squared
fractions, called light cone distribution amplitudes (LCDAs), and a deep
knowledge of their expression can provide further information about the form
factors shape. Using the $\mathcal{L}_{\text{QCD}}$ conformal symmetry [14],
the LCDAs are expanded on a polynomial basis, the most common choice being
represented by the orthonormal Appell polynomials, defined on the triangle
$T\left(x_{1},x_{3}\right)=\\{\left(x_{1},x_{3}\right)\in\mathbb{R}:x_{1}>0,x_{3}>0,x_{1}+x_{3}<1\\}$,
where $x_{i}=k_{i}^{+}/P^{+}$ are the quark’s light front momentum fractions
along is the $\left({}^{+}\right)$ direction and so the following relation
holds: $\sum_{i=1}^{3}x_{i}=1$. The only unknown quantities now are the
expansion coefficients, which have to be determined considering the
phenomenology of the problem. The nonperturbative coefficients admit an
evolution equation in the conformal symmetry framework, and their values can
be determined theoretically by QCD sum rules. On the other hand, we are
considering a center of mass energy of the system between 2.0 GeV and 3.5 GeV,
so we are not allowed to use perturbative methods. What we propose then is to
perform a truncated Laurent expansion of the non-perturbative coefficients
over the negative powers of the four momentum squared, subsequently performing
a fit over the recent BESIII experimental data to determine these
coefficients. The final goal of this description is to find whether the
oscillations of the EMFFs can be described by the model functions.

##  II The microscopic model

One of the most effective ways to describe subnuclear processes is to work on
a light front framework, expanding the involved particle states in a free
particle state basis, commonly known as Fock states. For a baryon we have

$\displaystyle\left|{\rm
baryon}\right\rangle=\left|0\right\rangle+\left|qqq\right\rangle+\left|qqqg\right\rangle+\left|qqqq\bar{q}\right\rangle+\dots\,,$
(1)

where the three-quark state can be expanded in a Lorentz series of its matrix
element between the vacuum and the particle states. The expansion has already
been performed in Ref. [15], e.g., for the proton has the form

$\displaystyle\left\langle
0\right|\varepsilon^{ijk}u^{i}_{\alpha}\left(a_{1}z\right)u^{j}_{\beta}\left(a_{2}z\right)d^{k}_{\gamma}\left(a_{3}z\right)\left|P\right\rangle$
(2) $\displaystyle=$
$\displaystyle\frac{1}{4}\Big{[}S_{1}MC_{\alpha\beta}\left(\gamma_{5}N^{+}\right)_{\gamma}+V_{1}\left(\not{p}C\right)_{\alpha\beta}\left(\gamma_{5}N^{+}\right)_{\gamma}\Big{.}$
$\displaystyle+P_{1}M\left(\gamma_{5}C\right)_{\alpha\beta}N_{\gamma}^{+}+A_{1}\left(\not{p}\gamma_{5}C\right)_{\alpha\beta}N_{\gamma}^{+}$
$\displaystyle\Big{.}+T_{1}\left(i\sigma_{\bot
p}C\right)_{\alpha\beta}\left(\gamma^{\bot}\gamma_{5}N^{+}\right)_{\gamma}+\ldots\Big{]}\,,$

where the functions $S_{1}$, $V_{1}$, $P_{1}$, $A_{1}$ and $T_{1}$ are called
light cone distribution amplitudes, they are functions of the scalar product
$P\cdot z$, $z$ being a light like four-vector. The dots in Eq. (2) indicate
that the expansion has been written explicitly only for twist-3 LCDAs, while
the complete expansion includes 24 LCDAs.

Considering now the whole expansion, we can find some conditions for the LCDAs
imposing that the nucleon state isospin is 1/2. For example, for twist-3
LCDAs, the following equation holds

$\displaystyle 2T_{1}\left(x_{1},x_{2},x_{3}\right)$ $\displaystyle=$
$\displaystyle\left[V_{1}-A_{1}\right]\left(x_{1},x_{3},x_{2}\right)$
$\displaystyle+\left[V_{1}-A_{1}\right]\left(x_{2},x_{3},x_{1}\right)\,,$

which allows to restrict the study of twist-3 LCDAs to a single function,
which is chosen to be

$\varphi_{N}\left(\mathbf{x}\right)=V_{1}\left(\mathbf{x}\right)-A_{1}\left(\mathbf{x}\right)\,,$
(3)

where $\mathbf{x}$ is 3-vector $(x_{1},x_{2},x_{3})$. Taking advantage from
the conformal symmetry of the Lagrangian density $\mathcal{L}_{\text{QCD}}$,
the twist-3 LCDA $\varphi_{N}$ can be expanded over the orthonormalized Appell
polynomials set $\left\\{P_{n}\left(\mathbf{x}\right)\right\\}_{n}$ as follows

$\varphi_{N}\left(\mathbf{x},Q^{2}\right)=120x_{1}x_{2}x_{3}\sum_{n}B_{n}\left(Q^{2}\right)P_{n}\left(\mathbf{x}\right)\,.$
(4)

The set of non-perturbative coefficients $\left\\{B_{n}\right\\}_{n}$ is
unknown and contains all the information about the form factor for the leading
twist.

Each coefficient $B_{n}$ is linked to the $\varphi_{N}$ LCDA’s momenta. The
first coefficient $B_{0}$ is fixed being linked to the normalization of
$\varphi_{N}$, i.e.,

$B_{0}=\int_{0}^{1}dx_{1}\int_{0}^{1}dx_{2}\int_{0}^{1}dx_{3}\varphi_{N}\left(\mathbf{x},Q^{2}\right)=1\,.$
(5)

As already stated, since we are considering the non perturbative aspect of the
LCDAs, we can perform an expansion over the negative powers of the four-
momentum squared,

$B_{n}\left(Q^{2}\right)=\sum_{k=0}^{M_{n}}b_{k}^{\left(n\right)}Q^{-2k},$ (6)

where $Q^{2}=-q^{2}$, $\left\\{\\{b_{k}^{(n)}\\}_{k=0}^{M_{n}}\right\\}_{n}$
is the set of coefficients and $M_{n}$ is the maximum power of $Q^{-2}$ in the
expansion of the $n^{\rm th}$ parameter $B_{n}$.

##  III Leading order contributing diagrams

Taking the leading order into account, the minimum number of contributing
diagrams has been evaluated in Ref. [15], where fourteen diagrams have been
considered. The Sachs form factors which we are interested in are related to
the Pauli and Dirac ones by the relations

$\displaystyle G_{E}\left(Q^{2}\right)$ $\displaystyle=$ $\displaystyle
F_{1}\left(Q^{2}\right)-\tau F_{2}\left(Q^{2}\right)\,,$ $\displaystyle
G_{M}\left(Q^{2}\right)$ $\displaystyle=$ $\displaystyle
F_{1}\left(Q^{2}\right)+F_{2}\left(Q^{2}\right)\,,$

where $\tau=Q^{2}/\left(4M_{B}^{2}\right)$. We performed our fit over the
effective form factor data, which is linked to the Sachs form factors by the
relation

$\displaystyle\left|G_{\text{eff}}\right|=\sqrt{\frac{\left|G_{E}\right|^{2}+2\tau\left|G_{M}\right|^{2}}{1+2\tau}}.$

Following the works of Brodsky and Lepage [16], the light front EMFF can be
written as the convolution of three probabilities, namely the probability of
describing the baryon and antibaryon as a system of three collinear quarks,
$\varphi_{N}\left(\mathbf{x},Q^{2}\right)$, and the probability of finding a
certain strong interaction, known as hard scattering kernel $K_{H}$. At the
leading order, each of the fourteen contributing diagrams corresponds to a
hard scattering kernel $K_{i}\left(\mathbf{x},\mathbf{y}\right)$, here
$\varphi_{N}\left(\mathbf{x}\right)=V_{1}\left(\mathbf{x}\right)-A_{1}\left(\mathbf{x}\right)$,
and $T\left(\mathbf{x}\right)=T_{1}\left(\mathbf{x}\right)$ are the light cone
distribution amplitudes involved in the calculation.

Table 1: Diagrams and leading twist baryon wave functions [15]. Index $i$ | Diagram | $K_{i}(\mathbf{x},\mathbf{y})$
---|---|---
1 |  | $\displaystyle\frac{\varphi_{N}(\mathbf{x})\varphi_{N}(\mathbf{y})+4T(\mathbf{x})T(\mathbf{y})}{(1-x_{1})^{2}x_{3}(1-y_{1})^{2}y_{3}}$
2 |  | 0
3 |  | $\displaystyle\frac{-4T(\mathbf{x})T(\mathbf{y})}{x_{1}x_{3}(1-x_{2})y_{1}y_{3}(1-y_{1})}$
4 |  | $\displaystyle\frac{\varphi_{N}(\mathbf{x})\varphi_{N}(\mathbf{y})}{x_{1}x_{3}(1-x_{3})y_{1}y_{3}(1-y_{1})}$
5 |  | $\displaystyle\frac{-\varphi_{N}(\mathbf{x})\varphi_{N}(\mathbf{y})}{x_{2}x_{3}(1-x_{3})y_{2}y_{3}(1-y_{1})}$
6 |  | 0
7 |  | $\displaystyle\left(\frac{1}{x_{1}y_{1}}+\frac{1}{x_{2}y_{2}}\right)\frac{\varphi_{N}(\mathbf{x})\varphi_{N}(\mathbf{y})}{(1-x_{3})^{2}(1-y_{3})^{2}}$
8 |  | 0
9 |  | $\displaystyle\frac{\varphi_{N}(\mathbf{x})\varphi_{N}(\mathbf{y})+4T(\mathbf{x})T(\mathbf{y})}{(1-x_{1})^{2}x_{2}(1-y_{1})^{2}y_{2}}$
10 |  | $\displaystyle\frac{\varphi_{N}(\mathbf{x})\varphi_{N}(\mathbf{y})+4T(\mathbf{x})T(\mathbf{y})}{(1-x_{1})^{2}x_{2}(1-y_{1})^{2}y_{2}}$
11 |  | 0
12 |  | $\displaystyle\frac{-\varphi_{N}(\mathbf{x})\varphi_{N}(\mathbf{y})}{x_{1}x_{2}(1-x_{3})y_{1}y_{2}(1-y_{1})}$
13 |  | $\displaystyle\frac{4T(\mathbf{x})T(\mathbf{y})}{x_{1}x_{2}(1-x_{1})y_{1}y_{2}(1-y_{2})}$
14 |  | $\displaystyle\frac{-\varphi_{N}(\mathbf{x})\varphi_{N}(\mathbf{y})}{x_{1}x_{2}(1-x_{1})y_{1}y_{2}(1-y_{3})}$

In order to evaluate the form factor $G_{M}$, we used the Chernyak-Zhitnitsky
[15] asymptotic formula

$\displaystyle\displaystyle
q^{4}G_{M}\left(q^{2}\right)\rightarrow\frac{\left(4\pi\bar{\alpha}_{s}\right)^{2}}{54}\left|f_{N}\right|^{2}$
$\displaystyle\displaystyle\int\\!\\!\left[dx\right]\\!\\!\int\\!\\!\left[dy\right]\left(2\sum_{i=1}^{7}e_{i}K_{i}\left(\mathbf{x},\mathbf{y}\right)+\\!\sum_{i=8}^{14}e_{i}K_{i}\left(\mathbf{x},\mathbf{y}\right)\right)\,,$
(7)

where
$\left[dx\right]=\delta\left(1-\sum_{i=1}^{3}x_{i}\right)dx_{1}dx_{2}dx_{3}$
and $\bar{\alpha}_{s}$ is the modified coupling constant. For the evaluation
of the modified coupling constant $\bar{\alpha}_{s}^{2}$, we follow the
procedure proposed by Chernyak and Zhitnitsky in Ref. [15]. The value of
$\bar{\alpha}_{s}^{2}$ is given by the product of the coupling constants for
the two subprocesses, namely the two gluon exchanges which appear in the tree-
level diagrams shown in Table 1. The average virtuality ${\bar{q}}_{1}^{2}$ of
the lightest gluon is ${\bar{x}}_{3}{\bar{y}}_{3}q^{2}$, while the rest gluon
has an averaged virtuality
${\bar{q}}_{2}^{2}=\overline{(1-x_{1})}\,\,\overline{(1-y_{1})}q^{2}$. The
typical values of a realistic nucleon wave function for the ${\bar{x}}_{i}$
are ${\bar{x}}_{1}\simeq 2/3$, ${\bar{x}}_{2}\simeq{\bar{x}}_{3}\simeq 1/6$.
Therefore
$\bar{\alpha}_{s}^{2}\left(q^{2}\right)=\alpha_{s}\left(\frac{1}{36}q^{2}\right)\alpha_{s}\left(\frac{1}{9}q^{2}\right)$.

The integrals in Eq. (III) are (weakly) convergent, it is possible to solve
them analytically. The results have already been obtained in Ref. [14], and
are reported in Appendix C. Non-perturbative-QCD effects are accounted for by
considering the $Q^{2}$-dependence of the $B_{n}$ defined in Eq. (6) as
truncated expansions in powers of $Q^{-2}$. The first parameter $B_{0}$ is
fixed to 1 and is related to the zero order of the LCDA $\varphi_{N}$, so it
is considered a constant. As for the other parameters, we are for now limiting
our discussion to the LCDA’s second order momenta, so we are only interested
in the first six parameters, namely those of the set $\\{B_{n}\\}_{n=0}^{5}$.

For the truncated expansion, we propose $M_{0}=0$, $M_{1}=M_{2}=1$,
$M_{3}=M_{4}=M_{5}=2$, so that

$\displaystyle B_{0}$ $\displaystyle=$ $\displaystyle 1\,,$ $\displaystyle
B_{n}(Q^{2})$ $\displaystyle=$ $\displaystyle
b_{0}^{(n)}+\frac{b_{1}^{(n)}}{Q^{2}}\,,\ \ \ n=1,2\,,$ $\displaystyle
B_{m}(Q^{2})$ $\displaystyle=$ $\displaystyle
b_{0}^{(m)}+\frac{b_{1}^{(m)}}{Q^{2}}+\frac{b_{2}^{(m)}}{(Q^{2})^{2}}\,,\ \ \
m=3,4,5\,.$ (8)

We obtain a closed expression of the form factor $G_{M}$ depending only on the
non-perturbative parameters for the proton and the neutron. Since the nucleons
are related by the isospin symmetry, their parameters are the same we can
performed a simultaneous fit to the recent BESIII data on proton and neutron
cross sections to determine the sets of coefficients
$\left\\{\\{b_{k}^{(n)}\\}_{k=0}^{M_{n}}\right\\}_{n=0}^{5}$. The fit is
performed using the ROOT Data Analysis Framework [17] program.

##  IV Results and discussion

Figure 1 shows the results for the effective proton (upper panel) and neutron
(lower panel) form factors in comparison with the experimental values measured
by the BESIII experiment [2, 3, 4]. The fit functions depend on 13 free
parameters, which are the coefficients of the expressions of Eq. (8). The
normalized minimum $\chi^{2}$ is

$\displaystyle\frac{\chi^{2}}{n_{\rm DoF}}=\frac{79.93}{43}\simeq 1.79\,.$

It has been obtained by using 48 data points of the proton cross section and
18 of the neutron one. The error bands have been determined by considering
both the errors of data, and the theoretical systematic error of the model,
which has ben estimated by using expressions for the $B_{n}(Q^{2})$ parameters
with the additional power $\left(Q^{-2}\right)^{M_{n}+1}$.

$G_{\rm eff}^{p}(q^{2})$

$\sqrt{q^{2}}$(GeV)

$G_{\rm eff}^{n}(q^{2})$

Figure 1: The bands represent the fit results for the proton (upper panel) and
neutron (lower panel) effective form factor. The data are from the BESIII
experiment [2, 3, 4].

$x_{1}$$x_{3}$

$\varphi_{N}(\mathbf{x})$

Figure 2: The nucleon distribution function $\varphi_{N}(\mathbf{x})$ obtained
at $\sqrt{-Q^{2}}=2.5$ GeV as a function of $x_{1}$ and $x_{3}$.

Figure 2 shows the twist-3 nucleon distribution amplitude
$\varphi_{N}(\mathbf{x})$ evaluated at $\sqrt{-Q^{2}}=2.5$ GeV. The maximum
value is reached at the light-cone momentum fractions

$x_{1}\simeq 0.49\,,\qquad x_{2}\simeq 0.24\,,\qquad x_{3}\simeq 0.27\,,$ (9)

which agree with the assumption made in the Chernyak-Zhitnitsky formula [15]
that the first quark has a momentum fraction about 50% larger than the other
two, which equally divide the remaining momentum fraction.

Summarizing, a coherent model has been developed to reproduce the data on
proton and neutron EMFFs, recently obtained by the BESIII collaboration. The
model is based on a parameterization of the light-cone distribution
amplitudes, and obeys conformal symmetry of the QCD Lagrangian.

The light-front quantization allows us to express the three-quark operator,
occurring in the expression of the hard scattering kernel, as the Fourier
transform of the light-cone distribution amplitudes, which, in turn, describe
the behaviour of the three valence quarks constituting the baryon in a light-
cone system. Always under the aegis of the QCD-Lagrangian conformal symmetry
we expanded the leading-twist baryon distribution amplitudes over a set of
Appell polynomials, which diagonalize the one gluon exchange kernel. The non-
perturbative nature of the baryon distribution amplitudes is implemented by
considering a $Q^{2}$-dependence of the expansion parameters
$\left\\{B_{n}\right\\}_{n}$. We have restricted our calculations to the first
six Appell polynomials.

These distribution amplitudes have been used to calculate the near threshold
behaviour of the nucleon effective form factors, extending the valence of an
asymptotic formula to the low-momentum transfer region. Such a near-threshold
extension has been obtained by considering a $Q^{2}$-dependence of the
expansion parameters $\left\\{B_{n}(Q^{2})\right\\}_{n}$, as polynomials of
zero, first and second degree of $Q^{-2}$.

Moreover, since the nucleon EMFFs are linked by the isospin symmetry, the same
set $\left\\{B_{n}(Q^{2})\right\\}_{n}$ can be used and the free coefficients
of their power series can be determined by means a simultaneous fit to the
proton and neutron data. The error bands of the near-threshold effective
nucleon form factors, have been obtained by considering the experimental
uncertainties of the data, used to determine the free parameters of the model
and the theoretical systematic uncertainties due to the particular
parametrization used.

This study aims to identify the origin or at least the dominant cause of the
oscillatory behaviour of the effective nucleon form factors. In particular, we
would like to distinguish between two hypotheses about the oscillation
phenomenon: the intrinsic dynamical origin and the final state interaction. In
the first case, the oscillations appear at the form factor level, while in the
second case, they are due to re-scattering reactions that happen when the
nucleons are already formed. A necessary condition in favour of the intrinsic
origin is that there exists a microscopic model of nucleons that can give
oscillatory form factors.

Even though for the neutron effective form factor the model reproduces quite
well the oscillatory behaviour, it seems to fail in the case of the proton.
Indeed, the obtained behaviour of the effective proton form factor, the orange
band shown in the upper panel of Fig. 1, is compatible with the so-called
regular background of Refs. [7, 8, 11]. It can be interpreted as the
contribution due to the short distance quark-level dynamics [18, 19], i. e.,
the $p\bar{p}$ final state is produced by the creation of quark-antiquark
pairs within a small volume, with a linear dimension much smaller than the
standard hadron size of about 1 fm.

Nevertheless, the model has the added value of proving that a unique
parametrization in all the kinematical ranges where data are present is
effective both for proton and neutron FFs. This is in contrast to previous
works, where a common fit could only be achieved either in a restricted
kinematical region, concluding in a change of the phase [3], or at the price
of three different models applicable in different kinematical regions [5].

##  Appendix A Contributing diagrams computation

Using the expressions [15]

$\displaystyle\varphi_{N}(\mathbf{x})$ $\displaystyle=$ $\displaystyle
120x_{1}x_{2}x_{3}\left(ax_{1}^{2}+bx_{2}^{2}+cx_{3}^{2}+dx_{3}+e\right)\,,$
$\displaystyle T(\mathbf{x})$ $\displaystyle=$ $\displaystyle
120x_{1}x_{2}x_{3}\left[\frac{a+c}{2}\left(x_{1}^{2}+x_{2}^{2}\right)+bp_{3}^{2}+\right.$
(10) $\displaystyle\left.\frac{d}{2}(1-x_{3})+e\right]\,,$

where the values of the coefficients $a$, $b$, $c$, $d$, $e$ and $f$ are given
in Ref. [15], the analytic solutions of the integrals of Eq. (III), namely the
ten non-vanishing expressions

$\displaystyle{\cal K}_{i}=\int[dx]\int[dy]K_{i}(\mathbf{x},\mathbf{y})\,,$
(11)

where the functions $K_{i}(\mathbf{x},\mathbf{y})$ are given in Table 1, with
$i\in\\{1,3,4,5,7,9,10,12,13,14\\}$, are:

$\displaystyle{\cal K}_{1}$ $\displaystyle=$
$\displaystyle\left(I_{1}^{\left(1\right)}\right)^{2}+4\left(I_{1}^{\left(2\right)}\right)^{2}\,,$
$\displaystyle I_{1}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{5}{12}(36a+6b+2c+12d+72e)\,,$ $\displaystyle
I_{1}^{\left(2\right)}$ $\displaystyle=$
$\displaystyle\frac{5}{12}(21a+20b+3c+32d+72e)\,;$ (12) $\displaystyle{\cal
K}_{3}$ $\displaystyle=$ $\displaystyle-4I_{3}^{(1)}I_{3}^{(2)}\,,$
$\displaystyle I^{(1)}_{3}$ $\displaystyle=$
$\displaystyle\frac{5}{3}(10a+2b+10c+15d+36e)\,,$ $\displaystyle I^{(2)}_{3}$
$\displaystyle=$ $\displaystyle\frac{5}{12}(15a+6b+15c+28d+72e)\,;$ (13)
$\displaystyle{\cal K}_{4}$ $\displaystyle=$ $\displaystyle
I_{4}^{(1)}I_{4}^{(2)}\,,$ $\displaystyle I_{4}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{5}{2}(a+3b+2(c+2d+6e))\,,$ $\displaystyle I^{(2)}_{4}$
$\displaystyle=$ $\displaystyle\frac{5}{6}(6a+9b+3c+8d+36e)\,;$ (14)
$\displaystyle{\cal K}_{5}$ $\displaystyle=$ $\displaystyle-
I_{5}^{(1)}I_{5}^{(2)}\,,$ $\displaystyle I_{5}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{5}{2}(3a+b+2(c+2d+6e))$ $\displaystyle I^{(2)}_{5}$
$\displaystyle=$ $\displaystyle\frac{10}{3}(9a+b+c+3d+18e)\,;$ (15)
$\displaystyle{\cal K}_{7}$ $\displaystyle=$
$\displaystyle\left(I^{(1)}_{7}\right)^{2}+\left(I^{(2)}_{7}\right)^{2}\,,$
$\displaystyle I^{(1)}_{7}$ $\displaystyle=$
$\displaystyle\frac{5}{6}(a+3(b+6c+8d+12e))\,,$ $\displaystyle I^{(2)}_{7}$
$\displaystyle=$ $\displaystyle\frac{5}{6}(3a+b+6(3c+4d+6e))\,;$ (16)
$\displaystyle{\cal K}_{9}$ $\displaystyle=$ $\displaystyle{\cal
K}_{10}=\left(I^{(1)}_{9}\right)^{2}+4\left(I^{(2)}_{9}\right)^{2}\,,$
$\displaystyle I^{(1)}_{9}$ $\displaystyle=$
$\displaystyle\frac{5}{6}(18a+b+3c+8d+36e)\,,$ $\displaystyle I^{(2)}_{9}$
$\displaystyle=$ $\displaystyle\frac{5}{12}(19a+6b+19c+28d+72e)\,;$ (17)
$\displaystyle{\cal K}_{12}$ $\displaystyle=$ $\displaystyle-
I_{12}^{(1)}I_{12}^{(2)}\,;$ $\displaystyle I_{12}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{10}{3}(a+b+9c+12d+18e)\,,$ $\displaystyle I_{12}^{(2)}$
$\displaystyle=$ $\displaystyle\frac{5}{6}(6a+3b+9c+16d+36e)\,;$ (18)
$\displaystyle{\cal K}_{13}$ $\displaystyle=$ $\displaystyle
4I_{13}^{(1)}I_{13}^{(2)}$ $\displaystyle I_{13}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{5}{12}(-3a+90b-3c+100d+360e)\,,$ $\displaystyle
I_{13}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{5}{12}(9a+18b+9c+20d+72e)\,;$ (19) $\displaystyle{\cal
K}_{14}$ $\displaystyle=$ $\displaystyle-I_{14}^{(1)}I_{14}^{(2)}\,,$
$\displaystyle I_{14}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{5}{6}(6a+3b+9c+16d+36e)\,,$ $\displaystyle I_{14}^{(2)}$
$\displaystyle=$ $\displaystyle\frac{10}{3}(a+b+9c+12d+18e)\,.$ (20)

##  Appendix B Distribution amplitudes computation

###  B.1 Quark distribution amplitudes complex conjugate computation

Omitting the color and current indices, there is a proportionality between the
hadronic current matrix element and the quark distribution amplitudes, given
by

$\left\langle P^{\prime}\right|\\!J\\!\left|P\right\rangle\propto\left\langle
P^{\prime}\right|\\!\bar{u}_{\mu}\bar{u}_{\nu}\bar{d}_{\rho}\\!\left|0\right\rangle\left\langle
0\right|\\!K_{H}^{\mu\nu\rho\alpha\beta\gamma}\\!\left|0\right\rangle\left\langle
0\right|\\!u_{\alpha}u_{\beta}d_{\gamma}\\!\left|P\right\rangle,$ (21)

so that the transition amplitude is proportional to the vacuum expectation
value of the hard scattering kernel. The Feynman diagram is depicted in Fig.
3.

$K_{H}^{\mu\nu\rho\alpha\beta\gamma}$$u_{\alpha}$$u_{\beta}$$d_{\gamma}$$\bar{u}_{\mu}$$\bar{u}_{\nu}$$\bar{d}_{\rho}$

Figure 3: Basic graph for the hard scattering kernel.

The hard scattering structure $K_{H}^{\mu\nu\rho\alpha\beta\gamma}$ is given
by a product of the three Lorentz structures which give the coupling of the
three fermion lines:

$K_{H}^{\mu\nu\rho\alpha\beta\gamma}=\Gamma_{1}^{\alpha\mu}\Gamma_{2}^{\beta\nu}\Gamma_{3}^{\gamma\rho}\,.$
(22)

First of all, we need to write the complex conjugate of the matrix element
representing the distribution amplitude $\varphi_{N}\left(\mathbf{x}\right)$.
We calculate this directly as we can take the complex conjugate of the three
parts composing the matrix element expansion:

$\displaystyle\left\langle
0\right|u_{\alpha}u_{\beta}d_{\gamma}\left|P\right\rangle^{*}$
$\displaystyle=$ $\displaystyle\left\langle
P\right|d^{\dagger}_{\gamma}u^{\dagger}_{\beta}u^{\dagger}_{\alpha}\left|0\right\rangle$
(23) $\displaystyle=$ $\displaystyle\left\langle
P\right|\bar{d}_{i}\bar{u}_{j}\bar{u}_{k}\left|0\right\rangle\gamma^{0}_{i\gamma}\gamma^{0}_{j\beta}\gamma^{0}_{k\alpha}\,.$

For the axial component we have

$\displaystyle\left\langle
P\right|\bar{d}_{i}\bar{u}_{j}\bar{u}_{k}\left|0\right\rangle_{A}$
$\displaystyle=$ $\displaystyle
A_{1}^{*}\left(\not{P}\gamma_{5}C\right)^{*}_{\alpha\beta}N\left(P\right)^{*}_{\gamma}\gamma^{0}_{\gamma
i}\gamma^{0}_{\beta j}\gamma^{0}_{\alpha k}$ (24) $\displaystyle=$
$\displaystyle
A_{1}^{*}\left(\left(\not{P}\gamma_{5}C\right)^{\dagger}\right)_{\beta\alpha}\bar{N}\left(P\right)_{i}\gamma^{0}_{\beta
j}\gamma^{0}_{\alpha k}$ $\displaystyle=$ $\displaystyle
A_{1}^{*}\left(C\gamma_{5}\not{P}\right)_{jk}\bar{N}\left(P\right)_{i}\,,$

for the vector one

$\displaystyle\left\langle
P\right|\bar{d}_{i}\bar{u}_{j}\bar{u}_{k}\left|0\right\rangle_{V}$
$\displaystyle=$ $\displaystyle
V_{1}^{*}\left(\not{P}C\right)^{*}_{\alpha\beta}\left(\gamma_{5}N\left(P\right)\right)^{*}_{\gamma}\gamma^{0}_{\gamma
i}\\!\gamma^{0}_{\beta j}\gamma^{0}_{\alpha k}$ (25) $\displaystyle=$
$\displaystyle-
V_{1}^{*}\left(\left(\not{P}C\right)^{\dagger}\right)_{\beta\alpha}\left(\bar{N}\left(P\right)\gamma_{5}\right)_{i}\gamma^{0}_{\beta
j}\gamma^{0}_{\alpha k}$ $\displaystyle=$ $\displaystyle-
V_{1}^{*}\left(C\not{P}\right)_{jk}\left(\bar{N}\left(P\right)\gamma_{5}\right)_{i}\,,$

and finally for the tensorial part

$\displaystyle\left\langle
P\right|\bar{d}_{i}\bar{u}_{j}\bar{u}_{k}\left|0\right\rangle_{T}$
$\displaystyle=$ $\displaystyle
T_{1}^{*}\left(\tensor{\sigma}^{\mu\nu}P_{\nu}C\right)^{*}_{\alpha\beta}\left(\gamma_{\mu}\gamma_{5}N\left(P\right)\right)^{*}_{\gamma}$
$\displaystyle\gamma^{0}_{\gamma i}\gamma^{0}_{\beta j}\gamma^{0}_{\alpha k}$
$\displaystyle=$ $\displaystyle-
T_{1}^{*}\left(\left(\tensor{\sigma}^{\mu\nu}P_{\nu}C\right)^{\dagger}\right)_{\beta\alpha}\left(\bar{N}\left(P\right)\gamma_{5}\gamma_{\mu}\right)_{i}$
$\displaystyle\gamma^{0}_{\beta j}\gamma^{0}_{\alpha k}$ $\displaystyle=$
$\displaystyle-
T_{1}^{*}\left(CP_{\nu}\tensor{\sigma}^{\mu\nu}\right)_{jk}\left(\bar{N}\left(P\right)\gamma_{5}\gamma_{\mu}\right)_{i}\,.$

###  B.2 Distribution amplitudes convolutions

Here we compute the convolution of the distribution amplitudes used for the
evaluation of the contributing diagrams. In each of the following calculations
we use the anti-commutative properties of the gamma matrices, in particular
the fact that the $\gamma_{5}$ matrix anti commutes with every component of
the matrix four-vector $\gamma^{\mu}$, while for the charge conjugation matrix
$C$ the following rule holds:

$C\gamma^{\mu}C=-\left(\gamma^{\mu}\right)^{T}\,.$ (27)

The axial, vector, tensor, and axial-vector components of the hadronic-current
matrix element are, respectively,

$\displaystyle\left\langle
P^{\prime}\right|\\!J\\!\left|P\right\rangle_{A}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\left\langle
P^{\prime}\right|\\!\bar{u}_{\mu}\bar{u}_{\nu}\bar{d}_{\rho}\\!\left|0\right\rangle_{A}\\!\left\langle
0\right|\\!K_{H}^{\mu\nu\rho\alpha\beta\gamma}\\!\left|0\right\rangle\\!\left\langle
0\right|\\!u_{\alpha}u_{\beta}d_{\gamma}\\!\left|P\right\rangle_{A}$ (28)
$\displaystyle=$ $\displaystyle
A_{1}^{*}\left(C\gamma_{5}\not{P}^{\prime}\right)_{\mu\nu}\left(\bar{N}P^{\prime}\right)_{\rho}\Gamma_{1}^{\mu\alpha}\Gamma_{2}^{\nu\beta}\Gamma_{3}^{\rho\gamma}A_{1}$
$\displaystyle\left(\not{P}\gamma_{5}C\right)_{\alpha\beta}N\left(P\right)_{\gamma}$
$\displaystyle=$
$\displaystyle\left|A_{1}\right|^{2}\bar{N}\left(P^{\prime}\right)\Gamma_{3}N\left(P\right)$
$\displaystyle\left[-\Gamma_{1}^{\mu\alpha}\left(\not{P}\gamma_{5}C\right)_{\alpha\beta}\left(\Gamma_{2}^{T}\right)^{\beta\nu}\left(C\gamma_{5}\not{P}^{\prime}\right)_{\nu\mu}\right]$
$\displaystyle=$
$\displaystyle-\left|A_{1}\right|^{2}\bar{N}\left(P^{\prime}\right)\Gamma_{3}N\left(P\right)$
$\displaystyle\text{Tr}\left[\Gamma_{1}\not{P}\gamma_{5}C\Gamma_{2}^{T}C\gamma_{5}\not{P}^{\prime}\right]$
$\displaystyle=$
$\displaystyle\left(-1\right)^{1+2n_{2}}\left|A_{1}\right|^{2}\bar{N}\left(P^{\prime}\right)\Gamma_{3}N\left(P\right)$
$\displaystyle\text{Tr}\left[\Gamma_{1}\not{P}\overset{\leftrightarrow}{\Gamma}_{2}\not{P}^{\prime}\right]\,;$
$\displaystyle\left\langle
P^{\prime}\right|\\!J\\!\left|P\right\rangle_{V}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\left\langle
P^{\prime}\right|\\!\bar{u}_{\mu}\bar{u}_{\nu}\bar{d}_{\rho}\\!\left|0\right\rangle_{V}\\!\left\langle
0\right|\\!K_{H}^{\mu\nu\rho\alpha\beta\gamma}\\!\left|0\right\rangle\left\langle
0\right|\\!u_{\alpha}u_{\beta}d_{\gamma}\\!\left|P\right\rangle_{V}$ (29)
$\displaystyle=$ $\displaystyle-
V_{1}^{*}\left(C\not{P}^{\prime}\right)_{\mu\nu}\left(\bar{N}\left(P^{\prime}\right)\gamma_{5}\right)_{\rho}\Gamma_{1}^{\mu\alpha}\Gamma_{2}^{\nu\beta}\Gamma_{3}^{\rho\gamma}$
$\displaystyle
V_{1}\left(\not{P}C\right)_{\alpha\beta}\left(N\left(P\right)\gamma_{5}\right)_{\gamma}$
$\displaystyle=$
$\displaystyle\left(-1\right)^{1+n_{3}}\left|V_{1}\right|^{2}\bar{N}\left(P^{\prime}\right)\Gamma_{3}N\left(P\right)$
$\displaystyle\left[\left(\Gamma_{1}\right)_{\mu\alpha}\left(\not{P}C\right)_{\alpha\beta}\left(\Gamma_{2}^{T}\right)_{\beta\nu}\left(C\not{P}^{\prime}\right)_{\nu\mu}\right]$
$\displaystyle=$
$\displaystyle\left(-1\right)^{1+n_{3}}\left|V_{1}\right|^{2}\bar{N}\left(P^{\prime}\right)\Gamma_{3}N\left(P\right)$
$\displaystyle\text{Tr}\left[\Gamma_{1}\not{P}C\Gamma_{2}^{T}C\not{P}^{\prime}\right]$
$\displaystyle=$
$\displaystyle\left(-1\right)^{1+n_{2}+n_{3}}\left|V_{1}\right|^{2}\bar{N}\left(P^{\prime}\right)\Gamma_{3}N\left(P\right)$
$\displaystyle\text{Tr}\left[\Gamma_{1}\not{P}\overset{\leftrightarrow}{\Gamma}_{2}\not{P}^{\prime}\right]$
$\displaystyle=$
$\displaystyle-\left|V_{1}\right|^{2}\bar{N}\left(P^{\prime}\right)\Gamma_{3}N\left(P\right)\text{Tr}\left[\Gamma_{1}\not{P}\overset{\leftrightarrow}{\Gamma}_{2}\not{P}^{\prime}\right]\,;$
$\displaystyle\left\langle
P^{\prime}\right|\\!J\\!\left|P\right\rangle_{T}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\left\langle
P^{\prime}\right|\\!\bar{u}_{\mu}\bar{u}_{\nu}\bar{d}_{\rho}\\!\left|0\right\rangle_{T}\\!\left\langle
0\right|\\!K_{H}^{\mu\nu\rho\alpha\beta\gamma}\\!\left|0\right\rangle\\!\left\langle
0\right|\\!u_{\alpha}u_{\beta}d_{\gamma}\\!\left|P\right\rangle_{T}$ (30)
$\displaystyle=$ $\displaystyle
T_{1}^{*}\left(CP^{\prime}_{\lambda}\tensor{\sigma}^{\delta\lambda}\right)_{\mu\nu}\left(\bar{N}\left(P^{\prime}\right)\gamma_{5}\gamma_{\delta}\right)_{\rho}\Gamma_{1}^{\mu\alpha}\Gamma_{2}^{\nu\beta}\Gamma_{3}^{\rho\gamma}$
$\displaystyle
T_{1}\left(\tensor{\sigma}^{\eta\kappa}P_{\kappa}C\right)_{\alpha\beta}\left(\gamma_{\eta}\gamma_{5}\left(N\left(P\right)\right)\right)_{\gamma}$
$\displaystyle=$
$\displaystyle\left|T_{1}\right|^{2}P_{\lambda}^{\prime}P_{\kappa}\bar{N}\left(P^{\prime}\right)\gamma_{\delta}\Gamma_{3}\gamma_{\eta}N\left(P\right)$
$\displaystyle\left[\Gamma_{1}^{\mu\alpha}\left(\tensor{\sigma}^{\eta\kappa}C\right)_{\alpha\beta}\left(\Gamma_{2}^{T}\right)_{\beta\nu}\left(C\tensor{\sigma}^{\delta\lambda}\right)_{\nu\mu}\right]$
$\displaystyle=$
$\displaystyle\left|T_{1}\right|^{2}P_{\lambda}^{\prime}P_{\kappa}\bar{N}\left(P^{\prime}\right)\gamma_{\delta}\Gamma_{3}\gamma_{\eta}N\left(P\right)$
$\displaystyle\text{Tr}\left[\Gamma_{1}\tensor{\sigma}^{\eta\kappa}C\Gamma_{2}^{T}C\tensor{\sigma}^{\delta\lambda}\right]$
$\displaystyle=$
$\displaystyle\left|T_{1}\right|^{2}P_{\lambda}^{\prime}P_{\kappa}\bar{N}\left(P^{\prime}\right)\gamma_{\delta}\Gamma_{3}\gamma_{\eta}N\left(P\right)$
$\displaystyle\text{Tr}\left[\Gamma_{1}\tensor{\sigma}^{\eta\kappa}\overset{\leftrightarrow}{\Gamma}_{2}\tensor{\sigma}^{\delta\lambda}\right]\,;$
$\displaystyle\left\langle
P^{\prime}\right|\bar{u}_{\mu}\bar{u}_{\nu}\bar{d}_{\rho}\left|0\right\rangle_{A}\left\langle
0\right|K_{H}^{\mu\nu\rho\alpha\beta\gamma}\left|0\right\rangle\left\langle
0\right|u_{\alpha}u_{\beta}d_{\gamma}\left|P\right\rangle_{V}$
$\displaystyle+\left\langle
P^{\prime}\right|\bar{u}_{\mu}\bar{u}_{\nu}\bar{d}_{\rho}\left|0\right\rangle_{V}\left\langle
0\right|K_{H}^{\mu\nu\rho\alpha\beta\gamma}\left|0\right\rangle\left\langle
0\right|u_{\alpha}u_{\beta}d_{\gamma}\left|P\right\rangle_{A}$
$\displaystyle=$
$\displaystyle\Big{[}A_{1}^{*}\left(C\gamma_{5}\not{P}^{\prime}\right)_{\mu\nu}\left(\bar{N}\left(P^{\prime}\right)\right)_{\rho}V_{1}\left(\not{P}C\right)_{\alpha\beta}\left(N\left(P\right)\gamma_{5}\right)_{\gamma}\Big{.}$
$\displaystyle\Big{.}-V_{1}^{*}\left(C\not{P}^{\prime}\right)_{\mu\nu}\\!\left(\bar{N}\left(P^{\prime}\right)\gamma_{5}\right)_{\rho}A_{1}\\!\\!\left(\not{P}\gamma_{5}C\right)_{\alpha\beta}N\left(P\right)_{\gamma}\\!\Big{]}$
$\displaystyle\Gamma_{1}^{\mu\alpha}\Gamma_{2}^{\nu\beta}\Gamma_{3}^{\rho\gamma}$
$\displaystyle=$
$\displaystyle-\left(V_{1}A_{1}^{*}+V_{1}^{*}A_{1}\right)\bar{N}\left(P^{\prime}\right)\Gamma_{3}\gamma_{5}N\left(P\right)\text{Tr}\left[\Gamma_{1}\not{P}\overset{\leftrightarrow}{\Gamma}_{2}\not{P}^{\prime}\right]\,.$

##  Appendix C Contributing diagrams in terms of the non-perturbative
parameters

Here we write the integrals used for the calculation of the nucleon form
factors $G_{M}^{p,n}$. The Chernyak-Zhitnitsky formula is

$\displaystyle
q^{4}G_{M}^{p,n}\left(q^{2}\right)\rightarrow\frac{\left(4\pi\bar{\alpha}_{s}\right)^{2}}{54}\left|f_{N}\right|^{2}I^{p,n}\,,$
(32)

where the integrals $I^{p,n}$ up to the second degree polynomials are

$\displaystyle I^{p}$ $\displaystyle=$ $\displaystyle
1400B_{0}B_{1}+\frac{2000}{9}B_{1}^{2}+1800B_{0}B_{2}+\frac{2800}{3}B_{1}B_{2}$
$\displaystyle+1200B_{2}^{2}+6600B_{0}B_{3}+\frac{22000}{9}B_{1}B_{3}+4800B_{2}B_{3}$
$\displaystyle+\frac{18800}{3}B_{3}^{2}-\frac{1000}{3}B_{0}B_{4}-\frac{2600}{27}B_{1}B_{4}-\frac{2200}{9}B_{2}B_{4}$
$\displaystyle-\frac{4600}{9}B_{3}B_{5}+\frac{2600}{243}B_{4}^{2}-1400B_{0}B_{5}-\frac{3500}{9}B_{1}B_{5}$
$\displaystyle-\frac{1100}{3}B_{2}B_{5}-\frac{5900}{3}B_{3}B_{5}+\frac{4100}{81}B_{4}B_{5}+\frac{7700}{27}B_{5}^{2}\,;$
$\displaystyle I^{n}$ $\displaystyle=$ $\displaystyle
1800B_{0}^{2}-1400B_{0}B_{1}+\frac{2000}{9}B_{1}^{2}-1800B_{0}B_{2}$ (33)
$\displaystyle-\frac{2800}{3}B_{1}B_{2}-200B_{2}^{2}-1000B_{0}B_{3}-\frac{22000}{9}B_{1}B_{3}$
$\displaystyle-2000B_{2}B_{3}-\frac{17000}{9}B_{3}^{2}+\frac{1000}{3}B_{0}B_{4}-\frac{2600}{27}B_{1}B_{4}$
$\displaystyle+\frac{2200}{9}B_{2}B_{4}+\frac{4600}{9}B_{3}B_{5}+\frac{2600}{243}B_{4}^{2}+\frac{2000}{3}B_{0}B_{5}$
$\displaystyle+\frac{3500}{9}B_{1}B_{5}+\frac{500}{9}B_{2}B_{5}+\frac{6500}{9}B_{3}B_{5}-\frac{4100}{81}B_{4}B_{5}$
$\displaystyle-\frac{7250}{81}B_{5}^{2}\,,$

where the coefficients $\left\\{B_{n}\right\\}_{n=0}^{5}$ are the ones
appearing in the expansion of the LCDA
$\varphi_{N}\left(\mathbf{x},Q^{2}\right)$, and their operative form is given
in Eq. (8).

## References

  * [1] Institute of high energy physics chinese academy of sciences, https://web.archive.org/web/20160316234203/http://www.ihep.ac.cn/english/E-Bepc/.
  * Ablikim  _et al._ [2019] M. Ablikim, M. N. Achasov, P. Adlarson, and N. Amhed (BESIII Collaboration), Study of the process ${e}^{+}{e}^{-}\rightarrow p\overline{p}$ via initial state radiation at besiii, Phys. Rev. D 99, 092002 (2019).
  * Ablikim  _et al._ [2021a] M. Ablikim, M. N. Achasov, P. Adlarson, S. Ahmed, and M. Albrecht, Oscillating features in the electromagnetic structure of the neutron, Nature Physics 17, 1200 (2021a).
  * Ablikim  _et al._ [2020] M. Ablikim, M. N. Achasov, P. Adlarson, and N. Amhed (BESIII Collaboration), Measurement of proton electromagnetic form factors in ${e}^{+}{e}^{-}\rightarrow p\overline{p}$ in the energy region 2.00–3.08 gev, Phys. Rev. Lett. 124, 042001 (2020).
  * Yang  _et al._ [2023] Q.-H. Yang, D. Guo, L.-Y. Dai, J. Haidenbauer, X.-W. Kang, and U.-G. Meißner, New insights into the oscillations of the nucleon electromagnetic form factors, Science Bulletin 68, 2729 (2023).
  * Ablikim  _et al._ [2021b] M. Ablikim, M. Achasov, P. Adlarson, S. Ahmed, and M. Albrecht, Measurement of proton electromagnetic form factors in the time-like region using initial state radiation at besiii, Physics Letters B 817, 136328 (2021b).
  * Bianconi and Tomasi-Gustafsson [2015] A. Bianconi and E. Tomasi-Gustafsson, Periodic interference structures in the timelike proton form factor, Phys. Rev. Lett. 114, 232301 (2015), arXiv:1503.02140 [nucl-th] .
  * Bianconi and Tomasi-Gustafsson [2016] A. Bianconi and E. Tomasi-Gustafsson, Phenomenological analysis of near threshold periodic modulations of the proton timelike form factor, Phys. Rev. C 93, 035201 (2016), arXiv:1510.06338 [nucl-th] .
  * Bianconi and Tomasi-Gustafsson [2017] A. Bianconi and E. Tomasi-Gustafsson, Fourth dimension of the nucleon structure: Spacetime analysis of the timelike electromagnetic proton form factors, Phys. Rev. C 95, 015204 (2017), arXiv:1611.02149 [hep-ph] .
  * Bianconi and Tomasi-Gustafsson [2018] A. Bianconi and E. Tomasi-Gustafsson, Soft rescattering in the timelike proton form factor within a spacetime scheme, Phys. Rev. C 98, 055204 (2018), arXiv:1809.05709 [nucl-th] .
  * Tomasi-Gustafsson  _et al._ [2021a] E. Tomasi-Gustafsson, A. Bianconi, and S. Pacetti, New fit of timelike proton electromagnetic form factors from $e^{+}e^{-}$ colliders, Phys. Rev. C 103, 035203 (2021a), arXiv:2012.14656 [hep-ph] .
  * Tomasi-Gustafsson  _et al._ [2021b] E. Tomasi-Gustafsson, A. Bianconi, and S. Pacetti, Dynamical Properties of Baryons, Symmetry 13, 1480 (2021b).
  * Brodsky  _et al._ [1998] S. J. Brodsky, H.-C. Pauli, and S. S. Pinsky, Quantum chromodynamics and other field theories on the light cone, Physics Reports 301, 299 (1998).
  * Braun  _et al._ [1999] V. Braun, S. Derkachov, G. Korchemsky, and A. Manashov, Baryon distribution amplitudes in QCD, Nuclear Physics B 553, 355 (1999).
  * Chernyak and Zhitnitsky [1984] V. Chernyak and I. Zhitnitsky, Nucleon wave function and nucleon form factors in qcd, Nuclear Physics B 246, 52 (1984).
  * Lepage and Brodsky [1980] G. P. Lepage and S. J. Brodsky, Exclusive processes in perturbative quantum chromodynamics, Phys. Rev. D 22, 2157 (1980).
  * Brun and Rademakers [1997] R. Brun and F. Rademakers, Root — an object oriented data analysis framework, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 389, 81 (1997), new Computing Techniques in Physics Research V.
  * Matveev  _et al._ [1973] V. A. Matveev, R. M. Muradian, and A. N. Tavkhelidze, Automodellism in the large - angle elastic scattering and structure of hadrons, Lett. Nuovo Cim. 7, 719 (1973).
  * Brodsky and Farrar [1973] S. J. Brodsky and G. R. Farrar, Scaling Laws at Large Transverse Momentum, Phys. Rev. Lett. 31, 1153 (1973).