# LASSO Estimation in Exponential Random Graph models

Sergio Buttazzo
Göran Kauermann

###### Abstract

The paper demonstrates the use of LASSO-based estimation in network models.
Taking the Exponential Random Graph Model (ERGM) as a flexible and widely used
model for network data analysis, the paper focuses on the question of how to
specify the (sufficient) statistics, that define the model structure. This
includes both, endogenous network statistics (e.g. two-stars, triangles, etc.)
as well as statistics involving exogenous covariates; on the node as well as
on the edge level. LASSO estimation is a penalized estimation that shrinks
some of the parameter estimates to be equal to zero. As such it allows for
model selection by modifying the amount of penalty. The concept is well
established in standard regression and we demonstrate its usage in network
data analysis, with the advantage of ’automatically’ providing a model
selection framework.

##  1 Introduction

The analysis of network data with exponential random graph models (ERGM) has
become very popular in recent years. Originally proposed by [19] and [11], the
model classed gained interest with stable fitting procedures proposed and
discussed in [41], [24] and [20], and implemented ready to use software, see
[23] and [17]. We refer to [27] and [28] for a general introduction to network
data analysis, see also [33] or [18] for a focus on ERGMs or more generally
[5] and [7]. The model class has been extended in various ways by including
random components ([25] or [44]), temporal aspects ([31] or [32]) or mixtures
of models ([40]), to mention a few ideas. Generally, the core of an ERGM is
the specification of sufficient network statistics that define the model
structure and allow for interpretations regarding the network formation.

A network consists of a set of nodes $\mathcal{N}$ and a set of edges
$E\subset\mathcal{N}\times\mathcal{N}$, which can be expressed as a network
related adjacency matrix $Y\in\\{0,1\\}^{N\times N}$ where
$Y_{ij}=1\Leftrightarrow(i,j)\in E$, i.e. if there is an edge from node $i$ to
node $j$. Assuming undirected edges, meaning that network ties do not contain
information about the orientation of the tie, we have $Y_{ij}=Y_{ji}$, i.e.
the adjacency matrix $Y$ is symmetric. Self-loops are typically left
undefined, meaning that the diagonal of $Y$ is not specified. With
$N=|\mathcal{N}|$ we define the number of nodes in the network. In this paper,
we will restrict ourselves to undirected networks, even though all proposed
methods in this paper can be extended readily to directed ties.

The general structure of an ERGM is of exponential family type, i.e. the
distribution of the random matrix $Y$ is given by:

$\displaystyle
P(Y=y;\theta)=\exp\\{\textbf{s}(y)^{\top}\theta-\kappa(\theta)\\}.$ (1)

Here $\theta\in\Theta$ is the parameter of interest and $\kappa(\theta)$ is a
log-normalization constant. The vector
$\textbf{s}(y)=(s_{0}(y),s_{1}(y),s_{2}(y),\cdots,s_{p}(y))^{\top}$ contains
the model statistics expressing endogenous network effects. The first
component $s_{0}(y)=\sum_{i=1}^{n}\sum_{j>i}y_{ij}$ takes the role of the
“intercept” in ERGMs. The remaining components take the role of model
variables. A number of different reasonable structural statistics have been
proposed in [42] and [35], see also [43]. The model can easily be extended to
also include exogenous covariates $x$, which results by constructing network
statistics $\textbf{s}(y,x)$ that depend on both, the network as well as tie
or node-specific quantities $x$. For ease of notation, we omit this
generalization for the moment.

In general, the concrete choice of the statistics typically depends on the
research questions one wants to focus on. In other words, quite commonly in
network data analysis, one selects the set of statistics based on substance
matter research questions and then estimates the specified model and
interprets the resulting parameter estimates. This approach has a couple of
pitfalls. First, many of the commonly used statistics are heavily correlated,
as they measure similar properties of the network. This leads to excessively
high running times and poor model performances when used in combination with
each other. Secondly, the model statistics need to be specified in advance,
which requires some level of expertise by the researcher. Third, overall,
model evaluation of the fitted model is necessary but it might be expedient to
have model selection and model estimation somewhat combined.

In this paper, we propose a data-driven approach to select the statistics in
$\textbf{s}(y)$. To do so we make use of LASSO, where the acronym stands for
Least Absolute Shrinkage to obtain a model Selection Operator, which allows
for data-driven model selection. To do so we start with a rich selection, that
is we include numerous statistics in the statistics vector $\textbf{s}(y)$ so
that its dimension $p$ is large. This in turn requires inducing a prior
structure on the parameter in terms of penalization to make the model
estimable. Originally proposed in [45], see also [39], the method has gained
wide and still increasing interest in statistics, data analysis, and machine
learning, see exemplary [49]. A Bayesian formulation of the LASSO approach is
provided in [36]. General properties of the LASSO when applied to OLS
regression have been extensively studied. [12] show consistency and asymptotic
normality for the LASSO estimator, provided $\lambda$ grows as $o(n)$, where
$n$ is the size of the sample. [48] provide a necessary and sufficient
condition for the LASSO estimator to achieve consistent selection (i.e.
selecting the true parameters).

Interestingly enough, variable selection and the development of model quality
measures in general, even though having huge importance when dealing with
ERGMs, have received little academic attention compared to other issues
connected with the statistical modeling of network data. An example is given
by [46], who introduce and discuss a bottom-up, stepwise variable selection
algorithm based on relative AIC change. We also refer to [22] who propose to
assess an ERGM’s performance by directly comparing the observed network with a
sample from the model, using a set of ad-hoc chosen statistics as a judgment
meter. LASSO itself has been applied to ERGMs by [15], even though their work
is focused on an entirely Bayesian point of view. Generally, though, while
LASSO estimation is widely used in a classical regression setup, its
application in ERGMs is certainly less developed. Some first ideas have been
proposed for instance in [3] or in recent dissertations, see [37], [34],[2].
In this paper, we elaborate this further and focus on the use of LASSO
estimation for model selection of ERGMs. It is noteworthy that the parameters
given by the LASSO approach are obtained through the introduction of a bias,
that leads to decreasing performances in terms of model fit. Instead, we
pursue to make use of the LASSO to obtain a ranking of importance of the model
variables and use this ranking as a guideline regarding which variables to
include in a new model, fit classically.

The paper is structured as follows: in Section 2 we define the mathematical
foundations upon which the LASSO variable selector is built. First we
explicitly define the ERGMs by characterizing the probability measure that
defines this class of models. Additionally, we present the set of network
statistics that are commonly used as model variables and discuss parameter
estimation. Subsequently, we switch the focus on a penalized version of the
ERGM likelihood function and derive a variable selection algorithm based on
penalized parameter estimation. In Section 3 we discuss applications of the
above-discussed method, by presenting results regarding a number of simulated
networks and two real-life (undirected, binary) networks. Finally, in Section
4, we provide a discussion on potential future research questions linked with
the topic at hand.

##  2 LASSO Estimation in ERGMs

###  2.1 Parameter estimation in ERGMs

We start our presentation with a short review of parameter estimation in the
ERGM (1). This will ease the presentation of LASSO estimation subsequently.

The problematic component in (1 is given by the normalizing constant

$\displaystyle\exp\left(\kappa\left(\theta\right)\right)=\sum_{y\in\mathcal{Y}}\textbf{s}(y)^{\top}\theta$
(2)

where the sum is over all possible network formations with $N$ nodes.
Calculation of (2) is numerically infeasible unless for miniature networks,
but the obstacle can be overcome by Markov Chain Monte Carlo (MCMC) methods
([13]). The idea of the approach is to estimate the parameters iteratively,
relying on samples from a temporary version of the model at each step. Markov
chain-based methods are particularly useful in the case of ERGMs: since the
likelihood function is only needed up to a normalization constant, a sample
can be obtained without having to directly evaluate the constant
$\kappa(\theta)$. See [41], [9] or [24] for further details.

If the normalization constant were tractable, the MLE could be directly
computed by the Newton-Raphson algorithm with the iteration step:

$\theta^{(n+1)}=\theta^{(n)}-\Sigma_{\theta^{(n)}}\left(\textbf{s}(Y)\right)^{-1}\left(\mathbb{E}_{\theta^{(n)}}[\textbf{s}(Y)]-\textbf{s}(y)\right)$

where $\Sigma_{\theta^{(n)}}\left(\textbf{s}(Y)\right)$ is the covariance
matrix of the statistics vector $\textbf{s}(Y)$ at iteration step $(n)$. The
quantities $\Sigma_{\theta^{(n)}}\left(\textbf{s}(Y)\right)$ and
$\mathbb{E}_{\theta^{(n)}}[\textbf{s}(Y)]$ are not directly computable but,
assuming an effective sampling algorithm, many variations are possible for
parameter estimation. The approach of [8], based on [13], is to use a sample
from the distribution with parameter $\theta^{(t)}$ to get an approximation of
the moment equation in a neighbourhood of $\theta^{(t)}$ and solve it to get
the new estimate $\theta^{(t+1)}$. An alternative approach, given by [41], is
to make use of the [38] stochastic iterative algorithm, and obtain an
approximate value for $\Sigma_{\theta^{(n)}}\left(\textbf{s}(Y)\right)$ and
$\mathbb{E}_{\theta^{(n)}}[\textbf{s}(Y)]$ through simulation:

$\displaystyle\mathbb{E}_{\theta^{(n)}}[\textbf{s}(Y)]$
$\displaystyle\approx\frac{1}{M}\sum_{i=1}^{M}\textbf{s}(y_{i})=\bar{\textbf{s}}$
$\displaystyle\Sigma_{\theta^{(n)}}\left(\textbf{s}(Y)\right)$
$\displaystyle\approx\frac{1}{M}\sum_{i=1}^{M}\left(\textbf{s}(y_{i})-\bar{\textbf{s}}\right)\left(\textbf{s}(y_{i})-\bar{\textbf{s}}\right)^{\top}$

The approximation for the covariance, however, can lead to singularities in
the iteration process when $M$ is small. On the other hand, due to the
sampling process being time-consuming, it would be preferable to not have to
rely on a large number of simulations in every step. [41] tackles the issue by
estimating $D=\Sigma_{\theta^{(1)}}\left(\textbf{s}(Y)\right)$ only for
$\theta^{(1)}$ and keeping it constant throughout the entire iteration. Modern
works on stochastic optimization (see [6], [26], [47]) show that the
evaluation of second-order derivatives of the objective function is not
needed, and convergence can be achieved by choosing an appropriate scalar
learning rate sequence. We direct the reader to [6] for additional details and
limit ourselves to using the suggested practical sequence of $\eta_{t}\sim
1/\sqrt{t}$. The resulting algorithm for ML estimation in ERGMs is sketched in
Algorithm 1 below.

Algorithm 1 Simulation-based SGD maximum likelihood estimation

  * •

Assume an observed network $y$ and an ERGM with chosen model statistics
$\mathbf{s}(\hat{y})$. Set $\theta=\theta_{0}$, $M\geq 0$ and a sequence of
learning rates $\eta_{t}>0$, such that $\eta_{t}=O(t^{-\frac{1}{2}})$.

  * •

For $t\geq 0$, cycle until convergence:

    * –

Using algorithm 1, draw an $M$-size sample $\\{y_{i}^{t},\dots,y_{M}^{t}\\}$
from the model, with parameters $\theta=\theta_{t}$;

    * –

Evaluate
$\Delta_{t}=\frac{1}{M}\sum_{i=1}^{M}\mathbf{s}(\hat{y})-\mathbf{s}(y_{i}^{t})$;

    * –

Update $\theta_{t+1}=\theta_{t}+\eta_{t}\cdot\Delta_{t}$

###  2.2 Penalization

Let us now extend the above to LASSO estimation. To incorporate LASSO approach
to penalized maximum likelihood estimator in ERGMs we add a penalty term to
the log likelihood function and obtain the penalized log likelihood

$\ell_{pen}(Y=y;\theta,\lambda)=\textbf{s(y)}^{\top}\theta-\kappa(\theta)-\lambda\cdot\|\theta^{\prime}\|_{1}.$
(3)

It needs to be noted that the edges statistic takes the role of the intercept
in the case of ERGMs and as such its parameter should remain unpenalized. With
this goal in mind, in (3), $\theta^{\prime}_{j}=\theta_{j}$ for each $j\neq
1$, while $\theta^{\prime}_{1}=0$. Due to the penalizing term having a
singularity at $\theta=0$, simply applying the Stochastic Gradient Descent
(SGD) iteration from Algorithm 1 to the log-likelihood in (3) would result in
excessive oscillations when approaching any $\theta$ having one or more
components $\theta_{j}=0$. Since the goal of the LASSO approach is exactly to
select a parameter vector containing null components, the update rule needs to
be adapted.

The research output on gradient-based optimization methods applied to
penalized functions is very vast. An overview is given in [1] and a selection
of elaborate approaches can be found, for example, in [10] and [29]. The
classic approach, which we follow here, is the method of subgradients. A
subgradient of a function $f(\textbf{x})$ at point $x_{0}$ is defined as a
vector v such that

$f(\mathbf{x})-f(\mathbf{x}_{0})\geq\mathbf{v}\left(\mathbf{x}-\mathbf{x}_{0}\right).$

This definition is a generalization of the concept of gradient and is
applicable in cases where the latter is not defined, such as the LASSO
penalizing term. In the specific case we are interested in, a suitable choice
of subgradient for the function $f(\theta)=\lambda\cdot\|\theta\|_{1}$ is

$\partial f(\theta)=\lambda\cdot sign(\theta)$

where $sign(\theta)$ is equal to $+1$ for positive components of $\theta$,
$-1$ for negative components and $0$ for null components.

The penalized maximum likelihood estimation algorithm we propose can be summed
as shown in Algorithm 2. Given an observed network, we choose a large set of
candidate model statistics, a penalizing factor, a starting parameter vector
and an appropriate learning rate sequence. At each step, a sample of networks
is drawn from the penalized likelihood with the current parameter, which is
used to obtain an estimation of the subgradient of the log-likelihood.
Finally, the parameter vector is updated by SGD, using the appropriate
learning rate and the subgradient estimation previously obtained.

Algorithm 2 LASSO variable selector

  * •

Assume an observed network $y$ and an ERGM with chosen model statistics
$\mathbf{s}(y)$. Set $\theta=\theta_{0}$, $M\geq 0$, penalization term
$\lambda>0$ and a sequence of learning rates $\eta_{t}>0$, such that
$\eta_{t}=O(t^{-\frac{1}{2}})$.

  * •

For $t\geq 0$, cycle until convergence:

    * –

Draw an $M$-size sample $\\{y_{i}^{t},\dots,y_{M}^{t}\\}$ from the model, with
parameters $\theta=\theta_{t}$;

    * –

Evaluate
$\Delta_{t}=\frac{1}{M}\sum_{i=1}^{M}\mathbf{s}(y)-\mathbf{s}(y_{i}^{t})$;

    * –

Update $\theta_{t+1}=\theta_{t}+\eta_{t}\cdot\left(\Delta_{t}-\lambda\cdot
sign\left(\theta_{t}\right)\right)$

  * •

Obtain converged vector of parameters $\hat{\theta}_{pen}(\lambda)$

  * •

Select the vector of statistics $\\{s_{i}\\}$ with
$\hat{\theta}_{pen}^{(i)}(\lambda)\neq 0$.

The vector of statistics $s_{i}(y)$ obtained through Algorithm 2 depends
strongly on the choice of the penalizing term $\lambda$. The higher the degree
of penalization, the lower the number of variables selected.

###  2.3 Variable Importance

Since the set of parameter estimates obtained by Algorithm 2 is penalized
(and, as a consequence, biased), it is not recommended to use the output as
final model estimates. However, we argue that relevant information is
contained in the penalizing term itself and it can be used as a factor that
determines the importance of a model statistic. Along this line of arguments
we define the importance score for the statistic $s_{r}(y)$ as the maximum
value of $\lambda$ for which the penalized estimated parameter
$\hat{\theta}_{pen}^{(i)}(\lambda)$ is non-zero:

$R_{r}=\max\\{\lambda\ |\ \hat{\theta}_{pen}^{(i)}(\lambda)\neq 0\\}.$

The larger the value $R_{r}$, the more penalization is needed to set the
parameter estimate to zero and hence the more important is the corresponding
statistic $s_{r}(y)$. The suggested use for the LASSO selector is therefore
the following:

  * •

Run the algorithm for different values of $\lambda$, in order to obtain a
ranking of variables based on their importance score;

  * •

Choose a threshold value $\lambda_{threshold}$ for the importance score, so
that all variable with $R_{t}>\lambda_{threshold}$ are included in the model.

  * •

Fit an unpenalized ERGM, using the variables selected in the previous step.

It is worth noting that above, no indication is given on how to choose the
threshold value $\lambda_{threshold}$. One possibility is to make use of the
Akaike Information Criterion (AIC). To do so we consider $R_{t}$ as variable
ranking. The choice of the threshold that dictates which variables to include
is then obtained such that adding the next variable in the ranked order of the
statistics would result in a lower AIC value. Another approach is to follow
p-values and choose the threshold such that adding a single new statistic
based on the ranking obtained through $R_{r}$ would result in non-significant
estimates. Whatever method one uses to determine the threshold
$\lambda_{threshold}$, the importance ranking $R_{t}$ based on the LASSO
estimate reduces the set of candidate models for model selection considerably.
We can conclude that the main appeal of the LASSO selector is to provide a
ranking of statistics importance.

###  2.4 Standardization of network features

Even in classical i.i.d. models, the use of penalized estimation often
requires to standardize the component variables, due to their different range
of possible values: the size of the effect of a given variable, namely its
parameter, depends on the range of values taken by the variable itself. When
penalization is applied to a set of parameters, it is necessary for their
numerical values to be comparable. A commonly used way to deal with this
problem is to standardize the variables, i.e. rescale the range of values
taken by each variable to a fixed interval.

When it comes to ERGMs the issue becomes slightly more tricky. Since we are
provided with a single network observation, for variables such as structural
statistics we are not provided with the sample standard deviation in the usual
way. Instead, a larger sample is to be generated from a model of the observed
network. In our case, we opted for a simple Erdős-Rényi model, as the goal
here is not precision in the fit, but rather an estimation of the magnitude of
structural variables per se. The estimation is provided by the observed values
in the sample: given a component $s_{r}(y)$ of the network statistics vector,
and a sample $\\{y_{1},\dots,y_{M}\\}$ of networks simulated from an Erdős-
Rényi model, the scaled statistic will be:

$s_{r,scaled}(y)=SD\left(s_{r}(y_{1}),\dots,s_{r}(y_{M})\right)^{-1}s_{r}(y)$
(4)

where the standard deviation is defined as usual for an $M$-sized sample as

$SD\left(s_{r}(y_{1}),\dots,s_{r}(y_{M})\right)=\sqrt{\frac{1}{M-1}\sum_{j=1}^{M}\left(s_{r}(y_{j})-\bar{s}_{r}\right)^{2}}.$

In Section 2.1 we discussed the problems with variance estimation when it
comes to simulations: in particular we discussed how, to avoid singularities,
a large sample size is needed, and how a large sample size can put a strain on
running times when samples need to be drawn in each step of an iterative
algorithm. In this case this is not a limiting issue: the variance estimation
used for standardization is performed one-off at the start of the cycle and
not repeated in each step of the iteration.

##  3 Simulation Studies

###  3.1 Endogeneous Network Statistics

Before presenting the estimation and model selection strategy on real-world
networks we want to explore the procedure in simulations. These simulations
allow us to study whether LASSO can be used for model selection in ERGMs. To
do so, we run simulations with known model setups which are described
subsequently. Traditionally the triangle count has been used as a meaningful
statistic to measure triadic closure for networks, while the 2-star count as a
meaningful statistic in order to measure triadic non-closure. However,
literature has shown ([16]) that these two particular statistics increase the
chances for the final model to be near-degenerate. This means that the
distribution resulting from the model will likely have a double peak shape and
assign unreasonably high probabilities to near-empty or near-complete
networks. A more favourable choice of variables uses a geometrically weighted
sum of edgewise shared partners (gwesp) as a substitute for the triangle
count, and a geometrically weighted sum of non-edgewise shared partners
(gwnsp) as a substitute for the 2-star count.

These two statistics are formally defined as follows:

$\displaystyle
s_{gwesp}(y,\alpha)=e^{\alpha}\sum_{1}^{n-2}\left(1-\left(1-e^{-\alpha}\right)^{i}\right)esp(i)$
$\displaystyle
s_{gwnsp}(y,\alpha)=e^{\alpha}\sum_{1}^{n-2}\left(1-\left(1-e^{-\alpha}\right)^{i}\right)nsp(i).$

Here $esp(i)$ and $nsp(i)$ count the number of edges (or non-edges,
respectively) that share exactly $i$ nodes as neighbours and $\alpha$ is a
decay parameter which can be set manually, as in our case, or estimated
through the use of a curved exponential family distribution ([21]).

This simulation study involves 6 different setups: in each setup a sample of
networks is drawn from a known ERGM distribution and a run through the LASSO
variable selector, in order to test if the procedure is able to select the
variable used in the original distribution. In every setup, the size of the
sample is 20.

The first three setups are performed on networks featuring $N=50$ nodes, the
last three on networks featuring $N=150$ nodes. In setup number 1 networks are
drawn from an ERGM distribution using edges and gwesp as variables, with
coefficients $\theta_{edges}=-3.5$ and $\theta_{gwesp}=1.0$. In setup number 2
the sample is drawn from an ERGM distribution using edges and gwnsp as
variables, with coefficients $\theta_{edges}=-2.0$ and $\theta_{gwesp}=0.5$.
In setup number 3 the sample is drawn from an Erdős-Rényi model, equivalent to
an ERGM using only edges as a variable, with coefficient
$\theta_{edges}=-1.5$. Setups 4 to 6 are a repetition of the first three,
where only the number of nodes is changed. In all setups the starting set for
the selecting procedure contains gwesp, gwnsp, and gwdegree as potential
candidates.

The study is conducted as follows: in each setup and for each sampled network,
the $\theta$ coefficients are evaluated for different magnitudes of the
penalization parameter $\lambda$. Afterwards, importance score $R_{r}$, i.e.
the highest $\lambda$ for which the coefficient is non-zero is recorded. It
must be noted that the importance score obtained depends on the sequence of
penalization that is used. We use the R package ERGM from the Statnet project
for all simulations ([17]).

In Figure 1 we take an example network from setup 1 and 2, respectively, and
show how the coefficients of each variable change according to the magnitude
of the penalization applied.

Setup 1 uses networks generated from an ERGM using gwesp as model statistic
and, as such, gwesp is selected in the first example, with an importance score
of 20. On the other hand, setup 2 uses networks generated from an ERGM using
gwnsp as model statistic. This is again accurately reflected in the second
example, where gwnsp is selected with an importance score of 50.

Figure 1: Example of coefficients path in setup 1 (left), and setup 2 (right).

A summary of the entire simulation setup including a visualisation of the
variability shown in Figure 2. Each plot refers to a simulation setup, where
different colors represent the coefficients of each variable evaluated for a
given penalization $\lambda$. The middle colored line represents the average
of each coefficient, over the sample of the 20 networks used in the setup,
while the shaded area spans between maximum and minimum.

From Figure 2 it is clear that whenever a variable is used in the original
distribution, it is consistently selected as first by the procedure based on
the importance score $R_{t}$, by a very large margin in the case of gwnsp or
by a smaller one in the case of gwesp. In the Erdős-Rényi setups, a
significantly lower penalization is necessary to obtain non-zero parameters.
This indicates a lower importance for all of the statistics, as expected given
that the networks were sampled from a distribution assuming independent edges.

Figure 2: Coefficient path in the different simulation setups, for decreasing
values of the penalization parameter. The colored line represents the average
over the sample, the coloured area covers the maximum and minimum over the
sample.

###  3.2 Endogeneous and Exogeneous Network Statistics

The above simulation study concerns endogenous network statistics. If
additional exogenous quantities are available, a typical research question
when analysing real life social networks is whether the reason behind the tie
formation process is driven endogenously or exogenously, i.e. depending on the
local structure of the ties surrounding the actors or on the actors’ intrinsic
characteristics. Our second simulation study focuses on the capability of
LASSO as a variable selection routine to distinguish between endogenous and
exogenous variables. As done before, we differentiate between 2 setups, one
with $N=50$ sized networks, and another with $N=150$. In each setup a sample
of 20 networks is drawn: each node of each network possesses an intrinsic
(binary) nodal attribute $x(i)$, being the only influencing factor in edge
formation.

This means that all edges are drawn independently with probabilities
$P(Y_{i,j})=0.05,\ 0.15,\ 0.30$ depending if $\chi(i)=\chi(j)=0$, if only one
of $\chi(i)$ and $\chi(j)$ is equal to $1$, or if $\chi(i)=\chi(j)=1$.

The same selection of endogenous statistics is used, with the addition of an
exogenous statistic summarizing the effect of the attribute. The exogenous
statistic, implemented through the nodefactor term in the ERGM package, counts
the number of occurrences of a node with a given attribute value in the
network edges and is formally defined as

$s_{attr}(y)=\sum_{i}^{n}\sum_{j<i}y_{i,j}(x_{attr}(i)+x_{attr}(j))$

where $x_{attr}(i)$ is equal to 1 if node $i$ has its binary attribute set to
1 and equal to 0 otherwise.

Figure 4 once again shows the coefficient path for different values of
$\lambda$: the exogenous attribute as the main driver for tie formation is
clearly visible from the coefficient path of nodefactor.Attribute.

Figure 3: Example network from the second simulation study, with N=50 nodes.
Nodes with the attribute yes are set to have a higher probability of forming
ties. Figure 4: Coefficient path in the two simulation setups using a nodal
attribute, for decreasing values of the penalization parameter. Network size
N=150 on the left, N=50 on the right.

##  4 Real data applications

###  4.1 London street gang

The first real network we investigate is the London Street Gang network (5),
first studied by [14] and re-used by [4] in the context of ERGM literature.
The 54 nodes in the network represent each a known member of a criminal gang
based in London between 2006 and 2009. An undirected tie between two nodes is
present if the two gang members have been arrested together while committing a
crime. The set of nodal covariates includes age and birthplace of the subject,
total number of arrests and convictions, and binary variables representing
whether or not the subject has been incarcerated before, is holding a legal
residence permit and has been active in the music industry.

All gang members are ethnically identifiable as black, but their birthplace
splits them into 4 different cultural subgroups: West African, Jamaican, UK
and Somali. Part of the goal of the original study by [14] was to show that
the cultural group had a major influence on how ties are formed between gang
members.

Our goal is to utilize our variable selection procedure to obtain a ranking of
importance between the different variables and get an insight into whether
this particular network is exogenously or endogenously driven.

Nodal attributes are encoded using the standard terms used in the ERGM R
package: for categorical attributes one network statistic is added for each
level minus one, counting the number of times a node with the set level
appears in an edge; for numerical attributes, one single statistic is added,
equal to the sum of the attribute value of the connected nodes, overall
network edges; the node matching statistics of a given categorical attribute
simply count all edges of the network connecting nodes with the same level.
More precisely:

$\displaystyle s_{numerical}(y,x)$
$\displaystyle=\sum_{i,j}y_{i,j}(x(i)+x(j))$ $\displaystyle
s_{categorical}(y,x=1)$
$\displaystyle=\sum_{i,j}y_{i,j}(\mathds{1}(x_{i}=1)+\mathds{1}(x_{j}=1))$
$\displaystyle s_{matching}(y,x)$
$\displaystyle=\sum_{i,j}y_{i,j}(\mathds{1}(x_{i}=x_{j}))$

Figure 5: London gangs network. Birthplace encoding: 1 = West Africa, 2 =
Jamaica, 3 = UK, 4 = Somalia.

The nodal covariates we include are all the aforementioned individual
information on the subjects, while the set of structural statistics, following
what was done in the simulation studies, includes the gwesp and gwnsp
statistic to measure triadic effects, the gwdegree statistic to measure
degree-based effects.

Figure 6: Coefficient values given different penalization parameters for the
London Street Gang network.

We show the coefficient paths in Figure 6. The results give a valuable insight
on tie formation: the network is mostly structure-driven, as the edge- and
non-edge-wise shared partners statistics are selected early on. The positive
and negative signs for the gwesp and gwnsp coefficients, respectively, are a
strong indicator of triadic closure, namely the tendency for the edges to form
closed triangular structures.

Of all the nodal covariates, the number of arrests is the first being selected
based on the importance score $R_{r}$ and the only one with an importance
value comparable to the structural variables. Other nodal covariates require a
lower penalization to be included in the model and thus can be considered of
lower importance to the tie formation process. Among the dummy variables
encoding the birthplace covariate, it’s worthwhile noting the higher
importance value for Birthplace.3 (UK born). This shows a higher relevance in
the tie forming process of being born in the UK, compared to different
birthplaces. In addition to the higher relevance, the negative value for the
coefficient highlights a negative influence or, in other words, a higher
probability for UK-born gang members to occupy marginal roles in the gang. The
last two variables to get selected are nodematch.Birthplace, indicating a
tendency towards associating with gang members born in the same region, and
nodefactor.Music.

The last step is to fit a non-penalized ERGM making use of the variables
mentioned above. Since adding the next ranked variable to the model causes an
increase of AIC score, our final (unpenalized) model is the one shown in Table
1.

Variable |  Coefficient |  SE |  P-value
---|---|---|---
edges |  -3.7072 |  0.4003 |  < 0.0001
gwnsp.fixed.0.5 |  -0.1179 |  0.0420 |  0.0049
gwesp.fixed.0.5 |  1.2011 |  0.1693 |  < 0.0001
nodefactor.UK.1 |  -0.2365 |  0.0711 |  0.0009
nodecov.Arrests |  0.1741 |  0.0525 |  0.0009
nodematch.Birthplace |  0.8331 |  0.1409 |  < 0.0001
nodefactor.Music.1 |  0.2533 |  0.1176 |  0.0313
Table 1: Estimated coefficents, standard errors and P-values of the final
proposed model for the London Street Gang network.

###  4.2 Coworker relationship in a corporate law partnership

The second dataset we apply our method to is the strong coworker network used
in [30] and later on featured in [42]. The 71 nodes in the network represent
lawyers of Northeastern US corporate law firm, referred to as SG&R, operative
in New England between 1988-1991. A tie between two members of the firm exists
if in the past year, the two lawyers either worked together on a case or an
administrative task or exchanged work products. A number of nodal attributes
are available, such as gender, age, years spent within the firm, status in the
firm (i.e. partner or associate), the office they are based in (Boston,
Hartford or Providence), the type of their practice (corporate or litigation)
and the law school they went to.

Once again, our goal is to use variable selection to rank the different
structural and exogenous variables. The coefficient paths are depicted in
Figure 7.

Figure 7: Coefficient values given different penalization parameters for the
Law Firm Coworkers network.

As opposed to the previous example, in this case, two node attributes, namely
the match between offices and the match between practices, are strongly
detected as the main drivers for tie formation. As expected, two workers
sharing the same area of practice or the same office have a very high chance
to collaborate.

Next selected is a group of variables with approximately the same importance
score. The structural variable gwesp, as well as the Status of the lawyer and
the match between Status. Both, the negative value for the nodematch.Status
coefficient and the positive value for nodefactor.Status.Partner highlights
the tendency for high-profile partners to not collaborate with each other, but
at the same time to form ties with lower-ranked lawyers. The nodal factor
office also belongs to this cluster.

The rest of the variables in the importance ranking are closely clustered
together with very little difference in score, the highest scoring one being
the match between law schools. We decided to draw the cutting line here, as
adding another variable ends up increasing the AIC score of the final model.

The proposed model, making use of the selected variables, is summarized in
Table 2.

Variable |  Coefficient |  SE |  P-value
---|---|---|---
edges |  -6.2482 |  0.3115 |  < 0.0001
nodematch.Practice |  1.3862 |  0.1317 |  < 0.0001
nodematch.Office |  1.8086 |  0.1646 |  < 0.0001
nodefactor.Status.Partner |  0.5655 |  0.1001 |  < 0.0001
gwesp.fixed.0.5 |  0.9608 |  0.1654 |  < 0.0001
nodematch.Status |  -0.7862 |  0.1495 |  < 0.0001
nodefactor.Office.Hartford |  0.4277 |  0.0937 |  < 0.0001
nodefactor.Office.Providence |  0.9954 |  0.1870 |  < 0.0001
nodematch.Law school |  0.3249 |  0.1277 |  0.011
Table 2: Estimated coefficents, standard errors and P-values of the final
proposed model for the Law Firm Coworkers network.

##  5 Discussion

The usage of a penalized likelihood provides for an effective variable
selection procedure in the framework of Exponential Random Graph models.
Refitting the penalized model for decreasing values of the penalization
parameter gives an order of importance for the initial set of variables, where
one is less important than another if a smaller penalization is needed to
obtain a non-zero coefficient. In other words, the value of the smallest
penalization parameter for which a given variable has a non-zero coefficient
can be used as an importance score for the variable. As noted before, the
coefficients obtained are biased and for this reason, we do not use them for
the final model but, instead, refit an unpenalized model using the selected
variables.

The importance score can be derived for both structural statistics and nodal
attributes making it possible to deduce if the main driver for tie formation
in the network is endogenous, in case the structural statistics rank high, or
exogenous, if nodal covariates rank higher.

Model selection is achieved through setting a cutoff value for the importance
score, and choosing the variables with higher scores than the cutoff. No
indication is given on the choice of the cutoff, which can be tailored to the
specific needs of the case, or be done following other model quality metrics,
such as AIC or variable p-values.

In the present paper we only deal with static, undirected networks, but as
ERGMs or ERG-derived models are widely used to analyse both directed networks
and networks that evolve through time. A possible direction for additional
development and future work will be to extend the usability of the LASSO based
variable selector to both directed and temporal ERGMs.

## References

  * [1] Francis R. Bach, Rodolphe Jenatton, Julien Mairal and Guillaume Obozinski  “Optimization with Sparsity-Inducing Penalties”  In _CoRR_ abs/1108.0775, 2011  arXiv: http://arxiv.org/abs/1108.0775
  * [2] Brenda Betancourt  “Modeling and Prediction of Time Series of Directed Binary Networks”, 2015
  * [3] Brenda Betancourt, Abel Rodríguez and Naomi Boyd  “Bayesian fused lasso regression for dynamic binary networks”  In _Journal of Computational and Graphical Statistics_ 26.4 Taylor & Francis, 2017, pp. 840–850
  * [4] Bart Blackburn and Mark S. Handcock  “Practical Network Modeling via Tapered Exponential-Family Random Graph Models”  In _Journal of Computational and Graphical Statistics_ 32.2 Taylor & Francis, 2023, pp. 388–401  DOI: 10.1080/10618600.2022.2116444
  * [5] Stephen P Borgatti, Martin G Everett and Jeffrey C Johnson  “Analyzing social networks”  Sage, 2018
  * [6] Léon Bottou  “Stochastic gradient descent tricks”  In _Neural networks: Tricks of the trade_ , 2012, pp. 421–436
  * [7] Skyler J Cranmer, Bruce A Desmarais and Jason W Morgan  “Inferential network analysis”  Cambridge University Press, 2020
  * [8] B. Crouch, S. Wasserman and F. Trachtenberg  “Markov Chain Monte Carlo maximum likelihood estimation for $p^{*}$ social network models.”  In _Paper presented at the Sunbelt XVIII and Fifth European International Social Networks Conference, Sitges (Spain)_ , 1998
  * [9] Mark S. David R.  “Inference in Curved Exponential Family Models for Networks”  In _Journal of Computational and Graphical Statistics_ 15:3, 2006, pp. 565–583
  * [10] John Duchi and Yoram Singer  “Efficient online and batch learning using forward backward splitting.”  In _Journal of Machine Learning Research_ , 2009, pp. 2899–2934
  * [11] Ove Frank and David Strauss  “Markov Graphs”  In _Journal of the American Statistical Association_ 81.395, 1986, pp. 832–842
  * [12] W. Fu and K. Knight  “Asymptotics for lasso-type estimators”  In _Ann. Statist._ 28 (5), 2000, pp. 1356–1378
  * [13] C.J. Geyer and E.A. Thompson  “Constrained Monte Carlo Maximum Likelihood for Dependent Data”  In _Journal of the Royal Statistical Society_ B 54, 1992, pp. 657–699
  * [14] Thomas Grund and James Densley  “Ethnic Homophily and Triad Closure: Mapping Internal Gang Structure Using Exponential Random Graph Models”  In _Journal of Contemporary Criminal Justice_ 31, 2015, pp. 354–370  DOI: 10.1177/1043986214553377
  * [15] Dan Han, Vicky Modisette, Melinda Forthofer and Rajib Paul  “Hierarchical Bayesian Adaptive Lasso Methods on Exponential Random Graph Models”  In _PREPRINT (Version 1) available at Research Square_ , December 2023  DOI: https://doi.org/10.21203/rs.3.rs-3669311/v1
  * [16] Mark Handcock  “Assessing degeneracy in statistical models of social networks”  In _Journal of the American Statistical Association_ 76, 2003, pp. 33–50
  * [17] Mark S. Handcock et al.  “ergm: Fit, Simulate and Diagnose Exponential-Family Models for Networks” R package version 3.11.0, 2021
  * [18] Jenine K Harris  “An introduction to exponential random graph modeling”  Sage Publications, 2013
  * [19] Paul W Holland and Samuel Leinhardt  “An exponential family of probability distributions for directed graphs”  In _Journal of the American Statistical Association_ 76.373 Taylor & Francis, 1981, pp. 33–50
  * [20] R.. Hummel, D.. Hunter and M.. Handcock  “Improving simulation based algorithms for fitting ERGMs”  In _Journal of Comp. and Graph. Statistics_ 21.4, 2012, pp. 920–939
  * [21] D.. Hunter  “Curved Exponential Family Models for Social Networks.”  In _Social Networks._ 29.2, 2007, pp. 216–230
  * [22] David R Hunter, Steven M Goodreau and Mark S Handcock  “Goodness of fit of social network models”  In _Journal of the american statistical association_ 103.481 Taylor & Francis, 2008, pp. 248–258
  * [23] David R Hunter et al.  “ergm: A package to fit, simulate and diagnose exponential-family models for networks”  In _Journal of statistical software_ 24.3 NIH Public Access, 2008, pp. nihpa54860
  * [24] David R. Hunter, P.N. Krivitsky and Michael Schweinberger  “Computational statistical methods for social network analysis”  In _Journal of Comp. and Graph. Statistics_ 21.4, 2012, pp. 856–882
  * [25] Sevag Kevork and Göran Kauermann  “Bipartite exponential random graph models with nodal random effects”  In _Social Networks_ 70 Elsevier, 2022, pp. 90–99
  * [26] Diederik P. Kingma and Jimmy Ba  “Adam: A Method for Stochastic Optimization”, 2017  arXiv:1412.6980 [cs.LG]
  * [27] E.. Kolaczyk  “Statistical Analysis of Network Data”  Springer, 2009
  * [28] E.. Kolaczyk and G. Csardi  “Statistical Analysis of Network Data with R”, 2014
  * [29] Lihong Li Langford John and Tong Zhang  “Sparse online learning via truncated gradient”  In _Journal of Machine Learning Research_ , 2009, pp. 777–801
  * [30] Emmanuel Lazega  “The Collegial Phenomenon: The Social Mechanisms of Cooperation Among Peers in a Corporate Law Partnership”  In _Oxford University Press_ , 2001
  * [31] Philip Leifeld and Skyler J Cranmer  “A theoretical and empirical comparison of the temporal exponential random graph model and the stochastic actor-oriented model”  In _Network science_ 7.1 Cambridge University Press, 2019, pp. 20–51
  * [32] Philip Leifeld, Skyler J Cranmer and Bruce A Desmarais  “Temporal exponential random graph models with btergm: Estimation and bootstrap confidence intervals”  In _Journal of Statistical Software_ 83.6 University of California, Los Angeles, 2018
  * [33] “Exponential Random Graph Models for Social Networks: Theory, Methods, and Applications”, Structural Analysis in the Social Sciences  Cambridge University Press, 2012  DOI: 10.1017/CBO9780511894701
  * [34] Vicki Modisette  “Penalized Bayesian exponential random graph models”  University of Louisville, Electronic ThesesDissertations. Paper 4139, 2023
  * [35] Martina Morris, Mark S. Handcock and David R. Hunter  “Specification of Exponential-Family Random Graph Models: Terms and Computational Aspects”  In _Journal of Statistical Software, 24(4), 1–24_ 24(4), 2008, pp. 1–24
  * [36] Trevor Park and George Casella  “The bayesian lasso”  In _Journal of the American Statistical Association_ 103.482 Taylor & Francis, 2008, pp. 681–686
  * [37] Sa Ren  “Inferring network structures using hierarchical exponential random graph models”, 2023
  * [38] H. Robbins and S. Monro  “A stochastic approximation method”  In _Annals of Mathematical Statistics_ 22, 1951, pp. 400–407
  * [39] Fadil Santosa and William W Symes  “Linear inversion of band-limited reflection seismograms”  In _SIAM journal on scientific and statistical computing_ 7.4 SIAM, 1986, pp. 1307–1330
  * [40] Michael Schweinberger and Pamela Luna  “HERGM: Hierarchical exponential-family random graph models”  In _Journal of Statistical Software_ 85, 2018, pp. 1–39
  * [41] Tom A.. Snijders  “Markov chain Monte Carlo estimation of exponential random graph models”  In _Journal of Social Structure_ 3.2, 2002, pp. 1–40
  * [42] Tom A.. Snijders, Philippa E. Pattison, Garry L. Robins and Mark S. Handcock  “New Specifications for Exponential Random Graph Models”  In _Sociological Methodology_ 36.1, 2006, pp. 99–153
  * [43] Tom AB Snijders  “Statistical models for social networks”  In _Annual review of sociology_ 37 Annual Reviews, 2011, pp. 131–153
  * [44] Stephanie Thiemichen, Nial Friel, Alberto Caimo and Göran Kauermann  “Bayesian exponential random graph models with nodal random effects”  In _Social Networks_ 46 Elsevier, 2016, pp. 11–28
  * [45] Robert Tibshirani  “Regression shrinkage and selection via the lasso”  In _Journal of the Royal Statistical Society Series B: Statistical Methodology_ 58.1 Oxford University Press, 1996, pp. 267–288
  * [46] Helal El-Zaatari, Fei Yu and Michael R. Kosorok  “Stochastic Step-wise Feature Selection for Exponential Random Graph Models (ERGMs)”  In _https://doi.org/10.48550/arXiv.2307.12862_ , 2023
  * [47] Matthew D. Zeiler  “An adaptive learning rate method”  In _arXiv:1212.5701, https://doi.org/10.48550/arXiv.1212.5701_ , 2012
  * [48] Peng Zhao and Bin Yu  “On Model Selection Consistency of Lasso”  In _Journal of Machine Learning Research_ 7.90, 2006, pp. 2541–2563  URL: http://jmlr.org/papers/v7/zhao06a.html
  * [49] Peng Zhao and Bin Yu  “On model selection consistency of Lasso”  In _The Journal of Machine Learning Research_ 7 JMLR. org, 2006, pp. 2541–2563


The LHCb collaboration

# First observation of $C\\!P$ violation in the decays of $B^{0}_{s}$ mesons

R. Aaij40, C. Abellan Beteta35,n, B. Adeva36, M. Adinolfi45, C. Adrover6, A.
Affolder51, Z. Ajaltouni5, J. Albrecht9, F. Alessio37, M. Alexander50, S.
Ali40, G. Alkhazov29, P. Alvarez Cartelle36, A.A. Alves Jr24,37, S. Amato2, S.
Amerio21, Y. Amhis7, L. Anderlini17,f, J. Anderson39, R. Andreassen56, R.B.
Appleby53, O. Aquines Gutierrez10, F. Archilli18, A. Artamonov 34, M.
Artuso58, E. Aslanides6, G. Auriemma24,m, S. Bachmann11, J.J. Back47, C.
Baesso59, V. Balagura30, W. Baldini16, R.J. Barlow53, C. Barschel37, S.
Barsuk7, W. Barter46, Th. Bauer40, A. Bay38, J. Beddow50, F. Bedeschi22, I.
Bediaga1, S. Belogurov30, K. Belous34, I. Belyaev30, E. Ben-Haim8, G.
Bencivenni18, S. Benson49, J. Benton45, A. Berezhnoy31, R. Bernet39, M.-O.
Bettler46, M. van Beuzekom40, A. Bien11, S. Bifani44, T. Bird53, A.
Bizzeti17,h, P.M. Bjørnstad53, T. Blake37, F. Blanc38, J. Blouw11, S. Blusk58,
V. Bocci24, A. Bondar33, N. Bondar29, W. Bonivento15, S. Borghi53, A.
Borgia58, T.J.V. Bowcock51, E. Bowen39, C. Bozzi16, T. Brambach9, J. van den
Brand41, J. Bressieux38, D. Brett53, M. Britsch10, T. Britton58, N.H. Brook45,
H. Brown51, I. Burducea28, A. Bursche39, G. Busetto21,q, J. Buytaert37, S.
Cadeddu15, O. Callot7, M. Calvi20,j, M. Calvo Gomez35,n, A. Camboni35, P.
Campana18,37, D. Campora Perez37, A. Carbone14,c, G. Carboni23,k, R.
Cardinale19,i, A. Cardini15, H. Carranza-Mejia49, L. Carson52, K. Carvalho
Akiba2, G. Casse51, L. Castillo Garcia37, M. Cattaneo37, Ch. Cauet9, M.
Charles54, Ph. Charpentier37, P. Chen3,38, N. Chiapolini39, M. Chrzaszcz 25,
K. Ciba37, X. Cid Vidal37, G. Ciezarek52, P.E.L. Clarke49, M. Clemencic37,
H.V. Cliff46, J. Closier37, C. Coca28, V. Coco40, J. Cogan6, E. Cogneras5, P.
Collins37, A. Comerma-Montells35, A. Contu15,37, A. Cook45, M. Coombes45, S.
Coquereau8, G. Corti37, B. Couturier37, G.A. Cowan49, D.C. Craik47, S.
Cunliffe52, R. Currie49, C. D’Ambrosio37, P. David8, P.N.Y. David40, A.
Davis56, I. De Bonis4, K. De Bruyn40, S. De Capua53, M. De Cian39, J.M. De
Miranda1, L. De Paula2, W. De Silva56, P. De Simone18, D. Decamp4, M.
Deckenhoff9, L. Del Buono8, N. Déléage4, D. Derkach14, O. Deschamps5, F.
Dettori41, A. Di Canto11, F. Di Ruscio23,k, H. Dijkstra37, M. Dogaru28, S.
Donleavy51, F. Dordei11, A. Dosil Suárez36, D. Dossett47, A. Dovbnya42, F.
Dupertuis38, R. Dzhelyadin34, A. Dziurda25, A. Dzyuba29, S. Easo48,37, U.
Egede52, V. Egorychev30, S. Eidelman33, D. van Eijk40, S. Eisenhardt49, U.
Eitschberger9, R. Ekelhof9, L. Eklund50,37, I. El Rifai5, Ch. Elsasser39, D.
Elsby44, A. Falabella14,e, C. Färber11, G. Fardell49, C. Farinelli40, S.
Farry12, V. Fave38, D. Ferguson49, V. Fernandez Albor36, F. Ferreira
Rodrigues1, M. Ferro-Luzzi37, S. Filippov32, M. Fiore16, C. Fitzpatrick37, M.
Fontana10, F. Fontanelli19,i, R. Forty37, O. Francisco2, M. Frank37, C.
Frei37, M. Frosini17,f, S. Furcas20, E. Furfaro23,k, A. Gallas Torreira36, D.
Galli14,c, M. Gandelman2, P. Gandini58, Y. Gao3, J. Garofoli58, P. Garosi53,
J. Garra Tico46, L. Garrido35, C. Gaspar37, R. Gauld54, E. Gersabeck11, M.
Gersabeck53, T. Gershon47,37, Ph. Ghez4, V. Gibson46, V.V. Gligorov37, C.
Göbel59, D. Golubkov30, A. Golutvin52,30,37, A. Gomes2, H. Gordon54, M.
Grabalosa Gándara5, R. Graciani Diaz35, L.A. Granado Cardoso37, E. Graugés35,
G. Graziani17, A. Grecu28, E. Greening54, S. Gregson46, P. Griffith44, O.
Grünberg60, B. Gui58, E. Gushchin32, Yu. Guz34,37, T. Gys37, C.
Hadjivasiliou58, G. Haefeli38, C. Haen37, S.C. Haines46, S. Hall52, T.
Hampson45, S. Hansmann-Menzemer11, N. Harnew54, S.T. Harnew45, J. Harrison53,
T. Hartmann60, J. He37, V. Heijne40, K. Hennessy51, P. Henrard5, J.A. Hernando
Morata36, E. van Herwijnen37, E. Hicks51, D. Hill54, M. Hoballah5, C.
Hombach53, P. Hopchev4, W. Hulsbergen40, P. Hunt54, T. Huse51, N. Hussain54,
D. Hutchcroft51, D. Hynds50, V. Iakovenko43, M. Idzik26, P. Ilten12, R.
Jacobsson37, A. Jaeger11, E. Jans40, P. Jaton38, A. Jawahery57, F. Jing3, M.
John54, D. Johnson54, C.R. Jones46, C. Joram37, B. Jost37, M. Kaballo9, S.
Kandybei42, M. Karacson37, T.M. Karbach37, I.R. Kenyon44, U. Kerzel37, T.
Ketel41, A. Keune38, B. Khanji20, O. Kochebina7, I. Komarov38, R.F. Koopman41,
P. Koppenburg40, M. Korolev31, A. Kozlinskiy40, L. Kravchuk32, K. Kreplin11,
M. Kreps47, G. Krocker11, P. Krokovny33, F. Kruse9, M. Kucharczyk20,25,j, V.
Kudryavtsev33, T. Kvaratskheliya30,37, V.N. La Thi38, D. Lacarrere37, G.
Lafferty53, A. Lai15, D. Lambert49, R.W. Lambert41, E. Lanciotti37, G.
Lanfranchi18, C. Langenbruch37, T. Latham47, C. Lazzeroni44, R. Le Gac6, J.
van Leerdam40, J.-P. Lees4, R. Lefèvre5, A. Leflat31, J. Lefrançois7, S.
Leo22, O. Leroy6, T. Lesiak25, B. Leverington11, Y. Li3, L. Li Gioi5, M.
Liles51, R. Lindner37, C. Linn11, B. Liu3, G. Liu37, S. Lohn37, I.
Longstaff50, J.H. Lopes2, E. Lopez Asamar35, N. Lopez-March38, H. Lu3, D.
Lucchesi21,q, J. Luisier38, H. Luo49, F. Machefert7, I.V. Machikhiliyan4,30,
F. Maciuc28, O. Maev29,37, S. Malde54, G. Manca15,d, G. Mancinelli6, U.
Marconi14, R. Märki38, J. Marks11, G. Martellotti24, A. Martens8, L. Martin54,
A. Martín Sánchez7, M. Martinelli40, D. Martinez Santos41, D. Martins Tostes2,
A. Massafferri1, R. Matev37, Z. Mathe37, C. Matteuzzi20, E. Maurice6, A.
Mazurov16,32,37,e, J. McCarthy44, A. McNab53, R. McNulty12, B. Meadows56,54,
F. Meier9, M. Meissner11, M. Merk40, D.A. Milanes8, M.-N. Minard4, J. Molina
Rodriguez59, S. Monteil5, D. Moran53, P. Morawski25, M.J. Morello22,s, R.
Mountain58, I. Mous40, F. Muheim49, K. Müller39, R. Muresan28, B. Muryn26, B.
Muster38, P. Naik45, T. Nakada38, R. Nandakumar48, I. Nasteva1, M. Needham49,
N. Neufeld37, A.D. Nguyen38, T.D. Nguyen38, C. Nguyen-Mau38,p, M. Nicol7, V.
Niess5, R. Niet9, N. Nikitin31, T. Nikodem11, A. Nomerotski54, A. Novoselov34,
A. Oblakowska-Mucha26, V. Obraztsov34, S. Oggero40, S. Ogilvy50, O.
Okhrimenko43, R. Oldeman15,d, M. Orlandea28, J.M. Otalora Goicochea2, P.
Owen52, A. Oyanguren 35,o, B.K. Pal58, A. Palano13,b, M. Palutan18, J.
Panman37, A. Papanestis48, M. Pappagallo50, C. Parkes53, C.J. Parkinson52, G.
Passaleva17, G.D. Patel51, M. Patel52, G.N. Patrick48, C. Patrignani19,i, C.
Pavel-Nicorescu28, A. Pazos Alvarez36, A. Pellegrino40, G. Penso24,l, M. Pepe
Altarelli37, S. Perazzini14,c, D.L. Perego20,j, E. Perez Trigo36, A. Pérez-
Calero Yzquierdo35, P. Perret5, M. Perrin-Terrin6, G. Pessina20, K.
Petridis52, A. Petrolini19,i, A. Phan58, E. Picatoste Olloqui35, B. Pietrzyk4,
T. Pilař47, D. Pinci24, S. Playfer49, M. Plo Casasus36, F. Polci8, G. Polok25,
A. Poluektov47,33, E. Polycarpo2, A. Popov34, D. Popov10, B. Popovici28, C.
Potterat35, A. Powell54, J. Prisciandaro38, V. Pugatch43, A. Puig Navarro38,
G. Punzi22,r, W. Qian4, J.H. Rademacker45, B. Rakotomiaramanana38, M.S.
Rangel2, I. Raniuk42, N. Rauschmayr37, G. Raven41, S. Redford54, M.M. Reid47,
A.C. dos Reis1, S. Ricciardi48, A. Richards52, K. Rinnert51, V. Rives
Molina35, D.A. Roa Romero5, P. Robbe7, E. Rodrigues53, P. Rodriguez Perez36,
S. Roiser37, V. Romanovsky34, A. Romero Vidal36, J. Rouvinet38, T. Ruf37, F.
Ruffini22, H. Ruiz35, P. Ruiz Valls35,o, G. Sabatino24,k, J.J. Saborido
Silva36, N. Sagidova29, P. Sail50, B. Saitta15,d, V. Salustino Guimaraes2, C.
Salzmann39, B. Sanmartin Sedes36, M. Sannino19,i, R. Santacesaria24, C.
Santamarina Rios36, E. Santovetti23,k, M. Sapunov6, A. Sarti18,l, C.
Satriano24,m, A. Satta23, M. Savrie16,e, D. Savrina30,31, P. Schaack52, M.
Schiller41, H. Schindler37, M. Schlupp9, M. Schmelling10, B. Schmidt37, O.
Schneider38, A. Schopper37, M.-H. Schune7, R. Schwemmer37, B. Sciascia18, A.
Sciubba24, M. Seco36, A. Semennikov30, K. Senderowska26, I. Sepp52, N.
Serra39, J. Serrano6, P. Seyfert11, M. Shapkin34, I. Shapoval16,42, P.
Shatalov30, Y. Shcheglov29, T. Shears51,37, L. Shekhtman33, O. Shevchenko42,
V. Shevchenko30, A. Shires52, R. Silva Coutinho47, T. Skwarnicki58, N.A.
Smith51, E. Smith54,48, M. Smith53, M.D. Sokoloff56, F.J.P. Soler50, F.
Soomro18, D. Souza45, B. Souza De Paula2, B. Spaan9, A. Sparkes49, P.
Spradlin50, F. Stagni37, S. Stahl11, O. Steinkamp39, S. Stoica28, S. Stone58,
B. Storaci39, M. Straticiuc28, U. Straumann39, V.K. Subbiah37, L. Sun56, S.
Swientek9, V. Syropoulos41, M. Szczekowski27, P. Szczypka38,37, T. Szumlak26,
S. T’Jampens4, M. Teklishyn7, E. Teodorescu28, F. Teubert37, C. Thomas54, E.
Thomas37, J. van Tilburg11, V. Tisserand4, M. Tobin38, S. Tolk41, D.
Tonelli37, S. Topp-Joergensen54, N. Torr54, E. Tournefier4,52, S. Tourneur38,
M.T. Tran38, M. Tresch39, A. Tsaregorodtsev6, P. Tsopelas40, N. Tuning40, M.
Ubeda Garcia37, A. Ukleja27, D. Urner53, U. Uwer11, V. Vagnoni14, G.
Valenti14, R. Vazquez Gomez35, P. Vazquez Regueiro36, S. Vecchi16, J.J.
Velthuis45, M. Veltri17,g, G. Veneziano38, M. Vesterinen37, B. Viaud7, D.
Vieira2, X. Vilasis-Cardona35,n, A. Vollhardt39, D. Volyanskyy10, D. Voong45,
A. Vorobyev29, V. Vorobyev33, C. Voß60, H. Voss10, R. Waldi60, R. Wallace12,
S. Wandernoth11, J. Wang58, D.R. Ward46, N.K. Watson44, A.D. Webber53, D.
Websdale52, M. Whitehead47, J. Wicht37, J. Wiechczynski25, D. Wiedner11, L.
Wiggers40, G. Wilkinson54, M.P. Williams47,48, M. Williams55, F.F. Wilson48,
J. Wishahi9, M. Witek25, S.A. Wotton46, S. Wright46, S. Wu3, K. Wyllie37, Y.
Xie49,37, F. Xing54, Z. Xing58, Z. Yang3, R. Young49, X. Yuan3, O.
Yushchenko34, M. Zangoli14, M. Zavertyaev10,a, F. Zhang3, L. Zhang58, W.C.
Zhang12, Y. Zhang3, A. Zhelezov11, A. Zhokhov30, L. Zhong3, A. Zvyagin37.
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é de Savoie, CNRS/IN2P3, Annecy-Le-Vieux, France
5Clermont Université, Université Blaise Pascal, CNRS/IN2P3, LPC, Clermont-
Ferrand, France
6CPPM, Aix-Marseille Université, CNRS/IN2P3, Marseille, France
7LAL, Université Paris-Sud, CNRS/IN2P3, Orsay, France
8LPNHE, Université Pierre et Marie Curie, Université Paris Diderot,
CNRS/IN2P3, Paris, France
9Fakultät Physik, Technische Universität Dortmund, Dortmund, Germany
10Max-Planck-Institut für Kernphysik (MPIK), Heidelberg, Germany
11Physikalisches Institut, Ruprecht-Karls-Universität Heidelberg, Heidelberg,
Germany
12School of Physics, University College Dublin, Dublin, Ireland
13Sezione INFN di Bari, Bari, Italy
14Sezione INFN di Bologna, Bologna, Italy
15Sezione INFN di Cagliari, Cagliari, Italy
16Sezione INFN di Ferrara, Ferrara, Italy
17Sezione INFN di Firenze, Firenze, Italy
18Laboratori Nazionali dell’INFN di Frascati, Frascati, Italy
19Sezione INFN di Genova, Genova, Italy
20Sezione INFN di Milano Bicocca, Milano, Italy
21Sezione INFN di Padova, Padova, Italy
22Sezione INFN di Pisa, Pisa, Italy
23Sezione INFN di Roma Tor Vergata, Roma, Italy
24Sezione INFN di Roma La Sapienza, Roma, Italy
25Henryk Niewodniczanski Institute of Nuclear Physics Polish Academy of
Sciences, Kraków, Poland
26AGH - University of Science and Technology, Faculty of Physics and Applied
Computer Science, Kraków, Poland
27National Center for Nuclear Research (NCBJ), Warsaw, Poland
28Horia Hulubei National Institute of Physics and Nuclear Engineering,
Bucharest-Magurele, Romania
29Petersburg Nuclear Physics Institute (PNPI), Gatchina, Russia
30Institute of Theoretical and Experimental Physics (ITEP), Moscow, Russia
31Institute of Nuclear Physics, Moscow State University (SINP MSU), Moscow,
Russia
32Institute for Nuclear Research of the Russian Academy of Sciences (INR RAN),
Moscow, Russia
33Budker Institute of Nuclear Physics (SB RAS) and Novosibirsk State
University, Novosibirsk, Russia
34Institute for High Energy Physics (IHEP), Protvino, Russia
35Universitat de Barcelona, Barcelona, Spain
36Universidad de Santiago de Compostela, Santiago de Compostela, Spain
37European Organization for Nuclear Research (CERN), Geneva, Switzerland
38Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
39Physik-Institut, Universität Zürich, Zürich, Switzerland
40Nikhef National Institute for Subatomic Physics, Amsterdam, The Netherlands
41Nikhef National Institute for Subatomic Physics and VU University Amsterdam,
Amsterdam, The Netherlands
42NSC Kharkiv Institute of Physics and Technology (NSC KIPT), Kharkiv, Ukraine
43Institute for Nuclear Research of the National Academy of Sciences (KINR),
Kyiv, Ukraine
44University of Birmingham, Birmingham, United Kingdom
45H.H. Wills Physics Laboratory, University of Bristol, Bristol, United
Kingdom
46Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom
47Department of Physics, University of Warwick, Coventry, United Kingdom
48STFC Rutherford Appleton Laboratory, Didcot, United Kingdom
49School of Physics and Astronomy, University of Edinburgh, Edinburgh, United
Kingdom
50School of Physics and Astronomy, University of Glasgow, Glasgow, United
Kingdom
51Oliver Lodge Laboratory, University of Liverpool, Liverpool, United Kingdom
52Imperial College London, London, United Kingdom
53School of Physics and Astronomy, University of Manchester, Manchester,
United Kingdom
54Department of Physics, University of Oxford, Oxford, United Kingdom
55Massachusetts Institute of Technology, Cambridge, MA, United States
56University of Cincinnati, Cincinnati, OH, United States
57University of Maryland, College Park, MD, United States
58Syracuse University, Syracuse, NY, United States
59Pontifícia Universidade Católica do Rio de Janeiro (PUC-Rio), Rio de
Janeiro, Brazil, associated to 2
60Institut für Physik, Universität Rostock, Rostock, Germany, associated to 11

aP.N. Lebedev Physical Institute, Russian Academy of Science (LPI RAS),
Moscow, Russia
bUniversità di Bari, Bari, Italy
cUniversità di Bologna, Bologna, Italy
dUniversità di Cagliari, Cagliari, Italy
eUniversità di Ferrara, Ferrara, Italy
fUniversità di Firenze, Firenze, Italy
gUniversità di Urbino, Urbino, Italy
hUniversità di Modena e Reggio Emilia, Modena, Italy
iUniversità di Genova, Genova, Italy
jUniversità di Milano Bicocca, Milano, Italy
kUniversità di Roma Tor Vergata, Roma, Italy
lUniversità di Roma La Sapienza, Roma, Italy
mUniversità della Basilicata, Potenza, Italy
nLIFAELS, La Salle, Universitat Ramon Llull, Barcelona, Spain
oIFIC, Universitat de Valencia-CSIC, Valencia, Spain
pHanoi University of Science, Hanoi, Viet Nam
qUniversità di Padova, Padova, Italy
rUniversità di Pisa, Pisa, Italy
sScuola Normale Superiore, Pisa, Italy

###### Abstract

Using $pp$ collision data, corresponding to an integrated luminosity of
$1.0~{}\mathrm{fb}^{-1}$, collected by LHCb in 2011 at a center-of-mass energy
of $7\\!$ $\mathrm{\,Te\kern-1.00006ptV}$, we report the measurement of direct
$C\\!P$ violation in $B^{0}_{s}\rightarrow K^{-}\pi^{+}$ decays,
$A_{C\\!P}(B^{0}_{s}\rightarrow K^{-}\pi^{+})=0.27\pm 0.04\,\mathrm{(stat)}\pm
0.01\,\mathrm{(syst)}$, with significance exceeding five standard deviations.
This is the first observation of $C\\!P$ violation in the decays of
$B^{0}_{s}$ mesons. Furthermore, we provide an improved determination of
direct $C\\!P$ violation in $B^{0}\rightarrow K^{+}\pi^{-}$ decays,
$A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})=-0.080\pm 0.007\,\mathrm{(stat)}\pm
0.003\,\mathrm{(syst)}$, which is the most precise measurement of this
quantity to date.

###### pacs:

Valid PACS appear here

EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

 |  |
---|---|---
|  | CERN-PH-EP-2013-068
|  | LHCb-PAPER-2013-018
|  | June 7, 2013

© CERN on behalf of the LHCb collaboration, license CC-BY-3.0.

Submitted to Phys. Rev. Lett.

The non-invariance of fundamental interactions under the combined action of
the charge conjugation ($C$) and parity ($P$) transformations is
experimentally well established in the $K^{0}$ and $B^{0}$ meson systems
Christenson:1964fg ; Aubert:2001nu ; Abe:2001xe ; PDG2012 . The Standard Model
(SM) description of $C\\!P$ violation, as given by the Cabibbo-Kobayashi-
Maskawa (CKM) theory of quark-flavor mixing Cabibbo:1963yz ; Kobayashi:1973fv
, has been very successful in describing existing data. However, the source of
$C\\!P$ violation in the SM is known to be too small to account for the
matter-dominated universe Cohen:1993nk ; Riotto:1999yt ; Hou:2008xd .

The study of $C\\!P$ violation in charmless charged two-body decays of neutral
$B$ mesons provides stringent tests of the CKM picture in the SM, and is a
sensitive probe to search for the presence of non-SM physics Deshpande:1994ii
; He:1998rq ; Fleischer:1999pa ; Gronau:2000md ; Lipkin:2005pb ;
Fleischer:2007hj ; Fleischer:2010ib . However, quantitative SM predictions for
$C\\!P$ violation in these decays are challenging because of the presence of
hadronic factors in the decay amplitudes, which cannot be accurately
calculated from quantum chromodynamics (QCD) at present. It is crucial to
combine several measurements from such two-body decays, exploiting approximate
flavor symmetries in order to cancel the unknown parameters. An experimental
program for measuring the properties of these decays has been carried out
during the last decade at the $B$ factories Lees:2013bb ; PhysRevD.87.031103
and at the Tevatron Aaltonen:2011qt , and is now continued by LHCb with
increased sensitivity. The discovery of direct $C\\!P$ violation in the
$B^{0}\rightarrow K^{+}\pi^{-}$ decay dates back to 2004 Aubert:2004qm ;
Chao:2004jy . This observation raised the question of whether the effect could
be accommodated by the SM or was due to non-SM physics. A simple but powerful
model-independent test was proposed in Refs. He:1998rq ; Lipkin:2005pb , which
required the measurement of direct $C\\!P$ violation in the
$B^{0}_{s}\rightarrow K^{-}\pi^{+}$ decay. However, $C\\!P$ violation has
never been observed with significance exceeding five Gaussian standard
deviations ($\sigma$) in any $B^{0}_{s}$ meson decay so far.

In this Letter we report measurements of direct $C\\!P$-violating asymmetries
in $B^{0}\rightarrow K^{+}\pi^{-}$ and $B^{0}_{s}\rightarrow K^{-}\pi^{+}$
decays using $pp$ collision data, corresponding to an integrated luminosity of
$1.0~{}\mathrm{fb}^{-1}$, collected with the LHCb detector in 2011 at a
center-of-mass energy of $7\\!$ $\mathrm{\,Te\kern-1.00006ptV}$. The present
results supersede those given in Ref. LHCb-PAPER-2011-029 . The inclusion of
charge-conjugate decay modes is implied except in the asymmetry definitions.
The direct $C\\!P$ asymmetry in the $B^{0}_{(s)}$ decay rate to the final
state $f_{(s)}$, with $f=K^{+}\pi^{-}$ and $f_{s}=K^{-}\pi^{+}$, is defined as

$A_{C\\!P}\\!\\!\left(\\!B^{0}_{(s)}\\!\\!\rightarrow\\!\\!f_{(s)}\\!\right)\\!\\!=\\!\Phi\\!\left[\Gamma\\!\left(\\!\overline{B}^{0}_{(s)}\\!\\!\rightarrow\\!\\!\bar{f}_{(s)}\\!\right)\\!\\!,\,\Gamma\\!\left(\\!B^{0}_{(s)}\\!\\!\rightarrow\\!\\!f_{(s)}\\!\right)\right]\\!\\!,$
(1)

where $\Phi[X,\,Y]=(X-Y)/(X+Y)$ and $\bar{f}_{(s)}$ denotes the charge-
conjugate of $f_{(s)}$.

The LHCb detector Alves:2008zz  is a single-arm forward spectrometer covering
the pseudorapidity range $2<\eta<5$, designed for the study of particles
containing $b$ or $c$ quarks. The trigger LHCb-DP-2012-004  consists of a
hardware stage, based on information from the calorimeter and muon systems,
followed by a software stage that applies a full event reconstruction. The
hadronic hardware trigger selects large transverse energy clusters in the
hadronic calorimeter. The software trigger requires a two-, three-, or four-
track secondary vertex with a large sum of the transverse momenta
($p_{\mathrm{T}}$) of the tracks and a significant displacement from the
primary $pp$ interaction vertices (PVs). At least one track should have
$p_{\mathrm{T}}$ and impact parameter (IP) $\chi^{2}$ with respect to all PVs
exceeding given thresholds. The IP is defined as the distance between the
reconstructed trajectory of a particle and a given $pp$ collision vertex, and
the IP $\chi^{2}$ is the difference between the $\chi^{2}$ of the PV
reconstructed with and without the considered track. A multivariate algorithm
is used for the identification of secondary vertices consistent with the decay
of a $b$ hadron. In order to improve the trigger efficiency on hadronic two-
body decays, a dedicated two-body software trigger is also used. This trigger
imposes requirements on the following quantities: the quality of the online-
reconstructed tracks, their $p_{\mathrm{T}}$ and IP; the distance of closest
approach of the decay products of the $B$ meson candidate, its
$p_{\mathrm{T}}$, IP and the decay time in its rest frame.

Figure 1: Invariant mass spectra obtained using the event selection adopted
for the best sensitivity on (a, b) $A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})$
and (c, d) $A_{C\\!P}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$. Panels (a) and (c)
represent the $K^{+}\pi^{-}$ invariant mass, whereas panels (b) and (d)
represent the $K^{-}\pi^{+}$ invariant mass. The results of the unbinned
maximum likelihood fits are overlaid. The main components contributing to the
fit model are also shown.

More selective requirements are applied offline. Two sets of criteria have
been optimized with the aim of minimizing the expected statistical uncertainty
either on $A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})$ or on
$A_{C\\!P}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$. In addition to the
requirements on the kinematic variables already used in the trigger,
requirements on the largest $p_{\mathrm{T}}$ and IP of the $B$ daughter
particles are applied. In the case of $B^{0}_{s}\rightarrow K^{-}\pi^{+}$
decays, a tighter selection is needed to achieve stronger rejection of
combinatorial background. For example, the decay time is required to exceed
$1.5$ ps, whereas in the $B^{0}\rightarrow K^{+}\pi^{-}$ selection a lower
threshold of $0.9$ ps is applied. This is because the probability for a $b$
quark to form a $B^{0}_{s}$ meson, which subsequently decays to the
$K^{-}\pi^{+}$ final state, is one order of magnitude smaller than that to
form a $B^{0}$ meson decaying to $K^{+}\pi^{-}$ LHCb-PAPER-2012-002 . The two
samples are then subdivided according to the various final states using the
particle identification (PID) provided by the two ring-imaging Cherenkov
(RICH) detectors LHCb-DP-2012-003 . Two sets of PID selection criteria are
applied: a loose set optimized for the measurement of
$A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})$ and a tight set for that of
$A_{C\\!P}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$. More details on the event
selection can be found in Ref. LHCb-PAPER-2011-029 .

To determine the amount of background events from other two-body $b$-hadron
decays with a misidentified pion or kaon (cross-feed background), the relative
efficiencies of the RICH PID selection criteria must be determined. This is
achieved by means of a data-driven method that uses $D^{*+}\rightarrow
D^{0}(K^{-}\pi^{+})\pi^{+}$ and $\mathchar 28931\relax\rightarrow p\pi^{-}$
decays as control samples. The production and decay kinematic properties of
the $D^{0}\rightarrow K^{-}\pi^{+}$ and $\mathchar 28931\relax\rightarrow
p\pi^{-}$ channels differ from those of the $b$-hadron decays under study.
Since the RICH PID information is momentum dependent, a calibration procedure
is performed by reweighting the distributions of the PID variables obtained
from the calibration samples, in order to match the momentum distributions of
signal final-state particles observed in data.

Unbinned maximum likelihood fits to the mass spectra of the selected events
are performed. The $B^{0}\rightarrow K^{+}\pi^{-}$ and $B^{0}_{s}\rightarrow
K^{-}\pi^{+}$ signal components are described by double Gaussian functions
convolved with a function that describes the effect of final-state radiation
Baracchini:2005wp . The background due to partially reconstructed three-body
$B$ decays is parameterized by means of two ARGUS functions Albrecht:1989ga
convolved with a Gaussian resolution function. The combinatorial background is
modeled by an exponential function and the shapes of the cross-feed
backgrounds, mainly due to $B^{0}\rightarrow\pi^{+}\pi^{-}$ and
$B^{0}_{s}\rightarrow K^{+}K^{-}$ decays with one misidentified particle in
the final state, are obtained from simulation. The cross-feed background
yields are determined from the $\pi^{+}\pi^{-}$, $K^{+}K^{-}$, $p\pi^{-}$ and
$pK^{-}$ mass spectra, using events passing the same selection as the signal
and taking into account the appropriate PID efficiency factors. The
$K^{+}\pi^{-}$ and $K^{-}\pi^{+}$ mass spectra for the events passing the two
selections are shown in Fig. 1. The average invariant mass resolution is about
$22~{}{\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$.

From the two mass fits we determine the signal yields $N(B^{0}\rightarrow
K^{+}\pi^{-})=41\hskip 1.42262pt420\pm 300$ and $N(B^{0}_{s}\rightarrow
K^{-}\pi^{+})=1065\pm 55$, as well as the raw asymmetries
$A_{\mathrm{raw}}(B^{0}\rightarrow K^{+}\pi^{-})=-0.091\pm 0.006$ and
$A_{\mathrm{raw}}(B^{0}_{s}\rightarrow K^{-}\pi^{+})=0.28\pm 0.04$, where the
uncertainties are statistical only. In order to derive the $C\\!P$ asymmetries
from the observed raw asymmetries, effects induced by the detector acceptance
and event reconstruction, as well as due to interactions of final-state
particles with the detector material, must be accounted for. Furthermore, the
possible presence of a $B^{0}_{(s)}-\overline{B}^{0}_{(s)}$ production
asymmetry must also be considered.

The $C\\!P$ asymmetry is related to the raw asymmetry by
$A_{C\\!P}=A_{\mathrm{raw}}-A_{\Delta}$, where the correction $A_{\Delta}$ is
defined as

$A_{\Delta}(B^{0}_{(s)}\rightarrow
K\pi)=\zeta_{d(s)}A_{\mathrm{D}}(K\pi)+\kappa_{d(s)}A_{\mathrm{P}}(B^{0}_{(s)}),$
(2)

with $\zeta_{d}=1$ and $\zeta_{s}=-1$. The instrumental asymmetry
$A_{\mathrm{D}}(K\pi)$ is given in terms of the detection efficiencies
$\varepsilon_{\mathrm{D}}$ of the charge-conjugate final states by
$A_{\mathrm{D}}(K\pi)=\Phi[\varepsilon_{\mathrm{D}}(K^{-}\pi^{+}),\,\varepsilon_{\mathrm{D}}(K^{+}\pi^{-})]$,
and the production asymmetry $A_{\mathrm{P}}(B^{0}_{(s)})$ is defined in terms
of the $\overline{B}^{0}_{(s)}$ and $B^{0}_{(s)}$ production rates,
$R(\overline{B}^{0}_{(s)})$ and $R(B^{0}_{(s)})$, as
$A_{\mathrm{P}}(B^{0}_{(s)})=\Phi[R(\overline{B}^{0}_{(s)}),\,R(B^{0}_{(s)})]$.
The factors $\kappa_{d}$ and $\kappa_{s}$ take into account dilutions due to
$B^{0}$ and $B^{0}_{s}$ meson mixing, respectively. Their values also depend
on event reconstruction and selection, and are $\kappa_{d}=0.303\pm 0.005$ and
$\kappa_{s}=-0.033\pm 0.003$ LHCb-PAPER-2011-029 . The factor $\kappa_{s}$ is
ten times smaller than $\kappa_{d}$, owing to the large $B^{0}_{s}$
oscillation frequency.

The instrumental charge asymmetry $A_{\mathrm{D}}(K\pi)$ is measured from data
using $D^{*+}\rightarrow D^{0}(K^{-}\pi^{+})\pi^{+}$ and $D^{*+}\rightarrow
D^{0}(K^{-}K^{+})\pi^{+}$ decays. The combination of the time-integrated raw
asymmetries of these two decay modes is used to disentangle the various
contributions to each raw asymmetry. The presence of open charm production
asymmetries arising from the primary $pp$ interaction constitutes an
additional complication. We write the following equations relating the
observed raw asymmetries to the physical $C\\!P$ asymmetries

$\displaystyle A^{*}_{\mathrm{raw}}(K\pi)$ $\displaystyle=$ $\displaystyle
A^{*}_{\mathrm{D}}(\pi_{s})+A^{*}_{\mathrm{D}}(K\pi)+A_{\mathrm{P}}(D^{*}),$
(3) $\displaystyle A^{*}_{\mathrm{raw}}(KK)$ $\displaystyle=$ $\displaystyle
A_{C\\!P}(KK)+A^{*}_{\mathrm{D}}(\pi_{s})+A_{\mathrm{P}}(D^{*}),$ (4)

where $A^{*}_{\mathrm{raw}}(K\pi)$ and $A^{*}_{\mathrm{raw}}(KK)$ are the
time-integrated raw asymmetries in $D^{*}$-tagged $D^{0}\rightarrow
K^{-}\pi^{+}$ and $D^{0}\rightarrow K^{-}K^{+}$ decays, respectively;
$A_{C\\!P}(KK)$ is the $D^{0}\rightarrow K^{-}K^{+}$ $C\\!P$ asymmetry;
$A^{*}_{\mathrm{D}}(K\pi)$ is the detection asymmetry in reconstructing
$D^{0}\rightarrow K^{-}\pi^{+}$ and $\overline{D}^{0}\rightarrow K^{+}\pi^{-}$
decays; $A^{*}_{\mathrm{D}}(\pi_{s})$ is the detection asymmetry in
reconstructing positively- and negatively-charged pions originated from
$D^{*}$ decays; and $A_{\mathrm{P}}(D^{*})$ is the production asymmetry for
prompt charged $D^{*}$ mesons. In Eq. (3) any possible $C\\!P$ asymmetry in
the Cabibbo-favored $D^{0}\rightarrow K^{-}\pi^{+}$ decay is neglected
Bianco:2003vb . By subtracting Eqs. (3) and (4), one obtains

$A^{*}_{\mathrm{raw}}(K\pi)-A^{*}_{\mathrm{raw}}(KK)=A^{*}_{\mathrm{D}}(K\pi)-A_{C\\!P}(KK).$
(5)

Once the raw asymmetries are measured, this equation determines unambiguously
the detection asymmetry $A^{*}_{\mathrm{D}}(K\pi)$, using the world average
for the $C\\!P$ asymmetry of the $D^{0}\rightarrow K^{-}K^{+}$ decay. Since
the measured value of the time-integrated asymmetry depends on the decay-time
acceptance, the existing measurements of $A_{C\\!P}(KK)$ Staric:2008rx ;
Aubert:2007if ; Aaltonen:2011se  are corrected for the difference in
acceptance with respect to LHCb LHCb-PAPER-2011-023 . This leads to the value
$A_{C\\!P}(KK)=(-0.24\pm 0.18)\%$. Furthermore, $B$ meson production and decay
kinematic properties differ from those of the $D$ decays being considered, and
different trigger and selection algorithms are applied. In order to correct
the raw asymmetries of $B$ decays, using the detection asymmetry
$A^{*}_{\mathrm{D}}(K\pi)$ derived from $D$ decays, a reweighting procedure is
needed. We reweight the $D^{0}$ momentum, transverse momentum and azimuthal
angle in $D^{0}\rightarrow K^{-}\pi^{+}$ and $D^{0}\rightarrow K^{-}K^{+}$
decays, to match the respective $B^{0}_{(s)}$ distributions in
$B^{0}\rightarrow K^{+}\pi^{-}$ and $B^{0}_{s}\rightarrow K^{-}\pi^{+}$
decays. The raw asymmetries are determined by means of $\chi^{2}$ fits to the
reweighted $\delta m=M_{D^{*}}-M_{D^{0}}$ distributions, where $M_{D^{*}}$ and
$M_{D^{0}}$ are the reconstructed $D^{*}$ and $D^{0}$ candidate invariant
masses, respectively.

From the raw asymmetries, values for the quantity $\Delta
A=A_{\mathrm{D}}(K\pi)-A_{C\\!P}(KK)$ are determined. We obtain the values
$\Delta A=(-0.91\pm 0.15)\%$ and $\Delta A=(-0.98\pm 0.11)\%$, using as target
kinematic distributions those of $B$ candidates passing the event selection
optimized for $A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})$ and for
$A_{C\\!P}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$, respectively. Using these two
values of $\Delta A$ and the value of $A_{C\\!P}(KK)$, we obtain the
instrumental asymmetries $A_{\mathrm{D}}(K\pi)=(-1.15\pm 0.23)\%$ for the
$B^{0}\rightarrow K^{+}\pi^{-}$ decay and $A_{\mathrm{D}}(K\pi)=(-1.22\pm
0.21)\%$ for the $B^{0}_{s}\rightarrow K^{-}\pi^{+}$ decay.

Assuming negligible $C\\!P$ violation in the mixing, as expected in the SM and
confirmed by current experimental determinations bib:hfagbase , the decay rate
of a $B^{0}_{(s)}$ meson with production asymmetry $A_{\mathrm{P}}$, decaying
into a flavor-specific final state $f_{(s)}$ with $C\\!P$ asymmetry
$A_{C\\!P}$ and detection asymmetry $A_{\mathrm{D}}$, can be written as

$\mathcal{R}(t;\,p)\propto\left(1\\!-\\!pA_{C\\!P}\right)\left(1\\!-\\!pA_{\mathrm{D}}\right)\left[H_{+}\\!\left(t\right)\\!-\\!pA_{\mathrm{P}}H_{-}\\!\left(t\right)\right]\\!,$
(6)

where $t$ is the reconstructed decay time of the $B$ meson and $p$ assumes the
values $p=+1$ for the final state $f_{(s)}$ and $p=-1$ for the final state
$\bar{f}_{(s)}$. The functions $H_{+}\left(t\right)$ and $H_{-}\left(t\right)$
are defined as

$\displaystyle H_{+}\\!\left(t\right)\\!$ $\displaystyle=$
$\displaystyle\\!\\!\left[e^{-\Gamma_{d(s)}t^{\prime}}\\!\\!\\!\cosh\\!{\left(\\!\\!\frac{\Delta\Gamma_{d(s)}}{2}t^{\prime}\\!\\!\right)}\\!\\!\otimes\\!R\\!\left(t,\,t^{\prime}\right)\right]\\!\varepsilon_{d(s)}\\!\left(t\right)\\!,$
(7) $\displaystyle H_{-}\\!\left(t\right)\\!$ $\displaystyle=$
$\displaystyle\\!\\!\left[e^{-\Gamma_{d(s)}t^{\prime}}\\!\\!\\!\cos\\!{\left(\Delta
m_{d(s)}t^{\prime}\right)}\\!\otimes\\!R\\!\left(t,\,t^{\prime}\right)\right]\\!\varepsilon_{d(s)}\\!\left(t\right)\\!,$
(8)

where $\Gamma_{d(s)}$ is the average decay width of the $B^{0}_{(s)}$ meson,
$\Delta\Gamma_{d(s)}$ and $\Delta m_{d(s)}$ are the decay width and mass
differences between the two $B^{0}_{(s)}$ mass eigenstates respectively,
$R\left(t,\,t^{\prime}\right)$ is the decay time resolution ($\sigma\simeq
50~{}\mathrm{fs}$ in our case) and the symbol $\otimes$ stands for
convolution. Finally $\varepsilon_{d(s)}\left(t\right)$ is the acceptance as a
function of the $B^{0}_{(s)}$ decay time. Using Eq. (6) we obtain the
following expression for the time-dependent asymmetry

$\displaystyle\mathcal{A}\left(t\right)$ $\displaystyle=$
$\displaystyle\Phi\\!\left[\mathcal{R}\left(t;\,-1\right)\\!,\,\mathcal{R}\left(t;\,+1\right)\right]$
$\displaystyle=$
$\displaystyle\frac{\left(A_{C\\!P}\\!+\\!A_{\mathrm{D}}\right)\\!H_{+}\\!\left(t\right)\\!+\\!A_{\mathrm{P}}\left(1\\!+\\!A_{C\\!P}A_{\mathrm{D}}\right)\\!H_{-}\\!\left(t\right)}{\left(1\\!+\\!A_{C\\!P}A_{\mathrm{D}}\right)\\!H_{+}\\!\left(t\right)\\!+\\!A_{\mathrm{P}}\left(A_{C\\!P}\\!+\\!A_{\mathrm{D}}\right)\\!H_{-}\\!\left(t\right)}.$

For illustrative purposes only, we consider the case of perfect decay time
resolution and negligible $\Delta\Gamma$, retaining only first-order terms in
$A_{C\\!P}$, $A_{\mathrm{P}}$ and $A_{\mathrm{D}}$. In this case, Eq. (First
observation of $C\\!P$ violation in the decays of $B^{0}_{s}$ mesons) reduces
to the expression

$\mathcal{A}\left(t\right)\approx
A_{C\\!P}+A_{\mathrm{D}}+A_{\mathrm{P}}\cos{\left(\Delta m_{d(s)}t\right)},$
(10)

i.e., the time-dependent asymmetry has an oscillatory term with amplitude
equal to the production asymmetry $A_{\mathrm{P}}$. By studying the full time-
dependent decay rate it is then possible to determine $A_{\mathrm{P}}$
unambiguously.

In order to measure the production asymmetry $A_{\mathrm{P}}$ for $B^{0}$ and
$B^{0}_{s}$ mesons, we perform fits to the decay time spectra of the $B$
candidates, separately for the events passing the two selections. The $B^{0}$
production asymmetry is determined from the sample obtained applying the
selection optimized for the measurement of $A_{C\\!P}(B^{0}\rightarrow
K^{+}\pi^{-})$, whereas the $B^{0}_{s}$ production asymmetry is determined
from the sample obtained applying the selection optimized for the measurement
of $A_{C\\!P}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$. We obtain
$A_{\mathrm{P}}(B^{0})=(0.1\pm 1.0)\%$ and $A_{\mathrm{P}}(B^{0}_{s})=(4\pm
8)\%$. Figure 2 shows the raw asymmetries as a function of the decay time,
obtained by performing fits to the invariant mass distributions of events
restricted to independent intervals of the $B$ candidate decay times.

Figure 2: Raw asymmetries as a function of the decay time for (a)
$B^{0}\rightarrow K^{+}\pi^{-}$ and (b) $B^{0}_{s}\rightarrow K^{-}\pi^{+}$
decays. In (b), the offset $t_{0}=1.5$${\rm\,ps}$ corresponds to the minimum
value of the decay time required by the $B^{0}_{s}\rightarrow K^{-}\pi^{+}$
event selection. The curves represent the asymmetry projections of fits to the
decay time spectra.

By using the values of the detection and production asymmetries, the
correction factors to the raw asymmetries $A_{\Delta}(B^{0}\rightarrow
K^{+}\pi^{-})=\left(-1.12\pm 0.23\pm 0.30\right)\\!\%$ and
$A_{\Delta}(B^{0}_{s}\rightarrow K^{-}\pi^{+})=\left(1.09\pm 0.21\pm
0.26\right)\\!\%$ are obtained, where the first uncertainties are due to the
detection asymmetry and the second to the production asymmetry.

Systematic uncertainties on the asymmetries are related to PID calibration,
modeling of the signal and background components in the maximum likelihood
fits and instrumental charge asymmetries. In order to estimate the impact of
imperfect PID calibration, we perform mass fits to determine raw asymmetries
using altered numbers of cross-feed background events, according to the
systematic uncertainties affecting the PID efficiencies. An estimate of the
uncertainty due to possible mismodeling of the final-state radiation is
determined by varying the amount of emitted radiation Baracchini:2005wp  in
the signal shape parameterization, according to studies performed on fully
simulated events, in which final state radiation is generated using Photos
Golonka:2005pn . The possibility of an incorrect description of the signal
mass model is investigated by replacing the double Gaussian function with the
sum of three Gaussian functions, where the third component has fixed fraction
($5\%$) and width ($50\,{\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$), and is
aimed at describing long tails, as observed in simulation. To assess a
systematic uncertainty on the shape of the partially reconstructed
backgrounds, we remove the second ARGUS function. For the modeling of the
combinatorial background component, the fit is repeated using a straight line.
Finally, for the case of the cross-feed backgrounds, two distinct systematic
uncertainties are estimated: one due to a relative bias in the mass scale of
the simulated distributions with respect to the signal distributions in data,
and another accounting for the difference in mass resolution between
simulation and data. All shifts from the relevant baseline values are
accounted for as systematic uncertainties. Systematic uncertainties related to
the determination of detection asymmetries are calculated by summing in
quadrature the respective uncertainties on $A_{\Delta}(B^{0}\rightarrow
K^{+}\pi^{-})$ and $A_{\Delta}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$ with an
additional uncertainty of $0.10\%$, accounting for residual differences in the
trigger composition between signal and calibration samples.

The systematic uncertainties for $A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})$
and $A_{C\\!P}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$ are summarized in Table 1.
Since the production asymmetries are obtained from the fitted decay time
spectra of $B^{0}\rightarrow K^{+}\pi^{-}$ and $B^{0}_{s}\rightarrow
K^{-}\pi^{+}$ decays, their uncertainties are statistical in nature and are
then propagated to the statistical uncertainties on
$A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})$ and $A_{C\\!P}(B^{0}_{s}\rightarrow
K^{-}\pi^{+})$.

Table 1: Systematic uncertainties on $A_{C\\!P}(B^{0}\rightarrow
K^{+}\pi^{-})$ and $A_{C\\!P}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$. The total
systematic uncertainties are obtained by summing the individual contributions
in quadrature.

Systematic uncertainty |  $A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})$ | $A_{C\\!P}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$
---|---|---
PID calibration | 0.0006 | 0.0012
Final state radiation | 0.0008 | 0.0020
Signal model | 0.0001 | 0.0064
Combinatorial background | 0.0004 | 0.0042
Three-body background | 0.0005 | 0.0027
Cross-feed background | 0.0010 | 0.0033
Detection asymmetry | 0.0025 | 0.0023
Total | 0.0029 | 0.0094

In conclusion, the parameters of $C\\!P$ violation in $B^{0}\rightarrow
K^{+}\pi^{-}$ and $B^{0}_{s}\rightarrow K^{-}\pi^{+}$ decays have been
measured to be

$\displaystyle A_{C\\!P}(B^{0}\\!\rightarrow\\!K^{+}\pi^{-})$ $\displaystyle=$
$\displaystyle-0.080\pm 0.007\,\mathrm{(stat)}\pm 0.003\,\mathrm{(syst)},$
$\displaystyle A_{C\\!P}(B^{0}_{s}\\!\rightarrow\\!K^{-}\pi^{+})$
$\displaystyle=$ $\displaystyle 0.27\pm 0.04\,\mathrm{(stat)}\pm
0.01\,\mathrm{(syst)}.$

Dividing the central values by the sum in quadrature of statistical and
systematic uncertainties, the significances of the measured deviations from
zero are $10.5\sigma$ and $6.5\sigma$, respectively. The former is the most
precise measurement of $A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})$ to date,
whereas the latter represents the first observation of $C\\!P$ violation in
decays of $B^{0}_{s}$ mesons with significance exceeding $5\sigma$. Both
measurements are in good agreement with world averages bib:hfagbase  and
previous LHCb results LHCb-PAPER-2011-029 .

These results allow a stringent test of the validity of the relation between
$A_{C\\!P}(B^{0}\rightarrow K^{+}\pi^{-})$ and $A_{C\\!P}(B^{0}_{s}\rightarrow
K^{-}\pi^{+})$ in the SM given in Ref. Lipkin:2005pb  as

$\Delta=\frac{A_{C\\!P}(B^{0}\\!\rightarrow\\!K^{+}\pi^{-})}{A_{C\\!P}(B^{0}_{s}\\!\rightarrow\\!K^{-}\pi^{+})}+\frac{\mathcal{B}(B^{0}_{s}\\!\rightarrow\\!K^{-}\pi^{+})}{\mathcal{B}(B^{0}\\!\rightarrow\\!K^{+}\pi^{-})}\frac{\tau_{d}}{\tau_{s}}=0,$
(11)

where $\mathcal{B}(B^{0}\rightarrow K^{+}\pi^{-})$ and
$\mathcal{B}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$ are $C\\!P$-averaged
branching fractions, and $\tau_{d}$ and $\tau_{s}$ are the $B^{0}$ and
$B^{0}_{s}$ mean lifetimes, respectively. Using additional results for
$\mathcal{B}(B^{0}\rightarrow K^{+}\pi^{-})$ and
$\mathcal{B}(B^{0}_{s}\rightarrow K^{-}\pi^{+})$ LHCb-PAPER-2012-002  and the
world averages for $\tau_{d}$ and $\tau_{s}$ bib:hfagbase , we obtain
$\Delta=-0.02\pm 0.05\pm 0.04$, where the first uncertainty is from the
measurements of the $C\\!P$ asymmetries and the second is from the input
values of the branching fractions and the lifetimes. No evidence for a
deviation from zero of $\Delta$ is observed with the present experimental
precision.

## Acknowledgements

We express our gratitude to our colleagues in the CERN accelerator departments
for the excellent performance of the LHC. We thank the technical and
administrative staff at the LHCb institutes. We acknowledge support from CERN
and from the national agencies: CAPES, CNPq, FAPERJ and FINEP (Brazil); NSFC
(China); CNRS/IN2P3 and Region Auvergne (France); BMBF, DFG, HGF and MPG
(Germany); SFI (Ireland); INFN (Italy); FOM and NWO (The Netherlands); SCSR
(Poland); ANCS/IFA (Romania); MinES, Rosatom, RFBR and NRC “Kurchatov
Institute” (Russia); MinECo, XuntaGal and GENCAT (Spain); SNSF and SER
(Switzerland); NAS Ukraine (Ukraine); STFC (United Kingdom); NSF (USA). We
also acknowledge the support received from the ERC under FP7. The Tier1
computing centres are supported by IN2P3 (France), KIT and BMBF (Germany),
INFN (Italy), NWO and SURF (The Netherlands), PIC (Spain), GridPP (United
Kingdom). We are thankful for the computing resources put at our disposal by
Yandex LLC (Russia), as well as to the communities behind the multiple open
source software packages that we depend on.

## References

  * (1) J. H. Christenson, J. W. Cronin, V. L. Fitch, and R. Turlay, Evidence for the $2\pi$ decay of the $K_{2}^{0}$ meson, Phys. Rev. Lett. 13 (1964) 138
  * (2) BaBar collaboration, B. Aubert et al., Observation of CP violation in the $B^{0}$ meson system, Phys. Rev. Lett. 87 (2001) 091801, arXiv:hep-ex/0107013
  * (3) Belle collaboration, K. Abe et al., Observation of large CP violation in the neutral B meson system, Phys. Rev. Lett. 87 (2001) 091802, arXiv:hep-ex/0107061
  * (4) Particle Data Group, J. Beringer et al., Review of particle physics, Phys. Rev. D86 (2012) 010001
  * (5) N. Cabibbo, Unitary symmetry and leptonic decays, Phys. Rev. Lett. 10 (1963) 531
  * (6) M. Kobayashi and T. Maskawa, CP violation in the renormalizable theory of weak interaction, Prog. Theor. Phys. 49 (1973) 652
  * (7) A. G. Cohen, D. Kaplan, and A. Nelson, Progress in electroweak baryogenesis, Ann. Rev. Nucl. Part. Sci. 43 (1993) 27, arXiv:hep-ph/9302210
  * (8) A. Riotto and M. Trodden, Recent progress in baryogenesis, Ann. Rev. Nucl. Part. Sci. 49 (1999) 35, arXiv:hep-ph/9901362
  * (9) W.-S. Hou, Source of CP violation for the baryon asymmetry of the universe, Chin. J. Phys. 47 (2009) 134, arXiv:0803.1234
  * (10) N. G. Deshpande and X.-G. He, CP asymmetry relations between $\bar{B}^{0}\rightarrow\pi\pi$ and $\bar{B}^{0}\rightarrow\pi K$ rates, Phys. Rev. Lett. 75 (1995) 1703, arXiv:hep-ph/9412393
  * (11) X.-G. He, SU(3) analysis of annihilation contributions and CP violating relations in $B\rightarrow PP$ decays, Eur. Phys. J. C9 (1999) 443, arXiv:hep-ph/9810397
  * (12) R. Fleischer, New strategies to extract $\beta$ and $\gamma$ from $B_{d}\rightarrow\pi^{+}\pi^{-}$ and $B_{s}\rightarrow K^{+}K^{-}$, Phys. Lett. B459 (1999) 306, arXiv:hep-ph/9903456
  * (13) M. Gronau and J. L. Rosner, The role of $B_{s}\rightarrow K\pi$ in determining the weak phase $\gamma$, Phys. Lett. B482 (2000) 71, arXiv:hep-ph/0003119
  * (14) H. J. Lipkin, Is observed direct CP violation in $B_{d}\rightarrow K^{+}\pi^{-}$ due to new physics? Check standard model prediction of equal violation in $B_{s}\rightarrow K^{-}\pi^{+}$, Phys. Lett. B621 (2005) 126, arXiv:hep-ph/0503022
  * (15) R. Fleischer, $B_{s,d}\rightarrow\pi\pi,\,\pi K,\,KK$: status and prospects, Eur. Phys. J. C52 (2007) 267, arXiv:0705.1121
  * (16) R. Fleischer and R. Knegjens, In pursuit of new physics with $B^{0}_{s}\rightarrow K^{+}K^{-}$, Eur. Phys. J. C71 (2011) 1532, arXiv:1011.1096
  * (17) BaBar collaboration, J. P. Lees et al., Measurement of CP asymmetries and branching fractions in charmless two-body B-meson decays to pions and kaons, Phys. Rev. D87 (2013) 052009, arXiv:1206.3525
  * (18) Belle collaboration, Y.-T. Duh et al., Measurements of branching fractions and direct $CP$ asymmetries for $B\rightarrow K\pi$, $B\rightarrow\pi\pi$ and $B\rightarrow KK$ decays, Phys. Rev. D87 (2013) 031103, arXiv:1210.1348
  * (19) CDF collaboration, T. Aaltonen et al., Measurements of direct CP violating asymmetries in charmless decays of strange bottom mesons and bottom baryons, Phys. Rev. Lett. 106 (2011) 181802, arXiv:1103.5762
  * (20) BaBar collaboration, B. Aubert et al., Direct CP violating asymmetry in $B^{0}\rightarrow K^{+}\pi^{-}$ decays, Phys. Rev. Lett. 93 (2004) 131801, arXiv:hep-ex/0407057
  * (21) Belle collaboration, Y. Chao et al., Improved measurements of partial rate asymmetry in $B\rightarrow hh$ decays, Phys. Rev. D71 (2005) 031502, arXiv:hep-ex/0407025
  * (22) LHCb collaboration, R. Aaij et al., First evidence of direct $C\\!P$ violation in charmless two-body decays of $B^{0}_{s}$ mesons, Phys. Rev. Lett. 108 (2012) 201601, arXiv:1202.6251
  * (23) LHCb collaboration, A. A. Alves Jr. et al., The LHCb detector at the LHC, JINST 3 (2008) S08005
  * (24) R. Aaij et al., The LHCb trigger and its performance in 2011, JINST 8 (2013) P04022, arXiv:1211.3055
  * (25) LHCb collaboration, R. Aaij et al., Measurement of $b$-hadron branching fractions for two-body decays into charmless charged hadrons, JHEP 10 (2012) 037, arXiv:1206.2794
  * (26) M. Adinolfi et al., Performance of the LHCb RICH detector at the LHC, arXiv:1211.6759, to appear in Eur. Phys. J. C
  * (27) E. Baracchini and G. Isidori, Electromagnetic corrections to non-leptonic two-body B and D decays, Phys. Lett. B633 (2006) 309, arXiv:hep-ph/0508071
  * (28) ARGUS collaboration, H. Albrecht et al., Search for $b\rightarrow s\gamma$ in exclusive decays of B mesons, Phys. Lett. B229 (1989) 304
  * (29) S. Bianco, F. L. Fabbri, D. Benson, and I. Bigi, A Cicerone for the physics of charm, Riv. Nuovo Cim. 26N7 (2003) 1, arXiv:hep-ex/0309021
  * (30) Belle collaboration, M. Staric et al., Search for a CP asymmetry in Cabibbo-suppressed $D^{0}$ decays, Phys. Lett. B670 (2008) 190, arXiv:0807.0148
  * (31) BaBar collaboration, B. Aubert et al., Search for CP violation in the decays $D^{0}\rightarrow K^{-}K^{+}$ and $D^{0}\rightarrow\pi^{-}\pi^{+}$, Phys. Rev. Lett. 100 (2008) 061803, arXiv:0709.2715
  * (32) CDF collaboration, T. Aaltonen et al., Measurement of CP–violating asymmetries in $D^{0}\rightarrow\pi^{+}\pi^{-}$ and $D^{0}\rightarrow K^{+}K^{-}$ decays at CDF, Phys. Rev. D85 (2012) 012009, arXiv:1111.5023
  * (33) LHCb collaboration, R. Aaij et al., Evidence for $C\\!P$ violation in time-integrated $D^{0}\rightarrow h^{-}h^{+}$ decay rates, Phys. Rev. Lett. 108 (2012) 111602, arXiv:1112.0938
  * (34) Heavy Flavor Averaging Group, Y. Amhis et al., Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of early 2012, arXiv:1207.1158, update available online at http://www.slac.stanford.edu/xorg/hfag
  * (35) 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