EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-PH-EP-2015-123 LHCb-PAPER-2015-011 September 21, 2015

Measurement of the exclusive $\Upsilon$ production cross-section in $pp$
collisions at $\sqrt{s}=7\mathrm{\,Te\kern-2.07413ptV}$ and
8$\mathrm{\,Te\kern-2.07413ptV}$

The LHCb collaboration†††Authors are listed at the end of this paper.

A study is presented of central exclusive production of $\Upsilon(nS)$ states,
where the $\Upsilon(nS)$ resonances decay to the ${\mu^{+}}{\mu^{-}}$ final
state, using $pp$ collision data recorded by the LHCb experiment. The cross-
section is measured in the rapidity range $2<y(\Upsilon)<4.5$ where the muons
are reconstructed in the pseudorapidity range $2<\eta(\mu^{\pm})<4.5$. The
data sample corresponds to an integrated luminosity of 2.9$\mbox{\,fb}^{-1}$
and was collected at centre-of-mass energies of
$7\mathrm{\,Te\kern-1.00006ptV}$ and $8\mathrm{\,Te\kern-1.00006ptV}$. The
measured $\Upsilon(1S)$ and $\Upsilon(2S)$ production cross-sections are

$\displaystyle\sigma(pp\rightarrow p\Upsilon(1S)p)$ $\displaystyle=$
$\displaystyle 9.0\pm 2.1\pm 1.7\textrm{ pb and}$
$\displaystyle\sigma(pp\rightarrow p\Upsilon(2S)p)$ $\displaystyle=$
$\displaystyle 1.3\pm 0.8\pm 0.3\textrm{ pb},$

where the first uncertainties are statistical and the second are systematic.
The $\Upsilon(1S)$ cross-section is also measured as a function of rapidity
and is found to be in good agreement with Standard Model predictions. An upper
limit is set at 3.4 pb at the 95% confidence level for the exclusive
$\Upsilon(3S)$ production cross-section, including possible contamination from
$\chi_{b}(3P)\rightarrow\Upsilon(3S)\gamma$ decays.

Published in JHEP

© CERN on behalf of the LHCb collaboration, licence CC-BY-4.0.

##  1 Introduction

Central exclusive production (CEP) of $\Upsilon(nS)$ $(n=1,2,3)$ resonances in
$pp$ collisions is thought to occur by photoproduction through the exchange of
a photon and a pomeron (a colour-singlet system) between two protons, as
illustrated in Fig. 1. Since the protons do not dissociate, typically only a
small component of momentum transverse to the beam direction ($p_{\rm T}$) is
exchanged in the interaction. The photoproduction of $\Upsilon$ resonances at
LHCb can be computed using perturbative quantum chromodynamics (QCD), given
the high photon-proton centre-of-mass energy, $W$, and the cross-section
depends on the square of the gluon parton-density function, $g(x)$, where
Bjorken-$x$ is the fraction of the proton’s momentum carried by the gluon [1].
Measurements of the production cross-sections for the $\Upsilon(nS)$
resonances in the forward region covered by the LHCb detector are sensitive to
$g(x)$ in the region of small $x$ down to approximately $1.5\times 10^{-5}$,
where the knowledge of $g(x)$ is limited. Furthermore, predictions for the
$\Upsilon(nS)$ cross-sections at leading order (LO) and next-to-leading order
(NLO) in the strong-interaction coupling differ greatly for the values of $W$
probed in $\Upsilon(nS)$ resonance production, and there are significant
variations depending on the models used to describe the $\Upsilon$ wave
function and the $t$-channel exchange [1, 2, 3].

Quarkonia photoproduction has been studied in exclusive production at HERA [4,
5, 6, 7, 8, 9], the Tevatron [10] and the LHC [11, 12, 13]. At LHCb, exclusive
production is associated with the absence of significant detector activity
apart from that associated with the exclusive candidate. The background from
proton dissociation occurring outside the detector acceptance is characterised
as having a value of $\Upsilon$ candidate $p_{T}$ which is larger than that
for exclusive production.

Figure 1: Leading Feynman diagram for photoproduction of $\Upsilon(nS)$
states, where the photon-pomeron interaction is indicated by the shaded grey
circle.

In this article, the exclusive production cross-section of $\Upsilon(nS)$
resonances is measured in the ${\mu^{+}}{\mu^{-}}$ final state where both
muons lie in the pseudorapidity ($\eta$) range $2<\eta(\mu^{\pm})<4.5$ and the
$\Upsilon(nS)$ candidate is reconstructed in the rapidity ($y$) range
$2<y(\Upsilon(nS))<4.5$. The $pp$ data correspond to an integrated luminosity
of 0.9$\mbox{\,fb}^{-1}$ at a $pp$ centre-of-mass energy of
$\sqrt{s}=7\mathrm{\,Te\kern-1.00006ptV}$ and 2.0$\mbox{\,fb}^{-1}$ at
$\sqrt{s}=8\mathrm{\,Te\kern-1.00006ptV}$. Given the limited statistical
precision, the data sets are combined to measure the production cross-
sections. The LHCb detector and the simulated event samples are outlined in
Sect. 2. In Sect. 3 selection criteria are discussed, which exploit the
absence of detector activity other than that associated with the
$\Upsilon(nS)$ candidate. The signal efficiency and the various sources of
background are also described. In Sect. 4 two fits are described, which allow
the determination of the exclusive signal yield: by fitting the $\Upsilon(nS)$
invariant mass spectrum in order to separate $\Upsilon$ resonances from dimuon
continuum background; and by fitting the $\Upsilon(nS)$ candidate $p^{2}_{\rm
T}$ distribution to distinguish exclusively produced $\Upsilon$ resonances
from those originating in hard interactions. Systematic uncertainties are
summarised in Sect. 5, and the measurements of the cross-sections are
discussed in Sect. 6. Finally, for the $\Upsilon(1S)$ the differential cross-
section, as a function of $\Upsilon(1S)$ candidate rapidity, is presented.

##  2 Detector and simulation

The LHCb detector [14, 15] is a single-arm forward spectrometer, designed for
the study of particles containing $b$ or $c$ quarks. It is fully instrumented
in the pseudorapidity range $2<\eta<5$ and has tracking capability in the
backward direction, in the range $-3.5<\eta<-1.5$. The detector includes a
high-precision tracking system consisting of a silicon-strip vertex detector
(VELO) surrounding the $pp$ interaction region, a large-area silicon-strip
detector located upstream of a dipole magnet with a bending power of about
$4{\rm\,Tm}$, and three stations of silicon-strip detectors and straw drift
tubes placed downstream of the magnet. The tracking system provides a
measurement of momentum, $p$, of charged particles with a relative uncertainty
that varies from 0.5% at low momentum to 1.0% at
200${\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$. Photons, electrons and hadrons are
identified by a calorimeter system consisting of scintillating-pad (SPD) and
preshower detectors, an electromagnetic calorimeter and a hadronic
calorimeter. Muons are identified by a system composed of alternating layers
of iron and multiwire proportional chambers.

The trigger consists of a hardware stage, based on information from the
calorimeter and muon systems, followed by a software stage, which applies a
full event reconstruction. The hardware trigger requires events to contain at
least one muon with a $p_{\rm T}$ greater than
200${\mathrm{\,Me\kern-1.00006ptV\\!/}c}$. Low-multiplicity events are
selected by requiring that fewer than ten hits should be detected in the
scintillating pad detector, positioned just upstream of the electromagnetic
calorimeter. In the subsequent software trigger, both of the final-state muons
are required to have $\mbox{$p_{\rm
T}$}>400{\mathrm{\,Me\kern-1.00006ptV\\!/}c}$.

The exclusive production of $\Upsilon(nS)$ resonances is simulated using the
SuperChiC software package [16], which provides the four-momentum of a single,
transversely polarised $\Upsilon(nS)$ resonance in each event. The decay of
the $\Upsilon(nS)$ candidate is described by EvtGen [17], in which final-state
radiation is generated using Photos [18]. The interaction of the generated
particles with the detector, and its response, are implemented using the
Geant4 toolkit [19, *Agostinelli:2002hh] as described in Ref. [21]. Samples,
each containing one million events, are prepared for $\Upsilon(nS)$ resonances
decaying to ${\mu^{+}}{\mu^{-}}$. In the same way, nine background samples of
a similar size are prepared containing events where an exclusively produced
$\chi_{b0,1,2}(1P,2P,3P)$ meson decays to the $\Upsilon(nS)\gamma$ final
states with $\Upsilon\rightarrow{\mu^{+}}{\mu^{-}}$. Separate samples are
prepared for every $\chi_{b}(mP)\rightarrow\Upsilon(nS)\gamma$
$(m,n=1,2,3;n\leq m)$ decay.

##  3 Candidate selection

Selection criteria are applied offline to events that pass the trigger
requirements, to select well-reconstructed $\Upsilon(nS)$ candidates and to
ensure the absence of unrelated detector activity. The latter set of
requirements favours events containing a single $pp$ interaction per bunch
crossing.

The final-state tracks, which must be associated with hits in the muon
chambers, are required to lie in the pseudorapidity range
$2<\eta(\mu^{\pm})<4.5$ and to be of good quality. In extracting the
differential cross-section for the $\Upsilon(1S)$, the following intervals in
$\Upsilon(1S)$ rapidity are considered: $2<y<3$, $3<y<3.5$ and $3.5<y<4.5$.
Dimuon candidates are selected if the invariant mass falls in the range
between 9${\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ and
20${\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$, and the candidate $p^{2}_{\rm
T}$ is less than 2$\mathrm{\,Ge\kern-1.00006ptV}^{2}/c^{2}$. The latter
requirement favours photoproduction candidates, which have a
characteristically low-$p_{\rm T}$. Events are rejected if one or more tracks
are reconstructed in the backward direction. In the forward region exactly two
tracks, corresponding to the muon candidates, are required, and these must be
reconstructed both in the VELO and in the downstream tracking detectors.

The selection criteria affect not only the $\Upsilon(nS)$ candidate but also
the level of activity in the rest of the event, specifically through the
requirements that there should be exactly two forward tracks, no backward
tracks and fewer than ten SPD hits. The event is excluded if more than one
proton-proton interaction occurs, causing a larger number of additional SPD
hits or extra tracks to be reconstructed. The probability for an exclusive
$\Upsilon$ event not to feature additional activity from another $pp$
interaction in the same beam crossing is determined as the fraction of events
containing no activity, according to these criteria, in a randomly accepted,
hence unbiased, sample. After subtracting the contribution of the dimuon
candidate, an event may contain fewer than eight SPD hits and no reconstructed
tracks in the backward direction or tracks in the forward region. The fraction
of randomly triggered events passing these criteria, $f_{\rm SI}$, is found to
be $(23.63\pm 0.04)\%$ for the 7$\mathrm{\,Te\kern-1.00006ptV}$ data and
$(18.48\pm 0.02)\%$ for the 8$\mathrm{\,Te\kern-1.00006ptV}$ data. The
difference arises because of the different beam conditions in the two data-
taking periods, leading to a different average number of proton-proton
interactions per event.

The reconstruction, trigger and offline selection efficiencies are determined
using simulated samples, and the combined efficiency varies between 77% and
84%. For signal candidates that pass the trigger and reconstruction stages,
the offline selection criteria are more than 99% efficient.

##  4 Determining the exclusive yield

Candidates reconstructed in the 7$\mathrm{\,Te\kern-1.00006ptV}$ and
8$\mathrm{\,Te\kern-1.00006ptV}$ data sets are combined in a single sample,
and two unbinned, extended, maximum-likelihood fits are carried out. A first
fit is performed to the dimuon invariant mass spectrum, between
9${\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ and
20${\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$. The fit contains a non-resonant
background component and three resonant components. The three resonant
components each receive contributions from exclusive signal, inelastic
background and $\chi_{b}\rightarrow\Upsilon\gamma$ feed-down decays. These
contributions are indistinguishable in the invariant mass distribution.

The probability density function (PDF) used to model each $\Upsilon(nS)$
signal peak is a Gaussian function with modified tails (a double-sided crystal
ball function [22]). The mass differences for the
$\Upsilon(2S){\rm-}\Upsilon(1S)$ and $\Upsilon(3S){\rm-}\Upsilon(1S)$
resonances are taken from Ref. [23]. The ratios of the $\Upsilon(2S)$ and
$\Upsilon(3S)$ resolutions with respect to the $\Upsilon(1S)$ are fixed to the
ratio of their masses with respect to the mass of the $\Upsilon(1S)$,
following the procedure used in previous $\Upsilon$ measurements using LHCb
data [24]. The parameters that govern the shapes of the tails are taken from
simulation, as is the resolution of the $\Upsilon(1S)$ resonance, which varies
from 35${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ to
57${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ in the different rapidity ranges.
The yields of the signal components are all free to vary independently.

A background PDF accounting for the non-resonant background is modelled using
an exponential shape where the slope and normalisation are allowed to vary.

The data are fitted in the whole rapidity range and in bins of rapidity. The
fit results are given in Table 1 and the fit in the full rapidity range is
shown between $9{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ and
$12{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$ in Fig. 2.

Figure 2: Invariant dimuon mass spectrum for 7 TeV and 8 TeV data in the rapidity range $2<y(\Upsilon)<4.5$ (black points). The fit PDF is superimposed (solid blue line). The $\Upsilon(1S,2S,3S)$ signal components, used to derive weights, are indicated with a long-dashed (red) line, and the non-resonant background is marked with a short-dashed (grey) line.  Table 1: Results of the invariant mass fits, within each rapidity interval.  Parameter | $2<y<4.5$ | $2<y<3$ | $3<y<3.5$ | $3.5<y<4.5$
---|---|---|---|---
$\Upsilon(1S,2S,3S)$ yield | 382 $\pm$ | 26 | 146 $\pm$ | 16 | 133 $\pm$ | 16 | 94 $\pm$ | 14
$\Upsilon(1S)$ fraction | 0.71 $\pm$ | 0.03 | 0.74 $\pm$ | 0.05 | 0.72 $\pm$ | 0.06 | 0.68 $\pm$ | 0.07
$\Upsilon(2S)$ fraction | 0.18 $\pm$ | 0.03 | 0.16 $\pm$ | 0.04 | 0.15 $\pm$ | 0.05 | 0.26 $\pm$ | 0.06
$\Upsilon(1S)$ mass (${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$) | 9452.5 $\pm$ | 3.3 | 9453.2 $\pm$ | 4.3 | 9452.4 $\pm$ | 5.6 | 9452.0 $\pm$ | 9.0

Two sources of background contribute to the fitted signal: feed-down from
$\chi_{b}\rightarrow\Upsilon\gamma$ decays, and inelastic interactions that
involve the undetected products of proton dissociation or additional gluon
radiation.

The feed-down background is estimated using a combination of data and
simulation, considering $\chi_{b}(mP)\rightarrow\Upsilon(nS)\gamma$ decays.
Events are considered in the data set if exactly one photon is found in
addition to the $\Upsilon$ candidate. Regions in the $\Upsilon\gamma$
invariant mass spectrum are defined, corresponding to the $\chi_{b}(1P,2P,3P)$
states, and the number of $\chi_{b}$ candidates, $N_{\chi_{b}}$, for each
decay $\chi_{b}(mP)\rightarrow\Upsilon(nS)\gamma$ is counted. An estimate of
the total feed-down content of the $\Upsilon$ data sample from each $\chi_{b}$
state is found using the expression:

$N_{\textrm{feed-down,
}\chi_{b}(mP)\rightarrow\Upsilon(nS)\gamma}=\frac{N_{\chi_{b}}\times\mathcal{F}}{\epsilon_{\gamma}\times\epsilon_{\textrm{mass-
range}}}.$ (1)

Here $\mathcal{F}$ is the purity of the $\Upsilon(nS)$ in the corresponding
mass window with respect to the non-resonant ${\mu^{+}}{\mu^{-}}\gamma$
background, determined by fitting the dimuon mass spectrum for events with
exactly one reconstructed photon; $\epsilon_{\gamma}$ is the efficiency for
reconstructing the photon produced in each $\chi_{b}(mP)$ decay, determined
using simulated exclusive $\chi_{b}(mP)\rightarrow\Upsilon(nS)\gamma$ decays;
and $\epsilon_{\textrm{mass-range}}=0.9$ corrects for the fraction of signal
$\Upsilon$ candidates which are expected to fall outside the mass window.
There are too few $\Upsilon(3S)\gamma$ candidates to estimate the purity
precisely so it is assumed to be 100%. Because of limited mass resolution and
small sample sizes the $\chi_{b}$ spin states cannot be resolved, so equal
contributions from the $\chi_{b1}(mP)$ and $\chi_{b2}(mP)$ states are assumed.
The $\chi_{b0}$ radiative decay rate is expected to be relatively suppressed
and is therefore neglected [23]. The feed-down background yields are given in
Table 2.

Table 2: Estimated yields of feed-down background from $\chi_{b}(mP)\rightarrow\Upsilon(nS)\gamma$ decays in each $\Upsilon(nS)$ sample, where the uncertainties are statistical only.  Signal window |  $\Upsilon$ sample | Estimated contamination yield
---|---|---
|  | $\chi_{b}(1P)$ | $\chi_{b}(2P)$ | $\chi_{b}(3P)$
$2<y(\Upsilon)<4.5$ | $\Upsilon(1S)$ | 63 $\pm$ | 10 | 14 $\pm$ | 5 | 3 $\pm$ | 2
| $\Upsilon(2S)$ | $-$ | 43 $\pm$ | 12 | 5 $\pm$ | 3
| $\Upsilon(3S)$ | $-$ | $-$ | 21 $\pm$ | 21
$2<y(\Upsilon)<3$ | $\Upsilon(1S)$ | 31 $\pm$ | 8 | 2 $\pm$ | 2 | 0 $\pm$ | 2
$3<y(\Upsilon)<3.5$ | $\Upsilon(1S)$ | 22 $\pm$ | 6 | 10 $\pm$ | 4 | 0 $\pm$ | 2
$3.5<y(\Upsilon)<4.5$ | $\Upsilon(1S)$ | 8 $\pm$ | 4 | 0 $\pm$ | 2 | 3 $\pm$ | 2

Since the mass shapes for signal and background do not significantly depend on
$p_{\rm T}$ over the $p_{\rm T}$ range considered, the $p^{2}_{\rm T}$
distribution of the $\Upsilon$ candidates is determined using the sPlot
technique [25]. A fit is then performed to the $p^{2}_{\rm T}$ distribution,
shown in Fig. 3, using candidates in the full rapidity range
$2.0<y(\Upsilon)<4.5$, with fit components corresponding to the $\Upsilon$
signal, inelastic background and feed-down background. The fraction of
exclusive signal calculated from this fit is assumed to be the same for each
rapidity bin.

The $p^{2}_{\rm T}$ distribution for the exclusive signal is derived from the
simulated sample. At HERA, the distribution of the exclusive charmonium signal
as a function of the candidate $p^{2}_{\rm T}$ was well described by an
exponential function, $\exp(-b\mbox{$p^{2}_{\rm T}$})$ [9, 6]. The $p^{2}_{\rm
T}$ distribution provides discrimination among various production sources
because it approximates the squared four-momentum transfer, $|t|$, which
depends on the production mechanism. Following Ref. [12], Regge phenomenology
is used to extrapolate the slope measured by HERA up to LHC energies according
to the expression

$b(W)=b_{0}+4\alpha^{\prime}\log\left(\frac{W}{W_{0}}\right),$ (2)

where $\alpha^{\prime}$ describes the slope of the exchange Regge trajectory
and the constant $b_{0}$ is specific for interactions at a given photon-proton
centre-of-mass energy, $W_{0}$ [1]. The SuperChiC generator [16] models the
pomeron-photon exchange and performs this extrapolation. Since the only
published measurement of $b_{\Upsilon}$ has very low precision [8], the
generator is tuned to reproduce the LHCb measurement of $b_{J/\psi}=5.7\pm
0.1$ $\mathrm{\,Ge\kern-1.00006ptV}^{-2}c^{2}$ [12] in exclusive $J/\psi$
production. The input values to SuperChiC are $b_{0}=5.6$
$\mathrm{\,Ge\kern-1.00006ptV}^{-2}c^{2}$, $\alpha^{\prime}=0.2$
$\mathrm{\,Ge\kern-1.00006ptV}^{-2}c^{2}$ and
$W_{0}=90$$\mathrm{\,Ge\kern-1.00006ptV}$. The fit PDF is obtained from the
simulated samples with these inputs, and a kernel estimation is employed to
derive a shape to fit to data [26]. It is assumed that the inelastic
background component is distributed according to a single exponential function
[12]. The slope and yield of this background function are free to vary in the
fit.

The total contamination from $\chi_{b}(mP)\rightarrow\Upsilon(nS)\gamma$
decays is constrained to be the sum of the contributions in Table 2 and enters
the $p^{2}_{\rm T}$ fit by means of a Gaussian constraint. Given that no
analysis of exclusive $\chi_{b}$ production has been undertaken, and the
consequent lack of knowledge of the exclusive purity of the very small sample
of reconstructed $\chi_{b}$ candidates, it is assumed that inelastic processes
contribute half of this feed-down background and the same inelastic background
PDF is employed as that used to model inelastic $\Upsilon(nS)$ production. For
the exclusive component, the shape of the dimuon $p^{2}_{\rm T}$ PDF depends
on the $\chi_{b}(mP)$ meson source, and the PDFs for each source, determined
from simulation, are combined according to their relative contributions to the
total feed-down yield.

The exclusive purity, ${\cal P}$, is defined as the ratio of the exclusive
signal yield, $N_{\rm exclusive}$, to the number of candidates remaining in
the sample, $N_{\rm exclusive}+N_{\rm inelastic}$, after subtraction of the
feed-down yield. The fit to the $p^{2}_{\rm T}$ distribution in the full
rapidity range, shown in Fig. 3, gives

${\cal P}\equiv\frac{N_{\rm exclusive}}{N_{\rm exclusive}+N_{\rm
inelastic}}=(54\pm 11)\,\%,$

with the exponential slope of the inelastic background measured to be
$-0.21\pm 0.26\mathrm{\,Ge\kern-1.00006ptV}^{-2}c^{2}$. The results of fits to
the $p^{2}_{\rm T}$ distribution in each rapidity interval are consistent with
this value. In order to validate the fit procedure, a set of pseudoexperiments
is generated using the parameters obtained from the fit to the data, and the
same fit is applied to each pseudoexperiment. The uncertainty on the purity is
underestimated by 15% in the fit and the statistical uncertainty quoted takes
account of this.

Figure 3: Fit to the $p^{2}_{\rm T}$ distribution of the $\Upsilon$ candidates
in the full rapidity range.

##  5 Systematic uncertainties

The relative systematic uncertainties for the $\Upsilon(1S,2S,3S)$ cross-
sections in the various rapidity ranges are summarised in Table 3.

Contributions to the systematic uncertainty arising from the $p^{2}_{\rm T}$
fit are considered: the uncertainty in the signal $p^{2}_{\rm T}$ distribution
as modelled by the SuperChiC generator and the variation of the exclusive
signal PDF expected in the various rapidity bins. The SuperChiC generator is
tuned to reproduce measurements of exclusive $J/\psi$ meson production made by
LHCb [12]. As no sufficiently precise measurements of the $p^{2}_{\rm T}$
distribution in exclusive $\Upsilon(nS)$ resonance production exist, an
estimate is made following Ref. [1], where it is argued from Regge theory that
the slope $b_{0}$ of the proton should be reduced by
$4\alpha^{\prime}\log(m_{\Upsilon(nS)}/m_{J/\psi})$. A simulated sample is
generated accordingly and used to derive a signal $p^{2}_{\rm T}$ template.
Changing $b_{0}$ from 5.6 to 4.7 produces a relative decrease in the exclusive
yields of 6%, and this change is taken as the systematic uncertainty. For the
differential cross-section measurements, the dependence of the signal
$p^{2}_{\rm T}$ shape on rapidity is studied by replacing the exclusive signal
$p^{2}_{\rm T}$ PDF with those determined in the smaller rapidity ranges, and
the largest change in purity is taken as the uncertainty. Combining the
systematic uncertainties in quadrature yields a total uncertainty for the
exclusive purity, ${\cal P}$, between 7.2% and 8.2%. In addition, the
possibility for variation in the shape of the continuum dimuon background in
$p^{2}_{\rm T}$ as a function of mass is considered. The determination of the
exclusive purity, $\mathcal{P}$, is repeated in the dimuon invariant mass
range from 9 to $12{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$, and the
difference is taken as a conservative estimate of the systematic uncertainty.
In Table 3 these sources contribute to the uncertainty denoted ‘purity fit’.

The uncertainty arising from the $p^{2}_{\rm T}$ shape derived from simulation
and used to describe the feed-down background is considered separately. The
feed-down background PDF is constructed using only contributions from the
$\chi_{b1}(mP)$ and $\chi_{b2}(mP)$ background components, in equal parts,
assuming no contribution from $\chi_{b0}(mP)$ decays. Since it is not possible
to resolve the spin states in data, we consider a conservative change where
the nominal PDF is replaced with that for background originating from the
decay of a $\chi_{b0}(mP)$ meson. The fit to data is repeated and the change
in exclusive purity is taken as the associated uncertainty. In addition, there
is uncertainty associated with the exclusive fraction of the $\chi_{b}(mP)$
feed-down background, which in turn affects the overall shape of the feed-down
background used in the $p^{2}_{\rm T}$ spectrum. Since the size of the data
set is too small to allow a data-driven estimate of the $\chi_{b}(mP)$ sample
exclusive purity, the PDF for this purity is assumed to be uniform between 0%
and 100%, and the effect of changing it by $\pm 1$ standard deviation
($0.50\pm 0.29$) is therefore considered. The resulting change in the
exclusive $\Upsilon(nS)$ sample purity is taken as the systematic uncertainty.
The estimate of the yield of $\chi_{b}(mP)$ meson feed-down in data includes
the determination of the photon reconstruction efficiency using simulated
samples of $\chi_{b}(mP)\rightarrow\Upsilon(nS)\gamma$ decays. Samples of
$B^{\pm}\rightarrow J\psi K^{*\pm}(\rightarrow{{K}^{\pm}}{{\pi}^{0}})$ and
${{B}^{\pm}}\rightarrow J/\psi{{K}^{\pm}}$ decays are used to validate the
agreement between photon reconstruction efficiencies in data and simulation.
An uncertainty of 5% is taken as the systematic uncertainty on the photon
reconstruction efficiency to account for the small differences seen [27]. The
resulting uncertainty is very small for the $\Upsilon(1S)$ but larger for the
$\Upsilon(2S,3S)$ samples where the relative contamination from the
$\chi_{b}(mP)$ background is larger. These three systematic uncertainties on
the cross-section are combined in quadrature and are presented as the ‘feed-
down b.g.’ systematic uncertainty in Table 3.

An estimate of the contamination of the $\Upsilon(1S)$ and $\Upsilon(2S)$
samples from the decays
$\Upsilon(2S)\rightarrow\Upsilon(1S){{\pi}^{0}}{{\pi}^{0}}$,
$\Upsilon(3S)\rightarrow\Upsilon(2S)\\{{{\pi}^{0}}{{\pi}^{0}},\gamma\gamma\\}$,
is made using the observed $\Upsilon(nS)$ candidate yields in data and the
relevant $\Upsilon^{\prime}\rightarrow\Upsilon X$ and
$\Upsilon^{(^{\prime})}\rightarrow{\mu^{+}}{\mu^{-}}$ branching fractions
[23]. The estimated contaminations are taken as systematic uncertainties for
the $\Upsilon(1S)$ and $\Upsilon(2S)$ cross-sections.

Uncertainties in the PDFs used to fit the $\Upsilon(nS)$ candidate invariant
mass spectrum are considered. Alternative $\Upsilon(nS)$ signal PDFs are
produced, obtained using kernel estimation. The systematic uncertainties on
the exclusive purity and each of the yields are assessed using
pseudoexperiments generated with the nominal invariant mass PDFs and fitted
with a model where the signal PDFs are replaced by those obtained using kernel
estimation. The effect of replacing the exponential PDF used to model the non-
resonant background in the invariant mass fit with a second-order polynomial
function is found in the same way using pseudoexperiments. Combining these two
sources of uncertainty in quadrature leads to a relative uncertainty which is
less than 4% for all the cross-sections. This uncertainty is labelled ‘mass
fit’ in Table 3.

The LHCb integrated luminosity has been measured with a relative uncertainty
of 1.7% at 7$\mathrm{\,Te\kern-1.00006ptV}$ and 1.2% at
8$\mathrm{\,Te\kern-1.00006ptV}$ [28]. The integrated luminosity is multiplied
by the estimated, selection-dependent, fraction of events that contain no
interactions other than the one that produces the signal candidate, $f_{\rm
SI}$. The determination of $f_{\rm SI}$ from data depends upon the subtraction
of the $\Upsilon$ signal candidate’s SPD hits. The spread in the signal
candidate SPD hit multiplicity is estimated from data to be one hit, and the
fraction $f_{\rm SI}$ is therefore recomputed with the signal subtraction
increased from two to three and the change is taken to be the systematic
uncertainty. To account for variations as a function of data-taking time, the
variation of the estimated single-interaction fraction is evaluated in each
uninterrupted period of data-taking during which conditions are constant,
typically an hour long, instead of considering each year as a whole, and the
change with respect to the nominal fraction is taken as the systematic
uncertainty. Combining the uncertainties in quadrature yields an overall
relative uncertainty for each year of 2.3%. The systematic uncertainties on
the luminosity for each year are assumed to be 100% correlated.

The branching fractions, $\mathcal{B}$, for $\Upsilon(1S)$, $\Upsilon(2S)$ and
$\Upsilon(3S)$ to decay to the dimuon final state are accounted for to
determine the $\Upsilon(nS)$ production cross-section. These branching
fractions are taken from Ref. [23] and, for the $\Upsilon(1S)$, $\Upsilon(2S)$
and $\Upsilon(3S)$ states, carry relative uncertainties of 2%, 9% and 10%,
respectively. These are propagated to an uncertainty on the production cross-
section.

Table 3: Summary of the relative systematic uncertainties, in %. | $2<y<3$ | $3<y<3.5$ | $3.5<y<4.5$ | $2<y<4.5$
---|---|---|---|---
| $\Upsilon(1S)$ | $\Upsilon(1S)$ | $\Upsilon(1S)$ | $\Upsilon(1S)$ | $\Upsilon(2S)$ | $\Upsilon(3S)$
Purity fit | 14.2 | 14.2 | 14.2 | 13.7 | 13.7 | 13.7
Feed-down b.g. | 12.2 | 12.2 | 12.3 | 12.2 | 14.6 | 12.5
$\Upsilon^{\prime}$ feed-down |  04.0 |  04.3 |  05.4 |  04.5 | 11.1 |  0$-$
Mass fit |  02.2 |  02.8 |  02.9 |  02.1 |  02.8 |  03.6
Luminosity |  02.3 |  02.3 |  02.3 |  02.3 |  02.3 |  02.3
$\mathcal{B}(\Upsilon\rightarrow{\mu^{+}}{\mu^{-}})$ |  02.0 |  02.0 |  02.0 |  02.0 |  08.8 |  09.6
Total | 19.5 | 19.7 | 20.0 | 19.3 | 24.8 | 21.4

##  6 Cross-section

The cross-section is obtained using

$\sigma=\frac{N_{\rm exclusive}}{{\cal
L}\times\epsilon\times\mathcal{B}(\Upsilon(nS)\rightarrow\mu^{+}\mu^{-})}.$
(3)

The effective integrated luminosity, $\cal{L}$, is 580$\mbox{\,pb}^{-1}$,
taking into account the values of $f_{\textrm{SI}}$ given in Sect. 3 [28].

The quantity $\epsilon$ is the efficiency correction, which is obtained in
each rapidity bin and for each resonance, and which is averaged for
7$\mathrm{\,Te\kern-1.00006ptV}$ and 8$\mathrm{\,Te\kern-1.00006ptV}$ data
according to the luminosity in each year, and
$\mathcal{B}(\Upsilon(nS)\rightarrow\mu^{+}\mu^{-})$ is the
$\Upsilon(nS)\rightarrow{\mu^{+}}{\mu^{-}}$ branching fraction. The measured
exclusive production cross-sections in the LHCb acceptance are

$\displaystyle\sigma(pp\rightarrow p\Upsilon(1S)p)$ $\displaystyle=$
$\displaystyle 9.0\pm 2.1\pm 1.7\textrm{ pb},$
$\displaystyle\sigma(pp\rightarrow p\Upsilon(2S)p)$ $\displaystyle=$
$\displaystyle 1.3\pm 0.8\pm 0.3\textrm{ pb},\textrm{ and}$
$\displaystyle\sigma(pp\rightarrow p\Upsilon(3S)p)$ $\displaystyle<$
$\displaystyle 3.4\textrm{ pb}\textrm{ at the 95\% confidence level},$

where the first uncertainties are statistical and the second are systematic,
and where the limit on $\sigma(\Upsilon(3S))$ includes possible contamination
from $\chi_{b}$ feed-down. The limit is calculated using pseudo-experiments
and includes the effect of systematic uncertainties, where correlations are
assumed to be negligible. The $\Upsilon(1S)$ production cross-section is given
in smaller ranges of $\Upsilon(1S)$ rapidity in Table 4.

Table 4: Production cross-section for the $\Upsilon(1S)$ resonance in ranges of $\Upsilon(1S)$ rapidity, where the muons are required to lie in the pseudorapidity range $2<\eta(\mu^{\pm})<4.5$. The first uncertainties are statistical and the second are systematic. | $2<y<3$ | $3<y<3.5$ | $3.5<y<4.5$
---|---|---|---
$\sigma(\Upsilon(1S))$ (pb) | 3.4 $\pm$ 0.9 $\pm$ 0.7 | 2.9 $\pm$ 0.8 $\pm$ 0.6 | 2.6 $\pm$ 0.8 $\pm$ 0.5

After correction for the LHCb geometrical acceptance, the cross-sections in
Table 4 can be compared to theoretical predictions. The efficiency for an
$\Upsilon$ candidate to be produced in the range $2<y(\Upsilon)<4.5$, and to
decay to muons which lie inside the acceptance, $2<\eta(\mu^{\pm})<4.5$, is
45%. In the smaller ranges of $\Upsilon(1S)$ rapidity considered for the
differential cross-section measurement, the efficiency is lowest in the outer
ranges, at 39% ($2<y<3$) and 36% ($3.5<y<4.5$), and highest in the central
range, at around 74% ($3<y<3.5$). The correction depends on the rapidity
distribution in the simulated sample, which has a different shape to those of,
for example, the predictions in Fig. 4a. To estimate the systematic
uncertainty on the geometrical acceptance correction, the simulated samples
are reweighted to obtain a uniform rapidity distribution within each rapidity
bin, and the change in the geometrical acceptance is taken as the systematic
uncertainty. This corresponds to a relative change in the geometrical
acceptance of less than 6%. The differential cross-sections,
$d\sigma(\Upsilon(1S))/dy$ are given in Table 5 and are shown in Fig. 4a
compared to LO and NLO predictions [1]. The LHCb data are in good agreement
with the NLO prediction.

Table 5: Measured $d\sigma(\Upsilon(1S))/dy$, where the data have been corrected for the effect of the LHCb geometrical acceptance. The statistical and systematic uncertainties are combined in quadrature. | $2<y<3$ | $3<y<3.5$ | $3.5<y<4.5$
---|---|---|---
$d\sigma(\Upsilon(1S))/dy$ (pb) | 8.8 $\pm$ 3.0 | 7.8 $\pm$ 2.7 | 7.1 $\pm$ 2.6

Figure 4: Measurements of exclusive $\Upsilon(1S)$ photoproduction compared to
theoretical predictions. In (a), the $\Upsilon(1S)$ cross-section in bins of
rapidity is shown, compared to LO and NLO predictions. The LHCb measurements
are indicated by black points with error bars for uncorrelated errors, and
solid rectangles indicating the total uncertainty. In (b), the photon-proton
cross-sections extracted from the LHCb results are indicated by black points,
where the statistical and systematic uncertainties are combined in quadrature.
The entire $W$-region in which these LHCb measurements are sensitive is
indicated. Measurements made by H1 and ZEUS in the low-$W$ region are
indicated by red and blue markers, respectively [4, 5, 7]. Predictions from
Ref. [1] are included, resulting from LO and NLO fits to exclusive $J/\psi$
production data. The filled bands indicate the theoretical uncertainties on
the 7 TeV prediction and the solid lines indicate the central values of the
predictions for 8 TeV. In (b) predictions from Ref. [2] using different models
for the $\Upsilon(1S)$ wave function are included, indicated by ‘bCGC’.

The LHCb results are also compared to theoretical predictions for the
underlying photon-proton cross-section, as a function of the centre-of-mass
energy of the photon-proton system, $W$, as shown in Fig. 4b. There are two
contributions to the photoproduction of an $\Upsilon(nS)$ resonance, depending
on which proton emits the virtual photon. The $pp$ cross-section is given by

$\frac{d\sigma^{\textrm{th}}(pp\rightarrow
p\Upsilon(1S)p)}{dy}=S^{2}(W_{+})\left(k_{+}\frac{dn}{dk_{+}}\right)\sigma_{+}^{\textrm{th}}(\gamma
p)+S^{2}(W_{-})\left(k_{-}\frac{dn}{dk_{-}}\right)\sigma_{-}^{\textrm{th}}(\gamma
p),$ (4)

where the predictions for the photon-proton cross-section are weighted by
absorptive corrections $S^{2}(W_{\pm})$ and the photon fluxes
$\frac{dn}{dk_{\pm}}$ for photons of energy
$k_{\pm}\approx(M_{\Upsilon(nS)}/2)\exp(\pm|y|)$. The absorptive corrections
and photon fluxes are computed following Ref. [1].

The three bins of $\Upsilon(1S)$ rapidity chosen in this analysis correspond
to ranges of $W$ for the $W_{+}$ and $W_{-}$ solutions. The contribution to
the total cross-section from the $W_{-}$ solutions is expected to be small and
is therefore neglected. The dominant $W_{+}$ solutions are therefore estimated
assuming that they dominate the cross-section, and are shown in Fig. 4b. The
magnitude of the theoretical prediction for the $W_{-}$ solutions is added as
a systematic uncertainty. The good agreement with the NLO prediction seen in
Fig. 4a is reproduced. The LHCb measurements probe a new kinematic region
complementary to that studied at HERA [4, 5, 7], as seen in Fig. 4b, and
discriminate between LO and NLO predictions. In Fig. 4b, the LHCb data are
also compared to the predictions given in Ref. [2] using models conforming to
the colour glass condensate (CGC) formalism [29] that take into account the
$t$-dependence of the differential cross-section. All agree well with the
data. The solid (black) and dotted (blue) lines correspond to two different
models for the scalar part of the vector-meson wave function.

##  7 Conclusion

The first measurement of exclusive $\Upsilon(nS)$ production in $pp$
collisions at 7 and 8$\mathrm{\,Te\kern-1.00006ptV}$ is presented, and a
differential cross-section is extracted as a function of $\Upsilon(1S)$
candidate rapidity. The data probe a previously unexplored kinematic region in
photon-proton centre-of-mass energy. The results are compared to theoretical
predictions and a strong preference for those including next-to-leading order
calculations is seen. Exclusive production studies at LHCb will be improved
during LHC Run II following the installation of scintillators at high
$|\eta|$, which will allow for improved trigger efficiency for exclusive
production processes and additional suppression of the background from
inelastic interactions [30].

## Acknowledgements

The authors gratefully acknowledge discussions with Lucian Harland-Lang
concerning the implementation and configuration of SuperChiC to model the
exclusive production process. We express our gratitude to our colleagues in
the CERN accelerator departments for the excellent performance of the LHC. We
thank the technical and administrative staff at the LHCb institutes. We
acknowledge support from CERN and from the national agencies: CAPES, CNPq,
FAPERJ and FINEP (Brazil); NSFC (China); CNRS/IN2P3 (France); BMBF, DFG, HGF
and MPG (Germany); INFN (Italy); FOM and NWO (The Netherlands); MNiSW and NCN
(Poland); MEN/IFA (Romania); MinES and FANO (Russia); MinECo (Spain); SNSF and
SER (Switzerland); NASU (Ukraine); STFC (United Kingdom); NSF (USA). The Tier1
computing centres are supported by IN2P3 (France), KIT and BMBF (Germany),
INFN (Italy), NWO and SURF (The Netherlands), PIC (Spain), GridPP (United
Kingdom). We are indebted to the communities behind the multiple open source
software packages on which we depend. We are also thankful for the computing
resources and the access to software R&D tools provided by Yandex LLC
(Russia). Individual groups or members have received support from EPLANET,
Marie Skłodowska-Curie Actions and ERC (European Union), Conseil général de
Haute-Savoie, Labex ENIGMASS and OCEVU, Région Auvergne (France), RFBR
(Russia), XuntaGal and GENCAT (Spain), Royal Society and Royal Commission for
the Exhibition of 1851 (United Kingdom).

## References

  * [1] S. P. Jones, A. D. Martin, M. G. Ryskin, and T. Teubner, _Probes of the small $x$ gluon via exclusive $J/\psi$ and $\Upsilon$ production at HERA and the LHC_, JHEP 1311 (2013) 085, arXiv:1307.7099
  * [2] V. P. Gonçalves, B. D. Moreira, and F. S. Navarra, _Exclusive $\Upsilon$ photoproduction in hadronic collisions at CERN LHC energies_, Phys. Lett. B742 (2015) 172, arXiv:1408.1344
  * [3] H. G. Dosch and E. Ferreira, _Scale dependent pomeron intercept in electromagnetic diffractive processes_ , arXiv:1503.06649
  * [4] ZEUS collaboration, J. Breitweg et al., _Measurement of elastic $\Upsilon$ photoproduction at HERA_, Phys. Lett. B437 (1998) 432, arXiv:hep-ex/9807020
  * [5] H1 collaboration, C. Adloff et al., _Elastic photoproduction of $J/\psi$ and $\Upsilon$ mesons at HERA_, Phys. Lett. B483 (2000) 23, arXiv:hep-ex/0003020
  * [6] ZEUS collaboration, S. Chekanov et al., _Exclusive photoproduction of $J/\psi$ mesons at HERA_, Eur. Phys. J. C24 (2002) 345, arXiv:hep-ex/0201043
  * [7] ZEUS collaboration, S. Chekanov et al., _Exclusive photoproduction of $\Upsilon$ mesons at HERA_, Phys. Lett. B680 (2009) 4, arXiv:0903.4205
  * [8] ZEUS collaboration, H. Abramowicz et al., _Measurement of the $t$ dependence in exclusive photoproduction of $\Upsilon(1S)$ mesons at HERA_, Phys. Lett. B708 (2012) 14, arXiv:1111.2133
  * [9] H1 collaboration, C. Alexa et al., _Elastic and proton-dissociative photoproduction of $J/\psi$ mesons at HERA_, Eur. Phys. J. C73 (2013) 2466, arXiv:1304.5162
  * [10] CDF collaboration, T. Aaltonen et al., _Observation of exclusive charmonium production and $\gamma+\gamma\rightarrow\mu^{+}\mu^{-}$ in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV_, Phys. Rev. Lett. 102 (2009) 242001, arXiv:0902.1271
  * [11] ALICE, B. Abelev et al., _Coherent $J/\psi$ photoproduction in ultra-peripheral Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV_, Phys. Lett. B718 (2013) 1273, arXiv:1209.3715
  * [12] LHCb collaboration, R. Aaij et al., _Updated measurements of exclusive $J/\psi$ and $\psi$(2S) production cross-sections in pp collisions at $\sqrt{s}=7$ TeV_, J. Phys. G41 (2014) 055002, arXiv:1401.3288
  * [13] ALICE, B. B. Abelev et al., _Exclusive $\mathrm{J/}\psi$ photoproduction off protons in ultra-peripheral p-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV_, Phys. Rev. Lett. 113 (2014), no. 23 232504, arXiv:1406.7819
  * [14] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [15] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [16] L. A. Harland-Lang, V. A. Khoze, M. G. Ryskin, and W. J. Stirling, _Standard candle central exclusive processes at the Tevatron and LHC_ , Eur. Phys. J. C69 (2010) 179, arXiv:1005.0695
  * [17] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [18] P. Golonka and Z. Was, _PHOTOS Monte Carlo: A precision tool for QED corrections in $Z$ and $W$ decays_, Eur. Phys. J. C45 (2006) 97, arXiv:hep-ph/0506026
  * [19] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [20] Geant4 collaboration, S. Agostinelli et al., _Geant4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [21] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [22] T. Skwarnicki, A study of the radiative cascade transitions between the $\Upsilon^{\prime}$ and $\Upsilon$ resonances, PhD thesis, Institute of Nuclear Physics, Krakow, 1986, DESY-F31-86-02
  * [23] Particle Data Group, K. A. Olive et al., _Review of particle physics_ , Chin. Phys. C38 (2014) 090001
  * [24] LHCb collaboration, R. Aaij et al., _Measurement of $\Upsilon$ production in pp collisions at $\sqrt{s}$ = 7 TeV_, Eur. Phys. J. C72 (2012) 2025, arXiv:1202.6579
  * [25] M. Pivk and F. R. Le Diberder, _sPlot: A statistical tool to unfold data distributions_ , Nucl. Instrum. Meth. A555 (2005) 356, arXiv:physics/0402083
  * [26] K. S. Cranmer, _Kernel estimation in high-energy physics_ , Comput. Phys. Commun. 136 (2001) 198, arXiv:hep-ex/0011057
  * [27] E. Govorkova, _Study of photon reconstruction efficiency using $B^{+}\rightarrow J/\psi K^{(*)+}$ decays in the LHCb experiment_, arXiv:1505.02960
  * [28] LHCb collaboration, R. Aaij et al., _Precision luminosity measurements at LHCb_ , JINST 9 (2014) P12005, arXiv:1410.0149
  * [29] J. Jalilian-Marian, A. Kovner, L. D. McLerran, and H. Weigert, _The intrinsic glue distribution at very small x_ , Phys. Rev. D55 (1997) 5414, arXiv:hep-ph/9606337
  * [30] FSC Team (CMS), HERSCHEL Team (LHCb), M. Albrow, P. Collins, and A. Penzo, _Forward shower counters for diffractive physics at the LHC_ , Int. J. Mod. Phys. A29 (2014) 1446018

LHCb collaboration

R. Aaij38, B. Adeva37, M. Adinolfi46, A. Affolder52, Z. Ajaltouni5, S. Akar6,
J. Albrecht9, F. Alessio38, M. Alexander51, S. Ali41, G. Alkhazov30, P.
Alvarez Cartelle53, A.A. Alves Jr57, S. Amato2, S. Amerio22, Y. Amhis7, L.
An3, L. Anderlini17,g, J. Anderson40, M. Andreotti16,f, J.E. Andrews58, R.B.
Appleby54, O. Aquines Gutierrez10, F. Archilli38, P. d’Argent11, A.
Artamonov35, M. Artuso59, E. Aslanides6, G. Auriemma25,n, M. Baalouch5, S.
Bachmann11, J.J. Back48, A. Badalov36, C. Baesso60, W. Baldini16,38, R.J.
Barlow54, C. Barschel38, S. Barsuk7, W. Barter38, V. Batozskaya28, V.
Battista39, A. Bay39, L. Beaucourt4, J. Beddow51, F. Bedeschi23, I. Bediaga1,
L.J. Bel41, I. Belyaev31, E. Ben-Haim8, G. Bencivenni18, S. Benson38, J.
Benton46, A. Berezhnoy32, R. Bernet40, A. Bertolin22, M.-O. Bettler38, M. van
Beuzekom41, A. Bien11, S. Bifani45, T. Bird54, A. Birnkraut9, A. Bizzeti17,i,
T. Blake48, F. Blanc39, J. Blouw10, S. Blusk59, V. Bocci25, A. Bondar34, N.
Bondar30,38, W. Bonivento15, S. Borghi54, M. Borsato7, T.J.V. Bowcock52, E.
Bowen40, C. Bozzi16, S. Braun11, D. Brett54, M. Britsch10, T. Britton59, J.
Brodzicka54, N.H. Brook46, A. Bursche40, J. Buytaert38, S. Cadeddu15, R.
Calabrese16,f, M. Calvi20,k, M. Calvo Gomez36,p, P. Campana18, D. Campora
Perez38, L. Capriotti54, A. Carbone14,d, G. Carboni24,l, R. Cardinale19,j, A.
Cardini15, P. Carniti20, L. Carson50, K. Carvalho Akiba2,38, R. Casanova
Mohr36, G. Casse52, L. Cassina20,k, L. Castillo Garcia38, M. Cattaneo38, Ch.
Cauet9, G. Cavallero19, R. Cenci23,t, M. Charles8, Ph. Charpentier38, M.
Chefdeville4, S. Chen54, S.-F. Cheung55, N. Chiapolini40, M. Chrzaszcz40,26,
X. Cid Vidal38, G. Ciezarek41, P.E.L. Clarke50, M. Clemencic38, H.V. Cliff47,
J. Closier38, V. Coco38, J. Cogan6, E. Cogneras5, V. Cogoni15,e, L.
Cojocariu29, G. Collazuol22, P. Collins38, A. Comerma-Montells11, A.
Contu15,38, A. Cook46, M. Coombes46, S. Coquereau8, G. Corti38, M. Corvo16,f,
B. Couturier38, G.A. Cowan50, D.C. Craik48, A. Crocombe48, M. Cruz Torres60,
S. Cunliffe53, R. Currie53, C. D’Ambrosio38, J. Dalseno46, P.N.Y. David41, A.
Davis57, K. De Bruyn41, S. De Capua54, M. De Cian11, J.M. De Miranda1, L. De
Paula2, W. De Silva57, P. De Simone18, C.-T. Dean51, D. Decamp4, M.
Deckenhoff9, L. Del Buono8, N. Déléage4, D. Derkach55, O. Deschamps5, F.
Dettori38, B. Dey40, A. Di Canto38, F. Di Ruscio24, H. Dijkstra38, S.
Donleavy52, F. Dordei11, M. Dorigo39, A. Dosil Suárez37, D. Dossett48, A.
Dovbnya43, K. Dreimanis52, G. Dujany54, F. Dupertuis39, P. Durante38, R.
Dzhelyadin35, A. Dziurda26, A. Dzyuba30, S. Easo49,38, U. Egede53, V.
Egorychev31, S. Eidelman34, S. Eisenhardt50, U. Eitschberger9, R. Ekelhof9, L.
Eklund51, I. El Rifai5, Ch. Elsasser40, S. Ely59, S. Esen11, H.M. Evans47, T.
Evans55, A. Falabella14, C. Färber11, C. Farinelli41, N. Farley45, S. Farry52,
R. Fay52, D. Ferguson50, V. Fernandez Albor37, F. Ferrari14, F. Ferreira
Rodrigues1, M. Ferro-Luzzi38, S. Filippov33, M. Fiore16,38,f, M. Fiorini16,f,
M. Firlej27, C. Fitzpatrick39, T. Fiutowski27, P. Fol53, M. Fontana10, F.
Fontanelli19,j, R. Forty38, O. Francisco2, M. Frank38, C. Frei38, M.
Frosini17, J. Fu21, E. Furfaro24,l, A. Gallas Torreira37, D. Galli14,d, S.
Gallorini22,38, S. Gambetta19,j, M. Gandelman2, P. Gandini55, Y. Gao3, J.
García Pardiñas37, J. Garofoli59, J. Garra Tico47, L. Garrido36, D. Gascon36,
C. Gaspar38, U. Gastaldi16, R. Gauld55, L. Gavardi9, G. Gazzoni5, A.
Geraci21,v, D. Gerick11, E. Gersabeck11, M. Gersabeck54, T. Gershon48, Ph.
Ghez4, A. Gianelle22, S. Gianì39, V. Gibson47, L. Giubega29, V.V. Gligorov38,
C. Göbel60, D. Golubkov31, A. Golutvin53,31,38, A. Gomes1,a, C. Gotti20,k, M.
Grabalosa Gándara5, R. Graciani Diaz36, L.A. Granado Cardoso38, E. Graugés36,
E. Graverini40, G. Graziani17, A. Grecu29, E. Greening55, S. Gregson47, P.
Griffith45, L. Grillo11, O. Grünberg63, B. Gui59, E. Gushchin33, Yu. Guz35,38,
T. Gys38, C. Hadjivasiliou59, G. Haefeli39, C. Haen38, S.C. Haines47, S.
Hall53, B. Hamilton58, T. Hampson46, X. Han11, S. Hansmann-Menzemer11, N.
Harnew55, S.T. Harnew46, J. Harrison54, J. He38, T. Head39, V. Heijne41, K.
Hennessy52, P. Henrard5, L. Henry8, J.A. Hernando Morata37, E. van
Herwijnen38, M. Heß63, A. Hicheur2, D. Hill55, M. Hoballah5, C. Hombach54, W.
Hulsbergen41, T. Humair53, N. Hussain55, D. Hutchcroft52, D. Hynds51, M.
Idzik27, P. Ilten56, R. Jacobsson38, A. Jaeger11, J. Jalocha55, E. Jans41, A.
Jawahery58, F. Jing3, M. John55, D. Johnson38, C.R. Jones47, C. Joram38, B.
Jost38, N. Jurik59, S. Kandybei43, W. Kanso6, M. Karacson38, T.M. Karbach38,†,
S. Karodia51, M. Kelsey59, I.R. Kenyon45, M. Kenzie38, T. Ketel42, B.
Khanji20,38,k, C. Khurewathanakul39, S. Klaver54, K. Klimaszewski28, O.
Kochebina7, M. Kolpin11, I. Komarov39, R.F. Koopman42, P. Koppenburg41,38, M.
Korolev32, L. Kravchuk33, K. Kreplin11, M. Kreps48, G. Krocker11, P.
Krokovny34, F. Kruse9, W. Kucewicz26,o, M. Kucharczyk26, V. Kudryavtsev34, K.
Kurek28, T. Kvaratskheliya31, V.N. La Thi39, D. Lacarrere38, G. Lafferty54, A.
Lai15, D. Lambert50, R.W. Lambert42, G. Lanfranchi18, C. Langenbruch48, B.
Langhans38, T. Latham48, C. Lazzeroni45, R. Le Gac6, J. van Leerdam41, J.-P.
Lees4, R. Lefèvre5, A. Leflat32, J. Lefrançois7, O. Leroy6, T. Lesiak26, B.
Leverington11, Y. Li7, T. Likhomanenko65,64, M. Liles52, R. Lindner38, C.
Linn38, F. Lionetto40, B. Liu15, S. Lohn38, I. Longstaff51, J.H. Lopes2, D.
Lucchesi22,r, M. Lucio Martinez37, H. Luo50, A. Lupato22, E. Luppi16,f, O.
Lupton55, F. Machefert7, F. Maciuc29, O. Maev30, S. Malde55, A. Malinin64, G.
Manca15,e, G. Mancinelli6, P. Manning59, A. Mapelli38, J. Maratas5, J.F.
Marchand4, U. Marconi14, C. Marin Benito36, P. Marino23,38,t, R. Märki39, J.
Marks11, G. Martellotti25, M. Martinelli39, D. Martinez Santos42, F. Martinez
Vidal66, D. Martins Tostes2, A. Massafferri1, R. Matev38, A. Mathad48, Z.
Mathe38, C. Matteuzzi20, K. Matthieu11, A. Mauri40, B. Maurin39, A. Mazurov45,
M. McCann53, J. McCarthy45, A. McNab54, R. McNulty12, B. Meadows57, F. Meier9,
M. Meissner11, M. Merk41, D.A. Milanes62, M.-N. Minard4, D.S. Mitzel11, J.
Molina Rodriguez60, S. Monteil5, M. Morandin22, P. Morawski27, A. Mordà6, M.J.
Morello23,t, J. Moron27, A.B. Morris50, R. Mountain59, F. Muheim50, J.
Müller9, K. Müller40, V. Müller9, M. Mussini14, B. Muster39, P. Naik46, T.
Nakada39, R. Nandakumar49, I. Nasteva2, M. Needham50, N. Neri21, S. Neubert11,
N. Neufeld38, M. Neuner11, A.D. Nguyen39, T.D. Nguyen39, C. Nguyen-Mau39,q, V.
Niess5, R. Niet9, N. Nikitin32, T. Nikodem11, D. Ninci23, A. Novoselov35, D.P.
O’Hanlon48, A. Oblakowska-Mucha27, V. Obraztsov35, S. Ogilvy51, O.
Okhrimenko44, R. Oldeman15,e, C.J.G. Onderwater67, B. Osorio Rodrigues1, J.M.
Otalora Goicochea2, A. Otto38, P. Owen53, A. Oyanguren66, A. Palano13,c, F.
Palombo21,u, M. Palutan18, J. Panman38, A. Papanestis49, M. Pappagallo51, L.L.
Pappalardo16,f, C. Parkes54, G. Passaleva17, G.D. Patel52, M. Patel53, C.
Patrignani19,j, A. Pearce54,49, A. Pellegrino41, G. Penso25,m, M. Pepe
Altarelli38, S. Perazzini14,d, P. Perret5, L. Pescatore45, K. Petridis46, A.
Petrolini19,j, M. Petruzzo21, E. Picatoste Olloqui36, B. Pietrzyk4, T.
Pilař48, D. Pinci25, A. Pistone19, S. Playfer50, M. Plo Casasus37, T.
Poikela38, F. Polci8, A. Poluektov48,34, I. Polyakov31, E. Polycarpo2, A.
Popov35, D. Popov10, B. Popovici29, C. Potterat2, E. Price46, J.D. Price52, J.
Prisciandaro39, A. Pritchard52, C. Prouve46, V. Pugatch44, A. Puig Navarro39,
G. Punzi23,s, W. Qian4, R. Quagliani7,46, B. Rachwal26, J.H. Rademacker46, B.
Rakotomiaramanana39, M. Rama23, M.S. Rangel2, I. Raniuk43, N. Rauschmayr38, G.
Raven42, F. Redi53, S. Reichert54, M.M. Reid48, A.C. dos Reis1, S.
Ricciardi49, S. Richards46, M. Rihl38, K. Rinnert52, V. Rives Molina36, P.
Robbe7,38, A.B. Rodrigues1, E. Rodrigues54, J.A. Rodriguez Lopez62, P.
Rodriguez Perez54, S. Roiser38, V. Romanovsky35, A. Romero Vidal37, M.
Rotondo22, J. Rouvinet39, T. Ruf38, H. Ruiz36, P. Ruiz Valls66, J.J. Saborido
Silva37, N. Sagidova30, P. Sail51, B. Saitta15,e, V. Salustino Guimaraes2, C.
Sanchez Mayordomo66, B. Sanmartin Sedes37, R. Santacesaria25, C. Santamarina
Rios37, E. Santovetti24,l, A. Sarti18,m, C. Satriano25,n, A. Satta24, D.M.
Saunders46, D. Savrina31,32, M. Schiller38, H. Schindler38, M. Schlupp9, M.
Schmelling10, T. Schmelzer9, B. Schmidt38, O. Schneider39, A. Schopper38,
M.-H. Schune7, R. Schwemmer38, B. Sciascia18, A. Sciubba25,m, A. Semennikov31,
I. Sepp53, N. Serra40, J. Serrano6, L. Sestini22, P. Seyfert11, M. Shapkin35,
I. Shapoval16,43,f, Y. Shcheglov30, T. Shears52, L. Shekhtman34, V.
Shevchenko64, A. Shires9, R. Silva Coutinho48, G. Simi22, M. Sirendi47, N.
Skidmore46, I. Skillicorn51, T. Skwarnicki59, E. Smith55,49, E. Smith53, J.
Smith47, M. Smith54, H. Snoek41, M.D. Sokoloff57,38, F.J.P. Soler51, F.
Soomro39, D. Souza46, B. Souza De Paula2, B. Spaan9, P. Spradlin51, S.
Sridharan38, F. Stagni38, M. Stahl11, S. Stahl38, O. Steinkamp40, O.
Stenyakin35, F. Sterpka59, S. Stevenson55, S. Stoica29, S. Stone59, B.
Storaci40, S. Stracka23,t, M. Straticiuc29, U. Straumann40, R. Stroili22, L.
Sun57, W. Sutcliffe53, K. Swientek27, S. Swientek9, V. Syropoulos42, M.
Szczekowski28, P. Szczypka39,38, T. Szumlak27, S. T’Jampens4, T. Tekampe9, M.
Teklishyn7, G. Tellarini16,f, F. Teubert38, C. Thomas55, E. Thomas38, J. van
Tilburg41, V. Tisserand4, M. Tobin39, J. Todd57, S. Tolk42, L. Tomassetti16,f,
D. Tonelli38, S. Topp-Joergensen55, N. Torr55, E. Tournefier4, S. Tourneur39,
K. Trabelsi39, M.T. Tran39, M. Tresch40, A. Trisovic38, A. Tsaregorodtsev6, P.
Tsopelas41, N. Tuning41,38, A. Ukleja28, A. Ustyuzhanin65,64, U. Uwer11, C.
Vacca15,e, V. Vagnoni14, G. Valenti14, A. Vallier7, R. Vazquez Gomez18, P.
Vazquez Regueiro37, C. Vázquez Sierra37, S. Vecchi16, J.J. Velthuis46, M.
Veltri17,h, G. Veneziano39, M. Vesterinen11, B. Viaud7, D. Vieira2, M. Vieites
Diaz37, X. Vilasis-Cardona36,p, A. Vollhardt40, D. Volyanskyy10, D. Voong46,
A. Vorobyev30, V. Vorobyev34, C. Voß63, J.A. de Vries41, R. Waldi63, C.
Wallace48, R. Wallace12, J. Walsh23, S. Wandernoth11, J. Wang59, D.R. Ward47,
N.K. Watson45, D. Websdale53, A. Weiden40, M. Whitehead48, D. Wiedner11, G.
Wilkinson55,38, M. Wilkinson59, M. Williams38, M.P. Williams45, M. Williams56,
F.F. Wilson49, J. Wimberley58, J. Wishahi9, W. Wislicki28, M. Witek26, G.
Wormser7, S.A. Wotton47, S. Wright47, K. Wyllie38, Y. Xie61, Z. Xu39, Z.
Yang3, X. Yuan34, O. Yushchenko35, M. Zangoli14, M. Zavertyaev10,b, L. Zhang3,
Y. Zhang3, A. Zhelezov11, A. Zhokhov31, L. Zhong3.

1Centro Brasileiro de Pesquisas Físicas (CBPF), Rio de Janeiro, Brazil
2Universidade Federal do Rio de Janeiro (UFRJ), Rio de Janeiro, Brazil
3Center for High Energy Physics, Tsinghua University, Beijing, China
4LAPP, Université Savoie Mont-Blanc, CNRS/IN2P3, Annecy-Le-Vieux, France
5Clermont Université, Université Blaise Pascal, CNRS/IN2P3, LPC, Clermont-
Ferrand, France
6CPPM, Aix-Marseille Université, CNRS/IN2P3, Marseille, France
7LAL, Université Paris-Sud, CNRS/IN2P3, Orsay, France
8LPNHE, Université Pierre et Marie Curie, Université Paris Diderot,
CNRS/IN2P3, Paris, France
9Fakultät Physik, Technische Universität Dortmund, Dortmund, Germany
10Max-Planck-Institut für Kernphysik (MPIK), Heidelberg, Germany
11Physikalisches Institut, Ruprecht-Karls-Universität Heidelberg, Heidelberg,
Germany
12School of Physics, University College Dublin, Dublin, Ireland
13Sezione INFN di Bari, Bari, Italy
14Sezione INFN di Bologna, Bologna, Italy
15Sezione INFN di Cagliari, Cagliari, Italy
16Sezione INFN di Ferrara, Ferrara, Italy
17Sezione INFN di Firenze, Firenze, Italy
18Laboratori Nazionali dell’INFN di Frascati, Frascati, Italy
19Sezione INFN di Genova, Genova, Italy
20Sezione INFN di Milano Bicocca, Milano, Italy
21Sezione INFN di Milano, Milano, Italy
22Sezione INFN di Padova, Padova, Italy
23Sezione INFN di Pisa, Pisa, Italy
24Sezione INFN di Roma Tor Vergata, Roma, Italy
25Sezione INFN di Roma La Sapienza, Roma, Italy
26Henryk Niewodniczanski Institute of Nuclear Physics Polish Academy of
Sciences, Kraków, Poland
27AGH - University of Science and Technology, Faculty of Physics and Applied
Computer Science, Kraków, Poland
28National Center for Nuclear Research (NCBJ), Warsaw, Poland
29Horia Hulubei National Institute of Physics and Nuclear Engineering,
Bucharest-Magurele, Romania
30Petersburg Nuclear Physics Institute (PNPI), Gatchina, Russia
31Institute of Theoretical and Experimental Physics (ITEP), Moscow, Russia
32Institute of Nuclear Physics, Moscow State University (SINP MSU), Moscow,
Russia
33Institute for Nuclear Research of the Russian Academy of Sciences (INR RAN),
Moscow, Russia
34Budker Institute of Nuclear Physics (SB RAS) and Novosibirsk State
University, Novosibirsk, Russia
35Institute for High Energy Physics (IHEP), Protvino, Russia
36Universitat de Barcelona, Barcelona, Spain
37Universidad de Santiago de Compostela, Santiago de Compostela, Spain
38European Organization for Nuclear Research (CERN), Geneva, Switzerland
39Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
40Physik-Institut, Universität Zürich, Zürich, Switzerland
41Nikhef National Institute for Subatomic Physics, Amsterdam, The Netherlands
42Nikhef National Institute for Subatomic Physics and VU University Amsterdam,
Amsterdam, The Netherlands
43NSC Kharkiv Institute of Physics and Technology (NSC KIPT), Kharkiv, Ukraine
44Institute for Nuclear Research of the National Academy of Sciences (KINR),
Kyiv, Ukraine
45University of Birmingham, Birmingham, United Kingdom
46H.H. Wills Physics Laboratory, University of Bristol, Bristol, United
Kingdom
47Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom
48Department of Physics, University of Warwick, Coventry, United Kingdom
49STFC Rutherford Appleton Laboratory, Didcot, United Kingdom
50School of Physics and Astronomy, University of Edinburgh, Edinburgh, United
Kingdom
51School of Physics and Astronomy, University of Glasgow, Glasgow, United
Kingdom
52Oliver Lodge Laboratory, University of Liverpool, Liverpool, United Kingdom
53Imperial College London, London, United Kingdom
54School of Physics and Astronomy, University of Manchester, Manchester,
United Kingdom
55Department of Physics, University of Oxford, Oxford, United Kingdom
56Massachusetts Institute of Technology, Cambridge, MA, United States
57University of Cincinnati, Cincinnati, OH, United States
58University of Maryland, College Park, MD, United States
59Syracuse University, Syracuse, NY, United States
60Pontifícia Universidade Católica do Rio de Janeiro (PUC-Rio), Rio de
Janeiro, Brazil, associated to 2
61Institute of Particle Physics, Central China Normal University, Wuhan,
Hubei, China, associated to 3
62Departamento de Fisica , Universidad Nacional de Colombia, Bogota, Colombia,
associated to 8
63Institut für Physik, Universität Rostock, Rostock, Germany, associated to 11
64National Research Centre Kurchatov Institute, Moscow, Russia, associated to
31
65Yandex School of Data Analysis, Moscow, Russia, associated to 31
66Instituto de Fisica Corpuscular (IFIC), Universitat de Valencia-CSIC,
Valencia, Spain, associated to 36
67Van Swinderen Institute, University of Groningen, Groningen, The
Netherlands, associated to 41

aUniversidade Federal do Triângulo Mineiro (UFTM), Uberaba-MG, Brazil
bP.N. Lebedev Physical Institute, Russian Academy of Science (LPI RAS),
Moscow, Russia
cUniversità di Bari, Bari, Italy
dUniversità di Bologna, Bologna, Italy
eUniversità di Cagliari, Cagliari, Italy
fUniversità di Ferrara, Ferrara, Italy
gUniversità di Firenze, Firenze, Italy
hUniversità di Urbino, Urbino, Italy
iUniversità di Modena e Reggio Emilia, Modena, Italy
jUniversità di Genova, Genova, Italy
kUniversità di Milano Bicocca, Milano, Italy
lUniversità di Roma Tor Vergata, Roma, Italy
mUniversità di Roma La Sapienza, Roma, Italy
nUniversità della Basilicata, Potenza, Italy
oAGH - University of Science and Technology, Faculty of Computer Science,
Electronics and Telecommunications, Kraków, Poland
pLIFAELS, La Salle, Universitat Ramon Llull, Barcelona, Spain
qHanoi University of Science, Hanoi, Viet Nam
rUniversità di Padova, Padova, Italy
sUniversità di Pisa, Pisa, Italy
tScuola Normale Superiore, Pisa, Italy
uUniversità degli Studi di Milano, Milano, Italy
vPolitecnico di Milano, Milano, Italy
†Deceased


# Plug-and-Play Secondary Control for Safety of LTI Systems under Attacks

Yankai Lin  Michelle S. Chong  Carlos Murguia  Department of Mechanical
Engineering, Eindhoven University of Technology, the Netherlands
(e-mail: {y.lin2,m.s.t.chong,c.g.murguia}@ tue.nl).

###### Abstract

We consider the problem of controller design for linear time-invariant cyber-
physical systems (CPSs) controlled via networks. Specifically, we adopt the
set-up that a controller has already been designed to stabilize the plant.
However, the closed loop system may be subject to actuator and sensor attacks.
We first perform a reachability analysis to see the effect of potential
attacks. To further ensure the safety of the states of the system, we choose a
subset of sensors that can be locally secured and made free of attacks. Using
these limited resources, an extra controller is designed to enhance the safety
of the new closed loop. The safety of the system will be characterized by the
notion of safe sets. Lyapunov based analysis will be used to derive sufficient
conditions that ensure the states always stay in the safe set. The conditions
will then be stated as convex optimization problems which can be solved
efficiently. Lastly, our theoretical results are illustrated through numerical
simulations.

###### keywords:

Safety, Cyber-physical systems, Control systems security, Linear systems,
Linear matrix inequality.

##  1 Introduction

Many industrial systems are now remotely controlled over communication
networks which makes the operation and maintenance of these systems more
efficient. However, due to the involvement of the communication channel,
cyber-physical systems or CPS (systems that integrate the physical process
with communication networks) are often subject to potential attacks which can
lead to disastrous outcomes. Examples include the StuxNet malware incident and
many other incidents shown in Cárdenas et al. (2008). Thus, preventing
disastrous outcomes or in other words, the safety of CPS is of significant
importance.

Various cyber attacks have been modelled and analyzed in Teixeira et al.
(2015a) using tools from systems theory. A core idea presented in Teixeira et
al. (2015a) is that, the amount of resources needed for attackers to launch an
attack against a CPS varies according to the complexity of attacks. In most
cases, the attacker wants to keep attacking the system while remaining
stealthy. This can be accomplished by the adversary via well-designed attacks
with sufficient resources. In this case, normal fault detection and isolation
strategies may fail to work Dutta and Langbort (2017). As detailed later in
this paper, we work with attack signals with bounded energy to reflect the
fact that adversaries might aim to remain stealthy and that attacks tampering
with sensing, actuation, and networked data always have limited resources due
to physical constraints.

One way to deal with attacks on CPS is to design a control scheme that can
detect the presence of attacks and then mitigate it, see Fawzi et al. (2014);
Chong et al. (2015) for results on linear systems. This approach ensures
security using redundancy and typically requires a large number of observers
to completely nullify the effect of attacks. In this work, however, we follow
a different approach. We aim to ensure safety of the system when attacks are
present. Namely, by modifying the dynamics of the closed-loop control system
we want the states of interest to stay in a set where safe operation of the
system can be guaranteed in the presence of stealthy attacks. This idea has
been explored in Murguia et al. (2020), where the volume of the set of states
that attacks can induce in the system is used as a security metric. After
appropriate modification of system dynamics, reachable sets of the attack
signal should be made as small as possible and fully contained in the safe
set. As a result, to compromise the safety of the system, an attacker needs to
invest more resources. One possible way is to inject attack signals with
larger amplitudes which will be used as a measure of safety of the system
against potential attacks. For other security measures of CPS in different
settings, we refer interested readers to recent works Bai et al. (2015); Tang
et al. (2019); Teixeira et al. (2015b); Milošević et al. (2020) and references
therein.

The main contribution of this work is, via set-theoretic methods, a systematic
tool to verify the safety of a given linear time-invariant closed-loop control
system (a primary controller in connection with the plant) subject to sensor
and actuator attacks. In addition, we also provide a tractable method of
synthesizing an output feedback linear controller, which we call a secondary
controller, that is connected concurrently to the existing closed-loop control
system to ensure safety of the overall closed-loop control system. We depict
this setup in Figure 1. This is different from control barrier function based
methods for safety-critical systems in Ames et al. (2017), where state
feedback control laws are calculated .

Figure 1: Ensuring safety with a secondary controller.

The rest of this manuscript is organized as follows. The problem formulation
is given in Section 2. Section 3 presents the analysis tool to evaluate the
safety of the closed loop system given the pre-designed controller. Section 4
proposes optimization-based controller synthesis tools to modify the given
system aiming to enhance safety of the overall system. Section 5 illustrates
the main results via numerical simulations. We will conclude the paper in
Section 6 and provide some future research directions.

Notation: Let $\mathbb{R}$ be the set of real numbers and $\mathbb{R}^{n}$ be
the $n$-dimensional Euclidean space. The matrix $I_{n}$ is used to denote the
$n$-dimensional identity matrix and $n$ will be omitted when the dimension is
clear. Similarly, $\mathbf{0}$ denotes the zero matrix with appropriate
dimensions. For a given vector $x$, For a given square matrix $R$,
$\text{Tr}[R]$ denotes the trace of $R$. We use $A\succ 0$ ($A\prec 0$) and
$A\succeq 0$ ($A\preceq 0$) to denote the matrix $A$ is positive (negative)
definite and positive (negative) semidefinite, respectively.

##  2 Problem Formulation

We consider the setup shown in Figure 1, where the plant is modelled as a
linear time-invariant (LTI) system as follows

$\begin{split}\dot{x}_{p}&=A_{p}x_{p}+B_{p}u\\\ y&=C_{p}x_{p},\end{split}$ (1)

where $x_{p}\in\mathbb{R}^{n_{p}}$ is the state vector of the plant,
$u\in\mathbb{R}^{n_{u}}$ and $y\in\mathbb{R}^{n_{y}}$ are the input and the
output of the plant, respectively. We assume that there is a controller which
has already been designed to stabilize the plant (1),

$\begin{split}\dot{x}_{1}&=A_{1}x_{1}+B_{1}(y+a_{y})\\\
u_{P}&=C_{1}x_{1}+D_{1}(y+a_{y})+a_{u},\end{split}$ (2)

which we call the primary controller of system (1). It is pre-designed to
stabilize (1) with the input signal $u_{P}\in\mathbb{R}^{n_{u}}$ and subject
to potential cyber attacks denoted by $a=[a_{u}^{\top}\
a_{y}^{\top}]^{\top}\in\mathbb{R}^{n_{u}+n_{y}}$, where $a_{u}$ and $a_{y}$
denote actuator and sensor attacks respectively. Since the primary controller
(2) is pre-designed without being aware of the attacks, the security and
safety of the closed loop may be compromised. Therefore, we propose
introducing a secondary controller, that runs in conjunction with the primary
controller. The secondary controller uses a subset of the sensors and
actuators, which are either available locally or known to be safeguarded
against malicious manipulation (e.g., via encryption or watermarking). We
propose designing a secondary controller which takes the following form

$\begin{split}\dot{x}_{2}&=A_{2}x_{2}+B_{2}y_{S}\\\
u_{S}&=C_{2}x_{2}+D_{2}y_{S},\end{split}$ (3)

where $y_{S}=C_{S}y$ is the set of sensors that are available to the secondary
controller, and $u_{S}$ is the control input signal. The overall input is
given by

$u=u_{P}+E_{u}u_{S},$ (4)

where $E_{u}$ is the selection matrix that dictates the set of secondary
control inputs that will be fed back to the plant (1). Note that the secondary
controller uses uncompromised sensor(s) and actuator(s) only, as it is co-
located with the plant (i.e., not connected to the network) and we assume that
the local sensors and actuators are not subject to network attacks.
Consequently, no attack signals appear in (3). The goal of the secondary
controller (3) is to ensure that when the overall closed loop system (1)-(4)
is subject to cyber attacks, the safety of the closed loop can be ensured,
i.e., the trajectories of the states of interest of the closed-loop system
(1)-(4) remain within a given safe set for all $t\geq 0$.

We describe our given safe set by the following ellipsoid
$\mathcal{E}(R_{\zeta},\bar{\zeta})$:

$(\zeta-\bar{\zeta})^{\top}R_{\zeta}(\zeta-\bar{\zeta})\leq 1,$ (5)

where $\zeta=[x_{p}^{\top}\ x_{1}^{\top}\
x_{2}^{\top}]^{\top}:=[\zeta_{1}^{\top}\ x_{2}^{\top}]^{\top}$ and
$\bar{\zeta}$ is defined similarly. Moreover, we have $R_{\zeta}\succeq 0$.
Ellipsoidal safe set has also been used in existing literature for security
and safety, see Romagnoli et al. (2020) for example. If the safe set is of
other shapes, (5) can be used as an outer approximation of the true safe set.
Depending on different scenarios, the matrix $R_{\zeta}$ may be rank
deficient. In this work, our primary concern are the plant states $x_{p}$.

We aim to solve the following problems in this paper.

  1. (Problem I)

Assess the worst case attacks $a_{u}$ and $a_{y}$ that can be tolerated by
connecting only the primary controller (2) to the plant (1).

  2. (Problem II)

Synthesize the secondary controller (3) such that the attacker needs to invest
more resources to violate the safety condition (5).

We provide a solution to Problem I and II in Section 3 and 4, respectively.

##  3 Invariant Set Based Safety Analysis

The first step of our analysis is to assess the worst possible attacks that
can be tolerated while ensuring safety of the closed loop system. In general,
when dealing with cyber attacks, it is natural to assume that the attacks are
unbounded. However, depending on the purpose of the attacker, intelligent
attackers often seek to remain stealthy and undetected, see Teixeira et al.
(2015a). Thus, we use the following condition of norm bounded attack signals
to reflect the fact that the attacker is resource limited and wants to remain
undetected, i.e., the attack signals $a$ satisfy the following

$a^{\top}R_{a}a\leq 1,$ (6)

with a positive definite matrix $R_{a}$. Larger valued entries in $R_{a}$
imply smaller upper bound on the corresponding attack signal $a$.

######  Remark 1.

Our analysis is based on invariant set analysis. Theoretically, one should
attempt to find the reachable set of the closed loop system subject to attacks
and check if the reachable set is a subset of the safe set. However,
calculating the exact reachable set of a given system is difficult in general.
Thus, we seek an invariant set of the closed loop which can be used as an
outer approximation of the reachable set, see Blanchini (1999) and Escudero et
al. (2022b).

The closed loop system consisting of the plant (1), primary controller (2),
and the secondary controller (3) can be written as

$\dot{\zeta}=\mathcal{A}\zeta+\mathcal{B}a,$ (7)

where

$\begin{split}&\mathcal{A}=\\\
&\left[\begin{array}[]{ccc}A_{p}+B_{p}D_{1}C_{p}+B_{p}E_{u}D_{2}C_{S}C_{p}&B_{p}C_{p}&B_{p}E_{u}C_{2}\\\
B_{1}C_{p}&A_{1}&\mathbf{0}\\\
B_{2}C_{S}C_{p}&\mathbf{0}&A_{2}\end{array}\right],\end{split}$ (8)
$\mathcal{B}=\left[\begin{array}[]{cc}B_{p}&B_{p}D_{1}\\\ \mathbf{0}&B_{1}\\\
\mathbf{0}&\mathbf{0}\end{array}\right].$ (9)

Moreover, we define the following partitions,

$\mathcal{A}:=\left[\begin{array}[]{cc}\mathcal{A}_{1}&\mathcal{A}_{2}\\\
\mathcal{A}_{3}&\mathcal{A}_{4}\end{array}\right],$ (10)
$\mathcal{B}:=\left[\begin{array}[]{c}\mathcal{B}_{1}\\\
\mathcal{B}_{2}\end{array}\right],$ (11)

where

$\mathcal{A}_{1}=\left[\begin{array}[]{cc}A_{p}+B_{p}D_{1}C_{p}+B_{p}E_{u}D_{2}C_{S}C_{p}&B_{p}C_{p}\\\
B_{1}C_{p}&A_{1}\end{array}\right],$ $\mathcal{A}_{4}=A_{2},$

and the expressions of $\mathcal{A}_{2}$ and $\mathcal{A}_{3}$ follow from
(8). Similarly,

$\mathcal{B}_{1}=\left[\begin{array}[]{cc}B_{p}&B_{p}D_{1}\\\
\mathbf{0}&B_{1}\end{array}\right],$

with $\mathcal{B}_{2}$ being $\mathbf{0}$.

Since we want to first find the worst attack signals that can be rendered in
the safe set by the primary controller (2) only, we set $E_{u}=\mathbf{0}$ and
$C_{S}=\mathbf{0}$. Moreover, because the secondary controller (3) is not
connected in feedback with the plant (1) in this case, entries corresponding
to secondary controller states are set to $\mathbf{0}$.

We observe that the closed-loop system (7) is a LTI system driven by the
attack signal $a$. We aim to construct an ellipsoidal outer approximation
$\mathcal{E}_{a}$ such that $\mathcal{R}_{a}\subseteq\mathcal{E}_{a}$ where
$\mathcal{R}_{a}$ is the reachable set of (7). We will solve this problem via
Lyapunov analysis and optimization tools.

Consider the quadratic function $V(\zeta)=\zeta^{\top}P\zeta$ for a positive
semi-definite matrix $P\succeq 0$. If we can find a $P$ such that $\dot{V}\leq
0$ whenever $V\geq 1$ along the trajectories of (7), then the ellipsoid
$\mathcal{E}_{P}=\\{\zeta\ |\ \zeta^{\top}P\zeta\leq 1\\}$ is an invariant set
of (7) as the states starting inside $\mathcal{E}_{P}$ cannot leave
$\mathcal{E}_{P}$.

First, for a matrix $Q$ with appropriate dimension, we define the following
matrices,

$E_{1}=\left[\begin{array}[]{ccc}\mathcal{A}_{1}^{\top}Q+Q\mathcal{A}_{1}&\mathbf{0}&Q\mathcal{B}_{1}\\\
\star&\mathbf{0}&\mathbf{0}\\\ \star&\star&\mathbf{0}\end{array}\right],$ (12)
$S=\left[\begin{array}[]{ccc}\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\star&1&\mathbf{0}\\\ \star&\star&-R_{a}\end{array}\right].$ (13)

Since in this part, the secondary controller has no impact on the closed loop
system, it is reasonable to work with only the projection of sets on to the
$\zeta_{1}$-hyperplane. We denote the projection of
$\mathcal{E}(R_{\zeta},\bar{\zeta})$ onto the $\zeta_{1}$-hyperplane by
$\hat{\mathcal{E}}(R_{\zeta},\bar{\zeta})$.

Then, we have the following result.

######  Proposition 2.

Let $E_{u}=\mathbf{0}$ and $C_{S}=\mathbf{0}$. Given a positive definite
matrix $R_{a}$, if there exist a positive definite matrix $Q$ and non-negative
constants $\alpha,\beta$ satisfying the following inequalities:

$-E_{1}-\alpha F-\beta S\succeq 0,\ s.t.\
\mathcal{E}_{Q}\subseteq\mathcal{E}(R_{\zeta},\bar{\zeta}),$ (14)

where $\mathcal{E}_{Q}=\\{\zeta_{1}\ |\ \zeta_{1}^{\top}Q\zeta_{1}\leq 1\\}$,
then there exists a forward invariant set for (7) under the constraints (6)
and it is a subset of the safe set $\mathcal{E}(R_{\zeta},\bar{\zeta})$.

We consider the vector $[\zeta_{1}^{\top}\ 1\ a^{\top}]^{\top}$. The condition
$\dot{V}(\zeta_{1})\leq 0$ can be restated as

$E_{1}=\left[\begin{array}[]{ccc}\mathcal{A}_{1}^{\top}Q+Q\mathcal{A}_{1}&\mathbf{0}&Q\mathcal{B}_{1}\\\
\star&\mathbf{0}&\mathbf{0}\\\ \star&\star&\mathbf{0}\end{array}\right]\preceq
0.$ (15)

Similarly, we can restate (6) as

$S=\left[\begin{array}[]{ccc}\mathbf{0}&\mathbf{0}&\mathbf{0}\\\
\star&1&\mathbf{0}\\\ \star&\star&-R_{a}\end{array}\right]\succeq 0.$ (16)

Finally, the condition $V(\zeta_{1})\geq 1$ can be written as

$F=\left[\begin{array}[]{ccc}Q&\mathbf{0}&\mathbf{0}\\\ \star&-1&\mathbf{0}\\\
\star&\star&\mathbf{0}\end{array}\right]\succeq 0.$ (17)

Since we want to ensure that $\dot{V}(\zeta_{1})\leq 0$ holds when (16) and
(17) hold, we can apply the $\mathcal{S}$-procedure (see Section 2.6.3 of Boyd
et al. (1994) and Lemma 1 of Escudero et al. (2022a)), which states that there
should exist non-negative constants $\alpha$ and $\beta$ such that the
following holds

$-E_{1}-\alpha F-\beta S\succeq 0.$

The extra constraint
$\mathcal{E}_{Q}\subseteq\hat{\mathcal{E}}(R_{\zeta},\bar{\zeta})$ ensures
that the resulting invariant set $\mathcal{E}_{Q}$ is a subset of the safe set
$\mathcal{E}(R_{\zeta},\bar{\zeta})$ since the secondary controller states
$x_{2}$ have no impact on the overall system. The invariant set can be
constructed as $\\{\zeta\ |\ \zeta^{\top}P\zeta\leq 1\\}$ where $Q$ is the
projection of $P$ onto the $\zeta_{1}$-hyperplane.

######  Remark 3.

Note that in Proposition 2, the condition (14) is in fact not a linear matrix
inequality (LMI) since both the matrix $Q$ and the constant $\alpha$ are
variables. In practice, one can choose to solve (14) for a fixed $\alpha\geq
0$.

######  Remark 4.

The result given in Proposition 2 is only a sufficient condition for finding
an invariant set for (7). We are restricting the shape of the invariant set to
be an ellipsoid for easier analysis. As a result, if there does not exist a
$Q,\alpha$, and $\beta$ that satisfy the conditions in Proposition 2, there
may still exist an invariant set for (7) which is a subset of the safe set.

######  Remark 5.

The constraint $\mathcal{E}_{Q}\subseteq\mathcal{E}(R_{\zeta},\bar{\zeta})$ is
a special case of the problem of inner approximation of the intersection of
ellipsoids. It is shown in Section 3.7.3 of Boyd et al. (1994) that via the
$\mathcal{S}$-procedure, the problem can be equivalently written as an LMI
problem which is convex.

Proposition 2 gives a sufficient condition to check, whether it is possible to
keep the states of the closed loop system (7) are always inside the safe set,
given the resources which are available to the attacker (characterized by
$R_{a}$). To find out the worst possible attack that can be dealt with by the
primary controller alone, one way is to maximize the volume of the ellipsoid
induced by (6). Since $R_{a}$ is a positive definite matrix, this is
equivalent to minimizing the volume of the ellipsoid induced by
$a^{\top}R_{a}^{-1}a\leq 1$. It is shown in Kurzhanski and Vályi (1997) that
the volume of the latter ellipsoid is proportional to $\sqrt{\det[R_{a}]}$,
where $\det$ means the determinant of a matrix. The function is again shown to
share the same minimizers with the function $\log\det[R_{a}]$, see Section 3.7
of Boyd et al. (1994). Unfortunately, for positive definite $R_{a}$,
$\log\det[R_{a}]$ is non-convex. This makes us follow the approach used in
Murguia et al. (2020) by minimizing a convex upper bound of
$\sqrt{\det[R_{a}]}$.

######  Lemma 6 (Lemma 4, Murguia et al. (2020)).

Given any positive definite matrix $R\in\mathbb{R}^{n\times n}$, the following
inequality holds

$\det[R]^{\frac{1}{2}}\leq\frac{1}{\sqrt{n}^{n}}\text{Tr}[R]^{\frac{n}{2}}.$
(18)

Moreover, $\arg\min[\text{Tr}[R]^{\frac{n}{2}}]=\arg\min[\text{Tr}[R]]$.

See Lemma 4 of Murguia et al. (2020).

The first result in this paper follows from Proposition 2 and is presented
below.

######  Theorem 7.

Let $E_{u}=\mathbf{0}$ and $C_{S}=\mathbf{0}$. If there exist a positive
definite matrix $R_{a}$, a positive definite matrix $Q$, and non-negative
constants $\alpha,\beta$ that solve the following optimization problem:

$\begin{split}\underset{Q,R_{a},\alpha,\beta}{\min}&\ \text{Tr}[R_{a}]\\\
s.t.&\ -E_{1}-\alpha F-\beta S\succeq 0\\\ &\
\mathcal{E}_{Q}\subseteq\hat{\mathcal{E}}(R_{\zeta},\bar{\zeta}),\end{split}$
(19)

then for all attack signals satisfying (6), $\mathcal{E}_{P}$ is a subset of
the safe set $\mathcal{E}(R_{\zeta},\bar{\zeta})$ and is a forward invariant
set for (7), where $Q$ is the projection of $P$ on to the
$\zeta_{1}$-hyperplane.

######  Remark 8.

Note that, in Theorem 7, $R_{a}$ is also a variable while it is given in
Proposition 2. Consequently, to make the optimization problem (19) convex, it
should be solved for fixed constants $\alpha$ and $\beta$.

##  4 Secondary Controller Synthesis

We now address the problem of secondary controller synthesis such that the
attacker needs to invest more resource to violate the safety condition (5).
The sensor and input selection matrices $C_{S}$ and $E_{u}$ are assumed to be
pre-selected. That is, the designer first decides which set of sensors are
locally available such that they can be secured and which set of secondary
controller inputs are to be fed back to the plant. Specifically, given $C_{S}$
and $E_{u}$, we want to find $\kappa:=(A_{2},B_{2},C_{2},D_{2})$ such that
$\text{Tr}[R_{a}]$ is minimized while there exists an invariant ellipsoid that
is contained in the safe set represented by (5). Let us define

$E_{2}=\left[\begin{array}[]{ccc}\mathcal{A}^{\top}P+P\mathcal{A}&\mathbf{0}&P\mathcal{B}\\\
\star&\mathbf{0}&\mathbf{0}\\\ \star&\star&\mathbf{0}\end{array}\right].$ (20)

The synthesis problem can be formulated as follows.

$\begin{split}\underset{P,R_{a},\alpha,\beta,\kappa}{\min}&\
\text{Tr}[R_{a}]\\\ s.t.&\ -E_{2}-\alpha F-\beta S\succeq 0\\\ &\
\mathcal{E}_{P}\subseteq\mathcal{E}(R_{\zeta},\bar{\zeta}),\end{split}$ (21)

Note that, the problem (21) is non-convex since we have introduced the set of
new variables $\kappa$ and there exist products of matrix variables such as
$\mathcal{A}^{\top}P$. This is different from the conditions in Theorem 7
where the matrix $\mathcal{A}_{1}$ is given, making (15) linear in
optimization variables. Therefore, we follow the approach introduced in
Scherer et al. (1997) to find an invertible change of variables that convexify
the problem (21).

For a positive definite matrix $P$, let it take the following form

$P=\left[\begin{array}[]{cc}Y&N\\\ \star&\tilde{Y}\end{array}\right],\
P^{-1}=\left[\begin{array}[]{cc}X&M\\\ \star&\tilde{X}\end{array}\right].$
(22)

Let the dimensions of $\mathcal{A}_{1}$ and $\mathcal{A}_{4}$ be $n_{1}$ and
$n_{2}$, respectively. Then we assume $Y$ and $X$ are both in
$\mathbb{R}^{n_{1}\times n_{1}}$ and $\tilde{Y}$ and $\tilde{X}$ are both in
$\mathbb{R}^{n_{2}\times n_{2}}$. Moreover, they are positive definite. Via
some calculations detailed in Scherer et al. (1997), we have the following
equality

$P\Pi_{1}=\Pi_{2},$ (23)

where

$\Pi_{1}=\left[\begin{array}[]{cc}X&I\\\
M^{\top}&\mathbf{0}\end{array}\right],\
\Pi_{2}=\left[\begin{array}[]{cc}I&Y\\\
\mathbf{0}&N^{\top}\end{array}\right].$ (24)

To better present how the change of variable works, we define the following
matrices

$\begin{split}\hat{A}&=\left[\begin{array}[]{cc}A_{p}+B_{p}D_{1}C_{p}&B_{p}C_{p}\\\
B_{1}C_{p}&A_{1}\end{array}\right],\
\hat{B}=\left[\begin{array}[]{c}B_{p}E_{u}\\\ \mathbf{0}\end{array}\right],\\\
&\quad\quad\quad\quad\ \ \ \ \ \ \ \
\hat{C}=\left[\begin{array}[]{cc}C_{S}C_{p}&\mathbf{0}\end{array}\right].\end{split}$
(25)

It is obvious that the matrices defined in (25) are all given. Moreover, it
can be verified that

$\mathcal{A}_{1}=\hat{A}+\hat{B}D_{2}\hat{C},\
\mathcal{A}_{2}=\hat{B}C_{2},\mathcal{A}_{3}=B_{2}\hat{C},\
\mathcal{A}_{4}=A_{2}.$ (26)

Then, the change of variables are given as follows

$\begin{split}\mathbf{A}&=Y\mathcal{A}_{1}X+Y\mathcal{A}_{2}M^{\top}+N\mathcal{A}_{3}X+N\mathcal{A}_{4}M^{\top}\\\
\mathbf{B}&=NB_{2}+Y\hat{B}D_{2}\\\ \mathbf{C}&=C_{2}M^{\top}+D_{2}\hat{C}X\\\
\mathbf{D}&=D_{2}.\end{split}$ (27)

From (24) and (27), after some calculations it can be derived that

$\begin{split}\Pi_{1}^{\top}P\mathcal{A}\Pi_{1}=\Pi_{2}^{\top}\mathcal{A}\Pi_{1}&=\left[\begin{array}[]{cc}\hat{A}X+\hat{B}\mathbf{C}&\hat{A}+\hat{B}\mathbf{D}\hat{C}\\\
\mathbf{A}&Y\hat{A}+\mathbf{B}\hat{C}\end{array}\right]\\\
:&=A(\eta),\end{split}$ (28)
$\Pi_{1}^{\top}P\mathcal{B}=\Pi_{2}^{\top}\mathcal{B}=\left[\begin{array}[]{c}\mathcal{B}_{1}\\\
Y\mathcal{B}_{1}\end{array}\right]:=B(\eta).$ (29)

Note that, (28) and (29) are linear in the new variables
$\eta:=(X,Y,\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})$. Based on (28) and
(29), if we perform the congruence transformation with
$\text{diag}(\Pi_{1},I)$ on (20), we have

$\mathbf{E}_{2}=\left[\begin{array}[]{ccc}A(\eta)^{\top}+A(\eta)&\mathbf{0}&B(\eta)\\\
\star&\mathbf{0}&\mathbf{0}\\\ \star&\star&\mathbf{0}\end{array}\right].$ (30)

After applying the same transformation to $F$ and $S$ in (21), we have

$\mathbf{F}=\left[\begin{array}[]{ccc}P(\eta)&\mathbf{0}&\mathbf{0}\\\
\star&-1&\mathbf{0}\\\ \star&\star&\mathbf{0}\end{array}\right],\
\mathbf{S}=S,$ (31)

where

$P(\eta)=\Pi_{1}^{\top}P\Pi_{1}=\left[\begin{array}[]{cc}X&I\\\
I&Y\end{array}\right].$ (32)

Next, we handle the other constraint in (21). To simplify the analysis, we
make the following assumption.

######  Assumption 9.

Let $Q$ be the projection of $P$ on to the $\zeta_{1}$-hyperplane, if
$\mathcal{E}_{Q}\subseteq\hat{\mathcal{E}}(R_{\zeta},\bar{\zeta})$ holds, then
$\mathcal{E}_{P}\subseteq\mathcal{E}(R_{\zeta},\bar{\zeta})$ holds.

Assumption 9 essentially means that the states of the secondary controller
$\zeta_{2}$ is not taken into account by the safe set. We believe this
assumption is reasonable, since safety of the states of the plant are often
the primary concern in most applications.

Since $P$ has the block structure given in (22), from Corollary 2 of Murguia
et al. (2020), it can be verified that $Q=Y-N^{\top}\tilde{Y}^{-1}N=X^{-1}$
using matrix inversion lemma. From Assumption 9, the constraint
$\mathcal{E}_{P}\subseteq\mathcal{E}(R_{\zeta},\bar{\zeta})$ can be
equivalently written as $\zeta_{1}^{\top}Q\zeta_{1}\leq
1\implies(\zeta_{1}-\bar{\zeta}_{1})^{\top}\hat{R}_{\zeta}(\zeta_{1}-\bar{\zeta}_{1})\leq
1$, where $\hat{R}_{\zeta}$ is such that

$\hat{\mathcal{E}}(R_{\zeta},\bar{\zeta})=\\{\zeta_{1}\ |\
(\zeta_{1}-\bar{\zeta}_{1})^{\top}\hat{R}_{\zeta}(\zeta_{1}-\bar{\zeta}_{1})\leq
1\\}.$ (33)

Since we have $Q=X^{-1}$, the condition $\zeta_{1}^{\top}Q\zeta_{1}\leq 1$ can
be equivalently stated as

$[\zeta_{1}^{\top}\ 1]\left[\begin{array}[]{cc}X^{-1}&\mathbf{0}\\\
\star&-1\end{array}\right][\zeta_{1}^{\top}\ 1]^{\top}\leq 0.$

Similarly, the condition $\zeta_{1}^{\top}\hat{R}_{\zeta}\zeta_{1}\leq 1$ can
be written as

$[\zeta_{1}^{\top}\
1]\left[\begin{array}[]{cc}\hat{R}_{\zeta}&-\hat{R}_{\zeta}\bar{\zeta}_{1}\\\
\star&\bar{\zeta}_{1}^{\top}\hat{R}_{\zeta}\bar{\zeta}_{1}-1\end{array}\right][\zeta_{1}^{\top}\
1]^{\top}\leq 0.$

Using $\mathcal{S}$-procedure again, $\zeta_{1}^{\top}Q\zeta_{1}\leq
1\implies(\zeta_{1}-\bar{\zeta}_{1})^{\top}\hat{R}_{\zeta}(\zeta_{1}-\bar{\zeta}_{1})\leq
1$ holds when there exists a non-negative constant $\delta$ such that the
following holds

$\left[\begin{array}[]{cc}\hat{R}_{\zeta}&-\hat{R}_{\zeta}\bar{\zeta}_{1}\\\
\star&\bar{\zeta}_{1}^{\top}\hat{R}_{\zeta}\bar{\zeta}_{1}-1\end{array}\right]-\delta\left[\begin{array}[]{cc}X^{-1}&\mathbf{0}\\\
\star&-1\end{array}\right]\preceq 0.$ (34)

Note that (34) is nonlinear in the variable $X$. Thus, we apply another
congruence transformation with $\text{diag}(X,I)$ to (34) leading to

$\left[\begin{array}[]{cc}X\hat{R}_{\zeta}X&-X\hat{R}_{\zeta}\bar{\zeta}_{1}\\\
\star&\bar{\zeta}_{1}^{\top}\hat{R}_{\zeta}\bar{\zeta}_{1}-1\end{array}\right]-\delta\left[\begin{array}[]{cc}X&\mathbf{0}\\\
\star&-1\end{array}\right]\preceq 0.$ (35)

Note that the first term of (35) can be written as

$\left[\begin{array}[]{cc}\mathbf{0}&-X\hat{R}_{\zeta}\bar{\zeta}_{1}\\\
\star&\bar{\zeta}_{1}^{\top}\hat{R}_{\zeta}\bar{\zeta}_{1}-1\end{array}\right]+[X\
\mathbf{0}]^{\top}\hat{R}_{\zeta}[X\ \mathbf{0}].$

Then, via the Schur complement result (see Boyd et al. (1994), for example),
inequality (35) is equivalent to

$\left[\begin{array}[]{ccc}\mathbf{0}&-X\hat{R}_{\zeta}\bar{\zeta}_{1}&-X\\\
\star&\bar{\zeta}_{1}^{\top}\hat{R}_{\zeta}\bar{\zeta}_{1}-1&\mathbf{0}\\\
\star&\star&-\hat{R}_{\zeta}^{-1}\end{array}\right]-\delta\left[\begin{array}[]{ccc}X&\mathbf{0}&\mathbf{0}\\\
\star&-1&\mathbf{0}\\\ \star&\star&\mathbf{0}\end{array}\right]\preceq 0,$
(36)

which is linear in $X$.

To summarize, the problem (21) can be equivalently formulated as

$\begin{split}\underset{\eta,\alpha,\beta,\delta}{\min}&\ \text{Tr}[R_{a}]\\\
s.t.&\ \mathbf{E}_{2}+\alpha\mathbf{F}+\beta\mathbf{S}\preceq 0\\\ &\
P(\eta)\succ 0\\\ &\ \mathbf{J}-\delta\mathbf{L}\preceq 0,\end{split}$ (37)

where

$\mathbf{J}=\left[\begin{array}[]{ccc}\mathbf{0}&-X\hat{R}_{\zeta}\bar{\zeta}_{1}&-X\\\
\star&\bar{\zeta}_{1}^{\top}\hat{R}_{\zeta}\bar{\zeta}_{1}-1&\mathbf{0}\\\
\star&\star&-\hat{R}_{\zeta}^{-1}\end{array}\right],\
\mathbf{L}=\left[\begin{array}[]{ccc}X&\mathbf{0}&\mathbf{0}\\\
\star&-1&\mathbf{0}\\\ \star&\star&\mathbf{0}\end{array}\right].$

Notice that (37) can be efficiently solved for fixed non-negative constants
$\alpha,\beta$ and $\delta$. Since $P(\eta)\succ 0$, we have $Y\succ 0$ and
$X-Y^{-1}\succ 0$. Moreover, by the definition of $\Pi_{1}$ and $\Pi_{2}$, we
have $MN^{\top}=I-XY$ guaranteed to be non-singular. Therefore, once (37) is
solved, appropriate $M$ and $N$ can always be found. Lastly, the secondary
controller variable $\kappa$ can be found via (27) in the order of $D_{2}$,
$C_{2}$, $B_{2}$, and $A_{2}$, see Scherer et al. (1997).

We summarize the main result presented above in the following theorem.

######  Theorem 10.

Suppose Assumption 9 holds and $n_{1}=n_{2}$. Given $E_{u}$ and $C_{S}$. If
there exist a positive definite matrix $R_{a}$, a positive definite matrix
$P(\eta)$, and non-negative constants $\alpha,\beta$, and $\delta$ that solve
the optimization problem (37), then for all attack signals satisfying (6),
$\mathcal{E}_{P}$ is a subset of the safe set
$\mathcal{E}(R_{\zeta},\bar{\zeta})$ and is a forward invariant set for (7).

When $R_{a}$ is given, we can give a synthesis result similar to Proposition
2.

######  Corollary 11.

Suppose Assumption 9 holds and $n_{1}=n_{2}$. Given $R_{a}$, $E_{u}$, and
$C_{S}$. If there exists a positive definite matrix $P(\eta)$, and non-
negative constants $\alpha,\beta$, and $\delta$ that satisfy

$\begin{split}&\ \mathbf{E}_{2}+\alpha\mathbf{F}+\beta\mathbf{S}\preceq 0\\\
&\ P(\eta)\succ 0\\\ &\ \mathbf{J}-\delta\mathbf{L}\preceq 0,\end{split}$

then for all attack signals satisfying (6), $\mathcal{E}_{P}$ is a subset of
the safe set $\mathcal{E}(R_{\zeta},\bar{\zeta})$ and is a forward invariant
set for (7).

Corollary 11 gives a sufficient condition to verify the possibility of
designing an output feedback linear controller to ensure safety of the closed
loop system. For example, if the optimization problem in Proposition 2 is
infeasible meaning the states of the system with primary controller only may
be driven outside the safe set. However, if the problem in Corollary 11 is
feasible under the same $R_{a}$, then it can be guaranteed that with the
addition of the secondary controller, the closed loop system can be rendered
safe.

##  5 Numerical Case Study

We now illustrate the main result of this work via numerical simulations. We
will consider a linear time-invariant plant with two inputs and two outputs
that is stabilized by a linear primary controller. However, under potential
attacks that satisfy (6), the safety of the overall system can not be verified
via Proposition 2. Given the same $R_{a}$, we attempt to ensure safety of the
system via Corollary 11 by designing a secondary output feedback linear
controller.

For the secondary controller, we use one secured sensor in local feedback and
modify one of the outputs. This makes $E_{u}=[1\ 0]^{\top}$ and $C_{S}=[1\
0]$. The linear plant is stabilized by the primary controller such that we
have

$\hat{A}=\left[\begin{array}[]{cc}-1&0\\\
0&-1\end{array}\right],\hat{B}=\left[\begin{array}[]{c}1\\\
0\end{array}\right],\hat{C}=[1\ 0].$

The matrix $R_{a}$ that characterizes the bound on the attack signals is given
by $R_{a}=\text{diag}(0.5,1,0.4,0.7)$. All inputs and outputs are subject to
potential attacks. Lastly, the safe set is given by $\\{\zeta\ |\
\zeta_{1}^{\top}\hat{R}_{\zeta}\zeta_{1}\leq 1\\}$ with
$\hat{R}_{\zeta}=0.022I_{2}$ which is a sphere centred around the origin.

Via Proposition 2, we can find an invariant set only if we drop the constraint
that the invariant set is a subset of the safe set. The plots of the found
invariant set and the safe set is shown in Figure 2. It can be seen that the
found invariant set which serves as an outer approximation of the reachable
set111It is mentioned in Murguia et al. (2020) that the approximation is tight
for linear systems. is not a subset of the safe set. This means the system may
be driven to an unsafe region by attack signals.

Figure 2: Plots of the safe set and the found invariant set using primary
controller only.

With the secondary controller design, we aim to guarantee that the invariant
ellipsoid is a subset of the safe set. If this is feasible, then there exists
at least one controller that keeps the closed loop system safe under all
attack signals that satisfy the constraint (6).

We choose set $\alpha=0.25$ and $\delta=0.99$ to make the LMI conditions in
Corollary 11. We solve these LMIs to obtain

$X=\left[\begin{array}[]{cc}28.6109&0\\\
0&31.9965\end{array}\right],Y=\left[\begin{array}[]{cc}3.9840&0\\\
0&0.1164\end{array}\right],$
$\mathbf{A}=\left[\begin{array}[]{cc}-26.8308&0\\\
0&-0.6765\end{array}\right],\ \mathbf{B}=\left[\begin{array}[]{c}-200.3492\\\
0\end{array}\right],$ $\mathbf{C}=[-78.5049\ 0],\ \mathbf{D}=-26.8308.$

We then choose $M=I_{2}$. The fact that $MN^{\top}=I-XY$ yields
$N=\left[\begin{array}[]{cc}-112.9854&0\\\ 0&-2.7259\end{array}\right]$. The
controller parameters can then be solved as follows in the order of $D_{2}$,
$C_{2}$, $B_{2}$, and lastly $A_{2}$ using (27).

$A_{2}=\left[\begin{array}[]{cc}-27.2049&0\\\ 0&-1.1187\end{array}\right],\
B_{2}=\left[\begin{array}[]{c}0.8271\\\ 0\end{array}\right],$
$C_{2}=[689.1488\ 0],\ D_{2}=-26.8308.$

Figure 3: Plots of the safe set and the found invariant set with the addition
of the secondary controller.

It can be seen from Figure 3 that, the found invariant ellipsoid is now a
subset of the safe set. Therefore, by adding the secondary controller to the
system, we rule out the possibility of unsafe operations of the system subject
to attacks satisfying (6). One can also try to minimize $\text{Tr}[X]$ under
the constraints given in Corollary 11. Recall that via Lemma 6, we prove that
$\text{Tr}[X]$ can be used as a convex upper bound on the volume of
$\mathcal{E}_{Q}$, with $Q$ being the projection of $P$ on to the
$\zeta_{1}$-hyperplane. The following controller is synthesized if we minimize
$\text{Tr}[X]$,

$A_{2}=\left[\begin{array}[]{cc}-1.1745\times 10^{8}&0\\\
0&-1.8725\end{array}\right],\ B_{2}=\left[\begin{array}[]{c}5.0553\times
10^{6}\\\ 0\end{array}\right],$ $C_{2}=[-977.4890\ 0],\ D_{2}=-232.8737.$

The plot of the invariant ellipsoid is shown in Figure. 4. It can be seen
that, the volume of the invariant set in Figure. 4 is significantly reduced
compared to the one shown in Figure. 3.

Figure 4: Plots of the safe set and the found invariant set with the addition
of the secondary controller that minimizes $\text{Tr}[X]$.

##  6 Conclusions

We have used an invariant set based method to provide a framework for checking
the safety of LTI systems subject to sensor and actuator attacks by resource
limited adversaries. In addition, by using a subset of sensors and secure
feedback, a plug-and-play linear secondary controller synthesis problem is
also solved to recover the safety of the overall system. The initial synthesis
problem is not convex, but can be rendered so with a congruence
transformation, which yields an LMI condition that can be solved efficiently.
The effectiveness of the proposed design is illustrated via a numerical
example, where a secondary controller is designed to guarantee safety which
can not be guaranteed by the primary controller alone.

One possible future research direction is the investigation of how to smarly
choose the matrices $C_{S}$ and $E_{u}$. In our work, we assume it is given,
i.e., the designer first chooses which sensors are to be secured and used for
secondary controller design and how the output of the secondary controller is
fed back to the plant. When the size of the system is large, it is desired to
have a more efficient and intelligent approach to choosing which set of
sensors are to be used for the secondary controller.

## References

  * Ames et al. (2017) Ames, A.D., Xu, X., Grizzle, J.W., and Tabuada, P. (2017).  Control barrier function based quadratic programs for safety critical systems.  _IEEE Transactions on Automatic Control_ , 62(8), 3861–3876.
  * Bai et al. (2015) Bai, C.Z., Pasqualetti, F., and Gupta, V. (2015).  Security in stochastic control systems: Fundamental limitations and performance bounds.  In _2015 American Control Conference (ACC)_ , 195–200. IEEE.
  * Blanchini (1999) Blanchini, F. (1999).  Set invariance in control.  _Automatica_ , 35(11), 1747–1767.
  * Boyd et al. (1994) Boyd, S., El Ghaoui, L., Feron, E., and Balakrishnan, V. (1994).  _Linear matrix inequalities in system and control theory_.  SIAM.
  * Cárdenas et al. (2008) Cárdenas, A.A., Amin, S., and Sastry, S. (2008).  Research challenges for the security of control systems.  _HotSec_ , 5, 15.
  * Chong et al. (2015) Chong, M.S., Wakaiki, M., and Hespanha, J.P. (2015).  Observability of linear systems under adversarial attacks.  In _2015 American Control Conference (ACC)_ , 2439–2444. IEEE.
  * Dutta and Langbort (2017) Dutta, A. and Langbort, C. (2017).  Stealthy output injection attacks on control systems with bounded variables.  _International Journal of Control_ , 90(7), 1389–1402.
  * Escudero et al. (2022a) Escudero, C., Massioni, P., Zamaï, E., and Raison, B. (2022a).  Analysis, prevention, and feasibility assessment of stealthy ageing attacks on dynamical systems.  _IET Control Theory & Applications_, 16(4), 381–397.
  * Escudero et al. (2022b) Escudero, C., Murguia, C., Massioni, P., and Zamaï, E. (2022b).  Enforcing safety under actuator injection attacks through input filtering.  In _2022 European Control Conference (ECC)_ , 1521–1528.
  * Fawzi et al. (2014) Fawzi, H., Tabuada, P., and Diggavi, S. (2014).  Secure estimation and control for cyber-physical systems under adversarial attacks.  _IEEE Transactions on Automatic control_ , 59(6), 1454–1467.
  * Kurzhanski and Vályi (1997) Kurzhanski, A. and Vályi, I. (1997).  _Ellipsoidal calculus for estimation and control_.  Springer.
  * Milošević et al. (2020) Milošević, J., Sandberg, H., and Johansson, K.H. (2020).  Estimating the impact of cyber-attack strategies for stochastic networked control systems.  _IEEE Transactions on Control of Network Systems_ , 7(2), 747–757.
  * Murguia et al. (2020) Murguia, C., Shames, I., Ruths, J., and Nešić, D. (2020).  Security metrics and synthesis of secure control systems.  _Automatica_ , 115, 108757.
  * Romagnoli et al. (2020) Romagnoli, R., Griffioen, P., Krogh, B.H., and Sinopoli, B. (2020).  Software rejuvenation under persistent attacks in constrained environments.  _IFAC-PapersOnLine_ , 53(2), 4088–4094.
  * Scherer et al. (1997) Scherer, C., Gahinet, P., and Chilali, M. (1997).  Multiobjective output-feedback control via lmi optimization.  _IEEE Transactions on automatic control_ , 42(7), 896–911.
  * Tang et al. (2019) Tang, Z., Kuijper, M., Chong, M.S., Mareels, I., and Leckie, C. (2019).  Linear system security—detection and correction of adversarial sensor attacks in the noise-free case.  _Automatica_ , 101, 53–59.
  * Teixeira et al. (2015a) Teixeira, A., Shames, I., Sandberg, H., and Johansson, K.H. (2015a).  A secure control framework for resource-limited adversaries.  _Automatica_ , 51, 135–148.
  * Teixeira et al. (2015b) Teixeira, A., Sou, K.C., Sandberg, H., and Johansson, K.H. (2015b).  Secure control systems: A quantitative risk management approach.  _IEEE Control Systems Magazine_ , 35(1), 24–45.