# Cyber-Physical Attack Leveraging Subsynchronous Resonance

Bosong Li, Baosen Zhang, Daniel S. Kirschen  Department of Electrical and
Computer Engineering
University of Washington

###### Abstract

This paper discusses how a cyber attack could take advantage of torsional
resonances in the shaft of turbo-generators to inflict severe physical damage
to a power system. If attackers were able to take over the control of a
battery energy storage device, they could modulate the injection of this
device at a frequency that matches one of the sub-synchronous resonance
frequencies of a generator. Small changes in injection might be sufficient to
excite one of these mechanical resonances, resulting in metal fatigue and
ultimately a catastrophic failure in the shaft of the generator. Using a
state-space model of the electro-mechanical system, the paper develops
transfer functions linking the magnitude of the malicious injections to the
magnitude of oscillations in the speed and angle of the various masses
connected to the shaft. Numerical results from a two-area power system
demonstrate the existence of vulnerable frequencies and show that damaging
mechanical oscillations can be triggered without causing easily detectable
signals at the generator terminals.

###### Index Terms:

Cyber attack, cyber-physical attack, battery energy storage, sub-synchronous
resonance, state space analysis

##  I Introduction

While a cyber attack [1, 2] [3, 4][5, 6] cripples a power system until the
malicious software has been expurgated, a physical attack that destroys a
major piece of equipment can take months to repair. To carry out a physical
attack, the adversary must get in close proximity to the target equipment,
which is expensive and difficult to carry out without detection. On the other
hand, a cyber-physical attack does not require the malicious actor to get in
close proximity of the target equipment. Instead, the attacker infiltrates
some aspect of the control infrastructure of the power system and manipulates
it to create physical damage [7]. Idaho National Laboratory demonstrated the
feasibility of this type of attack by taking over the protection system of a
generator and manipulating its synchronization until the generator self-
destructed [8].

To avoid detection and countermeasures, malicious manipulations of the control
system should remain small. It is therefore important to explore how an
attacker could use resonances to amplify their effects. Because their output
can be modulated at high frequency, battery energy storage systems (BESS)
represent an ideal vector for this type of attack. Furthermore, because the
number of BESS deployed in power systems is increasing rapidly, the attack
surface is growing [9] [10] [11].

This paper explores how a malicious actor could damage large turbo-generators
by taking over the control of a BESS and using it to inject small amounts of
power at frequencies corresponding to the torsional sub-synchronous resonance
(SSR) frequencies in the shaft of some generators. SSR is a condition of the
electric power system where with a turbine generator exchanges energy with the
rest of the system at one or more of the natural frequencies of the combined
system below the synchronous frequency of the system [12]. While sub-
synchronous oscillations between the various masses connected to the shaft of
the generator typically remain small, they cause metal fatigue and can over
time lead to a catastrophic failure of this shaft [13].

Physical and control countermeasures can be taken to avoid SSR [14]. Capacitor
compensation can be added to the system, but the investment cost of such
infrastructure measures is high. Protective relays can also be used to detect
oscillations at the generator terminals. However, as this paper will show, it
is possible to induce mechanical oscillations that are hard to detect at the
generator terminals.

To demonstrate the feasibility of cyber-physical attacks that leverage sub-
synchronous resonances, this paper develops a state-space model of the
combined electromechanical system. From this state space model, we derive
transfer functions linking the magnitude of the disturbances created by the
subverted BESS and the angular frequency of the various masses connected to
the shaft of a generator. These transfer functions exhibit resonance
frequencies that an attacker could target. A relatively low-power BESS could
therefore trigger sub-synchronous resonances that could ultimately destroy
large generators while remaining hard to detect.

The remainder of the paper is organized as follows: Section II describes a
state-space model of the system dynamics that combines the swing equations of
the generator, the mechanical characteristics of the shaft system, as well as
the power flow equations. Section III derives transfer functions relating the
BESS injections to the angular frequency and position of the various masses
connected to the generator shaft. Section IV describes numerical studies that
illustrate the frequency-domain analysis and correlate it with time-domain
simulations. Section V discusses possible countermeasures. Section VI
concludes and discusses further work.

##  II Model of System Dynamics

###  II-A Notation

In this paper, the rated angular velocity in electrical rad/s is denoted by
$\omega_{0}=2\pi f_{0}$, where frequency $f_{0}=60$ Hz. Assuming the number of
field poles $p_{f}=2$, we employ $\omega_{0m}$ as the rated angular velocity
in mechanical rad/s, and $\omega_{0m}=(2/p_{f})\omega_{0}=377$ rad/s.

Each generator contains a five-mass torsional system. This paper denotes the
speed and angle deviation of each rotor from the steady-state values
respectively with $\omega_{g},\omega_{s1},\omega_{s2},\omega_{s3},\omega_{s4}$
and $\theta_{g},\theta_{s1},\theta_{s2},\theta_{s3},\theta_{s4}$. The
subscription $gi$ represents the terminal of generator $i$, while the
subscriptions $s1-s4$ represent the other turbine sections. We further define
$\Delta\omega$ and $\Delta\theta$ as the speed and angle difference between
two adjacent masses connecting to the same shaft. A torsional system of
generator $i$ is presented in Fig. 1 with the variables illustrated, where the
shaft between two adjacent masses is denoted using a double subscript as
$12,23,34$ and $45$.

In the following formulation, the variables and parameters are employed with
their per-unit values if not specially mentioned.

Figure 1: Five-Mass Torsional System.

###  II-B State-Space Model

We investigate the dynamics of generators considering the swing equations of
the generator rotor together with the torsional characteristics of the shaft
system shown below in (1)-(2) and (3)-(11).

$\displaystyle\dot{\omega}_{g}=\frac{\omega_{0m}}{2H_{g}}~{}(P_{M}-P_{e}-D_{g}\cdot\omega_{g})$
(1) $\displaystyle\dot{\theta}_{g}=\omega_{g}$ (2)

$\displaystyle\dot{\omega}_{g}=\frac{\omega_{0m}}{2H_{g}}~{}\big{(}K_{12}(\theta_{s2}-\theta_{g})-D_{g}\cdot\omega_{g}-P_{e}\big{)}$
(3)
$\displaystyle\dot{\omega}_{s1}=\frac{\omega_{0m}}{2H_{s1}}~{}\big{(}K_{23}(\theta_{s3}-\theta_{s2})-K_{12}(\theta_{s2}-\theta_{g})$
(4)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}-D_{s1}\cdot\omega_{s1}+P_{s1}\big{)}$
$\displaystyle\dot{\omega}_{s2}=\frac{\omega_{0m}}{2H_{s2}}~{}\big{(}K_{23}(\theta_{s3}-\theta_{s2})-K_{12}(\theta_{s2}-\theta_{s1})$
(5)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}-D_{s2}\cdot\omega_{s2}+P_{s2}\big{)}$
$\displaystyle\dot{\omega}_{s3}=\frac{\omega_{0m}}{2H_{s3}}~{}\big{(}K_{34}(\theta_{s4}-\theta_{s3})-K_{23}(\theta_{s3}-\theta_{s2})$
(6)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}-D_{s3}\cdot\omega_{s3}+P_{s3}\big{)}$
$\displaystyle\dot{\omega}_{s4}=\frac{\omega_{0m}}{2H_{s4}}~{}\big{(}-K_{34}(\theta_{s4}-\theta_{s3})-D_{s4}\cdot\omega_{s4}-P_{s4}\big{)}$
(7) $\displaystyle\dot{\theta}_{s1}=\omega_{s1}$ (8)
$\displaystyle\dot{\theta}_{s2}=\omega_{s2}$ (9)
$\displaystyle\dot{\theta}_{s3}=\omega_{s3}$ (10)
$\displaystyle\dot{\theta}_{s4}=\omega_{s4}$ (11)

In (3)-(7), parameters $D$ and $H$ are respectively the damping coefficient
and inertia constant of each rotor section. Parameter $K$ is the shaft
stiffness. Note that (1) and (3) illustrate the dynamics of the generator
rotor in different forms. Furthermore, the mechanical torque on the generator
rotor equals the electrical torque in the steady state, i.e. $P_{M}=P_{e}$. We
can therefore combine (1) and (3) as (12).

$\dot{\omega}_{g}=\frac{\omega_{0m}}{2H_{g}}~{}\big{(}K_{12}(\theta_{s2}-\theta_{g})-D_{g}\cdot\omega_{g}-P_{M}\big{)}$
(12)

Hence the dynamics of a generator considering the shaft system is fully
described by (2), (4)-(11), and (12). In a power system containing $n$
generators, we take the rotor speed $\boldsymbol{\omega}$ and angular
displacement $\boldsymbol{\theta}$ of each generator as state variable
$\mathbf{x}=[\boldsymbol{\omega}^{T},\boldsymbol{\theta}^{T}]^{T}$, where
$\boldsymbol{\omega}$ and $\boldsymbol{\theta}$ are both vectors of length
$5n$ illustrated in details as:

$\displaystyle\boldsymbol{\omega}=[\underbrace{\omega_{g1},\omega_{g2},...,\omega_{gn}}_{n},\underbrace{\omega_{g1s1},...,\omega_{g1s4},...,\omega_{gns1},...,\omega_{gns4}}_{4n}]^{T}$
(13)
$\displaystyle\boldsymbol{\theta}=[\underbrace{\theta_{g1},\theta_{g2},...,\theta_{gn}}_{n},\underbrace{\theta_{g1s1},...,\theta_{g1s4},...,\theta_{gns1},...,\theta_{gns4}}_{4n}]^{T}$
(14)

In the steady state, the mechanical power $P_{M}$ of the generator rotor in
(12) equals the electrical power. By substituting the load bus angular
displacement in the power flow and power balance equations with state
variables, we write the mechanical power using the angular displacement vector
$\boldsymbol{\theta}$ and load vector $\mathbf{L}$:

$\mathbf{P_{M}}=\mathbf{P_{e}}=\mathbf{A_{e}}\boldsymbol{\theta}+\mathbf{B_{e}}\mathbf{L}$
(15)

where $\mathbf{P_{M}}=[P_{Mg1},P_{Mg2},...,P_{Mgn}]^{T}$. The matrices
$\mathbf{A_{e}}$ and $\mathbf{B_{e}}$ depend on the system admittance matrix
and topology. Since the focus of this paper is on the oscillations between
different masses in the generator, we take a load bus to be the slack bus to
emphasis generators internal dynamics.

This paper further denotes the input power on each mass of the generator shaft
system with a 4-dimensional vector
$\mathbf{P_{gi}}=[P_{gis1},P_{gis2},P_{gis3},P_{gis4}]^{T}=P_{Mgi}\cdot\mathbf{B_{f}}_{gi}$,
where the coefficient vector $\mathbf{B_{f}}_{gi}$ contains the fraction of
the total turbine power generated by each turbine in the steady-state.
Following the pattern of state variables $\boldsymbol{\omega}$ and
$\boldsymbol{\theta}$, the input power vector
$\mathbf{P_{I}}=[\mathbf{P_{M}}^{T},\mathbf{P_{g1}}^{T},...,\mathbf{P_{gn}^{T}}]^{T}$
is defined as:

$\mathbf{P_{I}}=\mathbf{B_{I}}\cdot\mathbf{P_{M}}=\mathbf{B_{I}}\mathbf{A_{e}}\boldsymbol{\theta}+\mathbf{B_{I}}\mathbf{B_{e}}\mathbf{L}$
(16) $\mathbf{B_{I}}=\left[\begin{smallmatrix}\mathbf{I_{n\times n}}\\\
\mathbf{B_{F}}\end{smallmatrix}\right]$ (17)
$\mathbf{B_{F}}=\left[\begin{smallmatrix}\mathbf{B_{f}}_{g1}&&&\\\
&\mathbf{B_{f}}_{g2}&&\\\ &&\ddots&\\\
&&&\mathbf{B_{f}}_{gn}\end{smallmatrix}\right]$ (18)

where $\mathbf{I_{n\times n}}$ is an $n$ by $n$ identity matrix. The $4n\times
n$-dimensional coefficient matrix $\mathbf{B_{F}}$ is constructed by matrices
$\mathbf{B_{f}}_{g1},\mathbf{B_{f}}_{g2},...,\mathbf{B_{f}}_{gn}$ indicating
how much input power each turbine provides in a generator.

We therefore build the state-space model of the whole power system as follows:

$\left[\begin{array}[]{c}\dot{\boldsymbol{\omega}}\\\
\dot{\boldsymbol{\theta}}\end{array}\right]=\left[\begin{array}[]{cc}\mathbf{A_{11}}&\mathbf{A_{12}}\\\
\mathbf{A_{21}}&\mathbf{A_{22}}\end{array}\right]\left[\begin{array}[]{c}\boldsymbol{\omega}\\\
\boldsymbol{\theta}\end{array}\right]+\mathbf{B}\cdot\mathbf{L}$ (19)

The system matrix $\boldsymbol{A}$ consists of four parts, where:

$\displaystyle\mathbf{A_{11}}\\!=\\!-\frac{\omega_{0m}}{2}\mathbf{diag}{(-\frac{D_{g}}{H_{g}},\\!-\frac{D_{s1}}{H_{s1}},\\!-\frac{D_{s2}}{H_{s2}},\\!-\frac{D_{s3}}{H_{s3}},\\!-\frac{D_{s4}}{H_{s4}})}$
(20)
$\displaystyle\mathbf{A_{12}}\\!\\!=\\!\\!\mathbf{B_{I}}\\!\mathbf{A_{e}}\\!\\!+\\!\\!\frac{\omega_{0m}}{2}\\!\\!\\!\left[{\small\begin{matrix}\begin{smallmatrix}-\frac{K_{12}}{H_{g}}&\frac{K_{12}}{H_{g}}&&&\\\
\frac{K_{12}}{H_{s2}}&-\frac{K_{12}+K_{K23}}{H_{s2}}&\frac{K_{23}}{H_{s2}}&&\\\
&\frac{K_{23}}{H_{s3}}&-\frac{K_{23}+K_{34}}{H_{s3}}&\frac{K_{34}}{H_{s3}}&\\\
&&\frac{K_{34}}{H_{s4}}&-\frac{K_{34}+K_{45}}{H_{s4}}&\frac{K_{45}}{H_{s4}}\\\
&&&\frac{K_{45}}{H_{s5}}&-\frac{K_{45}}{H_{s5}}\\\
\end{smallmatrix}\end{matrix}}\right]$ (21)

The sub-matrix $\mathbf{A_{21}}\\!=\\!\mathbf{I_{5n\times 5n}}$ is an identity
matrix, while sub-matrix $\mathbf{A_{22}}\\!=\\!\mathbf{0}$ is a zero matrix.
The input matrix $\mathbf{B}=\mathbf{B_{I}}\mathbf{B_{e}}$ according to (16).

##  III Transfer Function Magnitude Analysis

With access to the BESS on a load bus, potential attackers can inject
oscillatory signals into the system, illustrated mathematically as (22), where
the attack signal $\Delta L$ is a square-wave signal.  This paper assumes a
square-wave attack signal because it is created by an energy storage device by
constantly changing between the charging and discharging states. Moreover, we
built a model predictive control (MPC) problem to maximize the oscillation of
the generator terminal variables, and the optimal solution to this problem was
a square wave.

$\dot{\mathbf{x}}=\left[\begin{array}[]{cc}\mathbf{A_{11}}&\mathbf{A_{12}}\\\
\mathbf{A_{21}}&\mathbf{A_{22}}\end{array}\right]\mathbf{x}+\mathbf{B}(\mathbf{L}+\Delta\mathbf{L})$
(22)

All the elements of the attack signal $\Delta\mathbf{L}$ are 0 except for the
element corresponding to the load bus where the malicious energy storage
device resides. We assume the oscillatory signal magnitude is 1 p.u. of the
base power considering the limited power and capacity of the attack energy
storage device. The only factor that affects the system oscillation is,
therefore, the frequency of input signal $\Delta\mathbf{L}$. To explore the
impact of the injected oscillatory signal frequency on each rotor of the
generators in the system, we investigate two output signals denoted as (25)
and (28). To explore the impact of the oscillatory signal injection on the
generator terminal and the torsional system, we define the following output:

$\displaystyle\mathbf{y_{1}}=\left[\begin{array}[]{c}\boldsymbol{\omega}\\\
\boldsymbol{\theta}\end{array}\right]=\mathbf{C_{1}}\mathbf{x}+\mathbf{D_{1}}(\mathbf{L}+\Delta\mathbf{L})$
(25)
$\displaystyle\mathbf{y_{2}}=\left[\begin{array}[]{c}\Delta\boldsymbol{\omega}\\\
\Delta\boldsymbol{\theta}\end{array}\right]=\mathbf{C_{2}}\mathbf{x}+\mathbf{D_{2}}(\mathbf{L}+\Delta\mathbf{L})$
(28)

where output $\mathbf{y_{1}}$ is the generator terminal variables. The output
$\mathbf{y_{2}}$ contains torsional variables showing the speed and angle
differences between the rotors connecting to the same shaft. Output
$\mathbf{y_{2}}$ is further explained in (29)-(32).

$\displaystyle\Delta\boldsymbol{\omega}=[\Delta\boldsymbol{\omega}_{g1},\Delta\boldsymbol{\omega}_{g2},...,\Delta\boldsymbol{\omega}_{gn}]^{T}$
(29)
$\displaystyle\Delta\boldsymbol{\theta}=[\Delta\boldsymbol{\theta}_{g1},\Delta\boldsymbol{\theta}_{g2},...,\Delta\boldsymbol{\theta}_{gn}]^{T}$
(30)

For generator $i$, $\Delta\boldsymbol{\omega}_{gi}$ and
$\Delta\boldsymbol{\theta}_{gi}$ are both 4-dimensional vectors:

$\displaystyle\Delta\boldsymbol{\omega_{gi}}\\!\\!=\\!\\![\omega_{gi}\\!\\!-\\!\\!\omega_{gis1},\omega_{gis1}\\!\\!-\\!\\!\omega_{gis2},\omega_{gis2}\\!\\!-\\!\\!\omega_{gis3},\omega_{gis3}\\!\\!-\\!\\!\omega_{gis4}]^{T}$
(31)
$\displaystyle\Delta\boldsymbol{\theta_{gi}}\\!\\!=\\!\\![\theta_{gi}\\!\\!-\\!\\!\theta_{gis1},\theta_{gis1}\\!\\!-\\!\\!\theta_{gis2},\theta_{gis2}\\!\\!-\\!\\!\theta_{gis3},\theta_{gis3}\\!\\!-\\!\\!\theta_{gis4}]^{T}$
(32)

We analyze the magnitudes of the transfer functions between the input signal
$\Delta\mathbf{L}$ and the output signals in (25) and (28) over the frequency
spectrum from 0Hz to 60Hz. The magnitudes of the transfer functions
$|\Gamma_{1}|_{j}$ and $|\Gamma_{2}|_{k}$ correspond respectively to the $jth$
generator terminal variable and the $kth$ torsional variable of a generator.

$\displaystyle|y_{1}(j\omega_{a})|_{j}=|\Gamma_{1}(j\omega_{a})|_{j}\cdot|\Delta
L(j\omega_{a})|$ (33)
$\displaystyle|y_{2}(j\omega_{a})|_{k}=|\Gamma_{2}(j\omega_{a})|_{k}\cdot|\Delta
L(j\omega_{a})|$ (34)

Equations (33) and (34) show that with the same oscillatory input signal, the
oscillation magnitudes of the generator terminal variables as well as the
torsional variables are determined by the magnitude of the corresponding
transfer function. With an oscillatory input signal at frequency $\omega_{a}$,
if the magnitude $|\Gamma_{2}(j\omega_{a})|_{k}$ is higher than
$|\Gamma_{1}(j\omega_{a})|_{j}$ as shown in (35), the internal torsional
system of the generator suffers a more severe oscillation compared to the
generator terminal.

$|\Gamma_{2}(j\omega_{a})|_{k}>|\Gamma_{1}(j\omega_{a})|_{j},j,k\in
the~{}same~{}generator$ (35)

Additionally, the protection system is usually designed to prevent the
generator terminal rotor from significant oscillation, i.e., damping the
frequencies that lead to high $|\Gamma_{1}(j\omega_{a})|_{i}$ values. The
input oscillation frequencies denoted in (35) therefore remain commonly
neglected. A ratio parameter $R_{M}$ is further defined in (36) as the ratio
between the transfer function magnitudes of the torsional output and the
terminal output.

$R_{M}(\omega_{a})=\frac{|\Gamma_{2}(j\omega_{a})|_{k}}{|\Gamma_{1}(j\omega_{a})|_{i}},i,k\in
the~{}same~{}generator$ (36)

The higher the ratio $R_{M}$ is, the more severe the oscillation inside the
generator shaft system occurs while the measurements at the generator terminal
remain close to the steady-state values. If the attack frequency $\omega_{a}$
in (36) is lower than the nominal frequency, i.e. $\omega_{a}<\omega_{0}$, the
input attack signal leads to subsynchronous oscillation of the generator. Such
oscillation is easily neglected in the power system, for the generator
terminal measurements stay nearly unchanged. However, the subsynchronous
oscillation will cause mechanical fatigue and fracture in the long term. In
addition to the theoretical analysis in the frequency domain, we also testify
the vulnerability leveraging SSR in the time domain, shown with the numerical
studies in the following section.

##  IV A Two-area System Example

To investigate the existence of the vulnerability presented in this paper in
the power system, we employ a commonly used two-area system for numerical
studies shown in Fig. 2.

Figure 2: Two-area System

In Fig. 2, the line reactances are shown in p.u. on 230 kW and 100 MW [14]. A
linear system is considered here with a single circuit tie line. And an energy
storage device is integrated to Bus 3 to inject malicious input signal with a
magnitude of 1 MW, which is 1 p.u.. The power transfer from Area 1 to Area 2
is 400 MW, where L1 = 970 MW and L2 = 1770 MW. Four identical generating units
Gen 1, Gen 2, Gen 11, and G 12, are respectively loaded to 700 MW, 670 MW, 670
MW, and 700 MW, with a five-mass torsional system each. The inertia constant,
power fraction, and shaft stiffness of the generator torsional system are
presented in Table I.

TABLE I: Generator Torsional System parameters

###  IV-A Frequency-Domain Analysis

While the damping factor of the generator terminal is commonly neglected, i.e.
$D_{g1}\\!=\\!D_{g2}\\!\\!=\\!D_{g11}\\!\\!=\\!D_{g12}\\!=\\!0$, the damping
factor of each mass in the torsional system strongly affects the transfer
function magnitudes according to our tests. The impact of the mass damping
factor on the transfer function magnitude is depicted in Fig. 3. Note that
such impact on each generator is similar, we take Gen 1 as an example to avoid
repeating.

Figure 3: Transfer Function Magnitude Curves of Gen 1.

In Fig. 3, the transfer function magnitude curves of rotor speed $\omega$ and
angular displacement $\theta$ are respectively normalized. Thus the y-axis
ranges from 0 to 1. The $|\Gamma|_{\theta}$ curves in Fig. 3 see their
maximums with the input signal frequencies close to 0Hz. The maximums of the
$|\Gamma|_{\theta}$ curves are significantly higher than the other
$|\Gamma|_{\theta}$ values, which indicates the most notable oscillation
occurs when the input signals approximate DC signals. With the increase of
mass damping factor, the $|\Gamma|_{\omega}$ curve becomes smoother and moves
upward. A more significant oscillation is therefore expected due to the higher
values of the magnitude $|\Gamma|_{\omega}$. Hence it is reasonable for us to
employ a low damping factor for the torsional mass in further analysis. A low
mass damping factor is also consistent with the commonly-used steel masses in
real power systems.

The values of the ratio $R_{M}(\omega_{a})$ proposed in (36) are shown in Fig.
4 for Gen 1:

Figure 4: $R_{M}$ Values of Gen 1

where Fig. 4 (a) presents the magnitude ratio between the torsional speed
difference and the terminal rotor speed, and Fig. 4 (b) illustrates the
magnitude ratio between the torsional angle difference and the terminal rotor
angel. The higher $R_{M}(\omega_{a})$ value is, the more severe the
oscillation of the torsional system will be compared to the generator
terminal. Due to the low magnitudes of their transfer functions, the generator
terminal variables stay closely to their steady-state values. The internal
oscillation of the torsional system will thus be difficult to recognize. The
system is also more vulnerable when exposed to the malicious input at
frequency $(\omega_{a})$.

In Fig. 5, the blue curves are the terminal speed and angle of Gen 1. The
orange curves represent the speed and angel differences between the third and
fourth mass, both connecting to shaft 34 of Gen 1.

Figure 5: Comparison of $|\Gamma|_{\omega}$ , $|\Gamma|_{\theta}$ and
$|\Gamma|_{\Delta\omega}$ , $|\Gamma|_{\Delta\theta}$ of Gen 1

As shown in Fig. 5(a), the $|\Gamma|_{\theta g1}$ values are close to 0 when
the input signal frequency is above 0Hz, which explains the reason why
$R_{M\theta}$ values in Fig. 4(b) are significantly high. Therefore the the
local maximums of $R_{M\theta}$ shown in 4(b) coincide with the local maximums
of $|\Gamma|_{\Delta\theta}$ shown in 5(a). In Fig. 5(b), when the input
signal frequency is close to 34.47Hz, both curves see a local maximum value.
This local maximum is also the global maximum of $|\Gamma|_{\theta g1}$,
meaning the oscillation at this frequency is usually damped by power system
stabilizers, which are usually designed according to the speed and angle of
the generator terminals. However, the oscillations corresponding to other
local maximums often remain neglected from protective measures due to the lack
of attention to the speed and angle of the torsional masses. The torsional
system is therefore vulnerable to attack signals at these frequencies.
Potential attacks can thus attack the system with injection of oscillatory
signals at such frequencies with an energy storage device either remotely
hacked or physically controlled by them. Considering the results shown in Fig.
3, the example system is at risk when input signal frequency lies within the
range of 25.42-25.47Hz or 26.73-26.81Hz.

Note that the input signal frequency ranges depicted in Fig. 3-Fig. 5 are
0-50Hz, which is the typical frequency range of SSR. From 50 to 60Hz, the
y-axis values in in Fig. 3-Fig. 5 approximate 0, and are left out for
concision.

###  IV-B Time-Domain Analysis

Assuming the system is exposed to a square-wave input signal with magnitude
$|\Delta L|=1MW$ at Bus 3, we analyze the time-domain response of the output
considering two attack frequencies: $f_{a1}=25.42Hz$, and $f_{a2}=26.81Hz$.
The severity of an oscillation is measured by the amplitude of a signal’s
deviation from its steady-state value. We define $R_{\omega
j}=\frac{max(|\Delta\omega_{j}-\Delta\omega_{j,0}|)}{max(|\omega_{i}-\Delta\omega_{i,0}|)}$,
where the numerator is the deviation of the speed difference between mass $j$
and mass $j+1$ from the initial steady-state difference. The denominator is
the deviation of the terminal rotor speed from its initial state. Thus if the
value of $R_{\omega j}$ is higher than 1, the oscillation inside the torsional
system is more significant than that at the terminal rotor, and the
significance increases with the rise of $R_{\omega j}$ value. The same
definition and characteristic are applied to $R_{\theta j}$, thus is not
repeated. Considering a 10s control horizon with a time step of $10^{-3}s$, an
attack is issued at $t=2s$, the values of $R_{\omega j}$ and $R_{\theta j}$
are presented in the following table.

TABLE II: $R_{\omega}$ and $R_{\theta}$ Values under Different Attack
Frequencies

The maximums of the $R_{\omega j}$ and $R_{\theta j}$ under each attack signal
frequency are highlighted in bold in Table II. The time-domain responses of
unit Gen 2 is further illustrated in Fig. 6, where Fig. 6(a) depicts output
$\Delta\omega_{g2s4}$ and $\omega_{g2}$, and Fig. 6(b) presents output
$\Delta\theta_{g2s4}$ and $\theta_{g2}$.

Figure 6: Time-domain Response under Attack Frequency 26.81Hz

The $R_{\omega j}$ and $R_{\theta j}$ values in Table II higher than 1
indicate that the torsional system of the corresponding generators oscillate
notably while the terminal measurements close to the initial steady-state
values. Considering attack signal frequencies $f_{a1}$ and $f_{a2}$ are both
below 60Hz, the vulnerability revealed here results in SSR of the generators.
And we prove the existence of such vulnerability that leads to SSR when the
system is exposed to the proposed cyber-physical attacks. This paper assumes
access to a complete knowledge of system structure and data, which is
difficult to achieve in practice. To explore the system vulnerability with
limited access to the system information remains our future work.

##  V Countermeasures

The cyber-physical attack discussed in this paper requires only a relatively
low-power BESS that can be located at some distance from the targeted
generators. The resulting SSR leads to metal fatigue of the generator shafts
through small oscillations and that could ultimately cause a fatal mechanical
failure. Further work is needed to identify generators that might be
particularly susceptible to this type of attack and to develop effective
countermeasures. These countermeasures can be divided into four categories:
prevention, detection, reaction and mitigation.

Prevention is the first line of defense. Using best cyber security practices,
it aims to prevent attackers from taking control of power devices.
Unfortunately, the attack surface is getting larger due to the increasing
number of BESS and other controllable components.

Detecting this form of cyber-physical attack and locating their source are
challenging issues. Because the malicious power injections are small, they are
likely to be lost in the noise in conventional SCADA measurements, especially
if the network is large and complex. More sensitive measurements at a higher
time resolution are likely to be needed to be able to pinpoint the source of
the attack, which may be in a remote part of the network.

Once an attack has been detected and its source identified, to ensure a quick
reaction, procedures must be in place to disconnect it from the system.
Alternatively, the targeted generators may need to be disconnected.

If prevention, detection and reaction are deemed insufficiently effective,
mitigation may be required. Such measures could be similar to those that were
developed to deal with naturally occurring sub-synchronous resonance.
Reference [15] groups these countermeasures into four categories: system
switching and generator tripping, generator and system modifications, relaying
and detecting devices, and filtering and damping devices. While some of these
measures have been shown to be effective, they may be harder to implement in
an adversarial context because the attack may be aimed at any generator rather
than at a particular generator identified as being susceptible to sub-
synchronous resonance through careful system studies. Protecting all
generators against attacks at any dangerous frequency may also be difficult or
very costly

##  VI Conclusion

This paper exposes a potential vulnerability of power systems to a cyber-
physical attack where a malicious actor could trigger a SSR in the shaft of
large generators by creating small oscillations in the active power injections
of a battery energy storage system. These small mechanical oscillations
between the various masses connected to the shaft would cause metal fatigue
and ultimately lead to a catastrophic failure of the generator. Because this
resonance takes place within the shaft, it does not create significant
perturbations at the terminal of the generator and might therefore be
difficult to detect. Gaining access to the control system of a relatively low
power BESS would therefore give an attacker an opportunity to inflict physical
damage on equipment of a considerably larger rating and importance to the
reliable operation of the system. Our further work will explore detection,
reaction and mitigation countermeasures.

## Acknowledgment

The work described in this paper was carried out with funding from the US
National Science Foundation under its CRISP - Critical Resilient
Interdependent Infrastructure Systems and Processes program, grant number
1832287.

## References

  * [1] J. Liu, Y. Gu, L. Zha, Y. Liu, and J. Cao, “Event-triggered $H_{\infty}$ load frequency control for multiarea power systems under hybrid cyber attacks,” IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 49, no. 8, pp. 1665–1678, 2019.
  * [2] W. Chen, D. Ding, H. Dong, and G. Wei, “Distributed resilient filtering for power systems subject to denial-of-service attacks,” IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 49, no. 8, pp. 1688–1697, 2019.
  * [3] Q. Wang, W. Tai, Y. Tang, and M. Ni, “Review of the false data injection attack against the cyber-physical power system,” IET Cyber-Physical Systems: Theory & Applications, vol. 4, no. 2, pp. 101–107, 2019.
  * [4] J. Zhao, L. Mili, and M.Wang, “A generalized false data injection attacks against power system nonlinear state estimator and countermeasures,” IEEE Transactions on Power Systems, vol. 33, no. 5, pp. 4868–4877, 2018\.
  * [5] L. Zheng, T. Gao, and X. Zhang, “Security protection and testing system for cyber-physical based smart power grid,” in Proceedings of PURPLE MOUNTAIN FORUM 2019-international forum on smart grid protection and control. Springer, 2020, pp. 847–857.
  * [6] S. Mousavian, M. Erol-Kantarci, L. Wu, and T. Ortmeyer, “A risk-based optimization model for electric vehicle infrastructure response to cyber attacks,” IEEE Transactions on Smart Grid, vol. 9, no. 6, pp. 6160–6169, 2017\.
  * [7] T. Nguyen, S. Wang, M. Alhazmi, M. Nazemi, A. Estebsari, and P. Dehghanian, “Electric power grid resilience to cyber adversaries: State of the art,” IEEE Access, vol. 8, pp. 87 592–87 608, 2020.
  * [8] M. Zeller, “Myth or reality—does the aurora vulnerability pose a risk to my generator?” in 2011 64th Annual Conference for Protective Relay Engineers. IEEE, 2011, pp. 130–136.
  * [9] G. Cauley, and M. Lauby, “High-impact low-frequency event risk to the north american bulk power system,” North American Electric Reliability Corporation (NERC), Atlanta, GA, Tech. Rep, 2010.
  * [10] D. A. Tejada-Arango, A. S. Siddiqui, S. Wogrin, and E. Centeno, “A review of energy storage system legislation in the us and the european union,” Current Sustainable/Renewable Energy Reports, vol. 6, no. 1, pp. 22–28, 2019.
  * [11] C. Glenn, D. Sterbentz, and A. Wright, “Cyber threat and vulnerability analysis of the us electric sector,” Idaho National Lab.(INL), Idaho Falls, ID (United States), Tech. Rep., 2016.
  * [12] IEEE SSR Working Group, “Proposed terms and definitions for subsynchronous resonance,” in IEEE Symposium on Countermeasures for Subsynchronous Resonance, IEEE Pub, 81TH0086-9-PWR, 1981, pp. p92–97.
  * [13] P. M. Anderson, B. L. Agrawal, and J. E. Van Ness, Subsynchronous resonance in power systems. John Wiley & Sons, 1999, vol. 9.
  * [14] P. Kundur, “Power system stability,” Power system stability and control,vol. 10, 2007.
  * [15] IEEE Subsynchronous Resonance Working Group, “Countermeasures to subsynchronous resonance problems,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-99, no. 5, pp. 1810–1818, 1980.


EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-EP-2016-047 LHCb-PAPER-2016-002 June 6, 2016

Observation of the ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{\mathchar
28931\relax}\phi$ decay

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

The ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{\mathchar 28931\relax}\phi$
decay is observed using data corresponding to an integrated luminosity of
3.0$\mbox{\,fb}^{-1}$ recorded by the LHCb experiment. The decay proceeds at
leading order via a ${b}\rightarrow{{s}{\overline{s}}}{s}$ loop transition and
is therefore sensitive to the possible presence of particles beyond the
Standard Model. A first observation is reported with a significance of $5.9$
standard deviations. The value of the branching fraction is measured to be
$(5.18\pm 1.04\pm 0.35\,^{+0.67}_{-0.62})\times 10^{-6}$, where the first
uncertainty is statistical, the second is systematic, and the third is related
to external inputs. Triple-product asymmetries are measured to be consistent
with zero.

Published in Phys. Lett. B759 (2016) 282

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

##  1 Introduction

In the Standard Model (SM), the flavour-changing neutral current decay
${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{\mathchar 28931\relax}\phi$
proceeds via a ${b}\rightarrow{{s}{\overline{s}}}{s}$ loop (penguin) process.
A Feynman diagram of the gluonic penguin that contributes to this decay at
leading order is displayed in Fig. 1. This transition has been the subject of
theoretical and experimental interest in ${B}^{0}_{s}$ and ${B}^{0}$ decays,
since possible beyond the SM particles in the loop could induce non-SM $C\\!P$
violation [1, 2, 3]. The process has been probed with decay-time-dependent
methods in the ${{B}^{0}_{s}}\rightarrow\phi\phi$ and
${{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi}$ decay
modes [4, 5, 6, 7], which test for $C\\!P$ violation in the interference
between mixing and decay. In addition, measurements of $C\\!P$ violation in
the decay have been performed with the flavour-specific
${{B}^{0}}\\!\rightarrow{{K}^{*0}}\phi$ channel [8]. The results to date are
consistent with $C\\!P$ conservation in the
${b}\rightarrow{{s}{\overline{s}}}{s}$ process. Model-independently, non-SM
physics contributions could appear differently in these decay modes, though
many models contain strong correlations [9].

Figure 1: Feynman diagram contributing to the ${{\mathchar
28931\relax}^{0}_{b}}\rightarrow{\mathchar 28931\relax}\phi$ decay.

Measurements with ${\mathchar 28931\relax}^{0}_{b}$ baryons offer the
possibility to look for $C\\!P$ violation in the decay, both by studying
$C\\!P$ asymmetries and by means of $T$-odd observables. These observables
have been studied in greater detail for ${B}^{0}_{s}$ and ${B}^{0}$ meson
decays than those for ${\mathchar 28931\relax}^{0}_{b}$ baryons [10, 11, 4,
8]. Proposed methods to study $T$-odd asymmetries of ${\mathchar
28931\relax}^{0}_{b}$ baryons [12] exploit the polarisation structure of
${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{\mathchar 28931\relax}V$ decays,
where $V$ denotes a vector resonance [12], and can be affected by the initial
${\mathchar 28931\relax}^{0}_{b}$ polarisation if non-zero. An LHCb
measurement of the initial polarisation in ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}{\mathchar 28931\relax}$ decays has yielded a value consistent with
zero, though polarisation at the level of $10\%$ is possible given statistical
uncertainties [13]. No SM prediction exists specifically for the $T$-odd
asymmetries in ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ decays, though no large asymmetries are expected given the
prediction of $C\\!P$ conservation in the decays of beauty mesons for the same
transition. Measurements of $C\\!P$ asymmetries have been performed by LHCb in
an inclusive analysis of ${{\mathchar
28931\relax}^{0}_{b}}\rightarrow{\mathchar 28931\relax}hh^{\prime}$ decays
[14], where $h(h^{\prime})$ refers to a kaon or pion, with corresponding
$C\\!P$ asymmetries measured to be consistent with zero.

In this paper, a measurement of the ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ branching
fraction is presented using the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ decay
as a normalisation channel, which has a measured branching fraction of
$(7.3^{+0.7}_{-0.6})\times 10^{-6}$ [15]. The selection requirements used to
isolate the ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ decay with well-understood efficiencies reject suitable
control channels for a $\Delta A_{{C\\!P}}$ measurement. The ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ sample is then
used to perform measurements of the $T$-odd triple-product asymmetries, which
do not require a control channel. The results are based on $pp$ collision data
corresponding to an integrated luminosity of $1.0\mbox{\,fb}^{-1}$ and
$2.0\mbox{\,fb}^{-1}$ collected by the LHCb experiment at centre-of-mass
energies of $\sqrt{s}=7\mathrm{\,Te\kern-1.00006ptV}$ in 2011 and
8$\mathrm{\,Te\kern-1.00006ptV}$ in 2012, respectively.

##  2 Detector and simulation

The LHCb detector [16, 17] is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$, designed for the study of particles
containing $b$ or $c$ quarks. The detector includes a high-precision tracking
system consisting of a silicon-strip vertex detector surrounding the $pp$
interaction region, a large-area silicon-strip detector located upstream of a
dipole magnet with a bending power of about $4{\mathrm{\,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}$. The minimum
distance of a track to a primary vertex, the impact parameter, is measured
with a resolution of $(15+29/\mbox{$p_{\mathrm{T}}$}){\,\upmu\mathrm{m}}$,
where $p_{\mathrm{T}}$ is the component of the momentum transverse to the
beam, in ${\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$. Different types of charged
hadrons are distinguished using information from two ring-imaging Cherenkov
detectors. Photons, electrons and hadrons are identified by a calorimeter
system consisting of scintillating-pad and preshower detectors, an
electromagnetic calorimeter and a hadronic calorimeter. The online event
selection is performed by a trigger, which 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. At the hardware trigger
stage, events are required to have a muon with high $p_{\mathrm{T}}$ or a
hadron, photon or electron with high transverse energy in the calorimeters.
For hadrons, the transverse energy threshold is
3.5$\mathrm{\,Ge\kern-1.00006ptV}$. In the subsequent software trigger, at
least one charged particle must have a transverse momentum
$\mbox{$p_{\mathrm{T}}$}>1.7$${\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$ and be
inconsistent with originating from a PV. Finally, the tracks of two or more of
the final-state particles are required to form a vertex that is significantly
displaced from the PVs. The final state particles that are identified as kaons
are required to have a combined invariant mass consistent with that of the
$\phi$ meson.

In the simulation, $pp$ collisions are generated using Pythia8 [18,
*Sjostrand:2007gs] with a specific LHCb configuration [20]. Decays of hadronic
particles are described by EvtGen [21], in which final-state radiation is
generated using Photos [22]. The interaction of the generated particles with
the detector, and its response, are implemented using the Geant4 toolkit [23,
*Agostinelli:2002hh] as described in Ref. [25]. The decays of ${\mathchar
28931\relax}^{0}_{b}$ baryons are modelled according to a phase-space
description. Differences in the efficiencies of protons and anti-protons, at
the sub-percent level, are accounted for with the Geant4 implementation of the
detector description.

##  3 Selection

The ${{{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi}$ and
${{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi}$
decays are reconstructed through the ${\mathchar
28931\relax}\rightarrow{p}{{\pi}^{-}}$, ${{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ and $\phi\rightarrow{{K}^{+}}{{K}^{-}}$
final states, where the inclusion of charge conjugate processes is implied
throughout the paper. Decays of ${\mathchar
28931\relax}\\!\rightarrow{p}{{\pi}^{-}}$ and
${{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ are reconstructed in two different
categories. The first category contains $\mathchar 28931\relax$
(${K}^{0}_{\mathrm{\scriptscriptstyle S}}$) hadrons that decay inside the
vertex detector acceptance and the second contains $\mathchar 28931\relax$
(${K}^{0}_{\mathrm{\scriptscriptstyle S}}$) hadrons that decay outside. These
categories are referred to as _long_ and _downstream_ , respectively. The high
resolution of the vertex detector leads to enhanced momentum, vertex, and mass
resolutions for candidates in the long category relative to downstream
candidates.

Boosted decision trees (BDTs) [26, 27] are used to separate signal from
background. Different BDTs are trained for decays where the daughter tracks of
the ${\mathchar 28931\relax}$ (${K}^{0}_{\mathrm{\scriptscriptstyle S}}$)
hadron are classified as long or downstream and according to whether the data
was collected in 2011 (7$\mathrm{\,Te\kern-1.00006ptV}$) or 2012
(8$\mathrm{\,Te\kern-1.00006ptV}$), yielding eight separate BDTs in total. The
set of input variables used to train the ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$
(${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$) BDTs
consists of the ${\mathchar 28931\relax}^{0}_{b}$ (${B}^{0}$) vertex fit
quality, $p_{\mathrm{T}}$, $\eta$, the difference in $\chi^{2}$ of the PV
reconstructed with and without the candidate ($\chi^{2}_{\rm IP}$), the flight
distance squared divided by the associated variance ($\chi^{2}_{\rm FD}$), the
angle between the momentum vector and the vector from the PV to the decay
vertex, the ${\mathchar 28931\relax}$ (${K}^{0}_{\mathrm{\scriptscriptstyle
S}}$) vertex fit quality, and the $p_{\mathrm{T}}$ and $\eta$ of the $\phi$
and the ${\mathchar 28931\relax}$ (${K}^{0}_{\mathrm{\scriptscriptstyle S}}$)
hadrons. The minimum and maximum values of the $p_{\mathrm{T}}$ and $\eta$
associated to the final state particles are also included. In addition, the
BDT trained on the long category uses the $\chi^{2}_{\rm IP}$ and
$\chi^{2}_{\rm FD}$ of the ${\mathchar 28931\relax}$
(${K}^{0}_{\mathrm{\scriptscriptstyle S}}$) with respect to the associated PV.
A PV is reconstructed by requiring a minimum of five good quality tracks that
are consistent with originating from the same location within the luminous
region. Before the BDTs are trained, initial loose requirements are imposed on
the input variables. The BDTs are trained using simulated candidates for the
signal and data sidebands for the background. For the training samples, the
signal region is defined as being within
150${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ of the known ${\mathchar
28931\relax}^{0}_{b}$ (${B}^{0}$) mass [28]. In addition, the
${{K}^{+}}{{K}^{-}}$ invariant mass is required to be within
20${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ of the known $\phi$ mass and the
${p}{{\pi}^{-}}$ invariant mass is required to be within
15${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ of the known $\mathchar
28931\relax$ mass [28]. The sidebands are defined to be within
500${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ of the known ${\mathchar
28931\relax}^{0}_{b}$ (${B}^{0}$) mass excluding the signal region.

The figure of merit used to determine the requirement imposed on the
${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$
BDT output is defined as $\varepsilon/(3/2+\sqrt{N_{\rm bkg}})$ [29], where
$\varepsilon$ is the signal efficiency, and $N_{\rm bkg}$ is the number of
background events. This figure of merit is optimised for detection at three
standard deviations of decay modes not previously observed. The signal
efficiency is obtained from simulated signal candidates and the number of
background events is calculated from fits to the data sidebands interpolated
to the signal region. This optimisation procedure is performed separately for
each BDT.

In contrast to the ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ BDTs, the optimum response requirement for the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ BDTs is
chosen based on a figure of merit defined as $N_{\rm sig}/\sqrt{N_{\rm
sig}+N_{\rm bkg}}$, where $N_{\rm sig}$ is the number of signal events,
estimated from the BDT efficiency on simulated datasets normalised using the
known branching fraction of the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ decay
[15], and $N_{\rm bkg}$ is the expected number of background candidates in the
signal region, extrapolated from the ${B}^{0}$ sidebands. This figure of merit
is chosen as the ${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\phi$ branching fraction is well measured and is optimised separately for
each classifier.

##  4 Mass fit model

For both the ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ and
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ decay
modes, a three-dimensional fit is employed to determine the signal candidate
yields. In the ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ case, the three dimensions are the
${p}{{\pi}^{-}}{{K}^{+}}{{K}^{-}}$, ${p}{{\pi}^{-}}$, and ${{K}^{+}}{{K}^{-}}$
invariant masses, while in the fit to determine the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$
candidate yield, the three dimensions are the
${{\pi}^{+}}{{\pi}^{-}}{{K}^{+}}{{K}^{-}}$, ${{\pi}^{+}}{{\pi}^{-}}$, and
${{K}^{+}}{{K}^{-}}$ invariant masses.

Four components are present in the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ mass
fit: the signal ${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\phi$ component, the
${{B}^{0}}\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}{{K}^{+}}{{K}^{-}}$ non-resonant contribution, a
${{{\pi}^{+}}{{\pi}^{-}}{{K}^{+}}{{K}^{-}}}$ combinatorial component, along
with a true ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ component combined with
two random kaons. The
${{B}^{0}}\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}{{K}^{+}}{{K}^{-}}$ non-resonant component has been observed by the BaBar
[30], Belle [6] and LHCb [31] collaborations. This is separated from the
signal decay through the different ${{K}^{+}}{{K}^{-}}$ invariant mass line
shapes. No significant partially reconstructed background, in which one or
more of the final state particles are missed, is found in the ${B}^{0}$ mass
region. Peaking backgrounds, from decays in which at least one of the final
state particles has been misidentified, are suppressed by the narrow
${{K}^{+}}{{K}^{-}}$ mass window around the $\phi$ meson and are treated as
systematic uncertainties.

The ${B}^{0}$ signal is modelled with the same modified Gaussian function as
used in Ref. [32]. The modified Gaussian gives extra degrees of freedom to
accommodate extended tails far from the mean. The $\phi$ signal is modelled
with a relativistic Breit-Wigner shape [33] convolved with a Gaussian
resolution function. The ${{K}^{0}_{\mathrm{\scriptscriptstyle S}}}$ signal is
parametrised by the sum of two Gaussian functions with a common mean. Decays
from real ${B}^{0}$ mesons to the ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$
${K}^{+}$ ${K}^{-}$ final state in which the ${K}^{+}$ ${K}^{-}$ pair is non-
resonant are described by the same ${B}^{0}$ and
${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ line shapes as the signal, but with
a phase-space factor to describe the non-resonant kaon pairs. The phase-space
factor is given by the expression $(m^{2}-(2m_{K})^{2})/m^{2}$, where $m$ is
the ${K}^{+}$ ${K}^{-}$ invariant mass and $m_{K}$ is fixed to the value of
the charged kaon mass. The use of a Flatté function [34] rather than a phase-
space factor to describe a possible scalar component under the $\phi$
resonance is found to have a negligible effect on the results and is therefore
not included. The combinatorial background is modelled by exponential
functions in all three mass dimensions.

A simultaneous fit to the long and downstream datasets is performed. The
${B}^{0}$ resolution, modified Gaussian tail parameters and resolutions and
fractions of the ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ Gaussian functions
are constrained to values obtained from a fit to simulated data, performed
separateley for long and downstream datasets. The total yield and fraction in
the downstream dataset are left as free parameters for each component.

The fit to the ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ channel uses the same fit model as the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ control
channel: a modified Gaussian function is used to describe the ${\mathchar
28931\relax}^{0}_{b}$ mass shape, a double Gaussian model to describe the
${\mathchar 28931\relax}$ shape, and a relativistic Breit-Wigner convolved
with a Gaussian resolution function to describe that of the $\phi$ resonance.
Due to the relatively unexplored mass spectra present in the ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ decay, the
background contributions have been identified using the data sidebands. In the
final fit, four components are present. These are the signal ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ component, the
${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{\mathchar
28931\relax}{{K}^{+}}{{K}^{-}}$ non-resonant component in which the
${{K}^{+}}{{K}^{-}}$ dimension is described using the phase-space factor
defined previously, combinatorial components with true $\phi$ or $\mathchar
28931\relax$ resonances, and a component that has a combinatorial origin in
all three mass dimensions. Combinatorial backgrounds are modelled by
exponential functions in each fit dimension. As for the case of the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ fit,
the total yield and fraction in the downstream dataset are left as free
parameters for each component. In addition, the same parameters are
constrained to simulated data as in the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ fit.

##  5 Branching fraction measurement

The ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ branching fraction is obtained from the relation

$\displaystyle{\mathcal{B}}({{{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi})=\frac{\epsilon^{\rm
tot}_{{{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\phi}}}{\epsilon^{\rm tot}_{{{{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi}}}$
$\displaystyle\cdot\frac{f_{{d}}}{f_{{{\mathchar
28931\relax}^{0}_{b}}}}\cdot\frac{N_{{{{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi}}}{N_{{{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\phi}}}\cdot\frac{{\mathcal{B}}({{B}^{0}}\rightarrow{{K}^{0}}\phi)}{2}\cdot\frac{{\mathcal{B}}({{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\rightarrow{{\pi}^{+}}{{\pi}^{-}})}{{\mathcal{B}}({\mathchar
28931\relax}\rightarrow{p}{{\pi}^{-}})},$ (1)

where $\epsilon^{\rm tot}$ denotes the combined efficiency of the candidate
reconstruction, the offline selection, the trigger requirements, and the
efficiency of detector acceptance; $f_{{d}({{\mathchar
28931\relax}^{0}_{b}})}$ denotes the fraction of $b$ quarks that hadronise to
${B}^{0}$ (${\mathchar 28931\relax}^{0}_{b}$) hadrons. The ratio is taken from
the LHCb measured value $f_{{{\mathchar
28931\relax}^{0}_{b}}}/f_{{d}}=0.387\pm 0.033$ [35]. The extra factor $1/2$ in
Eq. 1 accounts for the fact that only half of ${K}^{0}$ mesons will decay as
${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ mesons. The value of the
${{B}^{0}}\\!\rightarrow{{K}^{0}}\phi$ branching fraction is taken to be
$(7.3^{+0.7}_{-0.6})\times 10^{-6}$ [15], while the PDG values of the
$\mathchar 28931\relax$ and ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$
branching fractions are used [28].

The reconstruction, selection and software trigger efficiencies, as well as
the acceptance of the LHCb detector, are determined from simulated samples,
using data-driven correction factors where necessary. The different
interaction cross-sections of the final-state particles with the detector
material is accounted for using simulated datasets.

For the case of the hardware trigger, the efficiency of events triggered by
the signal candidate is determined from control samples of
${{D}^{0}}\\!\rightarrow{{K}^{-}}{{\pi}^{+}}$ and ${\mathchar
28931\relax}\\!\rightarrow{p}{{\pi}^{-}}$ decays. The efficiency of events
triggered independently of the signal candidate is determined from simulation.
The agreement between data and simulation for the distributions of the
variables used in the BDT is verified with the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ data.

Data-driven corrections for the reconstruction efficiency of tracks
corresponding to the long category are obtained from
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ samples using a tag-and-probe
method [36]. This is applied after a separate weighting to ensure agreement in
detector occupancy between data and simulation. For measurements of the
relative branching fraction of ${{\mathchar
28931\relax}^{0}_{b}}\rightarrow{\mathchar 28931\relax}\phi$ to
${{B}^{0}}\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$, the final
state differs by substituting the proton from the decay of the $\mathchar
28931\relax$ with a pion. However, due to the differences in the kinematics of
the pions from the ${\mathchar 28931\relax}$ and the
${{K}^{0}_{\mathrm{\scriptscriptstyle S}}}$ decays, the distinct correction
factors for both daughters of the $\mathchar 28931\relax$ and
${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ are considered. In addition to the
track reconstruction efficiency, the vertexing efficiency of long-lived
particles contains disagreement between data and simulation. The corresponding
correction factors for the long and downstream datasets are determined
separately from ${{D}^{0}}\rightarrow\phi{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}$ decays.

Figure 2: Fit projections to the ${{p}{{\pi}^{-}}{{K}^{+}}{{K}^{-}}}$invariant
mass in the (a) long and (b) downstream datasets, the
${{{K}^{+}}{{K}^{-}}}$invariant mass in the (c) long and (d) downstream
datasets, and the ${{p}{{\pi}^{-}}}$invariant mass in the (e) long and (f)
downstream datasets. The total fit projection is given by the blue solid line.
The blue and green dotted lines represent the $\phi$ \+ ${\mathchar
28931\relax}$ and pure combinatorial fit components, respectively. The red and
magenta dashed lines represent the ${{{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi}$ signal and
the ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{\mathchar
28931\relax}{{K}^{+}}{{K}^{-}}$ non-resonant components, respectively. Black
points represent the data. Data uncertainties are Poisson 68% confidence
intervals.

Figure 3: Fit projections to the
${{{\pi}^{+}}{{\pi}^{-}}{{K}^{+}}{{K}^{-}}}$invariant mass in the (a) long and
(b) downstream datasets, the ${{{K}^{+}}{{K}^{-}}}$invariant mass in the (c)
long and (d) downstream datasets, and the ${{{\pi}^{+}}{{\pi}^{-}}}$ invariant
mass in the (e) long and (f) downstream datasets. The total fit projection is
given by the blue solid line. The green and blue dotted lines represent the
combinatorial and ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ \+ random
${K}^{+}$ ${K}^{-}$ fit components, respectively. The red and magenta dashed
lines represent the
${{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi}$
signal and the ${{B}^{0}}\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}{{K}^{+}}{{K}^{-}}$ non-resonant components, respectively. Black points
represent the data. Data uncertainties are Poisson 68% confidence intervals.

The yields of the ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ signal and
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ control
mode are determined from simultaneous extended unbinned maximum likelihood
fits to the respective datasets divided according to the data-taking period
and also according to whether the ${\mathchar
28931\relax}\,({{K}^{0}_{\mathrm{\scriptscriptstyle S}}})$ decay products are
reconstructed as long or downstream tracks. Efficiencies are applied to each
dataset individually. The projections of the fit result to ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ data are shown
in Fig. 2. The fitted yields are $350\pm 24$ and $89\pm 13$ for the
${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$ and
${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$
decay modes, respectively. The statistical significance of the ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ decay,
determined according to Wilks’ theorem [37] from the difference in the
likelihood value of the fits with and without the ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ component, is
found to be $6.5$ standard deviations. With the systematic uncertainties
discussed below included, the significance of the observed ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ decay yield is
calculated to be $5.9$ standard deviations. The projections of the fit result
to the ${{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi$
data are shown in Fig. 3. The fit is found to describe the data well in all
three dimensions and a clear peak from the control mode is seen.

The systematic contributions to the branching fraction uncertainty budget are
summarised in Table 1. The largest contributions to the systematic
uncertainties result from data-driven corrections applied to simulated data
along with the mass model used to determine the signal yields.

Table 1: Systematic uncertainty contributions to the branching fraction ratio. Source | Uncertainty (%)
---|---
Mass model | 3.0
Simulation sample size | 2.2
Tracking efficiency | 0.5
Vertex efficiency | 2.6
Hardware trigger | 2.8
Selection efficiency | 4.1
Peaking background | 0.1
Total | 6.7

Signal mismodelling is accounted for using a one-dimensional kernel estimate
for the description of the simulated mass distributions [38]. Background
mismodelling is accounted for using a linear function. The kernel estimate is
used in both the signal and control channels to describe the ${\mathchar
28931\relax}^{0}_{b}$, ${B}^{0}$, ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$,
and ${\mathchar 28931\relax}$ line shapes. In order to determine the
systematic uncertainties, 1000 pseudoexperiments are generated with the
alternative model and are subsequently fitted with the nominal model. The
average difference between the generated and fitted yield values is taken as
the systematic uncertainty. This leads to uncertainties of 3.0% and 0.6% for
the signal and control mode yields, respectively.

Systematic uncertainties associated with the efficiency corrections from
simulated datasets are considered. The limited size of the simulated sample
gives rise to an uncertainty of 2.2%. The main uncertainties in the tracking
and vertexing correction factors arise from the limited size of the control
sample, which leads to uncertainties of 0.5% and 2.6%, respectively. For the
case of the trigger efficiency, uncertainties related to the software trigger
cancel between the signal and control modes, as the software trigger decision
is made only on the decay products of the $\phi$ meson. Uncertainties in the
efficiency of the hardware trigger selections are estimated using data-driven
methods, for which an uncertainty of 2.8% is applied. The BDTs used to select
signal and control modes use the same input variables. Biases could exist if
the simulation mismodels these variables differently for signal and control
modes. In order to quantify this effect, the control mode is selected with the
same classifier as the signal decay. The difference in the measured branching
fraction is found to be 4.1%.

The ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow\mathchar
28934\relax^{0}(\rightarrow{\mathchar 28931\relax}\gamma){{K}^{+}}{{K}^{-}}$
and ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{p}{{K}^{-}}\phi$ decay modes
are found to be the only significant peaking background contributions.
However, for the case of the ${{\mathchar
28931\relax}^{0}_{b}}\rightarrow{p}{{K}^{-}}\phi$ decay, the resulting
candidates are reconstructed in the long dataset only. With the assumption
that the branching fraction for this decay is the same size as for the signal,
the contribution is $<1\%$ compared to the ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ decay and far
from the ${\mathchar 28931\relax}^{0}_{b}$ signal region, and is therefore
ignored. In order to determine the shape in the
${p}{{\pi}^{-}}{{K}^{+}}{{K}^{-}}$ spectrum of the ${{\mathchar
28931\relax}^{0}_{b}}\rightarrow\mathchar 28934\relax^{0}{{K}^{+}}{{K}^{-}}$
decay, a sample of ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow\mathchar
28934\relax^{0}{{K}^{+}}{{K}^{-}}$ simulated events is used with a requirement
that the ${{K}^{+}}{{K}^{-}}$ invariant mass is within
30${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ of the nominal $\phi$ mass. The
inclusion of an additional fit component using the shape from simulation is
found to have a small effect on the signal yield at the level of 0.1%, which
is assigned as a systematic uncertainty. For the case of the
${{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle S}}}\phi}$
control mode, no peaking background contributions have been identified.

The branching fraction ratio is measured to be

$\displaystyle\frac{{\mathcal{B}}({{{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi})}{{\mathcal{B}}({{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\phi})}=0.55\pm 0.11\mathrm{\,(stat)}\pm 0.04\mathrm{\,(syst)}\pm
0.05\,(f_{d}/f_{{{\mathchar 28931\relax}^{0}_{b}}}).$

The use of the world average value of
${\mathcal{B}}({{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\phi})=(3.65\,^{+0.35}_{-0.30})\times 10^{-6}$ [15] gives the final result
of

$\displaystyle{\mathcal{B}}({{{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi})/10^{-6}$
$\displaystyle=$ $\displaystyle 5.18\pm 1.04\mathrm{\,(stat)}\pm
0.35\mathrm{\,(syst)}\,^{+0.50}_{-0.43}\,({\mathcal{B}}({{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\phi}))\pm 0.44\,(f_{d}/f_{{{\mathchar 28931\relax}^{0}_{b}}}).$

##  6 Triple-product asymmetries

The ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ decay is a spin-$1/2$ to spin-$1/2$ plus vector transition.
Five angles are needed to describe this decay since ${\mathchar
28931\relax}^{0}_{b}$ baryons may potentially be produced with a transverse
polarisation in proton-proton collisions [13], as shown in Fig. 4. The angle
$\theta$ is defined as the polar angle of the ${\mathchar 28931\relax}$ baryon
in the ${\mathchar 28931\relax}^{0}_{b}$ rest frame with respect to the normal
vector defined through

$\displaystyle\hat{n}=\frac{\vec{p}_{1}\times\vec{p}_{{{\mathchar
28931\relax}^{0}_{b}}}}{|\vec{p}_{1}\times\vec{p}_{{{\mathchar
28931\relax}^{0}_{b}}}|},$ (2)

where $\vec{p}_{1}$ is the momentum of an incoming proton and
$\vec{p}_{{{\mathchar 28931\relax}^{0}_{b}}}$ is the momentum of the
${\mathchar 28931\relax}^{0}_{b}$ baryon. The angles $\theta_{\mathchar
28931\relax}$ and $\Phi_{\mathchar 28931\relax}$ are defined as the polar and
azimuthal angles of the proton from the decay of the ${\mathchar 28931\relax}$
baryon in the ${\mathchar 28931\relax}$ rest frame. The angles $\theta_{\phi}$
and $\Phi_{\phi}$ are defined as the polar and azimuthal angles of the
${K}^{+}$ meson in the rest frame of the $\phi$ meson.

Figure 4: Decay angles for the ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ decay, where
the angles are defined in the text.

Triple-product asymmetries, which are odd under time-reversal, have been
proposed by Leitner and Ajaltouni using the azimuthal angles
$\Phi_{n_{i}},i\in\\{{\mathchar 28931\relax},\phi\\}$, defined as [12]

$\displaystyle\cos\Phi_{n_{i}}$ $\displaystyle=\vec{e}_{Y}\cdot\vec{u}_{i},$
(3) $\displaystyle\sin\Phi_{n_{i}}$
$\displaystyle=\vec{e}_{Z}\cdot(\vec{e}_{Y}\times\vec{u}_{i}),$ (4)

where

$\displaystyle\vec{u}_{i}=\frac{\vec{e}_{Z}\times\hat{n}_{i}}{|\vec{e}_{Z}\times\hat{n}_{i}|}.$
(5)

The basis $\\{\vec{e}_{X},\vec{e}_{Y},\vec{e}_{Z}\\}$ is defined in the
${\mathchar 28931\relax}^{0}_{b}$ rest frame, in which $\vec{e}_{Z}$ is
parallel to $\hat{n}$, $\vec{e}_{X}$ is chosen to be parallel to the momentum
of the incoming proton, and $\hat{n}_{{\mathchar 28931\relax}(\phi)}$ is the
normal vector to the ${\mathchar 28931\relax}(\phi)$ decay plane, defined
through

$\displaystyle\hat{n}_{\mathchar 28931\relax}$
$\displaystyle=\frac{\vec{p}_{p}\times\vec{p}_{\pi}}{|\vec{p}_{p}\times\vec{p}_{\pi}|},$
(6) $\displaystyle\hat{n}_{\phi}$
$\displaystyle=\frac{\vec{p}_{{{K}^{+}}}\times\vec{p}_{{{K}^{-}}}}{|\vec{p}_{{{K}^{+}}}\times\vec{p}_{{{K}^{-}}}|}.$
(7)

Asymmetries in $\cos\Phi_{n_{i}}$ and $\sin\Phi_{n_{i}}$, where
$i\in\\{{\mathchar 28931\relax},\phi\\}$, are defined as

$\displaystyle A_{i}^{c}$
$\displaystyle=\frac{N^{+,c}_{i}-N^{-,c}_{i}}{N^{+,c}_{i}+N^{-,c}_{i}},$ (8)
$\displaystyle A_{i}^{s}$
$\displaystyle=\frac{N^{+,s}_{i}-N^{-,s}_{i}}{N^{+,s}_{i}+N^{-,s}_{i}},$ (9)

where $N^{+(-),c}_{i}$ and $N^{+(-),s}_{i}$ denote the number of candidates
for which the $\cos\Phi_{n_{i}}$ and $\sin\Phi_{n_{i}}$ observables are
positive (negative), respectively.

The asymmetries $A_{\mathchar 28931\relax}^{c,s}$ and $A_{\phi}^{c,s}$ are
determined experimentally through a simultaneous unbinned maximum likelihood
fit to datasets in which the relevant observables are positive and negative.
The fit construction and observables are identical to that used for the
branching fraction measurement. However, the yields for each dataset are
parametrised in terms of the total yield, $N_{j}$, and the asymmetry, $A_{j}$,
for fit component $j$ as

$\displaystyle N^{+}_{j}$ $\displaystyle=\frac{N^{j}}{2}(1+A_{j}),$ (10)
$\displaystyle N^{-}_{j}$ $\displaystyle=\frac{N^{j}}{2}(1-A_{j}).$ (11)

Distributions of the $\sin\Phi_{n_{({\mathchar 28931\relax},\phi)}}$ and
$\cos\Phi_{n_{({\mathchar 28931\relax},\phi)}}$ observables from ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ data have been
extracted using the sPlot method [39] and are provided in Fig. 5. The
numerical values of the fitted asymmetries are given in Table 2.

Figure 5: Distributions of the angular observables: (a) $\sin\Phi_{n_{\mathchar 28931\relax}}$, (b) $\cos\Phi_{n_{\mathchar 28931\relax}}$, (c) $\sin\Phi_{n_{\phi}}$, (d) $\cos\Phi_{n_{\phi}}$ from weighted ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ data. Table 2: Asymmetries measured from ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi$ data events. Asymmetry | Fit value
---|---
$A^{c}_{\mathchar 28931\relax}$ | $-0.22\pm 0.12$
$A^{s}_{\mathchar 28931\relax}$ | $0.13\pm 0.12$
$A^{c}_{\phi}$ | $-0.01\pm 0.12$
$A^{s}_{\phi}$ | $-0.07\pm 0.12$

Mismodelling of the mass components could lead to background contamination in
the determination of the asymmetries. In the determination of the uncertainty
related to the mass model, two contributions are considered. These are the
line shape models and the background asymmetries. The effects of the line
shapes are quantified using the same method as the branching fraction
measurement, i.e. the generation of datasets with a one-dimensional kernel
estimate of the simulation mass distributions in addition to modification of
the background description. In the nominal fit, components that are not from
the ${{\mathchar 28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar
28931\relax}\phi$ signal have zero asymmetries. For background components this
is justified due to the uncorrelated kinematics of the ${K}^{+}$ ${K}^{-}$ and
$p$ ${\pi}^{-}$ systems. However, the non-resonant ${{\mathchar
28931\relax}^{0}_{b}}\rightarrow{\mathchar 28931\relax}{{K}^{+}}{{K}^{-}}$
contribution could have non-zero asymmetries. The systematic uncertainty due
to the assumption of zero background asymmetries is determined through
comparing the nominal fit against the fit with all possible asymmetries
allowed to vary freely.

Efficiencies are found to be independent of the $\sin\Phi_{n_{i}}$ and
$\cos\Phi_{n_{i}}$ observables. The systematic uncertainty due to the angular
acceptance is then taken from the statistical uncertainty in fits to the
simulated datasets, after the application of an appropriate weighting to
account for the differences between data and simulation. The resolutions of
the angular observables are found from simulated events to be
32.3$\mathrm{\,mrad}$ and 22.1$\mathrm{\,mrad}$ for the $\Phi_{n_{\mathchar
28931\relax}}$ and $\Phi_{n_{\phi}}$ angles, respectively. The uncertainty due
to bin migration is then assigned assuming maximal asymmetry and leads to
minor uncertainties of 0.007 for the $\Phi_{n_{\phi}}$ angle and 0.010 for the
$\Phi_{n_{\mathchar 28931\relax}}$ angle. Systematic contributions to the
triple-product uncertainty budget are summarised in Table 3.

Table 3: Systematic uncertainty contributions to the triple-product asymmetries. Source | $A_{\mathchar 28931\relax}^{c}$ | $A_{\mathchar 28931\relax}^{s}$ | $A_{\phi}^{c}$ | $A_{\phi}^{s}$
---|---|---|---|---
Mass model | 0.061 | 0.051 | 0.026 | 0.009
Angular acceptance | 0.010 | 0.010 | 0.010 | 0.010
Angular resolution | 0.008 | 0.008 | 0.005 | 0.005
Total | 0.062 | 0.053 | 0.028 | 0.014

##  7 Summary

A search for the ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{\mathchar
28931\relax}\phi$ decay is presented based on a dataset of
3.0$\mbox{\,fb}^{-1}$ collected by the LHCb experiment in 2011 and 2012. The
decay is observed for the first time with a significance of $5.9$ standard
deviations including systematic uncertainties. The branching fraction is found
to be

$\displaystyle{\mathcal{B}}({{{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}\phi})/10^{-6}$
$\displaystyle=$ $\displaystyle 5.18\pm 1.04\mathrm{\,(stat)}\pm
0.35\mathrm{\,(syst)}\,^{+0.50}_{-0.43}\,({\mathcal{B}}({{{B}^{0}}\\!\rightarrow{{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\phi}))\pm 0.44\,(f_{d}/f_{{{\mathchar 28931\relax}^{0}_{b}}}).$

Triple-product asymmetries are measured to be

$\displaystyle A^{c}_{\mathchar 28931\relax}$ $\displaystyle=-0.22\pm
0.12\mathrm{\,(stat)}\pm 0.06\mathrm{\,(syst)},$ $\displaystyle
A^{s}_{\mathchar 28931\relax}$ $\displaystyle=\phantom{-}0.13\pm
0.12\mathrm{\,(stat)}\pm 0.05\mathrm{\,(syst)},$ $\displaystyle A^{c}_{\phi}$
$\displaystyle=-0.01\pm 0.12\mathrm{\,(stat)}\pm 0.03\mathrm{\,(syst)},$
$\displaystyle A^{s}_{\phi}$ $\displaystyle=-0.07\pm 0.12\mathrm{\,(stat)}\pm
0.01\mathrm{\,(syst)},$

and are consistent with zero. Data collected by the LHCb experiment in the
forthcoming years will improve the statistical precision of these measurements
and enable the dynamics of ${b}\rightarrow{s}$ transitions in beauty baryons
to be probed in greater detail, which will greatly enhance the reach of
searches for physics beyond the SM.

## Acknowledgements

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 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). We acknowledge the computing resources that
are provided by CERN, IN2P3 (France), KIT and DESY (Germany), INFN (Italy),
SURF (The Netherlands), PIC (Spain), GridPP (United Kingdom), RRCKI and Yandex
LLC (Russia), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil), PL-GRID
(Poland) and OSC (USA). We are indebted to the communities behind the multiple
open source software packages on which we depend. Individual groups or members
have received support from AvH Foundation (Germany), 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 and Yandex
LLC (Russia), GVA, XuntaGal and GENCAT (Spain), Herchel Smith Fund, The Royal
Society, Royal Commission for the Exhibition of 1851 and the Leverhulme Trust
(United Kingdom).

## References

  * [1] Y. Shimizu, M. Tanimoto, and K. Yamamoto, _Supersymmetry contributions to CP violations in $b\rightarrow s$ and $b\rightarrow d$ transitions taking account of new data_, Phys. Rev. D87 (2013) 056004, arXiv:1212.6486
  * [2] A. Datta, M. Imbeault, and D. London, _$B\rightarrow\phi K^{*}$ , $B\rightarrow\phi K_{S}$ and new physics_, Phys. Lett. B671 (2009) 256, arXiv:0811.2957
  * [3] T. Moroi, _CP violation in $B^{0}\rightarrow\phi K_{S}$ in SUSY GUT with right-handed neutrinos_, Phys. Lett. B493 (2000) 366, arXiv:hep-ph/0007328
  * [4] LHCb collaboration, R. Aaij et al., _Measurement of $C\\!P$ violation in $B_{s}^{0}\rightarrow\phi\phi$ decays_, Phys. Rev. D90 (2014) 052011, arXiv:1407.2222
  * [5] LHCb collaboration, R. Aaij et al., _First measurement of the $C\\!P$-violating phase in $B^{0}_{s}\rightarrow\phi\phi$ decays_, Phys. Rev. Lett. 110 (2013) 241802, arXiv:1303.7125
  * [6] Belle collaboration, K. Abe et al., _Measurement of time-dependent CP-violating asymmetries in $B^{0}\rightarrow\phi K_{S}^{0}$, ${{K}^{+}}{{K}^{-}}K_{S}^{0}$, and $\eta^{\prime}K_{S}^{0}$ decays_, Phys. Rev. Lett. 91 (2003) 261602, arXiv:hep-ex/0308035
  * [7] BaBar collaboration, J. P. Lees et al., _Study of CP violation in Dalitz-plot analyses of $B^{0}\rightarrow{{K}^{+}}{{K}^{-}}K_{S}^{0}$, ${{B}^{+}}\rightarrow{{K}^{+}}{{K}^{-}}{{K}^{+}}$, and ${{B}^{+}}\rightarrow K_{S}^{0}K_{S}^{0}{{K}^{+}}$_, Phys. Rev. D85 (2012) 112010, arXiv:1201.5897
  * [8] LHCb collaboration, R. Aaij et al., _Measurement of polarization amplitudes and $C\\!P$ asymmetries in $B^{0}\rightarrow\phi K^{*}(892)^{0}$_, JHEP 05 (2014) 069, arXiv:1403.2888
  * [9] M. Raidal, _$C\\!P$ asymmetry in $B\rightarrow\phi K_{S}$ decays in left-right models and its implications on $B_{s}$ decays_, Phys. Rev. Lett. 89 (2002) 231803, arXiv:hep-ph/0208091
  * [10] M. Gronau and J. L. Rosner, _Triple product asymmetries in $K$, $D_{(s)}$ and $B_{(s)}$ decays_, Phys. Rev. D84 (2011) 096013, arXiv:1107.1232
  * [11] S. K. Patra and A. Kundu, _$CPT$ violation and triple-product correlations in B decays_, Phys. Rev. D87 (2013), no. 11 116005, arXiv:1305.1417
  * [12] O. Leitner and Z. J. Ajaltouni, _Testing CP and time reversal symmetries with ${\mathchar 28931\relax}_{b}\rightarrow{\mathchar 28931\relax}V(1-)$ decays_, Nucl. Phys. Proc. Suppl. 174 (2007) 169, arXiv:hep-ph/0610189
  * [13] LHCb collaboration, R. Aaij et al., _Measurements of the $\Lambda_{b}^{0}\rightarrow J/\psi\Lambda$ decay amplitudes and the $\Lambda_{b}^{0}$ polarisation in $pp$ collisions at $\sqrt{s}=7$ TeV_, Phys. Lett. B724 (2013) 27, arXiv:1302.5578
  * [14] LHCb collaboration, R. Aaij et al., _Observations of ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{\mathchar 28931\relax}K^{+}\pi^{-}$ and ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{\mathchar 28931\relax}K^{+}K^{-}$ decays and searches for other ${{\mathchar 28931\relax}^{0}_{b}}$ and $\Xi_{b}^{0}$ decays to ${\mathchar 28931\relax}h^{+}h^{-}$ final states_, JHEP 1605 (2016) 081, arXiv:1603.00413
  * [15] Heavy Flavor Averaging Group, Y. Amhis et al., _Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of summer 2014_, arXiv:1412.7515, updated results and plots available at http://www.slac.stanford.edu/xorg/hfag/
  * [16] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [17] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [18] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [19] T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * [20] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , J. Phys. Conf. Ser. 331 (2011) 032047
  * [21] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [22] 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
  * [23] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [24] Geant4 collaboration, S. Agostinelli et al., _Geant4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [25] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [26] L. Breiman, J. H. Friedman, R. A. Olshen, and C. J. Stone, Classification and regression trees, Wadsworth international group, Belmont, California, USA, 1984
  * [27] R. E. Schapire and Y. Freund, _A decision-theoretic generalization of on-line learning and an application to boosting_ , J. Comput. Syst. Sci. 55 (1997) 119
  * [28] Particle Data Group, K. A. Olive et al., _Review of Particle Physics_ , Chin. Phys. C38 (2014) 090001
  * [29] G. Punzi, _Sensitivity of searches for new signals and its optimization_ , in Statistical Problems in Particle Physics, Astrophysics, and Cosmology (L. Lyons, R. Mount, and R. Reitmeyer, eds.), p. 79, 2003.  arXiv:physics/0308063
  * [30] BaBar collaboration, B. Aubert et al., _Measurement of CP asymmetries in $B^{0}\rightarrow\phi K^{0}$ and $B^{0}\rightarrow K^{+}K^{-}K^{0}_{S}$ decays_, Phys. Rev. D71 (2005) 091102, arXiv:hep-ex/0502019
  * [31] LHCb collaboration, R. Aaij et al., _Study of $B_{(s)}^{0}\rightarrow K^{0}_{S}h^{+}h^{\prime-}$ decays with first observation of $B_{s}^{0}\rightarrow K^{0}_{S}K^{\pm}\pi^{\mp}$ and $B_{s}^{0}\rightarrow K^{0}_{S}\pi^{+}\pi^{-}$_, JHEP 10 (2013) 143, arXiv:1307.7648
  * [32] LHCb collaboration, R. Aaij et al., _Observation of $C\\!P$ violation in $B^{\pm}\rightarrow DK^{\pm}$ decays_, Phys. Lett. B712 (2012) 203, Erratum ibid. B713 (2012) 351, arXiv:1203.3662
  * [33] HERA-B collaboration, I. Abt et al., _${K}^{*0}$ and $\phi$ meson production in proton-nucleus interactions at $\sqrt{s}$ = 41.6 GeV_, Eur. Phys. J. C50 (2007) 315, arXiv:hep-ex/0606049
  * [34] S. M. Flatté, _On the nature of $0^{+}$ mesons_, Phys. Lett. B63 (1976) 228
  * [35] LHCb collaboration, R. Aaij et al., _Study of the kinematic dependences of $\Lambda_{b}^{0}$ production in $pp$ collisions and a measurement of the $\Lambda_{b}^{0}\rightarrow\Lambda_{c}^{+}\pi^{-}$ branching fraction_, JHEP 08 (2014) 143, arXiv:1405.6842
  * [36] LHCb collaboration, R. Aaij et al., _Measurement of the track reconstruction efficiency at LHCb_ , JINST 10 (2015) P02007, arXiv:1408.1251
  * [37] S. S. Wilks, _The large-sample distribution of the likelihood ratio for testing composite hypotheses_ , Annals Math. Statist. 9 (1938) 60
  * [38] K. S. Cranmer, _Kernel estimation in high-energy physics_ , Comput. Phys. Commun. 136 (2001) 198, arXiv:hep-ex/0011057
  * [39] M. Pivk and F. R. Le Diberder, _sPlot: A statistical tool to unfold data distributions_ , Nucl. Instrum. Meth. A555 (2005) 356, arXiv:physics/0402083

LHCb collaboration

R. Aaij39, C. Abellán Beteta41, B. Adeva38, M. Adinolfi47, Z. Ajaltouni5, S.
Akar6, J. Albrecht10, F. Alessio39, M. Alexander52, S. Ali42, G. Alkhazov31,
P. Alvarez Cartelle54, A.A. Alves Jr58, S. Amato2, S. Amerio23, Y. Amhis7, L.
An3,40, L. Anderlini18, G. Andreassi40, M. Andreotti17,g, J.E. Andrews59, R.B.
Appleby55, O. Aquines Gutierrez11, F. Archilli39, P. d’Argent12, A.
Artamonov36, M. Artuso60, E. Aslanides6, G. Auriemma26,n, M. Baalouch5, S.
Bachmann12, J.J. Back49, A. Badalov37, C. Baesso61, S. Baker54, W. Baldini17,
R.J. Barlow55, C. Barschel39, S. Barsuk7, W. Barter39, V. Batozskaya29, V.
Battista40, A. Bay40, L. Beaucourt4, J. Beddow52, F. Bedeschi24, I. Bediaga1,
L.J. Bel42, V. Bellee40, N. Belloli21,k, I. Belyaev32, E. Ben-Haim8, G.
Bencivenni19, S. Benson39, J. Benton47, A. Berezhnoy33, R. Bernet41, A.
Bertolin23, F. Betti15, M.-O. Bettler39, M. van Beuzekom42, S. Bifani46, P.
Billoir8, T. Bird55, A. Birnkraut10, A. Bizzeti18,i, T. Blake49, F. Blanc40,
J. Blouw11, S. Blusk60, V. Bocci26, A. Bondar35, N. Bondar31,39, W.
Bonivento16, A. Borgheresi21,k, S. Borghi55, M. Borisyak67, M. Borsato38, M.
Boubdir9, T.J.V. Bowcock53, E. Bowen41, C. Bozzi17,39, S. Braun12, M.
Britsch12, T. Britton60, J. Brodzicka55, E. Buchanan47, C. Burr55, A.
Bursche2, J. Buytaert39, S. Cadeddu16, R. Calabrese17,g, M. Calvi21,k, M.
Calvo Gomez37,p, P. Campana19, D. Campora Perez39, L. Capriotti55, A.
Carbone15,e, G. Carboni25,l, R. Cardinale20,j, A. Cardini16, P. Carniti21,k,
L. Carson51, K. Carvalho Akiba2, G. Casse53, L. Cassina21,k, L. Castillo
Garcia40, M. Cattaneo39, Ch. Cauet10, G. Cavallero20, R. Cenci24,t, M.
Charles8, Ph. Charpentier39, G. Chatzikonstantinidis46, M. Chefdeville4, S.
Chen55, S.-F. Cheung56, M. Chrzaszcz41,27, X. Cid Vidal39, G. Ciezarek42,
P.E.L. Clarke51, M. Clemencic39, H.V. Cliff48, J. Closier39, V. Coco58, J.
Cogan6, E. Cogneras5, V. Cogoni16,f, L. Cojocariu30, G. Collazuol23,r, P.
Collins39, A. Comerma-Montells12, A. Contu39, A. Cook47, M. Coombes47, S.
Coquereau8, G. Corti39, M. Corvo17,g, B. Couturier39, G.A. Cowan51, D.C.
Craik51, A. Crocombe49, M. Cruz Torres61, S. Cunliffe54, R. Currie54, C.
D’Ambrosio39, E. Dall’Occo42, J. Dalseno47, P.N.Y. David42, A. Davis58, O. De
Aguiar Francisco2, K. De Bruyn6, S. De Capua55, M. De Cian12, J.M. De
Miranda1, L. De Paula2, P. De Simone19, C.-T. Dean52, D. Decamp4, M.
Deckenhoff10, L. Del Buono8, N. Déléage4, M. Demmer10, D. Derkach67, O.
Deschamps5, F. Dettori39, B. Dey22, A. Di Canto39, F. Di Ruscio25, H.
Dijkstra39, F. Dordei39, M. Dorigo40, A. Dosil Suárez38, A. Dovbnya44, K.
Dreimanis53, L. Dufour42, G. Dujany55, K. Dungs39, P. Durante39, R.
Dzhelyadin36, A. Dziurda27, A. Dzyuba31, S. Easo50,39, U. Egede54, V.
Egorychev32, S. Eidelman35, S. Eisenhardt51, U. Eitschberger10, R. Ekelhof10,
L. Eklund52, I. El Rifai5, Ch. Elsasser41, S. Ely60, S. Esen12, H.M. Evans48,
T. Evans56, A. Falabella15, C. Färber39, N. Farley46, S. Farry53, R. Fay53, D.
Fazzini21,k, D. Ferguson51, V. Fernandez Albor38, F. Ferrari15, F. Ferreira
Rodrigues1, M. Ferro-Luzzi39, S. Filippov34, M. Fiore17,g, M. Fiorini17,g, M.
Firlej28, C. Fitzpatrick40, T. Fiutowski28, F. Fleuret7,b, K. Fohl39, M.
Fontana16, F. Fontanelli20,j, D. C. Forshaw60, R. Forty39, M. Frank39, C.
Frei39, M. Frosini18, J. Fu22, E. Furfaro25,l, A. Gallas Torreira38, D.
Galli15,e, S. Gallorini23, S. Gambetta51, M. Gandelman2, P. Gandini56, Y.
Gao3, J. García Pardiñas38, J. Garra Tico48, L. Garrido37, P.J. Garsed48, D.
Gascon37, C. Gaspar39, L. Gavardi10, G. Gazzoni5, D. Gerick12, E. Gersabeck12,
M. Gersabeck55, T. Gershon49, Ph. Ghez4, S. Gianì40, V. Gibson48, O.G.
Girard40, L. Giubega30, V.V. Gligorov39, C. Göbel61, D. Golubkov32, A.
Golutvin54,39, A. Gomes1,a, C. Gotti21,k, M. Grabalosa Gándara5, R. Graciani
Diaz37, L.A. Granado Cardoso39, E. Graugés37, E. Graverini41, G. Graziani18,
A. Grecu30, P. Griffith46, L. Grillo12, O. Grünberg65, B. Gui60, E.
Gushchin34, Yu. Guz36,39, T. Gys39, T. Hadavizadeh56, C. Hadjivasiliou60, G.
Haefeli40, C. Haen39, S.C. Haines48, S. Hall54, B. Hamilton59, X. Han12, S.
Hansmann-Menzemer12, N. Harnew56, S.T. Harnew47, J. Harrison55, J. He39, T.
Head40, A. Heister9, K. Hennessy53, P. Henrard5, L. Henry8, J.A. Hernando
Morata38, E. van Herwijnen39, M. Heß65, A. Hicheur2, D. Hill56, M. Hoballah5,
C. Hombach55, L. Hongming40, W. Hulsbergen42, T. Humair54, M. Hushchyn67, N.
Hussain56, D. Hutchcroft53, M. Idzik28, P. Ilten57, R. Jacobsson39, A.
Jaeger12, J. Jalocha56, E. Jans42, A. Jawahery59, M. John56, D. Johnson39,
C.R. Jones48, C. Joram39, B. Jost39, N. Jurik60, S. Kandybei44, W. Kanso6, M.
Karacson39, T.M. Karbach39,†, S. Karodia52, M. Kecke12, M. Kelsey60, I.R.
Kenyon46, M. Kenzie39, T. Ketel43, E. Khairullin67, B. Khanji21,39,k, C.
Khurewathanakul40, T. Kirn9, S. Klaver55, K. Klimaszewski29, M. Kolpin12, I.
Komarov40, R.F. Koopman43, P. Koppenburg42,39, M. Kozeiha5, L. Kravchuk34, K.
Kreplin12, M. Kreps49, P. Krokovny35, F. Kruse10, W. Krzemien29, W.
Kucewicz27,o, M. Kucharczyk27, V. Kudryavtsev35, A. K. Kuonen40, K. Kurek29,
T. Kvaratskheliya32, D. Lacarrere39, G. Lafferty55,39, A. Lai16, D. Lambert51,
G. Lanfranchi19, C. Langenbruch49, B. Langhans39, T. Latham49, C. Lazzeroni46,
R. Le Gac6, J. van Leerdam42, J.-P. Lees4, R. Lefèvre5, A. Leflat33,39, J.
Lefrançois7, E. Lemos Cid38, O. Leroy6, T. Lesiak27, B. Leverington12, Y. Li7,
T. Likhomanenko67,66, R. Lindner39, C. Linn39, F. Lionetto41, B. Liu16, X.
Liu3, D. Loh49, I. Longstaff52, J.H. Lopes2, D. Lucchesi23,r, M. Lucio
Martinez38, H. Luo51, A. Lupato23, E. Luppi17,g, O. Lupton56, N. Lusardi22, A.
Lusiani24, X. Lyu62, F. Machefert7, F. Maciuc30, O. Maev31, K. Maguire55, S.
Malde56, A. Malinin66, G. Manca7, G. Mancinelli6, P. Manning60, A. Mapelli39,
J. Maratas5, J.F. Marchand4, U. Marconi15, C. Marin Benito37, P. Marino24,t,
J. Marks12, G. Martellotti26, M. Martin6, M. Martinelli40, D. Martinez
Santos38, F. Martinez Vidal68, D. Martins Tostes2, L.M. Massacrier7, A.
Massafferri1, R. Matev39, A. Mathad49, Z. Mathe39, C. Matteuzzi21, A. Mauri41,
B. Maurin40, A. Mazurov46, M. McCann54, J. McCarthy46, A. McNab55, R.
McNulty13, B. Meadows58, F. Meier10, M. Meissner12, D. Melnychuk29, M. Merk42,
A Merli22,u, E Michielin23, D.A. Milanes64, M.-N. Minard4, D.S. Mitzel12, J.
Molina Rodriguez61, I.A. Monroy64, S. Monteil5, M. Morandin23, P. Morawski28,
A. Mordà6, M.J. Morello24,t, J. Moron28, A.B. Morris51, R. Mountain60, F.
Muheim51, D. Müller55, J. Müller10, K. Müller41, V. Müller10, M. Mussini15, B.
Muster40, P. Naik47, T. Nakada40, R. Nandakumar50, A. Nandi56, I. Nasteva2, M.
Needham51, N. Neri22, S. Neubert12, N. Neufeld39, M. Neuner12, A.D. Nguyen40,
C. Nguyen-Mau40,q, V. Niess5, S. Nieswand9, R. Niet10, N. Nikitin33, T.
Nikodem12, A. Novoselov36, D.P. O’Hanlon49, A. Oblakowska-Mucha28, V.
Obraztsov36, S. Ogilvy52, O. Okhrimenko45, R. Oldeman16,48,f, C.J.G.
Onderwater69, B. Osorio Rodrigues1, J.M. Otalora Goicochea2, A. Otto39, P.
Owen54, A. Oyanguren68, A. Palano14,d, F. Palombo22,u, M. Palutan19, J.
Panman39, A. Papanestis50, M. Pappagallo52, L.L. Pappalardo17,g, C.
Pappenheimer58, W. Parker59, C. Parkes55, G. Passaleva18, G.D. Patel53, M.
Patel54, C. Patrignani20,j, A. Pearce55,50, A. Pellegrino42, G. Penso26,m, M.
Pepe Altarelli39, S. Perazzini15,e, P. Perret5, L. Pescatore46, K. Petridis47,
A. Petrolini20,j, M. Petruzzo22, E. Picatoste Olloqui37, B. Pietrzyk4, M.
Pikies27, D. Pinci26, A. Pistone20, A. Piucci12, S. Playfer51, M. Plo
Casasus38, T. Poikela39, F. Polci8, A. Poluektov49,35, I. Polyakov32, E.
Polycarpo2, A. Popov36, D. Popov11,39, B. Popovici30, C. Potterat2, E.
Price47, J.D. Price53, J. Prisciandaro38, A. Pritchard53, C. Prouve47, V.
Pugatch45, A. Puig Navarro40, G. Punzi24,s, W. Qian56, R. Quagliani7,47, B.
Rachwal27, J.H. Rademacker47, M. Rama24, M. Ramos Pernas38, M.S. Rangel2, I.
Raniuk44, G. Raven43, F. Redi54, S. Reichert55, A.C. dos Reis1, V. Renaudin7,
S. Ricciardi50, S. Richards47, M. Rihl39, K. Rinnert53,39, V. Rives Molina37,
P. Robbe7, A.B. Rodrigues1, E. Rodrigues55, J.A. Rodriguez Lopez64, P.
Rodriguez Perez55, A. Rogozhnikov67, S. Roiser39, V. Romanovsky36, A. Romero
Vidal38, J. W. Ronayne13, M. Rotondo23, T. Ruf39, P. Ruiz Valls68, J.J.
Saborido Silva38, N. Sagidova31, B. Saitta16,f, V. Salustino Guimaraes2, C.
Sanchez Mayordomo68, B. Sanmartin Sedes38, R. Santacesaria26, C. Santamarina
Rios38, M. Santimaria19, E. Santovetti25,l, A. Sarti19,m, C. Satriano26,n, A.
Satta25, D.M. Saunders47, D. Savrina32,33, S. Schael9, M. Schiller39, H.
Schindler39, M. Schlupp10, M. Schmelling11, T. Schmelzer10, B. Schmidt39, O.
Schneider40, A. Schopper39, M. Schubiger40, M.-H. Schune7, R. Schwemmer39, B.
Sciascia19, A. Sciubba26,m, A. Semennikov32, A. Sergi46, N. Serra41, J.
Serrano6, L. Sestini23, P. Seyfert21, M. Shapkin36, I. Shapoval17,44,g, Y.
Shcheglov31, T. Shears53, L. Shekhtman35, V. Shevchenko66, A. Shires10, B.G.
Siddi17, R. Silva Coutinho41, L. Silva de Oliveira2, G. Simi23,s, M.
Sirendi48, N. Skidmore47, T. Skwarnicki60, E. Smith54, I.T. Smith51, J.
Smith48, M. Smith55, H. Snoek42, M.D. Sokoloff58, F.J.P. Soler52, F. Soomro40,
D. Souza47, B. Souza De Paula2, B. Spaan10, P. Spradlin52, S. Sridharan39, F.
Stagni39, M. Stahl12, S. Stahl39, S. Stefkova54, O. Steinkamp41, O.
Stenyakin36, S. Stevenson56, S. Stoica30, S. Stone60, B. Storaci41, S.
Stracka24,t, M. Straticiuc30, U. Straumann41, L. Sun58, W. Sutcliffe54, K.
Swientek28, S. Swientek10, V. Syropoulos43, M. Szczekowski29, T. Szumlak28, S.
T’Jampens4, A. Tayduganov6, T. Tekampe10, G. Tellarini17,g, F. Teubert39, C.
Thomas56, E. Thomas39, J. van Tilburg42, V. Tisserand4, M. Tobin40, S. Tolk43,
L. Tomassetti17,g, D. Tonelli39, S. Topp-Joergensen56, E. Tournefier4, S.
Tourneur40, K. Trabelsi40, M. Traill52, M.T. Tran40, M. Tresch41, A.
Trisovic39, A. Tsaregorodtsev6, P. Tsopelas42, N. Tuning42,39, A. Ukleja29, A.
Ustyuzhanin67,66, U. Uwer12, C. Vacca16,39,f, V. Vagnoni15, S. Valat39, G.
Valenti15, A. Vallier7, R. Vazquez Gomez19, P. Vazquez Regueiro38, C. Vázquez
Sierra38, S. Vecchi17, M. van Veghel42, J.J. Velthuis47, M. Veltri18,h, G.
Veneziano40, M. Vesterinen12, B. Viaud7, D. Vieira2, M. Vieites Diaz38, X.
Vilasis-Cardona37,p, V. Volkov33, A. Vollhardt41, D. Voong47, A. Vorobyev31,
V. Vorobyev35, C. Voß65, J.A. de Vries42, R. Waldi65, C. Wallace49, R.
Wallace13, J. Walsh24, J. Wang60, D.R. Ward48, N.K. Watson46, D. Websdale54,
A. Weiden41, M. Whitehead39, J. Wicht49, G. Wilkinson56,39, M. Wilkinson60, M.
Williams39, M.P. Williams46, M. Williams57, T. Williams46, F.F. Wilson50, J.
Wimberley59, J. Wishahi10, W. Wislicki29, M. Witek27, G. Wormser7, S.A.
Wotton48, K. Wraight52, S. Wright48, K. Wyllie39, Y. Xie63, Z. Xu40, Z. Yang3,
H. Yin63, J. Yu63, X. Yuan35, O. Yushchenko36, M. Zangoli15, M.
Zavertyaev11,c, L. Zhang3, Y. Zhang3, A. Zhelezov12, Y. Zheng62, A. Zhokhov32,
L. Zhong3, V. Zhukov9, S. Zucchelli15.

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

aUniversidade Federal do Triângulo Mineiro (UFTM), Uberaba-MG, Brazil
bLaboratoire Leprince-Ringuet, Palaiseau, France
cP.N. Lebedev Physical Institute, Russian Academy of Science (LPI RAS),
Moscow, Russia
dUniversità di Bari, Bari, Italy
eUniversità di Bologna, Bologna, Italy
fUniversità di Cagliari, Cagliari, Italy
gUniversità di Ferrara, Ferrara, 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
†Deceased