sectioning

TUM-HEP 1054/16

On the challenge of estimating diphoton backgrounds at large invariant mass

J.F. Kamenika,b, G. Perezc, M. Schlafferc, A. Weilerd

a Jožef Stefan Institute, Jamova 39, 1000 Ljubljana, Slovenia

b Faculty of Mathematics and Physics, University of Ljubljana,

Jadranska 19, 1000 Ljubljana, Slovenia

c Department of Particle Physics and Astrophysics, Weizmann Institute of
Science,

7610001 Rehovot, Israel

d Physik Department T75, Technische Universität München, James-Franck-Strasse
1,

85748 Garching, Germany

Diphoton searches at high invariant mass are an integral part of the
experimental high energy frontier. Using the analyses of the 750 GeV diphoton
resonance as a case study, we examine the methodology currently employed by
the experimental analyses in estimating the dominant standard model
backgrounds. We assess the dependence of the significance associated with the
excess on the background modeling. In particular, we show that close to the
high energy tails of the distributions, estimates of the jet faking
backgrounds relying on functional extrapolations or Monte Carlo estimates of
the challenging photon-jet contributions introduce a large uncertainty.
Analyses with loose photon cuts, low photon $p_{T}$ cuts and those susceptible
to high photon rapidity regions are especially affected. Given that diphoton-
based searches beyond 1 TeV are highly motivated as discovery modes, these
considerations are relevant and applicable in the future.

We first consider a physics-driven deformation of the shape of the diphoton
faking photon-jet spectrum by next-to-leading order effects combined with a
rapidity and transverse momentum dependent fake rate. We show that the
resulting local significance of the excess is reduced due to such a
deformation. Then using a simple, but more general, ansatz to modify the
contribution of the jet faking backgrounds at high invariant masses but
keeping the inclusive and differential sample purities within experimental
uncertainty estimates, we demonstrate that the originally reported local 750
GeV excess significances could have been overestimated by more than one
standard deviation. We furthermore cross-check our analysis by comparing fit
results based on the smaller 2015 and the larger 2016 LHC datasets. Finally we
employ our methodology on the complete 13 TeV LHC dataset to asses the
systematics involved in the current diphoton searches beyond the TeV region.

##  1 Introduction

Searches for new physics at the energy frontier often look for new phenomena
at the edge of distributions. In this kinematical region the knowledge of the
Standard Model (SM) background is typically limited and the challenge is to
look for a new resonance where only partial knowledge on the SM background is
available. In this paper we focus in particular on new physics probes based on
the high diphoton invariant mass spectrum. We examine, using the analyses of
the 750 GeV diphoton resonance as a case study, the strategy currently used by
the experimental collaborations in estimating the dominant SM backgrounds. We
employ our methodology on the 13 TeV LHC dataset to asses the systematics
involved in the current diphoton searches beyond the TeV region

In their 2015 data sets, both ATLAS and CMS observed an excess in the diphoton
spectrum near $m_{\gamma\gamma}=750\,\,\text{GeV}\,.$ The relevant details of
the ATLAS and CMS analyses are described in [1, 2]. At face value the local
significances for a broad resonance were given by

$\begin{aligned} p_{\textrm{ATLAS}}&=4\times 10^{-5}\,,\\\
\sigma_{\textrm{ATLAS}}&=3.9\,,\end{aligned}\qquad\begin{aligned}
p_{\textrm{CMS}}&=5\times 10^{-3}\,,\\\
\sigma_{\textrm{CMS}}&=2.6\,,\end{aligned}\qquad\begin{aligned}
p_{\textrm{comb}}&=1\times 10^{-6}\,,\\\
\sigma_{\textrm{comb}}&=4.7\,,\end{aligned}$ (1.1)

where $p_{\textrm{ATLAS, CMS, comb}}$ $(\sigma_{\textrm{ATLAS, CMS, comb}})$
correspond to the local $p$ value (confidence level) of ATLAS, CMS, and their
naive combination.111We do not discuss here the global significance as it
strongly depends on the lower value of $m_{\gamma\gamma}$ defined for the
search region. ATLAS (CMS) chose it to be about 200 GeV (400 GeV).
Furthermore, as discussed below, the region below 500 GeV is dominating the
fit to the functional form which is used to estimate the background. Thus, it
is not clear whether one should consider this region as a control region or as
the region of interest for the search itself.

The local results quoted in Eq. (1.1) are quite significant and captured the
attention of the high energy community. Interpreting them naively, one would
be lead to one of the following conclusions

  1. (i)

this excess is a result of a rare statistical fluctuation;

  2. (ii)

this excess implies a discovery of non-Standard Model dynamics.

As both conclusions are quite extraordinary (certainly the second one), they
motivate an investigation into their robustness. In particular, we raise a
third option, to be considered in conjunction with (i), namely, we ask how
unlikely is the possibility that

  1. (iii)

the significance of the excess is overestimated due to underestimating fake-
based backgrounds.

With the inclusion of more data in the analyses the excess eventually vanished
[3, 4], ruling out the new physics hypothesis ii. However, the possibility of
claim iii remains unclear, affecting all analyses which rely on a precise
knowledge of the photon faking background and use the same techniques to
estimate it.

While our conclusion is independent of the 750 GeV resonance we use it as an
example case to scrutinize the hypothesis of the underestimated background and
its implications. First, the main rationale behind our hypothesis is presented
in Sec. 2, followed by a detailed description of our approach to background
estimation (Sec. 3) and the statistical treatment of the data (Sec. 4). The
comparison with the full 2016 data set is presented in Sec. 5. Our main
conclusions are summarized in Sec. 6. For other relevant works, see Refs. [5,
6].

##  2 The rationale

Superficially, the experimental situation related to the diphoton excess was
fairly straightforward. The experiments had reported a relatively narrow
“bump”, $\Gamma/m\lesssim 6\%\ll 1$. Such a bump implies a rise in the
differential distribution while, due to the rapidly falling parton luminosity
functions, it is expected that any reasonable background-related distribution
should be a monotonically decreasing function of the invariant mass.
Consequently, the presence of a non-Standard Model feature seemed to have been
indicated by the measurements. While this was qualitatively correct the
challenge is to quantify the significance of the excess. To endow the bump
with a significance, one needs to control and quantify the background.

The following approaches can be used to constrain the form of the background:

I. Data-driven approach. Assuming $\Gamma/m\ll 1$ and a featureless monotonic
background, a robust way to constrain it is through _interpolation_ via a two-
sided side band analysis. However, this requires to have enough measured
events at invariant masses both below the resonance and above it. In the case
of the 750 GeV excess, there were less than 40 events in all of the analyses
measured with invariant masses above 850 GeV. Such a small number of events
does not allow one to use this method reliably.

II. ”First-principle”/Monte-Carlo approach. There is a rather narrow class of
observables for which the theory has reached an advanced enough level such
that we can fully trust our ability to correctly predict the shape of the
background distributions. We believe that the invariant mass distribution of
experimentally measured diphoton events does not (yet) belong to this selected
class of observables. Namely, the continuous diphoton distribution consist of
an admixture of two dominant components: (i) The first is made of two real
isolated hard photons. This diphoton distribution is currently known to next-
to-next-to-leading order (NNLO) accuracy [7, 8] in perturbative QCD and
imposing cuts similar to the ATLAS spin-0 analysis suggests an overall
uncertainty of about 5% for the invariant mass distribution [8]. (ii) An
additional important background component is due to fakes coming mostly from
processes involving a hard photon and a jet that passes the various photon
quality and isolation cuts [9]. In addition, depending on these cuts, also the
dijet background could play an important role. The prompt photon-jet cross
section is currently known at next-to-leading order (NLO) in QCD, and several
codes are available to produce the relevant distributions, including JetPhox
[10] and PeTeR [11]. In addition, QCD threshold resummation at next-to-next-
to-next-to-leading logarithmic (N3LL) order [12, 13] as well as electroweak
Sudakov effects are being included [14], resulting in theory uncertainties of
about 10-20% [14, 9, 15]. However, a comparison with the 8 TeV ATLAS
measurement [9] shows that at low photon $p_{T}\sim 50\,$GeV the data exhibits
some level of deviations from the theoretical predictions (a larger
uncertainty is found for the invariant mass distribution, see [16]). In
addition, it is important to note, that the fake rate strongly depends on the
quark/gluon “flavor” of the tagged jet (for some discussion on jet flavor
definitions see [17, 18, 19, 20]): intuitively one can understand the
difference through the quark and gluon fragmentation functions to pions. At
large $x$, as required to be able to pass photon isolation criteria, gluon
fragmentation to few pions is much more suppressed (see e.g. Chapter 20 in
Ref. [21]). Accordingly, a dedicated ATLAS study [22] found that there is a
probability of about $1:2\times 10^{3}$ for a quark jet to fake a photon, and
only $1:2\times 10^{4}$ for a gluon jet to fake a photon, for jets with
$E_{T}>40\,$GeV. Applying this to the photon-jet background, we also note that
subleading jets might become an important source of fakes if the leading jet
is predominantly gluon-initiated.

In order to theoretically predict the purity of the diphoton mass
distribution, an appropriate admixture of the diphoton and the photon-jet(s)
components needs to be constructed [23]. Furthermore, for the latter
component, one is required to convolve the photon-jet distribution with the
relevant fragmentation functions or at least tag the flavor of the jet(s). It
is also important to note that the purity is a highly phase space dependent
quantity. Not only does it depend on the ratio of the differential jet-photon
and photon-photon production but also on the jet-to-photon fake rate. The fake
rate may exhibit strong dependence on the differential quantities such as
$p_{T}$ and (pseudo)rapidity $\eta$. For instance, as discussed below, in the
CMS analyses purity is estimated to be better than 90% in the (central-
central) EBEB event category but only better than 80% in the (forward-central)
EBEE one. Both experiments consider the purity in an inclusive way. However,
in the relevant kinematical region the data is not sufficient to constrain
possibly large deviations from the inclusive purity estimation (see Fig. 4).

III. Functional-fit approach. Given the present practical limitations of the
methods I and II, one is lead to a more phenomenological approach in which the
background estimate is obtained by fitting an universal function to control
regions in the data and then _extrapolating_ into the signal regions using the
fitted functional form. This allows one to predict the background at
relatively high invariant masses in a straightforward manner. Consequently,
both experiments are essentially following this approach in most of their
analyses222An exception is the ATLAS spin-2 analysis which employs a Monte
Carlo approach (II) with a data-driven estimate of the photon-jet and jet-jet
background, see Section 4., although the functional forms used by ATLAS in the
spin-0 analysis and by CMS are slightly different. Thus, the significance of
the excess is mostly determined by comparing measured events to a background
estimate predicted by a fitting function.

While method III is very transparent and makes the search for bumps easy to
analyze, it is also rather susceptible to systematic effects, in particular a
lack of understanding of the physics modifying the tails of the distributions,
as we argue below. The fitting functions used by ATLAS and CMS are well suited
for describing rapidly falling distributions and are fitted to the available
data. With the amount of data in the 2015 data sets, the differentially
measured number of events is abundant in the low invariant mass region and is
spare in the high mass region. The extraction of the functions’ parameters is
thus dominantly controlled by the low $m_{\gamma\gamma}$ region and hardly
affected by modifications of the invariant mass distribution at diphoton
masses of above roughly 500 GeV. However, the significance of the excess with
respect to the fitting function is very much affected by such deformations. As
it is hard to directly test or predict the correct form of the diphoton mass
distribution this raises the following questions:
${\mathcal{I}}$. Is the experimental signal over background estimation robust
against the presence of deviations from the fitting function predictions at
large invariant masses?
${\cal II}$. If this is not the case, can one produce smoking-gun predictions
to show that indeed the significance of the excess is being overestimated?

Let us first focus on point ${\cal I}$. To examine the sensitivity of the
significance of the excess to the variation of the tails of the distributions.
We consider a family of background shapes that are formed by an admixture of
the diphoton and photon-jet distributions. We keep the overall inclusive
purity of the samples at 90% and 80%, respectively, in accordance with the
measured data at low invariant masses. More specifically, we use two classes
of deformations. The first is derived from a modification of the photon-jet
spectrum due to NLO and showering effects combined with an increased fake rate
for larger transverse momenta and pseudo-rapidities of the jets.

We then consider a simpler ansatz where we allow the distribution of the
$pp\to\gamma j$ component to be reweighted at invariant masses above 500 GeV
such that the purity of events with large invariant masses is reduced leading
to a controlled deviation from the functional fit. In the following section we
provide a detailed description of our approach. We also provide some tests of
our procedure to check that our method complies with public data (below and
above the resonance region) and is passing the relevant statistical tests. We
then report how the significance is affected by the amount of rescaling of the
distributions of fakes. Finally we can use our ansatz to address item ${\cal
II}$ and provide smoking guns to test our hypothesis on overestimating the
excess significance. With the full statistics of the 2016 data sets at hand it
would be fairly easy to eliminate our hypothesis.

##  3 Reducible and irreducible backgrounds

The main background to the diphoton signal is the irreducible
$pp\to\gamma\gamma$ background. We consider in the following: ATLAS spin-0 and
spin-2, and CMS $13\,\,\text{TeV}$ EBEB and EBEE categories with magnets on.
We generate the diphoton invariant mass spectrum at NNLO with MCFM version 8.0
[24, 25, 26, 27, 8] applying the cuts as described in the respective analyses,
see Table 1. The main contribution to the reducible background is the
$pp\to\gamma j$ production where the hard jet is a quark jet that is wrongly
reconstructed as a photon. We generate this background at leading order (LO)
with MadGraph5 version 5.2 [28]. We note that at LO the $pp\to\gamma j$ sample
is dominated by quark jets, which as already mentioned, lead to a much larger
fake rate than gluon jets.

Analysis | ATLAS spin-0 | ATLAS spin-2 | CMS EBEB | CMS EBEE
---|---|---|---|---
$m_{\gamma\gamma}$ | $>150\,\,\text{GeV}$ | $>200\,\,\text{GeV}$ | $>230\,\,\text{GeV}$ | $>330\,\,\text{GeV}$
$p_{T,1}$ | $>0.4\,m_{\gamma\gamma}$ | $>55\,\,\text{GeV}$ | $>75\,\,\text{GeV}$ | $>75\,\,\text{GeV}$
$p_{T,2}$ | $>0.3\,m_{\gamma\gamma}$ | $>55\,\,\text{GeV}$ | $>75\,\,\text{GeV}$ | $>75\,\,\text{GeV}$
$\left|\eta_{1}\right|$ | $<2.37$ | $<2.37$ | $<1.44$ | $<1.44$
$\left|\eta_{2}\right|$ | $<2.37$ | $<2.37$ | $<1.44$ | $1.57<\eta_{2}<2.5$
$\left|\eta\right|$ excluded | $1.37<\eta_{1,2}<1.52$ | $1.37<\eta_{1,2}<1.52$ | n.a. | n.a.
$\sigma_{\gamma\gamma}$ [pb] (NNLO) | 2.7 | 1.9 | 0.52 | 0.23
$\sigma_{\gamma j}$ [pb] (NLO) | 1400 | 1000 | 250 | 130
$\left.{\sigma_{\gamma j}}{}/{\sigma_{\gamma\gamma}}{}\right|_{m>500\,\,\text{GeV}}$ | 510 | 670 | 470 | 640
Table 1: Cuts of the analyses where the subscript refers to the hardest and
second hardest photon candidate, the cross section of the $pp\to\gamma\gamma$
sample passing these cuts (calculated at NNLO with MCFM) and of the
$pp\to\gamma j$ sample at hadron level (before applying any photon mistag
rate), calculated at NLO with MadGraph5_aMC@NLO, showered with Pythia [29] and
the jets clustered with anti-$k_{T}$, $R=0.4$ algorithm using FastJet [30]. In
the last line the ratios of the two distributions in the invariant mass region
above 500 GeV are given.

The reconstructed diphoton distribution is a mixture of $pp\to\gamma\gamma$
and $pp\to\gamma j$ invariant mass distributions. Let us define a short-hand
notation for the normalized invariant mass distribution

$w_{\gamma X}\equiv\frac{1}{\sigma_{\gamma
X}}\times\frac{\textnormal{d}\sigma_{\gamma X}}{{\textnormal{d}m_{\gamma
X}}}\,,$ (3.2)

with $X=\gamma,j$. The mixed distribution $w_{\textnormal{mix}}$ is a function
of the normalisation $N_{\textnormal{mix}}$ and a parameter $\mathcal{R}$,
that controls the shape modification of $w_{\gamma j}$ and will be defined Eq.
(3.5). We write $w_{\textnormal{mix}}$ as

$w_{\textnormal{mix}}(N_{\textrm{mix}})=N_{\textrm{mix}}\left[\mathcal{P}w_{\gamma\gamma}+(1-\mathcal{P})w_{\gamma
j}\right]\,,$ (3.3)

where $\mathcal{P}$ is the inclusive purity of the sample. We set
$\mathcal{P}=90\%$ (80%) for the ATLAS and CMS EBEB (EBEE) analyses, which is
within the reported error bands. We will assume, that $w_{\gamma\gamma}$ is
obtained by normalizing the MCFM diphoton invariant mass distribution. As for
$w_{\gamma j}$, following the rationale described in Section 2, we modify in
two different ways as we now describe in detail.

###  3.1 QCD and jet-fake dependence of the diphoton shape

First, we calculate a photon-jet mass dependent K-factor using
MadGraph5_aMC@NLO, showered and hadronized with Pythia [29] and jets clustered
with an anti-$k_{T}$, $R=0.4$ algorithm [31] using FastJet [30]. We note that
in the NLO distribution we only consider the hardest jet of the event and we
do not record its flavor. This step may be potentially improved by the use of
an IR-safe jet flavor definition, see [17, 18, 19, 20]. Next, in order to
model the dependence of the fake rate on the pseudo rapidity and the
transverse momentum we use the following simplified ansatz for the jet
rejection $r(p_{T},\eta)$:

$r(p_{T},\eta)={\textrm{max}}\left\\{\frac{r_{0}}{1+p_{T}/p_{T}^{0}+\eta/\eta^{0}},r_{\textrm{min}}\right\\}\,,$
(3.4)

where the functional form is motivated by the kinematical dependence of the
jet-rejection rates as estimated by ATLAS [22] and the parameter values
$p_{T}^{0}=30$ GeV, $\eta^{0}=4$ are chosen to reproduce the rejection rate
ratios between the lowest and highest lying $\eta$ and $p_{T}$ bins within
uncertainties. Finally, $r_{0}/r_{\textrm{min}}$ is fairly uncertain as
estimates of rejection rates at very high $p_{T}$ and $\eta$ are not publicly
available, but reproducing experimental purity estimates in the forward region
[2] leads to values in the wide range $r_{0}/r_{\textrm{min}}\in[3,12]$ . The
resulting reweighting factors compared to the LO partonic $m_{\gamma j}$
distribution obtained from MadGraph5, $w_{\gamma j}^{\textnormal{MG}}$, at
both steps applied successively ($w_{\gamma j}^{\textrm{NLO}}/w_{\gamma
j}^{\textrm{MG}}$ and $w_{\gamma
j}^{\textrm{NLO}\,\times\,\textrm{fakes}}/w_{\gamma j}^{\textrm{MG}}$) are
shown in Fig. 1. We observe that with our choice of fake rate parameters, the
largest reweighting factors close to 3 are obtained above $m_{\gamma j}>800$
GeV for the ATLAS spin-2 cuts. However, all experimental categories are
affected by a reweighting factor which is a combination of NLO, hadronization
and faking effects, and which increases with the photon-jet invariant mass
until it saturates at some point. This suggests a simple functional form for
the effective photon-jet spectrum deformation which we discuss next.

Figure 1: Reweighting factor for $w_{\gamma j}$ with respect to the LO Monte
Carlo distribution as a function of the invariant mass for ATLAS (left) and
CMS (right). The dashed lines show the reweighting factor to modify the LO
parton distribution $w_{\gamma j}^{\textnormal{MG}}$ to the NLO shape,
including the effects of hadronization. The solid lines include in addition
the reweighting due to a phase space dependent fake rate.

###  3.2 Effective shape deformation

In our effective ansatz for the deformation of the photon-jet spectrum, we
focus on the invariant mass region above $m_{\gamma\gamma}>500$ GeV below
which the experiments have sufficient statistics to control the mass
distributions and calibrate their analyses and choose a simple, linear form

$\displaystyle w_{\gamma j}(\mathcal{R})$ $\displaystyle=\frac{w_{\gamma
j}^{\textrm{MG}}}{N(\mathcal{R})}\left[1+\mathcal{R}\times\begin{cases}0&m_{\gamma
j}<500\,\,\text{GeV}\\\ \frac{\left.w_{\gamma
j}^{\textnormal{MG}}\right|_{m_{\gamma j}=500\,\,\text{GeV}}}{w_{\gamma
j}^{\textnormal{MG}}}-1&500\,\,\text{GeV}\leq m_{\gamma j}\leq
760\,\,\text{GeV}\\\ \frac{\left.w_{\gamma
j}^{\textnormal{MG}}\right|_{m_{\gamma j}=500\,\,\text{GeV}}}{\left.w_{\gamma
j}^{\textnormal{MG}}\right|_{m_{\gamma
j}=760\,\,\text{GeV}}}-1&760\,\,\text{GeV}<m_{\gamma j}\end{cases}\right]\,,$
(3.5)

to roughly account for an overall kinematic dependence of the fake rate,
hadronization and higher order effects $N(\mathcal{R})$ is an
$\mathcal{R}$-dependent normalisation factor. We define $w_{\gamma
j}(\mathcal{R})$ such that $w_{\gamma j}(\mathcal{R}=0)$ corresponds to the
partonic LO distribution. Choosing $\mathcal{R}=1$, the shape of $w_{\gamma
j}$ is unmodified for $m_{\gamma j}<500\,\,\text{GeV}$, then flat up to
$m_{\gamma j}=760\,\,\text{GeV}$, just above the observed peak of the apparent
excess, and finally is rescaled by the ratio of the differential cross
sections at $m_{\gamma j}=500\,\,\text{GeV}$ and $m_{\gamma
j}=760\,\,\text{GeV}$ for $m_{\gamma j}>760\,\,\text{GeV}$. The
$\mathcal{R}$-dependence of $N(\mathcal{R})$ is chosen such that the integral
over $w_{\gamma j}$ is always 1, independent of the value of $\mathcal{R}$.
For the ATLAS spin-2 and the two CMS analyses, the intervals in the above
equation are shifted by $10\,\,\text{GeV}$ to larger values due to the
different binning in these searches.

In the left panel of Fig. 2, we show the normalized differential $pp\to\gamma
j$ cross section for the ATLAS spin-0 analysis as a function of the invariant
mass for several choices of $\mathcal{R}$. In the right panel, the
$\mathcal{R}$-dependent reweighting factor of Eq. (3.5) is shown. Since the
spin-0 analysis applies the strongest cuts on the transverse momenta of the
photon candidates ($0.4\,m_{\gamma\gamma}$ and $0.3\,m_{\gamma\gamma}$,
respectively) its distribution is the steepest. Thus the reweighting factor of
this analysis is the largest being almost 7 above 770 GeV. The maximal
reweighting factors for the other analyses are just above 6. We verified that
increasing the flat region by 20 GeV has only a small impact on the reported
results.

(a)

(b)

Figure 2: Left: Normalized invariant mass distributions after the ATLAS spin-0
cuts for the $pp\to\gamma j$ background with several choices for the
interpolating parameter $\mathcal{R}$. Right: Corresponding reweighting
factor.

In Fig. 3(a) and Fig. 3(c), we show the resulting invariant mass distributions
$w_{\textnormal{mix}}$ for the ATLAS spin-0 and CMS EBEB analysis,
respectively, on top of the normalized distribution as measured in the 2015
data set.

(a) ATLAS spin-0

(b) ATLAS spin-0

(c) CMS EBEB

(d) CMS EBEB and EBEE combined

Figure 3: Top left: Combined distributions $w_{\textnormal{mix}}$ for the
ATLAS spin-0 cuts, all with an overall purity of 90%. The distribution
obtained from the 2015 data sets is shown in blue. Bottom left: Corresponding
distribution for the CMS EBEB analysis. Top right: $\Delta\chi^{2}$ of fit to
the ATLAS spin-0 distribution as a function of $\mathcal{R}$. The 1- and
2-$\sigma$ regions are indicated by the thin lines. Bottom right:
corresponding plot for the combined CMS EBEB and EBEE analyses.

In addition to the average purity of the full sample, ATLAS and CMS try to
estimate the purity as a function of the diphoton invariant mass. This local
purity is given by

$\mathcal{P}^{i}=\frac{\mathcal{P}w^{i}_{\gamma\gamma}}{\mathcal{P}w^{i}_{\gamma\gamma}+(1-\mathcal{P})w^{i}_{\gamma
j}}$ (3.6)

for the $i$-th bin. It can deviate significantly from the average purity
$\mathcal{P}$ of the full sample. In Fig. 4, we show the binned purities for
the mixed samples with several choices of $\mathcal{R}$ compared to the purity
determined by ATLAS with the $2\times 2$ sideband [32] and the matrix method
[33] and by CMS with a method described in [34], respectively. We choose the
same binning of the purity as is used in the respective analysis.

(a) ATLAS spin-0

(b) ATLAS spin-2

(c) CMS EBEB

(d) CMS EBEE

Figure 4: Purity of the combined distribution as a function of
$m_{\gamma\gamma}$ for several choices of $\mathcal{R}$. The dashed lines show
the central value for the purity as determined by the experiments and shaded
the areas show the corresponding error.

While the local purity is within the error band in most of the considered mass
range (even for $\mathcal{R}=1$), it does decrease for large invariant masses
and our ansatz predicts a deviation from the experimental value. Given the low
statistics in this range, we consider this as a way to falsify our proposal in
the future rather than a contradiction with the currently available data.

##  4 Statistical treatment

The experimental analyses estimate the background shape by fitting a function
$f(x)$ with $x=m_{\gamma\gamma}/\sqrt{s}$ to the measured data. In the ATLAS
spin-0 analysis, the following ansatz is used:

$\displaystyle
f(x)=N\left(1-x^{1/3}\right)^{b}x^{\sum_{j=0}^{k}a_{j}\left(\log
x\right)^{j}}$ (4.7)

where $k=0$ was chosen. For the CMS analyses as well as the ATLAS spin-2, we
use

$\displaystyle f(x)=Nx^{a+b\log x}\,.$ (4.8)

Note, that the ATLAS spin-2 analysis uses a mixture of Monte Carlo (for the
$pp\to\gamma\gamma$ background) and data driven distributions (for the
$pp\to\gamma j$ and $pp\to jj$ background), leading to similar results as the
fit function approach. In the data driven method, the shape of the different
backgrounds is obtained by extracting the corresponding events from control
samples and fitting their distribution with a function. The relative
contribution to the observed $pp\to\gamma\gamma$ sample is extracted from the
data between $200\,\,\text{GeV}<m_{\gamma\gamma}<500\,\,\text{GeV}$. For more
details on this method see [1]. Given the small statistics in the large
invariant mass bins this approach roughly corresponds to our LO MG
distribution.

In order to see how the significance of the 750 GeV excess changes with our
ansatz, we fit the distribution $w_{\textrm{mix}}$, defined in Eq. (3.3), once
with $w_{\gamma j}^{\textrm{MG}}$ corresponding to $\mathcal{R}=0$, and then
with $w_{\gamma j}$ as estimated at NLO with showering and hadronization,
including fakes and finally with $\mathcal{R}$ as a free fit parameter (as
well as the appropriate fit function $f(x)$) to the measured data. As an
additional template, one could extend the fit function $f(x)$ by a
modification similar to the one described in Eq. (3.5), which we will however
not do for the sake of simplicity. The fits are performed with two methods
which yield similar results.

Firstly, we maximize the likelihood

$\displaystyle
L=\prod_{i=1}^{N_{\textnormal{bins}}}P_{N_{\textnormal{e}}^{i}}(N_{\textnormal{m}}^{i})$
(4.9)

where the product goes over all bins and
$P_{N_{\textnormal{e}}}(N_{\textnormal{m}})$ is the Poisson probability to
measure $N_{\textnormal{m}}$ events when $N_{\textnormal{e}}$ events are
expected.

Secondly, we minimize

$\displaystyle\chi^{2}=\sum_{i=1}^{N_{\textnormal{bins}}}(N_{\textnormal{m}}^{i}-N_{\textnormal{e}}^{i})^{2}/N_{\textnormal{e}}^{i}$
(4.10)

where we rebin the data such that each bin contains at least 10 events in
order for the $\chi^{2}$ distribution to provide a reasonable description of
the statistical uncertainties, see e.g. [35]. In both fit methods, the
overflow of the experimental histograms is treated as one single bin. The
best-fit parameters determine the number of expected events
$N_{\textnormal{e}}$ in the signal region (SR).

Since we are mostly interested in the local significance of the 750 GeV
excess, the SR is chosen by eye from the measured distribution with the aim to
capture the excess. We obtain a $p$-value by comparing the number of measured
events in the SR $N_{\textnormal{m}}$ with $N_{\textnormal{e}}$ (more
precisely: calculating the Poisson probability to measure at least
$N_{\textnormal{m}}$ events):

$p=\sum\limits_{n=N_{\textnormal{m}}}^{\infty}\frac{N_{\textnormal{e}}^{n}}{n!}e^{-N_{\textnormal{e}}}\,.$
(4.11)

Clearly this simple approach which does not use any signal modeling and relies
on a discrete width and position of the SR is far from perfect. Consequently,
it is not surprising that the obtained significances of the excess are smaller
than the ones reported by the experiments, even when we use the same fit
functions. Instead of focusing on absolute values one should therefore rather
consider the reduction of the significance that results from modifying the
background. The results of the fits are shown in Table 2.

| Analysis | ATLAS spin-0 | ATLAS spin-2 | CMS EBEB | CMS EBEE
---|---|---|---|---|---
meas. | SR | 730-770 GeV | 720-760 GeV | 720-780 GeV | 710-770 GeV | 710-770 GeV
$\int\mathcal{L}\,\textnormal{d}t$ | $3.2\,\text{fb}^{-1}$ | $15.4\,\text{fb}^{-1}$ | $3.2\,\text{fb}^{-1}$ | $2.7\,\text{fb}^{-1}$ | $12.9\,\text{fb}^{-1}$ | $2.7\,\text{fb}^{-1}$ | $12.9\,\text{fb}^{-1}$
$N_{m}$ | 15 | 33 | 40 | 12 | 24 | 21 | 53
fitfunction | $-2\log L$ | 270 | - | 330 | 200 | - | 200 | -
$\sigma$ | 3.4 | - | 2.9 | 1.9 | - | 1.7 | -
$\chi^{2}/n$ | 1.6 | 0.75 | 1.2 | 0.80 | 1.0 | 1.2 | 0.98
$\sigma$ | $3.1_{-0.2}^{+0.2}$ | $1.2^{+0.2}_{-0.2}$ | $2.9_{-0.3}^{+0.3}$ | $1.8_{-0.2}^{+0.3}$ | $-1.5^{+0.2}_{-0.3}$ | $1.6_{-0.3}^{+0.3}$ | $-1.3^{+0.2}_{-0.2}$
distribution |  $\mathcal{R}$: $-2\log L$ | 270 | 360 | 330 | 210 | 310 | 210 | 300
$\mathcal{R}$ | 0.86 | 0.11 | 0.97 | -0.048 | -0.7 | 0.24 | -0.12
$\sigma$ | 2.2 | 0.0 | 2.0 | 1.5 | -1.2 | 1.4 | -0.23
MG: $\chi^{2}/n$ | 1.5 | 0.76 | 1.3 | 0.78 | 1.6 | 1.1 | 1.1
$\sigma$ | $3.0_{-0.0}^{+0.0}$ | $0.2^{+0.0}_{-0.0}$ | $3.4_{-0.1}^{+0.1}$ | $1.4_{-0.1}^{+0.1}$ | $-2.6^{+0.1}_{-0.1}$ | $1.9_{-0.2}^{+0.2}$ | $-0.86^{+0.14}_{-0.14}$
NLO: $\chi^{2}/n$ | 1.5 | 0.76 | 1.2 | 0.80 | 1.6 | 1.2 | 1.0
$\sigma$ | $3.0_{-0.0}^{+0.0}$ | $0.3^{+0.0}_{-0.0}$ | $3.2_{-0.1}^{+0.1}$ | $1.4_{-0.1}^{+0.1}$ | $-2.6^{+0.1}_{-0.1}$ | $1.7_{-0.2}^{+0.2}$ | $-1.2^{+0.1}_{-0.1}$
NLO$\times$fakes: $\chi^{2}/n$ | 1.5 | 1.4 | 1.1 | 0.92 | 2.0 | 1.2 | 1.2
$\sigma$ | $2.7_{-0.0}^{+0.0}$ | $-0.3^{+0.0}_{-0.0}$ | $2.9_{-0.1}^{+0.1}$ | $1.2_{-0.1}^{+0.1}$ | $-3.0^{+0.1}_{-0.1}$ | $1.4_{-0.2}^{+0.2}$ | $-1.7^{+0.1}_{-0.1}$
$\mathcal{R}$: $\chi^{2}/n$ | 1.2 | 0.75 | 1.0 | 0.81 | 1.3 | 1.1 | 1.1
$\mathcal{R}$ | $1.2_{-0.5}^{+0.6}$ | $0.2^{+0.2}_{-0.2}$ | $1.1_{-0.4}^{+0.4}$ | $-0.15_{-0.39}^{+0.51}$ | $-0.6^{+0.2}_{-0.2}$ | $0.30_{-0.22}^{+0.29}$ | $-0.091^{+0.084}_{-0.074}$
$\sigma$ | $2.0_{-0.4}^{+0.4}$ | $-0.2^{+0.4}_{-0.4}$ | $1.9_{-0.5}^{+0.5}$ | $1.5_{-0.4}^{+0.4}$ | $-1.3^{+0.4}_{-0.4}$ | $1.2_{-0.4}^{+0.4}$ | $-0.40^{+0.42}_{-0.42}$
$p_{\textnormal{F-test}}$ | 0.021 | 0.20 | 0.0027 | 0.73 | 0.014 | 0.19 | 0.30

Table 2: Results of the fits to the data of all four analyses. In the first
block from the top the signal region is defined and the number of measured
events in this region $N_{m}$ is given. The results of a likelihood- and a
$\chi^{2}$ fit of the fit function (value of the maximal likelihood and
minimal $\chi^{2}$, respectively, and the local significance of the
$750\,\,\text{GeV}$ excess) are given in the second block. Finally, the third
block contains the results of a likelihood and $\chi^{2}$ fit of the
background distributions described in Section 3 to the data. When
$\mathcal{R}$ was fitted its best-fit value and the corresponding local
significance of the excess are given, otherwise just the significance. In the
last line the result of the F-test, testing whether $\mathcal{R}$ should be
used as fit parameter, is given. For the minimized $\chi^{2}$ the parameter
$n$ is the difference of number of bins and fit parameters. The errors
indicate the 1-$\sigma$ interval of the systematic uncertainty of the fit.
Note that the results of the spin-0 analysis with $3.2\,\text{fb}^{-1}$ are
based on the analysis with looser photon identification as described in [1].

Finally, an F-test is performed to determine if the generalization of our
mixed distribution with $w_{\gamma j}^{\textrm{MG}}$ to the one with
$w_{\gamma j}(\mathcal{R})$ given in Eq. (3.5) is needed to describe the data.
This test investigates the improvement of a fit when the fit function is
extended by an additional parameter. For this purpose, a test statistic

$\displaystyle
F=\frac{(\chi^{2}_{1}-\chi^{2}_{2})/(n_{1}-n_{2})}{\chi^{2}_{2}/n_{2}}$ (4.12)

is calculated, where $\chi^{2}_{1,2}$ are the minimized $\chi^{2}$’s of the
two fit functions, $n_{1,2}$ are the numbers of bins (27 for the ATLAS spin-0)
minus the number of input parameters (3 vs. 4 for the fitting function and 1
vs. 2 for our distribution), and the subscripts refer to the two fit functions
with 2 signifying the extended function. Eventually the $p$-value is
determined as

$\displaystyle
p_{\textnormal{F-test}}=\int_{F}^{\infty}\digamma(x;n_{1}-n_{2},n_{2})\textnormal{d}x,$
(4.13)

with $\digamma$ being the Fisher distribution. An additional fit parameter is
warranted if $p_{\textnormal{F-test}}<5\%$, see [1].

We find that for the ATLAS searches, the F-test suggests that $\mathcal{R}$
should be included as a fitting parameter. The probability of an accidental
improvement due to $\mathcal{R}>0$ is only $2.1\%$ (ATLAS spin-0) and $0.27\%$
(ATLAS spin-2). On the other hand, the CMS categories do not prefer a
significant non-zero $\mathcal{R}$, see Table 2. Furthermore, as a consistency
check, we apply the F-test on the ATLAS fitting function for spin-0, Eq.
(4.7), and find that adding a $k=1$ component to the function does not pass
the test. Hence, as mentioned in [1] only the leading term of the function
with $k=0$ is retained. The above in conjunction with the results collected in
Table 2 suggest that it is possible that the basis of functions used in Eq.
(4.7) is not sufficient to accommodate the deformation of the distribution
proposed by us (or at least not the first term in the functional form).

In the $\chi^{2}$ fit of the constructed distribution with $w_{\gamma j}^{MG}$
we find similar results for the local significance as with the $\chi^{2}$ fit
of the functional approach. However, in particular in the two ATLAS analyses,
using $\mathcal{R}$ as an additional fit parameter reduces the local
significance of the 750 GeV excess by 1-1.5 units. The fact that the reduction
is stronger in the spin-2 analysis corroborates our working assumption that
the background description deteriorates in the forward region. This is further
supported by the observation that in the CMS EBEB analysis, which collects
only events with both photon candidates in the central region, no reduction in
significance is observed and the best fit value for $\mathcal{R}$ is even
slightly negative. Only in the CMS EBEE analysis where one photon candidate is
in the forward region the significance is reduced by fitting $\mathcal{R}$,
albeit less than in the ATLAS analyses.

Since $\mathcal{R}>0$ flattens the distribution one might worry that the
reduction in the local significance is obtained by overshooting the measured
distribution in the high invariant mass region. By verifying that both the
minimal $\chi^{2}$ and the maximal likelihood hardly change between the
functional and the distribution fit we show that this is not the case.

As a final exercise, we try to obtain a “combined” significance from the
analyses of the 2015 data set. Clearly a proper statistical combination cannot
be done, since we neglect correlations between the various analyses and also
fit for a single universal value of $\mathcal{R}$. Realistically,
$\mathcal{R}$ is expected to be somewhat different for the different analyses
since they cover different regions of phase space. Nevertheless, since the
naive combination in Eq. (1.1) suffers from similar issues we set them aside
and proceed as follows. We sum the $\chi^{2}$’s of the analyses included in
the combination and fit for a common $\mathcal{R}$ while keeping the
normalizations as separate variables. By combining the two CMS analyses we
obtain $\sigma=2.4$ (1.9) for $w_{\gamma j}^{\textrm{MG}}$ (with $w_{\gamma
j}(\mathcal{R})$, best fit $\mathcal{R}=0.22$) and $\sigma=1.9$ with
$w_{\gamma j}^{\textnormal{NLO}\times\textnormal{fakes}}$. A combination of
the ATLAS analyses is impossible since they are not independent. However, we
can combine each of them with the two CMS analyses and obtain for ATLAS spin-0
combined with CMS $\sigma=3.6$ (2.6 with $w_{\gamma j}(\mathcal{R})$, best fit
$\mathcal{R}=0.46$; 3.1 with $w_{\gamma
j}^{\textnormal{NLO}\times\textnormal{fakes}}$) and for ATLAS spin-2 combined
with CMS $\sigma=4.2$ (2.8 with $w_{\gamma j}(\mathcal{R})$, best fit
$\mathcal{R}=0.53$; 3.4 with $w_{\gamma
j}^{\textnormal{NLO}\times\textnormal{fakes}}$), where the significance
numbers before the brackets are obtained for $w_{\gamma j}^{\textrm{MG}}$.

##  5 The new energy frontier: searches beyond 1 TeV

Around ICHEP 2016, ATLAS and CMS updated their analyses, now based on 15.4
$\text{fb}^{-1}$ and 12.9 $\text{fb}^{-1}$, respectively. In the updated ATLAS
spin-0 analysis [3] and the CMS EBEB and EBEE analyses [4] the large excess
around $750\,\,\text{GeV}$ vanished and no other significant excesses were
found. An update of the ATLAS spin-2 analysis has not been presented. While
CMS processed the data exactly as before, ATLAS made some adjustments, perhaps
most importantly, using a tighter photon isolation. We repeat the fits and the
statistical treatment of the reported results with the same methods as
described above and report the results for the larger dataset in Table 2. Note
that there is a downwards fluctuation in the signal region in the full CMS
dataset which even leads to a slightly negative significance.

Figure 5: Upper plots and lower right plot: Comparison of the measured and
fitted distributions with the small and the full datasets. In addition, the
pure digamma spectrum as obtained from MCFM is shown. In the upper left plot
the comparison is between the smaller Moriond 2016 dataset with the old photon
isolation method and the full ICHEP 2016 dataset. The distributions and
functions are normalized to have the same value at the low $m_{\gamma\gamma}$
end of the histograms.
The lower left plot shows the ratios of the normalized fit functions $f_{X}$
fitted to the ATLAS spin-0 data set $X$ with the year 2015 (2016) in the
brackets indicating the old Moriond (updated ICHEP) photon isolation criteria.

Figure 6: Best fit point of the $\chi^{2}$ fit of the appropriate function to
the Moriond 2016 dataset in red with the 1- and 2-$\sigma$ contours. The best
fit point for the fit to the full ICHEP 2016 dataset is shown in blue and in
the left plot the best fit point for the ATLAS spin-0 $3.2\,\text{fb}^{-1}$
dataset with the new photon isolation is shown in green.

Figure 7: Plot of the normalized distributions with several choices of
$\mathcal{R}$ and normalized fit functions to the old and new datasets, all
divided by the distribution for $\mathcal{R}=0$. In the plot for the ATLAS
spin-0 analysis the dashed lines show the results obtained using the old
dataset with the new photon isolation criteria.

Comparing the new fit functions to the ones based on the previous small
datasets presented at Moriond 2016 we find a steeper functional fit in all
three analyses, see Fig. 5. While the new best fit parameters are within one
standard deviation for the two CMS fits, the ones for the ATLAS fit deviate by
almost two standard deviations after marginalizing over the normalization, see
Fig. 6. This might, however, be an effect of the changed photon isolation as
the fit to the $3.2\,\text{fb}^{-1}$ dataset with the updated photon
identification also deviates by more than one standard deviation from the
previous best fit point. A better understanding of the effect of the fake
photons could be obtained by investigating the result of changing the
isolation criteria with the full $15.4\,\text{fb}^{-1}$ dataset. The tighter
isolation criteria are also reflected in the better agreement between the
fitted distributions and the MCFM generated digamma spectrum.

In order to show the changes in the fits, the ratio of the normalized fit
functions for the ATLAS spin-0 analysis is shown in the lower left plot of
Fig. 5. A direct comparison of the data and the un-normalized fit functions,
even for the fits to the two different $3.2\,\text{fb}^{-1}$ sets, is
difficult since the binning of data has changed. The large change in the fit
parameters is reflected in the deviation of more than 5 % for the comparison
of the fits to the two $3.2\,\text{fb}^{-1}$ datasets and the even greater
deviation when comparing with the fit function to the full
$15.4\,\text{fb}^{-1}$ dataset. While in the previous signal region near 750
GeV the change is of the order of 10 % it is greater than 30 % near 1.6 TeV.
This shows that the actual shape of the digamma spectrum at high invariant
masses is hard to predict precisely by an extrapolation and is therefore very
much subject to systematic uncertainties.

Finally in Fig. 7 the ratios of several normalized distributions and fit
functions to the normalized distribution with $\mathcal{R}=0$ are shown. These
include the distributions with the best fit value for $\mathcal{R}$ based on
the Moriond 2016 dataset and the ICHEP 2016 dataset and also the NLO
distributions and the fit functions to the old and new datasets. In the case
of the ATLAS spin-0 analysis also the distribution and fit function to the
2015 dataset with the new photon isolation is shown. By comparing the curves
we find that a sizable systematic uncertainty can be inferred from the
differences between the fit functions.

##  6 Conclusions

This paper deals with a problem that often arises in searches for new physics
at the energy frontier. In this context the challenge is to look for a new
resonance at the upper end of a distribution where only limited knowledge on
the SM background is available. As a case study we focus on the 750 GeV
anomaly where we examine in particular the implications of the possibility
that the excess in the 2015 data set is not only due to a (malicious)
statistical fluctuation but also a result of a physical effect. We discuss
possible issues with the background: how much photon-jet contamination is
still allowed in the region of interest? How could it affect the significance
of the excess?

We study these questions using currently available theoretical tools for
computing the photon-jet mass distributions and apply them to the small set of
publicly available data. However, this approach is limited by our ability to
thoroughly disentangle the effects of the $(p_{T},\eta)$-dependent jet fake
rate and the theoretical uncertainty of the shape of the photon-jets
background. We therefore choose to model these combined effects by an
$m_{\gamma j}$ dependent reweighting of the invariant mass distribution,
keeping the overall purity within the quoted ranges. We first study a physics-
driven reweighting procedure: we convolve a mass dependent K-factor with a
rapidity and transverse momentum dependent photon fake rate for the jets. The
K-factor is extracted comparing the NLO leading-jet-photon to the LO quark-
photon spectrum, and the phase space dependence of the fake rate is estimated
from the experimental literature [22]. Both correction factors are
approximate, based on incomplete information, and should be taken with a grain
of salt. Motivated by this result, we then consider a more phenomenological
deformation of the $pp\to\gamma j$ spectrum. It allows us to study the
sensitivity of the significances on a single continuous quantity $\mathcal{R}$
(see Eq. (3.5)) which parametrizes an effective deformation.

To summarise our results for the 750 GeV case study based on the 2015 data we
focus on the simpler effective ansatz where we find the following:

  * •

For the ATLAS spin-0 analysis, the significance of the excess can be reduced
by $\Delta\sigma\simeq 1.1$ when comparing the fitting function defined in
(4.7) with our best-fit to the $\mathcal{R}$-modified distribution. A
comparable reduction is found for the ATLAS spin-2 measurement. Here however,
it is less straightforward to determine the reduction, since the estimation of
the background shape in our ansatz differs from that of the ATLAS analysis,
which is not reproducible since the required data is not publicly available.
Strictly speaking, $w_{\gamma j}^{\textrm{MG}}$ does therefore not correspond
to the ATLAS approach but is the best approximation we can get. Since, however
ATLAS claims to find comparable results with the corresponding fitting
function defined in (4.8), we can reduce the significance with the
$\mathcal{R}$-modified distribution with respect to the fitting function by
$\Delta\sigma\simeq 1.0$ as well as with respect to the distribution with
$w_{\gamma j}^{\textrm{MG}}$ by $\Delta\sigma\simeq 1.5$.

  * •

The effect is smaller for the CMS 13 TeV analyses with $\Delta\sigma\simeq
0.3-0.4$, depending on the category.

  * •

In a “combined fit” to independent ATLAS and CMS datasets, the significance
can be reduced by as much as $\Delta\sigma\approx 1.0\,(1.4)$ for the ATLAS
spin-0 (spin-2) combined with CMS.

The larger preference for an enhanced photon-jet contribution in the spin-2
sample could point to its higher sensitivity to the large rapidity region
where jet fakes are more difficult to reject. Finally, an F-test shows that
the ATLAS data support using a more complex distribution.

To summarise our results for the 750 GeV case study based on the 2016 data we
find the following:

  * •

For the ATLAS spin-0 analysis, we find that the new data prefers $\mathcal{R}$
in the range $0.2\pm 0.2$, eliminating the remaining significance of
$1.2\,\sigma$ in the full data set. As for the spin-2 case no data is
currently available.

  * •

The updated CMS analyses based on $12.9\,\text{fb}^{-1}$ even have a downwards
fluctuation with respect to the fit function near 750 GeV leading to negative
significances. Correspondingly the best fit values for $\mathcal{R}$ are
negative and ameliorate the situation.

We emphasize that our simplified ansatz for the effective modification of the
photon-jet background with $\mathcal{R}$ is not meant to necessarily represent
a new background source nor the exact shape of the background contamination in
the signal region. Rather its envelope (corresponding to the shaded area in
Fig. 2(a)) is expected to reflect a possible combination of higher order QCD
contributions, fragmentation, isolation and detector effects, which are
outside of theoretical control and, in the high invariant mass signal region,
also beyond direct experimental probes with currently available data. We
further note that a quark flavor tagged dijet sample might provide a high-
statistics measurement of the relevant photon fake rates (see for instance
[36]).

We have also employed our analysis to compare the difference between the
fitting functions used by ATLAS (with the new isolation criteria) given the
2015 and 2016 data sets. The fitting functions where extrapolated to invariant
masses beyond the TeV region. In summary we have found that:

  * •

A variation of about 30% in the extrapolated background near
$m_{\gamma\gamma}=1.6\,\,\text{TeV}$ is obtained.

To conclude, we have extensively examined the status of LHC diphoton searches.
We have compared the analyses performed on both 2015 and 2016 datasets in
order to scrutinize the current state of the art measurements for systematic
effects. Using our approach we have reevaluated the current experimental
sensitivity to beyond standard model physics, especially in the tails of the
diphoton invariant mass distributions, beyond the TeV range. We found that the
extrapolation of background shapes is subject to sizeable uncertainties,
potentially affecting the significance of possible future excesses near the
edge of the measured distributions. Furthermore, our analysis motivates
further Monte Carlo studies of the dominant diphoton backgrounds, based on jet
flavor tagging algorithms. Knowledge of whether a jet is of “quark” or “gluon”
origin would improve our estimation for the jet-photon faking backgrounds to
next to leading order QCD accuracy. It is important to note that diphoton-
based searches at even larger invariant masses, that are highly motivated, are
being performed at present and will continue to be an integral part of the LHC
experimental physics program at the high energy frontier.

## Acknowledgments

We would like to thank Rikkert Frederix for useful discussions. J.F.K. would
like to thank CERN for hospitality while this work was being completed and
acknowledges the financial support from the Slovenian Research Agency
(research core funding No. P1-0035). The work of GP is supported by grants
from the BSF, ISF and ERC and the Weizmann-UK Making Connections Programme. AW
is supported by the DFG cluster of excellence “Origin and Structure of the
Universe” and the European Commission (AMVA4NewPhysics, 2020-MSCA-ITN-2015).

## References

  * [1] ATLAS collaboration, M. Aaboud et al., _Search for resonances in diphoton events at $\sqrt{s}$=13 TeV with the ATLAS detector_, 1606.03833.
  * [2] CMS collaboration, V. Khachatryan et al., _Search for resonant production of high-mass photon pairs in proton-proton collisions at sqrt(s) = 8 and 13 TeV_ , 1606.04093.
  * [3] ATLAS collaboration, T. A. collaboration, _Search for scalar diphoton resonances with 15.4 fb -1 of data collected at $\sqrt{s}$=13 TeV in 2015 and 2016 with the ATLAS detector_, .
  * [4] CMS collaboration, C. Collaboration, _Search for resonant production of high mass photon pairs using $12.9\,\mathrm{fb^{-1}}$ of proton-proton collisions at $\sqrt{s}=13\leavevmode\nobreak\ \mathrm{TeV}$ and combined interpretation of searches at 8 and 13 TeV_, .
  * [5] J. H. Davis, M. Fairbairn, J. Heal and P. Tunney, _The Significance of the 750 GeV Fluctuation in the ATLAS Run 2 Diphoton Data_ , 1601.03153.
  * [6] K. Bondarenko, A. Boyarsky, O. Ruchayskiy and M. Shaposhnikov, _Features in the Standard Model diphoton background_ , 1606.0959.
  * [7] S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini, _Diphoton production at hadron colliders: a fully-differential QCD calculation at NNLO_ , _Phys. Rev. Lett._ 108 (2012) 072001, [1110.2375].
  * [8] J. M. Campbell, R. K. Ellis, Y. Li and C. Williams, _Predictions for diphoton production at the LHC through NNLO in QCD_ , 1603.02663.
  * [9] ATLAS collaboration, G. Aad et al., _Measurement of the inclusive isolated prompt photon cross section in $pp$ collisions at $\sqrt{s}=8$ TeV with the ATLAS detector_, 1605.03495.
  * [10] S. Catani, M. Fontannaz, J. P. Guillet and E. Pilon, _Cross-section of isolated prompt photons in hadron hadron collisions_ , _JHEP_ 05 (2002) 028, [hep-ph/0204023].
  * [11] T. Becher, G. Bell, C. Lorentzen and S. Marti, _Transverse-momentum spectra of electroweak bosons near threshold at NNLO_ , _JHEP_ 02 (2014) 004, [1309.3245].
  * [12] T. Becher, C. Lorentzen and M. D. Schwartz, _Precision Direct Photon and W-Boson Spectra at High $p_{T}$ and Comparison to LHC Data_, _Phys. Rev._ D86 (2012) 054026, [1206.6115].
  * [13] T. Becher and X. Garcia i Tormo, _Addendum: Electroweak Sudakov effects in W, Z and gamma production at large transverse momentum_ , _Phys. Rev._ D92 (2015) 073011, [1509.01961].
  * [14] M. D. Schwartz, _Precision direct photon spectra at high energy and comparison to the 8 TeV ATLAS data_ , 1606.02313.
  * [15] CMS collaboration, V. Khachatryan et al., _Comparison of the Z/ $\gamma$? \+ jets to $\gamma$ \+ jets cross sections in pp collisions at $\sqrt{s}=8$ TeV_, _JHEP_ 10 (2015) 128, [1505.06520].
  * [16] ATLAS collaboration, G. Aad et al., _Dynamics of isolated-photon plus jet production in pp collisions at $\sqrt{(}s)=7$ TeV with the ATLAS detector_, _Nucl. Phys._ B875 (2013) 483–535, [1307.6795].
  * [17] A. Banfi, G. P. Salam and G. Zanderighi, _Infrared safe definition of jet flavor_ , _Eur. Phys. J._ C47 (2006) 113–124, [hep-ph/0601139].
  * [18] A. Buckley and C. Pollard, _QCD-aware partonic jet clustering for truth-jet flavour labelling_ , _Eur. Phys. J._ C76 (2016) 71, [1507.00508].
  * [19] A. J. Larkoski, J. Thaler and W. J. Waalewijn, _Gaining (Mutual) Information about Quark/Gluon Discrimination_ , _JHEP_ 11 (2014) 129, [1408.3122].
  * [20] C. W. Bauer and E. Mereghetti, _Heavy Quark Fragmenting Jet Functions_ , _JHEP_ 04 (2014) 051, [1312.5605].
  * [21] Particle Data Group collaboration, K. A. Olive et al., _Review of Particle Physics_ , _Chin. Phys._ C38 (2014) 090001.
  * [22] ATLAS, _Expected photon performance in the ATLAS experiment_ , Tech. Rep. ATL-PHYS-PUB-2011-007, CERN, Geneva, Apr, 2011.
  * [23] R. B. Neufeld, I. Vitev and B. W. Zhang, _The Physics of $Z^{0}/\gamma^{*}$-tagged jets at the LHC_, _Phys. Rev._ C83 (2011) 034902, [1006.2389].
  * [24] J. M. Campbell and R. K. Ellis, _An Update on vector boson pair production at hadron colliders_ , _Phys. Rev._ D60 (1999) 113006, [hep-ph/9905386].
  * [25] J. M. Campbell, R. K. Ellis and C. Williams, _Vector boson pair production at the LHC_ , _JHEP_ 07 (2011) 018, [1105.0020].
  * [26] J. M. Campbell, R. K. Ellis and W. T. Giele, _A Multi-Threaded Version of MCFM_ , _Eur. Phys. J._ C75 (2015) 246, [1503.06182].
  * [27] R. Boughezal, J. M. Campbell, R. K. Ellis, C. Focke, W. Giele, X. Liu et al., _Color singlet production at NNLO in MCFM_ , 1605.08011.
  * [28] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ , _JHEP_ 07 (2014) 079, [1405.0301].
  * [29] T. Sjostrand, S. Mrenna and P. Z. Skands, _PYTHIA 6.4 Physics and Manual_ , _JHEP_ 05 (2006) 026, [hep-ph/0603175].
  * [30] M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J._ C72 (2012) 1896, [1111.6097].
  * [31] M. Cacciari, G. P. Salam and G. Soyez, _The Anti-k(t) jet clustering algorithm_ , _JHEP_ 04 (2008) 063, [0802.1189].
  * [32] ATLAS collaboration, G. Aad et al., _Measurement of isolated-photon pair production in $pp$ collisions at $\sqrt{s}=7$ TeV with the ATLAS detector_, _JHEP_ 01 (2013) 086, [1211.1913].
  * [33] ATLAS collaboration, G. Aad et al., _Measurement of the isolated di-photon cross-section in $pp$ collisions at $\sqrt{s}=7$ TeV with the ATLAS detector_, _Phys. Rev._ D85 (2012) 012003, [1107.0581].
  * [34] CMS collaboration, S. Chatrchyan et al., _Measurement of differential cross sections for the production of a pair of isolated photons in pp collisions at $\sqrt{s}=7\,\text{TeV}$_, _Eur. Phys. J._ C74 (2014) 3129, [1405.7225].
  * [35] G. Cowan, K. Cranmer, E. Gross and O. Vitells, _Asymptotic formulae for likelihood-based tests of new physics_ , _Eur. Phys. J._ C71 (2011) 1554, [1007.1727].
  * [36] ZEUS collaboration, S. Chekanov et al., _Substructure dependence of jet cross sections at HERA and determination of alpha(s)_ , _Nucl. Phys._ B700 (2004) 3–50, [hep-ex/0405065].


# Witness-based Approach for Scaling Distributed Ledgers to Massive IoT
Scenarios

Duc-Lam Nguyen, Israel Leyva-Mayorga, and Petar Popovski
Connectivity Section, Department of Electronic Systems, Aalborg University
Aalborg, Denmark
Email: {ndl, ilm<EMAIL_ADDRESS>
###### Abstract

Distributed Ledger Technologies (DLTs) are playing a major role in building
security and trust in Internet of Things (IoT) systems. However, IoT
deployments with a large number of devices, such as in environment monitoring
applications, generate and send massive amounts of data. This would generate
vast number of transactions that must be processed within the distributed
ledger. In this work, we first demonstrate that the Proof of Work (PoW)
blockchain fails to scale in a sizable IoT connectivity infrastructure. To
solve this problem, we present a lightweight distributed ledger scheme to
integrate PoW blockchain into IoT. In our scheme, we classify transactions
into two types: 1) global transactions, which must be processed by global
blockchain nodes and 2) local transactions, which can be processed locally by
entities called witnesses. Performance evaluation demonstrates that our
proposed scheme improves the scalability of integrated blockchain and IoT
monitoring systems by processing a fraction of the transactions, inversely
proportional to the number of witnesses, locally. Hence, reducing the number
of global transactions.

###### Index Terms:

Distributed Ledgers, Blockchain, IoT, Witness, Environment Monitoring,
scalability.

##  I Introduction

Distributed Ledger Technologies (DLTs) provide high levels of security,
accountability, tractability, and privacy to the transmitted data [1]. This is
achieved by enabling key functionalities, such as transparency, distributed
operation, and immutability [2]. The benefits of DLTs are particularly
appealing for Internet of Things (IoT) applications, where large amounts of
data are generated and the devices can only implement weak security mechanisms
[3].

The trust provided by DLTs is greatly valuable in IoT monitoring applications
with a large number of devices. As an example, consider an urban IoT
application that monitors the air quality and gas emissions. The data
generated by this application is critical, so it must be protected, tractable,
immutable, and transparent. Nevertheless, in a traditional monitoring system,
the inter-organization sharing the data may be untrusted, complex, unreliable,
and non-transparent. Besides, the current IoT-based monitoring systems are
centralized, which leads to a single point of failure, where data can be lost
or modified [4].

Figure 1: Overview of our Blockchain-enabled IoT system wiBlock. The IoT nodes
generate and send the transactions to the base stations, which in turn send
them to the witnesses. These decide which transactions must be sent to the GB
and process the rest.

The problems described above may be solved by integrating Blockchain into IoT
applications. However, Blockchain architectures were not designed to handle a
large number of transactions, which would be generated by naively integrating
Blockchain into IoT. Specifically, IoT deployments usually present a star
topology, in which the devices communicate directly to the base stations (BS),
which then redirects the gathered data to the destination [5] (e.g., from
Narrowband IoT (NB-IoT) or LoRa deployments to a cloud server), as shown in
the left part of Fig. 1. In the most Blockchain and IoT integration, this same
architecture would be used, and the BS would be in charge of communicating
with the Blockchain [6]. Thus, every packet generated by the IoT devices would
represent a transaction, which can easily overload the Blockchain.

Three main challenges must be overcome to achieve an efficient integration of
Blockchain into IoT. First, DLTs use diverse resource-intensive gating
functions, for example, Proof-of-Work (PoW) and Proof-of-Stake (PoS), while
IoT devices are resource-constrained. As a consequence, the processing time of
these functions in IoT devices would be restrictive. Second, the widely-used
Blockchain arrangement cannot handle the massive transactions generated by IoT
devices. For example, Bitcoin network produces $1$ MB blocks, roughly once
every 10 minutes, with an average size of transaction around 500 bytes, which
give 7 transactions per second (tps). In comparison, Visa system performs
$2000$ tps on average, and an average daily peak of $4000$ tps, with a maximum
capacity of $56000$ tps. Third, the power saving mechanisms of the IoT devices
can cause problems during knowledge dissemination and synchronization. For
instance, an update may be severely delayed or even fail to arrive if a device
is in sleep mode.

In this paper, we present a witness-based Blockchain system called wiBlock,
especially designed to integrate Blockchain into resource-constrained IoT
applications. It is aimed to solve three of the main problems of traditional
IoT monitoring systems, namely trust, scalability, and cost. This is achieved
by: 1) enabling the use of DLTs to store IoT data, 2) limiting the number of
transactions that must be processed at the Global Blockchain (GB), and 3)
eliminating the need for complex computations and supporting sleep-awake
mechanisms at the IoT devices, respectively.

The architecture of wiBlock is illustrated in Fig. 1, where the IoT devices
interact exclusively with the witness system, which then may process the
transactions locally or communicate directly with the GB. The transactions
that must be processed by the GB are called global transactions, whereas the
transactions that can be verified locally at the witness system are called
local transactions. In order to see the need for this differentiation,
consider a pollution monitoring system, in which a number of sensors in a
given local area are associated to the same witness. Then a local transaction
can be used to send local sensing data from a device associated with the same
witness. For instance, the alarm sensor periodically requests gas sensor which
collects the concentration of pollutants e.g., SO2, CO2, NO to detect the
abnormal condition in air. In order to see the need for a global transaction,
note that sensors may wish to store their sensing data to external storage
system e.g, IPFS [7] or control a thermostat sensor, which is located in a
different area and associated with a different witness to adapt temperature.
In this case there is a need to communicate via different heterogeneous
networks and record the transaction results to the GB via global transactions.
Thus, the witness system reduces the number of transactions that need to be
processed by the GB and the latency of transaction verification. Furthermore,
wiBlock allows each IoT device to communicate with several witnesses. This
avoids having a single point of failure (i.e., bridge) between the IoT device
and the GB, which in turn greatly increases the reliability of the IoT
application. For example, Blockchain witness models have been found to be
beneficial for Cloud Service Level Agreement [8].

The contributions of this work are as follows:

  1. 1.

We investigate the possibilities of naively integrating Blockchain directly
into resource-constrained IoT systems. We identify some of the major problems
that arise in this setup, which illustrate that Blockchain technology is not
directly applicable to massive IoT.

  2. 2.

We propose a new IoT-friendly distributed ledger system named wiBlock. It aims
to solve the scalability issues of Blockchain in massive IoT environment by
defining two types of transactions: global and local.

  3. 3.

We thoroughly compare the performance _wiBlock_ with that of a naive
Blockchain and IoT integrated architecture. Our results show that our proposed
system enhances the scalability of the GB network.

The remainder of this paper is organized as follows. In Section II, we present
the system model, followed by the design of our novel wiBlock system in
Section III. We present the analysis and performance evaluation of wiBlock in
Section IV and Section V, respectively. Finally, we conclude the paper in
Section VI.

##  II System model

We consider an IoT application with $k$ devices. These are deployed uniformly
at random in a squared area of interest $A\in\mathbb{R}^{2}$. The IoT devices
generate transactions with the data collected from the environment according
to a Poisson process with rate $\lambda$.

In the most simple Blockchain and IoT integrated architecture, the
transactions are sent to the BS, which then redirects them to the GB. In
_wiBlock_ , the transactions are sent to the _witness system_ instead. This is
a set of $v$ witnesses, which have the capacity to verify transactions locally
and to communicate with the GB. The time needed for a witness to perform these
operations determine its capacity and depend on numerous factors. However, it
is out of the scope of this paper to derive their precise values. Transactions
are grouped into blocks of size $b$. Therefore, a new block is created when
$b$ new transactions are received at a server.

Witnesses may be either physical or logical entities, hence, their
organization is flexible. For simplicity, throughout this paper we assume one
witness is deployed at each BS and use these terms interchangeably. The BSs
are distributed randomly within $A$. We denote the set of IoT devices and
witnesses as $\mathcal{D}=\\{1,2,\dotsc,k\\}$ and
$\mathcal{W}=\\{1,2,\dotsc,v\\}$, respectively.

The IoT devices and witnesses communicate through wireless links under a
standard path loss model and large-scale (slow) fading. Thus, a transaction is
transmitted successfully from IoT device $i$ to a witness $w\in\mathcal{W}$
with probability $p_{s}(i,w)$. The IoT device $i$ selects the witness $w$
according to a predefined strategy. If the transmission fails, $i$ attempts
the transmission to a different witness. This process is repeated until the
transaction is confirmed or until a given number of attempts is reached
without success.

We consider a simple shadowing propagation model for the communication between
IoT devices and witnesses where, for a given transmission power $P_{t}$ and
carrier frequency $f$, the received power at a distance $d$ is

$P_{r}(d)=10\log_{10}\left(\frac{P_{t}G_{t}G_{r}\,c^{2}}{(4\pi
f)^{2}d^{\beta}}\right)+N(0,\sigma_{\text{dB}})~{}\text{dB}$ (1)

where $G_{t}$ and $G_{r}$ are the transmitter and receiver antenna gains,
respectively, $c=3\cdot 10^{8}$ m/s is the speed of light,
$N(0,\sigma_{\text{dB}})$ is a zero-mean Gaussian random variable (RV) with
standard deviation $\sigma_{\text{dB}}$ dB, and $\beta$ is the path loss
exponent.

From there, the outage probability at a given distance and receiver
sensitivity $\gamma$ is

$p_{\text{out}}(d)=1-Q\left(\frac{1}{\sigma_{\text{dB}}}10\log_{10}\left(\frac{\gamma(4\pi
f)^{2}d^{\beta}}{P_{t}G_{t}G_{r}c^{2}}\right)\right)$ (2)

and $p_{s}(i,w)=1-p_{\text{out}}(d(i,w))$. Throughout this paper, we assume
that the wireless resources are sufficient to support the communication
between the IoT devices and the witness system and do not go into the details
of the access protocols. Therefore, collisions caused by simultaneous
transmissions from the IoT devices to a witness $w$ can be avoided or resolved
if the links toward $w$ are not in outage. Finally, no errors occur in the
communication between the witness system and the GB.

##  III WiBlock Design

This section presents the detailed description of wiBlock architectural
elements and operation.

###  III-A Witness-based Blockchain System

As illustrated in Fig. 1, the witness-based Blockchain System consists of
three main components: the GB, the witness system, and the physical IoT
devices. The first action performed by the IoT devices after deployment is
authentication. For this, each device $i\in\mathcal{D}$ performs a key
exchange procedure with a witness $w\in\mathcal{W}$ to gain the necessary
permissions and build secure channels to perform transactions. After
authentication, the tuple $(i,w)$ is added by $w$ to the shared registry of
the witness system $\mathcal{R}$. For this, $w$ shares the authentication
information of $i$ (i.e., credentials) with the rest of the witnesses. By
keeping a shared registry, device $i$ can communicate with any witness, even
though is registered with $w$. After authentication, IoT devices collect the
data and sign it by using a SecretKey $s_{\text{key}}(i)$ that is unique for
each $i$ as Sign(data, $s_{\text{key}}(i)$, timestamp $\tau$). Next, the
transaction is created and transmitted to a witness $w$. Note that this latter
witness may be different to the one which $i$ is registered with. Local
transactions, denoted as $L_{l}$, are exchanged exclusively between $w$ and
all the IoT devices registered with it
$\\{i\in\mathcal{D}:(i,w)\in\mathcal{R}\\}$, for which the managers implement
a consensus procedure. On the other hand, Global transactions, denoted as
$L_{G}$, must be sent from a witness $w$ to the GB when
$(i,w)\notin\mathcal{R}$. These two types of transactions are further
described in the following.

Figure 2: Transaction flow in wiBlock, from generation to confirmation.

####  III-A1 Local Transactions

These transactions are transmitted from IoT device $i$ to the witness $w$,
whose key management component confirms that $(i,w)\in\mathcal{R}$. Then, this
same component checks whether the PublicKey $p_{\text{key}}(i)$ of $i$ has
been associated with any block in the local ledger. If $p_{\text{key}}(i)$ has
not been associated with any block, the witness $w$ generates a new block for
the given $i$. Then, $w$ arranges the transactions in order, updates the local
ledger, and a notification feedback message is transmitted to the devices.

####  III-A2 Global Transactions

These transactions are transmitted from IoT device $i$ to the witness $w$,
whose key management component confirms that $(i,w)\notin\mathcal{R}$. Then,
this same component will clarify which witness $i$ is registered with. If
$\exists w^{\prime}\in\mathcal{W}$ s.t. $(i,w^{\prime})\in\mathcal{R}$, the
transaction is forwarded to the GB. In case the GB has a block associated with
given device $i$, the transaction will be validated based on the corresponding
signature Sign(data, $s_{\text{key}}(i)$, timestamp $\tau$) and, if the
signature is valid, the transaction is appended to the block and transmitted
back to the witness $w^{\prime}$. Note that this type of transactions will be
frequently generated when the IoT devices are mobile. For example, cargo,
supply chain, and car subsystem monitoring.

###  III-B Witness Selection

123$\cdots$vWitness system $\mathcal{W}$123$\cdots$kSet of IoT devices
$\mathcal{D}$ Figure 3: In _wiBlock_ , each IoT device has a list of eligible
witnesses. Transactions generated by the IoT devices are sent to a witness in
this list, according to the witness selection strategy.

Numerous witness selection strategies can be implemented at the IoT devices
and each one may offer different benefits. However, the focus of the present
work is to evaluate the benefits of the witness-based architecture, rather
than to identify an optimal witness selection strategy. Therefore, we consider
the following a heuristic witness selection strategies and evaluate the
performance of the witness system. As illustrated in Fig. 3, IoT devices
select one of the $v$ available witnesses with probability $1/v$ and transmit
the transaction. Then, if the link between IoT device $i$ and witness $w$ is
not in outage, the transaction is confirmed. Otherwise, $i$ selects a new
witness uniformly at random from $\mathcal{W}\setminus w$ and transmits the
transaction. This process is repeated until the transaction is confirmed or
until a given number of attempts $l\leq v$ is reached without success. This is
the simplest strategy and assumes the IoT devices have no information about
the state of the wireless channel toward each witness separately.

##  IV Analysis

###  IV-A Queuing model of the witness system

We consider a queuing model for witness-based Blockchain network as described
in Fig. 4. The witnesses and Blockchain are modelled as queuing nodes to
capture the number of transactions that must be i) processed locally by
witnesses and ii) forwarded to the GB to be processed. We assume that
transactions are generated by the IoT devices following a Poisson process.
Hence, we denote $\lambda(i)$ as the transaction generation rate at IoT device
$i$.

Let $p(i,w)$ be the probability that $i$ chooses witness $w$ and $p_{s}(i,w)$
be the probability that the link between $i$ and $w$ is not in outage.
Building on this, the average transaction arrival rate at the witness $w$ is

$\lambda_{w}=\sum_{i=1}^{k}p(i,w)p_{s}(i,w)\lambda(i).$ (3)

Hence, the transaction arrival rate of different witnesses depends on the
density and location of the deployed IoT devices and witnesses, but also on
the witness selection criteria.

The probability $p(i,w)$ depends on the witness selection strategy. For the
strategy 1, random selection, let $A(w,u)$ be the matrix of permutations of
$u$ elements taken from
$\\{1-p_{s}(i,w^{\prime})\\}_{w^{\prime}\in\mathcal{W}\setminus w}$ with
${}_{(v-1)}P_{u}$ rows and $u$ columns. The element in row $x$ and column
$y\leq u$ of $A(w,u)$ is denoted $a_{xy}(w,u)$. From there, we can calculate
$p(i,w)$ as:

$p(i,w)=\frac{1}{v}+\frac{1}{v!}\sum_{u=1}^{l-1}(v-u-1)!\sum_{x=1}^{{}_{(v-1)}P_{u}}\prod_{y=1}^{u}a_{xy}(w,u)$
(4)

Figure 4:  Witness-based Blockchain queuing model described in Section IV.

As mentioned above, generated transactions are either global $L_{G}$ or local
$L_{l}$. We define $p$ as the probability that a transaction sent to a witness
is Global. Hence, $1-p$ is the probability that a transaction is local. Please
observe that the value of $p$ only depends on the number of witnesses $v$ and
is $p:=\Pr\left[(i,w)\notin\mathcal{R}\right]=(v-1)/v$.

The transaction processing time is assumed to follow an exponential
distribution with service rates $\mu_{1}$ and $\mu_{2}$ for global and local
transactions, respectively, and transactions are served according to a first-
come first-served (FCFS) policy. Building on this, we model the operation of
each witness as an M/H2/1 queue, which means that transactions arrive at the
witness $w$ at a rate $\lambda_{w}$ and the service time is represented by a
two-phase hyper-exponential distribution. With probability $p$, the first
transaction in the queue receives service at rate $\mu_{1}$, while with
probability $1-p$, it receives service rate at rate $\mu_{2}$. That is, the
type of transaction is defined at the beginning of service.

The state of each witness is represented by a pair $(m,n)$, in which $m$ is
the total number of transactions in the witness and $n\in\\{1,2\\}$ is the
current service phase, which depends on the type of transaction being served.
The stationary distribution of this queue in the witness $w$ can be obtained
by Neuts’ Matrix Geometric Method [9]. We denote the stationary probability
vector as:

$\bm{\tau}^{(w)}=\left[\tau_{0}^{(w)},\tau_{1}^{(w)},\tau_{2}^{(w)},\dotsc,\tau_{k}^{(w)},\dotsc\right],$
(5)

where $\tau_{m}^{(w)}$ is the steady-state probability of $m$ transactions in
the witness $w$. Alternatively, the mean service rate is

$\mu=\left(\frac{p}{\mu_{1}}+\frac{1-p}{\mu_{2}}\right)^{-1},$ (6)

and the offered load to $w$ is $\rho_{w}=\lambda_{w}/\mu$. From there, the we
calculate the variance of the service time

$\sigma^{2}_{w}=2\left(\frac{p}{\mu_{1}^{2}}+\frac{1-p}{\mu_{2}^{2}}\right)-\frac{1}{\mu^{2}}$
(7)

and the coefficient of variation $C_{w}^{2}=\mu^{2}\sigma_{w}^{2}$. Then, the
average number of transactions in the queue of $w$ is

$L(w)=\sum_{m=0}^{\infty}m\tau_{m}^{(j)}=\rho_{w}+\left(\frac{1+C_{w}^{2}}{2}\right)\frac{\rho_{w}^{2}}{1-\rho_{w}}.$
(8)

Then, the number of local transactions and Global transactions handled by $w$
are, respectively,

$L_{g}(w)=pL(w)$ (9)

and

$L_{l}(w)=(1-p)L(w)=L(w)-L_{g}(w).$ (10)

###  IV-B Global Blockchain (GB) System

We model the GB as a modified $M/G^{B}/1$ queue as in [10]. Let $L_{g}$ and
$T_{g}$ be the RVs that define the number of transactions in the Blockchain
queue and the confirmation time. We are interested in finding their mean
values. For this, we define $b$ to be the maximum number of transactions in a
block (i.e., the maximum block size). Hence, transactions are grouped into
blocks and a new block is created when there are $b$ transactions in the
Blockchain server.

Given that $p$ is the probability that a transaction sent to a witness is
processed at the GB, the transaction arrival rate at the GB from the $v$
witnesses in the IoT deployment is

$\lambda_{B}=\sum_{w=1}^{v}\lambda_{w}p.$ (11)

We denote $U$ as of the block generation time (i.e., the time it takes to
generate a block) at the GB. Then,

We define $U$ to be the continuous RV of the processing (i.e., service) time
of a block at the GB. Hence, the system is stable and a limiting probability
exists if and only if $\lambda_{B}E[U]<b$.

The cumulative distribution function (CDF) and the probability density
function (pdf) of $U$ are denoted $G(x)$ and $g(x)$, respectively. We use
these to calculate the hazard rate of $U$ as

$\theta(x)=\frac{g(x)}{1-G(x)}.$ (12)

Next, we define $L_{g}^{s}(t)$ as the number of transactions in server at time
$t$, $L_{g}^{q}(t)$ as the number of transaction in the queue at time $t$, and
$X(t)$ as the elapsed service time of the current transaction at $t$. From
[11], we define

$P_{m,n}(x,t)dx=\Pr\left[L_{g}^{s}(t)=m,L_{g}^{q}(t)=n,\right.\\\
\left.x<X(t)\leq x+dx\right]$ (13)

to be the joint probability that, at time $t\geq 0$, there are
$m\in\\{0,1,2,...,b\\}$ and $n\in\\{0,1,2,...,x\\}$ transactions in server and
queue, respectively, and the elapsed service time lies between $x$ and $x+dx$.
Next, we denote $P_{m,n}(x)=\lim_{t\xrightarrow{}\infty}P_{m,n}(x,t)$ and
consider the two following cases. In the first one we have
$\frac{d}{dx}P_{m,n}(x)=-\left[\lambda_{B}+\theta(x)\right]P_{m,n}(x)+\lambda_{B}P_{m,n-1}(x),\text{
for }0\leq m\leq b\text{ and }n\geq 1,$ which shows that the number of
transactions in the server and the queue does not change during a small
interval. In the second one we have
$\frac{d}{dx}P_{m,0}(x)=-\left[\lambda_{B}+\theta(x)\right]P_{m,0}(x),\quad\text{for
}0\leq m\leq b,$ which occurs when a transaction arrives at the system with
$0$ transactions in the queue. For the purposes of our study, it is sufficient
to calculate the mean confirmation time as

$\displaystyle E[T_{g}]\\!$ $\displaystyle=$
$\displaystyle\\!\Bigg{[}\lambda_{B}^{2}E[U^{2}]-b(b-1)-2\left(b-\lambda_{B}E[U^{2}]\right)$
$\displaystyle+\sum_{n=0}^{b-1}\alpha_{n}\Big{(}\lambda_{B}E[U^{2}](b-n)+2bE[U](b-n)$
$\displaystyle+E[U]\left(b^{2}-b-n^{2}+n\right)\\!\Big{)}\\!\Bigg{]}\frac{1}{2\lambda_{B}(b-\lambda_{B}E[U])},$

where
$\alpha_{n}=\sum_{m=0}^{b}\int\limits_{0}^{\infty}P_{m,n}(x)\theta(x)dx$. The
interested reader is referred to [10] for the fully detailed Blockchain
queuing model.

##  V Performance Evaluation

In this section, we use the queuing models described in Section IV to evaluate
the performance of _wiBlock_ in terms of scalability. For this, we obtain the
maximum transaction generation rate, along with the mean confirmation time and
ledger size for both, the local and global Blockchain. We use the performance
of a naive Blockchain and IoT integrated architecture, where the IoT devices
communicate directly to the GB, as a benchmark. The mean results regarding the
connectivity of the IoT devices with the witness system are obtained by a
large number of Monte Carlo simulations and then used as an input to the
queuing models.

In our analysis, each device generates transactions at a rate
$\lambda(i)=\lambda$ for all $i\in\mathcal{D}$. The block generation time $U$
is exponentially distributed with parameter $\mu_{B}=1.8\cdot 10^{-3}$ blocks
per second. So we have $g(x)=\mu_{B}\exp(-\mu_{B}x)$,
$E\left[U\right]=1/\mu_{B}$, and $E\left[U^{2}\right]=1/\mu_{B}^{2}$.
Furthermore, we define the default block size to be $b=1000$ transactions. The
rest of relevant parameters are listed in Table I; these values are used
unless otherwise stated.

For the selected parameter settings, the GB system is stable when
$\lambda_{B}^{*}=b/E\left[U\right]=1.8$. Building on this, from (11) we have
that $\lambda<b/(kE\left[U\right])=1.8/k$ must hold for the GB to be stable in
a traditional Blockchain architecture with $k$ identical IoT devices.
Conversely, for _wiBlock_ with random witness selection, we have that only a
fraction $p=(v-1)/v$ of the transactions must be sent to the GB. Hence,
assuming no wireless channel errors occur and all the generated transactions
are sent to a witness (i.e., $\sum_{w=1}^{v}\lambda_{w}=k\lambda$), the
maximum load per IoT device that _wiBlock_ can handle is

$\lambda<\frac{bv}{kE\left[U\right](v-1)}=\frac{1.8v}{k(v-1)}=\lambda^{*}(v),$
(15)

as shown in Fig. 5 for $v=\\{2,4,8\\}$.

TABLE I: Parameter settings for the performance evaluation. Parameter | Symbol | Value
---|---|---
Area of deployment | $A$ | $100\times 100~{}\mathrm{m}^{2}$
Number of IoT devices | $k$ | $500$
Number of witnesses | $v$ | $\\{2,3,...,10\\}$
Carrier frequency | $f$ |  $914$ MHz
Transmission power | $P_{t}$ |  $0.28183815$ W
Antenna gains | $G_{t},G_{r}$ | $1$
Receiver sensitivity | $\gamma$ |  $3.652\cdot 10^{-10}$ W
Standard deviation of shadow fading | $\sigma_{\text{dB}}$ |  $6$ dB
Path loss exponent | $\beta$ | 3
Block size | $b$ |  $1000$ transactions
Figure 5: Maximum transaction generation rate per IoT device $\lambda^{*}$ for
traditional Blockchain IoT and wiBlock with random witness selection. Figure
6: Ratio of transactions processed at the GB and at the witness system as a
function of the number of witnesses $v$.

Hence, from the GB perspective, _wiBlock_ allows to deploy $1/p=v/(v-1)$ times
more IoT devices than the naively integrated approach, as illustrated by Fig.
6. Note that the greatest gains in the scalability are obtained when $v$ is
small, however, other factors such as the area coverage and processing
capacity of the witness system must be taken onto account to select adequate
values of $v$.

Next, we evaluate the mean transaction confirmation time at the GB
$E\left[T_{g}\right]$. Note that, in case a single witness is deployed in the
system, all the transactions generated by the IoT devices will be considered
as local transactions and processed locally. This can overload the witness,
depending on its capabilities. In particular, the witness is stable if and
only if the load offered to the witness is $\lambda_{1}<\mu_{2}$. Furthermore,
deploying a single witness does not provide the necessary wireless coverage.
That is, the more witnesses are deployed, the higher the probability of being
able to communicate to, at least, one of them. Hence, we consider the cases
where at least two witnesses are deployed, as shown in Fig. 7a for
$v=\\{2,3,4\\}$.

(a)

(b)

Figure 7: (a) Mean transaction confirmation time $E\left[T_{g}\right]$ as a
function of the transaction generation rate at the IoT devices $\lambda$ and
(b) ledger size at the GB and at each witness for $v=2$ as a function of the
block size $b$ for traditional Blockchain IoT and wiBlock.

As Fig. 7a shows, the witness system reduces the number of transactions sent
to the GB and, as a consequence, greatly reduces the transaction confirmation
time. Besides, the ledger size is considerably reduced, depending on the
number of witnesses. This can be seen in Fig. 7b for $v=2$, where the ledger
size of the GB is half of that with the traditional Blockchain and IoT
integration, and the local ledger size at each witness is $1/v^{2}=1/4$ of it.

##  VI Conclusion

In this paper, we presented and evaluated the performance of a novel witness-
based Blockchain system for IoT applications. As a starting point, we
described the benefits of integrating Blockchain into IoT and the main
challenges that must be overcome to achieve this integration. Building on
these, we designed _wiBlock_ , an IoT-friendly distributed system that
incorporates a witness system to address scalability issues of Blockchain. The
scalability gains provided by _wiBlock_ are achieved by processing some of the
transactions generated by the IoT devices locally, at the witness system. Our
results show that the witness system greatly reduces the number of
transactions transmitted to the Blockchain network and the transaction
confirmation time. Future work includes the design of witness selection
algorithms and implement wiBlock in real testbed to further exploit the
benefits provided by the witness system.

##  VII Acknowledgement

This work has been in part supported by the European Research Council (Horizon
2020 ERC Consolidator Grant Nr.648382 WILLOW).

## References

  * [1] P. Danzi, A. E. Kalor, R. B. Sorensen, A. K. Hagelskjær, L. D. Nguyen, C. Stefanovic, and P. Popovski, “Communication aspects of the integration of wireless iot devices with distributed ledger technology,” _IEEE Network_ , vol. 34, no. 1, pp. 47–53, 2020.
  * [2] M. Swan, _Blockchain: Blueprint for a new economy_. O’Reilly Media, Inc, 2015.
  * [3] A. Dorri, S. Kanhere, R. Jurdak, and P. Gauravaram, “Blockchain for IoT security and privacy: The case study of a smart home,” in _Proc. IEEE International Conference on Pervasive Computing and Communications Workshops (PerCom Workshops)_ , Mar. 2017, pp. 618–623.
  * [4] J. Zhou, Z. Cao, X. Dong, and A. V. Vasilakos, “Security and privacy for cloud-based IoT: Challenges,” _IEEE Communications Magazine_ , vol. 55, no. 1, pp. 26–33, Jan. 2017.
  * [5] T. Maksymyuk, S. Dumych, M. Brych, D. Satria, and M. Jo, “An IoT based monitoring framework for software defined 5G mobile networks,” in _Proc. of the 11th International Conference on Ubiquitous Information Management and Communication_ , 2017, pp. 105:1–105:4.
  * [6] P. Danzi, A. Kalør, C̆. Stefanović, and P. Popovski, “Delay and communication tradeoffs for blockchain systems with lightweight IoT clients,” _IEEE Internet of Things Journal_ , vol. 6, no. 2, pp. 2354–2365, Apr. 2019.
  * [7] M. S. Ali, K. Dolui, and F. Antonelli, “Iot data privacy via blockchains and ipfs,” in _Proceedings of the Seventh International Conference on the Internet of Things_. ACM, 2017, p. 14.
  * [8] H. Zhou, X. Ouyang, Z. Ren, J. Su, C. de Laat, and Z. Zhao, “A blockchain based witness model for trustworthy cloud service level agreement enforcement,” in _Proc. IEEE Conference on Computer Communications (INFOCOM)_ , Apr. 2019, pp. 1567–1575.
  * [9] W. Stewart, _Probability, Markov chains, queues, and simulation: the mathematical basis of performance modeling_. Princeton University Press, 2009.
  * [10] Y. Kawase and S. Kasahar, “Transaction-confirmation time for bitcoin: A queueing analytical approach to blockchain mechanism,” in _Proc. International Conference on Queueing Theory and Network Applications_ , 2017, pp. 75–88.
  * [11] M. Chaudhry and J. Templeton, “The queuing system M/GB/l and its ramifications,” _European Journal of Operational Research_ , vol. 6, no. 1, pp. 56–60, 1981.