Amplitude analysis of the 
B+→π+π+π−
 decay

The results of an amplitude analysis of the charmless three-body decay 
B 
+ 
→ 
π 
+ 
π 
+ 
π 
− 
, in which 
C 
P 
-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 
p 
p 
 collisions recorded with the LHCb detector. The most challenging aspect of the analysis is the description of the behavior 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 
f 
2 
( 
1270 
) 
 resonance in the 
π 
+ 
π 
− 
 
D 
-wave, and the 
ρ 
3 
( 
1690 
) 
0 
 resonance in the 
π 
+ 
π 
− 
 
F 
-wave. Significant 
C 
P 
-violation effects are observed in both 
S 
- and 
D 
-waves, as well as in the interference between the 
S 
- and 
P 
-waves. The results from all three approaches agree and provide new insight into the dynamics and the origin of 
C 
P 
-violation effects in 
B 
+ 
→ 
π 
+ 
π 
+ 
π 
− 
 decays.


I. 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 localized in regions of phase space of charmless three-body B-meson decays have been observed in model-independent analyses [10][11][12][13],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 [14] 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 two-dimensional phase space known as the Dalitz plot [15,16]. 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 [17][18][19][20]. The amplitude analysis performed by the BABAR Collaboration [20] additionally observed a large S-wave contribution; however, measurements of CPviolating 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 center-ofmass energy of 7 TeV and 2 fb −1 was collected in 2012 with a center-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 [36].T h r e e different state-of-the-art approaches to the modeling 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 organized as follows: Sec. II gives a brief description of the LHCb detector and the event reconstruction and simulation software; the signal candidate selection procedure is described in Sec. III; Sec. IV describes the procedure for estimating the signal and background yields that enter into the amplitude fit; Sec. V outlines the formalism used for the construction of the amplitude models, as well as a description of the mass line shapes used to parametrize the intermediate structures; Sec. VI describes the systematic uncertainties associated with the analysis procedure; Sec. VII 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 Sec. VIII; and a summary of the work as a whole can be found in Sec. IX.As h o r t e r description of the analysis, more focused on the first observations of different sources of CP-violation effects, can be found in a companion article [37].

II. DETECTOR AND SIMULATION
The LHCb detector [38,39] 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 silicon-strip vertex detector surrounding the pp interaction region [40], 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 [41] 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, 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 [42]. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillatingpad 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 [43]. 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 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 [44] 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 are reconstructed, for data collected in 2011 (2012), in events where at least one charged particle has transverse momentum p T > 1.7ð1.6Þ GeV=c and is inconsistent with originating from a PV. A multivariate algorithm [45] 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 [46] with a specific LHCb configuration [47]. Decays of unstable particles are described by EVTGEN [48], in which final-state radiation is generated using PHOTOS [49]. The interaction of the generated particles with the detector and its response are implemented using the GEANT4 toolkit [50] as described in Ref. [51].

III. SELECTION
The selection of signal candidates follows closely the procedure used in the model-independent analysis of the same data sample [12]. Signal B þ candidates are formed from three good-quality tracks that are consistent with originating from the same secondary vertex (SV), which 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 decision-tree classifier [52] 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 optimized to maximize the expected approximate signal significance, N s = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N s þ N b p , where N s is the expected signal yield within 40 MeV=c 2 of the known B þ mass [53], 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 particleidentification 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 intermediateD 0 mesons and decays ofD 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 , with a corresponding full width at half maximum of better than 7 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.

IV. B + CANDIDATE INVARIANT-MASS FIT
An extended, unbinned, maximum-likelihood fit is performed to the mðπ þ π þ π − Þ invariant-mass 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, to permit correction for possible differences in efficiency between subsamples (see Sec. VE).
The shape of the B þ → π þ π þ π − signal decay is parametrized by the sum of a core Gaussian with two Crystal Ball functions [54], with tails on opposite sides of the peak 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 modeled with an ARGUS function [55] convolved with a Gaussian resolution function. The smooth combinatorial background is modeled 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 [56], 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 [57]; the details of the model used do not impact this analysis. This shape is modeled 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 Sec. VE. The mass fit results are shown in Fig. 1, while Table I 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 Sec. VI.

V. DALITZ-PLOT MODEL
The B þ → π þ 1 π þ 2 π − 3 decay amplitude can be expressed fully in terms of the invariant mass squared 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 visualization 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 line shape 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 Figs. 2(a) and 2(b).

A. 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 parametrizes the intermediate resonant or nonresonant processes, Combinatorial background asymmetry þ0.005 AE 0.010 B þ → K þ π þ π − background asymmetry 0.000 AE 0.008 FIG. 2. Conventional Dalitz-plot distributions for (a) B þ and (b) B − , and square Dalitz-plot (defined in Sec. VA 1) distributions for (c) B þ and (d) B − candidate decays to π AE π þ π − . Depleted regions are due to theD 0 veto. 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 parametrized as where R describes the mass line shape, T describes the angular dependence, and X are Blatt-Weisskopf barrier factors [58] depending on a radius parameter r BW .H e r e 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 )a s ⃗p, where both momenta are evaluated in the rest frame of the dipion system.
Using the Zemach tensor formalism [59,60], 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 provides 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 ¼j⃗ qjr P BW for the B þ decay and z ¼j⃗pjr 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 of the B þ meson and the intermediate resonance is taken to be r P

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 [19], which provides improved resolution in these critical regions when using uniform binning, for example when modeling efficiency effects. Furthermore, the mapping to a square space aligns the bin boundaries to the kinematic boundaries of the phase space. As such, all efficiencies and backgrounds described in Secs. VEand VFare determined as a function of the square Dalitz-plot variables; however the data and fit projections in Sec. VII are displayed as a function of the squares of the invariant-mass pairs.
The square Dalitz plot is defined in terms of m 0 and θ 0 where mðπ þ π þ Þ max ¼ m B þ − m π − and mðπ þ π þ Þ min ¼ 2m π þ represent the kinematic limits permitted in the B þ → π þ π þ π − decay and θðπ þ π þ Þ is the angle between π þ and π − in the π þ π þ rest frame. The Bose symmetry of the final state requires that distributions are symmetric with respect to θ 0 ¼ 0.5. The square Dalitz-plot distributions of the selected candidates can be seen in Figs. 2(c) and 2(d).

B. Mass line shapes
Resonant contributions are mostly described by the relativistic Breit-Wigner line shape with a mass-dependent decay width where q 0 is the value of q ¼j⃗ qj 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 [61]. It takes the form with an additional mass dependence 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 parametrizes the interference between these two contributions following Ref. [62], where R ρ ðmÞ is the Gounaris-Sakurai ρð770Þ 0 line shape, R ω ðmÞ is the relativistic Breit-Wigner ωð782Þ line shape, jζj 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 [53]. This is equivalent to the parametrization described in Ref. [63] if the small Δ 2 term in the denominator is ignored. The value for δ is fixed in the fit to δ ¼ 0.00215 AE 0.00035 GeV [62].

C. 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 behavior. 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, unitaritypreserving model informed by historical scattering data (K-matrix), and (iii) a quasi-model-independent (QMI) binned approach. 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 σ,o rf 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 [36,64], parametrized as where s σ is the square of the pole position ffiffiffiffi ffi extracted from the fit. The concept of ππ ↔ KK rescattering was originally developed in 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,35], The amplitude f rescatt ðmÞ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 − ηðmÞ 2 p 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 analyticity with dispersion relation techniques [65]. The inelasticity is described by with where m K ¼ 0.494 GeV is the charged kaon mass, The phase shift is given by where Except for m K , values of all parameters are taken from Ref. [65].

K-matrix model
The coherent sum of resonant contributions modeled with Breit-Wigner line shapes 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), two-body unitarity is naturally conserved within the K-matrix approach [66]. This approach was originally developed for two-body scattering [67] and the study of resonances in nuclear reactions [68,69], but was extended to describe resonance production and n-body decays in a more general way [70]. This section provides a brief introduction to the K-matrix approach as applied in this analysis; for more details see Ref. [71].
From unitarity conservation, the form factor for twobody 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 [53] 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 ijρ uu j when the channel is below its mass threshold (provided it does not cross into another channel). For the coupled multimeson final states, the corresponding expression can be found in Ref. [71].
The scattering matrix,K, can be parametrized 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 [72]. 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 [73]. 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 resonancesmixtures of bare states.
An 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 the production vector,P, 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. [74].
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. [75] and given in Table II (for five channels n ¼ 5 and five poles N p ¼ 5), which are obtained from a global analysis of ππ scattering data [72].

Quasi-model-independent analysis
In the 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 modeled 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 optimization 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π [76][77][78][79], ππ [80] and Dπ [81] 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 the premise of smoothness, which is appropriate for goals such as the 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 Swave distributions obtained by rescattering experiments.

D. Measurement quantities
The primary outputs of the Dalitz-plot fit are the complex isobar coefficients c AE j defined in Eq. (1). However, since these depend on the choice of phase convention, amplitude formalism and normalization, 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  [75], which are obtained from a global analysis of ππ scattering data [72]. 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. [74]. 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 line shapes 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-twobody 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 AE S is the coherent sum of all contributions to the S-wave.

E. Efficiency model
The efficiency of selecting a signal decay is parametrized in the two-dimensional square Dalitz plot and determined separately for B þ and B − decays. Nonuniformities 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 dipolemagnet 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 [82]. The particle identification efficiencies for the background-subtracted pions and kaons are parametrized in terms of their total and transverse momentum, and the number of tracks in the event. This efficiency is assumed to factorize 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 [83], 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 twodimensional cubic spline to mitigate effects of discontinuity at the bin edges, with bins abutting kinematic boundaries repeated to ensure good behavior at the edge of the phase space. The signal PDF for B þ or B − decays is then given by where ϵ AE represents the Dalitz-plot-dependent efficiency for the B AE decay.

F. Background model
The dominant source of background in the signal region is combinatorial in nature. In the Dalitz-plot fit, the AMPLITUDE ANALYSIS OF THE B þ → π þ π þ π − … PHYS. REV. D 101, 012006 (2020) distribution of this background is modeled 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 0 and θ 0 . 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 modeled 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 Sec. IV.
A secondary source of background arises from misidentified B þ → K þ π þ π − decays. This background is modeled 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 [56]. 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. 5.

G. 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 [71], which interfaces to the MINUIT function minimization algorithm [84,85]. In contrast, the QMI approach relies on the MINT2 [86] amplitude-analysis interface to MINUIT2 [85]. The fundamental difference between these amplitude-analysis software packages is in the handling of the normalization. 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 optimized is of the form where N k is the yield for the candidate category k (given in Table I), N is equal to P 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 minimizing twice the negative loglikelihood, −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 randomized. 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, and solutions near to this are also inspected (see Appendix B1 for the K-matrix model fit).

H. 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 [17][18][19][20], additional resonances are examined iteratively, in the order that maximizes 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 [87], 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 line shapes for resonances already included in the model, and the investigation of additional, more speculative, states. These form the basis of several important systematic uncertainties listed in Sec. VI and are further discussed in Sec. VIII. Last, 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, 0.050, 0.100, 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 (see Appendix A for the distribution for the isobar model fit, Appendix B for the distribution for the K-matrix model fit, and Appendix C for the distribution for the QMI model fit). 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 Sec. VII, 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 mismodeling. Nevertheless, a discrepancy between all of the models and the data is observed in the region around the f 2 ð1270Þ resonance, which is investigated further in Sec. VIII C.
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 III. Furthermore, the mass and width of the dominant ρð770Þ 0 contribution are left free to vary, which results in a significantly better fit quality.

VI. 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 I, comprise a statistical component as well as systematic effects due to the invariant-mass fit procedure. The uncertainty arising from assumptions regarding the signal parametrization 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 [88,89]. 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 [56], 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,  Group (PDG) tables [53]. 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 mismodeling in the f 2 ð1270Þ region, discussed in Sec. VIII, 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 [90], 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 parametrization 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 parametrization 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 is assigned as a systematic uncertainty. The parameters of the σ contribution to the S-wave are also varied within the uncertainties on the worldaverage mass and width, and the effect on the results is taken as a systematic uncertainty.
For fits using the K-matrix approach, both the fourtĥ P-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 TABLE IV. Systematic uncertainties on the CP-averaged fit fractions, given in units of 10 −2 , for the isobar method. Uncertainties are given both for the total S-wave, and for the individual components due to the σ pole and the rescattering amplitude. For comparison, the statistical uncertainties are also listed at the bottom.

Category
ρð770Þ 0 ωð782Þ   extent to which the S-wave interferes with other partial waves in a particular bin, and the approximation of an analytic line shape 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 recognized 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 summarized in Tables IV and V for the isobar approach,  Tables VI and VII for the K-matrix approach, and  Tables VIII and IX 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.

VII. RESULTS
Numerical results for the fit fractions and quasi-twobody CP asymmetries are given in this section; Appendixes A-C also provide the correlation matrices for the CP-averaged fit fractions and CP asymmetries, while Appendix D presents the combined values for the three S-wave approaches. 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 Appendixes A-C for the isobar, K-matrix and QMI approaches, respectively. Throughout the tables in this section, separate parameters for the ρð770Þ 0 and ωð782Þ resonances are presented, which are extracted from the combined component described by Eq. (15). 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 focusing 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 Appendixes F-H.

A. Fit fractions
The CP-averaged fit fractions are given in Table XI for all three S-wave approaches. The fit fractions and interference fit fractions, separated by B AE charge for each S-wave approach, are given in Tables XII-XVII. 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 model sources have been combined in quadrature for brevity. The total sums of fit fractions for the B þ and B − amplitudes can be found in Table X.

B. CP asymmetries
Quasi-two-body CP asymmetries associated to each component are shown in Table XVIII. A detailed discussion of these results is given in Sec. VIII; however it should be stressed that CP-violation effects can manifest in Dalitzplot 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 phase-space-integrated CP asymmetry is consistent, in all three models, with the value previously determined through model-independent analysis [12].

C. S-wave projections
The squared amplitude and phase motion of the S-wave models as a function of mðπ þ π − Þ can be seen in Figs. 13(a) and 13(b) for the isobar approach, Figs. 13(c) and 13(d) for the K-matrix approach and Figs. 13(e) and 13(f) for the QMI approach. A comparison of all three models, for the CP-averaged S-wave projections, can be seen in Fig. 14.TheQMIS-wave is recorded in Table XIX, while the statistical and systematic correlation matrices obtained with this approach are given in the Supplemental Material [91].

D. ρð770Þ 0 mass and width
The ρð770Þ 0 mass and width are allowed to vary freely in each fit as mentioned in Sec. VH. The fitted results are consistent with the world-average values, and can be seen in Table XX.

E. Multiple fit solutions
A search for secondary solutions with negative loglikelihood 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 secondary minimum closest to the most-negative minimum, with a change in log-likelihood of 0.8, are given in Appendix B.

VIII. DISCUSSION
The results and figures presented in Sec. VII 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 sizable asymmetries in localized regions of the Dalitz plot.
A. 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 XVIII. A number of theoretical calculations of these quantities are available in the literature, with some authors [92][93][94][95] predicting values for A CP ðB þ → ρð770Þ 0 π þ Þ that are consistent with the measured result, albeit sometimes with large uncertainties. Other approaches [96][97][98][99] 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 [27][28][29].
A significant CP asymmetry in the ρð770Þ 0 region can, however, be seen in the cos θ hel projections shown in Figs. 10(c) and 10(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 sizable 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 modeling of the π þ π − S-wave describe this effect well.
B. 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 lowmass tail; the CP asymmetry in the S-wave in this region is also seen in all three approaches in Figs. 13(a), 13(c) and 13(e).ACP asymmetry in the S-wave below the inelastic (KK) threshold cannot be explained via ππ ↔ KK rescattering, and therefore has a different origin to effects seen elsewhere in the Dalitz plot.
The combined significance of CP violation in the Swave 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 modeling 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 jAj 2 , and all models capture similar behavior around 1 GeV=c 2 .H o w e v e r ,i n 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 behavior 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. C. 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 AE 0.8 MeV=c 2 [53]. 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 AE 4, 1252 AE 4 and 1260 AE 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 0 2 ð1525Þ, with m f 0 2 ð1525Þ ¼ 1525 AE 5 MeV=c 2 and Γ f 0 2 ð1525Þ ¼ 73 þ6 −5 MeV, and f 2 ð1565Þ, with m f 2 ð1565Þ ¼ 1562 AE 13 MeV=c 2 and Γ f 2 ð1565Þ ¼ 134 AE 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-wellestablished f 2 ð1430Þ resonance. The fit results for the mass are consistent between each S-wave approach, with m f 2 ðXÞ ¼ 1600 AE 60, 1541 AE 24 and 1560 AE 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 AE 97, 204 AE 78 and  115 AE 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 ππDwave 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. TABLE XI. The CP-averaged fit fractions in units of 10 −2 , for each approach, where the first uncertainty is statistical, the second is the experimental systematic and the third is the model systematic.

Component
Isobar 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 is the quadratic sum of systematic and model sources.   TABLE XIV. 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 is the quadratic sum of systematic and model sources. The effect of additional decay channels on the energydependent 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é line shape [100] 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 Sec. VI. 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 nonzero 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 XV. 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 is the quadratic sum of systematic and model sources. TABLE XVI. 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 is the quadratic sum of systematic and model sources.
ρð770Þ 0 -ωð782Þ f 2 ð1270Þ ρð1450Þ 0 ρ 3 ð1690Þ 0 S-wave TABLE XVII. 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 is the quadratic sum of systematic and model sources.  B þ and B − complex coefficients, i.e., the difference in distance from (0,0) to the centers 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 normalized, unlike the case for the isobar and K-matrix approaches. Since some systematic variations can change the overall scale of various line shapes, 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 [26,101,102] that, however, have large uncertainties. shows the phase motion. Discontinuities in the phase motion are due to presentation in the range ½−180°; 180°. 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 jAj 2 plot in each approach is proportional to its respective S-wave fit fraction, the overall scale of the Kmatrix and QMI plots are set relative to the isobar S-wave fit fraction in order to facilitate comparison between the three approaches. TABLE XIX. QMI S-wave fit results where the first uncertainty is statistical and the second is the quadratic sum of systematic and model sources.    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.

IX. CONCLUSIONS
An amplitude analysis of the B þ → π þ π þ π − decay was performed, using a data sample corresponding to 3 fb −1 collected by the LHCb experiment during Run 1. Three complementary approaches were used to describe the large S-wave contribution to this decay. Overall good agreement was 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 was 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 was 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 were 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 π þ π − Swave 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 [14].
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 three-body 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.

APPENDIX 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 XXI, with the fitted pole parameters given in Table XXII. The correlation matrices for the CP-averaged fit fractions and quasi-two-body decay CP asymmetries, corresponding to those presented in Tables XI and XVIII,  can be found in Tables XXI-XXVI. 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.

APPENDIX 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 XXVII. Correlation matrices for these parameters are reported in Ref. [91]. 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 XI and XVIII, can be found in Tables XXVIII and XXIX,a n dXXX and XXXI, 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.

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 negativelog-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 XXXII (to be compared to the nominal   results in Table XXVII). Projections of the S-wave amplitude on mðπ þ π − Þ can be seen in Fig. 19.

APPENDIX C: QMI MODEL TABLES
The results for the Cartesian coefficients from the fit with the QMI S-wave can be found in Table XXXIII. Correlation matrices for these parameters are reported in Ref. [91]. The statistical and systematic correlation matrices for the CPaveraged fit fractions are given in Tables XXXIV and  XXXV, respectively, while those for the quasi-two-body decay CP asymmetries can be found in Tables XXXVI and XXXVII. 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.     Results are presented throughout this paper for each of three different approaches to the modeling of the ππS-wave: the isobar, K-matrix and QMI models. As discussed in Secs. VII and VIII, 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 XXXVIII 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.

APPENDIX 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 XXXIX.   TABLE XXXVIII. Results with S-wave model variation included as a source of systematic uncertainty. The first uncertainty is statistical, the second is experimental systematic and the third is the adjusted model systematic uncertainty.