Measurement of resonant and CP components in B̄s → J=ψπþπ− decays

Structure of the decay B̄s → J=ψπþπ− is studied using data corresponding to 3 fb−1 of integrated luminosity from pp collisions produced by the LHC and collected by the LHCb detector. Five interfering πþπ− states are required to describe the decay: f0ð980Þ; f0ð1500Þ; f0ð1790Þ; f2ð1270Þ, and f2ð1525Þ. An alternative model including these states and a nonresonant J=ψπþπ− 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 a 95% confidence level. The f0ð500Þ state is not observed, allowing a limit to be set on the absolute value of the mixing angle with the f0ð980Þ of < 7.7° at a 90% confidence level, consistent with a tetraquark interpretation of the f0ð980Þ substructure.


I. INTRODUCTION
CP violation studies in theB 0 s → J=ψπ þ π − decay mode complement studies usingB 0 s → J=ψϕ and improve the final accuracy in the measurement of the CP-violating phase, ϕ s [1]. While the CP content was previously shown to be more than 97.7% CP odd at a 95% confidence level (C.L.), 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 we 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 theB 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 nonresonant amplitudes.
II. AMPLITUDE FORMULA FORB 0 s → J=ψh þ h − The decay ofB 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 theB 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 theB 0 s rest frame; and (iii) χ, the angle between the J=ψ and h þ h − decay planes in theB 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 time-integrated and flavor-averaged decay rate is proportional to the function Sðm hh ; θ hh ; θ J=ψ ; χÞ ¼ jAðm hh ; θ hh ; θ J=ψ ; χÞj 2 þ jĀðm hh ; θ hh ; θ J=ψ ; χÞj 2 − 2DRe ( q p A Ã ðm hh ; θ hh ; θ J=ψ ; χÞ ×Āðm hh ; θ hh ; θ J=ψ ; χÞ ) ; 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 nonresonant) 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 theB 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 AE 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 theB resonance decay properties to form the expression Here P B is the J=ψ momentum in theB 0 s rest frame, P R is the momentum of either of the two hadrons in the dihadron rest frame, m B is theB 0 s mass, J R is the spin of R, L B is the orbital angular momentum between the J=ψ and h þ h − system, and L R is the orbital angular momentum in the h þ h − decay, and thus is the same as the spin of the h þ h − resonance. F ðL B Þ B and F ðL R Þ R are the Blatt-Weisskopf barrier factors for theB 0 s and R resonance, respectively [3]. The factor ffiffiffiffiffiffiffiffiffiffiffi ffi P R P B p 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 d functions are listed in Ref. [6].
The decay rates, j A 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 Here a R Assuming no direct CP violation, as this has not been observed inB 0 s → J=ψϕ decays [1], the relation between thē B 0 s and B 0 s variables isā R τ ¼ η R τ a R τ , where η R τ is the CP eigenvalueoftheτ transversity componentfor theintermediate state R, where τ denotes the 0, ∥, or ⊥ component. The finalstate CP parities for S, P, and D waves are given in Table I for each resonance R and each transversity τ. For the τ ¼ ⊥ amplitude, the L B value of a spin-1 (or spin-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 spin-2) resonances. In this analysis, the lower one is used. It is verified that our results are insensitive to the L B choices.

III. 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 siliconstrip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes [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 an 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 aB 0 s meson by the use of either 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].

IV. EVENT SELECTION
Preselection criteria are implemented to preserve a large fraction of the signal events and are identical to those used in Ref. [19]. AB 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 theB 0 s candidate is required to have the track fit χ 2 =n.d.f. to be less than 4, where n.d.f. 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 a dimuon invariant mass between −48 MeV and þ43 MeV from the observed J=ψ mass peak are selected, and they are then constrained to the J=ψ mass [6] for subsequent use.
Pion candidates are required to each have a p T greater than 250 MeV, and the sum, p T ðπ þ Þ þ p T ðπ − Þ, must be larger than 900 MeV. Both pions must have χ 2 IP greater than 9 to reject particles produced from the PV. (The reconstruction procedure and the PV resolution are given in Ref. [20].) The χ 2 IP is computed as the difference between the χ 2 's of the PV reconstructed with and without the considered track. Both pions must also come from a common vertex with χ 2 =n:d:f: < 16 and form a vertex with the J=ψ with a χ 2 =n.d.f. less than 10 (here n.d.f. equals 5). 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.
TheB 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 [21]. 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 theB 0 s properties of vertex χ 2 , pointing angle, flight distance, p T , and χ 2 IP . The BDT is trained on a simulated sample ofB 0 s → J=ψπ þ π − signal events and a background data sample from the sideband  3 We work in units where c ¼ 1.
5566 < mðJ=ψπ þ π − Þ < 5616 MeV. Then the BDT is tested on independent samples. The distributions of the BDT classifier for signal and background samples are shown in Fig. 3. By maximizing the signal significance, we set the requirement that the classifier be 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 theB 0 s mass and a smaller one at theB 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 knownB 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 contributions from where the K − in the former, and both K − and p in the latter, are misidentified as pions. The shape of theB 0 → J=ψπ þ π − signal is taken to be the same as that of theB 0 s . The combinatorial background shape is taken from like-sign combinations that are the sum of π þ π þ and π − π − candidates, and it was found to be well described by an exponential function in previous studies [3,22]. 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 AE 201, which is obtained from study of the events in the control region of 5066 < mðJ=ψπ þ π − Þ < 5141 MeV. The mass fit gives 27396 AE 207 signal and 7075 AE 101 background candidates, leading to the signal fraction f sig ¼ ð79.5 AE 0.2Þ%, within AE20 MeV of theB 0 s mass peak. The effective rms mass resolution is 9.9 MeV.

V. PROBABILITY DENSITY FUNCTION CONSTRUCTION
The correlated distributions of four variables m hh , cos θ hh , cos θ J=ψ , and χ are fitted using the candidates within AE20 MeV of theB 0 s mass peak. To improve the resolution of these variables, we perform a kinematic fit constraining theB 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 Fðm hh ;θ hh ; θ J=ψ ;χÞ ¼ f sig N sig εðm hh ;θ hh ; θ J=ψ ;χÞ × Sðm hh ; θ hh ; θ J=ψ ;χÞ þ ð1 − f sig ÞBðm hh ;θ hh ;θ J=ψ ;χÞ; (7) where ε is the detection efficiency, and B is the background PDF discussed later in Sec. V C. The normalization factor for the signal is given by The signal function S is defined in Eq. (1), where D ¼ ð8.7 AE 1.5Þ%, taking into account the acceptance [23], and choosing a phase convention q=p ¼ e −iϕ s . The phase ϕ s is fixed to the standard model value of −0.04 The data have been fitted with a double Crystal Ball signal and several background functions. The (red) solid curve shows theB 0 s signal, the (brown) dotted line shows the combinatorial background, the (green) short-dashed line shows the B − background, the (purple) dot-dashed curve isB 0 → J=ψπ þ π − , the (light blue) long-dashed line is the sum ofB 0 s → J=ψη 0 ,B 0 s → J=ψϕ with ϕ → π þ π − π 0 backgrounds and the Λ 0 b → J=ψK − p reflection, the (black) dot-long dashed curve is theB 0 → J=ψK − π þ reflection, and the (blue) solid curve is the total. radians [24]. Our results are found to be insensitive to the value of ϕ s used within the 95% C.L. limits set by the LHCb measurement [1].

A. 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 [25,26], 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. Apart from a large signal peak due to the f 0 ð980Þ, there are visible structures at about 1450 MeV and 1800 MeV.
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 parametrized in 26 bins of m 2 hh as where The efficiency in cos θ J=ψ also depends on m hh ; we fit the cos θ J=ψ distributions of the 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 second-order polynomial function FIG. 6 (color online). Distributions of (a) mðJ=ψπ þ Þ and (b) mðπ þ π − Þ forB 0 s → J=ψπ þ π − candidate decays within AE20 MeV of thē 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.
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 parametrized as a symmetric fifth-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.

C. Background composition
The main background source is combinatorial and is taken from the like-sign combinations within AE20 MeV of theB 0 s mass peak. The like-sign combinations also contain the B − background, which is peaked at cos θ hh ¼ AE1. 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 the mðπ AE π AE Þ 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 for the other (B other ), given by 14) 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 AE 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 AE1. The function for the B − background is defined as þ p b1 cos χ þ p b2 cos 2χÞ; (15) 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; and 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 theB 0 s signal region, the J=ψπ AE π AE candidates are kinematically constrained to theB 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 the background is shown in Fig. 11; fitting with the function determines the parameter α ¼ −0.34 AE 0.03. A fit to the like-sign combinations added with additional ρ 0 background determines the parameters describing the m hh , θ hh , and χ distributions.

A. Resonance models
To study the resonant structures of the decaȳ B 0 s → J=ψπ þ π − , we use the 34 471 candidates with invariant mass lying within AE20 MeV of theB 0 s mass peak, which include 7075 AE 101 background events. The π þ π − resonance candidates that could contribute toB 0 s → J=ψπ þ π − decay are listed in Table II. 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 nonresonance (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 1, and the Blatt-Weisskopf barrier factors F ð1Þ B and F ð0Þ R are both set to 1.
In the previous analysis [23], we observed a resonant state at ð1475 AE 6Þ MeV with a width of ð113 AE 11Þ MeV. We identified it with the f 0 ð1370Þ, though its mass and width values agreed neither with the f 0 ð1500Þ nor with the f 0 ð1370Þ. W. Ochs [29] 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 [28] in J=ψ → ϕπ þ π − decays.
The masses and widths of the resonances are also listed in Table II. 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 [28].
As suggested by D. V. Bugg [30], the Flatté model [31] for f 0 ð980Þ is slightly modified and is parametrized 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 the π þ π − 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 AE 0.25Þ GeV −2 [30]. This parametrization 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 andB 0 s , the interference terms between opposite CP states are canceled out, making it not possible to measure the relative phase between CP-even and CP-odd states here, so one CP-even phase, ϕ f 2 ð1270Þ ⊥ , is also fixed to 0.

B. 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 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

5R (Solution I)
Interference between different spin-J states vanishes when integrated over angle, because the d J λ0 angular functions are orthogonal.

C. Fit results
In order to compare the different models quantitatively, an estimate of the goodness of fit is calculated from fourdimensional (4D) partitions of the four variables, mðπ þ π − Þ, cos θ hh , cos θ J=ψ , and χ. We use the Poisson likelihood χ 2 [32], 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 =n:d:f. and the negative of the logarithm of the likelihood, − ln L, of the fits are given in Table III, where n.d.f. is the number of degrees of freedom, given as 1845, subtracted by the 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 0 2 ð1525Þ, f 0 ð1500Þ, and f 0 ð1790Þ. In adding NR to the 5R model, two minima with similar likelihoods are found. One minimum is consistent with the 5R results and has a NR fit fraction of ð0.3 AE 0.3Þ%; we group any fit models that are consistent with this 5R fit into the "Solution I" category. Another minimum has a significant NR fit fraction of ð5.9 AE 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 deviations (σ) of 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 quote results for both of them. In both cases, the dominant contribution is S wave, including f 0 ð980Þ, f 0 ð1500Þ, and f 0 ð1790Þ. The D wave, f 2 ð1270Þ and f 0 2 ð1525Þ, is only 2.3% for both solutions. Table IV shows the fit fractions from the baseline fits of two solutions, where systematic uncertainties are included; they will be discussed in Sec. VII. 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 V shows the fit fractions of the interference terms defined in Eq. (23). In addition, the phases are listed in Table VI. The other fit results are listed in Table VII, 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 a 90% C.L. 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 AE 0.30 þ0.08 −0.14 Þ% in Solution I and ð1.02 AE 0.36 þ0.09 −0.15 Þ% in Solution II. The largest upper limit is obtained by Solution II, where the ρð770Þ fit fraction is less than 1.7% at 90% C.L. Our previous study [3] did not consider the f 0 ð1790Þ resonance; instead, the NR component filled in the highermass region near 1800 MeV. It is found that including f 0 ð1790Þ improves the fit significantly in both solutions. Inclusion of this state reduces −2 ln L by 276 (97) 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.

D. Angular moments
We define the moments of the cosine of the helicity angle θ ππ , hY 0 l ð cos θ ππ Þi, 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 moments [3] is that hY 0 0 i is the efficiency-corrected and background-subtracted event distribution, hY 0 1 i is the interference of the sum of S-wave and P-wave and P-wave and D-wave amplitudes, hY 0 2 i is the sum of the P-wave, D-wave and the interference of S-wave and D-wave amplitudes, hY 0 3 i is the interference between the P-wave and D-wave amplitudes, hY 0 4 i is D wave, and

Components
Solution I Solution II f 0 ð980Þ þ f 0 ð1500Þ 9.50 −1.57 The points with error bars are the data points, and the solid curves are derived from the model 5R Solution I.

VII. SYSTEMATIC UNCERTAINTIES
The sources of the systematic uncertainties on the results of the amplitude analysis are summarized in Table VIII for  Solution I and Table IX 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 with 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 II but not used in the baseline fit models, varying the hadron scale r parameters in the Blatt-Weisskopf barrier factors for both the B meson and R resonance from 5.0 GeV −1 and 1.5 GeV −1 , respectively, to 3.0 GeV −1 , and using F KK ¼ 1 in the Flatté function. Compared to the nominal Flatté function, the new one improves the likelihood fit −2 ln L by 6.8 and 14.0 units for Solution I and Solution II, respectively. The largest variation among those changes is assigned as the systematic uncertainty for modeling.
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Þ.
Finally, we have tested the entire procedure with simulated pseudoexperiments producing both the signal and The points with error bars are the data points, and the solid curves are derived from the model 5R þ NR Solution II.
backgrounds and have verified that the fit finds the correct resonant substructure with the correct uncertainties.

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

B. CP content
The only CP-even content arises from the ⊥ projections of the f 2 ð1270Þ and f 0 2 ð1525Þ resonances, in addition to the 0 and ∥ of any possible ρð770Þ resonance. The CP-even measured values are ð0.89 AE 0.38 þ0.59 −0.10 Þ% and ð0.86 AE 0.42 þ0.66 −0.10 Þ% for Solutions I and II, respectively (see Table IV), 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 a 95% C.L., 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 0 2 ð1525Þ contribution.

C. 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. [33]). The mixing is parametrized by a normal 2 × 2 rotation matrix characterized by the angle φ m , giving in our case where jnni ≡ 1 ffiffi ffi 2 p ðjuūi þ jddiÞ: In this case, only the jssi wave function contributes. Thus, we have [2] Φð980Þ Φð500Þ ; where the Φ's are phase-space factors. The phase space in this pseudoscalar-to-vector-pseudoscalar decay is proportional to the cube of the f 0 momenta. Taking the average of the momentum-dependent phase space over the resonant line shapes results in the ratio of phase-space factors Φð500Þ Φð980Þ ¼ 1.25.
This limit is the most constraining ever placed on this mixing angle [19]. The value of tan 2 φ m is consistent with the tetraquark model, which predicts zero within a few degrees [2,33].

IX. CONCLUSIONS
TheB 0 s → J=ψπ þ π − decay can be described by the interfering sum of five resonant components: f 0 ð980Þ, f 0 ð1500Þ, f 0 ð1790Þ, f 2 ð1270Þ, and f 0 2 ð1525Þ. In addition, we find that a second model including these states plus nonresonant J=ψπ þ π − also provides a good description of the data. In both models, the largest component of the decay is the f 0 ð980Þ, with the f 0 ð1500Þ being almost an order of magnitude smaller. We also find that including the f 0 ð1790Þ resonance improves the data fit significantly. The π þ π − system is mostly S wave, with the D-wave components totaling only 2.3% in either model. No significant B 0 s → J=ψρð770Þ decay is observed; a 90% C.L. upper limit on the fit fraction is set to be 1.7%.
The most important result of this analysis is that the CP content is consistent with being purely odd, with the CPeven component limited to 2.3% at 95% C.L. Also of importance is the limit on the absolute value of the mixing angle between the f 0 ð500Þ and f 0 ð980Þ resonances of 7.7°a t 90% C.L., the most stringent limit ever reported. This is also consistent with these states being tetraquarks.

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