Amplitude analysis of the B+ -> pi+ pi+ pi- decay

The results of an amplitude analysis of the charmless three-body decay B+→ π+π+π−, in which CP -violation effects are taken into account, are reported. The analysis is based on a data sample corresponding to an integrated luminosity of 3 fb−1 of pp collisions recorded with the LHCb detector. The most challenging aspect of the analysis is the description of the behaviour of the π+π− S-wave contribution, which is achieved by using three complementary approaches based on the isobar model, the K-matrix formalism, and a quasi-model-independent procedure. Additional resonant contributions for all three methods are described using a common isobar model, and include the ρ(770)0, ω(782) and ρ(1450)0 resonances in the π+π− P-wave, the f2(1270) resonance in the π +π− D-wave, and the ρ3(1690) 0 resonance in the π+π− F-wave. Significant CP -violation effects are observed in both Sand D-waves, as well as in the interference between the Sand P-waves. The results from all three approaches agree and provide new insight into the dynamics and the origin of CP -violation effects in B+→ π+π+π− decays. Submitted to Phys. Rev. D c © 2019 CERN for the benefit of the LHCb collaboration. CC-BY-4.0 licence. †Authors are listed at the end of this paper. ar X iv :1 90 9. 05 21 2v 1 [ he pex ] 1 1 Se p 20 19


Introduction
In the Standard Model (SM), CP violation originates from a single irreducible complex phase in the Cabibbo-Kobayashi-Maskawa (CKM) matrix [1,2]. Thus far, all measurements of CP violation in particle decays are consistent with this explanation. Nevertheless, the degree of CP violation permitted in the SM is inconsistent with the macroscopic matter-antimatter asymmetry observed in the Universe [3], motivating further studies and searches for sources of CP violation beyond the SM.
For the manifestation of CP violation in decay, at least two interfering amplitudes with different strong and weak phases are required. In the SM, weak phases are associated with the complex elements of the CKM matrix and have opposite sign between charge-conjugate processes, while strong phases are associated with hadronic final-state effects and do not change sign under CP conjugation. In decays of b hadrons to charmless hadronic final states, contributions from both tree and loop (so-called "penguin") diagrams, which can provide the relative weak phase that is necessary for CP violation to manifest, are possible with comparable magnitudes. Indeed, significant CP asymmetries have been observed in both B 0 → K + π − [4][5][6][7] and B 0 → π + π − decays [4,8,9]. In multibody decays, variation across the phase space of strong phases, caused by hadronic resonances, allows for further enhancement of CP violation effects and a richer phenomenology compared to two-body decays. Large CP asymmetries localised in regions of phase space of charmless three-body B-meson decays have been observed in model-independent analyses [10][11][12], but until recently there has been no description of these effects with an accurate model of the contributing resonances. An amplitude analysis of B + → K + K + π − decays [13] has shown that ππ ↔ KK rescattering plays an important role in the observed CP violation, and it is anticipated that similar effects will occur in other charmless three-body B-meson decays.
This paper documents an analysis of the B + → π + π + π − decay amplitude in the twodimensional phase space known as the Dalitz plot [14,15]. The inclusion of charge-conjugate processes is implied, except where asymmetries are discussed. Previous studies of this decay mode indicate that the amplitude contains a sizable ρ(770) 0 component [16][17][18][19]. The amplitude analysis performed by the BaBar collaboration [19] additionally observed a large S-wave contribution, however, measurements of CP -violating quantities were limited by statistical precision.
The present analysis is performed on data corresponding to 3 fb −1 collected by the LHCb experiment, of which 1 fb −1 was collected in 2011 with a pp collision centre-of-mass energy of 7 TeV and 2 fb −1 was collected in 2012 with a centre-of-mass energy of 8 TeV. A model of the Dalitz plot distribution is constructed in terms of the intermediate resonant and nonresonant structures. Due to its magnitude and potential importance to the observed CP violation in this decay, particular attention is given to the π + π − S-wave contribution, which is known to consist of numerous overlapping resonances and open decay channels [34]. Three different state-of-the-art approaches to the modelling of the S-wave are used to ensure that any inaccuracies in the description of this part of the amplitude do not impact the interpretation of the physical quantities reported. This paper is organised as follows: Section 2 gives a brief description of the LHCb detector and the event reconstruction and simulation software; the signal candidate selection procedure is described in Section 3; Section 4 describes the procedure for estimating the signal and background yields that enter into the amplitude fit; Section 5 outlines the formalism used for the construction of the amplitude models, as well as a description of the mass lineshapes used to parameterise the intermediate structures; Section 6 describes the systematic uncertainties associated with the analysis procedure; Section 7 documents the physics parameters of interest obtained from the amplitude models and presents projections of the fit models on the selected data; these results are then discussed in Section 8; and a summary of the work as a whole can be found in Section 9. A shorter description of the analysis, more focussed on the first observations of different sources of CP -violation effects, can be found in a companion article [35].

Detector and simulation
The LHCb detector [36,37] 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 [38], 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 [39] placed downstream of the magnet. The tracking system provides a measurement of the momentum p of charged particles with relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary vertex (PV), or 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 (in GeV/c). Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [40]. 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 [41]. 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 almost cancels the effect. Approximately 60% of 2011 data and 52% of 2012 data was collected in the "down" polarity configuration, and the rest in the "up" configuration.
The online event selection is performed by a trigger [42] which consists of a hardware stage followed by a software stage. The hardware stage is based on information from the calorimeter and muon systems in which events are required to contain a muon with high p T , or a hadron, photon or electron with high transverse energy in the calorimeters. The software trigger requires a two-or three-track secondary vertex with significant displacement from all primary pp interaction vertices. All charged particles with p T > 500 (300) MeV/c, in which at least one charged particle must have transverse momentum p T > 1.7 (1.6) GeV/c and be inconsistent with originating from a PV, are reconstructed for data collected in 2011 (2012). A multivariate algorithm [43] is used for the identification of secondary vertices consistent with the decay of a b hadron.
Simulated data samples are used to investigate backgrounds from other b-hadron decays and also to study the detection and reconstruction efficiency of the signal. In the simulation, pp collisions are generated using Pythia [44] with a specific LHCb configuration [45]. Decays of unstable particles are described by EvtGen [46], in which final-state radiation is generated using Photos [47]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [48] as described in Ref. [49].

Selection
The selection of signal candidates follows closely the procedure used in the modelindependent analysis of the same data sample [12]. Signal B + candidates are formed from three good-quality tracks that originate from a secondary decay vertex (SV) with a fit χ 2 < 12. The SV must also be at least 3 mm away from any PV. The reconstructed B + candidate is associated with the PV that is most consistent with its flight direction. A requirement is also imposed on the angle between the B + momentum and the vector between the PV and SV, that must be less than approximately 6 mrad.
To reject random associations of tracks (combinatorial background), a boosted decisiontree classifier [50] is trained to discriminate between simulated signal candidates and candidates in collision data residing in a region where this background dominates, 5.4 < m(π + π + π − ) < 5.8 GeV/c 2 . The variables that enter this classifier are B + and decay product kinematic properties, quantities based on the quality of the reconstructed tracks and decay vertices, as well as the B + displacement from the PV. The requirement on the output of this classifier is optimised to maximise the expected approximate signal significance, N s / √ N s + N b , where N s is the expected signal yield within 40 MeV/c 2 of the known B + mass [51], and N b is the corresponding combinatorial background level within the same region.
To suppress backgrounds that arise when any number of kaons are misidentified as pions, requirements are placed on the particle-identification information associated with each final-state track. Furthermore, tracks associated with a hit in the muon system are also removed, as are tracks that are outside the fiducial region of the particle-identification system.
A veto in both combinations of the opposite-sign dipion mass, where 1.740 < m(π + π − ) < 1.894 GeV/c 2 , is applied to remove B + → (D 0 → π + π − )π + decays, along with partially reconstructed decays involving intermediate D 0 mesons and decays of D 0 mesons where one or more kaons are misidentified as pions. Approximately 2% of events contain multiple B + decay candidates following the aforementioned selection procedure, of which one is chosen at random.
The Dalitz-plot variables are calculated following a kinematic mass constraint, fixing the B + candidate mass to the known value to improve resolution and to ensure that all decays remain within the Dalitz-plot boundary. In consequence, the experimental resolution in the region with the narrowest resonance considered in this analysis, the ω(782) state, is better than 3 MeV/c 2 . This is smaller than the ω(782) width, and therefore effects related to the finite resolution in the Dalitz plot are not considered further.

B + candidate invariant-mass fit
An extended, unbinned, maximum-likelihood fit is performed to the m(π + π + π − ) invariantmass spectrum to extract yields and charge asymmetries of the B + → π + π + π − signal and various contributing backgrounds. The fit is performed to candidates in the range 5.080 < m(π + π + π − ) < 5.580 GeV/c 2 , and its results are used to obtain signal and background yields in the signal region, 5.249 < m(π + π + π − ) < 5.317 GeV/c 2 , in which the subsequent Dalitz-plot fit is performed. All shape parameters of the probability density functions (PDFs) comprising the fit model are shared between B + and B − candidates and only the yields are permitted to vary between these categories, which are fitted simultaneously. The data are also subdivided by data-taking year, and whether the hardware trigger decision is due to hadronic calorimeter deposits associated with the signal candidate, or due to other particles in the event.
The shape of the B + → π + π + π − signal decay is parameterised by the sum of a core Gaussian with two Crystal Ball functions [52] in order to describe the asymmetric tails of the distribution due to detector resolution and final-state radiation. The tail mean and width parameters of the Crystal Ball functions are determined from simulation relative to the core mean and width, which are left free in the fit to collision data to account for small differences between simulation and data. All remaining parameters, apart from the total yield, are obtained from a fit to simulated events.
Partially reconstructed backgrounds, which predominantly arise from four-body B-meson decays where a charged hadron or neutral particle is not reconstructed, are modelled with an ARGUS function [53] convolved with a Gaussian resolution function. The smooth combinatorial background is modelled with a falling exponential function. The only significant source of cross-feed background, where one or more kaons are misidentified as pions, is the B + → K + π + π − decay. To obtain an accurate model for this background, simulated B + → K + π + π − decays are weighted according to the amplitude model obtained by the BaBar collaboration [54], also accounting for the probability to be reconstructed as a B + → π + π + π − candidate. A model for B + → K + π + π − decays based on a similar-sized data sample has also been obtained by the Belle collaboration [55]; the details of the model used do not impact on the analysis. This shape is modelled in the invariant-mass fit as a sum of a Gaussian with two Crystal Ball functions, where all parameters are determined from a fit to the weighted simulation. Furthermore, the yield of this component is constrained to the B + → π + π + π − signal yield multiplied by the product of the relative branching fractions of these decays and the inverse of the relative overall reconstruction and selection efficiencies, which are described in Section 5.5.
The mass fit results are shown in Fig. 1, while Table 1 quotes the component yields and phase-space-integrated raw detection asymmetries in the B + signal region, which are subsequently used in the Dalitz-plot fit. The quoted uncertainties account for systematic effects evaluated with the procedures outlined in Section 6.   Parameter Value Signal yield 20 600 ± 1 600 Combinatorial background yield 4 400 ± 1 600 B + → K + π + π − background yield 143 ± 11 Combinatorial background asymmetry +0.005 ± 0.010 B + → K + π + π − background asymmetry 0.000 ± 0.008

Dalitz-plot model
The B + → π + 1 π + 2 π − 3 decay amplitude can be expressed fully in terms of the invariant-masssquared of two pairs of the decay products m 2 13 and m 2 23 . As no resonances are expected to decay to π + π + , the squared invariant-masses of the two combinations of oppositely charged pions are used as these two pairs. Due to Bose symmetry, the amplitude is invariant under exchange of the two like-sign pions, A(m 2 13 , m 2 23 ) ≡ A(m 2 23 , m 2 13 ), meaning that the assignments of π + 1 and π + 2 are arbitrary. Due to this symmetry, a natural "folding" occurs in the Dalitz plot about the axis m 2 13 = m 2 23 . Since the majority of the resonant structure is expected to be at low mass m(π + π − ) < 2 GeV/c 2 , the data and its projections are presented with the two axes being the squares of the low-mass m low and high-mass m high combinations of the opposite-sign pion pairs, for visualisation purposes. Plots of this kind are therefore similar to those found in other analyses with Dalitz plots that are expected to contain resonances along only one axis. In this case, structure resulting predominantly from the mass lineshape appears in m low , while m high is influenced by the angular momentum eigenfunctions. The Dalitz-plot distributions of the selected candidates can be seen in Fig. 2

Amplitude analysis formalism
In general, the isobar model is used to define the total amplitude for B + decays as a coherent sum over N components, each described by a function F j that parameterises the intermediate resonant or nonresonant processes: and similarly for B − : where the complex coefficients c j represent the relative contribution of component j. These are expressed in the "Cartesian" CP -violating convention, for c + j (c − j ) coefficients corresponding to B + (B − ) decays. The function F contains only strong dynamics and, for a resonant or nonresonant contribution to the m 13 spectrum, is parameterised as where R describes the mass lineshape, T describes the angular dependence, and X are Blatt-Weisskopf barrier factors [56] depending on a radius parameter r BW . Here and in the following, the momentum of one of the m 13 decay products is denoted as q and the momentum of the third pion (π 2 ) as p, where both momenta are evaluated in the rest frame of the dipion system. Using the Zemach tensor formalism [57,58], the angular probability distribution terms T ( p, q) are given by These are related to the Legendre polynomials P L (cos θ hel ), where the helicity angle θ hel is the angle between p and q and provide a good visual indicator of the spin of the intermediate state.
The Blatt-Weisskopf barrier factors account for the finite size of the decaying hadron and correct for the unphysical increase in the amplitude above the angular momentum barrier introduced by the form of angular momentum distributions given in Eq. (5). They are expressed in terms of z = | q |r P BW for the B + decay and z = | p |r R BW for the intermediate state, where L is the relative angular momentum between the B + meson and the resonance, which is equal to the spin of the resonance since the B + meson is spinless. The variable z 0 is equal to the value of z when the mass of the propagator is equal to the mass of the resonance. Unless otherwise stated, the barrier radius is taken to be r P BW = r R BW = 4.0 GeV −1 (≈ 0.8 fm) for all resonances. (To simplify expressions, natural units with c = = 1 are used throughout Sections 5.1-5.3.)

Square Dalitz plot
Since resonances tend to populate the edges of the conventional Dalitz plot in charmless B decays, it is useful to define the so-called "square" Dalitz plot [18], which provides improved resolution in these critical regions when using uniform binning, for example when modelling efficiency effects. Furthermore, the mapping to a square space aligns the bin boundaries to the kinematic boundaries of the phase space.

Mass lineshapes
Resonant contributions are mostly described by the relativistic Breit-Wigner lineshape with a mass-dependent decay width where q 0 is the value of q = | q | when the invariant mass, m, is equal to the pole mass, m 0 , of the resonance. Here, the nominal resonance width is given by Γ 0 . For the broad ρ(770) 0 resonance, an analytic dispersive term is included to ensure unitarity far from the pole mass, known as the Gounaris-Sakurai model [59]. It takes the form with an additional mass dependence where and dh dm 2 The constant parameter D is given by Isospin-violating ω(782) decays to two charged pions can occur via ω(782) mixing with the ρ(770) 0 state. To account for such effects, a model is constructed that directly parameterises the interference between these two contributions following Ref. [60], where R ρ (m) is the Gounaris-Sakurai ρ(770) 0 lineshape, R ω (m) is the relativistic Breit-Wigner ω(782) lineshape, |ζ| and φ ζ are free parameters of the fit that denote the respective magnitude and phase of the production amplitude of ω(782) with respect to that for the ρ(770) 0 state, and ∆ ≡ δ (m ρ + m ω ), where δ governs the electromagnetic mixing of ρ(770) 0 and ω(782) and m ρ and m ω represent the known particle masses [51]. The value for δ is fixed in the fit to δ = 0.00215 ± 0.00035 GeV [60]. This is equivalent to the parameterisation described in Ref. [61] if the small ∆ 2 term in the denominator is ignored.

S-wave models
The S-wave (L = 0) component of the B + → π + π + π − amplitude is both large in magnitude and contains many overlapping resonances and decay channel thresholds, i.e., where increasing two-body invariant mass opens additional decay channels that were previously inaccessible, thereby modulating the observed intensity. This analysis includes three distinct treatments of the S-wave component in B + → π + π + π − in an attempt to better understand its behaviour. They also increase confidence that parameters reported for the non-S-wave contributions are robust and provide additional information for further study. As such, three sets of results are presented here, corresponding to the cases where the π + π − S-wave is described by: (i) a coherent sum of specific resonant contributions (isobar); (ii) a monolithic, unitarity-conserving model informed by historical scattering data (K-matrix); and (iii) a quasi-model-independent binned approach (QMI). All approaches contain identical contributions to higher partial waves, where L > 0.

Isobar model
The isobar model S-wave amplitude is represented by the coherent sum of contributions from the σ, or f 0 (500), meson and a ππ ↔ KK rescattering amplitude within the mass range 1.0 < m(π + π − ) < 1.5 GeV. The σ meson is represented as a simple pole [62], parameterised as with pole mass, m σ , and width, Γ σ , extracted from the fit. The concept of ππ ↔ KK rescattering was originally developed inside the context of two-body interactions. For three-body decays, rescattering means that a pair of mesons produced in one channel will appear in the final state of a coupled channel. Therefore, a model is used that describes the source of the rescattering [32,33], where ∆ 2 ππ = ∆ 2 KK = 1 GeV. The total rescattering amplitude in the three-body B + decay is then The amplitude f rescatt (m) = 1 − η(m) 2 e 2iδ(m) is described in terms of the inelasticity, η(m), and a phase shift, δ(m). Functional forms of these are used that combine constraints from unitarity and analiticity with dispersion relation techniques [63]. The inelasticity is described by with where m K = 0.494 GeV is the charged kaon mass, 1 = 2.4, 2 = −5.5, and M = 1.5 GeV.
The phase shift is given by where M f = 1.32 GeV, c 0 = 1.3 and M s = 0.92 GeV. Except for m K , values of all parameters are taken from Ref. [63].

K-matrix model
The coherent sum of resonant contributions modelled with Breit-Wigner lineshapes can be used to describe the dynamics of three-body decays when the quasi-two-body resonances are relatively narrow and isolated. However, when there are broad, overlapping resonances (with the same isospin and spin-parity quantum numbers) or structures that are near open decay channels, this model does not satisfy S-matrix unitarity, thereby violating the conservation of quantum mechanical probability current.
Assuming that the dynamics is dominated by two-body processes (i.e. that the S-wave does not interact with other decay products in the final state), then two-body unitarity is naturally conserved within the K-matrix approach [64]. This approach was originally developed for two-body scattering [65] and the study of resonances in nuclear reactions [66,67], but was extended to describe resonance production and n-body decays in a more general way [68]. This section provides a brief introduction to the K-matrix approach as applied in this analysis; for more detail see Ref. [69].
From unitarity conservation, the form-factor for two-body production is related to the scattering amplitude for the same channel, when including all coupled channels. In this way, the K-matrix model describes the amplitude, F u , of a channel u in terms of the initialP -vector preparation of channel states v, that has the same form asK, "scattering" into the final state u via the propagator term (Î − iKρ) −1 , The diagonal matrixρ accounts for phase space, where the element for the two-body channel u is given by [51] where s = m 2 π + π − and m 1 and m 2 are the rest masses of the two decay products. This expression is analytically continued by settingρ uu to be i|ρ uu | when the channel is below its mass threshold (provided it does not cross into another channel). For the coupled multi-meson final states, the corresponding expression can be found in Ref. [69].
The scattering matrix,K, can be parameterised as a combination of the sum of N p poles with real bare masses m α , together with nonresonant "slowly-varying" parts (SVPs). These slowly-varying parts are so-called as they have a 1/s dependence, and incorporate real coupling constants, f scatt uv [70]. The scattering matrix is symmetric in u and v, where g α u and g α v denote the real coupling constants of the pole m α to the channels u and v, respectively. The factor is the Adler zero term, which suppresses the false kinematic singularity due to left-hand cuts when s goes below the ππ production threshold [71]. The parameters m 2 0 , s scatt 0 and s A0 are real constants. Note that the masses m α are those of the poles, or the so-called bare states of the system, which do not correspond to the masses and widths of resonances -mixtures of bare states.
Extension to three-body decays is achieved by fitting for the complex coefficients β α and f prod v of the production pole and SVP terms in theP -vector, where these coefficients are different for B + and B − decays to allow for CP violation. The parameter s prod 0 is dependent on the production environment and is taken from Ref. [72]. Using the above expressions, the amplitude for each production pole α to the ππ channel (denoted by the subscript u = 1) is given by where the sum is over the intermediate channels, v, while the SVP production amplitudes are separated out for each individual channel as All of these contributions are then summed to give the total S-wave amplitude TheK matrix elements in Eq. (24) are completely defined using the values quoted in Ref. [73] and given in Table 2 (for five channels n = 5 and five poles N p = 5), which are obtained from a global analysis of ππ scattering data [70]. Table 2: K-matrix parameters quoted in Ref. [73], which are obtained from a global analysis of ππ scattering data [70]. Only f 1v parameters are listed here, since only the dipion final state is relevant to the analysis. Masses m α and couplings g α u are given in GeV, while units of GeV 2 for s-related quantities are implied; s prod 0 is taken from Ref. [72].

Quasi-model-independent analysis
In the quasi-model-independent (QMI) analysis, the amplitude for the ππ S-wave is described by individual magnitudes and phases within each of 17 bins in m(π + π − ). The QMI approach exploits the fact that the S-wave amplitude is constant in cos θ hel (Eq. (5)) to disentangle this component from other contributions to the Dalitz plot, assuming the higher-order waves to be well modelled by the isobar approach. The bins are defined identically for B + and B − : 13 below the charm veto and 4 above. The binning scheme is chosen ad hoc as optimisation requires a priori knowledge of the S-wave and total B + → π + π + π − amplitude model. The bin boundaries are chosen adaptively, by requiring initially roughly equal numbers of candidates in each region and using the isobar model subsequently to tune the bins in order to reduce intrinsic bias in the method's ability to reproduce a known S-wave.
Similar QMI approaches have previously been performed in amplitude analyses of various decays of B and D mesons to study the Kπ [74][75][76][77], ππ [78] and Dπ [79] S-waves. In contrast to these previous approaches, no interpolation between bins is performed in this analysis -the amplitudes are constant within each bin. This choice is made as interpolation is based on a premise of smoothness, which is appropriate for goals such as confirmation of a new resonance. However, interpolation is not optimal in the description of the ππ S-wave due to the opening of various decay channels that become kinematically allowed as the mass increases, and which cause sharp changes in the amplitude on scales less than the bin width. A key difference between this and the other two S-wave approaches is that the QMI absorbs the average contribution from potential interactions with the other decay products in the final state. These final state interactions may be quantified by comparing with similar dipion S-wave distributions obtained by rescattering experiments.

Measurement quantities
The primary outputs of the Dalitz-plot fit are the complex isobar coefficients c ± j defined in Eq. (1). However, since these depend on the choice of phase convention, amplitude formalism and normalisation, they are not straightforwardly comparable between analyses and have limited physical meaning. Instead, it is useful to compare fit fractions, defined as the integral of the absolute value of the amplitude squared for each intermediate component, j, divided by that of the coherent matrix-element squared for all intermediate contributions, These fit fractions will not sum to unity if there is net constructive or destructive interference. The interference fit fractions are given by where the sum of the fit fractions and interference fit fractions must be unity. The fit fractions are defined for the B + and B − decay amplitudes separately, however a measure of the CP -averaged contribution can be obtained from the CP -conserving fit fraction Since the ρ(770) 0 and ω(782) components cannot be completely decoupled in the mixing amplitude, their effective lineshapes are defined for the purpose of calculating fit fractions. This is achieved by separating Eq. (15) into the two respective terms with a common denominator.
Another important physical quantity is the quasi-two-body parameter of CP violation in decay associated with a particular intermediate contribution The CP asymmetry associated with the S-wave cannot be determined using Eq. (33) since this component involves several contributions. Instead, it is determined as the asymmetry of the relevant B − and B + decay rates, where A ± S is the coherent sum of all contributions to the S-wave.

Efficiency model
The efficiency of selecting a signal decay is parameterised in the two-dimensional square Dalitz plot and determined separately for B + and B − decays. Non-uniformities in these distributions arise as a result of the detector geometry, reconstruction and trigger algorithms, particle identification selections and other background rejection requirements such as that imposed on the boosted decision-tree classifier to discriminate against combinatorial background. The efficiency map is primarily obtained using simulated decays, however effects arising from the hardware trigger and particle identification efficiency are determined using data calibration samples. The hardware-trigger efficiency correction is calculated using pions from D 0 → K − π + decays, arising from promptly produced D * + → D 0 (K − π + )π + decays, and affects two disjoint subsets of the selected candidates: those where the trigger requirements were satisfied by hadronic calorimeter deposits as a result of the signal decay and those where the requirements were satisfied only by deposits from the rest of the event. In the first case, the probability to satisfy the trigger requirements is calculated using calibration data as a function of the transverse energy of each final-state particle of a given species, the dipole-magnet polarity, and the hadronic calorimeter region. In the second subset, a smaller correction is applied following the same procedure in order to account for the requirement that these tracks did not fire the hardware trigger. These corrections are combined according to the relative abundance of each category in data.
The particle identification efficiency is calculated from calibration data also corresponding to the D * + → D 0 (K − π + )π + decay, where pions and kaons can be identified without the use of the LHCb particle identification system [80]. The particle identification efficiencies for the background-subtracted pions and kaons are parameterised in terms of their total and transverse momentum, and the number of tracks in the event. This efficiency is assumed to factorise with respect to the final-state tracks and therefore the efficiency for each track is multiplied to form the overall efficiency.
The effect of an asymmetry between the production rates of B − and B + mesons is indistinguishable in this analysis from a global detection efficiency asymmetry. Therefore, the B + production asymmetry, as measured within the LHCb acceptance [81], is taken into account by introducing a global asymmetry of approximately −0.6% into the efficiency maps. This is obtained as an average of the measured production asymmetries, weighted by the relative integrated luminosity obtained in 2011 and 2012. The overall efficiency, as a function of square Dalitz-plot position, can be seen in Fig. 3 for B + and B − decays separately. These histograms are smoothed by a two-dimensional cubic spline to mitigate effects of discontinuity at the bin edges. The signal PDF for B + or B − decays is then given by where ± represents the Dalitz-plot dependent efficiency for the B ± decay.

Background model
The dominant source of background in the signal region is combinatorial in nature. In the Dalitz-plot fit, the distribution of this background is modelled separately for B + and B − decays using square Dalitz-plot histograms of upper sideband data, from the region 5.35 < m(π + π + π − ) < 5.68 GeV/c 2 , with a uniform 16 × 16 binning in m and θ . In this region, a feature is observed that can be identified as arising from real B 0 → π + π − decays combined with a random track from the rest of the event. However, this background does not enter the signal region due to kinematics. As such, this feature is modelled in the m(π + π − ) spectrum using a Gaussian function located at the known B 0 mass, and events are subtracted from the combinatorial background histograms accordingly. The corresponding combinatorial background distributions can be seen in Fig. 4. For use in the fit, these histograms are smoothed using a two-dimensional cubic spline to mitigate effects of discontinuity at the bin edges. In the Dalitz-plot fit, the charge asymmetry in the combinatorial background yield is fixed to that obtained in the B + invariant-mass fit described in Section 4.
A secondary source of background arises from misidentified B + → K + π + π − decays. This background is modelled using simulated B + → K + π + π − events, reconstructed under the B + → π + π + π − decay hypothesis. To account for the phase-space distribution of this background, the events are weighted according to the amplitude model obtained by the BaBar collaboration [54]. Similarly to the combinatorial background, this contribution is described in terms of a uniform 16×16 binned square Dalitz-plot histogram, smoothed with a two-dimensional cubic spline. The corresponding distribution, without the smoothing, can be seen in Fig Yield (candidates) Figure 5: Square Dalitz-plot distribution for the misidentified B + → K + π + π − background model, scaled to represent its yield in the signal region.

Fit procedure
Of the three approaches to the S-wave, the isobar and K-matrix fits are performed using the Laura++ Dalitz-plot fitter package v3.5 [69], which interfaces to the MINUIT function minimisation algorithm [82,83]. In contrast, the QMI approach relies on the Mint2 [84] amplitude-analysis interface to Minuit2 [83]. The fundamental difference between these amplitude-analysis software packages is in the handling of the normalisation. The former approximates the definite integral by employing a Gaussian quadrature approach, while the latter invokes a Monte Carlo technique. Additionally, due to the size of its parameter space, the QMI greatly benefits from the use of GPU-accelerated solutions. In all cases, the combined ρ-ω mixing component is set to be the "reference" amplitude. In practice, this means that the average magnitude of the B + and B − coefficients for this component is set to unity (in terms of Eq. (3), x = 1), while the δx parameter is left free to vary to allow for CP violation. Since there is no sensitivity to the phase difference between the B + and B − amplitudes, the imaginary part of the ρ-ω component is set to zero for both B + and B − (y = δy = 0), which means that all other contributions to the model are measured relative to this component.
The extended likelihood function that is optimised is of the form where N k is the yield for the candidate category k (given in Table 1), N is equal to k N k , N cand is the total number of candidates, and P i k is the probability density function for candidates in category k in terms of the Dalitz-plot coordinates. The optimal values of the fitted parameters are found by minimising twice the negative log-likelihood, −2 log L.
Since Dalitz-plot analyses involve multidimensional parameter spaces, depending upon the initial parameter values the results may correspond to a local, rather than global, minimum of the −2 log L function. To attempt to find the global minimum, a large number of fits are performed where the initial values of the complex isobar coefficients c j are randomised. The fit result with the smallest −2 log L value out of this ensemble is then taken to be the nominal result for each S-wave method.

Model selection
The inclusion or exclusion of various resonant contributions to the amplitude is studied using the isobar and K-matrix S-wave approaches. This is not practical with the QMI approach as the large S-wave parameter space requires a detailed search for the global minimum given each model hypothesis. Starting with resonant contributions identified during previous analyses of the B + → π + π + π − decay [16][17][18][19], additional resonances are examined iteratively, in the order that maximises the change in log-likelihood between the current and proposed model with respect to the data. This procedure is terminated when the log-likelihood gain from including any contribution not yet in the model is less than 10. Only resonances that have been observed by two or more experiments and have been seen to decay to two charged pions are considered initially. Scalar and vector nonresonant contributions, and possible virtual excited B * states, are then investigated as possible improvements to the model, however none are found to have a significant contribution.
After this initial iterative procedure, a second step is performed that involves ad hoc trials of alternative mass lineshapes for resonances already included in the model, and investigation of additional, more speculative, states. These form the basis of several important systematic uncertainties listed in Section 6 and are further discussed in Section 8. Lastly, tests are performed for "latent" resonant contributions up to spin 4, where a resonance is inserted as a relativistic Breit-Wigner shape with a width of 0.025 GeV, 0.050 GeV, 0.100 GeV, or 0.150 GeV, in mass steps of 0.2 GeV. No significant evidence of any resonant structure not captured by the previously established model is observed.
The goodness of fit is assessed by comparing the fit model with the data in square Dalitz-plot bins and determining an associated χ 2 value. The binning is chosen through an adaptive procedure that requires an approximately constant number of candidates from the data sample in each bin. For various values of the required number of candidates per bin, the ratio of the χ 2 to the number of bins is approximately 1.5 accounting for statistical uncertainties only. Given the impact of the systematic uncertainties on the results, as shown in Section 7, the agreement of the fit models with the data is reasonable. Smaller χ 2 values are obtained for the S-wave models with larger numbers of free parameters, such that all three approaches have comparable goodness-of-fit overall. The distribution in the square Dalitz plot of bins that contribute significantly to the χ 2 does not reveal any clear source of mismodelling.
Resonant contributions with spin greater than zero that were identified through the model selection procedure are common to all three approaches and are listed in Table 3. Furthermore, the mass and width of the dominant ρ(770) 0 contribution are left free to vary, which results in a significantly better fit quality. Table 3: Non-S-wave resonances and their default lineshapes as identified by the model selection procedure. These are common to all S-wave approaches.

Systematic uncertainties
Sources of systematic uncertainty are separated into two categories: those that arise from experimental effects and those from the inherent lack of knowledge on the amplitude models. The experimental systematic uncertainties comprise those from the uncertainty on the signal and background yields, the phase-space-dependent efficiency description, the combinatorial and B + → K + π + π − background models in the Dalitz plot, and the intrinsic fit bias. Model systematic uncertainties comprise those introduced by the uncertainty on the known resonance masses and widths, the radius parameter of the ad hoc Blatt-Weisskopf barrier factors and from potential additional resonant contributions to the amplitude. Furthermore, this latter category also includes sources of uncertainty that are specific to each S-wave approach. The effects in each category are considered to be uncorrelated and are therefore combined in quadrature to obtain the total systematic uncertainty.
The uncertainties on the signal yield and the background yields and asymmetries, given in Table 1, comprise a statistical component as well as systematic effects due to the invariant-mass fit procedure. The uncertainty arising from assumptions regarding the signal parameterisation is found by replacing the model with two Crystal Ball functions with a common mean and width, but independent tail parameters. Similarly, the model for the combinatorial background is replaced with a first-order polynomial. The uncertainty on the cross-feed B + → K + π + π − background shape in the three-body invariant-mass fit is negligible, however the yield of this component is varied by three times the nominal uncertainty on the expectation from the simulation to account for possible inaccuracies in the constraint. Additionally, effects associated with allowing different relative signal and partially reconstructed background yields in the data subcategories separated by source of hardware trigger decision are investigated by constraining them to be common in both subcategories. The combined statistical and systematic uncertainties on the signal and background yields and the background asymmetries are then propagated to the Dalitz-plot fit, where those variations causing the largest upward and downward deviations with respect to the nominal yield values are taken to assign the systematic uncertainty relating to the three-body invariant-mass fit.
To account for the statistical uncertainty on the efficiency description, an ensemble of efficiency maps is created by sampling bin-by-bin from the baseline description, according to uncorrelated Gaussian distributions with means corresponding to the central value of the nominal efficiency in each bin, and widths corresponding to the uncertainty. The standard deviation of the distribution of resulting Dalitz-plot fit parameters obtained when using this ensemble is then taken to be the associated systematic uncertainty. To account for potential biases in the method used to correct the hardware trigger efficiency, an alternative method using B 0 → J/ψ (µ + µ − )K + π − decays, requiring a positive trigger decision on the muons from the J/ψ decay, is used to apply corrections to the simulation [85,86]. The effect on the baseline results is assigned as an uncertainty.
Additionally, to account for potential variation of the efficiency within a nominal square Dalitz-plot bin, the efficiency map is constructed using a finer binning scheme, and the total deviation of the results is taken as the systematic uncertainty. The effect arising from the uncertainty on the measured B + production asymmetry is also considered, but is found to be negligible.
The statistical uncertainty on the combinatorial background distribution is propagated to the Dalitz-plot fit results in a procedure similar to that for the efficiency map. Uncertainty associated with the Dalitz-plot model of the B + → K + π + π − decay is also assigned. This is calculated by fluctuating the parameters obtained in the B + → K + π + π − fit according to their uncertainties [54], taking into account the reported correlations on the statistical uncertainties. The standard deviation in the variation of the subsequent Dalitz-plot fit results is taken to be the systematic uncertainty due to this effect.
Systematic uncertainties related to possible intrinsic fit bias are investigated using an ensemble of pseudoexperiments. Differences between the input and fitted values from the ensemble for the fit parameters are generally found to be small. Systematic uncertainties are assigned as the sum in quadrature of the difference between the input and output values and the uncertainty on the mean of the output value determined from a fit to the ensemble.
Sources of model uncertainty independent of the S-wave approach are those arising from the uncertainties on the masses and widths of resonances in the baseline model, the Blatt-Weisskopf barrier factors and contributions from additional resonances. The systematic uncertainty due to resonance masses and widths are again assigned with an ensemble technique, where the parameter values, excluding those that appear in the isobar S-wave model, are fluctuated according to the uncertainties listed in the Particle Data Group tables [51]. Where appropriate, these are taken to be those from combinations only considering decays to π + π − . The uncertainty arising from the lack of knowledge of the radius parameter of the Blatt-Weisskopf barrier factors is estimated by modifying the value of this between 3 and 5 GeV −1 , with the maximum deviation of the fit parameters taken to be the systematic uncertainty.
To account for mismodelling in the f 2 (1270) region, discussed in Section 8, an additional systematic uncertainty is assigned as the maximal variation in fit parameters when either an additional spin-2 component with mass and width parameters determined by the fit is included into the model, or when the f 2 (1270) resonance mass and width are permitted to vary in the fit. Furthermore, a possible contribution from the ρ(1700) 0 resonance cannot be excluded. Using perturbative QCD calculations, the branching fraction of the B + → ρ(1700) 0 π + decay has been calculated to be around 3 × 10 −7 [87], which is plausibly within the sensitivity of this analysis. Therefore a systematic uncertainty is assigned as the deviation of the fit parameters with respect to the nominal values, when the ρ(1700) 0 contribution is included.
For fits related to the isobar approach, the nominal rescattering parametrisation relies on a source term with two components as given in Eq. (17). The fits have little sensitivity to the values chosen for the ∆ 2 ππ and ∆ 2 KK parameters, so the robustness of this parametrisation is investigated by using instead a source term with only one component, A source = [1 + (m/∆ 2 KK )] −1 , and the difference in the results obtained assigned as a systematic uncertainty. The parameters of the σ contribution to the S-wave are also varied within the uncertainties on the world-average mass and width, and the effect on the results taken as a systematic uncertainty.
For fits using the K-matrix approach, both the fourthP -vector pole and the fourth slowly varying part result in a negligible change to the total likelihood when removed, and therefore a systematic uncertainty is assigned that corresponds to the maximum deviation of the parameters, with respect to the nominal values, when these components are removed from the K-matrix model. Furthermore, in the baseline fit the s prod 0 parameter appearing in the slowly varying parts of Eq. (28) is fixed to a value of −3 GeV 2 /c 4 . However as this comprises part of the production component of the K-matrix, this is not fixed by scattering data and can depend on the production environment. As such, this value is varied between −1 and −5 GeV 2 /c 4 based on the likelihood profile and the maximum deviation from the nominal fit results taken to be the systematic uncertainty due to this effect.
For the fits involving the QMI approach, an additional bias may arise from the intrinsic ability of the approach to reproduce the underlying analytic S-wave. Causes of such a bias can include the definition of the binning scheme, the extent to which the S-wave interferes with other partial waves in a particular bin, and the approximation of an analytic lineshape by a constant amplitude in each bin. This systematic uncertainty is evaluated reusing the ensemble of pseudoexperiments generated for estimating the K-matrix fit bias, fitting them with the QMI model, and determining the difference between the obtained and true bin-averaged values of the S-wave amplitude. The QMI intrinsic bias is by far the dominant systematic uncertainty on the S-wave magnitude and phase motion. Previous quasi-model-independent partial wave analyses have not recognised such an effect as a possible source of bias; an important conclusion of this study is that the associated systematic uncertainty must be accounted for in analyses in which quantitative results from binned partial-wave amplitude models are obtained.
The systematic uncertainties for the CP -averaged fit fractions and quasi-two-body CP asymmetries are summarised in Tables 4 and 5 for the isobar approach, Tables 6  and 7 for the K-matrix approach, and Tables 8 and 9 for the QMI approach. In general the largest sources of systematic uncertainty are due to variations in the model, which tend to dominate the total uncertainties for the CP -averaged fit fractions while the CP asymmetries for well established resonances are somewhat more robust against these effects. In particular, the inclusion of an additional tensor or vector resonance, i.e. the f 2 (1430) or ρ(1700) 0 states, can have a large effect on parameters associated with other resonances, particularly when they are in the same partial wave. With larger data samples it may be possible to clarify the contributions from these amplitudes and thereby reduce these uncertainties. Intrinsic fit bias is also an important source of uncertainty for several measurements, in particular those using the QMI description of the S-wave.

Results
Numerical results for the fit fractions and quasi-two-body CP asymmetries are given in this section; Appendices A-C also provide the correlation matrices for the CP -averaged fit fractions and CP asymmetries, while Appendix D presents the results in a way that may be more convenient for some purposes. The complex coefficients are fitted in terms of Cartesian parameters as shown in Eq. (3), but it can also be convenient to interpret them in terms of magnitudes and phases. A comparison of the phases of the non-S-wave contributions between the three approaches can be found in Appendix E. The fitted complex coefficients are recorded for completeness in Appendices A, B and C for the isobar, K-matrix and QMI approaches, respectively. In general, the statistical uncertainty is lowest for the model with the fewest parameters (isobar), and highest for the model with the largest number of parameters (QMI), as expected.
A comparison of the data and all three fit models, projected onto m low for a large dipion mass range, along with the asymmetries between B − and B + decays, can be seen in Fig. 6. In these and subsequent figures, the difference between the data and model expectation, divided by the uncertainty on this quantity (the "pull") is also shown below in the same binning scheme. Projections focussing on the low m low and the ρ(770) 0 regions are shown in Fig. 7, while the f 2 (1270) and high m high regions are displayed in Fig. 8. The fit result projected on the helicity angle in each of the ρ(770) 0 and f 2 (1270) regions, shown in Fig. 9, is also given separately above and below the ρ(770) 0 pole in Fig. 10. The projection on the helicity angle in the vicinity of the ρ 3 (1690) 0 resonance is shown in Fig. 11. Figure 12 shows the raw difference in the number of B − and B + candidates in the low m low region for negative and positive helicity angle cosines. Additional projections separating the contributions for various components of the amplitude model are shown in Appendices F, G and H.

Fit fractions
The CP -averaged fit fractions are given in Table 10 for all three S-wave approaches. The fit fractions and interference fit fractions, separated by B ± charge for each S-wave approach, are given in Tables 11-16. In all cases, statistical uncertainties are calculated using 68% confidence intervals obtained from the results of fits performed to data sets sampled from the nominal fit models. Throughout this paper, if three uncertainties are listed, they are separated into statistical, systematic and amplitude model sources, whereas if only two are listed, the systematic and models sources have been combined in quadrature for brevity. The total sums of fit fractions for the B + amplitudes are found to be (93.7 ± 2.6 ± 1.5 ± 4.5)% for the isobar approach, (99.2 ± 1.8 ± 4.1 ± 5.7)% for K-matrix and (92.2 ± 1.2 ± 7.7 ± 3.2)% for QMI. The corresponding quantities for the B − amplitudes are (100.7 ± 2.7 ± 1.7 ± 6.0)%, (108.3 ± 1.7 ± 3.3 ± 9.3)% and (108.0 ± 1.7 ± 3.7 ± 6.3)%.

CP asymmetries
Quasi-two-body CP asymmetries associated to each component are shown in Table 17. A detailed discussion of these results is given in Section 8, however it should be stressed that CP -violation effects can manifest in Dalitz-plot distributions through interference effects that leave values of the quasi-two-body CP asymmetries consistent with zero, and  indeed this occurs in B + → π + π + π − decays. The CP asymmetry integrated across the Dalitz plot is consistent, in all three models, with the value previously determined through model-independent analysis [12].

S-wave projections
The squared amplitude and phase motion of the S-wave models as a function of m(π + π − ) can be seen in Fig. 13(a) and (b) for the isobar approach, Fig. 13(c) and (d) for the K-matrix approach and Fig. 13(e) and (f) for the QMI approach. A comparison of all three models, for the CP -averaged S-wave projections, can be seen in Fig. 14. The QMI S-wave is recorded in Table 18, while the statistical and systematic correlation matrices

Multiple fit solutions
A search for secondary solutions with negative log-likelihood values worse than, but close to, that of the best fit is performed for each S-wave approach by setting the initial values of the complex coefficients of the model to random values and repeating the fit to data. In both the isobar and QMI approaches, no secondary solutions are found within 25 units of −2 log L. For the K-matrix approach, however, secondary solutions are found in which some of the pole or SVP amplitude coefficients are rotated in the Argand plane with respect to the best fit result. Studies using data sampled from the nominal model indicate that these could potentially be resolved with larger data samples, and further improvements may also be possible by fitting for the scattering parameters along with the amplitude coefficients. Isobar coefficients and S-wave projections corresponding to the    Table 11: Fit (diagonal) and interference (off-diagonal) fractions for B + decay in units of 10 −2 , between amplitude components in the isobar approach. The first uncertainty is statistical and the second the quadratic sum of systematic and model sources.     Table 13: Fit (diagonal) and interference (off-diagonal) fractions for B + decay in units of 10 −2 , between amplitude components in the K-matrix approach. The first uncertainty is statistical and the second the quadratic sum of systematic and model sources.
S-wave 23.1 ± 0.9 ± 3.7 Table 14: Fit (diagonal) and interference (off-diagonal) fractions for B − decay in units of 10 −2 , between amplitude components in the K-matrix approach. The first uncertainty is statistical and the second the quadratic sum of systematic and model sources.
S-wave 28.1 ± 0.7 ± 3.1 Table 15: Fit (diagonal) and interference (off-diagonal) fractions for B + decay in units of 10 −2 , between amplitude components in the QMI approach. The first uncertainty is statistical and the second the quadratic sum of systematic and model sources. Table 16: Fit (diagonal) and interference (off-diagonal) fractions for B − decay in units of 10 −2 , between amplitude components in the QMI approach. The first uncertainty is statistical and the second the quadratic sum of systematic and model sources.   LHCb QMI (e) (f)  The blue curve indicates the isobar S-wave, the amber curve indicates the K-matrix S-wave, and the green points with error bars represent the QMI S-wave. The band or error bars in each case represent the total uncertainty, incorporating the dominant systematic uncertainties. As the integral of the |A| 2 plot in each approach is proportional to its respective S-wave fit fraction, the overall scale of the K-matrix and QMI plots are set relative to the isobar S-wave fit fraction in order to facilitate comparison between the three approaches.  3.50 ≤ m(π + π − ) < 5.14 0.00 ± 0.01 ± 0.03 +100 ± 72 ± 173 0.02 ± 0.01 ± 0.07 −144 ± 15 ± 121

Discussion
The results and figures presented in Section 7 show that the three models exhibit good overall agreement with the data and with each other, both in CP -average projections and in the variation of the asymmetries across the phase space. In this section the main features observed in the data and in the models are discussed in more detail. Many of the interference fit fractions are close to zero, as expected, since interference effects between partial waves with even and odd values of the relative angular momentum cancel when integrated over the helicity angle. The largest interference fit fraction is between the combined ρ-ω component and the ρ(1450) 0 resonance; since each of these is spin-1, the interference does not vanish when integrated over the Dalitz plot. No significant interference fit-fraction asymmetries are observed, however this does not preclude sizeable asymmetries in localised regions of the Dalitz plot.

The ρ(770) 0 -ω(782) region
The interference between the spin-1 ρ(770) 0 and ω(782) resonances is well described by the models, as shown in Fig. 7(b). No significant asymmetry is observed in this region when integrating over cos θ hel as shown in Fig. 7(d), and also seen in the ρ(770) 0 and ω(782) quasi-two-body CP asymmetry parameters in Table 17. A number of theoretical calculations of these quantities are available in the literature, with some authors [88][89][90][91] predicting values for A CP (B + → ρ(770) 0 π + ) that are consistent with the measured result, albeit sometimes with large uncertainties. Other approaches [92][93][94][95] give predictions for this quantity which appear to now be ruled out. There is also no evident CP -violation effect associated with ρ-ω mixing, contrary to some theoretical predictions [25][26][27].
A significant CP asymmetry in the ρ(770) 0 region can, however, be seen in the cos θ hel projections shown in Figs. 10(c) and (d) when dipion masses below and above the known ρ(770) 0 mass are inspected separately. The same effect can be seen in Fig. 12 where the data are separated by the sign of the value of cos θ hel . This feature, previously observed through a model-independent analysis [12], is characteristic of CP violation originating from a sizeable interference between the spin-1 ρ(770) 0 resonance and the broad spin-0 contribution present in this region. This effect cancels when integrating over the helicity angle, since the interference term is proportional to cos θ hel . In addition, the change in the asymmetry below and above the ρ(770) 0 peak indicates that the effect is mediated by a strong phase difference dominated by the evolution of the ρ(770) 0 Breit-Wigner amplitude phase. All three approaches to the modelling of the π + π − S-wave describe this effect well.

8.2
The π + π − S-wave A notable feature in Fig. 7(c) is the small but approximately constant asymmetry at m(π + π − ) values below the ρ(770) 0 mass. This region is dominated by the S-wave component with a small contribution from the ρ(770) 0 low-mass tail; the CP asymmetry in the S-wave in this region is also seen in all three approaches in Fig. 13(a), (c) and (e). A CP asymmetry in the S-wave below the inelastic (KK) threshold cannot be explained via a final-state interaction mechanism, and therefore has a different origin to effects seen elsewhere in the Dalitz plot.
The combined significance of CP violation in the S-wave and in the interference between the S-and P-waves is evaluated from the change in log-likelihood between the baseline fit and in fits where all relevant CP -violation parameters are fixed to be zero. Since the ρ-ω component serves as the reference amplitude in the baseline fit, this means that all δx and δy parameters, defined in Eq. (3), associated with the S-wave are fixed to zero. This is done in fits with each of the approaches to the S-wave model, with the resulting change in log-likelihood converted into a p-value, and subsequently into the number of Gaussian standard deviations (σ), accounting for the number of fixed parameters. The values obtained are around 30σ in all cases, despite the very different number of degrees of freedom associated with the S-wave in the different approaches. While this method can only be considered to give an approximation to the significance, it is sufficient to establish the presence of CP violation far beyond any reasonable doubt.
In order to separate the effects of CP violation in the S-wave and in the interference between the S-and P-waves, additional fits are performed in which the reference amplitude is changed to the S-wave. The δx and δy parameters associated with the S-wave are then fixed to zero, while those associated with the P-wave are allowed to vary in the fits. In this case, CP violation in the interference between the S-and P-waves is allowed, while none is possible in the S-wave itself, and hence the significance of each effect individually can be assessed. The values obtained are above 10σ in each of the S-wave modelling approaches, thus establishing that both CP -violation effects are present.
At low m(π + π − ) values, the S-wave magnitude and phase motion of the three approaches broadly agree, particularly for the CP -averaged |A| 2 , and all models capture similar behaviour around 1 GeV/c 2 . However, in the KK threshold region shown in Fig. 12, the change in sign of the difference between the number of B + and B − candidates between positive and negative cos θ hel is captured only by the K-matrix model. It is worth noting that this is the only model with an explicit f 0 (980) term: the isobar model includes only ππ ↔ KK rescattering above the KK threshold and the QMI binning is not sufficiently fine in this region to resolve a narrow structure. At 1.5 GeV/c 2 , the K-matrix has a clear phase motion, seen in Fig. 13, that is associated with the f 0 (1500) contribution. Consistent behaviour is seen in the QMI approach, although the uncertainties preclude a definite corroboration of the presence of the f 0 (1500) state. The isobar model does not include this component explicitly, and therefore it is expected that the phase is broadly constant here, continuing to the upper kinematic boundary. Above 3 GeV/c 2 , the magnitude of the S-wave component in all three approaches is consistent with zero and therefore the phase values are dominated by statistical and systematic uncertainties.

The f 2 (1270) region
Despite broad consistency between the three fit models, a clear discrepancy with the data is apparent for all of them in the f 2 (1270) region shown in Fig. 8(a). All fit model projections lie under (over) the data below (above) the f 2 (1270) peak, which is set to the known value of 1275.5 ± 0.8 MeV/c 2 [51]. Better agreement with the data is obtained when the f 2 (1270) mass and width are allowed to vary in the fits, as shown in Fig. 15(a). However, the obtained masses, equal to 1256 ± 4 MeV/c 2 , 1252 ± 4 MeV/c 2 and 1260 ± 4 MeV/c 2 for the isobar, K-matrix and QMI S-wave approaches, respectively, where the uncertainties are statistical only, are at least four standard deviations away from the world average. The  values obtained for the width are, however, consistent with the world average. Moreover, if the f 2 (1270) mass and width are allowed to vary independently in the B − and B + subsamples, inconsistent values are obtained. Alternatively, the discrepancy between the data and the models can be reduced by adding another spin-2 resonance in the f 2 (1270) region, as shown in Fig. 15(b). The established states f 2 (1525), with m f 2 (1525) = 1525 ± 5 MeV/c 2 and Γ f 2 (1525) = 73 +6 −5 MeV, and f 2 (1565), with m f 2 (1565) = 1562 ± 13 MeV/c 2 and Γ f 2 (1565) = 134 ± 8 MeV, are too high in mass and too narrow to be likely to induce a significant effect in the f 2 (1270) region. Therefore, an additional spin-2 resonance with mass and width parameters free to vary in the fit is introduced, with initial values corresponding to those of the not-well-established f 2 (1430) resonance. The fit results for the mass are consistent between each S-wave approach, with m f 2 (X) = 1600 ± 60, 1541 ± 24 and 1560 ± 14 MeV/c 2 for the isobar, K-matrix and QMI fits, respectively. However, the obtained values for the width are inconsistent, varying between Γ f 2 (X) = 367 ± 97, 204 ± 78 and 115 ± 37 MeV, where the uncertainties are statistical only. Therefore the addition of a second spin-2 state does not appear in the baseline model, but rather is considered as a source of systematic uncertainty on the model composition.
Consequently, the baseline model includes in the ππ D-wave only the f 2 (1270) resonance, with its mass and width fixed to the known values. Analysis of larger data samples will be required to obtain a more detailed understanding of the ππ D-wave in B + → π + π + π − decays.
The effect of additional decay channels on the energy-dependent width of the f 2 (1270) is also considered as another possibility for the baseline fit model discrepancy. A global analysis is performed to express the known branching fractions of f 2 (1270) decays to ππ, KK, ηη, π + π − π 0 π 0 , π + π − π + π − and π 0 π 0 π 0 π 0 in terms of their respective couplings in f 2 (1270) decays to these final states. Subsequent fits accounting for the energy-dependent width in a similar way as in the Flatté lineshape [96] are found to have minimal impact on the model and therefore do not contribute to the systematic uncertainties.
Despite the considerations outlined above, the CP violation associated with the f 2 (1270) resonance is robust with respect to the experimental and model systematic uncertainties documented in Section 6. This can be seen by comparing the coefficients that describe the f 2 (1270) resonance with those obtained during systematic variations as shown in Fig. 16. The fact that the contours for B + and B − coefficients do not overlap is a visual representation of the significantly non-zero values of the δx and δy parameters of Eq. (3). The quasi-two-body CP asymmetry, defined in Eq. (33), is related to the difference in the magnitudes of the B + and B − complex coefficients, i.e. the difference in distance from (0, 0) to the centres of the corresponding ellipses in Fig. 16. In addition to this, CP violation can also be observed in the difference in the phases relative to the ρ(770) 0 -ω(782) reference component in the B + and B − amplitudes, which would manifest in the Dalitz plot as a difference between B + and B − decays in the interference pattern between ρ(770) 0 and f 2 (1270) resonances. This is indicated by the difference in the polar angle in the Argand plane of Fig. 16. It is evident that systematic variations do not significantly modify the distance between the solid and dashed ellipses in the (x, y) plane, meaning that the significant overall CP violation associated with the f 2 (1270) resonance is robust.
When interpreting Fig. 16, it should be noted that in the QMI approach the amplitude components are not individually normalised, unlike the case for the isobar and K-matrix approaches. Since some systematic variations can change the overall scale of various lineshapes, their respective deviations as depicted in this plot cannot be directly interpreted as entirely systematic in origin, and as such naturally manifest as the largest deviations from the nominal model.
The statistical significance of CP violation in B + → f 2 (1270)π + is estimated by finding, for each S-wave method, the difference of twice the log-likelihood between two fits: the nominal one and another where the CP -violating parameters of the f 2 (1270) are fixed to zero. Since this quantity behaves like a χ 2 distribution with two degrees of freedom, it is converted into a p-value from which confidence intervals are derived. The significance of CP violation is found to be 20σ, 15σ and 14σ for the isobar, K-matrix and QMI approaches, respectively. This corresponds to the first observation of CP violation in B + → f 2 (1270)π + decay, which is the first observation of CP violation in any process with a final state containing a tensor resonance. The measured central value of A CP (B + → f 2 (1270)π + ) is consistent with some theoretical predictions [24,97,98] that, however, have large uncertainties.

The ρ 3 (1690) 0 region
The contribution from a spin-3 ρ 3 (1690) 0 component is evident in Fig. 11, with a dip in intensity at cos θ hel = 0, characteristic of an odd-spin resonance, as well as multiple "lobes" associated with a spin greater than 1. The central value of the CP asymmetry of this component is large and positive, however its systematic uncertainty is also large, mostly driven by ambiguities in the amplitude model. These uncertainties preclude any conclusive statement about CP violation in B + → ρ 3 (1690) 0 π + decays; an analysis with a larger data sample will be necessary to clarify this point.

Conclusions
An amplitude analysis of the B + → π + π + π − decay is performed, using a data sample corresponding to 3 fb −1 collected by the LHCb experiment during Run 1. Three complementary approaches are used to describe the large S-wave contribution to this decay. Overall good agreement is found between all three models and the data, although some discrepancies in the region around the f 2 (1270) region persist in the baseline models.
Significant CP violation associated with the f 2 (1270) resonance is observed, the first observation of CP violation in any process containing a tensor resonance, which is robust with respect to systematic uncertainties related to both experimental effects and to the composition of the amplitude model. The quasi-two-body CP asymmetry in the B + → ρ(770) 0 π + decay is found to be compatible with zero. However, CP -violation effects that are characteristic of interference between the spin-1 ρ(770) 0 resonance and the spin-0 S-wave contribution are observed, and are well described in all three approaches to the S-wave. This is the first observation of CP violation mediated entirely by interference between hadronic resonances.
All three approaches to the description of the π + π − S-wave broadly agree on the variation of its magnitude and phase. One striking feature is the presence of a significant CP asymmetry in the S-wave that is not associated with the aforementioned interference effect, of which a substantial component is at low m(π + π − ). Further phenomenological and experimental investigations will be required to better understand the underlying dynamics of these and other effects in B + → π + π + π − decays, and to elucidate connections with CP -violation effects observed in B + → K + K + π − decays [13].
The results of this analysis provide valuable input to phenomenological work on the underlying nature of the remarkably large CP violation observed in the charmless threebody decays of the charged B meson, and in B-meson decays in general. Furthermore, the robust description of low mass S-wave achieved with the approaches documented here gives insight into the effects of low-energy QCD in B-meson decays.

Appendices A Isobar model tables
The results for the Cartesian coefficients obtained from the fit with the isobar description of the S-wave are reported in Table 19. The correlation matrices for the CP -averaged fit fractions and quasi-two-body decay CP asymmetries, corresponding to those presented in Tables 10 and 17, can be found in Tables from 19 to 23. As an indication of fit quality, signed χ 2 distributions in the square Dalitz plot, separated by charge, are produced with an adaptive binning procedure requiring at least 15 events per bin. These are shown in Fig. 17.

B K-matrix model tables
The results for the Cartesian coefficients, and other fitted parameters, obtained with the K-matrix S-wave approach can be found in Table 24. Correlation matrices for these parameters are reported in Ref. [99]. Here, as the reference amplitude is the ρ-ω mixing component, rather than a component only representing the ρ(770) 0 resonance, the magnitude of the positive isobar coefficient describing the ρ(770) 0 resonance is not unity, but is calculated accounting for the small ω(782) contribution. Furthermore, the statistical and systematic uncertainty correlation matrices for the K-matrix fit CP -averaged fit fractions and quasi-two-body decay CP asymmetries, corresponding to those presented in Tables 10 and 17, can be found in Tables 25 and 26, and 27 and 28, respectively. As an indication of fit quality, signed χ 2 distributions in the square Dalitz plot, separated by charge, are produced with an adaptive binning procedure requiring at least 15 events per bin. These are shown in Fig. 18.

B.1 Secondary minimum
A secondary minimum is observed in the maximum likelihood fit of the model with the K-matrix S-wave approach, located approximately 0.8 units of negative-log-likelihood away from the primary minimum. This minimum results in fit results that are statistically consistent with the best minimum, except for in the parameters of the individual K-matrix components (fit fractions and overall CP -violation parameters are otherwise consistent).
The parameters obtained from this secondary solution can be seen in Table 29 (to be compared to the nominal results in Table 24). Projections of the S-wave amplitude on m(π + π − ) can be seen in Fig. 19.

C QMI model tables
The results for the Cartesian coefficients from the fit with the QMI S-wave can be found in Table 30. Correlation matrices for these parameters are reported in Ref. [99]. The statisical and systematic correlation matrices for the CP -averaged fit fractions are given in Tables 31 and 32, respectively, while those for the quasi-two-body decay CP asymmetries can be found in Tables 33 and 34. As an indication of fit quality, signed χ 2 distributions in the square Dalitz plot, separated by charge, are produced with an adaptive binning procedure requiring at least 15 events per bin. These are shown in Fig. 20.   S-wave +1.00

D Results with S-wave model variation included as systematic uncertainty
Results are presented throughout this paper for each of three different approaches to the modelling of the ππ S-wave: the isobar, K-matrix and QMI models. As discussed in Sections 7 and 8, all three give good descriptions of the data, with each describing some regions of the Dalitz plot better than the others. Therefore, it is not possible to conclude that one is preferred to the others. Nonetheless, it is anticipated that for some purposes it will be more useful to have a single set of results rather than three sets. Therefore, Table 35 provides such a presentation. The central values, statistical and experimental systematic uncertainties are taken from the results with the isobar model, while the largest deviation in the central value between the isobar model and the other two S-wave approaches is combined in quadrature with the other sources of model uncertainty.

E Phase comparison
The presentation of the complex coefficients c j in the Cartesian convention makes it difficult to compare the relative phases of the components in the different models. To facilitate this, the relevant information is presented in Table 36.

Supplemental material
In addition to the results presented in the main body of the paper, we provide a supplementary collection of files containing correlation matrices that are impractical to publish either in the main text or appendices. Files are provided in the JSON (JavaScript Object Notation) format, which is both machine and human readable, and contain parameter names, two-dimensional arrays corresponding to the correlations between parameters in the order given in the parameter list, and maps from strings describing entries in the correlation matrix and the associated value. Three files obtained from the QMI approach, related to tables presented in Section 7, are given here. These are comprised of statistical and systematic correlation matrices of the parameters given in Tables 15 and 16 in file qmi_ff_corr.json, statistical uncertainties on the isobar coefficients in file qmi_params_corr.json, along with statistical and systematic correlation matrices for the parameters in Table 18 in file qmi_SWaveParams_corr.json.
For the K-matrix fit model results, statistical correlation matrices for the values given in Tables 10, 13, 14, and 17, along with for the isobar coefficients listed in Table 24, can be found in file kMatrix_stat_matrices.json.