# Generation and measurement of sub-picosecond electron bunch in photocathode
rf gun††thanks: Supported by National Natural Science Foundation of China
(11205152) and Science Foundation of Ministry of Education of China (“985
project”: 173123200402002)

LI Wei-Wei HE Zhi-Gang1) JIA Qi-Ka<EMAIL_ADDRESS>National Synchrotron
Radiation Laboratory, University of Science and Technology of China, Hefei,
230029, Anhui, China

###### Abstract

We consider a scheme to generate sub-picosecond electron bunch in the
photocathode rf gun by improving the acceleration gradient in the gun,
suitably tuning the bunch charge, the laser spot size and the acceleration
phase, and reducing the growth of transverse emittance by laser shaping. A
nondestructive technique is also reported to measure the electron bunch
length, by measuring the high-frequency spectrum of wakefield radiation which
is caused by the passage of a relativistic electron bunch through a channel
surrounded by a dielectric.

###### keywords:

photocathode rf gun, sub-picosecond electron bunch, bunch length measurement,
wakefield radiation

###### pacs:

2

9.25.Bx, 29.27.Ac, 29.27.Fh

##  1 Introduction

Sub-picosecond electron bunch is generally required in lots of applications,
such as the generation of coherent synchrotron radiation[2], coherent Smith-
Purcell radiation[3] and coherent cherenkov radiation[4] whose spectra are in
the THz gap; the production of sub-picosecond X-ray pulse by Compton
backscattering[5]. For the time-resolved MeV ultra-fast electron
diffraction[6], the length of electron bunch is also needed to be sub-
picosecond and as shorter as possible to achieve higher time resolution.

Magnetic compression with chicane and velocity bunching in a linear
accelerator are customary and effective ways used to get a short bunch.
However both of the two techniques take large space, while the table-top
experimental facility is needed in some applications. In this paper, we
analyse the impact factors of the electron bunch length, and consider a scheme
to directly generate the sub-picosecond electron bunch in the photocathode rf
gun by improving the acceleration gradient in the gun, suitably tuning the
bunch charge, the laser spot size and the acceleration phase, and reducing the
growth of transverse emittance by laser shaping.

In respect of the bunch length measurement, there are many methods, such as rf
deflecting cavity[7], rf zero phasing[8], electro-optical sampling[9],
coherent radiation[10, 11], and so on. Each diagnostic has its advantages and
disadvantages. Ideally, one would want a diagnostic that disturbs the bunch as
little as possible, can be single shot measurement, and uses instrumentation
that is inexpensive, easy to adjust, and routine to calibrate. A
nondestructive bunch length diagnostic technique was reported in
Reference[12], which can measure the rms bunch length by observing the
frequency spectrum of wakefield radiation as the bunch passes through a vacuum
channel in a hollow dielectric element. The design of a bunch length
measurement system based on this technique is reported in this paper.

##  2 Generation of sub-picosecond electron bunch

The length of electron bunch is affected by such factors as space charge
effect, beam energy and energy spread, and the coupling effect between the
transverse and longitudinal emittances.

The bunch lengthening due to the space charge effect can be estimated in a
drift space by[13]:

$\Delta\sigma_{z}=2qcL^{2}/I_{a}R\sigma_{z}\gamma^{4}$ (1)

where q is the charge of bunch, c is the speed of light, L is the drift
distance, $I_{a}=1.7~{}kA$ , R is the bunch radius, $\sigma_{z}$ is the bunch
length and $\gamma$ is the beam energy. In the photocathode rf gun, the energy
of electron beam is low, so the space charge effect plays the dominant role.
In order to decrease the bunch lengthening caused by the space charge effect,
the acceleration gradient should be as high as possible and the bunch charge
should be chosen appropriately. Furthermore, the bunch length can be
compressed in the gun by tuning the acceleration phase[14]. For our laser
pulse (the measured rms length is about 2.0 ps), we use the code ASTRA[15] to
simulate the bunch length evolution as a function of acceleration phase at
different charges and acceleration gradients, and the results are shown in
Fig. 1.

The rms bunch length vs. acceleration phase at different charges and
acceleration gradients.

The acceleration gradient of our photocathode rf gun, which is a second
generation gun machined by the Department of Engineering Physics of Tsinghua
University, is achieved at about 80 MV/m. The acceleration gradient of the
third generation gun is achieved at 120 MV/m at present[16], whose cathode
seal technique is improved by replacing the HELICOFLEX seal with a MATSUMOTO
gasket to eliminate the cathode gap as much as possible[17].

The bunch length can be compressed further in the drift space by a suitable
energy spread of the beam, if the length of drive laser pulse is
appropriate[18]. For our relatively short laser pulse, the bunch length
changing caused by energy spread is small because of the relatively short
drift space and the small energy spread at the exit of the gun.

The coupling between the transverse and longitudinal emittances is another
factor in lengthening the electron bunch. The laser shaping technique is an
effective way to restrain the growth of transverse emittance, which consists
of spatial and temporal shaping[19]. In this paper, we only consider the
spatial shaping. Although the uniform laser spot can be achieved by using a
spatial shaper, it is difficult to transport the shaped spot to the cathode.
So we plan to clip the laser spot by an aperture, as shown in Fig. 2.

Clipping shaping: sketch and shaping result (UV light)

For the 2 mm diameter gaussian spot and clipped spot (a gaussian spot with 0.8
mm rms size clipped by a 2 mm diameter aperture), the evolution of transverse
emittance and rms bunch length are shown as Fig. 3, where the acceleration
gradient is 120 MV/m, the bunch charge is 200 pC, the acceleration phase is 5
degrees, and the magnetic field of the solenoid is 2500 Gauss. The rms bunch
length at the focal point of the solenoid (around 0.9 m) is about 0.97 ps for
the clipped laser spot, while it is 1.26 ps for the gaussian laser spot.

Transverse emittance and rms bunch length vs longitudinal position for
different laser spots

Equation (1) shows that the bunch lengthening due to the space charge effect
is proportional to the diameter of laser spot. So the bunch length can be
modulated by tuning the diameter of laser spot, as shown in Fig. 4, and the
current distributions of bunches at the focal point of the solenoid are shown
in Fig. 5. However the transverse emittance will grow seriously when the
diameter of laser spot is too large. The optimal transverse emittances are
$0.7~{}mm\cdot mrad$ , $1.2~{}mm\cdot mrad$ and $1.75~{}mm\cdot mrad$
respectively. Nevertheless, the transverse emittance is not critical in all
the applications.

The rms bunch length vs longitudinal position for clipped laser spots with
different diameters

Current distribution of bunches at the focal point of the solenoid

In summary, to generate sub-picosecond electron bunch directly in the
photocathode rf gun, we need to improve the acceleration gradient of the gun
as high as possible, restrain the growth of transverse emittance through laser
shaping technique, and carefully tune the energy (bunch charge) and diameter
of laser spot according to the requirements of applications.

##  3 Bunch length measurement

As a relativistic electron bunch travels along the vacuum channel in the tube,
it drives coherent Cherenkov radiation wakefields[20] that are confined to a
discrete set of modes due to the waveguide boundaries. The power (energy) at
certain modes is correlative to the length of electron bunch. So the bunch
length can be measured by observing the frequency spectrum of the wakefield
radiation. Fig. 6 is a sectional drawing of hollow cylindrical dielectric tube
coated on the outer surface with metal to form a dielectric-lined waveguide.
We now present a brief summary of analysis for the fields set up within this
structure.

The sectional drawing of a beam-driven cylindrical dielectric-lined waveguide

In the analysis, Gauss system of units is used. Fourier expansion of the
longitudinal fields in a circular cylindrical waveguide takes the form:

$\begin{array}[]{c}\left({\begin{array}[]{*{20}{c}}{{E_{z}}(r,t)}\\\
{{H_{z}}(r,t)}\\\
\end{array}}\right)=\frac{1}{{{{(2\pi)}^{3}}}}\int_{-\infty}^{\infty}{d\omega
d}k\times\\\ \sum\limits_{l=-\infty}^{\infty}{\exp\left[{-i(\omega
t-kz-l\theta)}\right]}\times\left({\begin{array}[]{*{20}{c}}{{e_{z}}(r)}\\\
{-i{h_{z}}(r)}\\\ \end{array}}\right)\\\ \end{array}$ (2)

Then $e_{z}(r)$ and $h_{z}(r)$ satisfy the Bessel’s equation

$\left[{\frac{{{d^{2}}}}{{d{r^{2}}}}+\frac{1}{r}\frac{d}{{dr}}+\left({k_{\bot}^{2}-\frac{{{l^{2}}}}{{{r^{2}}}}}\right)}\right]\left({\begin{array}[]{*{20}{c}}{{e_{z}}(r)}\\\
{{h_{z}}(r)}\\\ \end{array}}\right)=0$ (3)

In our case, only $TM_{0n}$ waveguide modes are excited. For these modes,
$l=0$ , and longitudinal magnetic field is zero. In the vacuum hole (
$0<r<R_{1}$ ), where $\varepsilon=\mu=1$ , the fields must be regular at $r=0$
, thus

${e_{z}}(r)={e_{z}}(0){I_{0}}(x)$ (4)

where $I_{0}(x)$ is the modified Bessel function, and
$x\equiv\left|{k_{\bot}^{\left(1\right)}}\right|r$ , and
${\left({k_{\bot}^{\left(1\right)}}\right)^{2}}\equiv{\left({\frac{\omega}{c}}\right)^{2}}-{k^{2}}<0$
. In the dielectric region( $R_{1}<r<R_{2}$ ) with outer conducting boundary,
where $\varepsilon=\varepsilon_{2}=\varepsilon_{r}$ , and $\mu=\mu_{2}$ , one
has ${e_{z}}({R_{2}})={e_{\theta}}({R_{2}})={h_{r}}({R_{2}})=0$ , these are
equivalent to simply ${e_{z}}({R_{2}})=d{h_{z}}({R_{2}})/dr=0$ . Thus

${e_{z}}(r)={e_{z}}(0){I_{0}}\left({{x_{1}}}\right)\frac{{{E_{0}}(y)}}{{{E_{0}}({y_{1}})}}$
(5)

where $x_{1}\equiv\left|{k_{\bot}^{\left(1\right)}}\right|R_{1}$ ,
${E_{0}}(y)\equiv{J_{0}}(y){N_{0}}({y_{2}})-{N_{0}}(y){J_{0}}({y_{2}})$ with
$y\equiv\left|{k_{\bot}^{\left(2\right)}}\right|r$ ,
$y_{1}\equiv\left|{k_{\bot}^{\left(2\right)}}\right|R_{1}$ ,
$y_{2}\equiv\left|{k_{\bot}^{\left(2\right)}}\right|R_{2}$ and
${\left({k_{\bot}^{\left(2\right)}}\right)^{2}}\equiv\varepsilon\mu{\left({\omega/c}\right)^{2}}-{k^{2}}>0$
, the $J_{0}$ and $N_{0}$ are ordinary Bessel functions of the first and
second kinds. The boundary condition at $r=R_{1}$ is that $e_{z}$ , $h_{z}$ ,
$e_{\theta}$ , $h_{\theta}$ are continuous, then we can get the equation

$\frac{{I_{{}_{0}}^{\prime}({x_{1}})}}{{{x_{1}}{I_{0}}({x_{1}})}}+\varepsilon\frac{{E_{{}_{0}}^{\prime}({y_{1}})}}{{{y_{1}}{E_{0}}({y_{1}})}}=0$
(6)

The $n_{th}$ root of this equation is the $k_{n}$ , and the $f_{n}=c\cdot
k_{n}/2\pi$ is the frequency of the $n_{th}$ mode radiation excited in the
structure.

Then, we go to the power solution of these excited modes when a charge bunch
traverses the structure, and the theory in Reference[21] is used. For $N=2$
concentric dielectric layers in the uniform cylindrical waveguide, the
orthonormality relation between any two modes can be written as:

$\sum\limits_{i=1}^{N=2}{\int_{{R_{i-1}}}^{{R_{i}}}{dr\cdot
r}}\left[{{\varepsilon_{i}}{e_{z,m}}\left(r\right){e_{z,n}}\left(r\right)+{\mu_{i}}{h_{z,m}}\left(r\right){h_{z,n}}\left(r\right)}\right]={C_{n}}{\delta_{mn}}$
(7)

where $C_{n}$ is the normalization constant to be used when a moving charge
bunch is the source of the fields. The $\overline{{P_{0n}}}$ is the power
radiated into the $TM_{0n}$ mode:

$\overline{{P_{0n}}}=-cq_{0}^{2}\beta\frac{{e_{z,n}^{2}\left(0\right)}}{{{C_{n}}}}\Theta(-s)\cdot
g({\sigma_{z}})$ (8)

where $q_{0}$ is the charge, $c\beta$ is the velocity of the electron,
$\Theta(-s)$ means the radiation is excited behind the electron,
$g({\sigma_{z}})$ is the form factor. For gaussian shape,
$g({\sigma_{z}})=\exp(-4{\pi^{2}}\sigma_{z}^{2}/\lambda_{n}^{2})$ , where
$\sigma_{z}$ is the rms length of the bunch, ${\lambda_{n}}=2\pi/{k_{n}}$ is
the wavelength of $TM_{0n}$ . For uniform shape,
$g({\sigma_{z}})={\sin^{2}}({k_{n}}\cdot\sqrt{3}{\sigma_{z}})/{({k_{n}}\cdot\sqrt{3}{\sigma_{z}})^{2}}$
. Fig. 7 shows the power of wakefield radiation as a function of frequency for
electron bunches with different rms lengthes, where the charge is 200 pC, the
beam energy is 5.55 MeV, the inner and outer radii of the structure are
$R_{1}=1.2~{}mm$ , $R_{2}=6.5~{}mm$ respectively. The material of the
dielectric is fused silica $\varepsilon_{r}=3.8$ .

The power of wakefield radiation as a function of frquency, for bunches with
different rms length. The black line is for gaussian shape, and the red line
is for uniform shape

The sketch of the measurement setup is shown in Fig. 8. The dielectric length
should be several centimeters to make sure that the Cerenkov wakefield
radiation dominates the transition radiation which is emitted as the bunch
enters or leaves the structure[22].

Sketch of the measurement setup

The radiation emits from the structure, and then is reflected parallel by a
parabolic mirror. The aperture in the mirror is used to ensure the passage of
electron beam, and the transition radiation generated by the passage of
electron beam through the aperture is weak and easy to calibrate. The filters
filtrate out the other radiation except for the radiation at specified
frequency. The detector can be a Schottky barrier diode, a golay cell, or a
bolometer. The precision of the diode is relatively low compared with the
golay cell and bolometer, and it works in certain bandwidth. The bolometer is
expensive and needs a cryogenically cooled environment. The golay cell is a
good choice with high precision and portability. Three or four filters will be
used to reconstruct the spectrum of the radiation, and the bunch length can be
measured. After calibration, just one filter is needed, and single shot
measurement can be achieved. Besides, the charge of the bunch can also be
concluded.

To ensure the passage of electron bunch, the dielectric structure should be
installed around the focal point of the solenoid. The transverse beam size
evolution along the longitudinal position is shown in Fig. 9.

Evolution of the transverse beam size along the longitudinal position

Considering the inner radius of the dielectric structure is 1.2 mm, the
electron beam can pass the measurement setup in a long range without electron
loss.

##  4 Summary and discussion

In this paper, we analyse the impact factors of the electron bunch length, and
draw the conclusion that to generate sub-picosecond electron bunch directly in
the photocathode rf gun, we need to improve the acceleration gradient of gun
as high as possible, restrain the growth of transverse emittance through laser
shaping technique, and carefully tune the energy (bunch charge) and diameter
of the laser spot according to the requirements of applications. In order to
measure the electron bunch length, the coherent Cherenkov radiation wakefield
is also analysed, which is excited by the relativistic electron bunch
traveling through a hollow cylindrical dielectric element. Based on the
analysis, a nondestructive technique for the measurement of bunch length is
reported. The advantage of this technique is routine to calibrate, and can be
single shot measurement.

To replace the dielectric element used in the measurement setup with a
redesigned one, narrow band THz radiation (around 0.3 THz) with high peak
power (hundreds KW) can be excited, and this can be applied to table-top THz
source.

## References

  * [1]
  * [2] Carr G, Martin M, McKinney W, et al. Nature (London), 2002, 420: 153
  * [3] Korbly S E, Kesar A S, Sirigiri J R, et al. Phys. Rev. Lett, 2005, 94: 054803
  * [4] Cook A M, Tikhoplav R, Tochitsky S Y, et al. Phys. Rev. Lett, 2009, 103: 095003
  * [5] Huang W H, He X Z, HUANG G, et al. High Energy Physics and Nuclear Physics, 2004, 28(4): 446
  * [6] Li R K, Tang C X, Huang W H, et al. Chinese Physics C, 2009, 33(Suppl.II): 165
  * [7] Shi J R, Chen H B, Tang C X, et al. Chinese Physics C, 2009, 32(10): 837
  * [8] Wang D X, Krafft G A, Sinclair C K. Phys. Rev. E, 1998, 57: 2283
  * [9] Berden G, Jamison S P, Macleod A M, et al. Phys. Rev. Lett, 2004, 93: 114802
  * [10] Lai R, Happek U, Sievers A J. Phys. Rev. E, 1994, 50: R4294
  * [11] Xiang D, Yang X F, Huang W H, et al. Chin. Phys. Lett, 2008, 25(7): 2440
  * [12] Shchelkunov S V, Marshall T C, Hirshfield J L, et al. Phys. Rev. ST Accel. Beams, 2005, 8: 062801
  * [13] Clayton C E, Serafini L. IEEE Transactions on Plasma Science, 1996, 24(2): 400
  * [14] Serafini L, Zhang R, Pellegrini C. Nucl. Instr. and Meth. A, 1997, 387: 305
  * [15] Flottmann K. ASTRA user manual, http://www.desy.de/~mpyflo/
  * [16] Huang W H, private communication
  * [17] Qian H J, Chen H B, Du Y C, et al. Proceedings of IPAC2011, 2011, 3170
  * [18] Wang X J, Qiu x, Ben-Zvi I. Phys. Rev. E, 1996, 54(4): R3121
  * [19] He Z G, Jia Q K. Proceedings of PAC2011, 2011, 1196
  * [20] Cherenkov P A. Phys. Rev, 1937, 52: 378
  * [21] Park S Y, Hirshfield J L. Phys. Rev. E, 2000, 62(1): 1266
  * [22] Wang C, Hirshfield J L, Fang J M, et al. Phys. Rev. ST Accel. Beams, 2004, 7: 051301


# Pleobot: a modular robotic solution for metachronal swimming

Sara Oliveira Santos  Brown University, Center for Fluid Mechanics, School of
Engineering, Providence, 02912, USA  Nils Tack  Brown University, Center for
Fluid Mechanics, School of Engineering, Providence, 02912, USA  Yunxing Su
Brown University, Center for Fluid Mechanics, School of Engineering,
Providence, 02912, USA  Francisco Cuenca-Jiménez  Universidad Nacional
Autónoma de México, Circuito Interior s/n, Engineering, Coyoacán, 04510,
México  Oscar Morales-Lopez  Universidad Nacional Autónoma de México, Circuito
Interior s/n, Engineering, Coyoacán, 04510, México  P. Antonio Gomez-Valdez
Universidad Nacional Autónoma de México, Circuito Interior s/n, Engineering,
Coyoacán, 04510, México  Monica M. Wilhelmus  Brown University, Center for
Fluid Mechanics, School of Engineering, Providence, 02912, USA
<EMAIL_ADDRESS>
###### Abstract

Metachronal locomotion is a widespread swimming mode used by aquatic swarming
organisms to achieve performance and maneuverability in the intermediate
Reynolds number regime. Our understanding of the mechanisms driving these
abilities is limited due to the challenges of studying live organisms. Designs
inspired by nature present an approach for developing small and maneuverable
underwater self-propelled robots. Here, we present the design, manufacture,
and validation of the _Pleobot_ –a unique krill-inspired robotic swimming
appendage constituting the first platform to study metachronal propulsion
comprehensively. Our methods combine a multi-link 3D printed mechanism with
active and passive actuation of the joints to generate natural kinematics.
Using force and fluid flow measurements in parallel with biological data, we
show the link between the flow produced by the appendage and thrust. Further,
we provide the first account of a leading-edge suction effect that contributes
to lift during the power stroke. The repeatability and modularity of the
_Pleobot_ enable the independent manipulation of particular motions and traits
to test hypotheses central to understanding the relationship between form and
function. Lastly, we outline future directions for the _Pleobot_ , including
adapting morphological features. We foresee a broad appeal to a wide array of
scientific disciplines, from fundamental studies in ecology, biology, and
engineering, to developing new platforms for studying oceans across the solar
system.

## Introduction

Nature-inspired robotics have become valuable tools in biology and ecology as
they allow better-controlled comparative studies than traditional organismal
methods [1]. Independent manipulation of a particular trait or behavior of
interest enables testing hypotheses central to exploring the underlying
mechanisms driving maneuverability and performance. While several test
platforms have shed light on important aspects of biological underwater
propulsion [2, 3, 4, 5], they have also produced practical engineering
solutions for unmanned underwater vehicles (UUV). Previous robotic solutions
for engineering problems, including tunas [6], scallop [7], and dolphins [8]
demonstrate improved propulsive efficiency and optimal operation in a range of
environments. However, because of experimental or technical challenges, most
advances have been for relatively high Reynolds numbers (Re). Consequently,
despite the tremendous diversity of species inhabiting the moderate realm, we
still need a comprehensive understanding of aquatic propulsion and its derived
engineering applications at low-to-intermediate Re.

Metachronal swimming is a ubiquitous swimming mode among many of the most
abundant invertebrate species inhabiting the transitional flow regime (Re $<$
1000), such as copepods and krill. This propulsive gait is characterized by
the sequential beating of several closely spaced swimming appendages, phase-
shifted in time, in a tail-to-head traveling wave [9]. Metachrony enables
large-scale behaviors in krill, particularly the ability to perform diel
vertical migrations (DVM) from the sea surface to oxygen minimum zones located
up to 1 km at depth [10, 11]. Efficient swimming is achieved by leveraging
drag during properly tuned power-and-recovery strokes [9, 12]. During the
power stroke, the swimming appendages (pleopods) move opposite the swimming
direction while maximizing their surface area. This generates drag in each
pleopod, producing thrust to propel the organism. Conversely, during the
recovery stroke, the profile area decreases, thus reducing the drag on the
appendages and creating a net thrust force sufficient to overcome the drag on
the body. Previous works showed drag-based metachronal swimming is more
effective than lift-based propulsion at intermediate Re for accelerating,
braking, and turning because it generates significant thrust over short
periods [13, 9, 14, 15].

Fluid flow experiments on live krill have established the foundations of our
understanding of the effect of fluid-structure interactions on the far-field
flow [16, 17]. Murphy et al. characterized the swimming kinematics of live
krill and linked the success of the propulsion system to appendage morphology,
stroke kinematics, and the resulting hydrodynamic effects [16]. Tip vortices
forming on the pleopods of krill [17] and mantis shrimp [18] were shown to be
central in the production of thrust [19, 20]. Kim and Gharib extended these
results using idealized pleopod shapes with simplified kinematics. They
demonstrated that the area over the surface of the appendage enclosed by the
tip vortex, specifically, modulates thrust production [19]. However,
challenges in characterizing the near-field flow and measuring the forces
generated during swimming have hindered efforts to establish the relationship
between fluid dynamics and force distribution (i.e., thrust and lift) in
intermediate Re metachronal swimmers.

Simplified robotic models and numerical simulations have complemented these
efforts, shedding light on the role of varying Re, phase lag, and appendage
spacing on the hydrodynamics [21, 22, 18, 23, 24, 25]. For instance, Ford et
al. showed that at low Re, the individual jets produced by adjacent pleopods
do not interact due to viscous dissipation, but at Re around 800, they form a
near-steady jet [22]. The outcome of this difference is an increase in
vertical and horizontal momentum that enables krill to generate thrust and
lift forces required for locomotion and hovering. Phase lag also emerges as an
important factor contributing to metachronal swimming efficiency, as it yields
near-maximal efficiency and thrust, and achieves higher average body velocity
when compared to synchronous motions [21, 22, 24, 26]. Likewise, appendage
spacing is a critical morphological factor linked to stroke kinematics.
Previous works demonstrated that greater swimming performance is achieved when
the ratio of appendage spacing to appendage length is less than one [27]. This
was hypothesized to be the direct consequence of inter-pleopod fluid flow
interactions enhancing the vortex strength and circulation during the power
stroke [18].

While the insights gathered through simplified systems are invaluable, we
still need a unifying theory of drag-based metachronal swimming to explain the
relationship between appendage kinematics, fluid flows, and force production.
For instance, even the production of lift needed for hovering still needs to
be explored. Adding complexity to better match the natural biological and
hydrodynamic phenomena is necessary to elucidate the effects of morphological
and kinematic characteristics on the near- and far-field flow. We should then
adopt an integrative framework allowing for the comprehensive study of the
link between locomotor kinematics, form, and function.

Motivated by a robotics-inspired biology approach to address this problem, we
present the first fully articulated, multi-link, metachronal robotic appendage
reproducing the swimming kinematics of krill. Our novel design achieves active
control of both the proximal and distal appendages and the passive out-of-
plane actuation of the two rami of the distal appendages through the
interaction between the pleopod and the surrounding fluid. In addition, we
implemented the main morphological features of krill pleopods through the
assembly of 3D-printed modules. Using force and fluid flow measurements, we
quantify the effects of particular flow features on thrust and lift production
and disentangle their role during the power and recovery phases of a stroke.
We use biological data to validate our model as a suitable test platform. In
this study, we leverage the _Pleobot_ to study the near-field flow of a
metachronal pleopod for the first time and provide evidence for a previously
unreported leading-edge suction effect that contributes to lift during the
power stroke. The repeatability and modularity of the _Pleobot_ enables a
range of configurations and motion programs to investigate metachrony across
Re, taxa, and morphological attributes. As such, we anticipate the _Pleobot_
will broadly appeal to a wide array of scientific disciplines, such as
ecology, biology, and engineering.

## Results

### Pleobot kinematics

The swimming kinematics of the _Pleobot_ are characterized by three main
angles: 1) $\alpha$, the angle between the body axis and the protopodite; 2)
$\beta$, the angle between the protopodite and the biramous distal appendage
(including the endopodite and exopodite); and 3) $\gamma$, the angle between
the endopodite and exopodite (Fig. 1A-D). The cupping angle, $\zeta$, is also
implemented in the _Pleobot_ as seen in _Euphausia superba_ (Fig. 1E,F). The
resulting krill-inspired model (Fig. 1G) uses a train of gears to actuate both
$\alpha$ and $\beta$ kinematics (Fig. 1B). $\gamma$ is actuated passively due
to hydrodynamic forces during the power and recovery strokes (Fig. 1D), and
the $\zeta$ angle of the exopodite is fixed (Fig. 1F).

The kinematics of the _Pleobot_ were validated against organismal data [16] by
tracking the motion of the robotic system over several consecutive cycles (n =
5). The cycle-averaged $\alpha$ and $\beta$ angles over time were in line with
kinematics data for _E. superba_ (Fig. 2A,B) [16]. The slight deviation of the
$\beta$ angle from krill data was associated with the tolerance between gears,
which is directly proportional to the number of units in the train. The motion
of the passively actuated exopodite ranged from 14∘ to 95∘, with abduction
occurring at the end of the return stroke and adduction coinciding with the
start of the recovery strokes (Fig. 1). Comparative values over a complete
beat cycle have not been reported for _E. superba_ , but Murphy and
collaborators[16] described $\gamma$ oscillating between 0∘ and 77∘. The
kinematic analysis of the _Pleobot_ demonstrates it can perform the prescribed
motion program with minimal error and agrees with the values reported for live
krill. As such, _Pleobot_ is an accurate and highly repeatable robotic analog
of the swimming appendages of live metachronal organisms [16].

### Flow visualization

PIV measurements were performed in four different planes to evaluate the flow
field (Fig. 3A). At the beginning of the power stroke (t/T, the non-
dimensional stroke time, is 0), a vortex forms at the tip of the endopodite
accompanied by a bound vortex uniformly distributed along most of the length
of this segment (Fig. 3E). This pattern agrees with biological data acquired
using live ghost shrimp (Fig. 3F). Meanwhile, the exopodite abducts as
$\alpha$ increases and $\beta$ decreases (Fig. 2). During the first half of
the power stroke (before reaching t/T = 0.25), the exopodite maintains a low
angle of attack (AoA, measured relative to the body axis, which is horizontal,
in a clockwise direction), promoting the formation of an attached leading-edge
vortex (LEV) along the anterior face (in the swimming direction) of the
exopodite (Fig. 3B). The LEV remains attached until the exopodite reaches AoA
= 73o ± 1.9o at t/T = 0.3. As $\alpha$ angle increases, the AoA of the
exopodite also increases to become vertical relative to the swimming direction
(AoA = 90o) during the second half of the power stroke (Fig. 4I). Given that
the pleopod decelerates during this phase, the LEV separates from the
exopodite (Figs 3C, 4I) and the bound vortex of the endopodite sheds at the
tip (Figs 3E,F, 4E). The combined effect of the flow entrained by the two
counter-rotating vortices results in the production of a large posteriorly
oriented downward jet.

Upon initiating the recovery stroke, the exopodite adducts and overlaps with
the endopodite to effectively reduce the overall profile area of the pleopod
subject to the flow (Figs 2C,E, 3D, 4C). Contrary to the power stroke, this
configuration generates a pair of counter-rotating side edge vortices that
induces spanwise flow in the swimming direction (Figs 3D, 4F,G).

### Hydrodynamic forces

Force measurements were performed to evaluate the contribution of lift and
thrust forces during steady forward swimming (Fig. 4). The lift ($C_{L}$) and
thrust ($C_{T}$) coefficients were computed for a 20x-scaled endopodite-and-
exopodite model. At the beginning of the power stroke (t/T = 0), there is a
significant amount of thrust force, likely because of the incoming spanwise
flow produced behind the pleopod at the end of the previous beat (Fig. 4G).
During the power stroke, thrust and lift strongly correlate with the abduction
of the exopodite. The forces on the endopodite decompose in the thrust and
lift directions, and the LEV forming along the exopodite drives this
association. Throughout the first third of the power stroke (t/T $\sim$ 0.16),
the lift coefficient increases sharply to reach its global maximum ($C_{L}$ =
0.44) as the exopodite is abducting and produces an attached LEV (Fig.
4B,C,H). In contrast, while thrust is being produced, it remains mostly
constant during this phase, and only increases significantly after lift peaks
and decreases (Fig. 4A). This inverse relationship coincides with the AoA of
the exopodite shifting to a more vertical orientation that causes the LEV to
shed. It also results from the increased contribution of the endopodite to
thrust, given the overall vertical orientation of the entire pleopod (see Fig.
4A,B). During the second half of the power stroke, thrust decreases sharply as
the bound vortex of the endopodite moves toward the tip and forms a stopping
vortex (Fig. 4E). Given the increase of the effective angle of the pleopod
relative to the flow (angle $\Psi$, the overall angle between the body axis
and the frontal face of the endopodite; see Fig. 4), the endopodite and the
exopodite produce a downward force inducing minimum lift. The subsequent
increase in lift force is likely due to the occurrence of spanwise flow
entrained by the stopping vortex that pushes against the anterior face of the
decelerating endopodite and exopodite (Fig. 4E).

During the first 2/3 of the recovery stroke, thrust is negative, suggesting
this phase of a beat is mostly dominated by drag forces. The thrust
coefficient is lowest around the mid-recovery stroke (Fig. 4A) when the
pleopod accelerates. Thrust increases gradually to become positive, thus
indicating some thrust is produced even during the recovery stroke. After the
start of the recovery stroke, lift drops sharply to another minimum with
decreasing $\psi$ angle (Fig. 4B). Lift rises again at the end of the recovery
stroke when the appendage decelerates. Increasing thrust and lift coincide
with the occurrence of spanwise flow behind the overlapping endopodite and
exopodite (Fig. 4G).

## Discussion

In recent years, the growing need for maneuverable AUVs for underwater
exploration has galvanized approaches inspired by nature [28, 29, 30, 31, 32].
In particular, the discovery of extraterrestrial oceans motivates the
development of novel robotic platforms that will likely require the
efficiency, versatility, and maneuverability of metachronal swimming. As such,
the _Pleobot_ lays down the foundations for the upcoming work on metachronal
swimming to inform the design of underwater explorers.

_Pleobot_ represents the first fully articulated robot to study metachronal
swimming. The design constitutes a modular, 3D-printed platform incorporating
the kinematic and morphological characteristics of krill – a marine organism
observed to form large aggregations and migrate hundreds of meters into the
ocean. While krill kinematics have been explored in great detail, a
quantitative analysis of the fluid flow near individual appendages is needed
to link kinematics to hydrodynamics and understand force production. _Pleobot_
is the first system to disentangle the roles of kinematic parameters,
morphology, and near-field hydrodynamics on thrust and lift generation at the
scale of a single appendage.

Our unique study shows that _Pleobot_ serves as an analog system to
investigate a wide range of parameters and characteristics that collectively
contribute to the efficiency and flow characteristics of metachronal
propulsion. Although krill has been used as a model organism in biological,
computational, and robotic studies, the focus has mostly been on the flow
field of several appendages rather than the singular mechanisms for the
generation of thrust and lift. Conclusions regarding fluid-structure
interactions are thus solely possible by measuring the near field flow and by
introducing complexities in the model that enable independent testing of
hypotheses relevant to the ecology and behavior of the natural systems (i.e.,
maneuverability, migrations, fast swimming).

In our study, the analysis of force and PIV data allows us to evaluate the
contribution of specific flow features to the generation of lift and thrust in
_Pleobot_. At the beginning of the power stroke, the endopodite contributes to
thrust and lift through the formation of tip and bound vortices. Similarly,
the exopodite also contributes to the generation of both forces. Creation of
lift during steady forward swimming is necessary for krill because they are
negatively buoyant and must create lift to maintain their position in the
water column [17]. As the stroke progresses, the exopodite becomes more
vertically oriented and its contribution to lift decreases, while its
contribution to thrust increases (Fig. 4A,B). This represents a trade-off
between the generation of lift and thrust production. While the primary
purpose of thrust in steady swimming is to propel the organism forward, lift
is needed to guarantee the vertical position in the water column. Our data
suggest that rather than being determined by behavioral processes, thrust and
lift production are constrained to a relatively narrow range of kinematics and
morphological criteria. Behavioral changes in the gait and beat frequency
invariably modulate these forces. However, the dependency on the $\psi$ angle
of the passive abduction of the exopodite and its orientation in the flow
demonstrates that a limited set of system requirements are fundamentally
important for thrust and lift production in metachronal propulsion.

Notably, we reveal the existence of a LEV along the exopodite during the power
stroke. Force measurements show that lift is enhanced by the formation of this
attached LEV on the exopodite during the power stroke. Furthermore, this
vortex delays the loss in lift forces slightly to maintain positive lift
throughout part of the power stroke. This lift-based mechanism likely
supplements the effects of the drag-based mechanism stemming from the pleopod
pushing on water. These results in the near field illustrate the importance of
the _Pleobot_ , not only as a platform to engineer solutions for ocean
exploration in the intermediate Reynolds number regime but also as a framework
to understand nature. Force measurements on a biologically-inspired appendage
highlight the importance of understanding the vortex dynamics around the
pleopod and reveal that both lift and thrust are generated simultaneously.

The incorporation of natural kinematics and morphological traits into the
_Pleobot_ provides promising results and incentives for increasing the
complexity of our robotic model and elucidate additional mechanisms
contributing to efficiency in krill and other metachronal swimmers. Markedly,
the presence of hair-like setae around the endopodite and exopodite, and
flexibility of the propulsors have the most potential for altering the flow
field around the pleopods. Setae are defining features in metachronal
organisms in the low-to-intermediate Re that increase the surface area in
contact with water to enhance thrust [33, 9, 34]. Flexibility has been shown
to be an important characteristic of propulsors in smoothing out thrust peaks
to generate nearly constant thrust during the power stroke [19]. _Pleobot_
provides a platform to evaluate the importance of these additional
morphological features. Moreover, its modular design can be easily replicated
and grouped in several units to investigate mechanisms such as constructive
pleopod interactions (Fig. 5). For example, the proximity of several
appendages beating at a phase lag promotes fluid flow interactions that likely
contribute to enhanced thrust [35, 18]. The Pleobot emerges as an ideal
platform to quantify these effects.

Finally, while organismal investigations form the core of our knowledge about
metachronal swimming [13, 18, 16, 17], working with live animals provides only
a partial understanding of the biomechanical and hydrodynamic mechanisms
leading to efficiency and performance. Recent mechanical and robotics
approaches allowing repeatable measurements highlight the valuable
contribution of low-cost, controllable test platforms for biological studies
[27, 22, 35]. Such a comparative approach will be central to developing a
unifying theory of metachronal swimming. Manufacturing of the _Pleobot_ is
accessible to both institutions and individuals who want to experiment (see
`oliveira:zenodo:2022)`). Its modular design allows for the quick prototyping
of different appendage shapes and sizes. Furthermore, the control of both its
proximal and distal segments facilitates the study of a wide range of systems,
from small copepods, to krill, and larger species like mantis shrimp and
lobsters. The _Pleobot_ can thus be employed to facilitate the comprehensive
characterization of metachronal, drag-based propulsion and help establish a
unifying theory for this locomotor mechanism.

## Methods

### Modeling of krill appendage kinematics

To reproduce the kinematics of _Euphausia superba_ , as reported by Murphy et
al.[16], we use values for $\alpha$ and $\beta$ to find kinematic
relationships for Pleobot. The Pleobot is dynamically scaled based on the
Reynolds number. This dimensionless parameter is defined using the velocity
and the length scale of the appendages and is approximately 600 for live
krill:

$\displaystyle Re=\frac{U_{tip}*L}{\nu}=\frac{2\theta nL^{2}}{\nu}\;$ (1)

where $L$ is the length of the pleopod, $U_{tip}$ the average velocity at its
tip, and $\nu$ the kinematic viscosity of the fluid. $U_{tip}$ is calculated
using the stroke amplitude $\theta$, the frequency $n$, and the length $L$ as
$2\theta nL$.

The locomotor system is implemented by a multi-link mechanism using a
transmission gear box actuating both the proximal and distal segments
(corresponding to $\alpha$ and $\beta$) with independently controlled servo
motors. Angle $\gamma$ represents out-of-plane motion, posing challenges for
active actuation and is thus passively actuated via the hydrodynamic
interaction between the structure (pleopods) and the fluid.

The kinematic relationships governing the positions of the gears are[36, 37]:

$\displaystyle\frac{\Delta\phi_{k+1}-\Delta\phi_{k}}{\Delta\phi_{k-1}-\Delta\phi_{k}}=-\frac{N_{k-1}}{N_{k+1}}=-\dfrac{r_{k-1}}{r_{k+1}}\;$
(2)

This equation gives the relationship between the driving and driven gears,
$\Delta\phi_{k+1}$ and $\Delta\phi_{k-1}$, respectively, and the link
connecting both, $\Delta\phi_{k}$ (Supplementary figure 1), where $N$ and $r$
are the number of teeth and radius of the corresponding gear. The position of
the pleopod, $\mathbf{R}_{pi}$, is governed for a given appendage $i$ as
(Supplementary figure 2,B):

$\displaystyle\mathbf{R}_{pi}=\mathbf{R}_{1i}+\mathbf{R}_{2i}\;$ (3)

where

(6)

The local reference frame vectors for the protopodite and exopodite are
$\mathbf{r}_{1i}$ and $\mathbf{r}_{2i}$. Similarly, $\mathbf{R(\theta_{1i})}$
and $\mathbf{R(\theta_{2i})}$ are the global rotation matrices and $x_{1i}$
and $x_{2i}$ the lengths (Supplementary figure 2,A).

Finally, the angular displacements of the gears are obtained using the
relationships described in Equation 2:

$\frac{\Delta\psi_{2i}-\Delta\theta_{1i}}{\Delta\psi_{1i}-\Delta\theta_{1i}}=-r_{e1},\\\
\frac{\Delta\psi_{3i}-\Delta\theta_{1i}}{\Delta\psi_{2i}-\Delta\theta_{1i}}=-r_{e2},\\\
\frac{\Delta\psi_{4i}-\Delta\theta_{1i}}{\Delta\psi_{3i}-\Delta\theta_{1i}}=-r_{e3}$
(7)

where $\psi$ is the rotation of the gear along link 1 (endopodite), $\theta$
is the angle of the first link measured from the global reference frame
(Supplementary figure 2,A), and $r_{e}$ is the gear ratio. Equations (3) and
(5) allow for the solution of direct kinematics of the mechanism by using
$\alpha_{i}\leavevmode\nobreak\ $and $\beta_{i}$, to calculate $\psi_{1i}$,
which moves the second bar (endopodite)(Supplementary figure 2,C). The first
bar (protopodite) is moved by angle $\alpha_{1i}$ (Supplementary figure 2,A).
The metachronal trajectory at the end of pleopod 1, is given by
$\mathbf{R}_{p1}=\left[x_{p1},y_{p1}\right]^{T}$.

Manufacturing was completed by 3D printing, keeping small tolerances to reduce
vibrations and loss of movement of the gear train. The supports house the
servos and gears, and were designed to be an above-water structure. Bearings
are used to reduce friction between gears and the axes. The transmission has
an amplification of 2.5 to achieve the desired angular speed for the links.

### Exopodite and endopodite

As _E. superba_ propels forward, the endopodite and exopodite abduct and
adduct to effectively change the profile area of the appendages to generate
thrust and reduce drag. This motion is characterized by $\gamma$, the angle
between the exopodite and the endopodite. In the Pleobot, this process is
actuated on the horizontal plane, lateral to the protopodite.

In forward-swimming krill, exopodite abduction occurs at the beginning of the
power stroke, and reaches a maximum of 77 ∘ [16]. This motion induces cupping
of the appendages creating a V-shaped structure [17], reminiscent of those
observed in swimming fish that have been shown to produce greater thrust
compared to flat fins (e.g., see [2]). The exopodite and endopodite adduct
during the recovery stroke.

Appendage cupping forms angle $\zeta$ between the endopodite and the
exopodite. Photographic evidence was used to quantify $\zeta$ by measuring the
angle between the midplanes of the endopodite and the exopodite (Fig. 1E) .
The Pleobot was set to match the mean quantified value of 37∘.

### Robotic design

Modular CAD designs were printed with polylactic acid (PLA) using a Prusa i3
MKS3+ 3D printer (Prusa Research, Prague, Czech Republic) for fast
prototyping. Each pair of appendages is actuated by two servos (HS-5087MH,
Hitec RCD, San Diego, CA, USA), controlled by a microcontroller (ELEGOO Mega
2560, Elegoo Industries, Shenzhen, China) programmed using Simulink
(MathWorks, Natick, MA, USA) via two repeating sequence interpolated blocks,
one for $\alpha$ and one for $\beta$, prescribing the angles adapted from
[16]. The full CAD library and assembly as well as the list of purchased
components can be accessed in the open-access repository by Oliveira Santos et
al. (see `oliveira:zenodo:2022)`.

The swimming kinematics of the Pleobot were analyzed and compared to those
reported for live krill [16]. One robotic appendage was tethered to a traverse
beam and submerged in a glycerin-water mixture at room temperature (60$\%$
glycerin and 40$\%$ water). Appendage motion was recorded at 125 fps using a
scientific camera (FASTCAM MINI WX, Photron, 2048 pixels x 2048 pixels). Black
markers on the surface of the robotic pleopod were tracked, both on the
protopodite and the exopodite by digitizing the video recordings via DLTdv8
for MATLAB using automatic point tracking [38].

### Flow field measurements

The flow field around the beating appendage was measured using 2D PIV. The
experiments were carried out in a cubic tank of 30 cm by side using a high-
speed camera (FASTCAM MINI WX, Photron, 2048 pixels x 2048 pixels) at 125
frames per second, a Nikon lens (Nikon AF-S VR Micro-NIKKOR 50 mm), and a
continuous laser (Laserglow, 1 W at 532 nm) with a cylindrical lens to create
a laser sheet. The flow was seeded with 10 $\mu$m particles (Dantec Dynamics,
Skovlunde, Denmark) (Supplementary Figure 3). The field of view is 174 mm x
174 mm, resulting in a spatial resolution of 0.085 mm per pixel. The velocity
fields were calculated using DaVis 10 (LaVision) with decreasing interrogation
window sizes (initial 48 pixels x 48 pixels, 50$\%$ overlap, 1 pass, and final
32 pixels x 32 pixels, 50$\%$, 3 passes). Standard vector post-processing was
performed to remove outliers in the flow field. The vorticity field was
calculated from the velocity field.

### Force measurements

The force measurements were conducted on a scaled-up (20x) model of the distal
appendage by using a 6-axis force transducer (Nano 17 F/T transducer, ATI) set
up in a cubic tank of 30 cm by side (Supplementary Figure 4). The Reynolds
number was matched to ensure dynamically scaled experiments. To this end, the
beating frequency and the fluid viscosity were modified from the baseline. A
glycerin-water mixture (60$\%$ glycerin and 40$\%$ water) was used with a
dynamic viscosity of 9 cSt, measured with a standard rheometer (Ares-G2, TA
Instruments). The appendage was designed to be neutrally buoyant in the water-
glycerin mixture.

The force transducer was mounted in the axis of rotation of the bi-ramous
pleopod, such that the x-axis was aligned with the pleopod, and the y-axis was
perpendicular to it. The forces measured, parallel ($F_{\parallel}$) and
perpendicular ($F_{\bot}$) to the pleopod, are decomposed into lift and thrust
by taking into account the orientation of this section in the body via the
angle $\Psi$, the overall angle between the body axis and the endopodite,
defined as $\alpha+180^{o}-\beta$. The lift and thrust coefficients ($C_{L}$
and $C_{T}$, respectively) are calculated using the lift and thrust forces
measured from the decomposed components (Equation 6). Using the instantaneous
force along a cycle, averaged over five cycles, we calculate lift and thrust
as:

(10)

and $C_{L}$ and $C_{T}$ as

(12)

where $\rho$ is the density of the glycerin-water mixture, $U$ is the
appendage tip velocity, and A the area of the abducted appendage.

### Organism experiments

Live ghost shrimp _Palaemonetes paludosus_ were used to validate the
hydrodynamics and $\gamma$ angle kinematics of the _Pleobot_. Ghost shrimp are
not obligate swimmers like krill, but they are appropriate analogs for several
reasons. First, they have a ubiquitous reliance on metachronal swimming.
Second, they have comparable morphology and swimming kinematics (i.e.,
parameters related to their exopodites, endopodites, and setae) to krill.
Finally, they are similar in body size and anatomy.

Specimens (n = 4) were tethered to a vertical wire using cyanoacrylate applied
to the dorsal surface of their carapace and hung in a filming vessel (15 × 10
× 5 cm3) brought to 100% oxygen saturation. We performed bright-field PIV by
seeding the water with 10 $\mu$m particles (Dantec Dynamics, Skovlunde,
Denmark) and back-lighting with a fiber optic illuminator (Fiber-Lite MI-152,
Dolan-Jenner Industries, Boxborough, MA, USA) coupled with a collimating lens
(N-BK7 Plano-Convex Lens, Thorlabs, Newton, NJ, USA) [39]. We recorded the
lateral and posterior views using a high-speed digital video camera (Fastcam
Nova R2, Photron, Tokyo, Japan) at 2000 fps and a resolution of 2048 × 1472
squared pixels.

Velocity vectors were calculated using the DaVis 10 software package
(LaVision, Göttingen, Germany). Image pairs were analyzed with three passes of
overlapping interrogation windows (75%) with decreasing size (96 × 96 to 64 ×
64 squared pixels). All frames were used for analysis, yielding a time
separation between frames of 0.5 ms. Masking of the body and pleopods before
image interrogation confirmed the absence of surface artifacts in the PIV
measurements.

## References

  * [1] Lauder, G. V.  Robotics as a comparative method in ecology and evolutionary biology.  _Integrative and Comparative Biology_ 62, 721–734 (2022).
  * [2] Esposito, C. J., Tangorra, J. L., Flammang, B. E. & Lauder, G. V.  A robotic fish caudal fin: effects of stiffness and motor program on locomotor performance.  _Journal of Experimental Biology_ 215, 56–67 (2012).
  * [3] Tytell E, L. J.  Biorobotic insights into neuromechanical coordination of undulatory swimming.  _Science Robotics_ 2, eabk0620 (2021).
  * [4] Howe, S., Bryant, K., Duff, A. & Astley, H.  Testing the effects of body depth on fish maneuverability via robophysical models.  _Bioinspiration & Biomimetics_ 17, 016002 (2021).
  * [5] Lauder, G. V. & Tangorra, J. L.  Fish locomotion: biology and robotics of body and fin-based movements.  _Robot fish: bio-inspired fishlike underwater robots_ 25–49 (2015).
  * [6] Zhu, J. _et al._ Tuna robotics: A high-frequency experimental platform exploring the performance space of swimming fishes.  _Science Robotics_ 4, eaax4615 (2019).
  * [7] Robertson, M. A., Efremov, F. & Paik, J.  Roboscallop: a bivalve inspired swimming robot.  _IEEE Robotics and Automation Letters_ 4, 2078–2085 (2019).
  * [8] Wu, Z., Yu, J., Yuan, J. & Tan, M.  Towards a gliding robotic dolphin: design, modeling, and experiments.  _IEEE/ASME Transactions on Mechatronics_ 24, 260–270 (2019).
  * [9] Vogel, S.  Life in moving fluids: The physical biology of flow-revised and expanded second edition.  _Princeton University Press_ (2020).
  * [10] Bianchi, D., Stock, C., Galbraith, E. D. & Sarmiento, J. L.  Diel vertical migration: Ecological controls and impacts on the biological pump in a one-dimensional ocean model.  _Global Biogeochemical Cycles_ 27, 478–491 (2013).
  * [11] Bianchi, D. & Mislan, K.  Global patterns of diel vertical migration times and velocities from acoustic data.  _Limnology and Oceanography_ 61, 353–364 (2016).
  * [12] Kils, U.  The Swimming behavior, swimming performance and energy balance of antarctic krill, _Euphausia superba_.  _BIOMASS Scientific Series No. 3_ (1981).
  * [13] Byron, M. L. _et al._ Metachronal motion across scales: current challenges and future directions.  _Integrative and Comparative Biology_ 61, 1674–1688 (2021).
  * [14] Vogel, S.  Comparative biomechanics: life’s physical world princeton university press.  _Princeton, New Jersey, USA_ (2003).
  * [15] Walker, J. A. & Westneat, M. W.  Mechanical performance of aquatic rowing and flying.  _Proceedings of the Royal Society of London. Series B: Biological Sciences_ 267, 1875–1881 (2000).
  * [16] Murphy, D., Webster, D., Kawaguchi, S., King, R. & Yen, J.  Metachronal swimming in antarctic krill: Gait kinematics and system design.  _Marine biology_ 158, 2541–2554 (2011).
  * [17] Murphy, D. W., Webster, D. R. & Yen, J.  The hydrodynamics of hovering in antarctic krill.  _Limnology and Oceanography: Fluids and Environments_ 3, 240–255 (2013).
  * [18] Garayev, K. & Murphy, D. W.  Metachronal swimming of mantis shrimp: kinematics and interpleopod vortex interactions.  _Integrative and Comparative Biology_ 61, 1631–1643 (2021).
  * [19] Kim, D. & Gharib, M.  Characteristics of vortex formation and thrust performance in drag-based paddling propulsion.  _Journal of Experimental Biology_ 214, 2283–2291 (2011).
  * [20] DeVoria, A. C. & Ringuette, M. J.  Vortex formation and saturation for low-aspect-ratio rotating flat-plate fins.  _Experiments in fluids_ 52, 441–462 (2012).
  * [21] Alben, S., Spears, K., Garth, S., Murphy, D. & Yen, J.  Coordination of multiple appendages in drag-based swimming.  _Journal of the Royal Society Interface_ 7, 1545–1557 (2010).
  * [22] Ford, M. P., Lai, H. K., Samaee, M. & Santhanakrishnan, A.  Hydrodynamics of metachronal paddling: Effects of varying reynolds number and phase lag.  _Royal Society open science_ 6, 191387 (2019).
  * [23] Granzier-Nakajima, S., Guy, R. D. & Zhang-Molina, C.  A numerical study of metachronal propulsion at low to intermediate reynolds numbers.  _Fluids_ 5, 86 (2020).
  * [24] Hayashi, R. & Takagi, D.  Metachronal swimming with rigid arms near boundaries.  _Fluids_ 5, 24 (2020).
  * [25] Lim, J. L. & DeMont, M. E.  Kinematics, hydrodynamics and force production of pleopods suggest jet-assisted walking in the american lobster (homarus americanus).  _Journal of Experimental Biology_ 212, 2731–2745 (2009).
  * [26] Zhang, C., Guy, R. D., Mulloney, B., Zhang, Q. & Lewis, T. J.  Neural mechanism of optimal limb coordination in crustacean swimming.  _Proceedings of the National Academy of Sciences_ 111, 13840–13845 (2014).
  * [27] Ford, M. P. & Santhanakrishnan, A.  On the role of phase lag in multi-appendage metachronal swimming of euphausiids.  _Bioinspiration & Biomimetics_ 16, 066007 (2021).
  * [28] Aubin, C. _et al._ Electrolytic vascular systems for energy-dense robots.  _Nature_ 571, 51–57 (2019).
  * [29] Katzschmann, R., DelPreto, J., MacCurdy, R. & Rus, D.  Exploration of underwater life with an acoustically controlled soft robotic fish.  _Science Robotics_ 3, eaar3449 (2018).
  * [30] Neveln, I. _et al._ Biomimetic and bio-inspired robotics in electric fish research.  _Journal of Experimental Biology_ 216, 2501–2514 (2013).
  * [31] Wang, Y., Wang, R., Wang, S., Tan, M. & Yu, J.  Underwater bioinspired propulsion: From inspection to manipulation.  _IEEE Transactions on Industrial Electronics_ 67, 7629–7638 (2019).
  * [32] White, C., Lauder, G. & Bart-Smith, H.  Tunabot flex: a tuna-inspired robot with body flexibility improves high-performance swimming.  _Bioinspiration & Biomimetics_ 16, 026019 (2021).
  * [33] Koehl, M.  Biomechanics of microscopic appendages: functional shifts caused by changes in speed.  _Journal of Biomechanics_ 37, 789–795 (2004).
  * [34] Cheer, A. & Koehl, M.  Paddles and rakes: fluid flow through bristled appendages of small organisms.  _Journal of Theoretical Biology_ 129, 17–39 (1987).
  * [35] Ford, M. P. & Santhanakrishnan, A.  Closer appendage spacing augments metachronal swimming speed by promoting tip vortex interactions.  _Integrative and Comparative Biology_ 61, 1608–1618 (2021).
  * [36] Erdman, A. G. & Sandor, G. N.  Mechanism design analysis and synthesis (vol. 1).  _Prentice-Hall_ (1997).
  * [37] Sandor, G. N. & Erdman, A. G.  Advanced mechanism design v. 2: Analysis and synthesis.  _Prentice-Hall_ (1984).
  * [38] Hedrick, T. L.  Software techniques for two-and three-dimensional kinematic measurements of biological and biomimetic systems.  _Bioinspiration & biomimetics_ 3, 034001 (2008).
  * [39] Gemmell, B. J., Jiang, H. & Buskey, E. J.  A new approach to micro-scale particle image velocimetry ($\mu$PIV) for quantifying flows around free-swimming zooplankton.  _Journal of Plankton Research_ 36, 1396–1401 (2014).

## Acknowledgements

We gratefully acknowledge Valentina Di Santo for insightful discussions that
strengthened the foundation of this work. We also thank Yair Sanchez Juarez
and Pedro Enrique Ávila Hernández for their work on the preliminary versions
of the Pleobot. Funding for this research was provided by the University of
California Institute for Mexico and the United States (UC MEXUS), the Consejo
Nacional de Ciencia y Tecnologia (CONACYT) through the UC MEXUS-CONACYT
collaborative grants program (UCMEXUS CN-18-138), Brown University OVPR
Research Seed Grant, and NASA Rhode Island EPSCoR Seed Grant (80NSSC22M0040).

## Author contributions statement

M.M.W., S.O.S, F.C.J, P.A.J.V. and O.M.L developed the concept and designed
experiments; S.O.S., N.T, Y.S. performed experiments and analyzed data;
M.M.W., S.O.S, F.C.J, P.A.J.V., O.M.L., N.T and Y.S. wrote the manuscript;
M.M.W. supervised the project. M.M.W. administered the project. All authors
reviewed the manuscript.

## Additional information

Competing interests The authors declare no competing interests.
Data Availability Statement Data will be provided upon reasonable request to
the corresponding author via email.

Figure 1: Morphology and kinematic parameters of the pleopod. Panels A
(_Euphausia superba_), C (_Palaemonetes paludosus_) and E show the kinematic
parameters of free-swimming shrimp incorporated in the Pleobot presented in
panel B, D and F: $\alpha$ is defined as the angle between the axis of the
body and the proximal segment (protopodite), $\beta$ is the angle between the
protopodite and the distal biramous segment (formed by the endopodite and the
exopodite), $\gamma$ appears during the power stroke as the exopodite and
endopodite separate, and $\zeta$ characterizes the cupping formed between the
exopodite and the endopodite. The Pleobot (panel G) is engineered based on a
mechanical gear train to actively control $\alpha$ and $\beta$, while
passively integrating $\gamma$. Note that panel C was captured from the back
of the organism, while panel E shows a close up view of the pleopods from the
side (as in panel A). Figure 2: Kinematic measurements. Validation of the
Pleobot was achieved by tracking its motion and plotting each kinematic
parameter alongside biological data [16]. Panels A-C show the evolution of
kinematic angles $\alpha$, $\beta$, and $\gamma$. Note that the standard
deviation over 5 cycles is indicated as shaded regions. Panels D and E show
the evolution of these angles from the lateral (panel D) and frontal (panel E)
views. Figure 3: Flow field measurements. The planes at which PIV data were
acquired are presented in panel A. Representative instantaneous velocity and
vorticity fields are shown in panels: B for the vertical plane I, C for the
frontal plane II, D for the horizontal plane III, and E for the vertical plane
IV. Flow field measurements acquired from a tethered live shrimp are presented
in panel F for reference. The black regions in panels B, and D-E are the
cross-sectional profiles of the Pleobot corresponding to the planes in A and
the white area in panel E represents a shadow area of the laser. In C, the
Pleobot is out-of-plane and its position is superimposed for clarity. Note
that the size and direction of the arrows scale with the value of the velocity
at a given grid point, while the colors scale with the value of the vorticity
according to the colormap. The red scale arrow represents 20 cm s-1. The scale
bar represents 1 cm in (B–E) and 1 mm in (F).  Figure 4: Force measurements.
The thrust (A) and lift (B) coefficients are averaged over five consecutive
cycles. The orange dashed lines show the kinematics of the appendage,
characterized by the angle $\Psi$. Corresponding pleopod profiles in panel C
for the anterior view of the $\gamma$ angle emphasize the change in pleopod
surface area over a beat cycle . Panels D-I show PIV measurements at different
times of a cycle, corresponding to the black dashed lines in panels A and B.
Panels D and E correspond to the lateral view along the endopodite during the
power stroke, showing a tip and bound vortex (D) that are shed at the end of
the power stroke (E). Panels F and G show a bottom view of the exopodite and
endopodite during the recovery stroke, with strengthening spanwise flow behind
the appendages between two side edge vortices. Panels H and I correspond to a
lateral view along the midplane of the exopodite during the power stroke,
where we identify a leading-edge vortex (LEV) that initially contributes
primarily to lift and then thrust as the AoA changes. The scale bar and the
red scale arrow represent 1 cm and 20 cm s-1, respectively for all frames.
Figure 5: Scaled KRill Inspired Modular Pleobots (SKRIMP). This sketch
illustrates how the modularity and rapid prototyping of the _Pleobot_ enables
the investigation of metachronal swimming by employing several units. This
includes fluid flows and forces - in several taxa, body plans, and appendage
configurations. _SKRIMP_ constitutes the baseline to engineer a new generation
of AUVs operating in complex marine environments leveraging the swimming
characteristics of metachronal swimmers.


# Learning Feynman Diagrams using Graph Neural Networks

Harrison Mitchell
Cavendish Laboratory
University of Cambridge
<EMAIL_ADDRESS>
&Alexander Norcliffe
Computer Laboratory
University of Cambridge
<EMAIL_ADDRESS>&Pietro Liò
Computer Laboratory
University of Cambridge
<EMAIL_ADDRESS>
###### Abstract

In the wake of the growing popularity of machine learning in particle physics,
this work finds a new application of geometric deep learning on Feynman
diagrams to make accurate and fast matrix element predictions with the
potential to be used in analysis of quantum field theory. This research uses
the graph attention layer which makes matrix element predictions to 1
significant figure accuracy above 90% of the time. Peak performance was
achieved in making predictions to 3 significant figure accuracy over 10% of
the time with less than 200 epochs of training, serving as a proof of concept
on which future works can build upon for better performance. Finally, a
procedure is suggested, to use the network to make advancements in quantum
field theory by constructing Feynman diagrams with effective particles that
represent non-perturbative calculations.

##  1 Introduction and Related Work

Particle physics has recently seen a rise in the popularity of machine
learning [1], [2], [3]. Applications include, analysing particle jets [4] [5],
tracking particle paths [6] and detecting collider events [7]. This work
explores a new application of geometric deep learning [8][9] on graph neural
networks (GNNs) by predicting matrix elements from Feynman diagrams for simple
particle interactions. Analytically, matrix elements are computed from Feynman
diagrams, a pictorial representation of a particle interaction, using a set of
rules called the Feynman rules, which translate the diagram into a
mathematical expression. The matrix elements then find wide use in computing
probabilities in particle physics. Graph neural networks are a type of network
that operate on graph valued data, to make node, edge or graph level
predictions. In this work, a graph level prediction is used for matrix
elements from a graph representation of Feynman diagrams. The convolutional
layer used is the graph attention layer (GAT) [10], which previously has shown
application in text classification [11], material property predictions [12]
and classical physics [13]. This work also aims to extend the testing of the
performance of the GAT layer and the effectiveness of its attention head
mechanism on Feynman diagrams. The advantage of using machine learning for
physics is partly in faster computation time for complicated matrix elements,
but the envisioned future of this study is that the neural network learns all
Feynman rules so that it can be connected to a dynamic graph network (for
example, the temporal graph network [14]). Such a network would progressively
construct new nodes and edges, creating new Feynman diagrams which would be
passed to the learned GNN whose output would be used to fit to matrix elements
from experimental data. This may predict effective Feynman diagrams in non-
perturbative regimes, such as low energy QCD, where truncating the infinite
series of perturbations does not approximate the matrix elements. Many of the
matrix elements that were used extensively in the learning process for the
neural networks were adapted from the 1995 work by Borodulin [15].

##  2 Method

(a) An example tree-level graph, that the network can take as an input. The
work of this paper is on diagrams with this structure but with different edge
features. Blue vectors are node features - interactions. Red vectors are edge
features - particles.

(b) Top left image depicts the action of the GAT layer. Top right is a
portrayal of the output of the GAT layers being pooled into a graph
representation and momenta added. This is passed to an FCN to get the
predictions.

Figure 1: Sample input (left) and network architecture (right).

###  2.1 Dataset encoding

Feynman diagrams have a natural graph representation as seen in Figure 1(a),
with interaction vertices being represented by nodes, and particles being
represented by edges. The edge features contain all the information to
describe a particle, whilst the targets are the matrix elements, averaged over
the helicity combinations. In our graph encodings, the momenta are not
included in the edge features and are instead concatenated after the GNN has
produced a graph representation. The concatenation process is outlined in
Figure 1(b). Since the Feynman rules are independent of momentum in both QED
and QCD, the representation of each graph should be as well. This method
provides scalability by removing the need for multiple diagrams for each
possible helicity combination. This is a more efficient procedure when the
number of Feynman diagrams is large. The particle encoding is given in
Equation 1, where masses are on-shell and in mega electronvolts (MeV).

$\text{edge
vector}=[m,S,I_{L}^{3},Y_{L},I_{R}^{3},Y_{R},r,g,b,\bar{r},\bar{g},\bar{b}]$
(1)

All data is given in the centre of mass frame with azimuthal angle $\phi=0$
and particles aligned with the $z$-axis, so that only the polar angle,
$\theta$, and the momentum, $p$, are needed to categorise the interaction.
Each data point includes a graph, representing the Feynman diagram, and a
$(p,\theta)$ pair. A large dataset consisting of 80,000 such data points with
ultra-relativistic momenta in the range $[1,1000]$GeV, was then constructed
and training was carried out on subsets containing only the appropriate
Feynman diagrams for the experiment. Each subset consisted of 2000 data
points, sampled evenly over the range of target matrix elements. The network
is fed three random splits of the dataset: training, validation and test
datasets The split ratio was 6:2:2 respectively. The training dataset is what
the model sees and uses to conduct back-propagation and the order is shuffled
each epoch. The validation dataset is passed through the model at the end of
each epoch. The test dataset is passed through after training.

A global node was also added to improve the training speed and performance of
the network. The three interaction strengths for QED, QCD and the weak force
($\alpha_{QED},\alpha_{S},\alpha_{W}$) were stored in the global node vector,
each given at zero energy. The logarithmic running of these constants can be
learnt during propagation through subsequent node embeddings of the global
node in the forward pass. This addition reduced the number of required layers
and trainable parameters in the GNN since node information on opposite ends of
the graph only needs two layers to pass to each other via the global node. It
was seen that even when the number of attention heads and layers in the GNN
were reduced, the network performed similarly to before adding the global
node.

###  2.2 Network Architecture

The aim is that the GNN will create a representation that encodes information
for the Feynman rules. The network consists of a GNN connected to a fully
connected network (FCN). The GNN architecture is built using three key
components: a convolutional layer, a linear layer and a pooling layer. The
convolutional layer used is the graph attention layer [10]. The linear layer
invokes simple fully connected layer with a leaky ReLu activation function,
needed to condense the size of the graph, since the attention mechanism
generates a larger graph depending on the number of attention heads. The
pooling layer uses global sum pooling in order to capture information across
all nodes as well as distinguishing different sized graphs by the magnitude of
the values in the representation. The graph representation is then passed to
the FCN consisting of two hidden layers and a leaky ReLU activation function.
The final output layer has just one neuron whose value is passed through a
sigmoid activation function so that the output is squeezed between zero and
one. The architecture is summarised in the schematic in Figure 1(b). Further
implementation details are included in Appendix B.

##  3 Results and Discussion

The network was trained on high energy tree-level QED diagrams involving
electrons and muons on which the model accurately predicts matrix elements for
several processes. There were 2 main metrics used to evaluate performance:
$L^{1}$ loss and the accuracy. Accuracy is defined as the proportion of
results in the dataset that are the same as the target after rounding up to a
number of decimal places (d.p.), $a$, as shown in Equation 2 .

$\mathcal{A}(\\{y_{i}\\},\\{\hat{y}_{i}\\},a)=\frac{1}{N}\sum_{i}^{N}\left(1-\Theta\left[L^{1}\left(\lceil
y_{i}\rceil_{a},\lceil\hat{y}_{i}\rceil_{a}\right)\right]\right)$ (2)

Where $\lceil x\rceil_{a}$ means round up to the $a^{th}$ decimal place and
$\Theta(x)$ is the step function. The predicted matrix elements are plotted
against the polar angle in the range $[0,2\pi]$. At large momenta, the matrix
element is independent of momentum, and so the variation is in the angle.
These metrics are tabulated in each experiment’s data tables. A further
extension of the model is tested on QCD, involving quark diagrams. The results
are included in Appendix C.

###  3.1 Three diagram dataset

(a) Results for $e^{-}e^{+}\to\mu^{-}\mu^{+}$ from a GNN that has seen both
matrix elements

(b) Results for $e^{-}e^{+}\to e^{-}e^{+}$ Bhabha scattering matrix elements.

Figure 2: Results for training on three diagrams

The initial model was tested on two graphs consisting of three Feynman
diagrams: the electron-positron to muon-antimuon s-channel diagram; and the
two Bhabha electron-positron to electron-positron s and t-channel diagrams,
batched into one graph with ultra-relativistic matrix element. Since both
graphs are combined into a single dataset the GNN is forced to distinguish
between the two different inputs, testing its ability to construct a
sufficiently detailed graph representation to feed to the FCN to communicate
which matrix element is to be fitted.

Figures 2(a) and 2(b) exhibit the model’s performance in fitting over the
range of angles. The predictions are good with a small discrepancy for the
muon matrix element at large angles. More detail on the loss and accuracy
values are given in Table 1, where it can be seen that the accuracy to one
decimal place is high but just shy of 100% because the model predicts some
matrix elements very close to the target, but some much farther away so that
the $L^{1}$ loss is consistently low but allowing for a much tighter fit for
some of the data points reaching two and three decimal place accuracy over 81%
and 11% of the time, respectively.

Table 1: Loss and accuracy values for the three combined Feynman diagrams for $e^{-}e^{+}\to\mu^{-}\mu^{+}$ and $e^{-}e^{+}\to e^{-}e^{+}$ Test metric | Data
---|---
Test accuracy (1d.p.) | 99.00%
Test accuracy (2d.p.) | 81.50%
Test accuracy (3d.p.) | 11.50%
Test L1 Loss | 0.0049

###  3.2 Tree-level QED dataset for electrons and muons

(a) Results from the fully trained network on all tree-level QED with
electrons and muons, for the $e^{-}e^{+}\to\mu^{-}\mu^{+}$ matrix element.

(b) Results from the fully trained network, for the $e^{-}e^{+}\to e^{-}e^{+}$
matrix element.

Figure 3: Training on full dataset

The dataset was expanded to encompass six graphs with four high energy matrix
elements at tree-level QED, for each of the permutations of electrons and
muons in the initial and final states and their allowed s and t-channels.
There is a drop in performance but the overall target shape is maintained.
Figure 3(a) shows the predicted and theoretical matrix elements for the same
s-channel $e^{-}e^{+}\to\mu^{-}\mu^{+}$ interaction. This time, this is one of
many graphs that the network is able to predict, but despite the increase in
number of points to fit, we still see strong performance. For this interaction
the discrepancy at the extreme angles 0 and $2\pi$ is more noticeable, but
there is additional error at the minimum $\theta=\pi$. Since the model
performs best in the regions that are close to linear, this may be indicative
of a requirement for additional non-linearities.

Figure 3(b) shows the predicted and theoretical matrix element for the s and
t-channel $e^{-}e^{+}\to e^{-}e^{+}$ interactions, where we see much stronger
agreement. Increasing the number of dense neurons showed signs of improvement,
but significantly increased the training time, but since this target consists
of two largely linear regimes, the models performs well considering the
training time. Table 2 shows the metrics for the learned GNN. The most
noticeable difference is the much higher loss. This behaviour is expected for
a larger dataset as there are more points for the network to fit to, but
despite increasing the width of the FCN, the GNN was unable to fit as
accurately as before. This may have been, in part, due to a need for an
increased number of neurons in the FCN to be able to decode and fit multiple
different functions depending on the graph representation it sees.
Alternatively it may be that the GNN’s graph representations were not large
enough to discriminate between the graphs. The graph representation was kept
to a constant dimension in the hope that the strong fit to the smaller three
diagram dataset implied that the Feynman rules were sufficiently encoded.
Increasing this may help to encode further detail of the Feynman rules with
different starting and final states. Finally, the minima search may be further
improved by smoothing the landscape with stochastic weight averaging by
allowing the network to find lower minima and therefore reduce the loss.

However, the ability to identify different graphs and reconstruct different
functions suggests that the GNN has the capacity to encode the Feynman rules
into its graph representation. With some improvement in implementation it is
possible that the GNN can make more accurate predictions and for a wider range
of particle states.

Table 2: Loss and accuracy values for all tree-level QED Feynman diagrams involving all permutations of $e^{-}e^{+}$ and $\mu^{-}\mu^{+}$ pairs Test metric | Data
---|---
Test accuracy (1d.p.) | 96.75%
Test accuracy (2d.p.) | 44.00%
Test accuracy (3d.p.) | 0.75%
Test L1 Loss | 0.015

##  4 Conclusion

Graph neural networks have been shown to be able to make basic matrix element
predictions by fitting to analytical data. This provides improved computation
time for matrix element calculations and scales well when introducing new
matrix elements. This work suggests that it is possible for the GNN’s graph
representation to encode the elements of the Feynman rules, and for the FCN to
decode it into a matrix element. Following on from the research in the field
of geometric deep learning, graph attention networks have been shown to be
effective in regression problems for Feynman diagrams. This supports the use
of geometric deep learning in particle physics beyond experimental analysis.
The aspiration is that this work will be extended to have learnt a full
dataset with all fundamental forces that, in itself, would be useful in making
experimental predictions. To further build on this, if the dataset is extended
to loop diagrams, the network will learn its own procedure to deal with
renormalization. This could then be connected to a network for dynamic graphs,
which will learn the edges and nodes to build, which upon passing to the GNN,
correctly fit to experimental data. The constructed diagrams are then ideally
interpretable as Feynman diagrams, corresponding to physical processes, and if
not, then at least it may provide effective Feynman diagrams, whose matrix
elements represent an infinite series in non-perturbative QFT.

## Acknowledgments and Disclosure of Funding

We thank Bentley Carr for his time and insight during the research of this
project. Additionally we would like to take the chance to thank the reviewers
for raising challenging points that drove us to find new perspectives and
better the future direction of this work.

Finally, we also thank the Cavendish Laboratory and Computer Laboratory at the
University of Cambridge for facilitating the original research.

## References

  * [1] Matthew Feickert and Benjamin Nachman  “A Living Review of Machine Learning for Particle Physics”  arXiv, 2021  DOI: 10.48550/ARXIV.2102.02770
  * [2] Dan Guest, Kyle Cranmer and Daniel Whiteson  “Deep Learning and Its Application to LHC Physics”  In _Annual Review of Nuclear and Particle Science_ 68.1, 2018, pp. 161–181  DOI: 10.1146/annurev-nucl-101917-021019
  * [3] Jonathan Shlomi, Peter Battaglia and Jean-Roch Vlimant  “Graph neural networks in particle physics”  In _Machine Learning: Science and Technology_ 2.2 IOP Publishing, 2021, pp. 021001  DOI: 10.1088/2632-2153/abbf9a
  * [4] James Barnard, Edmund Noel Dawe, Matthew J. Dolan and Nina Rajcic  “Parton shower uncertainties in jet substructure analyses with deep neural networks”  In _Physical Review D_ 95.1 American Physical Society (APS), 2017  DOI: 10.1103/physrevd.95.014018
  * [5] Luke Oliveira et al.  “Jet-images — deep learning edition”  In _Journal of High Energy Physics_ 2016.7 Springer ScienceBusiness Media LLC, 2016  DOI: 10.1007/jhep07(2016)069
  * [6] Javier Duarte and Jean-Roch Vlimant  “Graph Neural Networks for Particle Tracking and Reconstruction”  In _Artificial Intelligence for High Energy Physics_ WORLD SCIENTIFIC, 2022, pp. 387–436  DOI: 10.1142/9789811234033_0012
  * [7] Yasir Alanazi et al.  “A Survey of Machine Learning-Based Physics Event Generation”  In _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence_ International Joint Conferences on Artificial Intelligence Organization, 2021  DOI: 10.24963/ijcai.2021/588
  * [8] Franco Scarselli et al.  “The Graph Neural Network Model”  In _IEEE Transactions on Neural Networks_ 20.1, 2009, pp. 61–80  DOI: 10.1109/TNN.2008.2005605
  * [9] Michael Bronstein et al.  “Geometric Deep Learning: Going beyond Euclidean data”  In _IEEE Signal Processing Magazine_ 34.4 Institute of ElectricalElectronics Engineers (IEEE), 2017, pp. 18–42  DOI: 10.1109/msp.2017.2693418
  * [10] Petar Velickovic et al.  “Graph Attention Networks”  arXiv, 2017  DOI: 10.48550/ARXIV.1710.10903
  * [11] Liang Yao, Chengsheng Mao and Yuan Luo  “Graph Convolutional Networks for Text Classification”  In _Proceedings of the AAAI Conference on Artificial Intelligence_ 33(01), 2019, pp. 7370–7377  URL: https://doi.org/10.1609/aaai.v33i01.33017370
  * [12] Steph-Yves Louis et al.  “Graph convolutional neural networks with global attention for improved materials property prediction”  In _Phys. Chem. Chem. Phys._ 22, 2020, pp. 18141–18148  DOI: 10.1039/D0CP01474E
  * [13] Rui Wang et al.  “Bridging Physics-based and Data-driven modeling for Learning Dynamical Systems”  In _Proceedings of the 3rd Conference on Learning for Dynamics and Control_ 144, Proceedings of Machine Learning Research  PMLR, 2021, pp. 385–398  URL: https://proceedings.mlr.press/v144/wang21a.html
  * [14] Emanuele Rossi et al.  “Temporal Graph Networks for Deep Learning on Dynamic Graphs”  arXiv, 2020  DOI: 10.48550/ARXIV.2006.10637
  * [15] V.. Borodulin, R.. Rogalyov and S.. Slabospitsky  “CORE 2.1 (COmpendium of RElations, Version 2.1)”  arXiv, 1995  DOI: 10.48550/ARXIV.HEP-PH/9507456
  * [16] Adam Paszke et al.  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  In _Advances in Neural Information Processing Systems 32_ Curran Associates, Inc., 2019, pp. 8024–8035  URL: http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf
  * [17] Matthias Fey and Jan Eric Lenssen  “Fast Graph Representation Learning with PyTorch Geometric”, https://github.com/pyg-team/pytorch_geometric, 2019
  * [18] William Falcon and The PyTorch Lightning team  “PyTorch Lightning”, https://github.com/PyTorchLightning/pytorch-lightning, 2019  DOI: 10.5281/zenodo.3828935
  * [19] DeepFindr  “A Graph Neural Network project on HIV data”  In _GitHub repository_ GitHub, https://github.com/deepfindr/gnn-project, 2021
  * [20] Philippe  “uvadlc notebooks”  In _GitHub repository_ GitHub, https://github.com/phlippe/uvadlc_notebooks/blob/master/docs/tutorial_notebooks/tutorial7/GNN_overview.ipynb, 2022
  * [21] Joan Bruna, Wojciech Zaremba, Arthur Szlam and Yann LeCun  “Spectral Networks and Locally Connected Networks on Graphs”  arXiv, 2013  DOI: 10.48550/ARXIV.1312.6203
  * [22] R. Mertig, M. Böhm and A. Denner  “Feyn Calc - Computer-algebraic calculation of Feynman amplitudes”  In _Computer Physics Communications_ 64.3, 1991, pp. 345–359  DOI: https://doi.org/10.1016/0010-4655(91)90130-D
  * [23] Vladyslav Shtabovenko, Rolf Mertig and Frederik Orellana  “New developments in FeynCalc 9.0”  In _Computer Physics Communications_ 207 Elsevier BV, 2016, pp. 432–444  DOI: 10.1016/j.cpc.2016.06.008
  * [24] Vladyslav Shtabovenko, Rolf Mertig and Frederik Orellana  “FeynCalc 9.3: New features and improvements”  In _Computer Physics Communications_ 256 Elsevier BV, 2020, pp. 107478  DOI: 10.1016/j.cpc.2020.107478
  * [25] William L. Hamilton, Rex Ying and Jure Leskovec  “Inductive Representation Learning on Large Graphs”  arXiv, 2017  DOI: 10.48550/ARXIV.1706.02216
  * [26] William L. Hamilton  “Graph Representation Learning”  In _Synthesis Lectures on Artificial Intelligence and Machine Learning_ 14.3 MorganClaypool, pp. 1–159  URL: https://www.cs.mcgill.ca/~wlh/grl_book/
  * [27] Justin Gilmer et al.  “Neural Message Passing for Quantum Chemistry”  arXiv, 2017  DOI: 10.48550/ARXIV.1704.01212
  * [28] Hanjun Dai, Bo Dai and Le Song  “Discriminative Embeddings of Latent Variable Models for Structured Data”  arXiv, 2016  DOI: 10.48550/ARXIV.1603.05629
  * [29] Anja Butter and Tilman Plehn  “Generative Networks for LHC events”  arXiv, 2020  DOI: 10.48550/ARXIV.2008.08558
  * [30] Pavel Izmailov et al.  “Averaging Weights Leads to Wider Optima and Better Generalization”  arXiv, 2018  DOI: 10.48550/ARXIV.1803.05407
  * [31] Yoshua Bengio  “Practical recommendations for gradient-based training of deep architectures”  arXiv, 2012  DOI: 10.48550/ARXIV.1206.5533

##  Appendix A Societal impacts

This work has had two fields of development: one in the use of machine
learning to further physics and the other to use physics to research machine
learning. The development of physical theory has limited short term societal
impact, since particle physics is many years from finding implementation in
technology and engineering. This paper contributes towards the effectiveness
of the graph attention layer but only in its application to Feynman diagrams
and has its limitations. It does not suggest it could be easily substituted
for unethical datasets and does not promote its use in immoral situations. The
datasets used were self-constructed and did not use any sensitive information.
Creating a trained network is also fast and not energy intensive with Google
Colab’s free GPUs being sufficient with 12.68GB RAM used to train the whole
dataset, so there is little environmental impact.

##  Appendix B Implementation

###  B.1 Dataset encoding

Figure 4 gives an example Feynman diagram and a graphical depiction of its
encoding to a graph.

Figure 4: Encoding a Feynman diagram as a graph. The blue vectors are the node
features and the red are the edge features.

The nodes were also encoded as to make the graph orientable by capturing
information about the time-like direction of the interaction. This was
achieved by specifying which nodes are initial state particles, final state
particles and virtual particles, using a one-hot encoding.

$\text{node vector}=[\text{initial},\text{virtual},\text{final}]$ (3)

###  B.2 Training details

The network was constructed in PyTorch [16], with the aid of PyTorch Geometric
[17] and PyTorch Lightning [18]. We used the help of several GitHub
repositories including the work by DeepFindr [19] and the tutorial from
Pytorch Geometric [20]. All of the code was written in Google Colab using
their free GPUs, with about 13GB of RAM total.

The loss was generated with a variety of loss functions including mean square
error (MSE), absolute error ($L^{1}$), and LogCosh with the loss of a batch
being calculated by the sum of the individual losses. Ultimately all losses
performed similarly, but $L^{1}$ loss was primarily used in this study as
slight improvement on accuracy was found. A stochastic gradient descent
optimizer with momentum 0.4 was used. Damping was also tested but did not
improve the performance. The learning rate was started high at 1 but a
scheduler that reduced the learning rate when the loss plateaued was
implemented to hone in on the minima. The mini-batch size used depended on the
dataset but was kept at powers of 2, between 16 and 64, to maximize storage
use. This was a good compromise of speed and performance. The graph
representation vector had dimension that depended on the complexity of the
dataset, but reached a size of at most 8. This was sufficient for the network
to encode the relevant information to pass to the FCN. The activation
functions used throughout are ReLU and Leaky ReLU, these were fast to train
and did not give stagnation problems that sigmoid and tanh do.

The number of dense neurons in the FCN was progressively increased if the
model under-fitted, but a value between 64 and 128 was enough to be able to
fit to the functions for this paper. The depth of the GNN was kept to only 2
or 3 layers, due to the global node.

###  B.3 Runtime

The GNN training was completed on the order of minutes. Figure 5(a) shows how
the FCN begins to fit to the targets after just one epoch. Figure 5(b) shows
how the training time increases as the complexity of the dataset grows. The
blue line represents training on just one graph, which is why it converges on
a minimum so quickly. The orange line represents the most complex dataset, and
takes longer than the pink line to reach the same level of loss, being
overtaken at around 5000 steps. The slow start is due to the stochastic nature
of the optimization algorithm. Once it finds a minimum, the pink line
converges more quickly.

Training time is largely due to the number of parameters to fit. As the
dataset gets more complex, more parameters are required and so the parameter
space to search is larger. With the large learning rates used, it is possible
that steep local minima are also missed.

The sudden jumps in each line is due to the learning rate scheduler. As the
model starts to oscillate around the minima, the speed is dropped so that the
model does not overshoot.

(a) Training on $e^{-}e^{+}\to e^{-}e^{+}$ Bhabha scattering diagrams after
only one epoch of training

(b) Graph showing the improvement of the network loss at the end of each epoch
on the validation dataset. Blue: one graph. Pink:
$e^{-}e^{+}\to\mu^{-}\mu^{+}$ and $e^{-}e^{+}\to e^{-}e^{+}$ graphs. Orange:
whole dataset.

Figure 5: Graphs showing the speed to train the network

###  B.4 Program code

https://github.com/Clearbloo/Feynman_GNN.git

##  Appendix C QCD

QCD was also trialed to test the models generalization to a different set of
Feynman rules. The test process used was $u\bar{u}\to t\bar{t}$ via a gluon
s-channel. Since top quarks are massive, and not in chirality eigenstates,
this also doubled as a test for the network to distinguish momentum
dependence. There was a significant drop in performance, as shown in Figure 6
and Table 3. However, the shape of the original QED predictions was retained
and the model was able to identify the different curves based on momenta. This
suggests the GNN can be scaled up to include more matrix elements and
different sections of the standard model interaction Hamiltonian, with an
increasingly large graph representation and improvements on implementation.
This is expected as different Feynman rules need additional encodings on the
graph representation.

Table 3: Loss and accuracy values for $u\bar{u}\to t\bar{t}$ via gluon exchange Test metric | Data
---|---
Test accuracy (1dp) | 94.0%
Test accuracy (2dp) | 42.0%
Test accuracy (3dp) | 3.75%
Test L1 Loss | 0.0265
Figure 6: QCD results