Search for $CP$ violation in $\Xi_b^- \to p K^- K^-$ decays

A search for $CP$ violation in charmless three-body $\Xi_b^- \to p K^- K^-$ decays is performed using $pp$ collision data recorded with the LHCb detector, corresponding to integrated luminosities of $1\,\text{fb}^{-1}$ at a centre-of-mass energy $\sqrt{s} = 7\,\text{TeV}$, $2\,\text{fb}^{-1}$ at $\sqrt{s} = 8\,\text{TeV}$ and $2\,\text{fb}^{-1}$ at $\sqrt{s} = 13\,\text{TeV}$. A good description of the phase-space distribution is obtained with an amplitude model containing contributions from $\Sigma(1385)$, $\Lambda(1405)$, $\Lambda(1520)$, $\Lambda(1670)$, $\Sigma(1775)$ and $\Sigma(1915)$ resonances. The model allows for $CP$-violation effects, which are found to be consistent with zero. The branching fractions of $\Xi_b^- \to \Sigma(1385) K^-$, $\Xi_b^- \to \Lambda(1405) K^-$, $\Xi_b^- \to \Lambda(1520) K^-$, $\Xi_b^- \to \Lambda(1670) K^-$, $\Xi_b^- \to \Sigma(1775) K^-$ and $\Xi_b^- \to \Sigma(1915) K^-$ decays are also reported. In addition, an upper limit is placed on the product of ratios of $\Omega_b^-$ and $\Xi_b^-$ fragmentation fractions and the $\Omega_b^- \to p K^- K^-$ and $\Xi_b^- \to p K^- K^-$ branching fractions.


Introduction
In the Standard Model (SM), CP violation, defined as the breaking of symmetry under the combined charge conjugation and parity operations, owes its origin to a single irreducible complex phase in the Cabibbo-Kobayashi-Maskawa (CKM) matrix [1,2]. All effects of CP violation in particle decays observed so far are consistent with this paradigm. However, the degree of CP violation permitted in the SM is inconsistent with the observed matter-antimatter asymmetry in the Universe [3,4]. This motivates further searches for sources of CP violation beyond the SM.
Interference between two amplitudes with different weak and strong phases leads to CP violation in decay, where weak phases are those that change sign under CP conjugation while strong phases do not. In the SM, weak phases are associated with the complex elements of the CKM matrix and strong phases are associated with hadronic final-state effects. Two such amplitudes are potentially present in decays of b hadrons to final states that do not contain charm quarks, which therefore provide fertile ground for studies of CP violation. Significant asymmetries have been observed between B and B partial widths in B 0 → K − π + [5-9] B 0 → π + π − [5, 6, 10] and B 0 s → K + π − [7,8] decays. Even larger CP -violation effects have been observed in regions of the phase space of B − decays to π + π − π − , K − π + π − , K + K − K − and K + K − π − final states [11][12][13][14][15][16].
Breaking of CP symmetry has not yet been observed in the properties of any baryon. Tests of this symmetry have been performed through studies of Λ 0 b baryon decays to pπ − , pK − [7, 17], K 0 S pπ − [18], ΛK + K − , ΛK + π − [19], pπ − π + π − , pπ − K + K − , pK − π + π − and pK − K + K − [20-22] final states, as well as Ξ 0 b decays to pK − π + π − and pK − π + K − [21,22]. No significant evidence of CP violation has been found in any of these studies, nor in measurements of the properties of charm baryon decays [35]. In light of the large CP -violation effects observed in three-body charmless decays of B mesons, it is of great interest to extend the range of searches in b-baryon decays. In particular, the recently observed Ξ − b → pK − K − decay [23] provides an interesting new opportunity to search for CP -violation effects.
In this paper, the first amplitude analysis of Ξ − b → pK − K − decays is reported. This is also the first amplitude analysis of any b-baryon decay mode allowing for CP -violation effects. A search for the previously unobserved Ω − b → pK − K − decay is also presented. The analysis reported here is performed using proton-proton (pp) collision data recorded with the LHCb detector, corresponding to integrated luminosities of 1 fb −1 at a centre-of-mass energy of This paper is organised as follows. Section 2 gives a brief description of the LHCb detector, trigger requirements and simulation software. The signal candidate selection procedure is set out in Sec. 3. In Sec. 4, the procedure for estimating the signal and background yields that enter the amplitude fit is explained. Section 5 covers the modelling of the distribution of decays across the phase space. Sections 6 and 7 contain a description of the systematic uncertainties associated with the analysis procedure and a presentation of the results, respectively. A brief summary of the analysis is given in Sec. 8.

Detector, trigger and simulation
The LHCb detector [24,25] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV. 1 The minimum distance of a track to a primary pp collision vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam direction, in GeV. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The magnetic field deflects oppositely charged particles in opposite directions and this can lead to detection asymmetries. Periodically reversing the magnetic field polarity throughout the data-taking reduces this effect to a negligible level. Approximately 60% of 2011 data, 50% of 2012 data, 61% of 2015 data and 53% of 2016 data were collected in the "down" polarity configuration, and the rest in the "up" configuration.
The online event selection is performed by a trigger, which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. During offline analysis, reconstructed candidates are associated with trigger decisions. Events considered in the analysis are required to have been triggered at the hardware level in one of two ways: either through one of the final-state tracks of the signal decay depositing sufficient energy in the calorimeter system, or by one of the other tracks in the event, not reconstructed as part of the signal candidate, fulfilling any hardware trigger requirement. At the software stage, it is required that at least one charged particle associated to the b-hadron candidate has high p T and high χ 2 IP , where χ 2 IP is defined as the difference in PV fit χ 2 with and without the inclusion of a specific particle. A multivariate algorithm [26] is used to identify secondary vertices consistent with being a two-or three-track b-hadron decay. The PVs are fitted with and without the tracks that comprise the b-baryon candidate, and the PV that gives the smallest χ 2 IP is associated with the candidate. Finally, the momentum scale for charged particles is calibrated using samples of J/ψ → µ + µ − , B + → J/ψK + and Λ → pπ − decays collected concurrently with the data sample used for this analysis [27,28].
Simulation samples are used to investigate background from other b-hadron decays and to study the detection and reconstruction efficiency of the signal. In the simulation, pp collisions are generated using Pythia [29] with a specific LHCb configuration [30]. Decays of unstable particles are described by EvtGen [31], in which final-state radiation is generated using Photos [32]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [33] as described in Ref. [34].

Offline selection
The offline selection consists of an initial filtering stage followed by a requirement on the output of a multivariate algorithm (MVA). Compared to the procedure applied to select the Ξ − b → pK − K − channel in Ref. [23], improvements in both stages lead to a significant increase in efficiency. In particular, the inclusion in the multivariate algorithm of particle identification (PID) variables that distinguish the final-state charged hadrons from misidentified particles, is found to separate signal from background effectively.
In the filtering stage, tracks are required to be of good quality, to satisfy p > 1500 MeV and p T > 250 MeV, and to be displaced from all PVs. Tracks associated to proton candidates must, at this stage, satisfy a loose PID requirement and all tracks are required to not be associated to hits in the muon system. Each b-hadron (henceforth denoted as X − b ) candidate must form a good-quality decay vertex that is separated significantly from any PV, and must be consistent with originating from its associated PV. Only X − b candidates with p T > 3500 MeV and invariant mass 5545 < m(pK − K − ) < 6470 MeV are retained for further analysis.
In the selected m(pK − K − ) range there are three main categories of background that contribute: combinatorial background that results from random association of unrelated tracks; partially reconstructed background due to b-hadron decays into final states similar to the signal, but with additional soft particles that are not reconstructed; and cross-feed background that results from misidentification of one or more final-state particles. The MVA classifier is designed primarily to reduce combinatorial background while retaining high signal efficiency, but also has some discriminating power against the other background sources. It is trained with a signal sample comprised of simulated Ξ − b → pK − K − decays generated uniformly across the phase space, and a background sample obtained from candidates in data in the sideband regions 5545.0 < m(pK − K − ) < 5634.4 MeV and 6209.0 < m(pK − K − ) < 6470.0 MeV. The latter of these regions is dominated by combinatorial background, as its lower threshold excludes possible Ω − b → pK − K − decays from the sample. The former region includes also contributions from sources of partially reconstructed background such as decays is removed by assigning the proton candidate the kaon mass and vetoing the m(K + K − K − ) region within ±45 MeV around the known B − mass [35]. This veto corresponds to approximately ±3 times the invariant mass resolution for Variables that exhibit good discriminating power between the signal and background samples are chosen as inputs to the MVA. These are: the angle between the X − b candidate's momentum vector and the line connecting its decay vertex to its associated PV; the scalar sum of the p T of all final-state tracks; the χ 2 IP of the highest p T final-state track and of the X − b candidate; the square of the significance of the distance between the X − b decay vertex and its associated PV; the vertex fit χ 2 per degree of freedom of the X − b candidate; the minimum change in the X − b candidate vertex fit χ 2 when including an additional track; variables that characterise the PID information of the proton and kaon candidates; and a variable that quantifies the isolation of the X − b candidate. The last of these is defined as the p T asymmetry between the X − b candidate and the tracks within a circle, centred on the X − b candidate (but excluding its decay products), with a radius δη 2 + δφ 2 < 1.7 in the space of pseudorapidity η and azimuthal angle φ (in radians) around the beam direction [36].
To describe accurately the proton and kaon PID variables, the quantities in simulation are resampled according to values obtained from data calibration samples of Λ 0 b → Λ + c π − , D + s → φπ + and D * + → D 0 π + decays [37]. The procedure accounts for correlations between the variables associated to a particular track, as well as the dependence of the PID response on p T , η, and the multiplicity of tracks in the event. All other MVA input variables show good agreement between simulation and data, as validated with a control sample of B − → ppK − decays. The MVA input variables are also found to not be correlated strongly either with X − b candidate mass or with position in the phase space of the decay. Several types of MVA classifier are investigated, with a gradient boosted decision tree algorithm giving the best performance [38]. Four classifiers are trained separately with samples separated by data-taking period (Run 1 or Run 2) and by even or odd event numbers. The event number identifies the proton-proton bunch crossing, from which the X − b candidate was recorded, in a certain operational period of the experiment. To avoid possible MVA overtraining, for each data-taking period the classifier trained on the sample with even event numbers is validated and employed on the sample with odd event numbers, and vice versa.
A threshold on the output of the MVA is chosen to maximise N S / √ N S + N B , where N S and N B represent the estimated numbers of Ξ − b → pK − K − signal and combinatorial background candidates, respectively, within a signal region of ±40 MeV around the Ξ − b mass from Ref. [39]. This range corresponds to approximately ±2.5 times the Ξ − b → pK − K − invariant mass resolution. The value of N S is estimated using the signal efficiency evaluated from simulation, multiplied by the Ξ − b → pK − K − branching fraction, the Ξ − b fragmentation fraction, the bb production cross-section [40] and the integrated luminosity for the relevant data-taking period. The product of the Ξ − b → pK − K − branching fraction and the Ξ − b fragmentation fraction is obtained from the results of Ref. [23], where the B − → K + K − K − channel is used for normalisation, by multiplying the B − fragmentation fraction in the relevant kinematic range [41] and the B − → K + K − K − branching fraction [35]. The value of N B is estimated from data by fitting the region 6125 < m(pK − K − ) < 6470 MeV with a linear function and extrapolating the result into the signal region. The MVA output requirements have efficiencies of about 52% and 61% for Run 1 and Run 2, respectively, with combinatorial background rejection of about 98% for both data-taking periods. The choice of N S / √ N S + N B as the figure of merit is intended to obtain a sufficiently large data sample to make an amplitude analysis viable. After all offline selection requirements are applied, each selected event contains a single X − b candidate. The variables describing the phase space of the decay, which are used in the amplitude analysis, are calculated following a kinematic fit in which the X − b candidate mass is fixed to the Ξ − b mass from Ref. [39]. This procedure improves resolution of these variables and ensures that all decays remain within the phase-space boundary. The difference between the Ξ − b mass value used in this fit and recent more precise results [42][43][44] has negligible impact on the analysis. The experimental resolution of the m(pK − ) invariant mass, in the region with the narrowest resonance considered in this analysis, the Λ(1520) state, is expected to be around 1.5 MeV. This is smaller than the Λ(1520) width, and therefore effects related to finite resolution in the phase-space variables are not considered further.
The expected Ξ − b → pK − K − signal efficiency, assuming uniform distribution of decays across the phase space and taking into account the LHCb detector acceptance, reconstruction and both online and offline selection criteria, is (1.159 ± 0.005) % for Run 1 and (1.748 ± 0.006) % for Run 2. The corresponding Ω − b → pK − K − signal efficiencies are (1.257 ± 0.005) % and (1.921 ± 0.006) %. The quoted uncertainties are due to the limited size of the simulation samples only.

X − b candidate mass fit
Distributions of m(pK − K − ) for selected X − b candidates are shown in Fig. 1 for Run 1 and Run 2 separately. The signal yields are obtained from unbinned extended maximumlikelihood fits to these distributions. The fit model is composed of signal and background components whose shape parameters are mostly obtained from fits to the corresponding simulation samples, after imposing the same selection requirements as on the data. One exception is the combinatorial background component, which is modelled by an exponential function with slope parameter allowed to vary freely in the fit to data.
Signal Ξ − b → pK − K − and Ω − b → pK − K − components are each modelled with the sum of two Crystal Ball (CB) functions [45] where the core width and peak position are shared and with independent power-law tails on both sides. The tail parameters and the relative normalisation of the CB functions are determined from simulation. The peak positions are fixed to the Ξ − b mass from Ref. [42] and the known Ω − b mass [35], and a scale factor relating the width in data to that in simulation is introduced.
A possible cross-feed background contribution from Ξ − b → pK − π − decays [23], where the pion is misidentified as a kaon, is modelled with the sum of two CB functions. All shape parameters of this function are fixed according to the values obtained from a fit to simulation but the width is scaled by the same factor as the signal components. The phase-space distribution of these decays is not known, and the simulation sample is weighted according to a model, inspired by the m(pK − ) and m(pπ − ) mass spectra observed in Λ 0 b → J/ψpK − [46] and Λ 0 b → J/ψpπ − [47] decays, which consists of the Λ(1405), Λ(1520), Λ(1690), N (1440), N (1520), N (1535) and N (1650) resonances. The yield of the Ξ − b → pK − π − cross-feed component is expressed relative to the Ξ − b → pK − K − signal yield and constrained within uncertainty according to the previous branching fraction ratio measurement [23] and relative selection efficiency. The expected relative yields are 0.15 ± 0.06 and 0.14 ± 0.03 for Run 1 and Run 2, respectively.
Partially reconstructed and combinatorial background contributions are also included in the fit model. It is found that the m(pK − K − ) distributions of various potential sources of partially reconstructed background, such as [46]. Therefore, the baseline fit model includes a single partially reconstructed background component, which is modelled from simulated Ξ − b → K * − (K − π 0 )K − p decays with an ARGUS function [48] convolved with a Gaussian function. The threshold of the ARGUS function is fixed to the known value of m Ξ − b − m π 0 [35,42], and the width parameter of the Gaussian function is taken from the fit to simulation and scaled by the same factor as the signal components. Negligible contributions are expected from partially reconstructed Ω − b decays, such as

Parameter
Run 1 15 ± 9 Partially reconstructed background yield 231 ± 34 442 ± 36 Combinatorial background yield 721 ± 50 775 ± 51 The results of the fits to Run 1 and Run 2 data are shown in Table 1 and Fig. 1. The free parameters of each fit are the two signal yields, the partially reconstructed and combinatorial background yields, the width scale factor and the exponential shape parameter of the combinatorial background, while the cross-feed background yield is constrained to its expectation relative to the Ξ − b → pK − K − signal yield. Only candidates in the m(pK − K − ) signal region of ±40 MeV around the Ξ − b mass from Ref. [39] are retained for the amplitude analysis. In this region, the yields of the signal, cross-feed and combinatorial components are N sig = 181 ± 20, N cf = 16 ± 7 and N comb = 90 ± 6 for Run 1, and N sig = 278 ± 21, N cf = 25 ± 6 and N comb = 95 ± 6 for Run 2, where the quoted uncertainties are statistical only. These correspond to signal purities of (63 ± 3)% and (70 ± 2)% for Run 1 and Run 2, respectively. The contribution from partially reconstructed background in the signal region is negligible.
No significant signal from the Ω − b → pK − K − decay is observed. The results of the fits are used to set limits on the product of its branching fraction with the fragmentation fraction for Ω − b production, normalised to the corresponding quantities for where N and denote yield and efficiency, respectively, for the indicated mode, while and Ω − b fragmentation fractions. Results for the ratio R are reported, both for Run 1 and Run 2 separately and combined, in Sec. 7.

Amplitude analysis
The phase space of the three-body decay of a, potentially polarised, b baryon has five degrees of freedom. A baseline assumption is made that Ξ − b baryons produced in pp collisions within the LHCb acceptance have negligible polarisation, as observed for Λ 0 b baryons [49,50]. As a result, the phase space of the Ξ − b → pK − 1 K − 2 decay is characterised by two independent kinematic variables (subscripts here distinguish the two kaons in the final state). Since no resonances are expected to decay to K − 1 K − 2 , these variables chosen are the squared invariant masses m 2 (pK − 1 ) and m 2 (pK − 2 ). The presence of two identical kaons in the final state imposes a Bose symmetry that the decay amplitudes must be invariant under the exchange of these two particles. As a result, the two-dimensional distribution of m 2 (pK − 1 ) and m 2 (pK − 2 ) has a symmetry under interchange of the variables. This motivates the use of the variables m 2 low and m 2 high , which denote the lower and higher of m 2 (pK − 1 ) and m 2 (pK − 2 ), respectively, effectively removing a duplicated half of the m 2 (pK − 1 ), m 2 (pK − 2 ) plane. References hereafter to the Dalitz plot (DP) of Ξ − b → pK − K − decays refer to the two-dimensional m 2 low , m 2 high distribution. The DP distributions of selected candidates in Run 1 and Run 2 are shown in the top row of Fig. 2.
It is common practice in amplitude analysis to use the so-called "square" Dalitz plot (SDP) variables [15,51], which in this case are defined as Here m min (K − K − ) = 2m K and m max ( is the angle between one K − direction in the K − K − centre-of-mass frame and the direction of the K − K − system in the Ξ − b centre-of-mass frame. The symmetry of the final state requires that distributions are symmetric with respect to θ = 0.5, so only the region θ ∈ [0, 0.5] is considered. These SDP variables provide improved granularity, when using uniform binning, in the regions close to the DP boundaries that tend to be populated most densely. This is beneficial for example in the modelling of the signal efficiency. Furthermore, the mapping to a square space aligns the bin boundaries to the kinematic boundaries of the phase space. As such, all efficiencies and background distributions in the analysis are obtained as functions of the SDP variables. The SDP distributions of selected candidates in Run 1 and Run 2 are shown in the bottom row of Fig. 2.

Modelling of the signal component
The probability density function (PDF) for the signal component is expressed as where Q = +1 for Ξ − b decays and Q = −1 for Ξ + b decays and Ω denotes the phase space in terms of the DP variables. The efficiency is denoted by Q (Ω), and can differ for Q = +1 and −1 to accommodate efficiency asymmetries; as described in Sec. 5.2, the efficiency maps are determined using SDP coordinates, denoted Ω , but at any point in the phase space Q (Ω) = Q (Ω ). The term dΓ Q /dΩ describes the differential decay densities for Ξ − b and Ξ + b decays, including both local and overall rate asymmetries and the normalisation factor Γ is Equations (3) and (4) assume no asymmetry in the production rates of Ξ − b and Ξ + b baryons produced within the LHCb acceptance from high-energy pp collisions, consistent with measurement [44]. The effect of such a production asymmetry would, in this analysis, mimic a global (i.e. phase-space independent) difference between Q=+1 and Q=−1 , and therefore the systematic uncertainty due to this assumption can be evaluated straightforwardly.
The differential decay density is expressed as where A Q R,M Ξ b ,λp denotes the symmetrised decay amplitude for a given intermediate state R, Ξ − b spin component along a chosen quantisation axis M Ξ b , and proton helicity λ p . The quantisation axis is chosen to be the direction opposite to the proton momentum in the Ξ − b rest frame, and the proton helicity is defined in the rest frame of the K − K − system to ensure explicit symmetry between the pK − low and pK − high decay chains. Here K − low is the kaon whose four-momentum is used in the definition of m 2 low and K − high denotes the other kaon. The amplitude in Eq. (5) has been summed incoherently over the spins of the initial and final states (corresponding to an average over initial states) and coherently over all contributing intermediate states.
The helicity formalism is used to parametrise the decay dynamics. A detailed description of this formalism can be found in Refs. [46,47,[52][53][54]. In particular, the Dalitz-plot decomposition procedure [54] is followed to express the symmetrised decay amplitude as The first term corresponds to the amplitude for the weak decay where R decays to pK − low via the strong interaction. This decay amplitude is expressed as where the amplitude is summed coherently over the allowed helicities of the intermediate state, λ R , and of the proton, λ p , defined in the rest frames of the Ξ − b and the pK − low systems, respectively. The Ξ − b , R and proton spins are denoted by J Ξ b , J R and J p , respectively. The three functions of the form d J λ,λ in Eq. (7) are the small Wigner d-matrix elements [55] that impose angular momentum conservation giving rise to the condition |λ R | ≤ 1/2. As a result, for intermediate states with any half-integer spin, only helicities corresponding to λ R = ±1/2 contribute to the amplitude. The three angles θ R , θ p and ζ are functions of the DP variables. The angle θ R , defined in the Ξ − b rest frame, is formed between the direction opposite to the proton momentum and the combined momentum of the pK − low system. The angle θ p is between the direction opposite to the K − high momentum in the Ξ − b rest frame and the proton momentum in the rest frame of the pK − low system. The angle ζ gives the Wigner rotation that is required to relate the proton helicity state, |λ p , defined in the pK − low rest frame to the proton helicity state, |λ p , defined in the K − K − rest frame. This angle, computed in the proton rest frame, is formed between the momenta of K − low and of the K − K − system. Mathematical definitions of these three angles, each defined in the range [0, π], are where high , the Källén function is given by K(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc) and the Ξ − b , K and p masses are denoted by m Ξ b , m K and m p , respectively.
The second term in Eq. (6) corresponds to the weak decay of The expression for this amplitude can be obtained by interchanging m 2 low ↔ m 2 high in Eqs. (7)-(10). The term η λ p , in Eq. (7), arises as a consequence of parity conservation in the strong decay of the intermediate state R. It is defined as (7), encapsulates the combined couplings of the weak decay of the initial state and the strong decay of the intermediate state. This coefficient, subsequently referred to as the helicity coupling, can be expressed as where x R,λ R and y R,λ R denote the real and imaginary components of the CP -conserving part of the coupling, while δx R,λ R and δy R,λ R are CP -violating parameters.
The term R(m 2 low ), in Eq. (7), describes the lineshape of each resonant or nonresonant contribution. Resonances are parametrised with relativistic Breit-Wigner (RBW) functions, F RBW , that are modified by Blatt-Weisskopf barrier factors, B L Ξ b and B L R , and are given by Here m 2 x is either m 2 low or m 2 high , while p is the magnitude of the resonance momentum in the Ξ − b centre-of-mass frame and q is the magnitude of the proton momentum in the resonance centre-of-mass frame. The symbols p 0 and q 0 denote the values of these quantities at the resonance peak, i.e. when m x = m 0 . The orbital angular momentum released in the Ξ − b decay is denoted L Ξ b while that in the resonance decay is denoted L R . Angular momentum conservation in the Ξ − b decay imposes the condition x ). Angular momentum conservation in the resonance decay limits L R to J R ± 1 2 , which is then uniquely defined by parity conservation in the decay, η R = (−1) L R +1 .
The Blatt-Weisskopf barrier functions are and account for suppression creating high values of the orbital angular momentum L, which depends on the momentum of one of the decay products, k, in the centre-of-mass frame of the decaying particle and on the size of the decaying particle given by the constant d. The value d = 5.0 GeV −1 is used for Ξ − b decays while 1.5 GeV −1 is used for resonances [53].
The relativistic Breit-Wigner amplitude is given by where Here m 0 and Γ 0 denote the pole mass and width of the resonance, respectively. In the case of the Λ(1405) resonance, which peaks below the pK − threshold, m 0 is replaced by an effective mass in the kinematically allowed region [56], where m max and m min are the upper and lower limits of the kinematically allowed range, respectively. In this case, the q 0 value in Eq. (18) is the value of q at m = m eff 0 . This parameterisation ensures that only the tail of the RBW function enters the fit model as a virtual contribution. Nonresonant components are modelled using an exponential lineshape, where α is a slope parameter that is determined from the fit and m 0 is fixed to be the midpoint of the m low range, i.e. 2.83 GeV. The primary outputs of the amplitude analysis are the CP -conserving and CP -violating components of the helicity couplings introduced in Eq. (11). However, since these depend on the choice of phase convention, amplitude formalism and normalisation, they can be difficult to compare between analyses. It is therefore more useful to report the fit fractions F i for each intermediate component i of the fit model, defined by where It is also useful to report the interference fit fractions I ij between the two intermediate components i and j, defined by where The parameters of CP violation, A CP i , associated with each component i of the model are also reported. These are defined as

Modelling of signal efficiency and background distributions
The detector geometry and the online and offline selection procedure can induce variation in the signal efficiency across the phase space of the decay. This is accounted for, as shown in Eq. (3), by determining the efficiency as a function of the SDP variables. The efficiency maps are obtained from simulation, but with effects related to PID calibrated using data as outlined in Sec. 3. The efficiency maps for Ξ − b and Ξ + b decays can be seen in Fig. 3 separately for Run 1 and Run 2. These maps are obtained by employing a uniform 10 × 10 binning scheme and smoothing with a two-dimensional cubic spline to mitigate effects of discontinuity at the bin edges. No significant detection asymmetry is observed.
Candidates selected in data in the sideband 5890 < m(pK − K − ) < 6470 MeV are used to model the SDP distribution of the combinatorial background, which dominates this region as discussed in Sec. 4. The effect of the Ξ − b mass constraint used when calculating the SDP variables causes a distortion of the distribution from that of combinatorial background in the signal region. This is accounted for using a method [57] in which the unnormalised function describing the m(pK − K − ) and Ω = (m , θ ) space is expressed as where β is a free parameter determined from a fit to the sideband candidates. The functions f 0 and f 1 are modelled using neural networks that are trained using candidates from the data sideband region. This model is then extrapolated to predict the PDF of the combinatorial background at the Ξ − b mass, i.e. P comb (Ω ) = F (m Ξ b , Ω )/N , where the normalisation factor N = Ω F (m Ξ b , Ω )dΩ . The PDF in terms of DP variables is obtained using P comb (Ω) = |J|P comb (Ω ), where |J| is the Jacobian determinant of the transformation between variables, dΩ = |J|dΩ . The SDP distributions of the combinatorial background for Run 1 and Run 2 are shown in the top row of Fig. 4.
After imposing the selection criteria and splitting the data sample by the initial state charge, too few candidates are available in the sideband region to train the neural networks. As a result, the neural networks are trained using the combined sample of Ξ − b and Ξ + b candidates, and no asymmetry in the shape of the combinatorial background SDP distribution is assumed in the baseline model.
The SDP distribution of cross-feed background from misidentified Ξ − b → pK − π − decays that enter the signal region is modelled using simulation, and is shown in the bottom row of Fig. 4 separately for Run 1 and Run 2. These distributions are described in terms of a uniform 10 × 10 binned SDP histogram, smoothed with a two-dimensional cubic spline. As described in Sec. 4, the simulation is weighted to reproduce resonance structures expected in the phase space of Ξ − b → pK − π − decays. Differences in selection requirements, together with the limited statistics of the Ξ − b → pK − π − simulation samples, cause the PDFs to differ between Run 1 and Run 2. In the baseline fit it is assumed that there is no asymmetry between Ξ − b and Ξ + b candidates in the cross-feed yields or SDP distributions.

Fitting procedure
The total PDF that is used to model the phase-space distributions of Ξ − b → pK − K − and its conjugate decay is The yields of signal, combinatorial background and cross-feed background components are denoted by N sig , N comb and N cf , respectively, and are obtained as described in Sec. 4, separately for Run 1 and Run 2. The quantity N tot = N sig + N comb + N cf is the total yield in the signal region. The PDFs for the signal, combinatorial background and cross-feed background components are denoted by P Q sig (Ω), P comb (Ω) and P cf (Ω), respectively, where the former is given in Eq. (3) and the latter two are displayed in Fig. 4 in terms of the SDP variables. In the baseline model, only the signal PDF can differ for Ξ − b and Ξ + b candidates, although a possible global combinatorial background asymmetry, A comb , is a free parameter of the model. An unbinned maximum-likelihood fit is performed to the combined sample of candidates for Ξ − b → pK − K − and its conjugate decay to determine the parameters of the model, which are the CP -conserving and CP -violating coefficients of the helicity couplings of Eq. (11). The fit is performed simultaneously to the Run 1 and Run 2 data samples, which have separate efficiency and background models as described above. The fit model is implemented in a fitting package based on TensorFlow [58], interfaced with the Minuit function minimisation algorithm [59,60]. The function that is minimised, twice the negative log-likelihood, is where the index i runs over the N r candidates in the data sample from run period r (Run 1 or Run 2), and Ω i denotes the DP coordinates of candidate i.

Model selection
The Ξ − b → pK − K − decay can proceed via intermediate pK − resonances. Various Λ * and Σ * resonances that are known to decay to pK − are considered as potential components of the signal model. 2 The Particle Data Group (PDG) [35] reports a large number of such states; those that are sufficiently well established are considered in this study and are shown in Table 2. Masses and widths of all resonance components are fixed to either the central value or the midpoint in the range of values quoted in Table 2. Nonresonant components, labelled NR(J P ), with spin-parity J P = 1 Resonances with spin J ≥ 7/2 are excluded from consideration as these would require L Ξ b ≥ 3 and are therefore expected to be significantly suppressed. Table 2: Summary of the considered Λ * and Σ * resonances, ranked either **** or *** by the PDG [35]. Note that the pK − threshold is at 1432 MeV. Resonances marked † are included in the baseline model, as described in the text. The spin-parity of the Σ(2250) is not known and is assumed to be 3 2 + . For many of these states, the PDG does not report masses and widths with central values and uncertainties, but rather gives real and imaginary parts of the pole position. This reflects the fact that a simple Breit-Wigner parameterisation of these resonances may not fully describe their lineshapes; however, more sophisticated parametrisations are beyond the scope of the current analysis. In this analysis, as the normalisation of the decay density is arbitrary, the Λ(1520) resonance is chosen as the reference component. This implies that the coupling with the positive helicity of the Λ(1520) resonance a Q R,λ R =+1/2 is real. Explicitly, for the reference Λ(1520) resonance, y R,λ R =+1/2 = δy R,λ R =+1/2 = 0 and x R,λ R =+1/2 = 1, while δx R,λ R =+1/2 is free to vary in the fit to allow for CP violation in the Λ(1520) amplitude. The analysis is found to be insensitive to the coupling of the Λ(1520) component with negative helicity, and therefore x R,λ R =−1/2 = δx R,λ R =−1/2 = y R,λ R =−1/2 = δy R,λ R =−1/2 = 0 for the reference resonance. The helicity couplings of all other resonant and nonresonant components are left free to vary in the fit.
To establish a baseline fit model, the Λ(1520) component alone is initially included in the model, with additional components added iteratively in the order that maximises the change in −2 ln L obtained from fits to the data with prospective models. Components with different spin and parity should have zero interference fit fractions due to the orthogonality relation satisfied by the small Wigner d-matrix elements. However, the symmetrisation of the Dalitz plot can lead to non-zero values for such interference fit fractions in this analysis. As a result, in establishing the baseline model, it is possible to encounter "unphysical" interference fit fractions (> 40%) between two components. When such a case occurs, the component that gives the minimal change in −2 ln L when removed from the fit model is discarded. The procedure is terminated when the change in −2 ln L from including any further contribution is less than 9 units, limiting the potential for the model to be influenced by statistical fluctuations. This approach leads to a model that contains Σ(1385), Λ(1405), Λ(1520), Λ(1670), Σ(1775) and Σ(1915) components. The potential for additional components to be present in the true underlying model is considered as a source of systematic uncertainty.

Fit to data
In an attempt to find the global minimum, a large number of fits to data are performed, where the initial values of the helicity couplings are randomised. The baseline results are obtained from the fit that returns the smallest −2 ln L value out of this ensemble. This procedure is found to converge successfully to the global minimum without any secondary minima.
The Dalitz-plot distribution of the combined Run 1 and Run 2 data sample is compared to the model obtained from the fit to data, separately for Ξ − b and Ξ + b candidates, in Fig. 5. Projections of the fit results onto m low and m high are compared to the data in Fig. 6. Further comparisons of the fit result and the data in regions of the phase space are presented in Appendix A. There is no indication of CP violation, i.e. no significant difference between Ξ − b and Ξ + b decays, in the distributions. The overall agreement between the data and the model is good, with unbinned goodnessof-fit tests using the mixed-sample and point-to-point dissimilarity approaches [61] giving p-values of 0.20 and 0.25, respectively. In Fig. 6, there is an apparent discrepancy between the model and the data at m high between 3.4 GeV to 3.7 GeV, predominantly in the Ξ + b sample. However, due to the symmetry of the final state, any structure that appears in m high should also appear in m low , where no such structure is observed. Addition of extra components to the fit model does not significantly improve the data description. Moreover, the apparent discrepancy in Fig. 6 does not take into account the systematic uncertainty in the mismodelling of the combinatorial background, which is the largest component at m high ∼ 3.7 GeV. Therefore, this feature is not considered to be significant and is not investigated further.

Systematic uncertainties
The outcomes of the analysis are the ratio R of Ω − b and Ξ − b branching fractions and fragmentation fractions (see Eq. (1)), as well as the fit fractions, interference fit fractions and CP -asymmetry parameters obtained from the amplitude analysis. Various sources   of systematic uncertainty can affect these measurements. These are discussed one by one in this section, concluding with a summary. Systematic uncertainties are evaluated separately for Run 1 and Run 2, where appropriate.

Invariant mass fits
The fits to the m(pK − K − ) invariant mass distributions determine the signal and background yields, which are used in both the calculation of R and the amplitude analysis. Four sources of systematic uncertainty arising from these fits are considered. The first is due to the limited size of the data sample. This enters the calculation of R as statistical uncertainty, but is a source of systematic uncertainty in the amplitude analysis where the signal and background yields are fixed parameters. To evaluate the associated systematic uncertainty, these yields are varied according to the covariance matrix obtained from the m(pK − K − ) fit and for each variation the fit to the phase-space distribution is repeated. The root mean square (RMS) of the distribution of the change in each result of the amplitude analysis is assigned as the corresponding systematic uncertainty. The second source relates to the m(pK − K − ) fit model. Models for each of the components are varied to evaluate associated systematic uncertainties: the signal model is replaced with a Hypatia function [62]; the combinatorial background model is replaced with a second-order Chebyshev polynomial; the cross-feed and partially reconstructed background models are replaced with kernel density estimates; additional partially recon- structed background components are included; the model used to describe the phase-space distribution of the Ξ − b → pK − π − cross-feed background is varied. In each case, the change in each result from its baseline value is taken as the systematic uncertainty.
The third source concerns fixed parameters in the m(pK − K − ) fit that are taken from simulation. An ensemble of pseudoexperiments is generated using the nominal values of these parameters, and each is then fitted many times with parameters fixed to alternative values obtained using the covariance matrices of the fits to the simulation samples. The standard deviation of the change in each result is evaluated for every pseudoexperiment, and its average value over the ensemble is assigned as the systematic uncertainty.
Finally, potential fit bias is investigated by generating multiple pseudoexperiments with yield and fit parameters obtained from the nominal m(pK − K − ) fit. The difference between the mean fit result of the ensemble and the nominal value is assigned as the associated systematic uncertainty.

Selection efficiency maps
The efficiency maps are altered to evaluate systematic uncertainties from six sources. In each case the difference between the results obtained using the alternative efficiency maps and that with the baseline efficiency maps is assigned as the systematic uncertainty.
The first source reflects uncertainties in the p T distribution of Ξ − b baryons produced in the LHCb acceptance. Alternative efficiency maps are obtained where the simulation samples are weighted so that the p T distribution matches that of the background-subtracted data. Since there is no significant signal of Ω − b baryons, they are assumed to have the same p T distribution as Ξ − b baryons and the Ω − b → pK − K − efficiency map is altered in the same way.
The second source is a possible mismatch of the hardware trigger efficiency between simulation and data, which could arise due to miscalibration of transverse energy measurements from the calorimeter. Alternative efficiency maps are obtained by applying corrections that are calculated, as a function of track p T , using control samples of kaons from D * + → D 0 (K − π + )π + decays and protons from Λ 0 b → Λ + c (pK − π + )π − decays. The third source is due to uncertainty arising from binning the phase space when evaluating the efficiency maps. Alternative efficiency maps are obtained employing a different SDP binning scheme. An additional systematic uncertainty is associated to the efficiency of Ω − b → pK − K − decays, and arises from their unknown phase-space distribution. The standard deviation of the variation of the efficiency across the binned SDP histogram is assigned as the corresponding uncertainty.
The remaining sources relate to particle identification. The PID variables used in the MVA are drawn from data calibration samples accounting for dependence on the signal kinematics. Systematic uncertainties in this procedure arise from the limited statistics of both the simulation and calibration samples, and the modelling of the PID variables in the calibration samples. The limitations due to both simulation and calibration sample size are evaluated by bootstrapping to create multiple samples, and repeating the procedure for each sample. The impact of potential mismodelling of the PID variables in the calibration samples is evaluated by describing the corresponding distributions using density estimates with different kernel widths. For each of these cases, alternative efficiency maps are produced to determine the associated uncertainties on the results of the analysis.
In principle, mismodelling of the proton and kaon reconstruction efficiencies, and associated asymmetries, could be a source of systematic uncertainty. However, such effects are known to be negligible at the level of precision achieved in this analysis [?, ?], and therefore are not accounted for explicitly.

Background shapes
The combinatorial background SDP distribution is obtained by extrapolating from an m(pK − K − ) sideband region, and has uncertainties related to the available yield in the sideband and the extrapolation procedure itself. The former is evaluated by bootstrapping to create multiple combinatorial background samples, and repeating the amplitude fit with each. The RMS of the distribution of the change in each result is taken as the systematic uncertainty. The latter is evaluated by changing the architecture of the neural network, with the change in each result with respect to its baseline value assigned as the associated systematic uncertainty.

Background asymmetry
In the baseline model, it is assumed that there is no local asymmetry in the combinatorial background as described in Sec. 5.2. The associated systematic uncertainty is evaluated by considering separate background distributions for Ξ − b and Ξ + b candidates. In order to obtain sufficiently large background samples to determine these separate distributions, the MVA output requirement for candidates in the sideband region is relaxed.
A possible global combinatorial background asymmetry is accounted for in the baseline fit, while cross-feed background is assumed to have no asymmetry. A fit allowing for a global cross-feed background asymmetry is performed, and the differences between the results in this fit and their nominal values is assigned as the systematic uncertainty arising from this assumption.

Production asymmetry
The baseline fit model assumes no asymmetry in the Ξ − b production rates in the LHCb acceptance, consistent with measurements [44]. To evaluate the associated systematic uncertainty, the model is adjusted to include production asymmetries within the experimentally allowed range by introducing a global asymmetry in the efficiency maps. An ensemble of fits with varied Ξ − b production asymmetries is performed, and the RMS of the distribution of the change in each result with respect to its baseline value is assigned as the systematic uncertainty.

Polarisation
The transverse polarisation of Ξ − b baryons produced in pp collisions is assumed to be consistent with zero, as observed for Λ 0 b baryons [49,50]. The distributions of m 2 low and m 2 high are independent of the Ξ − b polarisation. However, if the Ξ − b baryons produced in LHC collisions are polarised, efficiency variation across additional phase-space variables should be considered in the analysis. To evaluate the systematic uncertainty due to potential Ξ − b polarisation, two sets of pseudoexperiments are generated, with the Ξ − b polarisation set to 20% in one case and −8% in the other. This corresponds to the ±2σ range measured for the Λ 0 b baryon in Ref. [49,50], where σ indicates the Gaussian standard deviation. A conservatively broad range is taken to allow for differences between Ξ − b and Λ 0 b polarisation. The pseudoexperiments are generated using a signal model whose helicity couplings are set to the values from the nominal fit and where the measured efficiency variation over the additional phase-space observables is introduced. A fit to the Dalitz-plot variables is then performed using the baseline model. The largest deviation of the parameters from the nominal case is assigned as the systematic uncertainty.

Modelling of the lineshapes
Each resonant contribution has fixed parameters in the amplitude fit. These include masses, widths and Blatt-Weisskopf radius parameters. An ensemble of fits is obtained varying the masses and widths of all resonances within the range of values quoted by the PDG and given in Table 2. The Blatt-Weisskopf radius parameter associated with the Ξ − b baryon is varied in the range 3-7 GeV −1 and that associated with the resonances is varied in the range 0-3.5 GeV −1 . The RMS of the distribution of the change in each fitted parameter is taken as the systematic uncertainty.
For the Λ(1405) and Σ(1385) resonances that peak below the m low threshold, the effective RBW used in the baseline fit is replaced with a lineshape equivalent to the Flatté parameterisation as done in Ref. [46]. The total width is modified to account for the Σ + π − channel, i.e. Γ(m) = Γ pK − (m) + Γ Σ + π − (m), assuming equal couplings to both channels. For each of these lineshape variations, the differences in the results between fits with the alternative and baseline models are assigned as the associated systematic uncertainties.

Alternative fit model
The effect of including additional signal components in the fit model is examined to assign systematic uncertainty due to the composition of the baseline model. The Λ(1690), Λ(1820), Σ(1670) and NR( 3 2 + ) components are added to the nominal model individually.
These modifications of the model are chosen since they improve −2 ln L, although not by a significant amount. The largest deviation, among the four cases, from the nominal value of each measured quantity is taken as the associated systematic uncertainty.

Summary of systematic uncertainties
Ratio of fragmentation and branching fractions: Since separate fits are performed to the m(pK − K − ) distributions from the Run 1 and Run 2 samples, and the signal efficiencies are also determined separately, results for R in each of the two samples are obtained. Systematic uncertainties on R are considered as being either completely uncorrelated or 100% correlated between the two results. The systematic uncertainties that are uncorrelated between Run 1 and Run 2 are folded into their respective likelihood functions, by convolution with a Gaussian of appropriate width. The correlated systematic uncertainties are later folded into the combined likelihood that is obtained by multiplying the likelihood functions of the two samples.
The uncorrelated systematic uncertainties are those that are related to the fixed parameters in the fit model, the fit bias and the impact on the efficiency of the Ξ − b production kinematics, and the descriptions of the hardware trigger and particle identification response. The correlated systematic uncertainties are those related to knowledge of the phase-space distributions of the decays and the fit model choice. Slightly different procedures are used to obtain the total uncertainty for the two sources of correlated systematic uncertainty. That related to knowledge of the phase-space distribution affects the efficiency, and hence is a constant relative uncertainty. The method to evaluate the uncertainty due to fit model choice gives different relative uncertainties between Run 1 and Run 2. Since the two samples have approximately equal statistical weight in the combination, the average of the relative uncertainties is taken and assigned to the combined result. Table 3 summarises the systematic uncertainties on R.
Amplitude analysis: The results of the amplitude analysis are the CP asymmetry parameters A CP , defined in Eq. (25), the fit fractions F defined in Eq. (21) and the interference fit fractions I defined in Eq. (23). A summary of the systematic uncertainties on these quantities is shown in Table 4. The most precise results are those related to the Λ(1520) and Λ(1670) resonances. For Λ(1520), the dominant systematic uncertainty on A CP is due to ignoring the efficiency variation over angular variables if Ξ − b baryons are produced polarised, whereas for F it is due to the limited size of the sample used for combinatorial background modelling and the variation of the Blatt-Weisskopf radius parameters. For Λ(1670), the largest systematic uncertainty on A CP is due to variation of the Λ(1405) lineshape which has the same spin and parity as this component, whereas for F it is due to use of an alternate fit model. For Σ(1385), the dominant systematic uncertainty on A CP is due to use of an alternate fit model, whereas for F it is due to variation of the Blatt-Weisskopf radius parameters. For Λ(1405), the largest systematic uncertainty on A CP is due to use of an alternate fit model, whereas for F it is due to variation in its lineshape. For Σ(1775), the dominant systematic uncertainty on A CP is due to the limited size of the sample used for modelling of the combinatorial background, whereas for F it is due to use of an alternate fit model. For Σ(1915), the dominant systematic uncertainty on both A CP and F is due to use of an alternate fit model.
For interference fit fractions, the largest systematic uncertainties are mainly due to the use of an alternate fit model, the limited size of the sample used for modelling of the combinatorial background, variation of the resonance lineshapes and variation of Blatt-Weisskopf radius parameters.

Results
Ratio of fragmentation and branching fractions: The results for the ratio R of the relative fragmentation and branching fractions for Ω − b → pK − K − and Ξ − b → pK − K − decays are R = (−20 ± 30 (stat) ± 1 (uncorr syst)) × 10 −3 for Run 1 and R = ( 51 ± 28 (stat) ± 2 (uncorr syst)) × 10 −3 for Run 2 , where the first uncertainty is statistical and the second includes only the uncorrelated systematic effects presented in Table 3. The negative log-likelihood functions for these two results are added to obtain a combined result, where both uncorrelated and correlated systematic uncertainties are included. In the combined result, it is implied that , which may vary with centre-of-mass energy of the LHC pp collisions, is an effective value averaged over the Run 1 and Run 2 data samples. This result is found to be consistent with, and more precise than, the previous measurement [23]. No significant evidence of the Ω − b → pK − K − decay is found, and therefore an upper limit on R is calculated at 90 (95) % confidence level by integrating the likelihood in the physical region of non-negative branching fraction, Amplitude analysis: The results for the CP -asymmetry parameters for each component of the signal model are shown in Table 5. No significant CP asymmetry is observed. The fit fraction matrix is reported in Table 6. The diagonal elements correspond to the fit fractions of the respective components, and the off-diagonal elements are the interference fit fractions. These results are derived from the helicity couplings that are the free parameters of the amplitude fit. Their statistical uncertainties are evaluated from an ensemble of pseudoexperiments, while systematic uncertainties are obtained as described in Sec. 6. The significance of each component in the baseline model is evaluated using pseudoexperiments. These are generated, each with a sample size corresponding to the data, according to the best fit model with the component of interest removed from the model. They are then fitted both with the model used to generate and with the model including the component of interest. Twice the difference between the negative log-likelihood values obtained in these two fits (−2∆ ln L) is used as a test statistic. A p-value, corresponding to the probability of observing −2∆ ln L values as large or larger than that found in the fit to data, is found by extrapolating the tail of the distribution obtained from the ensemble of pseudoexperiments. In order to account for dominant systematic uncertainties, this procedure is performed for the alternative model that gives the smallest value of −2 ln L in fits to data. The outcome is that the Λ(1520) and Λ(1670) components have p-values corresponding to 12.0 σ and 6.1 σ, respectively. All other components have significance below 3.5 σ. Table 5: Results for the CP -asymmetry parameters. The statistical uncertainties are obtained from pseudoexperiments while the systematic uncertainties are obtained following the procedure described in Sec. 6.

Component
A The branching fraction of each quasi-two-body contribution to the Ξ − b → pK − K − decay, corresponding to an intermediate resonance R, can be obtained from its fit fraction F i , The branching fraction of Ξ − b → pK − K − has not been measured directly, but the ratio of fragmentation and branching fractions relative to the B − → K + K − K − decay is known [23]. This can be combined with the known values of B ( [65], assuming that f u = f d , to obtain where the dominant uncertainty is that due to possible SU(3)-breaking effects which affect 44]. Consequently, the values of the quasi-two-body branching fractions are found to be where the uncertainties are statistical, systematic and due to the knowledge of

Summary
The structure of Ξ − b → pK − K − decays has been studied through an amplitude analysis. This is the first amplitude analysis of any b-baryon decay mode allowing for CP -violation effects. The analysis uses pp collision data recorded with the LHCb detector, corresponding to integrated luminosities of 1 fb −1 at √ s = 7 TeV, 2 fb −1 at √ s = 8 TeV and 2 fb −1 at √ s = 13 TeV. Due to the inclusion of more data and significantly improving the selection procedure compared to the previous study of this channel [23], a yield of about 460 signal decays within the m(pK − K − ) signal region is obtained, with a signal to background ratio of about 2 : 1. A good description of the data is obtained with an amplitude model containing contributions from Σ(1385), Λ(1405), Λ(1520), Λ(1670), Σ(1775) and Σ(1915) resonances. The CP asymmetry for each contributing component is evaluated and no significant CP -violation effect is observed. The Ξ − b → Λ(1520)K − and Ξ − b → Λ(1670)K − decays are observed with significance greater than 5 σ, and their branching fractions measured, together with those of No significant signal for Ω − b → pK − K − decays is found and an upper limit on the ratio of fragmentation and branching fractions of With the substantially larger samples that are anticipated following the upgrade of LHCb [66,67], it will be possible to reduce both statistical and systematic uncertainties on CP -violation observables in three-body b-baryon decays, and thereby test the Standard Model using the methods pioneered in this study.     Figure 8: Distributions of (top) m low and (bottom) m high , with 1.6 < m low < 1.8 GeV, for (left) Ξ − b and (right) Ξ + b candidates, with results of the fits superimposed. The total fit result is shown as the blue solid curve, with contributions from individual signal components and from combinatorial (Comb) and cross-feed (Crsfd) background shown as indicated in the legend.     Figure 10: Distributions of (top) m low and (bottom) m high , with 2.2 < m low < 2.8 GeV, for (left) Ξ − b and (right) Ξ + b candidates, with results of the fits superimposed. The total fit result is shown as the blue solid curve, with contributions from individual signal components and from combinatorial (Comb) and cross-feed (Crsfd) background shown as indicated in the legend.  Entries / (0.10 GeV)  Figure 12: Distributions of (top) m low and (bottom) m high , with m high < 3.15 GeV , for (left) Ξ − b and (right) Ξ + b candidates, with results of the fits superimposed. The total fit result is shown as the blue solid curve, with contributions from individual signal components and from combinatorial (Comb) and cross-feed (Crsfd) background shown as indicated in the legend.   [9] LHCb collaboration, R. Aaij et al., Measurement of CP asymmetries in two-body B 0 (s) -meson decays to charged pions and kaons, Phys. Rev. D98 (2018) 032004, arXiv:1805.06759.