EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 CERN-EP-2017-067 LHCb-PAPER-2017-012 25 July 2017

First observation

of a baryonic ${B}^{0}_{s}$ decay

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

We report the first observation of a baryonic ${B}^{0}_{s}$ decay,
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$, using
proton-proton collision data recorded by the LHCb experiment at center-of-mass
energies of 7 and 8$\mathrm{\,Te\kern-1.00006ptV}$, corresponding to an
integrated luminosity of 3.0$\mbox{\,fb}^{-1}$. The branching fraction is
measured to be ${\mathcal{B}}({{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{K}^{-}})+{\mathcal{B}}({{B}^{0}_{s}}\\!\rightarrow{\overline{p}}{\mathchar
28931\relax}{{K}^{+}})=\left[\vphantom{\sum}{5.46}\pm{0.61}\pm 0.57\pm
0.50({\mathcal{B}})\pm 0.32(f_{s}/f_{d})\right]\times 10^{-6},$ where the
first uncertainty is statistical and the second systematic, the third
uncertainty accounts for the experimental uncertainty on the branching
fraction of the ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ decay
used for normalization, and the fourth uncertainty relates to the knowledge of
the ratio of $b$-quark hadronization probabilities $f_{s}/f_{d}$.

Published in Phys. Rev. Lett. 119 (2017) 041802

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

The experimental study of $B$-meson decays to baryonic final states has a long
history, starting with the first observation of baryonic $B$ decays by the
CLEO collaboration in 1997 [1]. The asymmetric $e^{+}e^{-}$ collider
experiments BaBar and Belle reported numerous searches and observations of
decays of ${B}^{0}$ and ${B}^{+}$ mesons to baryonic final states [2]. The
LHCb collaboration published the first observation of a baryonic ${B}_{c}^{+}$
decay in 2014 [3]. Until now, no baryonic ${B}^{0}_{s}$ decay has ever been
observed with a significance in excess of five standard deviations; the Belle
collaboration provided the only evidence for such a process in the study of
${{B}^{0}_{s}}\\!\rightarrow{{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{}^{-}_{c}}{\mathchar
28931\relax}{{\pi}^{+}}$ decays, with a significance of 4.4 standard
deviations [4].

Areas of particular interest in baryonic $B$ decays are the study of the
hierarchy of branching fractions and the threshold enhancement in the baryon-
antibaryon mass spectrum [5, 2]. Multi-body baryonic $B$ decays are expected
to have higher branching fractions than two-body decays [6, 7]. The
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ and ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ branching
fractions are predicted to be of the order of $10^{-6}$ [8]. The notation
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ is used
hereafter for the sum of both accessible final states
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ and
${{B}^{0}_{s}}\\!\rightarrow{\overline{p}}{\mathchar 28931\relax}{{K}^{+}}$.
As emphasized in Ref. [8], which studied the decays
${{B}^{0}_{s}}\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}h^{-}$, the decay ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ is a
unique baryonic $B$ decay in that it is the only presently known decay where
all four processes, namely the decays of a ${B}^{0}_{s}$ or a ${\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}$ meson to either the
${p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$
or the ${\overline{p}}{\mathchar 28931\relax}{{K}^{+}}$ final state, can
occur. A $B$-flavor-tagged decay-time-dependent study is required in order to
separate the two possible final states and measure their individual branching
fractions as well as $C\\!P$ violation observables.

The current experimental knowledge on the family of
$\mbox{${B}^{0}_{({s})}$}\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}h^{-}$ decays
($h={\pi},{K}$) and related modes such as
$\mbox{${B}^{0}_{({s})}$}\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}h^{-}$, with ${{\overline{\mathchar
28934\relax}}{}^{0}}\rightarrow{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}\gamma$, is rather
scarce. The ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ decay
has been studied by the BaBar [9] and Belle [10, 11] collaborations and the
Belle collaboration has reported the 90% confidence level upper limits
${\mathcal{B}}({{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}})<8.2\times
10^{-7}$ and ${\mathcal{B}}({{B}^{0}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{\pi}^{-}})<3.8\times 10^{-6}$ [10].

Manifestations of $C\\!P$ and $T$ violation in baryonic $B$ decays have been
studied from a theoretical viewpoint, see for example Ref. [12] and references
therein. A large $C\\!P$-violation asymmetry of order 10% is expected for the
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ decay mode [12], which further motivates the
experimental study of $\mbox{${B}^{0}_{({s})}$}\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}h^{-}$ decays.

This Letter presents the first observation of a charmless baryonic
${B}^{0}_{s}$ decay. The branching fraction of the
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ decay is
measured relative to that of the topologically identical
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ decay to suppress common systematic uncertainties:

${\mathcal{B}}({{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{K}^{-}})+{\mathcal{B}}({{B}^{0}_{s}}\\!\rightarrow{\overline{p}}{\mathchar
28931\relax}{{K}^{+}})=\frac{f_{d}}{f_{s}}\,\frac{N({{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{K}^{-}})}{N({{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}})}\,\frac{\epsilon_{{{B}^{0}}\\!\rightarrow{p}{\kern
0.70004pt\overline{\kern-0.70004pt\mathchar
28931\relax}}{{\pi}^{-}}}}{\epsilon_{{{B}^{0}_{s}}\\!\rightarrow{p}{\kern
0.70004pt\overline{\kern-0.70004pt\mathchar
28931\relax}}{{K}^{-}}}}\,{\mathcal{B}}({{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}})\,,$ (1)

where $N$ represents yields determined from mass fits, $f_{q}$ stands for the
${b}$ hadronization probability to the meson $B_{q}$, and $\epsilon$
represents the selection efficiencies. The inclusion of charge-conjugate
processes is implied, unless otherwise stated.

The data sample analyzed corresponds to an integrated luminosity of
1.0$\mbox{\,fb}^{-1}$ of proton-proton collision data collected by the LHCb
experiment at center-of-mass energies of 7$\mathrm{\,Te\kern-1.00006ptV}$ in
2011 and 2.0$\mbox{\,fb}^{-1}$ at 8$\mathrm{\,Te\kern-1.00006ptV}$ in 2012.
The LHCb detector is a single-arm forward spectrometer covering the
pseudorapidity range $2<\eta<5$, designed for the study of particles
containing $b$ or $c$ quarks [13, 14]. The pseudorapidity is defined as
$\eta=-{\rm ln}\left[\tan(\theta/2)\right]$, where $\theta$ is the polar angle
with respect to the proton in the positive $z$ direction. The detector
elements that are particularly relevant to this analysis are a silicon-strip
vertex detector surrounding the proton-proton interaction region that allows
heavy hadrons to be identified from their characteristically long flight
distance; a tracking system that provides a measurement of momentum, $p$, of
charged particles; two ring-imaging Cherenkov detectors that are able to
discriminate between different species of charged hadrons; a calorimeter
system for the measurement of photons and neutral hadrons; and multiwire
proportional chambers for the detection of muons. Simulated data samples,
produced as described in Refs. [15, *Sjostrand:2006za, 17, 18, 19, 20,
*Agostinelli:2002hh, 22], are used to evaluate the response of the detector
and to investigate and characterize possible sources of background.

Events are selected in a similar way for both the signal decay
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ and the
normalization channel ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$, where
${\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}\rightarrow{\overline{p}}{{\pi}^{+}}$. Real-time event selection
is performed by a trigger [23] consisting of a hardware stage, based on
information from the calorimeter and muon systems, followed by a software
stage, which performs a full event reconstruction. The hardware trigger stage
requires events to have a muon with high transverse momentum,
$p_{\mathrm{T}}$, or a hadron, photon or electron with high transverse energy
deposited in the calorimeters. For this analysis, the hardware trigger
decision can either be made on the signal candidates or on other particles in
the event. The software trigger requires a two- or three-track secondary
vertex with a significant displacement from all the primary $pp$ interaction
vertices (PVs). At least one charged particle must have high $p_{\mathrm{T}}$
and be inconsistent with originating from a PV. A multivariate algorithm [24]
is used for the identification of secondary vertices consistent with the decay
of a $b$ or $c$ hadron.

The $\mathchar 28931\relax$ decays are reconstructed in two different
categories: the first consists of $\mathchar 28931\relax$ baryons that decay
early enough for the proton and pion to be reconstructed in the vertex
detector, while the second contains those that decay later such that track
segments cannot be reconstructed in the vertex detector. These reconstruction
categories are referred to as _long_ and _downstream_ , respectively.

The selection of ${B}^{0}_{({s})}$ candidates, formed by combining a
$\mathchar 28931\relax$ candidate with a proton and a pion or kaon, is carried
out with a filtering stage, a requirement on the response of a multilayer
perceptron [25] (MLP) classifier, and particle identification (PID) criteria
discussed below. The proton and pion or kaon, of opposite charge, both decay
products of the $B$ meson, are hereafter referred to as the charged hadrons.
Unless stated otherwise, the terms proton and pion refer to the charged
hadrons from the $B$-meson decay, not to the $\mathchar 28931\relax$ decay
products. Both the ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ and the
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ decay
chains are refitted [26] employing a mass constraint on the $\mathchar
28931\relax$ candidates.

In the filtering stage the $\mathchar 28931\relax$ decay products are required
to have a minimum momentum, $p$, form a good quality vertex and satisfy
$|m({p}{{\pi}^{-}})-m_{\mathchar
28931\relax}|<20(15){\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ for downstream
(long) candidates, where $m_{\mathchar 28931\relax}$ is the $\mathchar
28931\relax$ mass [27]. They must have a large impact parameter (IP) with
respect to all PVs, where the IP is defined as the minimum distance of a track
to a PV. A minimum $\chi^{2}_{\text{IP}}$ with respect to any PV is imposed on
each $\mathchar 28931\relax$ decay product, where $\chi^{2}_{\text{IP}}$ is
defined as the difference between the vertex-fit $\chi^{2}$ of a PV
reconstructed with and without the particle in question. A loose PID
requirement, based primarily on information from the ring-imaging Cherenkov
detectors, is imposed to select the proton candidate from the $\mathchar
28931\relax$ baryon to remove background from
${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ decays. For downstream $\mathchar
28931\relax$ candidates a minimum momentum is also required.

A minimum requirement is imposed on the scalar sum of the $p_{\mathrm{T}}$ of
the $\mathchar 28931\relax$ candidate and the two charged hadrons. The
distance of closest approach among any pair from (p, $\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}$, $h^{-}$) divided by
its uncertainty must be small. The $B$ candidate must have a good quality
vertex, have a minimum $p_{\mathrm{T}}$ and a small $\chi^{2}_{\text{IP}}$
with respect to the associated PV as its reconstructed momentum vector should
point to its production vertex; the associated PV is the one with which it
forms the smallest $\chi^{2}_{\text{IP}}$. The pointing condition of the $B$
candidate is further reinforced by requiring that the angle between the
$B$-candidate momentum vector and the line connecting the associated PV and
the $B$-decay vertex ($B$ direction angle, $\theta_{B}$) is close to zero.

Backgrounds from the ${{B}^{0}}\\!\rightarrow{{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{}^{-}_{c}}{p}$ decay
with ${{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{}^{-}_{c}}\\!\rightarrow{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$
(${{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{}^{-}_{c}}\\!\rightarrow{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$) are
removed from the ${p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ (${p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$) samples
with a veto around the ${\mathchar 28931\relax}^{+}_{c}$ mass [27] of three
times the ${\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ (${\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{K}^{-}}$) invariant mass resolution of approximately
6${\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$. No veto is found to be necessary
to suppress backgrounds from $B$ decays to charmonia and a
${{\pi}^{+}}{{\pi}^{-}}$ pair final states.

Further separation between signal and combinatorial background candidates
relies on MLPs implemented with the TMVA toolkit [28]. The MLPs are trained
using simulated ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ samples,
generated according to a constant matrix element without intermediate
resonances, to represent the signal, and with data from the high-mass sideband
region $5400<m({p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}})<5600{\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ for
the background, to avoid partially reconstructed backgrounds. Separate MLPs
are trained and optimized for each year of data taking and for the two
$\mathchar 28931\relax$ reconstruction categories. Each MLP is used to select
both ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ and
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$
candidates.

The seventeen variables used in the MLP classifiers are properties of the $B$
candidate, the charged hadrons and the $\mathchar 28931\relax$ decay products.
The input variables are the following: the $\chi^{2}$ per degree of freedom of
the kinematic fit of the decay chain [26]; the IP for all particles calculated
with respect to the associated PV; the distance of closest approach between
the two charged hadrons and the sum of their corresponding
$\chi^{2}_{\text{IP}}$; the $\mathchar 28931\relax$ candidate decay-length
significance with respect to the $B$ vertex, i.e. the decay length divided by
its uncertainty; the angle between the $\mathchar 28931\relax$ momentum and
the spacial vector connecting the $B$ and $\mathchar 28931\relax$ decay
vertices in the $B$ rest frame; the $\mathchar 28931\relax$ decay time; the
$B$-meson $p_{\mathrm{T}}$, pseudorapidity, direction angle $\theta_{B}$,
decay-length significance and decay time; the $\mathchar 28931\relax$ helicity
angle defined by the $\mathchar 28931\relax$ momentum in the $B$ rest frame
and the boost axis of the $B$ meson, which is given by the $B$-meson momentum
in the laboratory frame; the pointing variable defined as $P=[\sum_{p,{\kern
0.70004pt\overline{\kern-0.70004pt\mathchar
28931\relax}},h^{-}}{p}\times\sin{\theta_{B}}]/[\sum_{p,{\kern
0.70004pt\overline{\kern-0.70004pt\mathchar
28931\relax}},h^{-}}{p}\times\sin{\theta_{B}}+\sum_{p,{\kern
0.70004pt\overline{\kern-0.70004pt\mathchar
28931\relax}},h^{-}}{\mbox{$p_{\mathrm{T}}$}}]$. The optimal MLP requirement
for each of the four subsamples is determined by maximizing the signal
significance of the ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$
normalization decay, with the variation of the signal efficiency with MLP cut
value determined from simulation.

A PID selection is applied to the charged hadrons after the MLP selection. No
additional PID requirement is applied to the proton from the $\mathchar
28931\relax$ candidate since no contamination from misidentified
${{K}^{0}_{\mathrm{\scriptscriptstyle
S}}}\\!\rightarrow{{\pi}^{+}}{{\pi}^{-}}$ decays is observed. The optimization
of the PID requirements follows the same procedure as the optimization of the
MLP selection. If more than one candidate is selected in any event of any
subsample, which occurs in about 5% of selected events, one is chosen at
random.

Large data control samples of ${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}$,
$\mathchar 28931\relax\\!\rightarrow{p}{{\pi}^{-}}$ and ${{\mathchar
28931\relax}^{+}_{c}}\\!\rightarrow{p}{{K}^{-}}{{\pi}^{+}}$ decays are
employed [29] to determine the efficiency of the PID requirements. All other
selection efficiencies are determined from simulation. It is necessary to
account for the distribution of signal candidates and the variation of the
efficiency over the phase space of the decay. The variation is well described
by the factorized efficiencies in the two-dimensional space of the variables
$m^{2}(p{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}})$ and
$m^{2}(ph^{-})$ defining the Dalitz plot. Simulated events are binned in
$m^{2}(p{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}})$ in
order to determine the selection efficiencies, the variation in
$m^{2}(ph^{-})$ being mild and therefore integrated out. The distribution of
signal decays in the phase space is obtained separately for each spectrum with
the sPlot technique [30] with the $B$-meson candidate invariant mass used as
the discriminating variable. The overall efficiencies of this analysis are of
order $10^{-4}$.

The efficiency of the software trigger selection on both decay modes varied
during the data-taking period. During the 2011 data taking, downstream tracks
were not reconstructed in the software trigger. Such tracks were included in
the trigger during the 2012 data taking and a further significant improvement
in the algorithms was implemented mid-year. The corresponding changes to the
trigger efficiency are taken into account.

Potential sources of background to the ${p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}h^{-}$ spectra are
investigated using simulation samples. Cross-feed between the
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ and ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ decay
modes is the dominant source of peaking background. The loop-mediated decays
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{K}^{-}}$ and ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ are
suppressed and estimated to be insignificant [8]. Pion-kaon misidentification
from $b$-baryon decays such as the recently observed decays ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}h^{+}h^{\prime-}$
[31] is found to be negligible. The influence of proton-pion misidentification
in the reconstruction and selection of the $\mathchar 28931\relax$ baryon
arising from ${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ cross-feed is checked
since the PID requirement on the proton from the $\mathchar 28931\relax$ is
rather loose. It is verified with Armenteros-Podolanski plots [32] that the
${K}^{0}_{\mathrm{\scriptscriptstyle S}}$ contamination can be ignored. Cross-
feed from the presently unobserved decay ${{\mathchar
28931\relax}^{0}_{b}}\\!\rightarrow{\mathchar 28931\relax}{p}{\overline{p}}$
due to proton-pion and proton-kaon misidentification is assumed to be
negligible considering that the proton misidentification rate is small.
Partially reconstructed decays such as the unobserved
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\rho}^{-}}$ and ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{*-}}$ modes are
treated as a source of systematic uncertainty. Decay modes containing a
${\overline{\mathchar 28934\relax}}{}^{0}$ baryon decaying into
${{\overline{\mathchar 28934\relax}}{}^{0}}\rightarrow{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}\gamma$, where the
$\gamma$ is not detected, can pollute the signal regions due to the small mass
difference $m({{\mathchar 28934\relax}^{0}})-m({\mathchar 28931\relax})\approx
77{\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ [27]. The decay
${{B}^{0}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{\pi}^{-}}$ is expected to have a branching fraction at
the level of $10^{-6}$ [33], though searches for the
$\mbox{${B}^{0}_{({s})}$}\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}h^{-}$ family of decays have found no signal [10]. The
decays ${{B}^{0}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{\pi}^{-}}$ and
${{B}^{0}_{s}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{K}^{-}}$ are expected to be the dominant members of the
family and are included in the fits to the data.

The yields of the signal and background candidates in eight subsamples are
determined from a simultaneous unbinned extended maximum likelihood fit to the
${p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}h^{-}$
invariant mass distributions. The eight subsamples correspond to the 2011 and
2012 data-taking periods, the two $\mathchar 28931\relax$ reconstruction
categories, and the ${p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ and ${p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ final
state hypotheses. This approach allows the use of common shape parameters, and
the level of cross-feed background can be better constrained by fitting all
subsamples simultaneously. The probability density function in each subsample
is defined as the sum of components accounting for the signal decay, the
cross-feed contribution, the ${{B}^{0}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{\pi}^{-}}$ and
${{B}^{0}_{s}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{K}^{-}}$ decays, and combinatorial background.

The signal and normalization modes are modeled with the sum of two Novosibirsk
functions [34]. All shape parameters are fixed to the values obtained
separately for each subsample from simulation samples. The
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ and ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ peak
positions are free parameters determined simultaneously in all subsamples. The
cross-feed ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$
(${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$) in the ${p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$
(${p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{K}^{-}}$) invariant mass distribution is modeled with the sum
of a Gaussian and a modified Fermi function defined as the product of an
exponential and a Fermi-Dirac function. The
${{B}^{0}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{\pi}^{-}}$ and
${{B}^{0}_{s}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{K}^{-}}$ decays are modeled differently according to the
$\mathchar 28931\relax$ reconstruction category and the ${p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}h^{-}$ invariant mass
hypothesis under which they are reconstructed. Depending on the category a
modified Fermi function, a sum of two Novosibirsk functions, the sum of a
Novosibirsk and a Gaussian function, or the sum of a Novosibirsk and a
modified Fermi function are used. A combinatorial background component
described by an exponential function is present for both ${p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}h^{-}$ final states.

The yields of the ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$
candidates are determined in the fit together with the ratio of the
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ to
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ branching fractions, which is determined
simultaneously across all subsamples accounting for differences in selection
efficiencies. These depend on the data-taking period, $\mathchar 28931\relax$
reconstruction category and mass hypothesis of the meson from the $B$ decay.
The uncertainties arising from the ratios of efficiencies are included in the
fit as Gaussian constraints. The yields of the
${{B}^{0}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{\pi}^{-}}$ and
${{B}^{0}_{s}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{K}^{-}}$ decays are defined relative to those of the
corresponding ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ and
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ decays,
respectively. These two ${{\mathchar 28934\relax}^{0}}$-to-${\mathchar
28931\relax}$ decay yield ratios are determined simultaneously in the fit
across all subsamples following the same procedure as for the
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ decay. The
combinatorial background yield and shape parameters are treated independently
in each subsample and are allowed to vary in the fit.

Figure 1 presents the fit to the ${p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}h^{-}$ invariant mass
distributions for all subsamples combined. Both
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ and ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ signals
are prominent. In particular, the ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ decay is
observed with a statistical significance above 15 standard deviations,
estimated from the change in log-likelihood between fits with and without the
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ signal
component [35]. It constitutes the first observation of a baryonic
${B}^{0}_{s}$ decay. The yields summed over all subsamples are
$N({{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}})=519\pm
28$ and $N({{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}})=234\pm
29$, where the uncertainties are statistical only.

Figure 1: Mass distributions for $b$-hadron candidates for (left) the
${p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ and (right) the ${p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ sample for
the combined long and downstream categories. The black points represent the
data, the solid blue curve the result of the fit, the red dashed curve the
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$
contribution, the black (magenta) dotted curve the
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$
(${{B}^{0}_{s}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{K}^{-}}$) and the green dash-dotted curve the
contribution from ${{B}^{0}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{\pi}^{-}}$ decays. The combinatorial background
distribution is indicated by the shaded histogram.

The sPlot technique is used to subtract the background and obtain the phase
space distribution of signal candidates. Figure 2 shows the $m({p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}})$ invariant mass
distributions for the ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ and
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ candidates
after correcting for the distribution selection efficiencies. Both
distributions show a pronounced enhancement at threshold in the baryon-
antibaryon invariant mass, first suggested in Ref. [5] and observed in several
baryonic $B$ decay modes.

Figure 2: Efficiency-corrected and background-subtracted $m({p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}})$ invariant mass
distributions for (left) ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ and
(right) ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$
candidates. The distributions are normalized to unity.

The sources of systematic uncertainty arise from the fit model, the knowledge
of the selection efficiencies, and the uncertainties on the
${{B}^{0}}\\!\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{\pi}^{-}}$ branching fraction and on the ratio of hadronization
probabilities $f_{s}/f_{d}$. Uncertainties on the selection efficiencies arise
from residual differences between data and simulation in the trigger,
reconstruction, selection and particle identification. Additional
uncertainties arise due to the limited size of the simulation samples and the
corresponding uncertainty on the distribution of the efficiencies across the
decay phase space. As the efficiencies depend on the signal decay-time
distribution, the effect coming from the different lifetimes of the
${B}^{0}_{s}$ mass eigenstates has been evaluated [36]. Pseudoexperiments are
used to estimate the effect of using alternative shapes for the fit
components, of including additional backgrounds in the fit such as partially
reconstructed decays, and of excluding the
${{B}^{0}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{\pi}^{-}}$ and
${{B}^{0}_{s}}\\!\rightarrow{p}{{\overline{\mathchar
28934\relax}}{}^{0}}{{K}^{-}}$ decays that show no significant contribution.
Intrinsic biases in the fitted signal yields are investigated with ensembles
of simulated pseudoexperiments. A small bias is found and added to the
systematic uncertainty on the fit model. The systematic uncertainty due to the
knowledge of the efficiencies involved in the definition of fit constraints is
negligible. The total systematic uncertainty on the
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ branching
fraction is given by the sum of all uncertainties added in quadrature and
amounts to 10.5%; it is dominated by the systematic uncertainty on the fit
model.

The uncertainty on the branching fraction of the normalization decay,
${\mathcal{B}}({{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}})=(3.14\pm
0.29)\times 10^{-6}$ [27], is taken as a systematic uncertainty from external
inputs. The 5.8% uncertainty on the latest $f_{s}/f_{d}$ combination from
LHCb, $f_{s}/f_{d}=0.259\pm 0.015$ [37], is taken as a second source of
systematic uncertainty from external inputs.

The ${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ branching
fraction, determined relative to that of the ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$
normalization channel according to Eq. 1, is measured to be

${\mathcal{B}}({{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar
28931\relax}}{{K}^{-}})+{\mathcal{B}}({{B}^{0}_{s}}\\!\rightarrow{\overline{p}}{\mathchar
28931\relax}{{K}^{+}})=\left[\vphantom{\sum}{5.46}\pm{0.61}\pm 0.57\pm
0.50({\mathcal{B}})\pm 0.32(f_{s}/f_{d})\right]\times 10^{-6}\;,$

where the first uncertainty is statistical and the second systematic, the
third uncertainty accounts for the experimental uncertainty on the branching
fraction of the ${{B}^{0}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$ decay,
and the fourth uncertainty relates to the knowledge of $f_{s}/f_{d}$.

In summary, the first observation of the three-body charmless baryonic decay
${{B}^{0}_{s}}\\!\rightarrow{p}{\kern
1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{K}^{-}}$ is
reported using a proton-proton collision data sample collected by the LHCb
experiment, corresponding to an integrated luminosity of
3.0$\mbox{\,fb}^{-1}$. The decay is observed with a statistical significance
above 15 standard deviations, which constitutes the first observation of a
baryonic ${B}^{0}_{s}$ decay.

Decays of $B$ mesons to final states containing baryons are now observed for
all $B$-meson species. Their study provides valuable information on the
dynamics of hadronic decays of $B$ mesons. The present analysis motivates
further theoretical studies of baryonic ${B}^{0}_{s}$ decays in addition to
those currently published [8, 38, 39, 6].

## 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); MOST
and NSFC (China); CNRS/IN2P3 (France); BMBF, DFG and MPG (Germany); INFN
(Italy); NWO (The Netherlands); MNiSW and NCN (Poland); MEN/IFA (Romania);
MinES and FASO (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] CLEO collaboration, X. Fu, et al., _Observation of exclusive $B$ decays to final states containing a charmed baryon_, Phys. Rev. Lett. 79 (1997) 3125
  * [2] BaBar and Belle collaborations, A. J. Bevan et al., _The physics of the B factories_ , Eur. Phys. J. C74 (2014) 3026, arXiv:1406.6311
  * [3] LHCb collaboration, R. Aaij et al., _First observation of a baryonic ${{B}_{c}^{+}}$ decay_, Phys. Rev. Lett. 113 (2014) 152003, arXiv:1408.0971
  * [4] Belle collaboration, E. Solovieva et al., _Evidence for ${{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}}\rightarrow{{\mathchar 28931\relax}^{+}_{c}}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$_, Phys. Lett. B 726 (2013) 206, arXiv:1304.6931
  * [5] W.-S. Hou and A. Soni, _Pathways to rare baryonic $B$ decays_, Phys. Rev. Lett. 86 (2001) 4247, arXiv:hep-ph/0008079
  * [6] Y. K. Hsiao and C. Q. Geng, _Violation of partial conservation of the axial-vector current and two-body baryonic $B$ and $D_{s}$ decays_, Phys. Rev. D91 (2015) 077501, arXiv:1407.7639
  * [7] H.-Y. Cheng and C.-K. Chua, _On the smallness of tree-dominated charmless two-body baryonic $B$ decay rates_, Phys. Rev. D91 (2015) 036003, arXiv:1412.8272
  * [8] C. Q. Geng, Y. K. Hsiao, and E. Rodrigues, _Three-body charmless baryonic ${\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}$ decays_, Phys. Lett. B 767 (2017) 205, arXiv:1612.08133
  * [9] BaBar collaboration, B. Aubert et al., _Measurement of the branching fraction and $\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}$ polarization in ${{B}^{0}}\rightarrow{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{p}{{\pi}^{-}}$_, Phys. Rev. D79 (2009) 112009, arXiv:0904.4724
  * [10] Belle collaboration, M. Z. Wang et al., _Observation of ${{B}^{0}}\rightarrow{p}{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$_, Phys. Rev. Lett. 90 (2003) 201802, arXiv:hep-ex/0302024
  * [11] Belle collaboration, M.-Z. Wang et al., _Study of ${{{B}^{+}}}\rightarrow p{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}\gamma,p{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{0}}$ and ${{B}^{0}}\rightarrow p{\kern 1.00006pt\overline{\kern-1.00006pt\mathchar 28931\relax}}{{\pi}^{-}}$_, Phys. Rev. D76 (2007) 052004, arXiv:0704.2672
  * [12] C. Q. Geng and Y. K. Hsiao, _Direct $C\\!P$ and $T$ violation in baryonic $B$ decays_, Int. J. Mod. Phys. A23 (2008) 3290, arXiv:0801.0022
  * [13] LHCb collaboration, A. A. Alves Jr. et al., _The LHCb detector at the LHC_, JINST 3 (2008) S08005
  * [14] LHCb collaboration, R. Aaij et al., _LHCb detector performance_ , Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [15] T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * [16] T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * [17] I. Belyaev et al., _Handling of the generation of primary events in Gauss, the LHCb simulation framework_ , J. Phys. Conf. Ser. 331 (2011) 032047
  * [18] D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * [19] 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
  * [20] Geant4 collaboration, J. Allison et al., _Geant4 developments and applications_ , IEEE Trans. Nucl. Sci. 53 (2006) 270
  * [21] Geant4 collaboration, S. Agostinelli et al., _Geant4: A simulation toolkit_ , Nucl. Instrum. Meth. A506 (2003) 250
  * [22] M. Clemencic et al., _The LHCb simulation application, Gauss: Design, evolution and experience_, J. Phys. Conf. Ser. 331 (2011) 032023
  * [23] R. Aaij et al., _The LHCb trigger and its performance in 2011_, JINST 8 (2013) P04022, arXiv:1211.3055
  * [24] V. V. Gligorov and M. Williams, _Efficient, reliable and fast high-level triggering using a bonsai boosted decision tree_ , JINST 8 (2013) P02013, arXiv:1210.6861
  * [25] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, Parallel distributed processing: explorations in the microstructure of cognition, vol. 1, MIT, Cambridge, USA, 1986
  * [26] W. D. Hulsbergen, _Decay chain fitting with a Kalman filter_ , Nucl. Instrum. Meth. A552 (2005) 566, arXiv:physics/0503191
  * [27] Particle Data Group, C. Patrignani et al., _Review of particle physics_ , Chin. Phys. C40 (2016) 100001
  * [28] P. Speckmayer, A. Hoecker, J. Stelzer, and H. Voss, _The toolkit for multivariate data analysis: TMVA 4_ , J. Phys. Conf. Ser. 219 (2010) 032057
  * [29] M. Adinolfi et al., _Performance of the LHCb RICH detector at the LHC_, Eur. Phys. J. C73 (2013) 2431, arXiv:1211.6759
  * [30] M. Pivk and F. R. Le Diberder, _sPlot: A statistical tool to unfold data distributions_ , Nucl. Instrum. Meth. A555 (2005) 356, arXiv:physics/0402083
  * [31] 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 ${{\mathchar 28932\relax}^{0}_{b}}$ decays to ${\mathchar 28931\relax}h^{+}h^{-}$ final states_, JHEP 05 (2016) 081, arXiv:1603.00413
  * [32] J. Podolanski and R. Armenteros, _III. Analysis of V-events_ , The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 45 (1954) 13
  * [33] C.-K. Chua, W.-S. Hou, and S.-Y. Tsai, _Charmless three-body baryonic B decays_ , Phys. Rev. D66 (2002) 054004, arXiv:hep-ph/0204185
  * [34] BaBar collaboration, B. Aubert et al., _Search for decays of ${B}^{0}$ mesons into pairs of leptons_, in Proceedings of the $31^{st}$ international conference on high energy physics, ICHEP 2002, Amsterdam, The Netherlands, July 25-31, 2002, 2002.  arXiv:hep-ex/0207083
  * [35] S. S. Wilks, _The large-sample distribution of the likelihood ratio for testing composite hypotheses_ , Ann. Math. Stat. 9 (1938) 60
  * [36] K. De Bruyn et al., _Branching ratio measurements of ${B}^{0}_{s}$ decays_, Phys. Rev. D86 (2012) 014027, arXiv:1204.1735
  * [37] LHCb collaboration, R. Aaij et al., _Measurement of the fragmentation fraction ratio $f_{s}/f_{d}$ and its dependence on $B$ meson kinematics_, JHEP 04 (2013) 001, arXiv:1301.5286, $f_{s}/f_{d}$ value updated in LHCb-CONF-2013-011
  * [38] Y. K. Hsiao and C. Q. Geng, _$f_{J}(2220)$ and hadronic ${\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}$ decays_, Eur. Phys. J. C75 (2015) 101, arXiv:1412.4900
  * [39] C.-K. Chua, _Charmless two-body baryonic ${B}_{u,d,s}$ decays revisited_, Phys. Rev. D89 (2014) 056003, arXiv:1312.2335

LHCb collaboration

R. Aaij40, B. Adeva39, M. Adinolfi48, Z. Ajaltouni5, S. Akar59, J. Albrecht10,
F. Alessio40, M. Alexander53, S. Ali43, G. Alkhazov31, P. Alvarez Cartelle55,
A.A. Alves Jr59, S. Amato2, S. Amerio23, Y. Amhis7, L. An3, L. Anderlini18, G.
Andreassi41, M. Andreotti17,g, J.E. Andrews60, R.B. Appleby56, F. Archilli43,
P. d’Argent12, J. Arnau Romeu6, A. Artamonov37, M. Artuso61, E. Aslanides6, G.
Auriemma26, M. Baalouch5, I. Babuschkin56, S. Bachmann12, J.J. Back50, A.
Badalov38, C. Baesso62, S. Baker55, V. Balagura7,c, W. Baldini17, A.
Baranov35, R.J. Barlow56, C. Barschel40, S. Barsuk7, W. Barter56, F.
Baryshnikov32, M. Baszczyk27,l, V. Batozskaya29, V. Battista41, A. Bay41, L.
Beaucourt4, J. Beddow53, F. Bedeschi24, I. Bediaga1, A. Beiter61, L.J. Bel43,
V. Bellee41, N. Belloli21,i, K. Belous37, I. Belyaev32, E. Ben-Haim8, G.
Bencivenni19, S. Benson43, S. Beranek9, A. Berezhnoy33, R. Bernet42, A.
Bertolin23, C. Betancourt42, F. Betti15, M.-O. Bettler40, M. van Beuzekom43,
Ia. Bezshyiko42, S. Bifani47, P. Billoir8, A. Birnkraut10, A. Bitadze56, A.
Bizzeti18,u, T. Blake50, F. Blanc41, J. Blouw11,†, S. Blusk61, V. Bocci26, T.
Boettcher58, A. Bondar36,w, N. Bondar31, W. Bonivento16, I. Bordyuzhin32, A.
Borgheresi21,i, S. Borghi56, M. Borisyak35, M. Borsato39, F. Bossu7, M.
Boubdir9, T.J.V. Bowcock54, E. Bowen42, C. Bozzi17,40, S. Braun12, T.
Britton61, J. Brodzicka56, E. Buchanan48, C. Burr56, A. Bursche16, J.
Buytaert40, S. Cadeddu16, R. Calabrese17,g, M. Calvi21,i, M. Calvo Gomez38,m,
A. Camboni38, P. Campana19, D.H. Campora Perez40, L. Capriotti56, A.
Carbone15,e, G. Carboni25,j, R. Cardinale20,h, A. Cardini16, P. Carniti21,i,
L. Carson52, K. Carvalho Akiba2, G. Casse54, L. Cassina21,i, L. Castillo
Garcia41, M. Cattaneo40, G. Cavallero20,40,h, R. Cenci24,t, D. Chamont7, M.
Charles8, Ph. Charpentier40, G. Chatzikonstantinidis47, M. Chefdeville4, S.
Chen56, S.F. Cheung57, V. Chobanova39, M. Chrzaszcz42,27, A. Chubykin31, X.
Cid Vidal39, G. Ciezarek43, P.E.L. Clarke52, M. Clemencic40, H.V. Cliff49, J.
Closier40, V. Coco59, J. Cogan6, E. Cogneras5, V. Cogoni16,f, L. Cojocariu30,
P. Collins40, A. Comerma-Montells12, A. Contu40, A. Cook48, G. Coombs40, S.
Coquereau38, G. Corti40, M. Corvo17,g, C.M. Costa Sobral50, B. Couturier40,
G.A. Cowan52, D.C. Craik52, A. Crocombe50, M. Cruz Torres62, S. Cunliffe55, R.
Currie52, C. D’Ambrosio40, F. Da Cunha Marinho2, E. Dall’Occo43, J. Dalseno48,
A. Davis3, O. De Aguiar Francisco54, K. De Bruyn6, S. De Capua56, M. De
Cian12, J.M. De Miranda1, L. De Paula2, M. De Serio14,d, P. De Simone19, C.T.
Dean53, D. Decamp4, M. Deckenhoff10, L. Del Buono8, H.-P. Dembinski11, M.
Demmer10, A. Dendek28, D. Derkach35, O. Deschamps5, F. Dettori54, B. Dey22, A.
Di Canto40, P. Di Nezza19, H. Dijkstra40, F. Dordei40, M. Dorigo41, A. Dosil
Suárez39, A. Dovbnya45, K. Dreimanis54, L. Dufour43, G. Dujany56, K. Dungs40,
P. Durante40, R. Dzhelyadin37, M. Dziewiecki12, A. Dziurda40, A. Dzyuba31, N.
Déléage4, S. Easo51, M. Ebert52, U. Egede55, V. Egorychev32, S. Eidelman36,w,
S. Eisenhardt52, U. Eitschberger10, R. Ekelhof10, L. Eklund53, S. Ely61, S.
Esen12, H.M. Evans49, T. Evans57, A. Falabella15, N. Farley47, S. Farry54, R.
Fay54, D. Fazzini21,i, D. Ferguson52, G. Fernandez38, A. Fernandez Prieto39,
F. Ferrari15, F. Ferreira Rodrigues2, M. Ferro-Luzzi40, S. Filippov34, R.A.
Fini14, M. Fiore17,g, M. Fiorini17,g, M. Firlej28, C. Fitzpatrick41, T.
Fiutowski28, F. Fleuret7,b, K. Fohl40, M. Fontana16,40, F. Fontanelli20,h,
D.C. Forshaw61, R. Forty40, V. Franco Lima54, M. Frank40, C. Frei40, J.
Fu22,q, W. Funk40, E. Furfaro25,j, C. Färber40, E. Gabriel52, A. Gallas
Torreira39, D. Galli15,e, S. Gallorini23, S. Gambetta52, M. Gandelman2, P.
Gandini57, Y. Gao3, L.M. Garcia Martin70, J. García Pardiñas39, J. Garra
Tico49, L. Garrido38, P.J. Garsed49, D. Gascon38, C. Gaspar40, L. Gavardi10,
G. Gazzoni5, D. Gerick12, E. Gersabeck12, M. Gersabeck56, T. Gershon50, Ph.
Ghez4, S. Gianì41, V. Gibson49, O.G. Girard41, L. Giubega30, K. Gizdov52, V.V.
Gligorov8, D. Golubkov32, A. Golutvin55,40, A. Gomes1,a, I.V. Gorelov33, C.
Gotti21,i, E. Govorkova43, R. Graciani Diaz38, L.A. Granado Cardoso40, E.
Graugés38, E. Graverini42, G. Graziani18, A. Grecu30, R. Greim9, P.
Griffith16, L. Grillo21,40,i, L. Gruber40, B.R. Gruberg Cazon57, O.
Grünberg67, E. Gushchin34, Yu. Guz37, T. Gys40, C. Göbel62, T. Hadavizadeh57,
C. Hadjivasiliou5, G. Haefeli41, C. Haen40, S.C. Haines49, B. Hamilton60, X.
Han12, S. Hansmann-Menzemer12, N. Harnew57, S.T. Harnew48, J. Harrison56, M.
Hatch40, J. He63, T. Head41, A. Heister9, K. Hennessy54, P. Henrard5, L.
Henry70, E. van Herwijnen40, M. Heß67, A. Hicheur2, D. Hill57, C. Hombach56,
P.H. Hopchev41, Z.-C. Huard59, W. Hulsbergen43, T. Humair55, M. Hushchyn35, D.
Hutchcroft54, M. Idzik28, P. Ilten58, R. Jacobsson40, J. Jalocha57, E. Jans43,
A. Jawahery60, F. Jiang3, M. John57, D. Johnson40, C.R. Jones49, C. Joram40,
B. Jost40, N. Jurik57, S. Kandybei45, M. Karacson40, J.M. Kariuki48, S.
Karodia53, M. Kecke12, M. Kelsey61, M. Kenzie49, T. Ketel44, E. Khairullin35,
B. Khanji12, C. Khurewathanakul41, T. Kirn9, S. Klaver56, K. Klimaszewski29,
T. Klimkovich11, S. Koliiev46, M. Kolpin12, I. Komarov41, R. Kopecna12, P.
Koppenburg43, A. Kosmyntseva32, S. Kotriakhova31, M. Kozeiha5, L. Kravchuk34,
M. Kreps50, P. Krokovny36,w, F. Kruse10, W. Krzemien29, W. Kucewicz27,l, M.
Kucharczyk27, V. Kudryavtsev36,w, A.K. Kuonen41, K. Kurek29, T.
Kvaratskheliya32,40, D. Lacarrere40, G. Lafferty56, A. Lai16, G. Lanfranchi19,
C. Langenbruch9, T. Latham50, C. Lazzeroni47, R. Le Gac6, J. van Leerdam43, A.
Leflat33,40, J. Lefran$c$cois7, R. Lefèvre5, F. Lemaitre40, E. Lemos Cid39, O.
Leroy6, T. Lesiak27, B. Leverington12, T. Li3, Y. Li7, Z. Li61, T.
Likhomanenko35,68, R. Lindner40, F. Lionetto42, X. Liu3, D. Loh50, I.
Longstaff53, J.H. Lopes2, D. Lucchesi23,o, M. Lucio Martinez39, H. Luo52, A.
Lupato23, E. Luppi17,g, O. Lupton40, A. Lusiani24, X. Lyu63, F. Machefert7, F.
Maciuc30, B. Maddock59, O. Maev31, K. Maguire56, S. Malde57, A. Malinin68, T.
Maltsev36, G. Manca16,f, G. Mancinelli6, P. Manning61, J. Maratas5,v, J.F.
Marchand4, U. Marconi15, C. Marin Benito38, M. Marinangeli41, P. Marino24,t,
J. Marks12, G. Martellotti26, M. Martin6, M. Martinelli41, D. Martinez
Santos39, F. Martinez Vidal70, D. Martins Tostes2, L.M. Massacrier7, A.
Massafferri1, R. Matev40, A. Mathad50, Z. Mathe40, C. Matteuzzi21, A. Mauri42,
E. Maurice7,b, B. Maurin41, A. Mazurov47, M. McCann55,40, A. McNab56, R.
McNulty13, B. Meadows59, F. Meier10, D. Melnychuk29, M. Merk43, A.
Merli22,40,q, E. Michielin23, D.A. Milanes66, M.-N. Minard4, D.S. Mitzel12, A.
Mogini8, J. Molina Rodriguez1, I.A. Monroy66, S. Monteil5, M. Morandin23, M.J.
Morello24,t, O. Morgunova68, J. Moron28, A.B. Morris52, R. Mountain61, F.
Muheim52, M. Mulder43, M. Mussini15, D. Müller56, J. Müller10, K. Müller42, V.
Müller10, P. Naik48, T. Nakada41, R. Nandakumar51, A. Nandi57, I. Nasteva2, M.
Needham52, N. Neri22,40, S. Neubert12, N. Neufeld40, M. Neuner12, T.D.
Nguyen41, C. Nguyen-Mau41,n, S. Nieswand9, R. Niet10, N. Nikitin33, T.
Nikodem12, A. Nogay68, D.P. O’Hanlon50, A. Oblakowska-Mucha28, V. Obraztsov37,
S. Ogilvy19, R. Oldeman16,f, C.J.G. Onderwater71, A. Ossowska27, J.M. Otalora
Goicochea2, P. Owen42, A. Oyanguren70, P.R. Pais41, A. Palano14,d, M.
Palutan19,40, A. Papanestis51, M. Pappagallo14,d, L.L. Pappalardo17,g, C.
Pappenheimer59, W. Parker60, C. Parkes56, G. Passaleva18, A. Pastore14,d, M.
Patel55, C. Patrignani15,e, A. Pearce40, A. Pellegrino43, G. Penso26, M. Pepe
Altarelli40, S. Perazzini40, P. Perret5, L. Pescatore41, K. Petridis48, A.
Petrolini20,h, A. Petrov68, M. Petruzzo22,q, E. Picatoste Olloqui38, B.
Pietrzyk4, M. Pikies27, D. Pinci26, A. Pistone20,h, A. Piucci12, V.
Placinta30, S. Playfer52, M. Plo Casasus39, T. Poikela40, F. Polci8, M. Poli
Lener19, A. Poluektov50,36, I. Polyakov61, E. Polycarpo2, G.J. Pomery48, S.
Ponce40, A. Popov37, D. Popov11,40, B. Popovici30, S. Poslavskii37, C.
Potterat2, E. Price48, J. Prisciandaro39, C. Prouve48, V. Pugatch46, A. Puig
Navarro42, G. Punzi24,p, C. Qian63, W. Qian50, R. Quagliani7,48, B. Rachwal28,
J.H. Rademacker48, M. Rama24, M. Ramos Pernas39, M.S. Rangel2, I. Raniuk45,†,
F. Ratnikov35, G. Raven44, M. Ravonel Salzgeber40, M. Reboud4, F. Redi55, S.
Reichert10, A.C. dos Reis1, C. Remon Alepuz70, V. Renaudin7, S. Ricciardi51,
S. Richards48, M. Rihl40, K. Rinnert54, V. Rives Molina38, P. Robbe7, A.B.
Rodrigues1, E. Rodrigues59, J.A. Rodriguez Lopez66, P. Rodriguez Perez56,†, A.
Rogozhnikov35, S. Roiser40, A. Rollings57, V. Romanovskiy37, A. Romero
Vidal39, J.W. Ronayne13, M. Rotondo19, M.S. Rudolph61, T. Ruf40, P. Ruiz
Valls70, J.J. Saborido Silva39, E. Sadykhov32, N. Sagidova31, B. Saitta16,f,
V. Salustino Guimaraes1, D. Sanchez Gonzalo38, C. Sanchez Mayordomo70, B.
Sanmartin Sedes39, R. Santacesaria26, C. Santamarina Rios39, M. Santimaria19,
E. Santovetti25,j, A. Sarti19,k, C. Satriano26,s, A. Satta25, D.M. Saunders48,
D. Savrina32,33, S. Schael9, M. Schellenberg10, M. Schiller53, H. Schindler40,
M. Schlupp10, M. Schmelling11, T. Schmelzer10, B. Schmidt40, O. Schneider41,
A. Schopper40, H.F. Schreiner59, K. Schubert10, M. Schubiger41, M.-H. Schune7,
R. Schwemmer40, B. Sciascia19, A. Sciubba26,k, A. Semennikov32, A. Sergi47, N.
Serra42, J. Serrano6, L. Sestini23, P. Seyfert21, M. Shapkin37, I. Shapoval45,
Y. Shcheglov31, T. Shears54, L. Shekhtman36,w, V. Shevchenko68, B.G.
Siddi17,40, R. Silva Coutinho42, L. Silva de Oliveira2, G. Simi23,o, S.
Simone14,d, M. Sirendi49, N. Skidmore48, T. Skwarnicki61, E. Smith55, I.T.
Smith52, J. Smith49, M. Smith55, l. Soares Lavra1, M.D. Sokoloff59, F.J.P.
Soler53, B. Souza De Paula2, B. Spaan10, P. Spradlin53, S. Sridharan40, F.
Stagni40, M. Stahl12, S. Stahl40, P. Stefko41, S. Stefkova55, O. Steinkamp42,
S. Stemmle12, O. Stenyakin37, H. Stevens10, S. Stoica30, S. Stone61, B.
Storaci42, S. Stracka24,p, M.E. Stramaglia41, M. Straticiuc30, U. Straumann42,
L. Sun64, W. Sutcliffe55, K. Swientek28, V. Syropoulos44, M. Szczekowski29, T.
Szumlak28, S. T’Jampens4, A. Tayduganov6, T. Tekampe10, G. Tellarini17,g, F.
Teubert40, E. Thomas40, J. van Tilburg43, M.J. Tilley55, V. Tisserand4, M.
Tobin41, S. Tolk49, L. Tomassetti17,g, D. Tonelli24, S. Topp-Joergensen57, F.
Toriello61, R. Tourinho Jadallah Aoude1, E. Tournefier4, S. Tourneur41, K.
Trabelsi41, M. Traill53, M.T. Tran41, M. Tresch42, A. Trisovic40, A.
Tsaregorodtsev6, P. Tsopelas43, A. Tully49, N. Tuning43, A. Ukleja29, A.
Ustyuzhanin35, U. Uwer12, C. Vacca16,f, A. Vagner69, V. Vagnoni15,40, A.
Valassi40, S. Valat40, G. Valenti15, R. Vazquez Gomez19, P. Vazquez
Regueiro39, S. Vecchi17, M. van Veghel43, J.J. Velthuis48, M. Veltri18,r, G.
Veneziano57, A. Venkateswaran61, T.A. Verlage9, M. Vernet5, M. Vesterinen12,
J.V. Viana Barbosa40, B. Viaud7, D. Vieira63, M. Vieites Diaz39, H. Viemann67,
X. Vilasis-Cardona38,m, M. Vitti49, V. Volkov33, A. Vollhardt42, B. Voneki40,
A. Vorobyev31, V. Vorobyev36,w, C. Voß9, J.A. de Vries43, C. Vázquez Sierra39,
R. Waldi67, C. Wallace50, R. Wallace13, J. Walsh24, J. Wang61, D.R. Ward49,
H.M. Wark54, N.K. Watson47, D. Websdale55, A. Weiden42, M. Whitehead40, J.
Wicht50, G. Wilkinson57,40, M. Wilkinson61, M. Williams40, M.P. Williams47, M.
Williams58, T. Williams47, F.F. Wilson51, J. Wimberley60, M.A. Winn7, J.
Wishahi10, W. Wislicki29, M. Witek27, G. Wormser7, S.A. Wotton49, K.
Wraight53, K. Wyllie40, Y. Xie65, Z. Xu4, Z. Yang3, Z. Yang60, Y. Yao61, H.
Yin65, J. Yu65, X. Yuan61, O. Yushchenko37, K.A. Zarebski47, M.
Zavertyaev11,c, L. Zhang3, Y. Zhang7, A. Zhelezov12, Y. Zheng63, X. Zhu3, V.
Zhukov33, J.B. Zonneveld52, 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
17Universita e INFN, Ferrara, Ferrara, Italy
18Sezione INFN di Firenze, Firenze, Italy
19Laboratori Nazionali dell’INFN di Frascati, Frascati, Italy
20Sezione INFN di Genova, Genova, Italy
21Universita & INFN, Milano-Bicocca, Milano, Italy
22Sezione 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
35Yandex School of Data Analysis, Moscow, Russia
36Budker Institute of Nuclear Physics (SB RAS), Novosibirsk, Russia
37Institute for High Energy Physics (IHEP), Protvino, Russia
38ICCUB, Universitat de Barcelona, Barcelona, Spain
39Universidad de Santiago de Compostela, Santiago de Compostela, Spain
40European Organization for Nuclear Research (CERN), Geneva, Switzerland
41Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL),
Lausanne, Switzerland
42Physik-Institut, Universität Zürich, Zürich, Switzerland
43Nikhef National Institute for Subatomic Physics, Amsterdam, The Netherlands
44Nikhef National Institute for Subatomic Physics and VU University Amsterdam,
Amsterdam, The Netherlands
45NSC Kharkiv Institute of Physics and Technology (NSC KIPT), Kharkiv, Ukraine
46Institute for Nuclear Research of the National Academy of Sciences (KINR),
Kyiv, Ukraine
47University of Birmingham, Birmingham, United Kingdom
48H.H. Wills Physics Laboratory, University of Bristol, Bristol, United
Kingdom
49Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom
50Department of Physics, University of Warwick, Coventry, United Kingdom
51STFC Rutherford Appleton Laboratory, Didcot, United Kingdom
52School of Physics and Astronomy, University of Edinburgh, Edinburgh, United
Kingdom
53School of Physics and Astronomy, University of Glasgow, Glasgow, United
Kingdom
54Oliver Lodge Laboratory, University of Liverpool, Liverpool, United Kingdom
55Imperial College London, London, United Kingdom
56School of Physics and Astronomy, University of Manchester, Manchester,
United Kingdom
57Department of Physics, University of Oxford, Oxford, United Kingdom
58Massachusetts Institute of Technology, Cambridge, MA, United States
59University of Cincinnati, Cincinnati, OH, United States
60University of Maryland, College Park, MD, United States
61Syracuse University, Syracuse, NY, United States
62Pontifícia Universidade Católica do Rio de Janeiro (PUC-Rio), Rio de
Janeiro, Brazil, associated to 2
63University of Chinese Academy of Sciences, Beijing, China, associated to 3
64School of Physics and Technology, Wuhan University, Wuhan, China, associated
to 3
65Institute of Particle Physics, Central China Normal University, Wuhan,
Hubei, China, associated to 3
66Departamento de Fisica , Universidad Nacional de Colombia, Bogota, Colombia,
associated to 8
67Institut für Physik, Universität Rostock, Rostock, Germany, associated to 12
68National Research Centre Kurchatov Institute, Moscow, Russia, associated to
32
69National Research Tomsk Polytechnic University, Tomsk, Russia, associated to
32
70Instituto de Fisica Corpuscular, Centro Mixto Universidad de Valencia -
CSIC, Valencia, Spain, associated to 38
71Van Swinderen Institute, University of Groningen, Groningen, The
Netherlands, associated to 43

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 Genova, Genova, Italy
iUniversità di Milano Bicocca, Milano, Italy
jUniversità di Roma Tor Vergata, Roma, Italy
kUniversità di Roma La Sapienza, Roma, Italy
lAGH - University of Science and Technology, Faculty of Computer Science,
Electronics and Telecommunications, Kraków, Poland
mLIFAELS, La Salle, Universitat Ramon Llull, Barcelona, Spain
nHanoi University of Science, Hanoi, Viet Nam
oUniversità di Padova, Padova, Italy
pUniversità di Pisa, Pisa, Italy
qUniversità degli Studi di Milano, Milano, Italy
rUniversità di Urbino, Urbino, Italy
sUniversità della Basilicata, Potenza, Italy
tScuola Normale Superiore, Pisa, Italy
uUniversità di Modena e Reggio Emilia, Modena, Italy
vIligan Institute of Technology (IIT), Iligan, Philippines
wNovosibirsk State University, Novosibirsk, Russia
†Deceased


# Sparse time-frequency representation via atomic norm minimization

###### Abstract

Nonstationary signals are commonly analyzed and processed in the time-
frequency (T-F) domain that is obtained by the discrete Gabor transform (DGT).
The T-F representation obtained by DGT is spread due to windowing, which may
degrade the performance of T-F domain analysis and processing. To obtain a
well-localized T-F representation, sparsity-aware methods using
$\ell_{1}$-norm have been studied. However, they need to discretize a
continuous parameter onto a grid, which causes a model mismatch. In this
paper, we propose a method of estimating a sparse T-F representation using
atomic norm. The atomic norm enables sparse optimization without
discretization of continuous parameters. Numerical experiments show that the
T-F representation obtained by the proposed method is sparser than the
conventional methods.

Index Terms—  Short-time Fourier transform (STFT), convex optimization, atomic
norm, basis pursuit, semidefinite programming.

##  1 Introduction

Nonstationary signals are commonly analyzed and processed in the time-
frequency (T-F) domain. For converting the time-domain signal into the T-F
domain, the short-time Fourier transform (STFT) or the discrete Gabor
transform (DGT) is usually utilized owing to its simplicity and well-
understood structure [1, 2]. However, the T-F representation obtained by it is
spread due to windowing of the analyzed signal. This spread may affect the
performance of T-F domain analysis and processing.

To achieve a well-localized T-F representation, many approaches have been
proposed [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Reassignment and
synchrosqueezing methods aim to relocate the spread components into the
original positions using phase derivative information [3, 4, 5, 6]. Their
performances are affected by mixing of the components due to windowing [15].
Sparsity-aware methods are powerful tools that are robust against such mixing
of the components and noises [7, 8, 9, 10, 11, 12, 13, 14]. Sparsity-aware
methods aim to find a sparse solution of an underdetermined system. However,
the typical formulation based on $\ell_{1}$-norm minimization involves
discretization of a continuous parameter onto a grid. It may degrade the
performance due to a model mismatch between the signal and the predefined grid
[16].

Recently, sparse optimization using atomic norm has been studied [17, 18, 19,
20] and applied to many applications such as line spectrum estimation [21,
22], direction of arrival estimation [23, 24, 25], and target localization in
radar [26]. Atomic norm does not require discretization of continuous
parameters. Thus, introducing the atomic norm into sparse T-F representation
should obtain a better-localized T-F representation.

In this paper, we propose an estimation method of a sparse T-F representation.
In the proposed method, the estimation problem is formulated as atomic norm
minimization under the condition that the analyzed time-domain signal can be
reconstructed. Numerical experiments confirmed that the proposed method
provides a sparser T-F representation than the conventional methods.

Throughout this paper, ${\mathbb{N}}$, ${\mathbb{R}}$ and ${\mathbb{C}}$
denote the sets of all natural, real and complex numbers, respectively.
${\mathbf{x}}[n]$ is the $n$th element of a vector ${\mathbf{x}}$, and the
${\mathbf{X}}[m,n]$ is the $(m,n)$th element of a matrix ${\mathbf{X}}$.
$\overline{(\cdot)}$, ${(\cdot)}^{\mathrm{T}}$ and ${(\cdot)}^{*}$ represent
the conjugate, the transpose and the conjugate transpose, respectively.
$\operatorname{Tr}({\mathbf{X}})$ is the trace of ${\mathbf{X}}$. The inner
product of vectors ${\mathbf{x}},{\mathbf{y}}$ is defined as
$\left\langle{\mathbf{x}},{\mathbf{y}}\right\rangle={\mathbf{x}}^{\mathrm{T}}\overline{\mathbf{y}}$.
The Frobenius inner product of two matrices ${\mathbf{X}},{\mathbf{Y}}$ is
defined as
$\left\langle{\mathbf{X}},{\mathbf{Y}}\right\rangle_{\mathrm{F}}=\operatorname{Tr}({\mathbf{X}}^{\mathrm{T}}\overline{\mathbf{Y}})$.
$\left\|\cdot\right\|_{\text{F}}$ denotes the Frobenius norm.
${\mathbf{I}}_{L}\in\mathbb{R}^{L\times L}$ is the identity matrix.
${\mathbf{X}}\succeq 0$ stands for ${\mathbf{X}}$ being positive semidefinite.

##  2 Preliminaries

###  2.1 Gabor system and discrete Gabor transform (DGT)

Let ${\mathbf{g}}\in\mathbb{R}^{L}$ denote a window. A Gabor system is defined
as [1, 2]

$\mathcal{G}({\mathbf{g}},a,M)=\left\\{{\mathbf{g}}_{m,n}\right\\}_{m=0,\ldots,M-1,\,n=0,\ldots,N-1},$
(1)

where

${\mathbf{g}}_{m,n}[l]=\mathrm{e}^{\mathrm{i}\frac{2\pi
m(l-an)}{M}}{\mathbf{g}}[l-an],$ (2)

$a\in\mathbb{N}$ is the time-shifting width, and $M\in\mathbb{N}$ is the
number of frequency channels. DGT and the inverse DGT with respect to the
Gabor system $\mathcal{G}({\mathbf{g}},a,M)$ are defined by

$({\mathbf{G}}_{\mathbf{g}}^{*}{\mathbf{f}})[m+nM]=\left\langle{\mathbf{f}},{\mathbf{g}}_{m,n}\right\rangle,\qquad{\mathbf{G}}_{\mathbf{g}}{\mathbf{c}}=\sum_{m,n}{\mathbf{c}}[m+nM]{\mathbf{g}}_{m,n},$

where ${\mathbf{c}}\in{\mathbb{C}}^{MN}$ is a collection of the coefficients
corresponding to a T-F representation. If $\mathcal{G}({\mathbf{g}},a,M)$ is a
frame [1, 2], there exist a dual frame $\mathcal{G}({\mathbf{h}},a,M)$ which
satisfies

${\mathbf{f}}=\sum_{m,n}\left\langle{\mathbf{f}},{\mathbf{g}}_{m,n}\right\rangle{\mathbf{h}}_{m,n}.$
(3)

That is, a T-F representation ${\mathbf{c}}$ satisfying
${\mathbf{f}}={\mathbf{G}}_{\mathbf{g}}{\mathbf{c}}$ can be obtained by DGT
with a dual window ${\mathbf{h}}$ associated with ${\mathbf{g}}$. A standard
constuction of the dual window is the canonical dual window:

$\tilde{\mathbf{h}}=({\mathbf{G}}_{\mathbf{g}}{\mathbf{G}}_{\mathbf{g}}^{*})^{-1}{\mathbf{g}}.$
(4)

###  2.2 Sparse T-F representation using $\ell_{1}$-norm

A T-F representation obtained by DGT ${\mathbf{G}}_{\mathbf{h}}^{*}$ is spread
due to multiplication of a dual window ${\mathbf{h}}$. If
${\mathcal{G}}({\mathbf{g}},a,M)$ is a frame, the T-F representation
${\mathbf{c}}$ is a redundant representation of a time-domain signal
${\mathbf{f}}$, i.e., the T-F representation ${\mathbf{c}}$ satisfying

${\mathbf{f}}={\mathbf{G}}_{\mathbf{g}}{\mathbf{c}}$ (5)

is not unique. The direct formulation for finding the sparsest solution of
this underdetermined system is to minimize the number of non-zero
coefficients, called the $\ell_{0}$-norm. Unfortunately, this problem is
usually an intractable combinatorial optimization problem, and its solution is
sensitive to noise. Instead of the $\ell_{0}$-norm, the $\ell_{1}$-norm has
been widely used as a cost function for promoting sparsity [7, 8, 9, 10, 11,
12, 13, 14]. It is formulated as

$\operatorname*{minimize}_{\mathbf{c}}\quad\left\|{\mathbf{c}}\right\|_{1}\quad\operatorname*{subject~{}to}\quad{\mathbf{f}}={\mathbf{G}}_{\mathbf{g}}{\mathbf{c}},$
(6)

where $\left\|\cdot\right\|_{1}$ is the $\ell_{1}$-norm, defined as
$\left\|{\mathbf{c}}\right\|_{1}=\sum_{j=0}^{NM-1}\left|{\mathbf{c}}[j]\right|$.
This problem is convex, which can be solved by convex optimization algorithms.
Nonetheless, ${\mathcal{G}}({\mathbf{g}},a,M)$ contains windowed sinusoids
whose frequency is discretized onto the grid
$\left\\{m/M\right\\}_{m=0,\ldots,M-1}$. Eq. (6) may provide a poor result
when the signal ${\mathbf{f}}$ has a component whose frequency is not included
in the grid.

###  2.3 Line spectrum estimation using atomic norm

To avoid the effects of grid mismatch, a method using atomic norm has been
studied as a gridless sparse optimization method [17, 18, 19, 20, 21, 22, 23,
24, 25, 26]. Here, line spectrum estimation is applied only to the $n$th
windowed signal. A windowed signal at time index $n$ is denoted as
${\mathbf{f}}_{n}={\mathbf{W}}_{n}{\mathbf{f}}$, where
${\mathbf{W}}_{n}\in\mathbb{R}^{L\times L}$ is a diagonal matrix whose
diagonal elements are given by ${\mathbf{W}}_{n}[l,l]={\mathbf{g}}[l-an]$. We
assume that the $n$th windowed signal ${\mathbf{f}}_{n}$ can be expressed as a
sum of complex sinusoids,

${\mathbf{f}}_{n}={\mathbf{W}}_{n}\sum_{k}c_{n,k}{\mathbf{a}}_{n,k},\qquad{\mathbf{a}}_{n,k}\in{\mathcal{A}},$
(7)

where ${\mathcal{A}}$ is a collection of complex sinusoids

$\mathcal{A}=\left\\{{\mathbf{a}}\in\mathbb{C}^{L}\,\middle|\,{\mathbf{a}}[l]=\mathrm{e}^{\mathrm{i}2\pi\omega
l},\omega\in[0,1)\right\\}.$ (8)

The atomic norm is used to express the $n$th windowed signal
${\mathbf{f}}_{n}$ with a few coefficient $c_{n,k}$. Let us denote
${\mathbf{x}}_{n}=\sum_{k}c_{n,k}{\mathbf{a}}_{n,k}$, then the atomic norm of
${\mathbf{x}}_{n}$ associated with a set of atoms ${\mathcal{A}}$ is given by
[18]

$\displaystyle\left\|{\mathbf{x}}_{n}\right\|_{\mathcal{A}}$
$\displaystyle=\inf\left\\{\nu_{n}\geq
0\,\middle|\,{\mathbf{x}}_{n}\in\nu_{n}\operatorname{conv}({\mathcal{A}})\right\\},$
(9)
$\displaystyle=\inf\left\\{\sum_{k}|c_{n,k}|\,\middle|\,{\mathbf{x}}_{n}=\sum_{k}c_{n,k}{\mathbf{a}}_{n,k},\,\,\,{\mathbf{a}}_{n,k}\in\mathcal{A}\right\\},$

where $\operatorname{conv}({\mathcal{A}})$ is the convex hull of
${\mathcal{A}}$. It corresponds to the infimum of the $\ell_{1}$-norm of
coefficients when ${\mathbf{x}}_{n}$ is represented by a linear combination of
elements in ${\mathcal{A}}$. That is, the atomic norm can be interpreted as an
extension of the $\ell_{1}$-norm to the continuous parameter $\omega\in[0,1)$.
The atomic norm in Eq. (9) is characterized by the following optimization
problem [19]:

$\displaystyle\left\|{\mathbf{x}}_{n}\right\|_{\mathcal{A}}=$
$\displaystyle\min_{{\mathbf{u}}_{n},\nu_{n}}\quad\frac{1}{2L}\operatorname{Tr}\left(T({\mathbf{u}}_{n})\right)+\frac{1}{2}\nu_{n}$
$\displaystyle\operatorname*{subject~{}to}\quad\left[\begin{array}[]{cc}T({\mathbf{u}}_{n})&{\mathbf{x}}_{n}\\\
{\mathbf{x}}_{n}^{*}&\nu_{n}\end{array}\right]\succeq 0,$ (12)

where $T:{\mathbb{C}}^{L}\to{\mathbb{C}}^{L\times L}$ is the Hermitian
Toeplitz operator:

$T({\mathbf{u}})=\left[\begin{array}[]{cccc}{\mathbf{u}}[0]&\overline{{\mathbf{u}}[1]}&\cdots&\overline{{\mathbf{u}}[L-1]}\\\
{\mathbf{u}}[1]&{\mathbf{u}}[0]&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\overline{{\mathbf{u}}[1]}\\\
{\mathbf{u}}[L-1]&\cdots&{\mathbf{u}}[1]&{\mathbf{u}}[0]\end{array}\right].$
(13)

If the Hermitian Toeplitz matrix $T({\mathbf{u}}_{n})$ is positive
semidefinite and singular, it can be uniquely decomposed as [27]

$T({\mathbf{u}}_{n})=\sum_{k=0}^{K-1}\left|c_{n,k}\right|{\mathbf{a}}_{n,k}{\mathbf{a}}^{*}_{n,k},\qquad{\mathbf{a}}_{n,k}\in\mathcal{A},$
(14)

where $K$ corresponds to the rank of $T({\mathbf{u}}_{n})$.
${\mathbf{a}}_{n,k}$ in Eq. (14) can be obtained by Prony’s method [28], the
matrix pencil method [29], or other linear prediction methods [30]. Then,
coefficients $c_{n,k}$ can be obtained by solving the linear equation

$\left[{\mathbf{a}}_{n,0},\ldots,{\mathbf{a}}_{n,K-1}\right]{\mathbf{c}}_{n}={\mathbf{x}_{n}},$
(15)

where ${\mathbf{c}}_{n}=\left[c_{n,0},\ldots,c_{n,K-1}\right]^{\mathrm{T}}$.
Therefore, estimating sinusoids in the windowed signal ${\mathbf{f}}_{n}$
using the atomic norm is formulated as

$\operatorname*{minimize}_{\mathbf{x}_{n}}\quad\left\|{\mathbf{x}}_{n}\right\|_{\mathcal{A}}\quad\operatorname*{subject~{}to}\quad{\mathbf{f}}_{n}={\mathbf{W}}_{n}{\mathbf{x}}_{n}.$
(16)

Substituting Eq. (12) into Eq. (16), Eq. (16) can be rewritten as the
following semidefinite programming:

$\displaystyle\operatorname*{minimize}_{\mathbf{x}_{n},\mathbf{u}_{n},\nu_{n}}$
$\displaystyle\quad\frac{1}{2L}\operatorname{Tr}(T({\mathbf{u}}_{n}))+\frac{1}{2}\nu_{n}$
$\displaystyle\operatorname*{subject~{}to}$
$\displaystyle\quad\left[\begin{array}[]{cc}T({\mathbf{u}}_{n})&{\mathbf{x}}_{n}\\\
{\mathbf{x}}_{n}^{*}&\nu_{n}\end{array}\right]\succeq
0,\quad{\mathbf{f}}_{n}={\mathbf{W}}_{n}{\mathbf{x}}_{n}.$ (19)

While the eigenvalue decomposition of the $(L+1)\times(L+1)$ matrix needs to
be iterated when solving Eq. (19), it can be reduced to the eigenvalue
decomposition of a $(J+1)\times(J+1)$ matrix if the window function is
supported on $[0,J-1]$.

##  3 Proposed method

The method of estimating the T-F representation based on the $\ell_{1}$-norm
suffers from degradation of performance due to discretization onto the grid as
shown in Fig. 1 (a). On the other hand, performing the line spectrum
estimation using the atomic norm for each time index, a sparse T-F
representation can be obtained without a grid in frequency direction. However,
Eq. (19) estimates sinusoids at each time index $n$ independently as Fig. 1
(b), which does not efficiently take advantage of the sparsity of T-F
representation (see Fig. 2).

In this paper, we propose a estimation method of a sparse T-F representation
using the atomic norm to avoid the effect of grid mismatch. In our
formulation, the atomic norm is minimized under the constraint of the perfect
reconstruction of the entire signal. It can be interpreted as an extension of
Eq. (6) with an infinite number of frequency channels as Fig. 1 (c).

Fig. 1: Estimating T-F representations using (a) $\ell_{1}$-norm, (b) the line
spectrum estimation by atomic norm, and (c) the proposed method.

###  3.1 Sparse T-F representation using atomic norm

Reconstruction of the $n$th windowed signal
${\mathbf{f}}_{n}={\mathbf{W}}_{n}{\mathbf{f}}_{n}$ is considered in Eq. (19).
In contrast, reconstructing the entire signal ${\mathbf{f}}$ by windowing and
summing ${\mathbf{x}}_{n}$ is considered in the proposed method, which is
written as

${\mathbf{A}}_{\mathbf{g}}{\mathbf{x}}=\sum_{n=0}^{N-1}{\mathbf{W}}_{n}{\mathbf{x}}_{n},$
(20)

where
${\mathbf{x}}=\left[{\mathbf{x}}_{0}^{\mathrm{T}},{\mathbf{x}}_{1}^{\mathrm{T}},\ldots,{\mathbf{x}}_{N-1}^{\mathrm{T}}\right]^{\mathrm{T}}$.
As a cost function, a sum of atomic norms for time index $n$ is chosen because
it corresponds to the grid-less version of the cost function in Eq. (6),

$\sum_{n=0}^{N-1}\left\|{\mathbf{x}}_{n}\right\|_{\mathcal{A}}=\inf\left\\{\sum_{n,k}|c_{n,k}|\,\middle|\,{\mathbf{x}}_{n}=\sum_{k}c_{n,k}{\mathbf{a}}_{n,k},\,\,\,{\mathbf{a}}_{n,k}\in\mathcal{A}\right\\}.$

Thus, the estimation of a sparse T-F representation is formulated as

$\operatorname*{minimize}_{{\mathbf{x}}}\quad\sum_{n=0}^{N-1}\left\|{\mathbf{x}}_{n}\right\|_{\mathcal{A}}\quad\operatorname*{subject~{}to}\quad{\mathbf{f}}={\mathbf{A}}_{\mathbf{g}}{\mathbf{x}}.$
(21)

Since ${\mathbf{x}}_{n}=\sum_{k}c_{n,k}{\mathbf{a}}_{n,k}$, the $n$th windowed
signal
${\mathbf{W}}_{n}{\mathbf{x}}_{n}=\sum_{k}c_{n,k}({\mathbf{W}}_{n}{\mathbf{a}}_{n,k})$
corresponds to the vertical sum of elements in Fig. 1. The signal
${\mathbf{f}}$ can be reconstructed by
$\sum_{n=0}^{N-1}{\mathbf{W}}_{n}{\mathbf{x}}_{n}={\mathbf{A}}_{\mathbf{g}}{\mathbf{x}}$.
After solving Eq. (21), the coefficients $c_{n,k}$ can be calculated as the
line spectrum estimation using atomic norm (in Sec. 2.3).

###  3.2 Algorithm for solving Eq. (21)

We firstly reformulate it as a semidefinite programming to solve Eq. (21).
Substituting Eq. (12) into Eq. (21), it can be rewritten as

$\displaystyle\operatorname*{minimize}_{{\mathbf{x}},{\mathbf{u}},{\boldsymbol{\nu}}}$
$\displaystyle\quad\sum_{n=0}^{N-1}\frac{1}{2L}\mathrm{Tr}(T({\mathbf{u}}_{n}))+\frac{1}{2}\nu_{n}$
$\displaystyle\mathop{\rm subject~{}to}$
$\displaystyle\quad\left[\begin{array}[]{cc}T({\mathbf{u}}_{n})&{\mathbf{x}}_{n}\\\
{\mathbf{x}}_{n}^{*}&\nu_{n}\end{array}\right]\succeq 0,\text{ for
}n=0,\ldots,N-1$ (24)
$\displaystyle\quad{\mathbf{f}}={\mathbf{A}}_{\mathbf{g}}{\mathbf{x}},$ (25)

where
${\mathbf{u}}=\left[{\mathbf{u}}_{0}^{\mathrm{T}},{\mathbf{u}}_{1}^{\mathrm{T}},\ldots,{\mathbf{u}}_{N-1}^{\mathrm{T}}\right]^{\mathrm{T}}$
and
${\boldsymbol{\nu}}=\left[\nu_{0},\nu_{1},\ldots,\nu_{N-1}\right]^{\mathrm{T}}$.
We adopt the alternating direction method of multipliers (ADMM) [31] to solve
Eq. (25). For applying ADMM to the proposed method, we introduce auxiliary
variables ${\mathbf{Z}}_{n}\in\mathbb{C}^{(L+1)\times(L+1)}$ for
$n=0,\ldots,N-1$ and a set corresponding to the reconstruction constraint
$C=\left\\{{\mathbf{x}}\middle|{\mathbf{A}}_{\mathbf{g}}{\mathbf{x}}={\mathbf{f}}\right\\}$.
Then, Eq. (25) is reformulated as

$\displaystyle\operatorname*{minimize}_{\begin{subarray}{c}{\mathbf{x}}\in
C,{\mathbf{u}},{\boldsymbol{\nu}}\\\ {\mathbf{Z}}_{n}\succeq 0\end{subarray}}$
$\displaystyle\quad\sum_{n=0}^{N-1}\frac{1}{2L}\mathrm{Tr}(T({\mathbf{u}}_{n}))+\frac{1}{2}\nu_{n}$
$\displaystyle\mathop{\rm subject~{}to}$
$\displaystyle\quad\left[\begin{array}[]{cc}T({\mathbf{u}}_{n})&{\mathbf{x}}_{n}\\\
{\mathbf{x}}_{n}^{*}&\nu_{n}\end{array}\right]={\mathbf{Z}}_{n},\text{ for
}n=0,\ldots,N-1.$ (28)

The augmented Lagrangian associated with this problem is given by

$\displaystyle{\mathcal{L}}({\mathbf{x}},{\mathbf{u}},{\boldsymbol{\nu}},{\mathbf{Z}}_{n},{\mathbf{\Lambda}}_{n})=$
$\displaystyle\sum_{n=0}^{N-1}\frac{1}{2L}\mathrm{Tr}(T({\mathbf{u}}_{n}))+\frac{1}{2}\nu_{n}$
$\displaystyle+\left\langle{\mathbf{\Lambda}}_{n},\left[\begin{array}[]{cc}T({\mathbf{u}}_{n})&{\mathbf{x}}_{n}\\\
{\mathbf{x}}_{n}^{*}&\nu_{n}\end{array}\right]-{\mathbf{Z}}_{n}\right\rangle_{\mathrm{F}}$
(31)
$\displaystyle+\frac{\rho}{2}\left\|\left[\begin{array}[]{cc}T({\mathbf{u}}_{n})&{\mathbf{x}}_{n}\\\
{\mathbf{x}}_{n}^{*}&\nu_{n}\end{array}\right]-{\mathbf{Z}}_{n}\right\|_{\mathrm{F}}^{2},$
(34)

where ${\mathbf{\Lambda}}_{n}\in\mathbb{C}^{(L+1)\times(L+1)}$ for
$n=0,\ldots,N-1$ are dual variables, and $\rho>0$ is the augmented Lagrangian
parameter. Then, ADMM consists of the following iterations:

$\displaystyle({\mathbf{x}}^{(i+1)},$
$\displaystyle{\mathbf{u}}^{(i+1)},{\boldsymbol{\nu}}^{(i+1)})=\operatorname*{argmin}_{{\mathbf{x}}\in
C,{\mathbf{u}},{\boldsymbol{\nu}}}\,\,{\mathcal{L}}({\mathbf{x}},{\mathbf{u}},{\boldsymbol{\nu}},{\mathbf{Z}}_{n}^{(i)},{\mathbf{\Lambda}}_{n}^{(i)}),$
(35) $\displaystyle{\mathbf{Z}}_{n}^{(i+1)}=$
$\displaystyle\operatorname*{argmin}_{{\mathbf{Z}}_{n}\succeq
0}\,\,{\mathcal{L}}({\mathbf{x}}^{(i+1)},{\mathbf{u}}^{(i+1)},{\boldsymbol{\nu}}^{(i+1)},{\mathbf{Z}}_{n},{\mathbf{\Lambda}}_{n}^{(i)}),$
(36) $\displaystyle{\mathbf{\Lambda}}_{n}^{(i+1)}=$
$\displaystyle{\mathbf{\Lambda}}_{n}^{(i)}+\rho\left(\left[\begin{array}[]{cc}T({\mathbf{u}}_{n}^{(i+1)})&{\mathbf{x}}_{n}^{(i+1)}\\\
({\mathbf{x}}_{n}^{(i+1)})^{*}&\nu_{n}^{(i+1)}\end{array}\right]-{\mathbf{Z}}_{n}^{(i+1)}\right).$
(39)

Introducing

${\mathbf{Z}}_{n}=\left[\begin{array}[]{cc}{{\mathbf{Z}}_{\mathrm{T}}}_{n}&{{\mathbf{z}}_{\mathrm{x}}}_{n}\\\
{{\mathbf{z}}_{\mathrm{x}}}_{n}^{*}&{z_{\nu}}_{n}\end{array}\right],\quad{\mathbf{\Lambda}}_{n}=\left[\begin{array}[]{cc}{{\mathbf{\Lambda}}_{\mathrm{T}}}_{n}&{{\boldsymbol{\lambda}}_{\mathrm{x}}}_{n}\\\
{{\boldsymbol{\lambda}}_{\mathrm{x}}}_{n}^{*}&{\lambda_{\nu}}_{n}\end{array}\right],$

${\mathbf{z}}_{\mathrm{x}}=\left[{{\mathbf{z}}_{\mathrm{x}}}_{0}^{\mathrm{T}},{{\mathbf{z}}_{\mathrm{x}}}_{1}^{\mathrm{T}},\ldots,{{\mathbf{z}}_{\mathrm{x}}}_{N-1}^{\mathrm{T}}\right]^{\mathrm{T}}$,
and
${\boldsymbol{\lambda}}_{\mathrm{x}}=\left[{{\boldsymbol{\lambda}}_{\mathrm{x}}}_{0}^{\mathrm{T}},{{\boldsymbol{\lambda}}_{\mathrm{x}}}_{1}^{\mathrm{T}},\ldots,{{\boldsymbol{\lambda}}_{\mathrm{x}}}_{N-1}^{\mathrm{T}}\right]^{\mathrm{T}}$,
Eq. (35) can be separately solved for ${\mathbf{x}}$, ${\mathbf{u}}$ and
${\boldsymbol{\nu}}$.

The proposed algorithm is summarized in Algorithm 1. The update for
${\mathbf{x}}$ is written as the projection onto the set
$C=\left\\{{\mathbf{x}}\middle|{\mathbf{A}}_{\mathbf{g}}{\mathbf{x}}={\mathbf{f}}\right\\}$,
denoted by $P_{C}(\cdot)$. It is given by

$P_{C}({\mathbf{v}})={\mathbf{v}}-{\mathbf{A}}_{\mathbf{g}}^{*}({\mathbf{A}}_{\mathbf{g}}{\mathbf{A}}_{\mathbf{g}}^{*})^{-1}({\mathbf{A}}_{\mathbf{g}}{\mathbf{v}}-{\mathbf{f}}).$
(40)

The computation of matrix inversion in Eq. (40) can be avoided by the
cannonical dual window $\tilde{\mathbf{h}}$ as

${\mathbf{A}}_{\mathbf{g}}^{*}\left({\mathbf{A}}_{\mathbf{g}}{\mathbf{A}}_{\mathbf{g}}^{*}\right)^{-1}={\mathbf{A}}_{\tilde{\mathbf{h}}}^{*}.$
(41)

The updates for ${\mathbf{u}}$, ${\boldsymbol{\nu}}$, ${\mathbf{Z}}_{n}$, and
${\mathbf{\Lambda}}_{n}$ can be computed in parallel for each $n$.
$T^{\dagger}:{\mathbb{C}}^{L\times L}\to{\mathbb{C}}^{L}$ in the update for
${\mathbf{u}}$ represents the pseudo-inverse operator of $T$,

$T^{\dagger}({\mathbf{X}})[n]=\frac{1}{2(L-n)}\sum_{k=0}^{L-n-1}\left({\mathbf{X}}[k,k+n]+\overline{\mathbf{X}[k+n,k]}\right).$

The update for ${\mathbf{Z}}_{n}$ can be calculated by a projection onto the
positive semidefinite cone ${\mathbb{S}}_{+}$, which is implemented by setting
the negative eigenvalues to 0. Since Eq. (25) is a convex optimization
problem, Algorithm 1 can obtain the global optimal solution regardless of
initialization for ${\mathbf{Z}}_{n}$ and ${\mathbf{\Lambda}}_{n}$.

Algorithm 1 ADMM for solving Eq. (21)

${\mathbf{A}}$, ${\mathbf{f}}$, $\rho$

${\mathbf{x}}$, ${\mathbf{u}}$, ${\boldsymbol{\nu}}$

Initialize ${\mathbf{Z}}_{n}$ and ${\mathbf{\Lambda}}_{n}$ for
$n=0,\ldots,N-1$

for $i=0,1,\cdots$ do

${\mathbf{x}}\leftarrow
P_{C}\left({\mathbf{z}}_{\mathrm{x}}-\frac{1}{\rho}{\boldsymbol{\lambda}}_{\mathrm{x}}\right)$

for $n=1,\cdots,N$ do

${\mathbf{u}}_{n}\leftarrow
T^{\dagger}\left({{\mathbf{Z}}_{\mathrm{T}}}_{n}-\frac{1}{\rho}\left({{\mathbf{\Lambda}}_{\mathrm{T}}}_{n}+\frac{1}{2}{\mathbf{I}}_{L}\right)\right)$

$\nu_{n}\leftarrow{z_{\nu}}_{n}-\frac{1}{\rho}\left({\lambda_{\nu}}_{n}+\frac{1}{2}\right)$

${\mathbf{Z}}_{n}\leftarrow
P_{{\mathbb{S}}_{+}}\left(\left[\begin{array}[]{cc}T({\mathbf{u}}_{n})&{\mathbf{x}}_{n}\\\
{\mathbf{x}}_{n}^{*}&\nu_{n}\end{array}\right]+\frac{1}{\rho}{\mathbf{\Lambda}}_{n}\right)$

${\mathbf{\Lambda}}_{n}\leftarrow{\mathbf{\Lambda}}_{n}+\rho\left(\left[\begin{array}[]{cc}T({\mathbf{u}}_{n})&{\mathbf{x}}_{n}\\\
{\mathbf{x}}_{n}^{*}&\nu_{n}\end{array}\right]-{\mathbf{Z}}_{n}\right)$

end for

end for

##  4 Numerical experiments

To evaluate the performance of the proposed method, firstly, the proposed
method was applied to an artificial signal that contains a sinusoid, a linear
chirp, and a quadratic chirp. The proposed method was compared to DGT with the
canonical dual window, the reassignment method [3], the $\ell_{1}$-norm
minimization [11], and the window-wise atomic norm minimization. The Slepian
window [32] was chosen as a window ${\mathbf{g}}$, whose length and bandwidth
were set to $2^{7}$ and $0.04$. The time-shifting width and the number of
frequency channels were set to $a=2^{4}$ and $M=2^{10}$. ${\mathbf{Z}}_{n}$
and ${\mathbf{\Lambda}}_{n}$ in Algorithm 1 were initialized to the zero
matrices. Prony’s method was used to estimate ${\mathbf{a}}_{n,k}$ from the
solution of Eq. (25).

The estimated T-F representations are shown in Fig. 2. The T-F representation
obtained using the $\ell_{1}$-norm was better-localized than DGT and the
reassignment, but it had multiple non-zero coefficients in each time index to
represent a sinusoid. The window-wise atomic norm minimization estimated a
sparse representation corresponding to the sinusoid. On the other hand, it
cannot express the chirps sparsely. The T-F representation obtained by the
proposed method was the most-localized among these T-F representations.

Then, to evaluate the sparseness of these T-F representations, the squared
absolute value of coefficients are plotted in Fig. 3. It can be seen that the
proposed method can represent the signal using the least coefficients. These
results indicate that the proposed method provides a sparse representation
using the atomic norm while taking into account the relationship among each
time index.

Finally, the proposed method was applied to a speech signal. The settings
associated with the Gabor system and the proposed algorithm were the same as
in the previous experiment. The estimated T-F representations are shown in
Fig. 4. It can be seen that the T-F representation obtained by the proposed
method was the sparsest among these T-F representations. The result suggests
that the proposed method performs well also for a real audio signal.

Fig. 2: T-F representations of an artificial signal. Each column shows (from
left to right) the T-F representations obtained by DGT, the reassignment
method, the $\ell_{1}$-norm minimization, the window-wise atomic norm
minimization, and the proposed method, respectively. The bottom row
illustrates these enlargements in the red box.

##  5 Conclusion

In this paper, we proposed the method of estimating sparse T-F representation
via the atomic norm minimization. The proposed method estimates a sparse T-F
representation without discretization of frequency using the atomic norm. The
experimental results show that the proposed method can estimate a sparser T-F
representation than the existing methods. Future work includes applications of
the proposed method to denoising and mode decomposition.

Fig. 3: Squared absolute value of coefficients in decending order. Fig. 4: T-F
representations of a speech signal. Each row and column represents the same as
Fig. 2.

## References

  * [1] H. G. Feichtinger and T. Strohmer, Eds.,  Gabor Analysis and Algorithms,  Birkhäuser Boston, Boston, MA, 1998.
  * [2] K. Gröchenig,  Foundations of Time-Frequency Analysis,  Birkhäuser Boston, Boston, MA, 2001.
  * [3] F. Auger and P. Flandrin,  “Improving the readability of time-frequency and time-scale representations by the reassignment method,”  IEEE Trans. Signal Process., vol. 43, no. 5, pp. 1068–1089, May 1995.
  * [4] I. Daubechies, J. Lu, and H. T. Wu,  “Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool,”  Appl. Comput. Harmon. Anal., vol. 30, no. 2, pp. 243–261, 2011\.
  * [5] S. Meignen, T. Oberlin, and S. McLaughlin,  “A new algorithm for multicomponent signals analysis based on synchrosqueezing: With an application to signal sampling and denoising,”  IEEE Trans. Signal Process., vol. 60, no. 11, pp. 5787–5798, Nov. 2012.
  * [6] F. Auger, P. Flandrin, Y.-T. Lin, S. McLaughlin, S. Meignen, T. Oberlin, and H.-T. Wu,  “Time-frequency reassignment and synchrosqueezing: An overview,”  IEEE Signal Process. Mag., vol. 30, no. 6, pp. 32–41, Nov. 2013\.
  * [7] G. E. Pfander and H. Rauhut,  “Sparsity in time-frequency representations,”  J. Fourier Anal. Appl., vol. 16, pp. 233–260, 2010.
  * [8] A. Gholami,  “Sparse time-frequency decomposition and some applications,”  IEEE Trans. Geosci. Remote Sens., vol. 51, no. 6, pp. 3598–3604, June 2013.
  * [9] M. D. Plumbley, T. Blumensath, L. Daudet, R. Gribonval, and M. E. Davies,  “Sparse representations in audio and music: From coding to source separation,”  Proc. IEEE, vol. 98, no. 6, pp. 995–1005, 2010.
  * [10] I. Bayram,  “Mixed norms with overlapping groups as signal priors,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), May 2011, pp. 4036–4039.
  * [11] P. Balazs, M. Doerfler, M. Kowalski, and B. Torresani,  “Adapted and adaptive linear time-frequency representations: A synthesis point of view,”  IEEE Signal Process. Mag., vol. 30, no. 6, pp. 20–31, 2013.
  * [12] P.-Y. Chen and I. W. Selesnick,  “Group-sparse signal denoising: Non-convex regularization, convex optimization,”  IEEE Trans. Signal Process., vol. 62, no. 13, pp. 3464–3478, July 2014.
  * [13] E. Sejdić, I. Orović, and S. Stanković,  “Compressive sensing meets time-frequency: An overview of recent advances in time-frequency processing of sparse signals,”  Digital Signal Process., vol. 77, pp. 22–35, 2018.
  * [14] M. Kowalski, A. Meynard, and H. T. Wu,  “Convex optimization approach to signals with fast varying instantaneous frequency,”  Appl. Comput. Harmon. Anal., vol. 44, no. 1, pp. 89–122, 2018.
  * [15] T. Kusano, K. Yatabe, and Y. Oikawa,  “Maximally energy-concentrated differential window for phase-aware signal processing using instantaneous frequency,”  in Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP), May 2020, pp. 5825–5829.
  * [16] Y. Chi, L. L. Scharf, A. Pezeshki, and A. R. Calderbank,  “Sensitivity to basis mismatch in compressed sensing,”  IEEE Trans. Signal Process., vol. 59, no. 5, pp. 2182–2195, 2011\.
  * [17] E. J. Candès and C. Fernandez-Granda,  “Towards a mathematical theory of super-resolution,”  Commun. Pure Appl. Math., vol. 67, no. 6, pp. 906–956, 2014.
  * [18] V. Chandrasekaran, B. Recht, P. A. Parrilo, and A. S. Willsky,  “The convex geometry of linear inverse problems,”  Found. Comput. Math., vol. 12, no. 6, pp. 805–849, Dec. 2012.
  * [19] G. Tang, B. N. Bhaskar, P. Shah, and B. Recht,  “Compressed sensing off the grid,”  IEEE Trans. Inf. Theory, vol. 59, no. 11, pp. 7465–7490, 2013.
  * [20] Y. Chi and M. Ferreira Da Costa,  “Harnessing sparsity over the continuum: Atomic norm minimization for superresolution,”  IEEE Signal Process. Mag., vol. 37, no. 2, pp. 39–57, 2020.
  * [21] Z. Yang and L. Xie,  “On gridless sparse methods for line spectral estimation from complete and incomplete data,”  IEEE Trans. Signal Process., vol. 63, no. 12, pp. 3139–3153, 2015\.
  * [22] Y. Li and Y. Chi,  “Off-the-grid line spectrum denoising and estimation with multiple measurement vectors,”  IEEE Trans. Signal Process., vol. 64, no. 5, pp. 1257–1269, Mar. 2016.
  * [23] A. Xenaki and P. Gerstoft,  “Grid-free compressive beamforming,”  J. Acoust. Soc. Am., vol. 137, no. 4, pp. 1923–1935, Apr. 2015\.
  * [24] Y. Yang, Z. Chu, Z. Xu, and G. Ping,  “Two-dimensional grid-free compressive beamforming,”  J. Acoust. Soc. Am., vol. 142, no. 2, pp. 618–629, 2017.
  * [25] K. Mahata and M. M. Hyder,  “Grid-less T.V minimization for DOA estimation,”  Signal Process., vol. 132, pp. 155–164, Mar. 2017.
  * [26] R. Heckel, V. I. Morgenshtern, and M. Soltanolkotabi,  “Super-resolution radar,”  Inf. Inference, vol. 5, no. 1, pp. 22–75, Mar. 2016.
  * [27] C. Caratheodory and L. Fejér,  “über den zusammenhang der extremen von harmonischen funktionen mit ihren koeffizienten und über den picard-landau’schen satz,”  Rendiconti Circolo Mat. Palermo, vol. 32, no. 1, pp. 218–239, Dec. 1911.
  * [28] G. Riche de Prony,  “Essai expérimental et analytique: Sur les lois de la dilatabilité de fluides élastiques et sur celles de la force expansive de la vapeur de l’eau et de la vapeur de l’alcool à différentes températures,”  J. L’École Polytech., vol. 1, no. 22, pp. 24–76, 1795.
  * [29] Y. Hua and T. Sarkar,  “Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise,”  IEEE Trans. Acoust. Speech Signal Process., vol. 38, no. 5, pp. 814–824, May 1990.
  * [30] D. Tufts and R. Kumaresan,  “Estimation of frequencies of multiple sinusoids: Making linear prediction perform like maximum likelihood,”  Proc. IEEE, vol. 70, no. 9, pp. 975–989, 1982.
  * [31] S. Boyd, N. Parikh, E. Chu, J. Eckstein, S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein,  “Distributed optimization and statistical learning via the alternating direction method of multipliers,”  Found. Trends Mach. Learn., vol. 3, no. 1, pp. 1–122, 2010.
  * [32] D. Slepian,  “Prolate spheroidal wave functions, Fourier analysis, and uncertainty-V: The discrete case,”  Bell Syst. Tech. J., vol. 57, no. 5, pp. 1371–1430, May 1978.