aainstitutetext: Konya Technical University,
Road 42250, Selcuklu/Konya, Turkeybbinstitutetext: Fermilab,
P.O. Box 500, Batavia, IL, 60510, U.S.A.

#  Mass distributions of dijet resonances from excited quarks at proton-proton
colliders

Emine Gurpinar Guler a Yalcin Guler b and Robert M. Harris
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study the expected experimental mass distributions of dijet resonances from
excited quarks in proton-proton collisions at energies $\sqrt{s}=$ 13, 14, 27,
100, 300, and 500 TeV. We explore in detail the expected shapes at both the
generator and experimental levels, and identify within the distributions the
effects of the excited quark natural width, parton momentum distributions of
the proton, radiation, and experimental resolution. We present both
differential and cumulative probability distributions as a function of dijet
mass, and the signal acceptance of a window in dijet mass centered on each
resonance. We find that for a range of resonance masses, between 10% and 50%
of $\sqrt{s}$, the dijet mass distributions and window acceptance are
practically universal, approximately invariant under changes in resonance mass
and $\sqrt{s}$. This work supports our Snowmass 2021 study on the sensitivity
to dijet resonances at proton-proton colliders.

FERMILAB-PUB-22-301-PPD-SCD

May 5, 2022

##  1 Introduction

A significant benchmark for discovery at a proton-proton ($pp$) collider is
the sensitivity to a dijet resonance Harris:2022kls , $X$, the intermediate
state of the s-channel process $pp\rightarrow X\rightarrow 2\mbox{ jets}$. To
search for new particles with the highest values of resonance mass ($M$),
hadron collider experiments have used the classic technique of searching for
bumps in the invariant mass spectrum of dijets Harris:2011bh , the two
individually resolved jets with the largest transverse momentum. The most
common benchmark model for dijet resonances is the excited quark, $q^{*}$
ref_qstar ; Baur:1989kv , explicitly sought and constrained by every LHC
search for resolved dijet resonances Aad:2011aj ; Khachatryan:2010jd ;
ATLAS2010 ; Chatrchyan2011123 ; ATLAS:2012pu ; CMS:2012yf ; CMS:2013qha ;
ATLAS:2014aqa ; CMS:2015sja ; CMS:2015dcf ; ATLAS:2015nsi ; CMS:2016ecr ;
CMS:2016iap ; ATLAS:2017yvp ; ATLAS:2018fzt ; CMS:2018xlo ; ATLAS:2019hjw ;
CMS:2019vgj . We expect the $q^{*}$ model will also be a benchmark of dijet
resonance searches at future $pp$ colliders.

In this paper we explore the expected shape at future $pp$ colliders of the
distribution in dijet mass ($m$) of the $q^{*}$ signal. This shape determines
the fraction of signal events that form a characteristic resonant peak, a
bump, that experiments search for on top of the large continuum background of
dijet events from quantum chromodynamics (QCD). The signal shape as a function
of dijet mass is used by the experiments to determine statistical significance
of potential signals, determine upper limits on the production cross section,
and set mass limits on the $q^{*}$ model. To provide a foundation to better
understand the experimental sensitivity of $pp$ colliders to excited quarks,
this paper will explore the generator and experiment level contributions to
the shapes of expected dijet mass distributions from excited quarks. We will
establish the approximate invariance of the resonance shapes as a function of
both $M$ and $\sqrt{s}$, and discuss the cases where that invariance is
violated. We expect this exploration and discussion will be of interest to
searches for dijet resonances at both current and future $pp$ colliders. In a
previous paper Harris:2022kls  we estimated the sensitivity of future $pp$
colliders to dijet resonances, including $q^{*}$. This study will use the
shapes of excited quarks at future $pp$ colliders to support that paper,
demonstrating the range of validity of the dijet mass window acceptance we
used.

We study existing, planned, proposed and envisioned $pp$ colliders at six
values of $\sqrt{s}$. The energy $\sqrt{s}=13$ TeV, used to collect data from
2015 to 2018 at the Large Hadron Collider (LHC) at CERN, connects our studies
with published searches to provide deeper understanding of the signal shapes
that were used. There is also a strong interest in $\sqrt{s}=14$ TeV, the
planned energy of the LHC in its final high luminosity operation (HL-LHC)
anticipated to begin in 2029 and continue for about a decade111See the LHC
schedule at https://lhc-commissioning.web.cern.ch/schedule/LHC-long-term.htm.
We present just a few results for $\sqrt{s}=27$ TeV, a previously studied high
energy option for the LHC (HE-LHC) FCC:2018bvk , an energy which also has been
considered for a $pp$ collider to fill the Fermilab site Bhat:2022hdi . Our
studies are more complete for $\sqrt{s}=100$ TeV, because it is a large next
step in energy to a proposed and widely studied $pp$ collider stage of a
Future Circular Collider (FCC-hh) at CERN FCC:2018vvp . One result is
presented for $\sqrt{s}=300$ TeV, because it was a high energy option
discussed at Snowmass 2021 and included in our previous paper Harris:2022kls .
Finally, to understand how the excited quark shapes evolve to the highest
possible collision energies, we include complete studies for $\sqrt{s}=500$
TeV, the envisioned energy of a Collider in the Sea McIntyre:2017ibd .

##  2 Generator samples

We have produced generator samples for the sub-process $qg\rightarrow
q^{*}\rightarrow qg$ in $pp$ collisions using PYTHIA8 Version 8.243
Sjostrand:2014zea . Each sample contains $20$ thousand events with $q^{*}$
mass $M$ equal to 10%, 20%, 30%, 40%, 50% and 60% of the considered $pp$
collision energy, $\sqrt{s}=13$, $14$, $27$, $100$, $300$ and $500$ TeV. The
PYTHIA8 settings for all the samples were identical to those used by
previously published CMS simulations CMS:2018xlo  for $\sqrt{s}=13$ TeV, and
to directly compare with those simulations, we also produced samples at $M=2$,
$4$, $6$ and $8$ TeV for $\sqrt{s}=13$ TeV. As in Ref. CMS:2018xlo , we used
the CTEQ6L1 parton distribution function (PDF) set refCTEQ , the CUETP8M1
event tune Khachatryan:2015pea ; Skands:2014pea , and the standard assumptions
about excited quarks: mass degenerate $u^{*}$ and $d^{*}$ only, compositeness
scale $\Lambda=M$, and the same couplings of excited quarks to quarks and
gauge bosons as in the standard model ($f=f^{\prime}=f_{s}=1$).

##  3 Dijet types

In this paper we consider dijet mass distributions for four different types of
dijets, representing different steps in the production, reconstruction and
potential observation of the excited quark.

First, we construct "di-parton" dijet mass distributions, the invariant mass
of the quark-gluon state within the sub-process $qg\rightarrow
q^{*}\rightarrow qg$. Di-parton mass is equivalent to the sub-process total
energy in the center-of momentum frame, often denoted as $\sqrt{\hat{s}}$.

Second, we construct "AK4 genjet" dijet mass distributions, the invariant mass
of the two AK4 genjets with highest transverse momentum ($p_{T}$) in the
event, also called the two leading AK4 genjets. Here, AK4 genjets are found by
clustering all the particles within the event using the anti-$k_{T}$ algorithm
Cacciari:2008gp ; Cacciari:2011ma  with a distance parameter $R=0.4$, the
effective radius of a narrow jet cone. This AK4 algorithm is a jet
reconstruction algorithm, chosen as the LHC standard because it is useful for
the widest range of analyses that make use of jets, but the narrow jet cone
size makes it suboptimal for dijet resonance event reconstruction.

Third, we construct "wide genjet" dijet mass distributions, where we implement
the wide jet algorithm discussed in Ref. CMS:2018xlo  using AK4 genjets as
input. In this algorithm the two leading AK4 genjets are used as seeds for
each of the two wide genjets, and every event is required to have exactly two
wide genjets. Each wide genjet contains exactly one leading AK4 genjet, as
well as all other AK4 genjets with pseudorapidity $|\eta|<2.5$ and
$p_{T}>(\sqrt{s}/13\ \mbox{TeV})\times 30$ GeV found within a distance $\Delta
R=\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}<1.1$ from the leading AK4 genjet.
The genjet $p_{T}$ requirement above is generalized for all pp collision
energies from the value in Ref. CMS:2018xlo , scaled with $\sqrt{s}$ from the
requirement $p_{T}>30$ GeV at $\sqrt{s}=13$ TeV. The wide jet algorithm
effectively widens the cone size to $R=1.1$, recovering the energy lost to
jets radiated outside that distance, without including low energy jets from
pileup. It is an event algorithm, forcing the dijet interpretation on each
event, and it has been optimized for dijet resonance event reconstruction
CMS:2018xlo .

Fourth, we construct "smeared wide genjet" dijet mass distributions, by
convoluting the dijet mass distribution from wide genjets with the anticipated
experimental dijet mass resolution of wide jets estimated in section 5. This
final set of dijet mass distributions are the ones with all estimated
experimental effects.

##  4 Generator-level resonances

Generator-level dijet mass distributions predict the shape of dijet resonances
before the effects of experimental resolution. These distributions are
valuable for understanding the physics of the shape of dijet resonances at
future $pp$ colliders. They are determined by the properties of the underlying
model and the parton shower within the generator, which are relatively well
known compared to the significant uncertainties in the experimental resolution
of future detectors. Generator level dijet mass distributions for excited
quarks are presented in Fig. 1, 2 and 3, and we discuss in this section what
we learn from them about the physics determining the shape of dijet
resonances.

Figure 1:  Dijet mass distributions for excited quarks with a resonance mass
equal to 50% of $\sqrt{s}$, from $pp$ collisions at $\sqrt{s}$ equal to 13,
100, and 500 TeV (top row) and 13, 14, and 27 TeV (bottom row), for the cases
where the dijet mass is calculated from the Lorentz vectors of the two final-
state partons (left column), the two wide genjets (middle column), and the two
leading AK4 genjets in the event (right column).

###  4.1 Collision energy invariance

Figure 1 demonstrates the approximate invariance of the dijet mass
distribution to the choice of $\sqrt{s}$. For a fixed value of the $q^{*}$
mass ratio, $M/\sqrt{s}=0.5$, Fig. 1 shows the dijet mass distributions of di-
partons, wide genjets and AK4 genjets. This choice of $q^{*}$ mass ratio was
shown to be a critical value in Ref. Harris:2022kls , the approximate mass
limit for excited quarks at $pp$ colliders at the baseline integrated
luminosity. All of the distributions in Fig. 1 are for the dimensionless
ratio, $m/\sqrt{s}$, which allows a fair comparison of the shapes for various
choices of collider $\sqrt{s}$. Variable width bins are used corresponding to
1% of the dijet mass and starting at a value $m/\sqrt{s}=0.05$: the bin edges
are $0.05$, $0.0505$, $0.051005$, etc.

Figure 1 (left) shows the di-parton mass distributions from the generator,
which are almost completely invariant with changes in $\sqrt{s}$. They
originate from the differential cross section which has the general form:

$\frac{d\sigma}{dm}\sim\int_{0}^{1}\int_{0}^{1}{dx_{1}dx_{2}(q(x_{1})g(x_{2})+q(x_{2})g(x_{1}))\left[\frac{\Gamma^{2}}{\left(m^{2}-M^{2}\right)^{2}+M^{2}\Gamma^{2}}\right]\delta(x_{1}x_{2}-m^{2}/s)}$
(1)

where $q(x)$ and $g(x)$ are the PDFs of the initial state quark and gluon,
respectively, and $x_{1}$ and $x_{2}$ are their fractional momenta within each
of the two colliding protons. The term in square brackets in Eq. 1 is the
relativistic Breit-Wigner distribution of the resonant sub-process with width
$\Gamma$. Further details of how it is treated within PYTHIA are given in
Eq.(7.47) of Ref. Sjostrand:2006za . For the $q^{*}$ model, $\Gamma$ can be
written as

$\Gamma/M=\frac{\alpha_{s}}{3}+k\alpha\approx.03$ (2)

where $\alpha_{s}$ is the strong gauge coupling of the $q^{*}\rightarrow qg$
decay that produces dijets. In Eq. 2, $\alpha$ is the fine structure constant,
the coupling of electroweak interactions. It is multiplied by the factor
$k\approx 0.81$ to obtain the electroweak partial width $k\alpha$, for the sum
of the decays to an electroweak gauge boson ($q^{*}\rightarrow qW,qZ,\mbox{and
}q\gamma$). The $q^{*}$ width evolves only very slightly as a function of
$\sqrt{s}$, due to the logarithmic running of the strong coupling $\alpha_{s}$
as a function of the excited quark mass $M$, as shown in table 1. Here
$\alpha\approx 1/125$ at $M=6.5$ TeV, and the running of $\alpha$ is very
small and has negligible affect on the $q^{*}$ width. As a result, the di-
parton resonance width narrows only slightly with increasing $\sqrt{s}$ in
Fig. 1 (left). Finally, notice the long tail at low dijet mass, which arises
at di-parton level from the larger values of the PDFs at lower values of
fractional momenta $x$, making these distributions slightly asymmetric. This
effect is also discussed in the paragraphs following Eq.(7.49) in Ref.
Sjostrand:2006za . We note it is fairly independent of $\sqrt{s}$. Throughout
this paper we will discuss similar PDF effects on the distributions, and we
will elaborate more on their cause and size.

$\sqrt{s}$ | $M$ | $\Gamma/M$ | $\frac{q^{*}\rightarrow qg}{q^{*}\rightarrow\mbox{all}}$
---|---|---|---
(TeV) | (TeV) |  |
13 | 6.5 | 0.032 | 0.80
100 | 50 | 0.028 | 0.77
500 | 250 | 0.026 | 0.75
Table 1:  As a function of $\sqrt{s}$, we list a corresponding choice of
resonance mass $M/\sqrt{s}=0.5$, the full $q^{*}$ width $\Gamma/M$, and the
branching fraction for the decay to dijets $q^{*}\rightarrow qg$.

In Fig. 1 (middle) the wide genjet dijet mass distributions again show almost
no variation with collider energy in the narrow range $13<\sqrt{s}<27$ TeV
(bottom plot), and a small variation with collider energy in a very wide range
$13<\sqrt{s}<500$ TeV (top plot). Comparing with di-partons (left), we see
that for wide genjets the width has grown and the tail to low mass has
increased noticeably. At this value of M, the tail is no longer predominantly
due to PDFs, it is now dominated by the energy loss from very wide angle QCD
radiation emitted by the final state quark and gluon. Here the radiation is
also decreasing as $\sqrt{s}$ increases, again this is due to the logarithmic
decrease of $\alpha_{s}$ with $M$.

In Fig. 1 (right) the AK4 genjet dijet mass distributions show a little
variation with collider energy. Comparing to wide genjets (middle), we see
that for AK4 genjets both the width and the tail to low mass has increased
significantly. The narrow effective cone size $R=0.4$ of AK4 genjets results
in a significant energy loss from moderate angle QCD radiation emitted by the
final state quark and gluon, greatly increasing the tail to low mass and the
width of the AK4 genjet distribution. The wide genjet algorithm, with its
larger effective cone size $R=1.1$, collects more final state radiation in the
genjet than the AK4 genjet algorithm, and therefore results in dijet mass
distributions that are significantly closer to the di-parton mass distribution
and with less variation as a function of $\sqrt{s}$

Figure 2:  Dijet mass distributions for excited quarks from $pp$ collisions at
$\sqrt{s}$ equal to 13, 100, and 500 TeV, for resonance mass equal to 10%,
30%, and 50% of $\sqrt{s}$ (top row) and 20%, 40%, and 60% of $\sqrt{s}$
(bottom row), for the cases where the dijet mass is calculated from the
Lorentz vectors of the two final-state partons (left column), the two wide
genjets (middle column), and the two leading AK4 genjets in the event (right
column).

Figure 2 also demonstrates the approximate insensitivity of the dijet mass
distributions to the choice of $\sqrt{s}$, but does so for all six values of
excited quark mass we consider, in the range $0.1\leq M/\sqrt{s}\leq 0.6$.
Once again the distributions are of the dimensionless ratio $m/\sqrt{s}$,
which spaces the dijet mass distributions across the plot, because the dijet
mass, $m$, peaks at the excited quark mass, $M$, for each distribution. We can
already begin to see that there is something significantly different about the
distribution for $M/\sqrt{s}=0.6$, which we will explore in more detail in the
next section using a more appropriate choice of dimensionless variable to
study this effect. For di-partons (left), the long tail to low mass we had
noted previously for $M/\sqrt{s}=0.5$, has evolved into a longer tail with
significantly larger probability for $M/\sqrt{s}=0.6$. The same is true for
wide genjets (middle) and AK4 genjets (right). For the case of $\sqrt{s}=13$
TeV the $M/\sqrt{s}=0.6$ distribution for AK4 genjets is practically flat as a
function of dijet mass without a significant peak all the way to the kinematic
edge at $m/\sqrt{s}=0.6$.

###  4.2 Resonance mass invariance

Figure 3 demonstrates the approximate invariance of the dijet mass
distribution to changes in resonance mass. Here we show distributions of
another dimensionless ratio, dijet mass divided by resonance mass ($m/M$),
which allows a more direct comparison of the resonance shapes at different
values of resonance mass. For each of the three types of dijets reconstructed
at generator-level (di-partons, wide genjets, and AK4 genjets), regardless of
$\sqrt{s}$, the four distributions within the resonance mass range
$0.1<M/\sqrt{s}<0.4$ are almost indistinguishable, and for practical purposes
the five distributions within the slightly larger range $0.1<M/\sqrt{s}<0.5$
are approximately invariant.

Figure 3:  Dijet mass distributions for excited quarks with resonance mass
equal to 10%, 20%, 30%, 40%, 50% and 60% of $\sqrt{s}$, from $pp$ collisions
at $\sqrt{s}$ equal to 13 TeV (top row), 14 TeV (2nd row), 100 TeV (3rd row),
and 500 TeV (bottom row), for the cases where the dijet mass is calculated
from the Lorentz vectors of the two final-state partons (left column), the two
wide genjets (middle column), and the two leading AK4 genjets in the event
(right column).

Here we discuss why the shape of the dijet mass distribution is different at
high values of resonance mass. The low mass tail begins to become noticeable
for $M/\sqrt{s}=0.5$, and is significant for $M/\sqrt{s}=0.6$. This is
consistent with it’s source, dominated by the product of the resonance Breit-
Wigner and the steeply falling PDFs of the proton, which decrease rapidly as a
function of fractional momentum, $x$, and are much steeper at high values of
$x$. To appreciate this, note that the pole mass of the excited quark
resonance occurs at the dijet mass value $m=\sqrt{x_{1}x_{2}s}=M$. In Eq. 1,
the PDFs are steeply falling with increasing $x$, so at the small values of
$x$, corresponding to the low dijet mass tail of the resonance, the PDFs are
orders of magnitude larger than at the relatively high value $x\approx 0.5$.
This multiplies the Breit-Wigner tail at low dijet mass, and begins to
compensate the fall in the Breit-Wigner at low mass, creating an observable
low mass tail for $M/\sqrt{s}=0.5$. That steepness of the PDFs increases at
higher resonance mass, so that the ratio of PDF values between the tail and
the peak in dijet mass is roughly an additional order of magnitude larger at
$x\approx 0.6$ than they are at $x\approx 0.5$, making the low mass tail
practically a second peak at low mass, which is the dominant feature of the
excited quark resonance shapes for $M/\sqrt{s}=0.6$.

The other source of the low mass tail is final-state radiation, which
dominates the tail for $0.1<M/\sqrt{s}<0.4$, and roughly makes the same amount
of low mass tail independent of resonance mass. However, for $M/\sqrt{s}=0.6$,
the effect of PDFs discussed above is the dominant source of the low mass tail
in Fig. 3, and final-state radiation does not contribute significantly in
comparison. To see this, note that the probability density values on the tail
at very low dijet mass for di-partons, wide genjets, and AK4 genjets are
approximately the same. This can also be seen even clearer when we compare di-
partons and wide genjets directly in Fig. 6. It is only the height and width
of the resonance peak that are different for the three generator-level types
of dijets. The dijet mass distribution of di-partons is not affected by final
state radiation, because the mass is reconstructed from partons before the
radiation is emitted. The wide genjet distributions are moderately affected by
final state radiation, and the AK4 genjet distributions are significantly
affected by final state radiation. The observation that the probability
density values on the tail at very low mass for both types of genjets are
approximately the same as for di-partons, for excited quarks with
$M/\sqrt{s}=0.6$, demonstrates that at this value of resonance mass the tail
at very low dijet mass is not significantly affected by radiation, it is
dominated by PDF effects. The excited quark mass value $M/\sqrt{s}=0.5$ is
roughly a transition point, where the tail at very low mass has significant
contributions from both PDFs and final-state radiation, and neither is
dominant.

##  5 Experimental resolution

The establishment at the generator-level of the approximate invariance of the
dijet mass distribution to changes in $\sqrt{s}$ and $M$, for the range
$0.1<M/\sqrt{s}<0.5$, accomplishes most of what is needed to support our
Snowmass 2021 study Harris:2022kls , which assumed this invariance.
Nevertheless, we include estimates of experimental resolution, both to
estimate the expected experimental distributions at future $pp$ colliders, and
to derive additional understanding of the previously published distributions
from CMS CMS:2018xlo . In short, we will extract the effective CMS
experimental resolution for dijet resonances from excited quarks, generalize
that experimental resolution for all $pp$ colliders, and smear our wide genjet
distributions with that resolution to obtain experimental distributions at all
$pp$ colliders. This procedure has the advantage that it is simple, quick, and
easy to understand. It scales from an existing experiment to obtain a
plausible resolution that is sufficient for our study of future pp colliders.

###  5.1 Extraction of Gaussian experimental resolution for CMS

Full simulations of dijet resonances in $pp$ collisions at $\sqrt{s}=13$ TeV
are shown in Fig. 4(left) for excited quarks with resonance mass 2, 4, 6 and 8
TeV. The reconstruction used is the wide jet algorithm, equivalent to our wide
genjet algorithm, but using reconstructed experiment-level AK4 jets as input,
instead of generator-level AK4 genjets. The CMS publication CMS:2018xlo
discussed these dijet mass distributions as having “ _… Gaussian cores from
jet energy resolution, and tails towards lower mass values primarily from QCD
radiation_ ”, and presented values for the dijet mass resolution within the
Gaussian core for quark-quark and gluon-gluon resonances. These came from
fitting the distributions with truncated Gaussians, to estimate the core
resolution. We find that the corresponding values of the RMS deviation for the
Gaussian core of quark-gluon resonances in Fig. 4(left) is well represented by
the function

$\sigma_{\mbox{core}}(\%)=3.6\%+\frac{2.0\%}{\sqrt{M\mbox{[TeV]}}-0.5}.$ (3)

This function, the pink curve labeled "Gaussian core" in Fig. 4(middle), gives
RMS values for the Gaussian cores of quark-gluon resonances that are the
average of the RMS values presented for quark-quark and gluon-gluon resonances
in Ref. CMS:2018xlo . We assume this Gaussian core RMS deviation is a
convolution of two resolutions: the resolution of the underlying wide genjets
convolved with the true Gaussian resolution of the CMS detector for
experimentally observing wide jets. The measured half width at half maximum of
the wide genjets dijet mass distribution, labeled "Genjets (half width)" in
Fig. 4(middle), increases from 2.5% to 3.5% as $M$ increases from 2 to 8 TeV.
Deconvolving the Gaussian core RMS, by subtracting in quadrature the genjet
RMS, where genjet $\mbox{RMS}\approx\mbox{(half width)}/1.18$, gives the curve
labeled "Gaussian resolution" in Fig. 4(middle). This estimate of the Gaussian
experimental resolution of the CMS experiment for wide jets, in Run 2 at
$\sqrt{s}=13$ TeV, is well modeled by the function

$\sigma_{\mbox{gauss}}(\%)=1.5\%+\frac{5.4\%}{\sqrt{M\mbox{[TeV]}}}.$ (4)

Figure 4:  Extraction and re-application of experimental resolution. Left)
Reproduced from Ref. CMS:2018xlo  are the dijet mass distributions of fully
simulated dijet resonances in the CMS experiment, where excited quarks were
used to model quark-gluon resonances (blue). Middle) the resolution of the
Gaussian core of the CMS fully simulated dijet resonances (violet), the half-
width of the wide genjet mass distributions (green), and the approximate
Gaussian detector resolution resulting from the deconvolution of these two
resolutions (black). Right) the Gaussian distribution as a function of dijet
mass from detector resolution alone (black), the excited quark dijet mass
distribution for wide genjets (green), and the convolution of this generator-
level distribution with Gaussian detector resolution (violet).

###  5.2 Application of resolution and check of experimental distribution for
CMS

In Fig. 4(right) we show Gaussian distributions (black histograms), with RMS
equal to the experimental resolution in Eq. 4, at the four different values of
excited quark mass that were simulated in the CMS publication CMS:2018xlo .
Also shown are the dijet mass distributions of wide genjets (green histograms)
at the same values of excited quark mass, which we convolve with the Gaussian
resolution to obtain smeared wide genjets (pink histograms). This re-
application of the Gaussian experimental resolution to our generator-level
distributions, produces experiment-level dijet mass distributions, the smeared
wide genjets in Fig. 4 (right). We check how well our procedure works in Fig.
5, where we compare the smeared wide genjets (pink curves) with the fully-
simulated CMS wide jet dijet mass distributions (blue curves).

The level of agreement between our smeared wide genjets and the fully-
simulated CMS wide jets in Fig. 5 is adequate for our purpose in this paper,
namely to produce excited quark dijet mass distributions containing
approximate experimental resolution. Nevertheless, we can learn a little more
from the details of the comparison. The width of the peaks are only slightly
overestimated by this procedure. The small difference may be due to biases in
directly fitting the fully simulated distribution using truncated Gaussians,
to try and obtain the experimental Gaussian core resolution, since final-state
radiation will slightly widen the core of the distribution in addition to
producing the low mass tail. However, the long tails to low mass of the two
distributions are in good agreement. This is likely because the low mass tail
is produced almost completely by final-state radiation and PDFs at the
generator-level, as discussed in section 4, and is not affected by the
Gaussian experimental resolution. This is especially true for the excited
quark mass value 8 TeV, where the overall agreement between the two
distributions is very good, because the main feature of this shape is the long
tail to low dijet mass. Finally, we note a very small difference in the 2 TeV
shapes. The published shape from CMS is truncated at a dijet mass of 3 TeV,
while the full tail to high mass is shown for smeared wide genjets.

Figure 5:  Comparison of the fully simulated dijet resonance shapes of an
excited quark in the CMS experiment from Ref. CMS:2018xlo  (blue) and
predictions of the dijet resonance shapes using wide genjets smeared with an
estimate of the Gaussian resolution of the CMS experiment (violet).

###  5.3 Experimental resolution at pp Colliders

The estimated experimental resolution at $\sqrt{s}=13$ TeV in Eq. 4 can be
scaled to any pp collision energy by replacing $M\mbox{[TeV]}$ with
$M(13/\sqrt{s})$, and rewriting in terms of the dimensionless variable
$M/\sqrt{s}$:

$\sigma_{\mbox{gauss}}(\%)=1.5\%+\frac{1.5\%}{\sqrt{M/\sqrt{s}}}.$ (5)

Equation 5 is a generalized experimental resolution for the Gaussian core of
dijet resonances reconstructed with wide jets from excited quarks at pp
colliders. This resolution only depends on $M/\sqrt{s}$, so we are assuming
that the experimental resolution is constant for a fixed value of
$M/\sqrt{s}$. This is the same as assuming that future detectors will be
designed to give the same experimental resolution in percent as current
detectors, at the values of dijet mass appropriate for high mass search and
discovery, which is likely a conservative assumption. We apply this resolution
to the generator-level distributions (wide genjets) to obtain experiment-level
distributions (smeared wide genjets).

###  5.4 From generator-level to experiment-level resonances

Fig. 6 compares three levels of the excited quark dijet mass distributions. We
make this comparison for resonance masses in the range $0.4<M/\sqrt{s}<0.6$,
because that is near the critical mass $M/\sqrt{s}=0.5$ where $pp$ colliders
can exclude or discover an excited quark Harris:2022kls , and is also the mass
range where the generator-level distributions show a significant change in
shape, as discussed in section 4.2.

Figure 6:  Dijet mass distributions for excited quarks at parton level (blue),
genjet level (red) and experimental level (green) from $pp$ collisions at
$\sqrt{s}$ equal to 13 TeV (top row), 14 TeV (2nd row), 100 TeV (3rd row), and
500 TeV (bottom row), for resonance mass equal to 40% (left column), 50%
(middle column), and 60% of $\sqrt{s}$ (right column).

First, we show parton level, the generator level closest to the collision sub-
process, the dijet mass distribution of di-partons shown in Fig. 6 in blue.
Naturally it is the narrowest of the distributions, and for $M/\sqrt{s}=0.4$
it is a truly narrow Breit-Wigner-like distribution, with a full width at half
maximum $\Gamma\approx.03M$. At larger masses the product of the steeply
falling PDFs with the underlying Breit-Wigner process makes the distribution
noticeably asymmetric, as discussed in section 4.2. For $M/\sqrt{s}=0.5$ the
di-parton distribution has developed a long tail to low dijet mass, and for
$M/\sqrt{s}=0.6$ that tail is the dominant feature of the distribution,
containing the majority of the probability.

Second, we show genjet level, the generator-level after initial-state and
final-state radiation, the parton shower, hadronization, and clustering of
particles into jets, which results in the dijet mass distribution of wide
genjets shown in Fig. 6 in red. Wide genjets include most of the energy from
excited quark decay products, but energy from wide-angle final-state radiation
can be lost, when the particles from the hadronization of that radiation are
outside the $\Delta R=1.1$ radius of the wide genjets. As discussed in section
4.1, this energy loss is the primary cause of the widening of the peak, a very
small shift of the peak, and a significant tail to low mass, which are the
majority of the differences between the di-parton and wide genjet
distributions shown. As discussed in section 4.2, for $M/\sqrt{s}\leq 0.4$
this radiation is the dominant source of the tail at low mass,
$M/\sqrt{s}=0.5$ is a transition point, where radiation and the product of
PDFs with the underlying Breit-Wigner are roughly the same size contributions
to the tail, and for $M/\sqrt{s}=0.6$ radiation is only a small component of
the tail.

Third, we show the experiment-level estimates of fully reconstructed wide
jets, the dijet mass distributions of smeared wide genjets shown in Fig. 6 in
green. We have used Gaussians with RMS given by Eq. 5, to smear the wide
genjet dijet mass distributions. Fig. 6 shows that this Gaussian experimental
resolution widens the peak of the resonance shape but does not affect the
tails. This is because the experimental resolution is wider than the peaks of
these narrow resonances, but significantly narrower than the long tail to low
mass. The experimental width of the peak is determined by the experimental
resolution, as expected for narrow resonances. However, the long tail at low
mass is completely independent of experimental resolution. As discussed in the
paragraphs above and in sections 4.1 and 4.2, these tails are determined
completely at generator-level by the dijet reconstruction algorithm, the
initial-state and final-state radiation, the underlying Breit-Wigner line
shape, and the PDFs.

Figure 6 also demonstrates that the comparisons of the different levels of
dijet mass shapes are essentially independent of $\sqrt{s}$.

##  6 Experiment-level Resonances

The remaining figures in this paper will all be made with smeared wide
genjets, to approximate the expected measured shapes, cumulative probability
functions, and efficiencies of actual resonances observed by an experiment at
a pp collider.

Figure 7:  Comparisons of experiment level dijet mass distributions of excited
quarks, with dijet mass plotted as a fraction of $\sqrt{s}$, from $pp$
collisions at $\sqrt{s}$ equal to 13 TeV (blue), 100 TeV (red), and 500 TeV
(green) for resonance masses of 10%, 20% and 30% of $\sqrt{s}$ (top row) and
40%, 50% and 60% of $\sqrt{s}$ (bottom row).

###  6.1 Collision energy and resonance mass invariances at experiment-level

The approximate invariance of the dijet mass distribution at the experiment-
level to the choice of $\sqrt{s}$ is presented in Fig. 7 and Fig. 8. These two
figures differ only in the choice of the dimensionless ratio plotted on the
horizontal axis: $m/\sqrt{s}$ and $m/M$. They both demonstrate the approximate
invariance of the experiment-level resonance shapes to variations in
$\sqrt{s}$, which naturally follows from the invariance at the generator-level
observed in Fig 2 and discussed in section 4.1, and the invariance of the
experimental resolution with $\sqrt{s}$, which was explicitly assumed in
section 5.3. The change in shape of the experiment-level resonances when
$\sqrt{s}$ increases from 13 to 100 TeV, although small, is significantly more
than the change in shape when $\sqrt{s}$ increases from 100 to 500 TeV. This
is because the decrease in width due to logarithmic evolution with
$\alpha_{s}$ of the generator-level half width is being damped by convolution
with a dominant and constant experimental resolution.

Figure 8:  Comparisons of experiment level dijet mass distributions of excited
quarks, with dijet mass plotted as a fraction of resonance mass, from $pp$
collisions at $\sqrt{s}$ equal to 13 TeV (blue), 100 TeV (red), and 500 TeV
(green) for resonance masses of 10%, 20% and 30% of $\sqrt{s}$ (top row) and
40%, 50% and 60% of $\sqrt{s}$ (bottom row).

The invariance of the dijet mass distribution at the experiment-level to the
choice of resonance mass $M$, approximately valid for the mass range
$0.1<M/\sqrt{s}<0.5$, is shown in Fig. 9. For all four values of $\sqrt{s}$
shown, only the resonance mass value $M/\sqrt{s}=0.6$ has a significantly
different dijet mass distribution, because the significant probability within
the long tail to low mass breaks the invariance. We discuss here the evolution
of the shape with $M/\sqrt{s}$ for both the peak and the tails.

On the peak of the dijet mass distribution, the probability density increases
with $M/\sqrt{s}$ for the resonance masses in the range $0.1<M/\sqrt{s}<0.4$.
This is because the experimental resolution in Eq. 5 improves as $M/\sqrt{s}$
increases. The highest peak is always for $M/\sqrt{s}=0.4$, because when the
mass reaches $M/\sqrt{s}=0.5$ noticeable numbers of events have migrated from
the peak towards the low mass tail, so that the peak probability density for
$M/\sqrt{s}=0.5$ always lies between the peaks for $M/\sqrt{s}=0.2$ and
$M/\sqrt{s}=0.3$. The migration is caused by the steepening PDFs, discussed in
section 4.2 and 5.4, that dominate the distribution for $M/\sqrt{s}=0.6$ which
therefore has the lowest peak probability density. Despite this issue with the
low mass tail, the width of the peak is steadily decreasing as $M/\sqrt{s}$
increases within the full range $0.1<M/\sqrt{s}<0.6$, because the peak is
dominated by the experimental resolution.

On the tail at high dijet mass, the probability density decreases with
increasing $M/\sqrt{s}$ for all the resonance masses in the range
$0.1<M/\sqrt{s}<0.6$. The high mass tail has the natural shape of the
underlying Breit-Wigner distribution, slightly enhanced by initial state
radiation and experimental resolution, both of which decrease with increasing
$M/\sqrt{s}$.

On the tail at low dijet mass, the probability density is virtually
indistinguishable for $0.1<M/\sqrt{s}<0.3$, because the process of final state
radiation which dominates the tail in this region is independent of
$M/\sqrt{s}$. The first distribution to show a very small increase in the low
mass tail due to the steepening PDFs is $M/\sqrt{s}=0.4$, which remains
dominated by final-state radiation. Once again, as discussed in section 4.2
and 5.4, the value $M/\sqrt{s}=0.5$ is the transition mass where the low mass
tail is roughly equally composed of the effects of final-state radiation and
the steepening PDFs, and for $M/\sqrt{s}=0.6$ the low mass tail is dominated
by the effect of the steepening PDFs. Nowhere within the low mass tail are
there any indications of the effects of experimental resolution, which are
only apparent in the height of the peak, the width of the peak, and possibly
the high-mass tail.

Figure 9:  Comparisons of experiment level dijet mass distributions of excited
quarks, with resonance mass equal to 10%, 20%, 30%, 40%, 50% and 60% of
$\sqrt{s}$, from $pp$ collisions at $\sqrt{s}$ equal to 13 and 14 TeV (top
row), and 100 and 500 TeV (bottom row).

###  6.2 Cumulative distributions of experiment-level resonances

Cumulative probability distributions allow easier visualization and extraction
of the acceptance of requirements made on dijet mass in a resonance search.
These distributions give the total probability, $p$, arising from integrating
the probability density distribution, $dp/dm$, between a dijet mass of $0$ and
$m$: $\int(dp/dm)dm$. They have the value $p=0$ for $m=0$ and the value $p=1$
for large enough $m$ to contain the entire distribution. Cumulative
probability distributions as a function of dijet mass, using the dimensionless
ratio $m/M$, are shown in Fig. 10 and Fig. 11.

Figure 10 shows the cumulative probability from integrating Fig. 9 between 0
and a given value of the dimensionless ratio $m/M$. To within a probability of
roughly 10%, the cumulative probability is approximately invariant as a
function of resonance mass for the mass range $0.1<M/\sqrt{s}<0.5$.

Similarly, Fig. 11 shows the cumulative probability from integrating Fig. 8
between 0 and a given value of the dimensionless ratio $m/M$. Again, to within
a probability of roughly 10%, the cumulative probability is approximately
invariant as a function of collision energy $\sqrt{s}$.

Figure 10:  Cumulative probability distributions from experiment level dijet
mass distributions of excited quarks, for resonance mass $M$ equal to 10%,
20%, 30%, 40%, 50% and 60% of $\sqrt{s}$, from $pp$ collisions at $\sqrt{s}$
equal to 13, 14, and 27 TeV (top row), and 100, 300 and 500 TeV (bottom row).

Figure 11:  Cumulative probability distributions for experiment level dijet
mass distributions of excited quarks, from $pp$ collisions at $\sqrt{s}$ equal
to 13 TeV (blue), 100 TeV (red), and 500 TeV (green), for resonance masses of
10%, 20% and 30% of $\sqrt{s}$ (top row) and 40%, 50% and 60% of $\sqrt{s}$
(bottom row).

###  6.3 Acceptance of a dijet mass window

A classic method to search for dijet resonance is to look for a bump in the
dijet mass distribution within a window in dijet mass centered on the
resonance peak. This is because the majority of the sensitivity comes from a
narrow interval in dijet mass surrounding the peak. In our previous paper
Harris:2022kls  we used a dijet mass window, $0.836<m/M<1.164$, to estimate
search sensitivity. In that paper, the acceptance of that mass window was
found to be $A=0.6$, for a dijet resonance from an excited quark at $M=6$ TeV
and $\sqrt{s}=13$ TeV. We argued that acceptance was independent of resonance
mass and collision energy, and used that acceptance to estimate the
sensitivity to dijet resonances from excited quarks at $pp$ colliders.

Figure 12 shows the $q^{*}$ acceptance of the dijet mass window,
$0.836<m/M<1.164$, centered on the experiment level dijet resonance shapes in
Fig. 9. The acceptance in Fig. 12 can also be easily extracted from Fig. 10 by
reading off the cumulative probabilities at the window boundaries in the ratio
of dijet mass to resonance mass, $m/M=1.164$ and $m/M=0.836$, and taking the
difference. Similarly, the acceptance for any choice of dijet mass window can
be derived from Fig. 10 by taking the difference between the cumulative
probabilities at the upper edge and lower edge of the mass window chosen.

Figure 12:  The experimental acceptance of the dijet mass window used in Ref.
Harris:2022kls , as a function of resonance mass, for excited quarks from $pp$
collisions at $\sqrt{s}$ equal to 13, 14, 100 and 500 TeV.

Figure 12 shows the acceptance at four values of $\sqrt{s}$ as a function of
resonance mass, which we will denote as $A_{\sqrt{s}}(M/\sqrt{s})$. There is
less than about 10% variation in window acceptance within the resonance mass
range $0.1<M/\sqrt{s}<0.5$, approximately invariant as a function of resonance
mass for the same mass interval where the resonance shapes and cumulative
probabilities are approximately invariant as a function of resonance mass. The
acceptance increases slowly with $\sqrt{s}$, $A_{13}<A_{14}<A_{100}<A_{500}$,
but again doesn’t change by more than about 10%, approximately invariant as a
function of $\sqrt{s}$. For the lowest $pp$ collision energy, LHC at
$\sqrt{s}=13$ TeV, almost indistinguishable from HL-LHC at $\sqrt{s}=14$ TeV,
the acceptance has a maximum value of $A_{13}(0.3)=0.68$ near $M/\sqrt{s}=0.3$
and is $A_{13}(0.5)=0.58$ for $M/\sqrt{s}=0.5$, which does not result in
measurably different sensitivities than using $A=0.6$ as we did in our paper
Harris:2022kls . For the highest $pp$ collision energy, $\sqrt{s}=500$ TeV,
the acceptance has a maximum value of $A_{500}(0.3)=0.76$ at $M/\sqrt{s}=0.3$
and is $A_{500}(0.5)=0.67$ for $M/\sqrt{s}=0.5$, which again would not result
in a measurable gain in sensitivity if we had used these values instead of the
more conservative single choice $A=0.6$. However, at $M\sqrt{s}=0.6$ the
window acceptance for $\sqrt{s}=13$ and $14$ TeV has fallen to
$A_{13}(0.6)=A_{14}(0.6)=0.32$, less than half its maximum value, and the
window acceptance for $\sqrt{s}=100$ and $500$ TeV has fallen to
$A_{100}(0.6)=A_{500}(0.6)=0.40$, just greater than half its maximum value. As
discussed in Ref. Harris:2022kls , this loss of acceptance means our quoted
sensitivities for excited quarks with resonance masses near $M/\sqrt{s}=0.6$
need to be reduced by roughly 6%, and our estimate of $q^{*}$ excluded mass at
95% confidence level for HL-LHC with 3 ab-1 in Refs. Harris:2022kls ;
Bernardi:2022hny  needs to be reduced from $7.9$ TeV to roughly $7.6$ TeV.

##  7 Summary and conclusions

We have presented and explained in detail the expected experimental dijet mass
distributions from excited quarks at present and future $pp$ colliders. The
natural width of the peak of the experimental distribution narrows only
slightly with the collider energy $\sqrt{s}$ due to the running of the strong
coupling constant. The expected width of the observed peak is consistent with
the natural width convolved with a Gaussian experimental resolution. We have
estimated and parameterized this resolution for wide jets from the CMS
experiment, and scaled with $\sqrt{s}$ to estimate the resolution of future
experiments. The expected experimental dijet mass distributions have a long
tail at low dijet mass which is formed from two components, both of which come
from physics at the generator-level. The first component of the low mass tail,
which is visible and approximately invariant for all resonance masses, is the
energy lost due to final-state radiation. That tail is significantly reduced
by using the wide jet algorithm to reconstruct dijets, as opposed to using
jets with a smaller cone size. This tail from final-state radiation interferes
with accurate extraction of the Gaussian experimental resolution, making it
approximate. The second component of the low mass tail, which is visible only
at high resonance masses, is from the multiplication of the tail of the Breit-
Wigner line shape by steeply falling parton momentum distributions of the
proton, which are significantly steeper at high resonance mass. This second
component contributes as much as final state radiation to the low dijet mass
tail for resonance masses equal to half of $\sqrt{s}$, and dominates the dijet
mass shape completely for resonance masses greater than 60% of $\sqrt{s}$.
Finally, the dijet mass distributions also have a short and small tail to high
mass which becomes longer and larger at lower resonance masses, consistent
with the effects of the natural line shape, initial state radiation and
experimental resolution.

We have also presented cumulative probability distributions as a function of
dijet mass, that can be used to extract the acceptance of requirements made on
dijet mass within an excited quark search. As an example of how to use these
cumulative probability distributions, we have presented the excited quark
signal acceptance of the dijet mass window which we used in our Snowmass 2021
study on the sensitivity to dijet resonances at $pp$ colliders Harris:2022kls
.

We find an approximate universality of the differential dijet mass
distributions, the cumulative probability distributions, and the window
acceptance, which are all approximately invariant to changes in $\sqrt{s}$ and
resonance mass for excited quark masses between 10% and 50% of $\sqrt{s}$.
This universality, arising from the relatively constant natural width of the
resonance and the approximate invariance of final-state radiation to changes
in $\sqrt{s}$ and resonance mass, is violated by the steepening falloff of the
parton momentum distributions of the proton when the resonance mass reaches
and exceeds 60% of $\sqrt{s}$.

###### Acknowledgements.

This work was supported by the Turkish Energy, Nuclear and Mineral Research
Agency. This work was also supported by the Fermi National Accelerator
Laboratory, managed and operated by Fermi Research Alliance, LLC under
Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy. The U.S.
Government retains and the publisher, by accepting the article for
publication, acknowledges that the U.S. Government retains a non-exclusive,
paid-up, irrevocable, world-wide license to publish or reproduce the published
form of this manuscript, or allow others to do so, for U.S. Government
purposes.

## References

  * (1) R.M. Harris, E. Gurpinar Guler and Y. Guler, _Sensitivity to Dijet Resonances at Proton-Proton Colliders_ , in _2022 Snowmass Summer Study_ [2202.03389].
  * (2) R.M. Harris and K. Kousouris, _Searches for dijet resonances at hadron colliders_ , _Int. J. Mod. Phys. A_ 26 (2011) 5005 [1110.5302].
  * (3) U. Baur, I. Hinchliffe and D. Zeppenfeld, _Excited quark production at hadron colliders_ , _Int. J. Mod. Phys. A_ 02 (1987) 1285.
  * (4) U. Baur, M. Spira and P.M. Zerwas, _Excited quark and lepton production at hadron colliders_ , _Phys. Rev. D_ 42 (1990) 815.
  * (5) ATLAS Collaboration, _Search for new physics in dijet mass and angular distributions in $pp$ collisions at $\sqrt{s}=7$ TeV measured with the ATLAS detector_, _New J. Phys._ 13 (2011) 053044 [1103.3864].
  * (6) CMS Collaboration, _Search for dijet resonances in 7 TeV pp collisions at CMS_ , _Phys. Rev. Lett._ 105 (2010) 211801 [1010.0203].
  * (7) ATLAS Collaboration, _Search for new particles in two-jet final states in 7 TeV proton-proton collisions with the ATLAS detector at the LHC_ , _Phys. Rev. Lett._ 105 (2010) 161801 [1008.2461].
  * (8) CMS Collaboration, _Search for resonances in the dijet mass spectrum from 7 TeV pp collisions at CMS_ , _Phys. Lett. B_ 704 (2011) 123 [1107.4771].
  * (9) ATLAS collaboration, _ATLAS search for new phenomena in dijet mass and angular distributions using $pp$ collisions at $\sqrt{s}=7$ TeV_, _JHEP_ 01 (2013) 029 [1210.1718].
  * (10) CMS collaboration, _Search for Narrow Resonances and Quantum Black Holes in Inclusive and $b$-Tagged Dijet Mass Spectra from $pp$ Collisions at $\sqrt{s}=7$ TeV_, _JHEP_ 01 (2013) 013 [1210.2387].
  * (11) CMS collaboration, _Search for Narrow Resonances Using the Dijet Mass Spectrum in $pp$ Collisions at $\sqrt{s}$=8 TeV_, _Phys. Rev._ D87 (2013) 114015 [1302.4794].
  * (12) ATLAS collaboration, _Search for new phenomena in the dijet mass distribution using $p-p$ collision data at $\sqrt{s}=8$ TeV with the ATLAS detector_, _Phys. Rev._ D91 (2015) 052007 [1407.1376].
  * (13) CMS collaboration, _Search for resonances and quantum black holes using dijet mass spectra in proton-proton collisions at $\sqrt{s}=$ 8 TeV_, _Phys. Rev._ D91 (2015) 052009 [1501.04198].
  * (14) CMS collaboration, _Search for narrow resonances decaying to dijets in proton-proton collisions at $\sqrt{(}s)=$ 13 TeV_, _Phys. Rev. Lett._ 116 (2016) 071801 [1512.01224].
  * (15) ATLAS collaboration, _Search for new phenomena in dijet mass and angular distributions from $pp$ collisions at $\sqrt{s}=$ 13 TeV with the ATLAS detector_, _Phys. Lett._ B754 (2016) 302 [1512.01530].
  * (16) CMS collaboration, _Search for narrow resonances in dijet final states at $\sqrt{(}s)=$ 8 TeV with the novel CMS technique of data scouting_, _Phys. Rev. Lett._ 117 (2016) 031802 [1604.08907].
  * (17) CMS collaboration, _Search for dijet resonances in proton-proton collisions at $\sqrt{s}$ = 13 TeV and constraints on dark matter and other models_, _Phys. Lett._ B769 (2017) 520 [1611.03568].
  * (18) ATLAS collaboration, _Search for new phenomena in dijet events using 37 fb -1 of $pp$ collision data collected at $\sqrt{s}=$13 TeV with the ATLAS detector_, _Phys. Rev._ D96 (2017) 052004 [1703.09127].
  * (19) ATLAS collaboration, _Search for low-mass dijet resonances using trigger-level jets with the ATLAS detector in $pp$ collisions at $\sqrt{s}=13$ TeV_, _Phys. Rev. Lett._ 121 (2018) 081801 [1804.03496].
  * (20) CMS collaboration, _Search for narrow and broad dijet resonances in proton-proton collisions at $\sqrt{s}=13$ TeV and constraints on dark matter mediators and other new particles_, _JHEP_ 08 (2018) 130 [1806.00843].
  * (21) ATLAS collaboration, _Search for new resonances in mass distributions of jet pairs using 139 fb -1 of $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _arXiv:1910.08447, submitted to JHEP_ (2019) [1910.08447].
  * (22) CMS collaboration, _Search for high mass dijet resonances with a new background prediction method in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _arXiv:1911.03947, submitted to JHEP_ (2019) [1911.03947].
  * (23) FCC collaboration, _HE-LHC: The High-Energy Large Hadron Collider: Future Circular Collider Conceptual Design Report Volume 4_ , _Eur. Phys. J. ST_ 228 (2019) 1109.
  * (24) P.C. Bhat et al., _Future Collider Options for the US_ , in _2022 Snowmass Summer Study_ [2203.08088].
  * (25) FCC collaboration, _FCC-hh: The Hadron Collider: Future Circular Collider Conceptual Design Report Volume 3_ , _Eur. Phys. J. ST_ 228 (2019) 755.
  * (26) P.M. McIntyre, S. Assadi, S. Bannert, J. Breitschopf, D. Chavez, J. Gerity et al., _Collider in the Sea: Vision for a 500 TeV World Laboratory_ , in _2nd North American Particle Accelerator Conference_ , p. MOB2CO03, 2017, DOI.
  * (27) T. Sjöstrand, S. Ask, J.R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159 [1410.3012].
  * (28) J. Pumplin, D.R. Stump, J. Huston, H.-L. Lai, P. Nadolsky and W.-K. Tung, _New generation of parton distributions with uncertainties from global QCD analysis_ , _JHEP_ 07 (2002) 012 [hep-ph/0201195].
  * (29) CMS Collaboration, _Event generator tunes obtained from underlying event and multiparton scattering measurements_ , _Eur. Phys. J. C_ 76 (2016) 155 [1512.00815].
  * (30) P. Skands, S. Carrazza and J. Rojo, _Tuning PYTHIA 8.1: the Monash 2013 tune_ , _Eur. Phys. J. C_ 74 (2014) 3024 [1404.5630].
  * (31) M. Cacciari, G.P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063 [0802.1189].
  * (32) M. Cacciari, G.P. Salam and G. Soyez, _FastJet user manual_ , _Eur. Phys. J. C_ 72 (2012) 1896 [1111.6097].
  * (33) T. Sjostrand, S. Mrenna and P.Z. Skands, _PYTHIA 6.4 Physics and Manual_ , _JHEP_ 05 (2006) 026 [hep-ph/0603175].
  * (34) G. Bernardi et al., _The Future Circular Collider: a Summary for the US 2021 Snowmass Process_ , 2203.06520.


# Data-Driven Encoding: A New Numerical Method for Computation of the Koopman
Operator

Jerry Ng1, H. Harry Asada2, This material is based upon work supported by
National Science Foundation Grant NSF-CMMI 2021625.1Jerry Ng is with
Department of Mechanical Engineering at the Massachusetts Institute of
Technology<EMAIL_ADDRESS>Harry Asada is with Department of Mechanical
Engineering at the Massachusetts Institute of Technology<EMAIL_ADDRESS>
###### Abstract

This paper presents a data-driven method for constructing a Koopman linear
model based on the Direct Encoding (DE) formula. The prevailing methods,
Dynamic Mode Decomposition (DMD) and its extensions are based on least squares
estimates that can be shown to be biased towards data that are densely
populated. The DE formula consisting of inner products of a nonlinear state
transition function with observable functions does not incur this biased
estimation problem and thus serves as a desirable alternative to DMD. However,
the original DE formula requires knowledge of the nonlinear state equation,
which is not available in many practical applications. In this paper, the DE
formula is extended to a data-driven method, Data-Driven Encoding (DDE) of
Koopman operator, in which the inner products are calculated from data taken
from a nonlinear dynamic system. An effective algorithm is presented for the
computation of the inner products, and their convergence to true values is
proven. Numerical experiments verify the effectiveness of DDE compared to
Extended DMD. The experiments demonstrate robustness to data distribution and
the convergent properties of DDE, guaranteeing accuracy improvements with
additional sample points. Furthermore, DDE is applied to deep learning of the
Koopman operator to further improve prediction accuracy.

##  I Introduction

Dynamic Mode Decomposition (DMD) was presented as a method to produce linear
models from data generated through nonlinear dynamical processes by using
Singular Value Decomposition (SVD) [1]. Later, this method was developed
further to create Extended Dynamic Mode Decomposition (EDMD), which introduced
the concept of using observable functions, nonlinear functions of state
variables, as a method of augmenting the state space [2]. EDMD referenced the
Koopman Operator as justification and a theoretical underpinning for lifting
the state space. In his seminal work, Bernard Koopman showed the existence of
this operator that transforms nonlinear systems into linear systems [3].
Another extension of DMD has shown the viability of using DMD for control on
non-autonomous systems [4]. This enabled complex nonlinear Model Predictive
Control (MPC) to be converted to linear MPC [5], leading to numerous studies
utilizing the methodology to real systems[6, 7, 8, 9].

To improve the accuracy of the models based on the Koopman Operator, two
avenues of research have formed. The first avenue regards the selection of the
observable functions used for constructing a lifted state space, as these
functions are a key ingredient in creating an accurate linear model. Various
methods have been developed, including deep neural networks for learning
effective observable functions [10, 11, 12] and optimization [13]. However, an
efficient selection of observables does not solve all the issues that arise
when attempting to construct an accurate linear model. For example, it is
known that unstable modes are involved in Koopman-based DMD models and their
extensions although the underlying nonlinear systems are known to be stable.
The second avenue of research involves the formulation of the linear
transition matrix. Extensive studies have been done to create stable linear
models to remedy the situations where an outright use of DMD would lead to the
creation of an unstable linear model [14, 15, 16, 17]. Recently, an extension
of DMD, called Robust Dynamic Mode Decomposition (RDMD), utilizes statistical
measures to suppress the effect of outliers on modeling the linear Koopman
matrix[18].

A fundamental difficulty in constructing a proper linear model is data
dependency. Least Squares Estimation (LSE), involved in all DMD based methods,
often produces a significant bias due to the distribution of the dataset. To
eliminate this dependency on data distributions, the current work takes an
alternative approach to LSE.

Recently, a new formulation of the Koopman Operator, termed Koopman Direct
Encoding (DE), was produced [19]. This method directly encodes the nonlinear
function of state transition using observables as basis functions to obtain a
Koopman linear model. Inner products of observable functions in composition
with the nonlinear state transition function are used to construct the state
transition matrix without use of LSE. While DE theoretically guarantees the
exact linear model, it requires access to the nonlinear state equations, which
are often not available in practical applications. The current work aims to
fill the gap between DE and data-driven approaches.

There are four significant contributions presented in this work. The first is
the conversion of the DE formula of the Koopman Operator to a data-driven
formula. The second is a computational algorithm and proof of its convergence
to the true inner products that constitute the DE formula. The third is
numerical experiments that provide evidence that the proposed method, unlike
EDMD, does not exhibit biases to data distribution, but can produce
consistently higher accuracy compared to EDMD. Finally, the DDE algorithm is
utilized in modeling a high order nonlinear system in combination with deep
learning.

##  II Koopman Operator and the Direct Encoding Method

In this section, we give a brief overview of the Koopman Operator and dynamic
mode decomposition, and introduce the direct encoding method for obtaining a
Koopman operator directly from nonlinear dynamics.

###  II-A Least Squares Estimation of the Koopman operator

Consider a discrete-time dynamical system, given by

$x_{t+1}=f(x_{t})$ (1)

where $x\in X\subset\mathbb{R}^{n}$ is the independent state variable vector
representing the dynamic state of the system, $f$ is a self-map, nonlinear
function $f:X\rightarrow X$, and $t$ is the current time step. Also consider a
real-valued observable function of the state variables
$g:X\rightarrow\mathbb{R}$. The Koopman Operator $\mathcal{K}$ is an infinite-
dimensional linear operator acting on the observable function $g$ :

$\mathcal{K}g=g\circ f$ (2)

where $g\circ f$ is the composition of function $g$ with function $f$:
$(g\circ f)(x)=g(f(x))$.

A common data-driven method for constructing the operator is Extended Dynamic
Mode Decomposition (EDMD) [2], where observables that are experimentally
obtained or simulated from the governing equation of the system are augmented
by including real-valued observable functions of the independent state vector
$x_{t}$. Collectively, a high-dimensional state vector is formed.

$z_{t}=\begin{bmatrix}g_{1}(x_{t})\\\ g_{2}(x_{t})\\\ \vdots\\\
g_{m}(x_{t})\end{bmatrix}$ (3)

where $m$ is the order of the lifted state corresponding to the number of
observable functions. Underpinned by the Koopman Operator theory, EDMD assumes
the existence of a linear state transition matrix $A$ relating $z_{t+1}$ to
$z_{t}$, and determine $A$ by solving a least squares regression that
minimizes the Sum of Squared Error (SSE).

$A=\arg\min_{A}\sum_{t}||z_{t+1}-Az_{t}||^{2}$ (4)

Singular Value Decomposition (SVD) is used for the least squares optimization.

###  II-B Direct Encoding of the Koopman Operator

An alternative to the least squares estimate and EDMD is to obtain the exact
$A$ matrix by directly encoding the self-map, nonlinear state transition
function $f(x)$ with an independent and complete set of observable functions
through inner product computations. This Direct Encoding method is introduced
next, while the full proof can be found in [19].

Let us first consider the case where $g_{1},g_{2},g_{3},...$ are orthonormal
basis functions spanning a Hilbert space $\mathcal{H}$. We assume that the
self-map nonlinear function $f(x)$ is continuous and that the composition of
$g_{j}$ with $f$ is also involved in the Hilbert space.

$g_{j}\circ f\in\mathcal{H}\quad\forall j$ (5)

This implies that the function $g_{j}\circ f$ can be expanded in
$[g_{1},g_{2},g_{3},...]$.

$g_{j}\circ f=\sum_{k=1}^{\infty}\langle g_{j}\circ f,g_{k}\rangle g_{k}$ (6)

Concatenating $g_{1},g_{2},g_{3},...$ and $g_{1}\circ f,g_{2}\circ
f,g_{3}\circ f...$ in infinite dimensional column vectors, respectively,

$\bar{z}_{t}=\begin{bmatrix}g_{1}(x_{t})\\\ g_{2}(x_{t})\\\
\vdots\end{bmatrix},\bar{z}_{t+1}=\begin{bmatrix}g_{1}[f(x_{t})]\\\
g_{2}[f(x_{t})]\\\ \vdots\end{bmatrix}$ (7)

eq. (6) can be written in matrix and vector form.

$\bar{z}_{t+1}=\bar{A}\bar{z}_{t}$ (8)

where $\bar{A}$ is an infinite dimensional matrix consisting of the inner
products involved in eq. (6),

$\bar{A}=\begin{bmatrix}\langle g_{1}\circ f,g_{1}\rangle&\langle g_{1}\circ
f,g_{2}\rangle&\dots\\\ \langle g_{2}\circ f,g_{1}\rangle&\langle g_{2}\circ
f,g_{2}\rangle&\dots\\\ \vdots&\vdots&\ddots\end{bmatrix}$ (9)

Eq. (8) manifests that the state lifted to the infinite dimensional space
$\bar{z_{t}}$ makes linear state transition with matrix $\bar{A}$.

The observables $g_{1},g_{2},g_{3},...$ were assumed to be orthonormal basis
functions in the above derivation. This assumption can be relaxed to an
independent and complete set of basis functions spanning the Hilbert space.
Hereafter, let $[g_{1},g_{2},g_{3},...]$ be an independent and complete set of
basis functions spanning the Hilbert space.

It can be shown that the time evolution of lifted state $z_{t}$ is given by a
constant matrix $A_{f}$ for the independent and complete set of basis
functions $[g_{1},g_{2},g_{3},...]$.

$z_{t+1}=A_{f}z_{t}$ (10)

The matrix $A_{f}$ can be computed directly from the self-map, state
transition function $f(x)$ and an independent and complete set of observables
$[g_{1},g_{2},g_{3},...]$ through inner product computations. Post-multiplying
the transpose of $z_{t}$ to both sides of eq. (10) and integrating them over
$X$ yield:

$\int_{X}z(f(x))z^{T}(x)dx=A_{f}\int_{X}z(x)z^{T}(x)dx$ (11)

which can be written as

$Q=A_{f}R$ (12)

where

$\displaystyle Q$ $\displaystyle=\begin{bmatrix}\langle g_{1}\circ
f,g_{1}\rangle&\langle g_{1}\circ f,g_{2}\rangle&\dots\\\ \langle g_{2}\circ
f,g_{1}\rangle&\langle g_{2}\circ f,g_{2}\rangle&\dots\\\
\vdots&\vdots&\ddots\end{bmatrix}$ (13) $\displaystyle R$
$\displaystyle=\begin{bmatrix}\langle g_{1},g_{1}\rangle&\langle
g_{1},g_{2}\rangle&\dots\\\ \langle g_{2},g_{1}\rangle&\langle
g_{2},g_{2}\rangle&\dots\\\ \vdots&\vdots&\ddots\end{bmatrix}$ (14)

Because the observables are independent, the matrix $R$ is non-singular.
Therefore, the matrix $A_{f}$ is given by

$A_{f}=QR^{-1}$ (15)

This formula for obtaining the matrix $A_{f}$ directly from the governing
nonlinear state equation with the function $f(x)$ and the independent
observables through inner products, which are guaranteed to exist in Hilbert
space $\mathcal{H}$, is the Direct Encoding method.

##  III Data-Driven Koopman Encoding

The prevailing method for construction of the Koopman Operator, EDMD, is based
on LSE and SVD. This method, however, cannot provide a consistent estimate;
the result is highly dependent on the distribution of data as the Koopman
Operator is being approximated [20]. This dependency on distribution occurs
because a core assumption of LSE is that the model structure is correct; when
this assumption is violated, LSE is unable to create an unbiased estimator
[21]. As the Koopman operator is truncated in practical use, this assumption
does not hold.

Non-uniform data distributions inevitably occur in practical applications. For
a nonlinear dynamical system with a stable equilibrium, for example, data
collected from experiments and/or simulation of the system tend to be dense in
the vicinity of the equilibrium, as all trajectories that begin within a
region of attraction converge to the equilibrium. Because LSE applies equal
weighting to all data points, the model is heavily tuned to the behavior of
the densely populated region.

The Direct Encoding method described previously enables us to obtain the exact
linear state transition matrix $A$ through inner product computations. As the
formulation is based on integration over the entire state space, there is no
bias towards particular parts of the domain.

However, the original form of the Direct Encoding method utilizes the
nonlinear state equation, i.e. the self-map $f(x)$, to compute the inner
products. In practical applications, such a nonlinear function is not always
available; only data are available. The objective of this section is to
establish a computational algorithm to obtain the $A$ matrix by numerically
computing the inner products, $\langle g_{i},g_{j}\rangle,\langle g_{i}\circ
f,g_{j}\rangle$, from a given set of data.

The method presented consists of three operations.

  * •

The integral of the inner products is reduced in range from the entire state
space to the dynamic range encapsulated by the data.

  * •

The dynamic range is discretized with data points.

  * •

The inner product integral is reduced to a weighted summation of the integrand
evaluated at each data point multiplied by the volume $\Delta v$ associated to
each point.

Naturally, if data are densely populated in a small region, the discretized
integral interval is small and thereby the volume also becomes small.
Similarly, the volume tends to be larger where the data are sparse. In the
summation, the integrand evaluated at individual data points are ”weighted” by
the size of the volume. This numerical inner product calculation prevents
overemphasis of clustered data.

###  III-A Inner Product Computation

We present the data-driven encoding method (DDE) as an alternative data-driven
method to DMD for calculating a finite order approximation of the Koopman
Operator. The objective of this method is to compute the matrices $R$ and $Q$
in (13) and (14) from data. This entails the computation of inner products:

$\displaystyle\langle g_{i},g_{j}\rangle=\int_{X}G_{ij}(\xi)d\xi$ (16)
$\displaystyle\langle g_{i}\circ f,g_{j}\rangle=\int_{X}F_{ij}(\xi)d\xi$ (17)

where

$G_{ij}(x)=g_{i}(x)\bar{g}_{j}(x),F_{ij}(x)=g_{i}[f(x)]\bar{g}_{j}(x)$ (18)

are assumed to be Riemann Integrable; the functions are bounded and continuous
[22].

There are two data sets used for the inner product computation. The first data
set is

$D_{N}=\\{x_{i}\mid i=1,\cdots,N;x_{i}\in X\\}$ (19)

Note that all the data values are finite, $|x_{i}|<\infty$. As such, the
integral interval of the inner products is finite in computing them from the
data. To define the integral interval, we consider the dynamic range of the
system, $X_{D}$, determined from the data set $D_{N}$. See Fig.1. The dynamic
range $X_{D}$ is defined to be the minimum domain in the space $X$ that
includes all the data points in $D_{N}$, $X_{D}\supset D_{N}$, and that is
convex. Namely, for any two states in $X_{D}$, $x,x^{\prime}\in X_{D}$,

$\xi=\alpha x+(1-\alpha)x^{\prime}\in X_{D}$ (20)

where $0\leq\alpha\leq 1$.

Figure 1: Illustration of the dynamic range of a dataset defined by the convex
hull containing all points in the set, partitioned using a triangulation
method. The data points are in black and the convex hull that encapsulates all
data points is in grey.

Each data point $x_{i}$ is mapped to $f(x_{i})$, following the state
transition law in eq.(1). We assume that the transferred state, too, stays
within the same dynamic range $X_{D}$. Collecting all the transferred states
yields the second data set.

$\displaystyle D_{N}^{f}=\\{f(x_{i})\mid i=1,\cdots,N;x_{i}\in D_{N}\\}$ (21)
$\displaystyle D_{N}^{f}\subset X_{D}$ (22)

This implies that the state space of the nonlinear system under consideration
is closed within the dynamic range $X_{D}$.

With this dynamic range, we redefine our objective to compute the inner
products over $X_{D}$.

$\displaystyle R_{ij}=\int_{X_{D}}G_{ij}(x)dx$ (23) $\displaystyle
Q_{ij}=\int_{X_{D}}F_{ij}(x)dx$ (24)

The integrals can be computed by partitioning the domain $X_{D}$ into many
segments $X_{1},\cdots X_{P}$, as shown in Fig. 1.

$X_{D}=\bigcup_{p=1}^{P}X_{p}$ (25)

We generate these segments by applying a meshing technique to the data set
$D_{N}$, where the $n$-dimensional coordinates of individual data points are
treated as nodes of a mesh. Delaunay Triangulation, for example, generates a
triangular mesh structure with desirable properties [23]. As illustrated in
Fig. 1, each triangular element is convex and has no internal node. The volume
of the dynamic range $V(X_{D})$ is the sum of the volumes of all the elements.

$V(X_{D})=\sum_{p=1}^{P}\Delta v_{p}$ (26)

Accordingly, the integral $R_{ij}$ in eq.(23) can be segmented to

$R_{ij}=\sum_{p=1}^{P}\int_{X_{p}}G_{ij}(x)dx$ (27)

Suppose that the $p$-th element has $K_{p}$ nodes, as shown in Fig.2.
Renumbering these nodes 1 through $K_{p}$,

$\\{x[k_{p}]\Big{|}x[k_{p}]\in X_{p};k_{p}=1,\cdots,K_{p}\\},p=1,\cdots,P$
(28)

The integrand $G_{ij}$ within the $p$-th element can be approximated to the
mean of the $K_{p}$ nodes involved in the $p$-th element.

$G_{ij}(x;p)\approx\bar{G}_{ij,p}=\frac{1}{K_{p}}\sum_{k_{p}=1}^{K_{p}}G_{ij}(x[k_{p}]),x[k_{p}]\in
X_{p}$ (29)

If Delaunay Triangulation is used, $K_{p}=n+1$. See Fig. 2. Substituting this
into (23) yields the approximate value of $R_{ij}$.

$\hat{R}_{ij}=\sum_{p=1}^{P}\frac{1}{K_{p}}\sum_{k_{p}=1}^{K_{p}}G_{ij}(x[k_{p}])\Delta
v_{p}$ (30)

where

$\Delta v_{p}=\int_{X_{p}}1\cdot dx$ (31)

Similarly, each component of the matrix $Q$ can be computed by using the same
meshing.

$\hat{Q}_{ij}=\sum_{p=1}^{P}\frac{1}{K_{p}}\sum_{k_{p}=1}^{K_{p}}F_{ij}(x[k_{p}])\Delta
v_{p}$ (32)

Note that $F_{ij}$ is evaluated by using the data points in both $D_{N}^{f}$
and $D_{N}$,

$F_{ij}(x[k_{p}])=g_{i}[f(x[k_{p}])]\bar{g}_{j}(x[k_{p}])$ (33)

where $f(x[k_{p}])\in D_{N}^{f}$.

Figure 2: Visualization of the integrand calculation process. The volume of
the partition is encapsulated by the data points is denoted as $\Delta v_{p}$.
With this grouping, the value of $G_{ij}$ is calculated for each point and the
average among this group is computed, $\bar{G}_{ij,p}$. In turn, this value,
weighted by the volume of this partition, is summed across other partitions
(not shown) to approximate the value of the element $R_{ij}$.

###  III-B Convergence

Consider the center of each partition, $\bar{x}_{p}=\int_{X_{p}}xdx/\Delta
v_{p}$, and the distance between $\bar{x}_{p}$ and each point, $x[k_{p}]$:

$\Delta x[k_{p}]=\bar{x}_{p}-x[k_{p}]$ (34)

See Fig.2. The maximum distance from the center of the partition to each point
that makes up the partition is

$|\Delta x_{p}|=\max\\{|\Delta x[1]|,\,...\,|\Delta x[k_{p}-1]|,\,|\Delta
x[k_{p}]|\\}$ (35)

Consider a sequence of refining the approximate inner product integral
$\hat{R}_{ij}$ by increasing data points $N$. We can show that, as the number
of partition $P$ tends infinity and the maximum subintervals $|\Delta x_{p}|$
approach zero, the approximate inner product integral $\hat{R}_{ij}$ converges
to its true integral.

$R_{ij}=\lim_{\begin{subarray}{c}P\to\infty\\\ |\Delta x_{p}|\to
0\end{subarray}}\
\sum_{p=1}^{P}\frac{1}{K_{p}}\sum_{k_{p}=1}^{K_{p}}G_{ij}(x[k_{p}])\Delta
v_{p}$ (36)

This formulation takes the form of weighted sums, specifically Riemann sums.
Given functions that are bounded and continuous over the subdomain of
interest, sequences of this form are known to have a common limit and thus
converge upon refinement to the Riemann integral value over that subdomain,
according to Numerical Integration theory [22, Section 1.5].

###  III-C Algorithm

In the prior section, integrals (30) and (32) are presented as summations over
partitions. This computation can be streamlined by converting the summations
over partitions to the one over nodes. Consider node 3 associated to data
point $x_{3}$ in Fig.1, for example. This node is an apex of the 5 surrounding
triangles. This implies that integrand $G_{ij}(x)$ is calculated or recalled 5
times in computing (30) and (32). This repetition can be eliminated by
computing volume $\Delta v_{k}$ associated to node $k$, rather than partition
$p:\Delta v_{p}$. Namely, we compute

$\Delta v_{k}=\sum_{p=1}^{P}\frac{\Delta v_{p}}{K_{p}}I(k,p)$ (37)

where $I(k,p)$ is a membership function that takes value 1 when node $k$ is an
apex of triangle $p$, that is, node $k$ is involved in partition $p$. Using
this volume as a new weight we can rewrite (30) and (32) to be

$\displaystyle\hat{R}_{ij}=\sum_{k=1}^{K}G_{ij}(x[k])\Delta v_{k}$ (38)
$\displaystyle\hat{Q}_{ij}=\sum_{k=1}^{K}F_{ij}(x[k])\Delta v_{k}$ (39)

Using this conversion, the computation can be streamlined and cleanly
separated into three steps, as shown by pseudo-code in Algorithm 1. The steps
are: (1) Graph Creation: data are connected to create partitions of the domain
using a mesh generator: lines 3 to 8, (2) Weighting Calculation: calculation
of the weights for each data point: lines 10 to 17, and (3) Matrix
Calculation: the calculation of the $R$ and $Q$ matrices to find the matrix
$A$, lines 19 to 21.

Input:

1 $D_{N}$, $D_{N}^{f}$

Output:

2 $P$: $p$, $K$: $k$, $R$, $Q$ $A$

3 Graph Creation:

4 for  _$x_{i}$ in $D_{N}$ and corresponding $f(x_{i})$ in $D_{N}^{f}$_ do

5 Create node, $k$

6 Assign node attributes for current data point $k.x_{t}=x_{t}$ and
$k.x_{t+1}=x_{t+1}$

7 Append node to node list $K$.

8 end for

9then

10 Weighting Calculation:

11 Create list of triangles, $P$ using Delaunay Triangulation on node list $K$
using $k.{x_{t}}$

12 for  _$p$ in $P$_ do

13 Calculate volume, $V$, of $p$

14 for  _$k$ corresponding to $K_{p}$_ do

15 Update volume of each node $k.\Delta v_{k}=k.\Delta v_{k}+\frac{V}{n+1}$

16 end for

17

18 end for

19then

20 Matrix Calculation:

21 Find $R$ and $Q$ via (38) and (39)

Find $A=QR^{-1}$

Algorithm 1 Algorithm for DDE in pseudocode

##  IV Experiments

In this section, the DDE algorithm is implemented for the sake of evaluating
its validity and comparing its modeling accuracy to EDMD. Consider a 2nd order
nonlinear system consisting of a pendulum with a nonlinear damper. See Fig. 3.
The pendulum also bounces against walls with nonlinear compliance. The state
variables for this system are $x=[\theta,\dot{\theta}]^{T}$, and the equation
of motion can be written as:

$\ddot{\theta}=-sin(\theta)+F_{k}+F_{c}$ (40)

where $F_{k}$ and $F_{c}$ are wall reaction moment and damping moment,
respectively,

$\displaystyle F_{k}$
$\displaystyle=\begin{cases}-\text{sign}(\theta)\>k(|\theta|-\frac{\pi}{4})^{2}&\text{if
}|\theta|\geq\frac{\pi}{4}\\\ 0&\text{if }|\theta|<\frac{\pi}{4}\end{cases}$
(41) $\displaystyle F_{c}$
$\displaystyle=-\text{sign}(\dot{\theta})\>c\>\dot{\theta}^{2}$ (42)

where $k=200$ and $c=1$. We present a two part numerical experiment for this
system. The first experiment regards variations in dataset size and
distribution, and the second experiment varies the usage of observable
functions.

Figure 3: Diagram of pendulum with with walls. (a) depicts the range of the
pendulum where the walls are equally angularly displaced from the vertical.
(b) depicts the forces exerted on the pendulum due to contact with walls. (c)
depicts the damping force exerted onto the pendulum.

###  IV-A Dataset Variations

The datasets tested are of three types:

  1. 1.

Uniform: These datasets are composed of a rectangular dynamic range which is
sampled uniformly, like an evenly divided grid. The range varies from
$\theta=[-0.8,0.8]$ and $\dot{\theta}=[-2,2]$, where the mass can hit the
walls and the damping can vary from 0 to a significant value. See Fig. 3-(b),
(c).

  2. 2.

Gaussian: Data points are sampled with a finite-support Gaussian distribution.
The data are distributed non-uniformly with their highest density at the peak
of the Gaussian placed at diverse locations. In addition, each dataset
contains 100 data points uniformly distributed along the boarder of the
dynamic range to guarantee the same dynamic range as the uniform datasets.
Samples outside the dynamic range are excluded.

  3. 3.

Trajectories: These datasets are composed of trajectories, beginning from 100
initial conditions that are simulated forward the same number of time steps.
The dynamic range of this dataset differs from the two other dataset types.

The models constructed for DDE and EDMD use the same observable functions. The
observable functions chosen are two dimensional radial basis functions (RBFs),
uniformly distributed between the maximum and minimum values of each state
variable in their respective dataset, and the state variables. The total order
of the system is 27th order with 25 RBFs and 2 state variables.

A trajectory dataset graph is generated using Delaunay Triangles in DDE, shown
in Fig. 4.

Figure 4: Graph of connections for a trajectory dataset composed of 10000 points formed when utilizing the data-driven direct encoding method. The lightly red shaded region denotes the dynamic range. In the zoomed-in image, the difference in volumes associated to different data points can be observed. Noting the difference in size of the purple triangle versus the green triangle. TABLE I: Sum of Squared Errors over dynamic range with varying dataset sizes. Dataset Size | Total SSE | SSE Variance
---|---|---
EDMD / DDE | EDMD / DDE
Uniform Datasets
900 | 19.470 / 17.167 |  0.0094 / 0.0097
2500 | 17.995 / 16.471 |  0.0095 / 0.0103
10000 | 17.010 / 16.184 |  0.0099 / 0.0105
22500 | 16.698 / 16.133 |  0.0101 / 0.0105
Trajectory Datasets
1000 | 56.532 / 33.392 | 0.0349 / 0.0193
2500 | 33.330 / 25.064 | 0.0200 / 0.0130
5000 | 31.690 / 25.099 | 0.0195 / 0.0129
10000 | 30.184 / 25.101 | 0.0186 / 0.0129
25000 | 29.380 / 25.106 | 0.0181 / 0.0129

The accuracy of the models is tested through calculating sum of squared errors
(SSE) for one-step ahead predictions over the dynamic range of the datasets.
These error values are calculated for a uniform grid of points, similar to
that used in the uniform datasets. A visualization of the SSE is plotted in
Fig. 5. The results of these calculations are shown in Table I and II. In the
computation, the dynamic range is discretized, and the SSE value of each point
is summed. For the Gaussian datasets, the test is run for eight iterations of
each dataset to account for randomness and the average result is noted.

Figure 5: Sum of Squared Error plots for various datasets. (a) and (b) are EDMD models for Gaussian datasets; (a) uses a dataset that is centered at $[0,0]$ and (b) uses a dataset that is centered at $[0,2]$. (c) DDE and (d) EDMD models using the trajectory dataset composed of 2500 data points and using 25 RBFs. Circled in red are the regions with greatest variation in SSE between models. Only the dynamic range is shown in all plots. TABLE II: Sum of Squared Errors over dynamic range for Gaussian distributions with different centers. Centers of the distribution are noted above each column. Dataset Size | Total SSE
---|---
Center: $[0,0]$ | Center: $[0.8,0]$ | Center: $[0,2]$
EDMD / DDE | EDMD / DDE | EDMD / DDE
Gaussian Datasets
1000 | 25.565 / 24.377 | 25.161 / 22.98 | 24.338 / 22.445
2500 | 24.991 / 23.739 | 25.421 / 22.747 | 24.221 / 21.590
5000 | 24.662 / 23.518 | 26.805 / 22.415 | 23.914 / 21.195
10000 | 24.395 / 23.085 | 28.424 / 21.825 | 25.770 / 21.052
25000 | 25.219 / 22.167 | 30.788 / 21.687 | 26.941 / 20.704
TABLE III: Sum of Squared Errors over dynamic range with varying order of lifted linear models. # Observables | Total SSE
---|---
EDMD | DDE
Trajectory Dataset, 5000 points
27 | 31.690 | 25.099
51 | 36.657 | 21.637
83 | 28.437 | 13.613

###  IV-B Observable Function Variation

The second experiment varies the number of observable functions selected, thus
increasing the order. In this experiment, the number of RBFs is varied through
uniformly increasing the density of the centers of the function over the range
of the dataset. The results are noted in Table III. In the table, the number
detailing number of observables is the number of functions including the state
variables.

###  IV-C Discussion

From these results we can make the following observations.

  * •

All the numerical experiments show that DDE outperforms EDMD in total SSE.

  * •

For uniform datasets, both models are nearly equivalent, though DDE has
slightly lower SSE in all cases, shown in Table I. This result is expected as
all data points are weighted equally in a uniform distribution.

  * •

EDMD models exhibit significantly different distributions of prediction error,
depending on dataset distribution. In Fig. 5-(a), the EDMD model learned from
a dataset with high density near the origin produces a prediction error
distribution that is low in accuracy in the top-right and the bottom-left
corners of the dynamic range. In contrast, Fig.5-(b) illustrates that when
using EDMD to learn from data with high density at $\theta=0,\dot{\theta}=2$,
the model results in low accuracy in the lower half of the dynamic range. DDE
does not exhibit this high distribution dependency and achieves lower total
SSE, as shown in Table II.

  * •

In the trajectory datasets in Table I, DDE is not only lower in total SSE than
EDMD but is also smaller in variance. Fig. 5-(c) shows that DDE has a
uniformly low error distribution across the dynamic range, while EDMD in
Fig.5-(d) has two regions, as circled in the figure, with significantly larger
error. These regions coincide with sparsity in the dataset, providing evidence
of EDMD’s bias towards regions of high data density and explaining the
difference in performance between models for the non-uniform datasets.

  * •

In the trajectory datasets, the total SSE converges for DDE beginning with
small dataset sizes. This result implies that the elements in the $R$ and $Q$
matrices of DDE, that is, the inner product integral computations, converge as
the data size and the data density increase. This convergence is confirmed in
Fig. 6, where several elements of the $Q$ matrix are plotted against the data
size.

  * •

The second experiment, regarding variations in observable function numbers
demonstrates the effect of DDE remains even for significant increases in
observable functions, referring to Table III. In all cases tested, DDE
significantly outperforms EDMD over the dynamic range, as expected for a non-
uniform dataset.

Figure 6: Change in values of diagonal elements in the $Q$ matrix of DDE with
respect to trajectory dataset size. The specific elements shown correspond to
a state variable, $Q_{0,0}$, and three RBF functions of varying distances away
from the equilibrium.

##  V Application to Neural Net Koopman Modeling

The use of deep neural networks for finding effective observable functions and
constructing a Koopman linear model has been reported by several groups [10,
11, 12]. This method, sometimes referred to as Deep Koopman, is effective for
approximating the Koopman operator to a low-order model, compared to the use
of locally activated functions, such as RBFs, which scale poorly for high-
order nonlinear systems. The proposed DDE method can be incorporated into Deep
Koopman, further improving approximation accuracy.

Fig. 7 shows the architecture of the neural network similar to the prior works
[10, 11, 12]. The input layer receives training data of independent state
variables. The successive hidden layers produce observable functions; these
functions feed into the output layer consisting of linear activation
functions. This linear output layer corresponds to the $A$ matrix that maps
the observables of the current time to those of the next time step, i.e. the
state transition in the lifted space. In the Deep Koopman approach, the output
layer, that is, the $A$ matrix, is trained together with observable functions
in the hidden layers. Using the observable functions learned from deep
learning, the output layer weights are replaced with the $A$ matrix obtained
from DDE, captioned in Fig. 7.

Figure 7: Feedforward neural network model used to generate observable
functions. The final layer of the neural network is a linear layer. In the
standard Deep Koopman model this remains the same after the model is fully
trained. However, with the DDE model, this final layer is recalculated using
DDE by taking in the dataset used to train the model.

The effectiveness of this Deep Koopman-DDE method is applied to a multi-cable
manipulation system [6]. A simplified single cable version is utilized
consisting of six independent state variables. The network is constructed
using PyTorch, trained with an Adam optimizer, to generate 40 observable
functions. These model concatenates the state variables of the nonlinear
system, resulting in a 46th order model. Table IV shows parameters used for
training the Deep Koopman model. The training dataset is composed of 3000 data
points drawn from trajectories. Table V compares the Deep Koopman model to the
proposed model that uses DDE. Results are in terms of sum of squared error
over a set of test trajectories. A significant improvement is achieved by
incorporating DDE into the deep learning method.

TABLE IV: Neural network parameters and characteristics. Parameters and Characteristics | Value
---|---
Number of Hidden layers | 3
Activation Functions, Both Hidden Layers | ReLU
Width of 1st Hidden Layer | 16
Width of 2nd Hidden Layer | 16
Width of 3rd Hidden Layer | 40
Learning Rate, $\alpha$ | 0.01
TABLE V: Average SSE prediction error for trajectories in set of test data for single winch system. Modeling Method | 1 Time Step | 20 Time Steps
---|---|---
Deep Koopman only | 0.2471 | 9.8610
Deep Koopman + DDE | 0.2350 | 4.1131

Practical concerns arose from the use of Delaunay Triangulation in this
experiment. Specifically, the method was found to be limiting and inconsistent
for systems that are eighth order or higher, due to computational cost and
numerical instability. Alternative numerical integration approaches or
alternative partitioning methods may be necessary when applying the method to
higher order nonlinear systems.

##  VI Conclusion

In this work, a new data-driven approach to generating a Koopman linear model
based on the direct encoding of Koopman Operator (DDE) is presented as an
alternative to dynamic mode decomposition (DMD) and other related methods
using least squares estimate (LSE). The major contributions include: 1) The
analytical formula of Direct Encoding is converted to a numerical formula for
computing the inner product integrals from given data; 2) An efficient
algorithm is developed and its convergence conditions to the true results are
analyzed; 3) Numerical experiments demonstrate a) greater accuracy compared to
EDMD, b) lower sensitivity to data distribution, and c) rapid convergence of
inner product computation. Furthermore, the DDE method is incorporated to Deep
Koopman, i.e. neural network based methods for construction of the Koopman
Operator, for improving prediction accuracy.

## References

  * [1] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” _Journal of fluid mechanics_ , vol. 656, pp. 5–28, 2010.
  * [2] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data–driven approximation of the koopman operator: Extending dynamic mode decomposition,” _Journal of Nonlinear Science_ , vol. 25, no. 6, pp. 1307–1346, 2015.
  * [3] B. O. Koopman, “Hamiltonian systems and transformation in hilbert space,” _Proceedings of the National Academy of Sciences_ , vol. 17, no. 5, pp. 315–318, 1931.
  * [4] J. L. Proctor, S. L. Brunton, and J. N. Kutz, “Dynamic mode decomposition with control,” _SIAM Journal on Applied Dynamical Systems_ , vol. 15, no. 1, pp. 142–161, 2016.
  * [5] M. Korda and I. Mezić, “Linear predictors for nonlinear dynamical systems: Koopman operator meets model predictive control,” _Automatica_ , vol. 93, pp. 149–160, jul 2018. [Online]. Available: https://doi.org/10.1016%2Fj.automatica.2018.03.046
  * [6] J. Ng and H. H. Asada, “Model predictive control and transfer learning of hybrid systems using lifting linearization applied to cable suspension systems,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 682–689, 2021.
  * [7] V. Cibulka, T. Haniš, M. Korda, and M. Hromčík, “Model predictive control of a vehicle using koopman operator,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 4228–4233, 2020, 21st IFAC World Congress. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2405896320331748
  * [8] H. M. Calderón, E. Schulz, T. Oehlschlägel, and H. Werner, “Koopman operator-based model predictive control with recursive online update,” in _2021 European Control Conference (ECC)_ , 2021, pp. 1543–1549.
  * [9] D. Bruder, X. Fu, R. B. Gillespie, C. D. Remy, and R. Vasudevan, “Data-driven control of soft robots using koopman operator theory,” _IEEE Transactions on Robotics_ , vol. 37, no. 3, pp. 948–961, 2021.
  * [10] B. Lusch, J. N. Kutz, and S. L. Brunton, “Deep learning for universal linear embeddings of nonlinear dynamics,” _Nature Communications_ , vol. 9, no. 1, Nov 2018. [Online]. Available: http://dx.doi.org/10.1038/s41467-018-07210-0
  * [11] E. Yeung, S. Kundu, and N. Hodas, “Learning deep neural network representations for koopman operators of nonlinear dynamical systems,” 2017. [Online]. Available: https://arxiv.org/abs/1708.06850
  * [12] N. S. Selby and H. H. Asada, “Learning of causal observable functions for koopman-DFL lifting linearization of nonlinear controlled systems and its application to excavation automation,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 4, pp. 6297–6304, oct 2021. [Online]. Available: https://doi.org/10.1109%2Flra.2021.3092256
  * [13] M. Korda and I. Mezić, “Optimal construction of koopman eigenfunctions for prediction and control,” 2020.
  * [14] F. Fan, B. Yi, D. Rye, G. Shi, and I. R. Manchester, “Learning stable koopman embeddings,” 2021. [Online]. Available: https://arxiv.org/abs/2110.06509
  * [15] G. Mamakoukas, I. Abraham, and T. D. Murphey, “Learning stable models for prediction and control,” _IEEE Trans Robot_ , 2020.
  * [16] P. Bevanda, M. Beier, S. Kerz, A. Lederer, S. Sosnowski, and S. Hirche, “Diffeomorphically learning stable koopman operators,” _IEEE Control Systems Letters_ , vol. 6, pp. 3427–3432, 2022.
  * [17] M. Han, J. Euler-Rolle, and R. K. Katzschmann, “Desko: Stability-assured robust control with a deep stochastic koopman operator,” in _International Conference on Learning Representations_ , 2021.
  * [18] A. H. Abolmasoumi, M. Netto, and L. Mili, “Robust dynamic mode decomposition,” _IEEE Access_ , vol. 10, pp. 65 473–65 484, 2022. [Online]. Available: https://doi.org/10.1109%2Faccess.2022.3183760
  * [19] H. Asada, “Global, unified representation of heterogenous robot dynamics using composition operators,” 2022. [Online]. Available: https://arxiv.org/abs/2208.00175
  * [20] L. Ljung, “System identification,” in _Signal analysis and prediction_. Springer, 1998, pp. 163–173.
  * [21] D. A. Freedman, _Statistical models: theory and practice_. Cambridge University Press, 2009.
  * [22] P. J. Davis and P. Rabinowitz, _Methods of numerical integration_. Courier Corporation, 2007.
  * [23] S.-W. Cheng, T. K. Dey, J. Shewchuk, and S. Sahni, _Delaunay mesh generation_. CRC Press Boca Raton, 2013\.