aainstitutetext: Department of Physics and Astronomy, 9500 Gilman Drive,
University of California, San Diegobbinstitutetext: Department of Physics,
Faculty of Science, Chulalongkorn University, Bangkok, Thailand

# Narrow Resonances Revisited - Simplifying Multidimensional Constraints

R. Sekhar Chivukula b Pawin Ittisamai a James Osborne a and Elizabeth H.
Simmons<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

As we amass more LHC data, we continue to search for new and improved methods
of visualizing search results, in ways that are as model-independent as
possible. The simplified limits framework is an approach developed to recast
limits on searches for narrow resonances in terms of products of branching
ratios (BRs) corresponding to the resonance’s production and decay modes. In
this work, we extend the simplified limits framework to a multidimensional
parameter space of BRs, which can be used to unfold an ambiguity in the
simplified parameter $\zeta$ introduced when more than one channel contributes
to the production of the resonance. It is also naturally applicable to
combining constraints from experimental searches with different observed final
states. Constraints can be visualized in a three-dimensional space of
branching ratios by employing ternary diagrams, triangle plots which utilize
the inherent unitarity of the sum of the resonance’s BRs. To demonstrate this
new methodology, we recast constraints from recent ATLAS searches in diboson
final states for spin-0, 1, and 2 narrow resonances into constraints on the
resonance’s width-to-mass ratio and display them in the space of relevant
branching ratios. We also demonstrate how to generalize the method to cases
where more than three branching ratios are relevant by using N-simplex
diagrams, and we suggest a broader application of the general method to
digital data sets.

##  1 Introduction

Narrow resonance searches have long been a staple of experimental efforts to
identify or constrain new physics beyond the standard model (BSM). Typically,
the invariant mass of the system is measured and limits are placed on the
production cross section times branching ratio,
$\sigma_{\textrm{prod}}\times\textrm{BR}$, of new resonances and interpreted
within the context of specific benchmark (BM) models. In refs.
Chivukula:2016hvp ; Chivukula:2017lyk , the authors explored recasting these
model-independent constraints in terms of a simplified parameter, $\zeta$,
which depends only on the product of BRs and the ratio of the resonance width
to its mass. This reparameterization of narrow resonance limits in terms of
partonic quantities can often simplify the task of interpreting constraints in
terms of many models of interest.

Combining narrow resonance searches from multiple channels can extend the
exclusion limits of current and future searches for BSM physics with more than
one dominant decay mode. This possibility has been explored by both ATLAS
Aaboud:2018bun  and CMS Sirunyan:2019vgt  in diboson and dilepton channels at
the LHC. The constraints from different channels were combined in the context
of specific BM models where the relationships between BRs are known, a
necessity in order to cast constraints in terms of
$\sigma_{\textrm{prod}}\times\textrm{BR}$. This, however, introduces an
additional level of model dependence to the results, compared to monochannel
searches, which makes it difficult to apply such constraints to models whose
BRs vary from the BM choices.

In this work, we extend the concept of simplified limits, demonstrating a
methodology by which one may display constraints on the masses or widths of
narrow resonances within the parameter space of up to three BRs without
introducing new sources of model dependence. This method can be applied to
single channel searches where the production of a BSM resonance is predicted
to receive contributions from multiple channels. It can also be readily
applied to combined constraints on BSM scenarios with a single production mode
and several decay modes.

In sec. 2, we review the foundations of simplified limits for a single
experimental search. In sec. 3, we introduce ternary diagrams as a method of
unfolding constraints coming from multiple production modes or final states.
In sec. 4, we apply these methods to experimental searches for spin-0, spin-1,
and spin-2 resonances. Section 5 discusses the use of N-simplex diagrams as a
generalization of our method to cases where more than three branching ratios
are sources of experimental constraints on the properties of new resonances.
We conclude in sec. 6 and highlight the value of using the digital data record
to enable exploration of an arbirtrary-dimensional parameter space of
branching ratios.

##  2 Simplified Limits on Narrow Resonances

In this section, we review the formulation, presented in Refs.
Chivukula:2016hvp ; Chivukula:2017lyk , leading to the model-independent
simplified parameter $\zeta$ allowing one to compare bounds from data with the
easily calculated product of BRs corresponding to production and decay of a
narrow resonance as well as its width-to-mass ratio. We begin with the
resonance’s partonic cross section.

The tree-level cross section for resonant production of a state $R$ from
initial state partons $i,\,j$ and decaying to final state $x,\,y$ can be
written in the Breit-Wigner form as

$\displaystyle\hat{\sigma}_{ij\rightarrow R\rightarrow
xy}(\hat{s})=16\pi\mathcal{N}_{ij}\left(1+\delta_{ij}\right)\frac{\Gamma_{xy}\,\Gamma_{ij}}{(\hat{s}-m_{R}^{2})^{2}+m_{R}^{2}\Gamma_{R}^{2}}\,,$
(1)

with $\Gamma_{ab}\equiv\Gamma(R\rightarrow ab)$ the resonance’s partial decay
width, $\Gamma_{R}$ its total width, and $m_{R}$ its mass. Here, $\hat{s}$ is
the partonic center of mass energy of the system, while $\mathcal{N}_{ij}$ is
a ratio of spin and color factors,

$\displaystyle\mathcal{N}_{ij}\equiv\frac{N_{S_{R}}}{N_{S_{i}}N_{S_{j}}}\,\frac{C_{R}}{C_{i}C_{j}}\,,$
(2)

with $N_{S}$ and $C$, respectively, counting the number of spin and color
states of the incoming partons and resonance. In the narrow width
approximation (NWA), $\Gamma_{R}/m_{R}\ll 1$ and

$\displaystyle\frac{1}{(\hat{s}-m_{R}^{2})^{2}+m_{R}^{2}\Gamma_{R}^{2}}\approx\frac{\pi}{m_{R}\Gamma_{R}}\delta(\hat{s}-m_{R}^{2})\,.$
(3)

Thus, the tree-level cross section in the NWA is given by

$\displaystyle\hat{\sigma}_{ij\rightarrow R\rightarrow
xy}(\hat{s})=16\pi^{2}\mathcal{N}_{ij}\left(1+\delta_{ij}\right)\textrm{BR}_{xy}\,\textrm{BR}_{ij}\,\frac{\Gamma_{R}}{m_{R}}\,\delta(\hat{s}-m_{R}^{2})\,,$
(4)

with $\textrm{BR}_{ab}\equiv\Gamma_{ab}/\Gamma_{R}$ the BR of the resonance.

For hadron colliders, the partonic cross section is related to the
experimentally observable cross section by convolving it with the hadrons’
parton distribution functions (PDFs). For proton-proton colliders like the
LHC, we have

$\displaystyle\sigma_{pp\rightarrow R\rightarrow
xy}(s)=16\pi^{2}\mathcal{N}_{ij}\left(1+\delta_{ij}\right)\textrm{BR}_{ij}\,\textrm{BR}_{xy}\,\frac{\Gamma_{R}}{m_{R}}\left[\frac{1}{s}\frac{dL_{ij}}{d\tau}\right]_{\tau=\frac{m_{R}^{2}}{s}}\,,$
(5)

with $s$ the proton-proton center of mass energy. Here, $dL_{ij}/d\tau$ is the
parton luminosity function,

$\displaystyle\frac{dL_{ij}}{d\tau}\equiv\frac{1}{1+\delta_{ij}}\int_{\tau}^{1}\frac{dx}{x}\left[f_{i}(x,\,\mu_{F}^{2})f_{j}\left(\frac{\tau}{x},\,\mu_{F}^{2}\right)+f_{j}(x,\,\mu_{F}^{2})f_{i}\left(\frac{\tau}{x},\,\mu_{F}^{2}\right)\right]\,,$
(6)

with $f_{i}$ the PDF for parton $i$, $x$ the fraction of the proton’s momentum
carried by the parton, and $\mu_{F}$ the factorization scale. If multiple
partons contribute to the same experimental signal (e.g. light quark
production or decay), this can be extended to a sum over initial and/or final
state partons,

$\displaystyle\sigma_{pp\rightarrow R\rightarrow
XY}(s)=16\pi^{2}\sum_{i^{\prime}j^{\prime}}\textrm{BR}_{i^{\prime}j^{\prime}}\sum_{xy\in
XY}\textrm{BR}_{xy}\,\frac{\Gamma_{R}}{m_{R}}\left[\sum_{ij}\omega_{ij}\mathcal{N}_{ij}\frac{1+\delta_{ij}}{s}\frac{dL_{ij}}{d\tau}\right]_{\tau=\frac{m_{R}^{2}}{s}}\,,$
(7)

with $XY$ the observable final state. The weight function $\omega_{ij}$,

$\displaystyle\omega_{ij}\equiv\frac{\textrm{BR}_{ij}}{\sum_{i^{\prime}j^{\prime}}\textrm{BR}_{i^{\prime}j^{\prime}}}\,,$
(8)

lies between 0 and 1 such that $\sum\omega_{ij}=1$ by construction;
$\omega_{ij}$ represents the fraction of the resonance’s total production rate
due to each individual partonic channel.

We are now ready to define the $\zeta$ parameter,111The arrangement of
kroenecker deltas in this definition of the $\zeta$ parameter differs slightly
from that of Refs. Chivukula:2016hvp ; Chivukula:2017lyk . This definition of
$\zeta$ is more convenient because it has two distinct upper limits instead of
four, depending only on whether or not the production modes also contribute to
the observed final state.

$\displaystyle\zeta\equiv\sum_{ij}\textrm{BR}_{ij}\sum_{xy\in
XY}\textrm{BR}_{xy}\,\frac{\Gamma_{R}}{m_{R}}=\frac{\sigma_{pp\rightarrow
R\rightarrow
XY}(s)}{16\pi^{2}}\left[\sum_{ij}\omega_{ij}\mathcal{N}_{ij}\frac{1+\delta_{ij}}{s}\frac{dL_{ij}}{d\tau}\right]^{-1}_{\tau=\frac{m_{R}^{2}}{s}}\,.$
(9)

Note that $\zeta$ retains much of the model-independence of the experimental
search, which depends predominantly on the spin and helicity of the resonance,
while translating the constraint into purely partonic quantities. As defined
in eq. (9), $\zeta$ is bounded from above by $\Gamma_{R}/m_{R}$
($\Gamma_{R}/4m_{R}$) when there is overlap (no overlap) between initial and
final states. As we are working in the NWA where $\Gamma_{R}/m_{R}\lesssim
10\%$, this corresponds to an upper bound on $\zeta$ of $1/10$ ($1/40$) when
there is overlap (no overlap) between initial and final states.

##  3 Ternary Diagrams

In this section, we introduce a prescription for extending the simplified
limits framework to situations which necessitate considering constraints on
more than a one dimensional parameter. We begin by discussing the situations
considered in this paper: the limitations of the original simplified limits
formulation and combining constraints from multiple independent experimental
final states.

The simplified parameter $\zeta$ defined in eq. 9 offers a relatively model-
independent conversion from the limits on
$\sigma_{\textrm{prod}}\times\textrm{BR}$ to partonic quantities which more
directly describe the properties of the resonance: its mass, width, and BRs.
There are, however, several situations that strain the applicability of these
one-dimensional limits, which we will explore here. The first exception is in
situations where there may be more than one production mode. In these
situations, deconvolving the hadronic PDFs introduces an ambiguity in
constraints due to the weight factors $\omega_{ij}$, which are unknown a
priori without the introduction of additional model-dependent assumptions. An
example of this would be Drell-Yan (DY) production of a heavy $Z^{\prime}$,
which leads to a band in limits on $\zeta$. See the right panel of fig. 5 in
sec. 4.2 for an example of this scenario.

The second scenario we consider involves combining searches from multiple
experimentally distinguishable final states. Combining the results of multiple
searches offers increased statistics and is particularly effective when the
experiments in question are similarly sensitive to two or more distinguishable
final states. While monochannel searches place constraints on
$\sigma_{\textrm{prod}}\times\textrm{BR}$, doing so for combined searches
requires making an additional model-dependent assumption about the
relationship between the various channels’ BRs. The question we therefore wish
to address is how best to represent model-independent constraints within the
multidimensional space of parameters necessitated by the above two situations.

In the framework of simplified limits, the relevant model-independent
quantities are the resonance’s BRs, its mass, and its total width. In
combining channels, a model-independent scenario would be one where
constraints on the properties of the resonance could be displayed within the
space of BRs. Consider first the scenario of a resonance with one production
mode and two decay modes. Here, we are naively forced to place limits within a
three-dimensional space of BRs. However, we also have the sum rule

$\displaystyle\sum_{i=1}^{3}\textrm{BR}_{i}=1\,,$ (10)

reducing the space of independent BRs by one. Ternary diagrams,222Perhaps the
most familiar example of a ternary diagram in particle physics is the Dalitz
plot for three-body decays Dalitz:1953cp , where the sum of the two-body
invariant masses of the final state is constrained by the kinematics of the
system, e.g. for particle 0 decaying into particles 1, 2, and 3, the
constraint is $m_{12}^{2}+m_{13}^{2}+m_{23}^{2}=\sum_{i=0}^{3}m_{i}^{2}$.
representations of the space of three variables which sum to a constant, are
ideally suited to displaying constraints within this parameter space.

Fig. 1 shows an example of a ternary diagram, where at each point in the
diagram the sum of BRs is one. The tick marks on the axes are skewed to
indicate the lines of constant $\textrm{BR}_{i}$. Lines of constant
$\textrm{BR}_{2}$ run parallel to the $\textrm{BR}_{1}$ axis, lines of
constant $\textrm{BR}_{3}$ are parallel to the $\textrm{BR}_{2}$ axis, and
lines of constant $\textrm{BR}_{1}$ are parallel to the $\textrm{BR}_{3}$
axis. The point $P$ is labeled as an example, corresponding to
$\\{\textrm{BR}_{1},\,\textrm{BR}_{2},\,\textrm{BR}_{3}\\}=\\{0.5,\,0.4,\,0.1\\}$
with $\sum_{i=1}^{3}\textrm{BR}_{i}=1$ by construction. Within this parameter
space, one may then plot contours of upper limits on the remaining model
parameters. In the simplified limits framework, one may therefore use ternary
diagrams to display constraints on $m_{R}$ for fixed values of
$\Gamma_{R}/m_{R}$ or on $\Gamma_{R}/m_{R}$ for fixed values of $m_{R}$. We
present examples of ternary diagrams, displaying upper limits on
$\Gamma_{R}/m_{R}$, in sec. 4.

Figure 1: Example of a ternary diagram representing the space of BRs, with
each point within the diagram satisfying the constraint
$\sum_{i=1}^{3}\textrm{BR}_{i}=1$. The dashed lines represent hypothetical
experimental constraints from combining limits, showing contours of constant
$m_{R}$ or $\Gamma_{R}/m_{R}$. The point $P$ corresponds to
$\\{\textrm{BR}_{1},\,\textrm{BR}_{2},\,\textrm{BR}_{3}\\}=\\{0.5,\,0.4,\,0.1\\}$.

Above, we mentioned the scenario where one production mode and two decay modes
saturates all possible BRs. Interpreting the ternary diagrams for models of
resonances with more than three BRs turns out to be straightforward. For a
resonance with $n$ BRs, we instead have
$\sum_{i=1}^{3}\textrm{BR}_{i}=1-\sum_{i=4}^{n}\textrm{BR}_{i}$. We can thus
define “effective” BRs,

$\displaystyle\widetilde{\textrm{BR}}_{i}\equiv\textrm{BR}_{i}\left(1-\sum_{j=4}^{n}\textrm{BR}_{j}\right)^{-1}\,,$
(11)

which satisfy the unitary sum rule,
$\sum_{i=1}^{3}\widetilde{\textrm{BR}}_{i}=1$, implicit in the construction of
ternary diagrams. To interpret constraints displayed on the ternary diagrams
in the framework of simplified limits, one must then also define an effective
width,

$\displaystyle\widetilde{\Gamma}_{R}\equiv\Gamma_{R}\left(1-\sum_{j=4}^{n}\textrm{BR}_{j}\right)^{2}\,.$
(12)

A ternary diagram with sides spanning the range $[0,\,1]$ in this context
generically displays the space of effective BRs, where
$\widetilde{\textrm{BR}}\geq\textrm{BR}$ and
$\widetilde{\Gamma}_{R}\leq\Gamma_{R}$ with the equalities saturated only when
the resonance has exactly three decay modes. In many cases it is expected that
only a few channels will have similar experimental sensitivity over a wide
range of the parameter space, and ternary diagrams then capture the most
interesting region of the possible decay modes. However, a discussion of the
generalization of this method to N-simplexes is presented in sec. 5.

With this framework in place, we have introduced the means of extending
searches for narrow resonances into the multidimensional parameter space of
BRs. Within the simplified limits framework, this allows us to unfold the
ambiguity introduced by deconvolving the hadronic PDFs when there is more than
one production mode. In the case of multiple experimentally distinguishable
final states, it allows us to combine search results without introducing
further model-dependent assumptions about the relationship between BRs. In
what follows, we explore specific examples of each of these situations.

##  4 Applications

Searches for resonances in the diboson final state have been an important part
of the search for BSM physics at the LHC, having been studied in detail by
both ATLAS Aad:2020tps ; Aad:2020ddw ; Aad:2019fbh  and CMS Sirunyan:2019jbg .
Models usually considered for such searches include composite and little Higgs
models Marzocca:2012zn ; Bellazzini:2014yua ; Schmaltz:2005ky , heavy vector
triplet (HVT) models deBlas:2012qp ; Pappadopulo:2014qza , and models of
gravity in warped extra dimensions Randall:1999ee ; Agashe:2007zd . The
efficiency of a given search depends primarily on the spin and helicity of the
resonance while remaining relatively model-independent otherwise. In this
section, we explore the implications of searches for spin-0, 1, and 2 narrow
resonances, using diboson final states as examples. We apply the simplified
limits framework discussed above to existing searches by the ATLAS
collaboration.

For the conversion of experimental constraints and the calculation of the
simplified limits parameter $\zeta$ we use the CT18 NLO central PDF set
Hou:2019efy . The RS radion BM is calculated using leading-order formulae
given in ref. Barger:2011qn . Possible large higher-order corrections to the
radion’s $gg$ coupling, parameterized by a K-factor, and heavy quark couplings
are not included. Such corrections can be comparable to those of heavy Higgs
production via gluon fusion where $\textrm{K}>1$, leading to conservative
estimates of the experimental constraints presented here. The HVT BMs are
calculated at leading-order using CalcHEP 3.8.7 Belyaev:2012qa  with the model
file provided by ref. Pappadopulo:2014qza . The RS graviton BM is calculated
at leading-order using formulae given in refs. Agashe:2007zd ; Bijnens:2001gh
.

Combining the statistics of multiple BSM physics searches requires delicate
attention to the details of each experiment. In what follows, we do not
attempt to reproduce a complete statistical analysis of the experiments
considered. Instead, we make the conservative assumption that the combination
of constraints from multiple channels is given by the strongest of the bounds,

$\displaystyle\sigma_{\textrm{prod}}^{95}=\textrm{min}\\!\left[\sigma_{1}^{95},\,\sigma_{2}^{95},\,\dots\right]\,.$
(13)

A detailed combination of the statistics of each search will in general
produce stronger results, however such an analysis is beyond the scope of this
paper. For more thorough discussions of the statistics involved in such
searches, see e.g. refs. Read:2002hq ; Cowan:2010js .

###  4.1 Spin-0 Resonance

We first consider a neutral spin-0 resonance ($\phi$) produced via gluon
fusion. In ref. Aad:2020ddw , ATLAS reports constraints on the production of a
Randall-Sundrum (RS) radion Goldberger:1999uk  in the combined $WW+ZZ$
channels, with the ratio $\textrm{BR}(\phi\rightarrow
WW)/\textrm{BR}(\phi\rightarrow ZZ)$ fixed by the model. The neutral scalar
radion is a feature of extra-dimensional models which is predicted to
stabilize the size of the extra dimension. The radion coupling to SM fields is
inversely proportional to the vacuum expectation value of the radion field,
$\Lambda_{\phi}$ and proportional to the mass (mass squared) of the SM
fermions (bosons) it couples to. As the light fermion couplings to the radion
are suppressed by their masses, the radion is predominantly produced via gluon
fusion. We use $\Lambda_{\phi}=3$ TeV and $kL=35$ as a BM, where $kL$ is the
size of the extra dimension. The radion’s BRs are roughly constant for masses
above a few TeV, with sizable BRs of

$\displaystyle\begin{array}[]{l l
l}\textrm{BR}_{WW}=47\%\,,&\textrm{BR}_{hh}=23\%\,,&\textrm{BR}_{gg}=5.8\%\,,\\\
\textrm{BR}_{ZZ}=23\%\,,&\textrm{BR}_{t\overline{t}}=0.9\%\,,&\end{array}$
(16)

at $M_{\phi}=3$ TeV.

Figure 2: Constraints on narrow scalar resonance production from gluon fusion
in the $WW$ and $ZZ$ channels Aad:2020ddw . Constraints from ATLAS are
converted to upper limits on the $\zeta$ parameter, eq. (9), represented by
solid red and blue lines for the $WW$ and $ZZ$ channels, respectively. The
gray shaded region corresponds to the upper limit on the product of branching
ratios times $\Gamma_{\phi}/M_{\phi}=10\%$, approximately where the NWA breaks
down. The dot-dashed lines show the predictions for the radion in our BM
model.

Fig. 2 shows the ATLAS constraints from the individual channels converted into
the language of simplified limits. Generally the constraints from each channel
are quite competitive, except below $M_{\phi}\sim 1$ TeV and between
approximately 2.0 – 2.6 TeV where the $ZZ$ channel is significantly more
constraining. Above $M_{\phi}\sim 4.5$ TeV, the search does not constrain any
model which satisfies the NWA assumption. Also shown are the predictions for
our BM radion, which sets limits on the mass of the radion of $M_{\phi}\gtrsim
3.3~{}(2.9)$ TeV in the $WW$ ($ZZ$) channel. Other models would be represented
by different curves in $\zeta$ vs. $M_{\phi}$, corresponding to different
limits on the mass of the resonance.

Figure 3: Ternary diagram showing constraints on a narrow scalar resonance
produced via gluon fusion and decaying to $WW$ and $ZZ$ Aad:2020ddw  for a
resonance mass of 3.3 TeV. The colors depict constraints on
$\widetilde{\Gamma}_{\phi}/M_{\phi}$ (labeled right of the legend), which have
been translated to the corresponding constraint on $\Lambda_{\phi}$ for the BM
value $kL=35$ (labeled left of the legend). The black dot represents the
location in the parameter space of BRs corresponding to the radion in our BM
model, and the black line through the legend labels the value of the
experimental constraint on the BM radion. The BM prediction
$\widetilde{\Gamma_{\phi}}/M_{\phi}\sim 10^{-1.5}$ (corresponding to
$\Lambda_{\phi}=3$ TeV), being slightly above the the experimental limit,
demonstrates that the model is just slightly excluded.

Applying the principles of sec. 3, fig. 3 shows constraints on
$\widetilde{\Gamma}_{\phi}/M_{\phi}$ for $M_{\phi}=3.3$ TeV, which is near the
experimental limit for our BM radion. The radion BM values are labeled by a
point in the plane of $\widetilde{\textrm{BR}}$s and by a line on the legend
of $\widetilde{\Gamma}_{\phi}/M_{\phi}$. In general, one can see that as the
BRs for either the production mode or both decay modes decrease, the
constraints unsurprisingly also weaken. Conversely, constraints are strongest
when $\widetilde{\textrm{BR}}(\phi\rightarrow
gg)\sim\widetilde{\textrm{BR}}(\phi\rightarrow VV)$ for a single channel, with
the other channel’s BR negligible. Contours of constant
$\log_{10}\widetilde{\Gamma}_{\phi}/M_{\phi}$ are shown in increments of
$10^{-1}$, and the BM value $\widetilde{\Gamma}_{\phi}/M_{\phi}\sim 10^{-1.5}$
roughly corresponds to the experimental limit, which is labeled by a solid
black line through the legend. Sharp edges in the contours of constant
$\widetilde{\Gamma}_{\phi}/M_{\phi}$ occur when the constraints from both
channels are equal, and highlight the region of parameter space where one
would expect to gain the most from a proper statistical combination of the
independent searches. As the BM radion model also predicts a sizable
$\textrm{BR}(\phi\rightarrow hh)$, we also see the application of effective
BRs, where $\widetilde{\textrm{BR}}_{WW}\;(\textrm{BR}_{WW})=62\%\;(47\%)$,
$\widetilde{\textrm{BR}}_{ZZ}\;(\textrm{BR}_{ZZ})=31\%\;(23\%)$, and
$\widetilde{\textrm{BR}}_{gg}\;(\textrm{BR}_{gg})=7.7\%\;(5.8\%)$.

Figure 4: Ternary diagrams showing constraints on scalar narrow resonance
production from gluon fusion in the $WW$ and $ZZ$ channels Aad:2020ddw  for a
selection of resonance masses. The colors depict constraints on
$\widetilde{\Gamma}_{\phi}/M_{\phi}$ (labeled right of the legends), which
have been translated to the corresponding constraint on $\Lambda_{\phi}$ for
the BM value $kL=35$ (labeled left of the legends). The black dot represents
the location in the parameter space corresponding to the radion in our BM
model, and the black line through the legend labels the value of the
experimental constraint on the BM radion.

To demonstrate how the diagrams change as a function of mass, fig. 4 shows
constraints on $\widetilde{\Gamma}_{\phi}/M_{\phi}$ for several choices of
resonance mass between $1.0~{}\textrm{TeV}\leq M_{\phi}\leq
4.0~{}\textrm{TeV}$. Here one can clearly see how the constraints weaken as
the mass of the resonance increases, while within the space of
$\widetilde{\textrm{BR}}$s the relative trends remain similar to those
discussed above for fig. 3. Because the radion couplings are all proportional
to $\Lambda_{\phi}^{-1}$, increasing $\Lambda_{\phi}$ decreases the radion’s
total width and vice versa, while keeping individual BRs constant at leading
order. Therefore, each point on the ternary diagram also translates directly
into a constraint on $\Lambda_{\phi}$ for a fixed value of $kL$. Constraints
on $\Lambda_{\phi}$ for the BM value $kL=35$ are also labeled to the left of
the legend of each diagram, and the experimental constraint on the BM point is
labeled by a solid black line through the legend. For the BM value
$\Lambda_{\phi}=3$ TeV, the radion is shown to be excluded in the diagrams
corresponding to $M_{\phi}\leq 3.0$ TeV, while it is unconstrained in the
$M_{\phi}=4.0$ TeV diagram.

###  4.2 Spin-1 Resonance

For a spin-1 resonance, we consider the production and decay of either a
charged ($W^{\prime\,\pm}$) or neutral ($Z^{\prime}$) narrow vector resonance.
Here gluon fusion production is forbidden by charge conservation or Yang’s
theorem, so production occurs via either Drell-Yan (DY) or vector boson fusion
(VBF) processes. For DY production of a charged vector resonance, production
occurs via $q\overline{q}^{\prime}$ annihilation and in most cases is
dominated by the valence quark combination $u\overline{d}$ or $d\overline{u}$.
For a neutral vector resonance, however, production occurs via $q\overline{q}$
and is dominated by a combination of valence quarks,
$u\overline{u}+d\overline{d}$. As the ratio of couplings of the neutral
resonance to up and down quarks is not known a priori, deconvolving the proton
PDFs introduces an ambiguity to the $\zeta$ parameter, as discussed in sec. 3.

The production and decay of $W^{\prime\,\pm}$ and $Z^{\prime}$ can be
parameterized in terms of an HVT model, which is a phenomenological framework
proposed by the authors of ref. Pappadopulo:2014qza  to cover a variety of
explicit BSM models. In the HVT model, an $SU(2)_{L}$ vector triplet
$V^{\prime}$ is introduced, and its interactions with the SM fields are
parameterized by a variety of couplings. The parameter $g_{V}$ characterizes
the typical interaction strength of the new triplet, and the parameters
$c_{H}$ and $c_{F}$ characterize deviations from this strength in coupling to
the Higgs and fermion currents, respectively. A factor of $g^{2}/g_{V}^{2}$ is
inserted in the coupling of $V^{\prime}$ to the SM fermions to make contact
with many specific extended gauge models in the literature, where $g$ denotes
the $SU(2)_{L}$ gauge coupling. Therefore, the interaction of $V^{\prime}$
with the Higgs doublet current (and therefore with the longitudinal components
of the SM $W$ and $Z$) is parameterized by $(c_{H}g_{v})$. Likewise, the
coupling of $V^{\prime}$ to the SM fermions is controlled by the combination
$(g^{2}/g_{V})c_{f}$.

BM values are typically chosen to represent the range of various specific BSM
extensions. Model A, with $g_{V}=1$, is representative of a weakly coupled
scenario such as in theories with an extended gauge symmetry. Model B, with
$g_{V}=3$, is chosen to represent a strongly-coupled composite Higgs
scenario.333Our BM values for $c_{F}$ and $c_{H}$ are calculated from the
relations in appendix A of ref. Pappadopulo:2014qza  with
$\widetilde{c}_{VW}=-1$, $\widetilde{c}_{H}=\widetilde{c}_{F}=0$ for model A
and $\widetilde{c}_{VW}=-\widetilde{c}_{H}=1$, $\widetilde{c}_{F}=0$ for model
B. In both BM scenarios we assume a universal coupling of $V^{\prime}$ to all
SM fermions; all interactions not mentioned above, which only contribute to
the decay of $V^{\prime}$ via the small mixing between $V^{\prime}$ and the SM
weak gauge bosons, are turned off. The two BM models predict dramatically
different decays for the heavy spin-one resonance. Model A predicts a BR of a
few percent to bosons, with dijets representing the dominant decay mode, while
model B predicts decays predominantly into dibosons.

Figure 5: Constraints on charged (left pane) and neutral (right pane) narrow
vector resonances from DY production in the $WV$ Aad:2020ddw  and $VH$
Aad:2020tps  channels. Constraints from ATLAS are converted to upper limits on
the $\zeta$ parameter, eq. (9), represented by solid red and blue lines for
the $WV$ and $VH$ channels, respectively. The shaded bands in the constraints
of the right plot are due to the model-dependent relationship between the
production mode couplings of the $Z^{\prime}$, with the lower constraint
corresponding to purely $u\overline{u}$ production and the upper constraint
corresponding to purely $d\overline{d}$ production. The gray shaded region
corresponds to the upper limit on the product of branching ratios times
$\Gamma_{V^{\prime}}/M_{V^{\prime}}=10\%$, where the NWA is no longer valid.
Also shown are the HVT $pp\rightarrow V^{\prime}\rightarrow WV$ theory
predictions for BM models A (black dot-dashed line) and B (black dotted line).

We recast constraints from ATLAS narrow resonance searches in the $WV$
Aad:2020ddw  and $VH$ Aad:2020tps  channels in terms of the simplified limits
framework. The left plot of fig. 5 shows constraints on the production of a
narrow $W^{\prime}$ resonance in the $WZ$ and $WH$ channels, as well as the
$W^{\prime}\rightarrow WZ$ predictions for the HVT BM models A and B. Although
not shown, predictions for the $W^{\prime}\rightarrow WH$ channel are quite
similar. We see that the $WZ$ channel dominates the constraints, with a narrow
$W^{\prime}$ excluded for masses below 3.8 (4.0) TeV for model A (model B).
The right plot of fig. 5 shows constraints on the production of a narrow
$Z^{\prime}$ resonance in the $WW$ and $ZH$ channels, as well as the
$Z^{\prime}\rightarrow WW$ predictions for the HVT BM models A and B. Here,
the band in the constraints is due to the model-dependent relationship between
the up and down quark couplings to $Z^{\prime}$. The lower limit of $\zeta$
corresponds to $\omega_{u\overline{u}}=1$ and $\omega_{d\overline{d}}=0$,
while the upper limit of $\zeta$ corresponds to $\omega_{d\overline{d}}=1$ and
$\omega_{u\overline{u}}=0$ with other values of $\omega$ falling between these
two extremes. The HVT $Z^{\prime}$ is seen to be excluded for masses below
$3.1~{}(3.3)~{}\textrm{TeV}\lesssim M_{Z^{\prime}}\lesssim
3.7~{}(4.0)~{}\textrm{TeV}$ for model A (model B).

Figure 6: Ternary diagrams showing constraints on a charged vector narrow
resonance from DY production in the $WZ$ Aad:2020ddw  and $WH$ channels
Aad:2020tps . The colors depict constraints on
$\widetilde{\Gamma}_{W^{\prime}}/m_{W^{\prime}}$ for a selection of resonance
masses. Points labeled “HVT A” and “HVT B” mark the locations in the parameter
space corresponding to the BMs for model A and model B, respectively. The
black lines through the legend (labeled “EXP A” and “EXP B”) show the values
of the experimental constraints on the BM models, while the tick marks to the
left of the legend (labeled “HVT A” and “HVT B”) show their predicted values.

As the production of a charged vector resonance suffers no ambiguity in the
initial state, we may apply the same principles used in sec. 4.1 to combine
limits from the $WZ$ and $WH$ channels. Fig. 6 shows the combined constraints
on $\widetilde{\Gamma}_{W^{\prime}}/M_{W^{\prime}}$ for several choices of
resonance mass. At $M_{W^{\prime}}=1.0$ TeV, there are no constraints from the
$WH$ search, so the ternary diagram displays only constraints from the $WZ$
channel. On the other hand, the constraints at $M_{W^{\prime}}=2.0$ TeV are
similar for both channels which is reflected in the diagram. The remaining
diagrams are again dominated by the $WZ$ channel, in agreement with the left
panel of fig. 5. The HVT BMs are also displayed on the ternary diagrams, both
by points labeling their predicted branching ratios and by tick marks on the
left of legend labeling the predicted values of
$\log_{10}(\widetilde{\Gamma}_{W^{\prime}}/M_{W^{\prime}})$. The solid lines
through the legends label the experimental constraints on the BM points, so
that one can more easily see that both BMs are excluded for
$M_{W^{\prime}}\leq 3.0$ TeV, while model A is unconstrained at
$M_{W^{\prime}}=4.0$ TeV and model B is very close to the experimental
constraint.

Figure 7: Ternary diagrams showing constraints on a neutral vector narrow
resonance from DY production in the $WW$ Aad:2020ddw  and $ZH$ channels
Aad:2020tps . The colors depict constraints on
$\widetilde{\Gamma}_{Z^{\prime}}/m_{Z^{\prime}}$ for a selection of resonance
masses. Points labeled “HVT A” and “HVT B” mark the locations in the parameter
space corresponding to the BMs for model A and model B, respectively. The
black lines through the legend (labeled “EXP A” and “EXP B”) show the values
of the experimental constraints on the BM HVT models, while the tick marks to
the left of the legend (labeled “HVT A” and “HVT B”) show their predicted
values.

Conversely, while searches for a neutral vector resonance via DY production
may not be ideal candidates for model-independent combined searches, one may
instead probe its coupling to light quarks via ternary diagrams. We take for
example the search for a DY-produced $Z^{\prime}$ resonance in the $WW$
channel, which dominates the constraints when compared to the $ZH$ channel for
a majority of the HVT parameter space. The main production modes are via
either $u\overline{u}$ or $d\overline{d}$ annihilation, and fig. 7 shows
constraints on $\widetilde{\Gamma}_{Z^{\prime}}/M_{Z^{\prime}}$ over a range
of masses. The HVT BMs are also displayed on the ternary diagrams, both by
points labeling their predicted branching ratios and by tick marks on the left
of the legend labeling the predicted values of
$\log_{10}(\widetilde{\Gamma}_{Z^{\prime}}/M_{Z^{\prime}})$. The solid lines
through the legends label the experimental constraints on the BM points, so
that one can more easily see that both BMs are excluded for
$M_{Z^{\prime}}\leq 3.0$ TeV, while they are both unconstrained at
$M_{Z^{\prime}}=4.0$ TeV. In these examples, the tilt in the contours of
constant $\log_{10}(\widetilde{\Gamma}_{Z^{\prime}}/M_{Z^{\prime}})$ follows
from the roughly 2-to-1 luminosity ratio of $u$-to-$d$ quarks in the proton’s
PDF.444For a discussion of $Z^{\prime}$ properties characterized in part by
their couplings to quarks, see e.g. Ref. Carena:2004xs .

###  4.3 Spin-2 Resonance

For a spin-2 resonance, we consider the production of a narrow resonance via
gluon fusion in the $WW$ and $ZZ$ channels. Heavy spin-2 resonances are a
generic feature of models of quantum gravity in extra dimensions, where
Kaluza-Klein (KK) towers of heavy gravitons ($G_{\textrm{KK}}$) are predicted.
Typical BM models for considering the lightest graviton mode are the original
RS model, referred to as RS1 Randall:1999ee , where all SM fields are
localized on the IR brane, and the bulk RS model Agashe:2007zd , where the SM
fields are allowed to propagate in the bulk. In RS1, the localization of the
graviton in the warped bulk near the IR brane induces couplings to all SM
fields which are only TeV suppressed, so the graviton has significant BRs to
light fermions. On the other hand, in the bulk RS model light fields are
localized toward the UV brane which greatly suppresses their couplings to
gravitons. Instead, the graviton is mainly produced via gluon fusion or VBF
and has significant BRs to $t\overline{t}$, $WW$, and $ZZ$. In what follows we
will therefore consider the bulk RS model.

Figure 8: Constraints on spin-2 narrow resonance production from gluon fusion
in the $WW$ and $ZZ$ channels Aad:2020ddw . Constraints from ATLAS are
converted to upper limits on the $\zeta$ parameter, eq. (9), represented by
solid red and blue lines for the $WW$ and $ZZ$ channels, respectively. The
gray shaded region corresponds to the upper limit on the product of branching
ratios times $\Gamma_{\phi}/M_{\phi}=10\%$, approximately where the NWA breaks
down. The dot-dashed lines show the predictions from the Bulk RS BM model.

The graviton’s couplings to SM fields are determined by the parameter
$k/\overline{M}_{\textrm{Pl}}$ where $k$ is the warped curvature scale and
$\overline{M}_{\textrm{Pl}}$ is the reduced Planck mass. For our BM model of
graviton production, we assume $k/\overline{M}_{\textrm{Pl}}=1.0$. Fig. 8
shows constraints from ATLAS searches in the $WW$ and $ZZ$ channels
Aad:2020ddw , converted into the language of simplified limits. The
constraints from each channel are competitive except for the region below
$M_{G_{\textrm{KK}}}\lesssim 1.2$ TeV, where the $ZZ$ channel dominates. Also
shown is the prediction from the bulk RS BM model, which sets a limit on the
mass of the graviton of $M_{G_{\textrm{KK}}}\gtrsim 1.7~{}(1.5)$ TeV in the
$WW$ ($ZZ$) channel.

Figure 9: Ternary diagrams showing constraints on spin-2 narrow resonance
production from gluon fusion in the $WW$ and $ZZ$ channels Aad:2020ddw . The
colors depict constraints on $\widetilde{\Gamma}_{R}/m_{R}$ for a selection of
resonance masses. The black dot represents the location in the parameter space
corresponding to the bulk RS BM model. The black line through the legend
labeled “EXP” shows the value of the experimental constraints on the bulk RS
BM model, while the tick mark to the left of the legend labeled “RS” shows its
predicted value.

Fig. 9 shows ternary diagrams associated with these searches, displaying
constraints on $\widetilde{\Gamma}_{G_{\textrm{KK}}}/M_{G_{\textrm{KK}}}$ for
a selection of graviton masses. The predictions from the bulk RS BM are also
shown, both by a point labeling their predicted branching ratios and by a tick
mark on the left of the legend labeling the predicted value of
$\log_{10}(\widetilde{\Gamma}_{G_{\textrm{KK}}}/M_{G_{\textrm{KK}}})$. The
solid line through the legend labels the experimental constraint on the BM
point, so that one can more easily see that the BM model is excluded for
$M_{G_{\textrm{KK}}}\leq 1.5$ TeV and unconstrained for
$M_{G_{\textrm{KK}}}\geq 2.0$ TeV, in agreement with fig. 8.

##  5 Further Generalization: N-Simplexes

In this section, we indicate how to further extend the simplified limits
framework to situations in which more than a few branching ratios are
important sources of experimental constraints on new resonances.

In sec. 3, we introduced a means of extending searches for narrow resonances
into the multidimensional parameter space of BRs, focusing on situations where
only three BR’s were of greatest importance: either two production modes and
one primary decay mode, or one production mode and a pair of experimentally
distinguishable final states. Section 4 explored a number of specific
examples, showing how ternary diagrams can provide insight.

Let us now consider the situation where a new resonance has N+1 branching
ratios, all able to provide significant experimental constraints. Then we can
generalize the work of sec. 3, by noting that the branching ratios now obey
the sum rule:

$\displaystyle\sum_{i=1}^{N+1}\textrm{BR}_{i}=1\,,$ (17)

which tells us that the dimension of the space of independent BR’s is N. Just
as compositional data in three variables can be represented by a ternary plot,
where the ratios summing to 1 are plotted within a two-dimensional equilateral
triangle, so compositional data in N+1 variables can be represented in a
simplicial sample space where the ratios summing to 1 are located within an
N-simplex.555A simplex generalizes a triangle to arbitrary dimensions; it is
the minimal polytope in the space of a given number of dimensions. A 0-simplex
is a point; a 1-simplex is a line segment; a 2-simplex is a triangle; a
3-simplex is a tetrahedron; a 4-simplex is a 5-cell, and so forth.

Moreover, if the new resonance has M+1 branching ratios, of which only N+1 are
able to provide significant experimental constraints, we can generalize the
notion of ”effective” branching ratios from sec. 3 as well. Now, we have
$\sum_{i=1}^{N+1}\textrm{BR}_{i}=1-\sum_{i=N+2}^{M+1}\textrm{BR}_{i}$. We can
thus define “effective” BRs,

$\displaystyle\widetilde{\textrm{BR}}_{i}\equiv\textrm{BR}_{i}\left(1-\sum_{j=N+2}^{M+1}\textrm{BR}_{j}\right)^{-1}\,,$
(18)

which satisfy the unitary sum rule,
$\sum_{i=1}^{N+1}\widetilde{\textrm{BR}}_{i}=1$, implicit in the construction
of simplex diagrams. To interpret constraints displayed on the N-simplex
diagrams in the framework of simplified limits, one must then also define an
effective width,

$\displaystyle\widetilde{\Gamma}_{R}\equiv\Gamma_{R}\left(1-\sum_{j=N+2}^{M+1}\textrm{BR}_{j}\right)^{2}\,.$
(19)

An N-simplex diagram with sides spanning the range $[0,\,1]$ in this context
generically displays the space of effective BRs, where
$\widetilde{\textrm{BR}}\geq\textrm{BR}$ and
$\widetilde{\Gamma}_{R}\leq\Gamma_{R}$ with the equalities saturated only when
$N=M$.

Again, within the simplified limits framework, this construction allows us to
handle a wider variety of scenarios. On the one hand, we can unfold the
ambiguity introduced by deconvolving the hadronic PDFs when there is more than
one production mode. On the other hand, we can combine search results for
multiple experimentally distinguishable final states without introducing
further model-dependent assumptions about the relationship between BRs.

While one cannot easily plot a higher-dimensional N-simplex as a two-
dimensional image, one can nonetheless still perform a statistical analysis to
discern how the experimental constraints shape the allowed region of the
N-simplex and check whether a given new resonance’s location within the
N-simplex is in the allowed region. One might even illustrate this in a
journal article by displaying the ternary diagram sub-space of the full
N-simplex that provides the strongest limit on the model in question.

##  6 Discussion

In this article we have introduced a more integrative method of presenting
constraints on production and decay of narrow resonances when multiple
branching ratios yield valuable experimental information about the properties
of new resonances. The method utilizes the NWA to parameterize constraints in
terms of products of BRs, the mass of the resonance, and the total width of
the resonance. We have seen that representing the results of searches for
narrow resonances in terms of the parameterization of the simplified limits
framework—$\widetilde{\textrm{BR}}$s, $\widetilde{\Gamma}_{R}$, and
$m_{R}$—provides a natural context for combining the statistics from multiple
search channels for a common resonance.

We have largely focused on cases where only three channels (two production and
one decay mode, or vice versa) are relevant, by employing ternary diagrams to
display the combined constraints on the expanded space of BRs. We have
illustrated applications to resonances of spin 0, spin 1, and spin 2, arising
in a variety of beyond-the-standard-model scenarios. Our approach clearly
offers a more model-independent method of interpreting constraints from
multiple channels compared to the traditional product of production cross
section times BR, which would require making assumptions about the
relationship between decay BRs. It is also applicable to situations with
multiple production modes by unfolding the uncertainty inherent in the one-
dimensional simplified limits parameter $\zeta$. This method is complementary
to traditional limits, with $\sigma\times\textrm{BR}$ offering the cleanest
display of constraints at the expense of sometimes introducing specific model
assumptions, while ternary diagrams can encompass more a more model-
independent parameter space, easily translatable to a variety of disparate
models.

While the applications we considered in detail had three-channel parameter
spaces, we have also discussed how the use of ternary diagrams can be readily
generalized to the use of N-simplex diagrams for situations where additional
branching ratios also offer valuable experimental constraints on the
properties of new resonances. Making use of this more generalized method will,
of course, rely on the availability of data about the multiple branching
ratios. Indeed, having access to digital data sets for searches combining
multiple experimental channels would be ideal for documenting and leveraging
limits on a many-dimensional parameter space, capable of encompassing all of
the initial and final state BRs relevant for the searches being considered.
Recently, there has been a tremendous effort to make digitalized data
available to all researchers with the introduction of the HepData repository
Maguire:2017ypu .

As a closing thought, we would like to advocate for experimental
collaborations to provide these larger digital data sets (parameterized in
terms of the BRs, total width, and mass of the new resonance) to supplement
the information presentable in the traditional article format. This will
enable the data to be most fully leveraged to explore the widest possible
range of models in detail, enabling constraints to be quickly understood for a
plethora of interesting theories.

##  7 Acknowledgments

We thank J. Duarte for suggesting the application of this work to digital data
repositories. We also thank D. Foren, K.A. Mohan, D. Sengupta, and X. Wang for
useful discussions and comments. This material is based upon work supported by
the National Science Foundation under Grant No. PHY-1915147. P.I. was
supported by the CUniverse research promotion project of Chulalongkorn
University in Bangkok, Thailand, under Grant No. CUAASC.

## References

  * (1) R.S. Chivukula, P. Ittisamai, K. Mohan and E.H. Simmons, _Simplified Limits on Resonances at the LHC_ , _Phys. Rev. D_ 94 (2016) 094029 [1607.05525].
  * (2) R.S. Chivukula, P. Ittisamai, K. Mohan and E.H. Simmons, _Broadening the Reach of Simplified Limits on Resonances at the LHC_ , _Phys. Rev. D_ 96 (2017) 055043 [1707.01080].
  * (3) ATLAS collaboration, _Combination of searches for heavy resonances decaying into bosonic and leptonic final states using 36 fb -1 of proton-proton collision data at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Phys. Rev. D_ 98 (2018) 052008 [1808.02380].
  * (4) CMS collaboration, _Combination of CMS searches for heavy resonances decaying to pairs of bosons or leptons_ , _Phys. Lett. B_ 798 (2019) 134952 [1906.00057].
  * (5) R. Dalitz, _On the analysis of tau-meson data and the nature of the tau-meson_ , _Phil. Mag. Ser. 7_ 44 (1953) 1068.
  * (6) ATLAS collaboration, _Search for resonances decaying into a weak vector boson and a Higgs boson in the fully hadronic final state produced in proton $-$proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, 2007.05293.
  * (7) ATLAS collaboration, _Search for heavy diboson resonances in semileptonic final states in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, 2004.14636.
  * (8) ATLAS collaboration, _Search for diboson resonances in hadronic final states in 139 fb -1 of $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _JHEP_ 09 (2019) 091 [1906.08589].
  * (9) CMS collaboration, _A multi-dimensional search for new heavy resonances decaying to boosted WW, WZ, or ZZ boson pairs in the dijet final state at 13 TeV_ , _Eur. Phys. J. C_ 80 (2020) 237 [1906.05977].
  * (10) D. Marzocca, M. Serone and J. Shu, _General Composite Higgs Models_ , _JHEP_ 08 (2012) 013 [1205.0770].
  * (11) B. Bellazzini, C. Csáki and J. Serra, _Composite Higgses_ , _Eur. Phys. J. C_ 74 (2014) 2766 [1401.2457].
  * (12) M. Schmaltz and D. Tucker-Smith, _Little Higgs review_ , _Ann. Rev. Nucl. Part. Sci._ 55 (2005) 229 [hep-ph/0502182].
  * (13) J. de Blas, J. Lizana and M. Perez-Victoria, _Combining searches of Z’ and W’ bosons_ , _JHEP_ 01 (2013) 166 [1211.2229].
  * (14) D. Pappadopulo, A. Thamm, R. Torre and A. Wulzer, _Heavy Vector Triplets: Bridging Theory and Data_ , _JHEP_ 09 (2014) 060 [1402.4431].
  * (15) L. Randall and R. Sundrum, _A Large mass hierarchy from a small extra dimension_ , _Phys. Rev. Lett._ 83 (1999) 3370 [hep-ph/9905221].
  * (16) K. Agashe, H. Davoudiasl, G. Perez and A. Soni, _Warped Gravitons at the LHC and Beyond_ , _Phys. Rev. D_ 76 (2007) 036006 [hep-ph/0701186].
  * (17) T.-J. Hou et al., _New CTEQ global analysis of quantum chromodynamics with high-precision data from the LHC_ , 1912.10053.
  * (18) V. Barger and M. Ishida, _Randall-Sundrum Reality at the LHC_ , _Phys. Lett. B_ 709 (2012) 185 [1110.6452].
  * (19) A. Belyaev, N.D. Christensen and A. Pukhov, _CalcHEP 3.4 for collider physics within and beyond the Standard Model_ , _Comput. Phys. Commun._ 184 (2013) 1729 [1207.6082].
  * (20) J. Bijnens, P. Eerola, M. Maul, A. Mansson and T. Sjostrand, _QCD signatures of narrow graviton resonances in hadron colliders_ , _Phys. Lett. B_ 503 (2001) 341 [hep-ph/0101316].
  * (21) A.L. Read, _Presentation of search results: The CL(s) technique_ , _J. Phys. G_ 28 (2002) 2693.
  * (22) G. Cowan, K. Cranmer, E. Gross and O. Vitells, _Asymptotic formulae for likelihood-based tests of new physics_ , _Eur. Phys. J. C_ 71 (2011) 1554 [1007.1727].
  * (23) W.D. Goldberger and M.B. Wise, _Modulus stabilization with bulk fields_ , _Phys. Rev. Lett._ 83 (1999) 4922 [hep-ph/9907447].
  * (24) M. Carena, A. Daleo, B.A. Dobrescu and T.M.P. Tait, _$Z^{\prime}$ gauge bosons at the Tevatron_, _Phys. Rev. D_ 70 (2004) 093009 [hep-ph/0408098].
  * (25) E. Maguire, L. Heinrich and G. Watt, _HEPData: a repository for high energy physics data_ , _J. Phys. Conf. Ser._ 898 (2017) 102006 [1704.05473].


# Training neural networks with end-to-end optical backpropagation

James Spall  Clarendon Laboratory, University of Oxford, Parks Road, Oxford,
OX1 3PU, UK  Xianxin Guo<EMAIL_ADDRESS>Clarendon Laboratory,
University of Oxford, Parks Road, Oxford, OX1 3PU, UK  Lumai Ltd., Wood Centre
for Innovation, Quarry Road, Headington, Oxford, OX3 8SB, UK  A. I. Lvovsky
<EMAIL_ADDRESS>Clarendon Laboratory, University of Oxford,
Parks Road, Oxford, OX1 3PU, UK  Lumai Ltd., Wood Centre for Innovation,
Quarry Road, Headington, Oxford, OX3 8SB, UK

###### Abstract

Optics is an exciting route for the next generation of computing hardware for
machine learning, promising several orders of magnitude enhancement in both
computational speed and energy efficiency. However, to reach the full capacity
of an optical neural network it is necessary that the computing not only for
the inference, but also for the training be implemented optically. The primary
algorithm for training a neural network is backpropagation, in which the
calculation is performed in the order opposite to the information flow for
inference. While straightforward in a digital computer, optical implementation
of backpropagation has so far remained elusive, particularly because of the
conflicting requirements for the optical element that implements the nonlinear
activation function. In this work, we address this challenge for the first
time with a surprisingly simple and generic scheme. Saturable absorbers are
employed for the role of the activation units, and the required properties are
achieved through a pump-probe process, in which the forward propagating signal
acts as the pump and backward as the probe. Our approach is adaptable to
various analog platforms, materials, and network structures, and it
demonstrates the possibility of constructing neural networks entirely reliant
on analog optical processes for both training and inference tasks.

## Introduction

Machine learning, one of the most revolutionary scientific breakthroughs in
the past decades, has completely transformed the technology landscape,
enabling innovative applications in fields ranging from natural language
processing to drug discovery. As the demand for increasingly sophisticated
machine learning models continues to escalate, there is a pressing need for
faster and more energy-efficient computing solutions. In this context, analog
computing has emerged as a promising alternative to traditional digital
electronics [1, 2, 3, 4, 5, 6, 7]. A particularly exciting platform for analog
neural networks (NNs) is optics, in which the interference and diffraction of
light during propagation implements the linear part of every computational
layer [8, 9].

Most of the current analog computing research and development is aimed at
using the NN for inference [10, 8]. Training such NNs, on the other hand, is a
challenge. This is because the backpropagation algorithm [11], the workhorse
of training in digital NNs, requires the calculation to be performed in the
order opposite to the information flow for inference, which is difficult to
implement on an analog physical platform. Hence analog models are typically
trained offline (_in silico_), on a separate digital simulator, after which
the parameters are transferred to the analog hardware. In addition to being
slow and inefficient, this approach can lead to errors arising from imperfect
simulation and systematic errors (‘reality gap’). In optics, for example,
these effects may result from dust, aberrations, spurious reflections and
inaccurate calibration [12].

To enable learning in analog NNs, different approaches have been proposed and
realized [13]. Several groups explored various ‘hardware-in-the loop’ schemes,
in which, while the backpropagation was done _in silico_ , the signal acquired
from the analog NN operating in the inference regime was incorporated into the
calculation of the feedback for optimizing the NN parameters [14, 12, 15, 16,
17, 18, 19, 20]. This has partially reduced the training error, but has not
addressed the low speed and inefficiency of _in silico_ training.

Recently, several optical neural networks (ONNs) were reported that were
trained online (_in situ_) using methods alternative to backpropagation.
Bandyopadhyay _et al._ trained an ONN based on integrated photonic circuits
using simultaneous perturbation stochastic approximation, i.e. randomly
perturbing all ONN parameters and using the observed change of the loss
function to approximate its gradient [21]. Filipovich  _et al._ applied direct
feedback alignment, wherein the error calculated at the output of the ONN is
used to update the parameters of all layers [22]. However both these methods
are inferior to backpropagation as they take much longer to converge,
especially for sufficiently deep ONNs [23].

An optical implementation of the backpropagation algorithm was proposed by
Hughes _et al._ [24], and recently demonstrated experimentally [25], showing
that the training methods of current digital NNs can be applied to analog
hardware. However, their scheme omitted a crucial step for optical
implementation: backpropagation through nonlinear activation layers. Their
method requires digital nonlinear activation and multiple opto-electronic
inter-conversions inside the network, complicating the training process.
Furthermore, the method applies only to a specific type of ONN that uses
interferometer meshes for the linear layer, and does not generalise to other
ONN architectures. Complete implementation of the backpropagation algorithm in
optics, through all the linear and nonlinear layers, that can generalise to
many ONN systems, remains a highly challenging goal.

In this work, we address this long-standing challenge and present the first
complete optical implementation of the backpropagation algorithm in a two-
layer ONN. The gradients of the loss function with respect to the NN
parameters are calculated by light travelling through the system in the
reverse direction. The main difficulty of all-optical training lies in the
requirement that the nonlinear optical element used for the activation
function needs to exhibit different properties for the forward and backward
propagating signals. Fortunately, as demonstrated in our earlier theoretical
work [26] and explained below, there does exist a group of nonlinear
phenomena, which exhibits the required set of properties with sufficient
precision.

We optically train our ONNs to perform classification tasks, and our results
surpass those trained with a conventional in silico method. Our optical
training scheme can be further generalized to other platforms using different
linear layers and analog activation functions, making it an ideal tool for
exploring the vast potential of analog computing for training neural networks.

Figure 1: Illustration of optical training. (a) Network architecture of the
ONN used in this work, which consists of two fully-connected linear layers and
a hidden layer. (b) Simplified experimental schematic of the ONN. Each linear
layer performs optical MVM with a cylindrical lens and a spatial light
modulator (SLM) that encodes the weight matrix. Hidden layer activations are
computed using SA in an atomic vapour cell. Light propagates in both
directions during optical training. (c) Working principle of SA activation.
The forward beam (pump) is shown by solid red arrows, backward (probe) by
purple wavy arrows. The probe transmission depends on the strength of the pump
and approximates the gradient of the SA function. For high forward intensity
(top panel), a large portion of the atoms are excited to the upper level.
Stimulated emission produced by these atoms largely compensates the absorption
due to the atoms in the ground level. For weak pump (bottom panel), the
excited level population is low and the absorption is significant. (d) Neural
network training procedure. (e) Optical training procedure. Both signal and
error propagation in the two directions are fully implemented optically. Loss
function calculation and parameter update are left for electronics without
interrupting the optical information flow.

## Optical training algorithm

We consider a multilayer perceptron — a common type of NN which consists of
multiple linear layers that establish weighted connections between neurons,
interlaid by activation functions that enable the network to learn complex
nonlinear functions. To train the NN, one presents it with a training set of
labeled examples and iteratively adjusts the NN parameters (weights and
biases) to find the correct mapping between the inputs and outputs.

The training steps are summarised in Fig. 1(d). The weight matrices, denoted
$W^{(i)}$ for the $i$-th layer, are first initialized with random values. Each
iteration of training starts by entering the input examples from the training
set as input vectors $x=a^{(0)}$ into the NN, and forward propagating through
all of its layers. In every layer $i$, one performs a matrix-vector
multiplication (MVM) of the weight matrix and the activation vector,

$z^{(i)}=W^{(i)}\times a^{(i-1)},$ (1)

followed by element-wise application of the activation function $g(\cdot)$ to
the resulting vector:

$a^{(i)}=g\left(z^{(i)}\right).$ (2)

The output $y=a^{(L)}$ of an $L$-layer NN allows one to compute the loss
function $\mathcal{L}(y,t)$ that determines the difference between the network
predictions $y$ and ground truth labels $t$ from the training set. The
backpropagation algorithm helps calculating the gradient of this loss function
with respect to all the parameters in the network, through what is essentially
an application of the chain rule of calculus. The network parameters are then
updated using these gradients and optimization algorithms such as stochastic
gradient descent. The training process is repeated until convergence.

The gradients we require are given by [11] as

$\frac{\partial\mathcal{L}}{\partial W^{(i)}}=\delta^{(i)}\otimes a^{(i-1)},$
(3)

where $\delta^{(i)}$ is referred to as the “error vector” at the $i$th layer
and $\otimes$ denotes the outer product. The error vector is calculated as

$\displaystyle\delta^{(i-1)}$
$\displaystyle=\left({W^{(i)}}^{T}\times\delta^{(i)}\right)g^{\prime}\left(z^{(i-1)}\right),$
(4)

going through layers in reverse sequence. The expression for the error vector
$\delta^{(L)}$ in the last layer depends on the choice of the loss function,
but for the common loss functions of mean-squared error and cross-entropy
(with an appropriate choice of activation function) it is simply the
difference between the NN output and the label: $\delta^{(L)}=y-t$.

Therefore, to calculate the gradients at each layer we need one vector from
the forward pass through the network (the activations) and one vector from the
backward pass (the errors).

We see from Eq. (4) that the error backpropagation consists of two operations.
First we must perform an MVM, mirroring the feedforward linear operation (1).
In an ONN, this can be done by light that propagates backwards through the
same linear optical arrangement [27]. The second operation consists in
modulation of the MVM output by the activation function derivative and poses a
significant challenge for optical implementation. This is because most optical
media exhibit similar properties for forward and backward propagation. On the
other hand, our application requires an optical element that is (1) nonlinear
in the forward direction, (2) linear in the backward direction and (3)
modulates the backward light amplitude by the derivative of the forward
activation function.

We have solved this challenge with our optical backpropagation protocol, which
calculates the right-hand side of Eq. (4) entirely optically with no opto-
electronic conversion or digital processing. The first component of our
solution is the observation that many optical media exhibit nonlinear
properties for strong optical fields, but are approximately linear for weak
fields. Hence, we can satisfy the conditions (1) and (2) by maintaining the
back-injected beam at a much lower intensity level than the forward.
Furthermore, there exists a set of nonlinear phenomena that also addresses the
requirement (3). An example is saturable absorption (SA). The transmissivity
of SA medium in the backward direction turns out to approximate the derivative
of its intensity-dependent transmission in the forward direction
$g^{\prime}\left(z^{(i-1)}\right)$. This approximation is valid up to a
certain numerical factor and only for small values of $z^{(i-1)}$; however, as
shown in our prior work [26], this is sufficient for successful training.

## Multilayer ONN

Our ONN as shown in Fig. 1(a,b) is implemented in a free-space tabletop
setting. The neuron values are encoded in the transverse spatial structure of
the propagating light field amplitude. Spatial light modulators are used to
encode the input vectors and weight matrices. The NN consists of two fully-
connected linear layers implemented with optical MVM following our previously
demonstrated experimental design [28].

This design has a few characteristics that make it suitable for use in a deep
neural network. First, it is reconfigurable, so that both neuron values and
network weights can be arbitrarily changed. Second, multiple MVM blocks can be
cascaded to form a multilayer network, as the output of one MVM naturally
forms the input of the next MVM. Using a coherent beam also allows us to
encode both positive- and negative-valued weights. Finally, the MVM works in
both directions, meaning the inputs and outputs are reversible, which is
critical for the implementation of our optical backpropagation algorithm. The
hidden layer activation between the two layers is implemented optically by
means of SA in a rubidium atomic vapour cell [Fig. 1(c)].

## Results

Figure 2: Multi-layer ONN characterisation. (a) Scatter plots of measured
against theory results for MVM-1 (first layer forwards), MVM-2a (second layer
forwards) and MVM-2b (second layer backwards). All three MVM results are taken
simultaneously. Histograms of the signal and noise error for each MVM are
displayed underneath. (b) First-layer activations $a^{(1)}_{\rm meas}$
measured after the vapor cell, plotted against the theoretically expected
linear MVM-1 output $z^{(1)}_{\rm theory}$ before the cell. The green line is
a best fit curve of the theoretical SA nonlinear function. (c) The amplitude
of a weak constant probe passed backwards through the vapor cell as a function
of the pump $z^{(1)}_{\rm theory}$, with constant input probe. Measurements
for both forward and backward beams are taken simultaneously.

### Linear layers

We first set up the linear layers that serve as the backbone of our ONN, and
we make sure that they work accurately and simultaneously in both directions —
a highly challenging task that has never been achieved before to our best
knowledge.

This involves three MVMs: first layer in the forward direction (MVM-1), second
layer in both forward (MVM-2a) and backward (MVM-2b) directions. To
characterise these MVMs, we apply random vectors and matrices and
simultaneously measure the output of all three: the results for 300 random
MVMs are presented in Fig. 2(a). To quantify the MVM performance, we define
the signal-to-noise ratio (SNR, see Methods for details). As illustrated by
the histograms, MVM-1 has the greatest SNR of 14.9, and MVM-2a has a lower SNR
of 7.1, as a result of noise accumulation from both layers and the reduced
signal range. MVM-2b has a slightly lower SNR of 6.7, because the optical
system is optimized for the forward direction. Comparing these experimental
results with a simple numerical model, we estimate $1.3\%$ multiplicative
noise in our MVMs, which is small enough not to degrade the ONN performance
[12].

### Nonlinearity

With the linear layers fully characterized, we now measure the response of the
activation units in both directions. With the vapor cell placed in the setup
and the laser tuned to resonance with the atomic transition, we pass the
output of MVM-1 through the vapor cell in the forward direction. The response
as presented in Fig. 2(b) shows strong nonlinearity. We fit the data with the
theoretically expected SA transmissivity (see Supplementary for details),
thereby finding the optical depth to be $\alpha_{0}=7.3$, which is sufficient
to achieve high accuracy in ONNs [26]. The optical depth and the associated
nonlinearity can be easily tuned to fit different network requirements by
controlling the temperature of the vapor cell. In the backward direction, we
pass weak probe beams through the vapor cell and measure the output. Both the
forward and backward beams are simultaneously present in the vapor cell during
the measurement.

In Fig. 2(c) we measure the effect of the forward amplitude $z^{(1)}$ on the
transmission of the backward beam through the SA. The theoretical fit for
these data — the expected backward transmissivity calculated from the physical
properties of SA — is shown by the red curve. For comparison, the orange curve
shows the rescaled exact derivative $g^{\prime}\left(z^{(1)}\right)$ of the SA
function, which is the dependence required for the calculation (4) of the
training signal. Although the two curves are not identical, they both match
the experimental data for a broad range of neuron values generated from the
random MVM, hence the setting is appropriate for training.

### All optical classification

Figure 3: Optical training performance. (a) Decision boundary charts of the
ONN inference output for three different classification tasks, after the ONN
has been trained optically (top) or _in-silico_ (bottom). (b) Learning curves
of the ONN for classification of the ‘rings’ dataset, showing mean and
standard deviation of the validation loss and accuracy averaged over 5
repeated training runs. Shown above are decision boundary charts of the ONN
output for the test set, after different epochs. (c) Evolution of output
neuron values, and of output errors, for the training set inputs of the two
classes.

After setting up the two-layer ONN, we perform end-to-end optical training and
inference on classification tasks: distinguishing two classes of data points
on a two-dimensional plane (Fig. 3). We implement a fully-connected feed-
forward architecture, with three input neurons, five hidden layer neurons and
two output neurons (Fig. 1). Two input neurons are used to encode the input
data point coordinates $(x_{1},x_{2})$, and the third input neuron of constant
value is used to set the first layer bias. The class label is encoded by a
‘one-hot’ vector $(0,1)$ or $(1,0)$, and we use categorical cross-entropy as
the loss function.

We optically train the ONN on three 400-element datasets with different
nonlinear boundary shapes, which we refer to as ‘rings’, ‘XOR’ and ‘arches’
[Fig. 3(a)]. Another 200 similar elements of each set are used for validation,
i.e. to measure the loss and accuracy after each epoch of training. The test
set consists of a uniform grid of equally-spaced $(x_{1},x_{2})$ values. The
optical inference results for the test set are displayed in Fig. 3(a) by light
purple and orange backgrounds, whereas the blue circles and orange triangles
show the training set elements.

For all three datasets, each epoch consists of 20 mini-batches, with a mini-
batch size of 20, and we use the Adam optimizer to update the weights and
biases from the gradients. We tune hyperparamters such as learning rate and
number of epochs to maximise network performance. Table 1 summarises the
network architecture and hyperparameters used for each dataset.

Dataset |  | Input
---
neurons
| Hidden
---
neurons
| Output
---
neurons
| Learning
---
rate
Epochs |  | Batches
---
per epoch
Batch size
Rings | 2 | 5 | 2 | 0.01 | 16 | 20 | 20
XOR | 0.005 | 30
Arches | 0.01 | 25
Table 1: Summary of network architecture and hyperparameters used in optical
and digital training.

Figure 3(b) shows the optical training performance on the ‘rings’ dataset. We
perform five repeated training runs, and plot the loss and accuracy for the
validation set after each epoch of training. To visualise how the network is
learning the boundary between the two classes, we also run a test dataset
after each epoch. Examples of the network output after 1, 3, 6 and 10 epochs
are shown. We see that the ONN quickly learns the nonlinear boundary and
gradually improves the accuracy to $100\%$. This indicates a strong optical
nonlinearity in the system and a good gradient approximation in optical
backpropagation. Details of the training procedure are provided in the Methods
section, and results for the other two datasets in Supplementary Note 3.

To better understand the optical training process, we explore the evolution of
the output neuron and error vector values in Fig. 3(c). First, we plot the
mini-batch mean value of each output neuron, $\bar{z}_{j}^{(2)}$, for the
inputs from the two classes separately in the upper and lower panels, over the
course of the training iterations. We see the output neuron values diverge in
opposite ways for the two classes, such that the inputs can be distinguished
and correctly classified.

Second, we similarly plot the evolution of the mini-batch mean output error,
$\bar{\delta}^{(2)}$, for each neuron. This is calculated as the difference
between the network output vector $a^{(2)}$ and the ground truth label $y$,
averaged over each mini-batch. As expected, we see the output errors converge
towards zero as the system learns the correct boundary.

### Optical training vs in-silico training

To demonstrate the optical training advantage, we perform in-silico training
of our ONN as a comparison. We digitally model our system with a neural
network of the equivalent architecture, including identical learning rate,
number of epochs and all other hyperparamters. The hidden layer nonlinearity
and the associated gradient are given by the best fit curve and theoretical
probe response of Fig. 2(b). The trained weights are subsequently used for
inference with our ONN. The top and bottom rows in Fig. 3(a) plot the network
output of the test boundary set, after the system has been trained optically
and digitally, respectively, for all three datasets. In all cases, the
optically trained network achieves almost perfect accuracy, whilst the
digitally trained network is clearly not optimised, with the network
prediction not matching the data. This is further evidence for the already
well-documented advantages of hardware-in-the-loop training schemes.

## Discussion

Network layer | Function | Implementation examples
---|---|---
| MVM | Free-space optical multiplier, photonic crossbar array
Linear layer | Diffraction | Programmable optical mask
| Convolution | Lens Fourier transform
| Saturable absorption | Atomic vapour cell, semiconductor absorber, Graphene
Nonlinear layer | Saturable gain | EDFA, SOA, Raman amplifier
| Intensity-dependent phase modulation | Optical Kerr medium, thermal nonlinear material
Table 2: Generalization of the optical training scheme. EDFA: Erbium-doped
fiber amplifier. SOA: semiconductor optical amplifier.

Our optical training scheme is surprisingly simple and effective. It adds
minimal computational overhead to the network, since it doesn’t require a
priori simulation, or intricate mapping of network parameters to physical
device settings. It also imposes minimal hardware complexity on the system, as
it requires only a few additional beam splitters and detectors to measure the
activation and error values for parameter updates.

Our scheme can be generalised and applied to many other analog neural networks
with different physical implementations of the linear and nonlinear layers. We
list a few examples in Table. 2. Common optical linear operations include MVM,
diffraction and convolution. Compatible optical MVM examples include our free-
space multiplier and photonic crossbar array [29], as they are both
bidirectional, in the sense that optical field propagating backwards through
these arrangements gets multiplied by the transpose of the weight matrix.
Diffraction naturally works in both directions, hence diffractive neural
networks constructed using different programmable amplitude and phase masks
also satisfy the requirements [30]. Optical convolution, achieved with the
Fourier transform by means of a lens, and mean pooling, achieved through an
optical low pass filter, also work in both directions. Therefore, a
convolutional neural network can be optically trained as well. Detailed
analysis on the generalization to these linear layers can be found in
Supplementary Note 4.

Regarding the generalization to other nonlinearity choices, the critical
requirement is the ability to acquire gradients during backpropagation. Our
pump-probe method is compatible with multiple types of optical nonlinearities:
saturable absorption, saturable gain and intensity-dependent phase modulation
[31]. Using saturable gain as the nonlinearity offers the added advantage of
loss compensation in a deep network, and using intensity-dependent phase
modulation nonlinearity, such as self-lensing, allows one to build complex-
valued optical neural networks with potentially stronger learning capabilities
[32, 26, 33].

In our ONN training implementation, some computational operations remain
digital, specifically the calculation of the last layer error $\delta^{(2)}$
and the outer product (3) between the activation and error vectors. Both these
operations can be done optically if needed [26]. The error vector can in many
cases be obtained by subtracting the ONN output from the label vector by way
of destructive interference [12]. Interference can also be utilized to compute
the outer product by fanning out the two vectors and overlapping them in a
criss-cross fashion onto a pixelated photosensor.

Our optically trained ONN can be scaled up to improve computing performance.
Previously, using a similar experimental setup, we have demonstrated an ONN
with 100 neurons per layer and high multiplier accuracy [12], and 1000 neurons
can be supported by commercial SLMs or liquid crystal displays. Optical
encoders and detectors can work at speeds up to 100 GHz using off-the-shelf
components, enabling ultra-fast end-to-end optical computing at low latency.
Therefore computational speeds up to $10^{17}$ operations per second are
within reach, and our optical training method is compatible with this
productivity rate.

## Methods

### Multilayer ONN

To construct the multi-layer ONN, we connect two optical multipliers in
series. For the first layer (MVM-1), the input neuron vector $x$ is encoded
into the field amplitude of a coherent beam using a digital micromirror device
(DMD), DMD-1. This is a binary amplitude modulator, and every physical pixel
is a micromirror that can reflect at two angles representing 0 or 1. By
grouping 128 physical pixels as a block, we are able to represent 7-bit
positive-valued inputs on DMD-1, with the input value proportional to the
number of binary pixels ‘turned on’ in each block.

Since MVM requires performing dot products of the input vector with every row
of the matrix, we create multiple copies of the input vector on DMD-1, and
image them onto the $W^{(1)}$ matrix mask — a phase-only liquid-crystal
spatial light modulator (LC-SLM), SLM-1 — for element-wise multiplication. The
MVM-1 result $z^{(1)}$ is obtained by summing the element-wise products using
a cylindrical lens (first optical ‘fan-in’), and passing the beam through a
narrow adjustable slit to select the zero spatial frequency component. The
weights in our ONN are real-valued, encoded by LC-SLMs with 8-bit resolution
using a phase grating modulation method that enables arbitrary and accurate
field control [34].

The beam next passes through a rubidium vapor cell to apply the activation
function, such that immediately after the cell the beam encodes the hidden
layer activation vector, $a^{(1)}$. The beam continues to propagate and
becomes the input for the second linear layer. Another cylindrical lens is
used to expand the beam (first optical ‘fan-out’), before modulation by the
second weight matrix mask SLM-2. Finally, summation by a third cylindrical
lens (second optical ‘fan-in’) completes the second MVM in the forward
direction (MVM-2a), and the final beam profile encodes $z^{(2)}$.

To read out the activation vectors required for the optical training, we
insert beam splitters at the output of each MVM to tap-off a small portion of
the beam. The real-valued vectors are measured by high-speed cameras, using
coherent detection techniques detailed in Supplementary Note 2.

At the output layer of the ONN we use a digital softmax function to convert
the output values into probabilities, and calculate the loss function and
output error vector, which initiates the optical backpropagation.

### Optical backpropagation

The output error vector, $\delta^{(2)}$ is encoded in the backward beam by
using DMD-2 to modulate a beam obtained from the same laser as the forward
propagating beam. The backward beam is introduced to the system through one of
the arms of the beam splitter placed at the output of MVM-2a, and carefully
aligned so as to overlap with the forward beam. SLM-2 performs element-wise
multiplication by the transpose of the second weight matrix.

The cylindrical lens that performs ‘fan-out’ for the forward beam, performs
‘fan-in’ for the backward beam into a slit, completing the second layer
backwards MVM (MVM-2b). Passing through the vapor cell modulates the beam by
the derivative of the activation function, after which the beam encodes the
hidden layer error vector $\delta^{(1)}$. Another beam splitter and camera are
used to tap off the backward beam and measure the result.

In practice, two halves of the same DMD act as DMD-1 and DMD-2, and a portion
of SLM-1 is used to encode the sign of the error vector. A full experiment
diagram is provided in Supplementary Note 1.

Each training iteration consists of optically measuring all of $a^{(1)}$,
$z^{(2)}$ and $\delta^{(1)}$. These vectors are used, along with the inputs
$x=a^{(0)}$, to calculate the weight gradients according to Eq. (3) and weight
updates, which are then applied to the LC-SLMs. This process is repeated for
all the mini-batches until the network converges.

### SA activation

The cell with atomic rubidium vapor is heated to 70 degrees by a simple
heating jacket and temperature controller. The laser wavelength is locked to
the $D_{2}$ transition at $780$nm.

The power of the forward propagating beam is adjusted to ensure the beam at
the vapor cell is intense enough to saturate the absorption, whilst the
maximum power of the backward propagating beam is attenuated to approximately
$2\%$ of the maximum forward beam power, to ensure a linear response when
passing through the cell in the presence of a uniform pump.

In the experiment, the backward probe response does not match perfectly with
the simple two-level atomic model, due to two factors.

First, the probe does not undergo $100\%$ absorption even with the pump turnes
off. Second, a strong pump beam causes the atoms to fluoresce in all
directions, including along the backward probe path. Therefore, the backward
signal has a background offset proportional to the forward signal. To
compensate for these issues, three measurements are taken to determine the
probe response $\delta^{(1)}$ for each training iteration: pump only; probe
only; and both pump and probe. In this way, the background terms due to pump
fluorescence and unabsorbed probe could be negated.

Acknowledgements This work is supported by Innovate UK Smart Grant 10043476.
X.G. acknowledges support from the Royal Commission for the Exhibition of 1851
Research Fellowship.

Author contributions X.G. and A.L. conceived the experiment. J.S. carried out
the experiment and performed the data analysis. All the authors jointly
prepared the manuscript. This work was done under the supervision of A.L.

Competing interests The authors declare no competing interests in this work.

## References

  * Roy  _et al._ [2019] K. Roy, A. Jaiswal, and P. Panda, Towards spike-based machine intelligence with neuromorphic computing, Nature 575, 607 (2019).
  * Yao  _et al._ [2020] P. Yao, H. Wu, B. Gao, J. Tang, Q. Zhang, W. Zhang, J. J. Yang, and H. Qian, Fully hardware-implemented memristor convolutional neural network, Nature 577, 641 (2020).
  * Xu  _et al._ [2021] X. Xu, M. Tan, B. Corcoran, J. Wu, A. Boes, T. G. Nguyen, S. T. Chu, B. E. Little, D. G. Hicks, R. Morandotti, _et al._ , 11 tops photonic convolutional accelerator for optical neural networks, Nature 589, 44 (2021).
  * Feldmann  _et al._ [2021] J. Feldmann, N. Youngblood, M. Karpov, H. Gehring, X. Li, M. Stappers, M. Le Gallo, X. Fu, A. Lukashchuk, A. S. Raja, _et al._ , Parallel convolutional processing using an integrated photonic tensor core, Nature 589, 52 (2021).
  * Sludds  _et al._ [2022] A. Sludds, S. Bandyopadhyay, Z. Chen, Z. Zhong, J. Cochrane, L. Bernstein, D. Bunandar, P. B. Dixon, S. A. Hamilton, M. Streshinsky, _et al._ , Delocalized photonic deep learning on the internet’s edge, Science 378, 270 (2022).
  * Lin  _et al._ [2018] X. Lin, Y. Rivenson, N. T. Yardimci, M. Veli, Y. Luo, M. Jarrahi, and A. Ozcan, All-optical machine learning using diffractive deep neural networks, Science 361, 1004 (2018).
  * Shen  _et al._ [2017] Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund, _et al._ , Deep learning with coherent nanophotonic circuits, Nature Photonics 11, 441 (2017).
  * Shastri  _et al._ [2021] B. J. Shastri, A. N. Tait, T. F. de Lima, W. H. Pernice, H. Bhaskaran, C. D. Wright, and P. R. Prucnal, Photonics for artificial intelligence and neuromorphic computing, Nature Photonics 15, 102 (2021).
  * Wu  _et al._ [2022] J. Wu, X. Lin, Y. Guo, J. Liu, L. Fang, S. Jiao, and Q. Dai, Analog Optical Computing for Artificial Intelligence, Engineering 10, 133 (2022).
  * Wetzstein  _et al._ [2020] G. Wetzstein, A. Ozcan, S. Gigan, S. Fan, D. Englund, M. Soljačić, C. Denz, D. A. Miller, and D. Psaltis, Inference in artificial intelligence with deep optics and photonics, Nature 588, 39 (2020).
  * LeCun  _et al._ [2015] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, nature 521, 436 (2015).
  * Spall  _et al._ [2022] J. Spall, X. Guo, and A. I. Lvovsky, Hybrid training of optical neural networks, Optica 9, 803 (2022).
  * Buckley  _et al._ [2023] S. M. Buckley, A. N. Tait, A. N. McCaughan, and B. J. Shastri, Photonic online learning: a perspective, Nanophotonics 12, 833 (2023).
  * Wright  _et al._ [2022] L. G. Wright, T. Onodera, M. M. Stein, T. Wang, D. T. Schachter, Z. Hu, and P. L. McMahon, Deep physical neural networks trained with backpropagation, Nature 601, 549 (2022).
  * Zhou  _et al._ [2021] T. Zhou, X. Lin, J. Wu, Y. Chen, H. Xie, Y. Li, J. Fan, H. Wu, L. Fang, and Q. Dai, Large-scale neuromorphic optoelectronic computing with a reconfigurable diffractive processing unit, Nature Photonics 15, 367 (2021).
  * Tam  _et al._ [1990] S. M. Tam, B. Gupta, H. A. Castro, and M. Holler, Learning on an analog vlsi neural network chip, in  _1990 IEEE International Conference on Systems, Man, and Cybernetics Conference Proceedings_ (IEEE, 1990) pp. 701–703.
  * Mennel  _et al._ [2020] L. Mennel, J. Symonowicz, S. Wachter, D. K. Polyushkin, A. J. Molina-Mendoza, and T. Mueller, Ultrafast machine vision with 2d material neural network image sensors, Nature 579, 62 (2020).
  * Li  _et al._ [2018] C. Li, D. Belkin, Y. Li, P. Yan, M. Hu, N. Ge, H. Jiang, E. Montgomery, P. Lin, Z. Wang, _et al._ , Efficient and self-adaptive in-situ learning in multilayer memristor neural networks, Nature communications 9, 1 (2018).
  * Wang  _et al._ [2019] Z. Wang, C. Li, P. Lin, M. Rao, Y. Nie, W. Song, Q. Qiu, Y. Li, P. Yan, J. P. Strachan, _et al._ , In situ training of feed-forward and recurrent convolutional memristor networks, Nature Machine Intelligence 1, 434 (2019).
  * Feldmann  _et al._ [2019] J. Feldmann, N. Youngblood, C. D. Wright, H. Bhaskaran, and W. H. Pernice, All-optical spiking neurosynaptic networks with self-learning capabilities, Nature 569, 208 (2019).
  * Bandyopadhyay  _et al._ [2022] S. Bandyopadhyay, A. Sludds, S. Krastanov, R. Hamerly, N. Harris, D. Bunandar, M. Streshinsky, M. Hochberg, and D. Englund, Single chip photonic deep neural network with accelerated training, arXiv preprint arXiv:2208.01623 (2022).
  * Filipovich  _et al._ [2022] M. J. Filipovich, Z. Guo, M. Al-Qadasi, B. A. Marquez, H. D. Morison, V. J. Sorger, P. R. Prucnal, S. Shekhar, and B. J. Shastri, Silicon photonic architecture for training deep neural networks with direct feedback alignment, Optica 9, 1323 (2022).
  * Bartunov  _et al._ [2018] S. Bartunov, A. Santoro, G. E. Hinton, B. A. Richards, L. Marris, and T. P. Lillicrap, Assessing the scalability of biologically-motivated deep learning algorithms and architectures, Advances in Neural Information Processing Systems 2018-Decem, 9368 (2018), arXiv:1807.04587 .
  * Hughes  _et al._ [2018] T. W. Hughes, M. Minkov, Y. Shi, and S. Fan, Training of photonic neural networks through in situ backpropagation and gradient measurement, Optica 5, 864 (2018).
  * Pai  _et al._ [2023] S. Pai, Z. Sun, T. W. Hughes, T. Park, B. Bartlett, I. A. Williamson, M. Minkov, M. Milanizadeh, N. Abebe, F. Morichetti, _et al._ , Experimentally realized in situ backpropagation for deep learning in photonic neural networks, Science 380, 398 (2023).
  * Guo  _et al._ [2021] X. Guo, T. D. Barrett, Z. M. Wang, and A. Lvovsky, Backpropagation through nonlinear units for the all-optical training of neural networks, Photonics Research 9, B71 (2021).
  * Wagner and Psaltis [1987] K. Wagner and D. Psaltis, Multilayer optical learning networks, Applied Optics 26, 5061 (1987).
  * Spall  _et al._ [2020] J. Spall, X. Guo, T. D. Barrett, and A. Lvovsky, Fully reconfigurable coherent optical vector–matrix multiplication, Optics Letters 45, 5752 (2020).
  * Ohno  _et al._ [2022] S. Ohno, R. Tang, K. Toprasertpong, S. Takagi, and M. Takenaka, Si microring resonator crossbar array for on-chip inference and training of the optical neural network, ACS Photonics 9, 2614 (2022).
  * Zhou  _et al._ [2020] T. Zhou, L. Fang, T. Yan, J. Wu, Y. Li, J. Fan, H. Wu, X. Lin, and Q. Dai, In situ optical backpropagation training of diffractive optical neural networks, Photonics Research 8, 940 (2020).
  * Boyd [2020] R. W. Boyd,  _Nonlinear optics_ (Academic press, 2020).
  * Zhang  _et al._ [2021] H. Zhang, M. Gu, X. Jiang, J. Thompson, H. Cai, S. Paesani, R. Santagati, A. Laing, Y. Zhang, M. Yung, _et al._ , An optical neural chip for implementing complex-valued neural network, Nature communications 12, 457 (2021).
  * Cruz-Cabrera  _et al._ [2000] A. A. Cruz-Cabrera, M. Yang, G. Cui, E. C. Behrman, J. E. Steck, and S. R. Skinner, Reinforcement and backpropagation training for an optical neural network using self-lensing effects, IEEE transactions on neural networks 11, 1450 (2000).
  * Arrizón  _et al._ [2007] V. Arrizón, U. Ruiz, R. Carrada, and L. A. González, Pixelated phase computer holograms for the accurate encoding of scalar complex fields, JOSA A 24, 3500 (2007).


11institutetext:  Heriot-Watt University

# CNN based IoT Device Identification

Kahraman Kostas

###### Abstract

While the use of the Internet of Things is becoming more and more popular,
many security vulnerabilities are emerging with the large number of devices
being introduced to the market. In this environment, IoT device identification
methods provide a preventive security measure as an important factor in
identifying these devices and detecting the vulnerabilities they suffer from.
In this study, we present a method that identifies devices in the Aalto
dataset using the convolutional neural network (CNN).

While the use of the Internet of Things is becoming more and more popular,
many security vulnerabilities are emerging with the large number of devices
being introduced to the market. In this environment, IoT device identification
methods provide a preventive security measure as an important factor in
identifying these devices and detecting the vulnerabilities they suffer from.
In this study, we present a method that identifies devices in the Aalto
dataset[1, 2] using the convolutional neural network (CNN) .

The CNN is used in network security or classification of network data, while
the network packet information is converted to pseudo-images to present to the
algorithm. In this part, it is mentioned how the internet flow is translated
into pseudo-pictures. Some of the paths followed in these conversion processes
are talk about in the studies below.

Lim et al. [3] converted the network data into pseudo-images for use on CNN.
For this process, firstly, the payload portion of the network packet was
converted to binary numbers to obtain 4-bit groups (two-part for every byte).
These nibbles were converted to decimal numbers, each of which acted as pixels
of a picture. Pseudo-pictures of 36,64,256,1024 pixels were produced according
to their size in payloads (Please see Figure 1).

Figure 1: Representation of the operation of the algorithm simply for the
study of Lim et al. [3]

For this operation, the length of the series of numbers to create the image is
rounded to the nearest image size (36,64,256 or 1024). For this, the interrupt
or 0 padding method is applied. For example, the algorithm adds 26 zeros to a
series of 10 numbers to reach the closest length option (36). On the other
hand, in a series with a length of 500 numbers, it is ignored by cutting after
256. Some of the pseudo pictures obtained from this process (pictures from the
first package for 10 different devices) are shared below (Please see Figure
2).

Figure 2: Pictures from the first package for 10 different devices for the
study of Lim et al. [3]

Lotfollahi et al. [4] use pcap files for CNN and SAE. These operations can be
listed as follows. The Ethernet header is deleted. Standard packets such as
DNS or tree-way handshake packets are deleted because they do not carry
payloads and are not effective in determining the class of data. Adding 0 to
the end of UDP packets and IP payloads less than 1480 Bytes ensures that all
data arrays are the same size. Then normalization is performed by dividing
each byte to 255. The result is a pseudo-image of 1480 pixels (37x40). The
representation of this operation can be seen from Figure 3.

Figure 3: Representation of the operation of the algorithm simply for the
study of Lotfollahi et al. [4]

Some of the pseudo pictures obtained from this process (pictures from the
first package for 10 different devices) are shared below (See Figure 4).

Figure 4: Pictures from the first package for 10 different devices for the
study of the study of Lotfollahi et al. [4]

In their study, Wang et al. [5] convert 4 different data groups (Session +
All, Session + L7, Flow + All, Flow + L7) to pseudo-pictures of 784 (28 * 28)
pixels to be used in CNN. In this process, Session and flow refer to the
5-tuples such as source IP, source port, destination IP, destination port and
transport-level protocol. However, while the session is bi-directional, flow
is one-way. L7 stands for application layer and ALL stands for all protocol
layers. The representation of this operation can be seen from Figure 5.

Figure 5: Representation of the operation of the algorithm simply for the
study of Wang et al. [5]

We have performed only one (Session + L7) of these options and shared the
methods and results below. The first 13 pixels of the generated image
symbolize the session (Source IP:4 Byte, Destination IP:4 Byte, Source Port:2
Byte, Destination Port:2 Byte, UDP/TCP:1 Byte) and the remaining 771 pixels
symbolize the payload. Some of the pseudo pictures obtained from this process
(pictures from the first package for 10 different devices) are shared below
(See Figure 6).

Figure 6: Pictures from the first package for 10 different devices for the
study of the study of Wang et al. [5]

All three of the above-mentioned studies have been done on traditional
networks. In the application below, using these methods, devices will be
classified through payload from the Pcap file obtained from IoT appliance.

In the approach of Lim et al. [3], It is seen that the data has been put into
more than one form, and this approach will require training of different
models. However, we preferred to give this data a standard shape and evaluate
it with a single model. For this, we chose to get the first 784 bytes like
Wang et al. This is because the size of the packets produced by the majority
of the network packets produced by the devices ranges from 0 to 800 (Please
see Figure 6).

Figure 7: The distribution of the size of the network packets

On the other hand, unlike Wang et al. [5], we did not use the features
extracted from the layer headers because we see this process as some kind of
feature extraction. Whereas with this application, we want to focus on to
classify the devices using only the payload.

In this step, we will compare the performance of pseudo-images obtained from
packet payloads and pseudo-images obtained from the feature sets from
IoTDevIDv1 [6] using CNN. The process of creating pseudo-images from payload
can be summarized briefly as follows111The script available at:
github.com/kahramankostas/CNN-based-IoT-Device-Identification: During this
process, all headers of the internet packages are deleted and the data
obtained is converted into a two-dimensional matrix (28X28). In this process,
packages that do not contain any loads such as 3-handshake packages are
ignored. If the package contains more than 784 bytes of data, the excess
portion is ignored. However, if the package size contains less than 784 data,
it is brought to this size by adding ”0” to the end of the payload. In order
to select the optimum hyperparameters of CNNs in these two processes, the
ImageClassifier feature of the library named Autokeras[7], which scans
parameters in ANNs, is used. Selected parameters and architecture can be seen
in Appendix K. The data obtained as a result of this process is given in Table
1 below. It can be seen the distribution of accuracy scores in Figure 8.

| CNN with Payload | CNN with Fingerprint
---|---|---
Accuracy | 0.631 | 0.625
Time (seconds) | 209.926 | 20.809
Utest (CNN with Payload, CNN with Fingerprint) | 5.84e-09 |

Table 1: Comparison of CNN Results with 2 approaches.
The significance level is 0.05 for the Mann Whitney U test.  Figure 8:
Distribution of Accuracy Values

When Table 1 is looked at, it is seen that the results obtained with the
fingerprint method are quite close to the results obtained by examining the
payloads (accuracy score is approximately 0.006 lower). However, if the speeds
of these two methods are compared, it is seen that the fingerprint method is
about 10 times faster than the other method.

## References

  * [1] S. Marchal, “IoT devices captures, aalto university,” 2017.  Accessed: 2021-08-25.
  * [2] M. Miettinen, S. Marchal, I. Hafeez, N. Asokan, A.-R. Sadeghi, and S. Tarkoma, “Iot sentinel: Automated device-type identification for security enforcement in iot,” in 2017 IEEE 37th International Conference on Distributed Computing Systems (ICDCS), pp. 2177–2184, IEEE, 2017.
  * [3] H.-K. Lim, J.-B. Kim, J.-S. Heo, K. Kim, Y.-G. Hong, and Y.-H. Han, “Packet-based network traffic classification using deep learning,” in 2019 International Conference on Artificial Intelligence in Information and Communication (ICAIIC), pp. 046–051, IEEE, 2019.
  * [4] M. Lotfollahi, M. J. Siavoshani, R. S. H. Zade, and M. Saberian, “Deep packet: A novel approach for encrypted traffic classification using deep learning,” Soft Computing, vol. 24, no. 3, pp. 1999–2012, 2020.
  * [5] W. Wang, M. Zhu, J. Wang, X. Zeng, and Z. Yang, “End-to-end encrypted traffic classification with one-dimensional convolution neural networks,” in 2017 IEEE International Conference on Intelligence and Security Informatics (ISI), pp. 43–48, IEEE, 2017.
  * [6] K. Kostas, M. Just, and M. A. Lones, “IoTDevID: A behaviour-based fingerprinting method for device identification in the IoT,” 2021.
  * [7] H. Jin, Q. Song, and X. Hu, “Auto-keras: An efficient neural architecture search system,” in Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 1946–1956, ACM, 2019.

Appendix

Appendix K-1
CNN Model Architecture for the Fingerprint dataset.

Figure 9: CNN Model Architecture for the Fingerprint dataset.

Appendix K-2
CNN Model Architecture for Payload.

Figure 10: CNN Model Architecture for Payload.