Measurement of resonant and $CP$ components in $\overline{B}_s^0\rightarrow J/\psi \pi^+\pi^-$ decays

The resonant structure of the decay $\overline{B}_s^0\to J/\psi\pi^+\pi^-$ is studied using data corresponding to 3 fb$^{-1}$ of integrated luminosity from $pp$ collisions by the LHC and collected by the LHCb detector. Five interfering $\pi^+\pi^-$ states are required to describe the decay: $f_0(980),f_0(1500),f_0(1790),f_2(1270)$, and $f_2^{\prime}(1525)$. An alternative model including these states and a non-resonant $J/\psi \pi^+\pi^-$ component also provides a good description of the data. Based on the different transversity components measured for the spin-2 intermediate states, the final state is found to be compatible with being entirely $CP$-odd. The $CP$-even part is found to be $<2.3$% at 95% confidence level. The $f_0(500)$ state is not observed, allowing a limit to be set on the absolute value of the mixing angle with the $f_0(980)$ of $<7.7^{\circ}$ at 90% confidence level, consistent with a tetraquark interpretation of the $f_0(980)$ substructure.


Introduction
CP violation studies in the B 0 s → J/ψ π + π − decay mode complement studies using B 0 s → J/ψ φ and improve the final accuracy in the CP -violating phase, φ s , measurement [1]. While the CP content was previously shown to be more than 97.7% CP -odd at 95% confidence level (CL), it is important to determine the size of any CP -even components as these could ultimately affect the uncertainty on the final result for φ s . Since the π + π − system can form light scalar mesons, such as the f 0 (500) and f 0 (980), we can investigate if these states have a quark-antiquark or tetraquark structure, and determine the mixing angle between these states [2]. The tree-level Feynman diagram for the process is shown in Fig. 1. We have previously studied the resonance structure in B 0 s → J/ψ π + π − decays using data corresponding to an integrated luminosity of 1 fb −1 [3]. 1 In this paper we use 3 fb −1 of luminosity, and also change the analysis technique substantially. Here the π + π − mass, and all three decay angular distributions are used to determine the resonant and nonresonant components. Previously the angle between the decay planes of J/ψ → µ + µ − and π + π − in the B 0 s rest frame, χ, was integrated over. This simplified the analysis, but sacrificed some precision and also prohibited us from measuring separately the helicity +1 and −1 components of any π + π − resonance, knowledge of which would permit us to evaluate the CP composition of resonances with spin greater than or equal to 1. Since one of the particles in the final state, the J/ψ, has spin-1 its three decay amplitudes must be considered, while the π + π − system is described as the coherent sum of resonant and possibly non-resonant amplitudes.
2 Amplitude formula for B 0 s → J/ψ h + h − The decay of B 0 s → J/ψ h + h − , where h denotes a pseudoscalar meson, followed by J/ψ → µ + µ − can be described by four variables. We take the invariant mass of h + h − (m hh ) and three helicity angles defined as (i) θ J/ψ , the angle between the µ + direction in the J/ψ rest frame with respect to the J/ψ direction in the B 0 s rest frame; (ii) θ hh , the angle between the h + direction in the h + h − rest frame with respect to the h + h − direction in the B 0 s 1 Charged conjugated modes are also used when appropriate. rest frame, and (iii) χ, the angle between the J/ψ and h + h − decay planes in the B 0 s rest frame. Figure 2 shows these angles pictorially 2 . In this paper hh is equivalent to π + π − .
From the time-dependent decay rate of ( -) B 0 s → J/ψ h + h − derived in Ref. [4], the timeintegrated and flavor-averaged decay rate is proportional to the function S(m hh , θ hh , θ J/ψ , χ) =|A(m hh , θ hh , θ J/ψ , χ)| 2 + |A(m hh , θ hh , θ J/ψ , χ)| 2 − 2D Re q p A * (m hh , θ hh , θ J/ψ , χ)A(m hh , θ hh , θ J/ψ , χ) , (1) where ( -) A, the amplitude of ( -) B 0 s → J/ψ h + h − at proper time t = 0, is a function of m hh , θ J/ψ , θ hh , χ, and is summed over all resonant (and possibly non-resonant) components; q and p are complex parameters that describe the relation between mass and flavor eigenstates [5]. The interference term arises because we must sum the B 0 s and B 0 s amplitudes before squaring. Even when integrating over proper time, the terms proportional to sinh (∆Γ s t/2) do not vanish because of the finite ∆Γ s in the B 0 s system, where ∆Γ s is the width difference between the light and the heavy mass eigenstates. The factor D is where Γ s is the average B 0 s decay width, and ε(t) is the detection efficiency as a function of t. For a uniform efficiency, D = ∆Γ s /(2Γ s ) and is (6.2 ± 0.9)% [6].
The amplitude, A R (m hh ), is used to describe the mass line-shape of the resonance R, that in most cases is a Breit-Wigner function. It is combined with the B resonance decay properties to form the expression Here P B is the J/ψ momentum in the B 0 s rest frame, P R is the momentum of either of the two hadrons in the dihadron rest frame, m B is the B 0 s mass, J R is the spin of R, 2 These definitions are the same for B 0 s and B 0 s , namely, µ + and h + are used to define the angles in both cases. L B is the orbital angular momentum between the J/ψ and h + h − system, and L R the orbital angular momentum in the h + h − decay, and thus is the same as the spin of the are the Blatt-Weisskopf barrier factors for the B 0 s and R resonance, respectively [3]. The factor √ P R P B results from converting the phase space of the natural Dalitz-plot variables m 2 hh and m 2 J/ψ h + to that of m hh and cos θ hh [7]. We must sum over all final states, R, so for each J/ψ helicity, denoted by λ, equal to 0, +1, and −1 we have where ( -) h R λ are the complex coefficients for each helicity amplitude and the Wigner dfunctions are listed in Ref. [6].
The decay rates, | ( -) A(m hh , θ hh , θ J/ψ , χ)| 2 , and the interference term, A * (m hh , θ hh , θ J/ψ , χ)A(m hh , θ hh , θ J/ψ , χ), can be written as functions of ( -) H λ (m hh , θ hh ), θ J/ψ and χ. These relationships are given in Ref. [4]. In order to use the CP relations, it is convenient to replace the helicity complex coefficients ( -) h R λ by the complex transversity coefficients ( -) a R τ using the relations Here ( -) a R 0 corresponds to longitudinal polarization of the J/ψ meson, and the other two coefficients correspond to polarizations of the J/ψ meson and h + h − system transverse to the decay axis: ( -) a R for parallel polarization of the J/ψ and h + h − , and ( -) a R ⊥ for perpendicular polarization.
Assuming no direct CP violation, as this has not been observed in B 0 s → J/ψ φ decays [1], the relation between the B 0 s and B 0 s variables isā R τ = η R τ a R τ , where η R τ is CP eigenvalue of the τ transversity component for the intermediate state R, where τ denotes 0, , or ⊥ component. The final state CP parities for S, P, and D-waves are given in Table 1. Spin In this analysis a fit determines the amplitude strength a R τ and the phase φ R τ of the amplitude for each resonance R and each transversity τ . For the τ =⊥ amplitude, the L B value of a spin-1 (or -2) resonance is 1 (or 2); the other transversity components have two possible L B values of 0 and 2 (or 1 and 3) for spin-1 (or -2) resonances. In this analysis the lower one is used. It is verified that our results are insensitive to the L B choices.

Data sample and detector
The data sample corresponds to an integrated luminosity of 3 fb −1 collected with the LHCb detector [8] using pp collisions. One-third of the data was acquired at a center-ofmass energy of 7 TeV, and the remainder at 8 TeV. The detector 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, 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 [9] placed downstream. The combined tracking system provides a momentum 3 measurement with relative uncertainty that varies from 0.4% at 5 GeV to 0.6% at 100 GeV, and impact parameter (IP) resolution of 20 µm for tracks with large transverse momentum (p T ). Different types of charged hadrons are distinguished by information from two ring-imaging Cherenkov detectors (RICH) [10]. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [11]. The trigger consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage that applies a full event reconstruction [12]. Events selected for this analysis are triggered by a J/ψ → µ + µ − decay, where the J/ψ is required at the software level to be consistent with coming from the decay of a B 0 s meson by use either of IP requirements or detachment of the J/ψ from the primary vertex (PV). In the simulation, pp collisions are generated using Pythia [13] with a specific LHCb configuration [14]. Decays of hadronic particles are described by EvtGen [15], in which final state radiation is generated using Photos [16]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [17] as described in Ref. [18].
Preselection criteria are implemented to preserve a large fraction of the signal events, and are identical to those used in Ref. [19]. A B 0 s → J/ψ π + π − candidate is reconstructed by combining a J/ψ → µ + µ − candidate with two pions of opposite charge. To ensure good track reconstruction, each of the four particles in the B 0 s candidate is required to have the track fit χ 2 /ndf to be less than 4, where ndf is the number of degrees of freedom of the fit. The J/ψ → µ + µ − candidate is formed by two identified muons of opposite charge, having p T greater than 500 MeV, and with a geometrical fit vertex χ 2 less than 16. Only candidates with dimuon invariant mass between −48 MeV and +43 MeV from the observed J/ψ mass peak are selected, and are then constrained to the J/ψ mass [6] for subsequent use.
Pion candidates are required to each have p T greater than 250 MeV, and the sum, p T (π + ) + p T (π − ) larger than 900 MeV. Both pions must have χ 2 IP greater than 9 to reject particles produced from the PV. The χ 2 IP is computed as the difference between the χ 2 of the PV reconstructed with and without the considered track. Both pions must also come from a common vertex with χ 2 /ndf < 16, and form a vertex with the J/ψ with a χ 2 /ndf less than 10 (here ndf equals five). Pion candidates are identified using the RICH and muon systems. The particle identification makes use of the logarithm of the likelihood ratio comparing two particle hypotheses (DLL). For pion selection we require DLL(π − K) > −10 and DLL(π − µ) > −10.
The B 0 s candidate must have a flight distance of more than 1.5 mm. The angle between the combined momentum vector of the decay products and the vector formed from the positions of the PV and the decay vertex (pointing angle) is required to be less than 2.5 • .
Events satisfying this preselection are then further filtered using a multivariate analyzer based on a boosted decision tree (BDT) technique [20]. The BDT uses eight variables that are chosen to provide separation between signal and background. These are the minimum of DLL(µ − π) of the µ + and µ − , p T (π + ) + p T (π − ), the minimum of χ 2 IP of the π + and π − , and the B 0 s properties of vertex χ 2 , pointing angle, flight distance, p T and χ 2 IP . The BDT is trained on a simulated sample of B 0 s → J/ψ π + π − signal events and a background data sample from the sideband 5566 < m(J/ψπ + π − ) < 5616 MeV. Then the BDT is tested on independent samples. The distributions of BDT classifier for signal and background samples are shown in Fig. 3. By maximizing the signal significance we set the requirement that the classifier is greater than zero, which has a signal efficiency of 95% and rejects 90% of the background.
The invariant mass of the selected J/ψπ + π − combinations is shown in Fig. 4. There is a large peak at the B 0 s mass and a smaller one at the B 0 mass on top of a background. A double Crystal Ball function with common means models the radiative tails and is used to fit each of the signals. The known B 0 s − B 0 mass difference [6] is used to constrain the difference in mean values. Other components in the fit model take into account where the K − in the former, and both K − and p in the latter, are misidentified as pions. The shape of the B 0 → J/ψπ + π − signal is taken to be the same as that of the B 0 s . The combinatorial background shape is taken from like-sign combinations that are the sum of π + π + and π − π − candidates, and was found to be well described by an exponential function in previous studies [3,21]. The shapes of the other components are taken from simulation with their yields allowed to vary. The Λ 0 b → J/ψ K − p reflection yield in the fit region is constrained to the expected number 2145 ± 201, which is obtained from study of the events in the control region of 5066 < m(J/ψπ + π − ) < 5141 MeV. The mass fit gives 27396 ± 207 signal and 7075 ± 101 background candidates, leading to the signal fraction f sig = (79.5 ± 0.2)%, within ±20 MeV of the B 0 s mass peak. The effective r.m.s. mass resolution is 9.9 MeV.

Probability density function construction
The correlated distributions of four variables m hh , cos θ hh , cos θ J/ψ , and χ are fitted using the candidates within ±20 MeV of the B 0 s mass peak. To improve the resolution of these variables we perform a kinematic fit constraining the B 0 s and J/ψ masses to their world average mass values [6], and recompute the final state momenta.
The overall PDF given by the sum of signal, S, and background functions is where ε is the detection efficiency, and B is the background PDF discussed later in Sec. 5.3.
dot-long dashed curve is the B 0 → J/ψK − π + reflection and the (blue) solid curve is the total.
The normalization factor for signal is given by The signal function S is defined in Eq. (1), where D = (8.7 ± 1.5)%, taking into account the acceptance [22], and choosing a phase convention q/p = e −iφs . The phase φ s is fixed to the standard model value of −0.04 radians [23]. Our results are found to be insensitive to the value of φ s used within the 95% CL limits set by the LHCb measurement [1].

Data distributions of the Dalitz-plot
The event distribution for m 2 (π + π − ) versus m 2 (J/ψπ + ) in Fig. 5 shows clear structures in m 2 (π + π − ). The presence of possible exotic structures in the J/ψπ + system, as claimed in similar decays [24,25], is investigated by examining the J/ψπ + mass distribution shown in Fig. 6 (a). No resonant effects are evident. Figure 6 (b) shows the π + π − mass distribution.  Figure 6: Distributions of (a) m(J/ψπ + ) and (b) m(π + π − ) for B 0 s → J/ψπ + π − candidate decays within ±20 MeV of the B 0 s mass. The (red) points with error bars show the background contribution determined from m(J/ψπ + π − ) fits performed in each bin of the plotted variables. 8

Detection efficiency
The detection efficiency is determined from a phase space simulation sample containing 4 × 10 6 B 0 s → J/ψπ + π − events with J/ψ → µ + µ − . The efficiency can be parameterized in terms of analysis variables as where s 12 ≡ m 2 (J/ψ π + ) and s 13 ≡ m 2 (J/ψ π − ) are functions of (m hh , θ hh ); such parameter transformations in ε 1 are implemented in order to use the Dalitz-plot based efficiency model developed in previous publications [3,19]. The efficiency functions take into account correlations between m hh and each of the three angles as determined by the simulation. The efficiency as a function of the angle χ is shown in Fig. 7. To simplify the normalization of the PDF, the efficiency as a function of χ is parameterized in 26 bins of m 2 where The efficiency in cos θ J/ψ also depends on m hh ; we fit the cos θ J/ψ distributions of J/ψπ + π − simulation sample with the function giving 26 values of a as a function of m 2 hh . The resulting distribution in a is shown in Fig. 8 and is best described by a 2nd order polynomial function  The function ε 1 (s 12 , s 13 ) can be determined from the simulation after integrating over cos θ J/ψ and χ, because the functions ε 2 and ε 3 are normalized in cos θ J/ψ and χ, respectively. It is parameterized as a symmetric 5th order polynomial function given by where x = s 12 /GeV 2 − 18.9, and y = s 13 /GeV 2 − 18.9. The phase space simulation is generated uniformly in the two-dimensional distribution of (s 12 , s 13 ), therefore the distribution of selected events reflects the efficiency and is fit to determine the efficiency parameters ε i . The projections of the fit are shown in Fig. 9, giving the efficiency as a function of cos θ π + π − versus m(π + π − ) in Fig. 10.

Background composition
The main background source is combinatorial and is taken from the like-sign combinations within ±20 MeV of the B 0 s mass peak. The like-sign combinations also contain the B − background which is peaked at cos θ hh = ±1. The like-sign combinations cannot contain any ρ 0 , which is measured to be 3.5% of the total background. To obtain the ρ 0 contribution, the background m(π + π − ) distribution shown in Fig. 6 (b), found by fitting the m(J/ψ π + π − ) distribution in bins of m(π + π − ), is compared to m(π ± π ± ) distribution from the like-sign combinations. In this way simulated ρ 0 background is added into the like-sign candidates. The background PDF B is the sum of functions for B − (B B − ) and  Figure 10: Parameterization of the detection efficiency as a function of cos θ π + π − and m(π + π − ). The scale is arbitrary.
for the other (B other ), given by where N other and N B − are normalization factors, and f B − is the fraction of the B − background in the total background. The J/ψπ + π − mass fit gives f B − = (1.7 ± 0.2)%.
The B − background is separated because its invariant mass is very close to the highest allowed limit, resulting in its cos θ hh distribution peaking at ±1. The function for the B − background is defined as where G is the Gaussian function, and the parameters m 0 , σ m , σ θ , p b1 , and p b2 are determined by the fit. The last term is the same function for χ.
The function for the other background is where the function Here where m min and m max give the fit boundaries of m hh , B 2 (ζ) is a fifth-order Chebychev polynomial; q(ζ) and p(ζ) are both second-order Chebychev polynomials with the coefficients c 1 and c 2 being free parameters. In order to better approximate the real background in the B 0 s signal region, the J/ψ π ± π ± candidates are kinematically constrained to the B 0 s mass, and µ + µ − to the J/ψ mass. The second part 1 + α cos 2 θ J/ψ is a function of the J/ψ helicity angle. The cos θ J/ψ distribution of background is shown in Fig. 11; fitting with the function determines the parameter α = −0.34 ± 0.03. A fit to the like-sign combinations added with additional ρ 0 background determines the parameters describing the m hh , θ hh , and χ distributions. Figures 12 and 13 show the projections of cos θ hh and m hh , and of χ of the total background, respectively.

Resonance models
To study the resonant structures of the decay B 0 s → J/ψπ + π − we use the 34 471 candidates with invariant mass lying within ±20 MeV of the B 0 s mass peak which include 7075±101 background events. The π + π − resonance candidates that could contribute to B 0 s → J/ψπ + π − decay are listed in Table 2. The resonances that decay into a π + π − pair must be isoscalar (I = 0), because the ss system forming the resonances in Fig. 1 has I = 0. To test the isoscalar argument, the isospin-1 ρ(770) meson is also added to the baseline fit. The non-resonance (NR) is assumed to be S-wave, its shape is defined by Eq. (3) where the amplitude function A R (m hh ) is set to be equal to one, and the Blatt-Weisskopf In the previous analysis [22], we observed a resonant state at (1475 ± 6) MeV with a width of (113 ± 11) MeV. We identified it with the f 0 (1370) though its mass and width values agreed neither with the f 0 (1500) or the f 0 (1370). W. Ochs [26] argues that the better assignment is f 0 (1500); we follow his suggestion. In addition, a structure is clearly visible in the 1800 MeV region (see Fig. 6 (b)), which was not the case in our previous analysis [3]. This could be the f 0 (1790) resonance observed by BES [27] in J/ψ → φπ + π − decays.
From the measured ratios B B 0 , using the measured π + π − and K + K − branching fractions [6], the expected f 2 (1525) fit fraction for the transversity 0 component is (0.45 ± 0.13)%, and the ratio of helicity λ = 0 to |λ| = 1 components, which is equal to the ratio of transversity 0 to the sum of ⊥ and components, is 1.9 ± 0.8, where the uncertainties are dominated by that on f 2 (1525) fit fractions in B 0 s → J/ψ K + K − decays. This information is used as constraints in the fit.
The masses and widths of the resonances are also listed in Table 2. When used in the fit they are fixed to these central values, except for the parameters of f 0 (980) and f 0 (1500) that are determined by the fit. In addition, the parameters of f 0 (1790) are constrained to those determined by the BES measurement [27].
As suggested by D. V. Bugg [29], the Flatté model [30] for f 0 (980) is slightly modified, and is parameterized as where m R is the f 0 (980) pole mass, the parameters g ππ and g KK are the f 0 (980) coupling constants to π + π − and K + K − final states, respectively, and the phase space ρ factors are given by Lorentz-invariant phase spaces as Compared to the normal Flatté function, a form factor F KK = exp(−αk 2 ) is introduced above the KK threshold and serves to reduce the ρ KK factor as m 2 π + π − increases, where k is momentum of each kaon in the KK rest frame, and α = (2.0 ± 0.25) GeV −2 [29]. This parameterization slightly decreases the f 0 (980) width above the KK threshold. The parameter α is fixed to 2.0 GeV −2 as it is not very sensitive to the fit.
To determine the complex amplitudes in a specific model, the data are fitted maximizing the unbinned likelihood given as where N is the total number of candidates, and F is the total PDF defined in Eq. (7). In order to converge properly in a maximum likelihood method, the PDFs of the signal and background need to be normalized. This is accomplished by first normalizing the χ and cos θ J/ψ dependent parts analytically, and then normalizing the m hh and cos θ hh dependent parts using a numerical integration over 1000×200 bins. The fit determines amplitude magnitudes a R i i and phases φ R i i defined in Eq. (6). The a f 0 (980) 0 amplitude is fixed to 1, since the overall normalization is related to the signal yield. As only relative phases are physically meaningful, φ f 0 (980) 0 is fixed to 0. In addition, due to the averaging of B 0 s and B 0 s , the interference terms between opposite CP states are cancelled out, making it not possible to measure the relative phase between CP -even and odd states here, so one CP -even phase, φ f 2 (1270) ⊥ , is also fixed to 0.

Fit fraction
Knowledge of the contribution of each component can be expressed by defining a fit fraction for each transversity τ , F R τ , which is the squared amplitude of R integrated over the phase space divided by the entire amplitude squared over the same area. To determine F R τ one needs to integrate over all the four fitted observables in the analysis. The interference terms between different helicity components vanish, after integrating over the two variables of cos θ J/ψ and χ. Thus we define the transversity fit fraction as 15 where λ = 0 in the d-function for τ = 0, and λ = 1 for τ =⊥ or . Note that the sum of the fit fractions is not necessarily unity due to the potential presence of interference between two resonances. Interference term fractions are given by and Interference between different spin-J states vanishes, when integrated over angle, because the d J λ0 angular functions are orthogonal.

Fit results
In order to compare the different models quantitatively, an estimate of the goodness of fit is calculated from four-dimensional (4D) partitions of the four variables, m(π + π − ), cos θ hh , cos θ J/ψ and χ. We use the Poisson likelihood χ 2 [31] defined as where n i is the number of events in the four-dimensional bin i and x i is the expected number of events in that bin according to the fitted likelihood function. A total of 1845 bins are used to calculate the χ 2 , where 41(m hh ) × 5(cos θ hh ) × 3(cos θ J/ψ ) × 3(χ) equal size bins are used, and m hh is required to be between 0.25 and 2.30 GeV. The χ 2 /ndf, and the negative of the logarithm of the likelihood, −lnL, of the fits are given in Table 3, where ndf is the number of degree of freedom given as 1845 subtracted by number of fitting parameters and 1. The nomenclature describing the models gives the base model first and then "+" for any additions. The 5R model contains the resonances f 0 (980), f 2 (1270), f 2 (1525), f 0 (1500), and f 0 (1790). If adding NR to 5R model, two minima with similar likelihoods are found. One minimum is consistent with the 5R results and has NR fit fraction of (0.3 ± 0.3)%; we group any fit models that are consistent with this 5R fit into the "Solution I" category. Another minimum has significant NR fit fraction of (5.9 ± 1.4)%, this model and other consistent models are classified in the "Solution II" category. Among these resonance models, we select the baseline model by requiring each resonance in the model to have more than 3 standard deviation (σ) significance evaluated by the fit fraction divided by its uncertainty. The baseline fits are 5R in Solution I and 5R+NR in Solution II. No additional components are significant when added to these baseline fits. Unfortunately, we cannot distinguish between these two solutions and will     Table 4 shows the fit fractions from the baseline fits of two solutions, where systematic uncertainties are included; they will be discussed in Sec. 7. Figures 14 and 15 show the fit projections of m(π + π − ), cos θ ππ , cos θ J/ψ and χ from 5R Solution I and 5R+NR Solution II, respectively. Also shown in Figs. 16 and 17 are the contributions of each resonance as a function of m(π + π − ) from the baseline Solution I and II fits, respectively. Table 5 shows the fit fractions of the interference terms defined in Eq. (23). In addition, the phases are listed in Table 6. The other fit results are listed in Table 7 including the f 0 (980) mass, the Flatté function parameters g ππ , g KK /g ππ , and masses and widths of f 0 (1500) and f 0 (1790) resonances.
In both solutions the f 0 (500) state does not have a significant fit fraction. We set an upper limit for the fit fraction ratio between f 0 (500) and f 0 (980) of 0.3% from Solution I and 3.4% from Solution II, both at 90% CL. A similar situation is found for the ρ(770) state. When including it in the fit, the fit fraction of ρ(770) is measured to be (0.60 ± 0.30 +0.08 −0.14 )% in Solution I and (1.02 ± 0.36 +0.09 −0.15 )% from Solution II. The largest upper limit is obtained by Solution II, where the ρ(770) fit fraction is less than 1.7% at 90% CL.
Our previous study [3] did not consider the f 0 (1790) resonance, instead the NR component filled in the higher mass region near 1800 MeV. It is found that including  Figure 18 compares the total S-wave amplitude strength and phase as a function of m(π + π − ) between the two solutions, showing consistent amplitude strength but distinct phase. The total S-wave amplitude is calculated as Eq. (4) summing over all spin-0 component R with λ = 0, where the d-function is equal to 1. The amplitude strength can be well measured from the m(π + π − ) distribution, but this is not the case for the phase, which is determined from the interference with the small fraction of higher spin resonances.

Angular moments
We define the moments of the cosine of the helicity angle θ ππ , Y 0 l (cos θ ππ ) as the efficiency corrected and background subtracted π + π − invariant mass distributions, weighted by spherical harmonic functions. The moment distributions provide an additional way of visualizing the presence of different resonances and their interferences, similar to a partial wave analysis. Figures 19 and 20 show the distributions of the angular moments for 5R Solution I and 5R+NR Solution II, respectively. In general the interpretation of these

Systematic uncertainties
The sources of the systematic uncertainties on the results of the amplitude analysis are summarized in Table 8 for Solution I and Table 9 for Solution II. The contributions to the systematic error due to φ s , the function ε(t), Γ s and ∆Γ s [6] uncertainties, and L B choices for transversity 0 and of spin ≥ 1 resonances, are negligible. The systematic errors associated to the acceptance or background modeling are estimated by repeating the fit to the data 100 times. In each fit the parameters in the acceptance or background function are randomly generated according to the corresponding error matrix. The uncertainties due to the fit model include possible contributions from each resonance listed in Table  2 but not used in the baseline fit models, varying the hadron scale r parameters in the Blatt-Weisskopf barrier factors for the B meson and R resonance from 5.0 GeV −1 and 1.5 GeV −1 , respectively, to both 3.0 GeV −1 , and using F KK = 1 in the Flatté function. among those changes is assigned as the systematic uncertainties for modeling.
Finally, we repeat the data fit by varying the mass and width of resonances within their errors one at a time, and add the changes in quadrature. To assign a systematic uncertainty from the possible presence of the f 0 (500) or ρ(770), we repeat the above procedures using the model that has the baseline resonances plus f 0 (500) or ρ(770).

Fit fraction intervals
The fit fractions shown in Table 4 differ considerably for some of the states between the two solutions. Table 10 lists the 1σ regions for the fit fractions taking into account the differences between the solutions and including systematic uncertainties. The regions covers both 1σ intervals of the two solutions.

CP content
The only CP -even content arises from the ⊥ projections of the f 2 (1270) and f 2 (1525) resonances, in addition to the 0 and of any possible ρ(770) resonance. The CP -even measured values are (0.89 ± 0.38 +0.59 −0.10 )% and (0.86 ± 0.42 +0.66 −0.10 )% for Solutions I and II, respectively (see Table 4), where the systematic uncertainty is dominated by the forbidden ρ(770) transversity 0 and components added in quadrature. To obtain the corresponding upper limit, the covariance matrix and parameter values from the fit are used to generate 2000 sample parameter sets. For each set, the CP -even fraction is calculated and is then smeared by the systematic uncertainty. The integral of 95% of the area of the distribution yields an upper limit on the CP -even component of 2.3% at 95% CL, where the larger value given by Solution II is used. The upper limit is the same as our previous measurement [3], while the current measurement also adds in a possible f 2 (1525) contribution.

Mixing angle and interpretation of light scalars
The I = 0 resonanances, f 0 (500) and f 0 (980), are thought to be mixtures of underlying states whose mixing angle has been estimated previously (see references cited in Ref. [32]). The mixing is parameterized by a normal 2×2 rotation matrix characterized by the angle