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

The resonant structure of the reaction $\overline{B}^0\rightarrow J/\psi \pi^+\pi^-$ is studied using data from 3 fb$^{-1}$ of integrated luminosity collected by the LHCb experiment, one-third at 7 Tev center-of-mass energy and the remainder at 8 Tev. The invariant mass of the $\pi^+\pi^-$ pair and three decay angular distributions are used to determine the fractions of the resonant and non-resonant components. Six interfering $\pi^+\pi^-$ states: $\rho(770)$, $f_0(500)$, $f_2(1270)$, $\rho(1450)$, $\omega(782)$ and $\rho(1700)$ are required to give a good description of invariant mass spectra and decay angular distributions. The positive and negative CP fractions of each of the resonant final states are determined. The $f_0(980)$ meson is not seen and the upper limit on its presence, compared with the observed $f_0(500)$ rate, is inconsistent with a model of tetraquark substructure for these scalar mesons at the eight standard deviation level. In the $q\overline{q}$ model, the absolute value of the mixing angle between the $f_0(980)$ and the $f_0(500)$ scalar mesons is limited to be less than $17^{\circ}$ at 90% confidence level.


Introduction
The decay mode B 0 → J/ψπ + π − is of particular interest in the study of CP violation in the B system. 1 The decay can proceed either via a tree level process, shown in Fig. 1(a), or via the penguin mechanisms shown in Fig. 1(b). The ratio of penguin to tree amplitudes is enhanced in this decay relative to B 0 → J/ψ K 0 S [1]. Thus the effects of penguin topologies can be investigated by using the J/ψπ + π − decay and comparing different measurements of the CP violating phase, β, in J/ψ K 0 S , and individual channels such as B 0 → J/ψρ 0 . The B 0 → J/ψπ + π − decay is also useful for the study of the substructure of light mesons that decay into π + π − . Tests have been proposed to ascertain if the scalar f 0 (500) and f 0 (980) mesons are formed of qq or tetraquarks. In the model of Ref. [2], if these scalar states are tetraquarks, the ratio of decay widths is predicted to be 1/2. If instead these are qq states, they can be mixtures of two base states; in this scenario the width ratio can be any value and is determined principally by the mixing angle between the base states. The B 0 → J/ψπ + π − decay was first observed by the BaBar collaboration [?]. It has been previously studied by LHCb using data from 1 fb −1 of integrated luminosity [3]. The branching fraction was measured to be (3.97 ± 0.22) × 10 −5 . The mass and angular distributions were used to measure the resonant substructure. That analysis, however, did not use the angle between the J/ψ and π + π − decay planes, due to limited statistics.
A new theoretical approach [4] now allows us to include all the angular information and measure the fraction of CP -even and CP -odd states. This information is vital to any subsequent CP violation measurements.

Data sample and detector
In this paper, we measure the resonant substructure and CP content of the B 0 → J/ψπ + π − decay from data corresponding to 3 fb −1 of integrated luminosity collected with the LHCb detector [5] using pp collisions. One-third of the data was acquired at a center-of-mass energy of 7 TeV, and the remainder at 8 TeV. The detector is a singlearm 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 [6] placed downstream. The combined tracking system provides a momentum measurement with relative uncertainty that varies from 0.4% at 5 GeV to 0.6% at 100 GeV, 2 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 (RICH) detectors [7]. 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 [8].
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 [9]. Events selected for this analysis are triggered by a J/ψ → µ + µ − decay, where the J/ψ meson is required at the software level to be consistent with coming from the decay of a B 0 meson by use either of IP requirements or detachment of the J/ψ meson decay vertex from the primary vertex (PV). In the simulation, pp collisions are generated using Pythia [10] with a specific LHCb configuration [11]. Decays of hadronic particles are described by EvtGen [12], in which final state radiation is generated using Photos [13]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [14,15] as described in Ref. [16].

Observables used in the analysis
The B 0 → J/ψπ + π − decay with J/ψ → µ + µ − can be described by the invariant mass of the π + π − (m hh ) pair, and three angles: (i) the angle between the µ + direction in the J/ψ rest frame with respect to the J/ψ direction in the B 0 rest frame, θ J/ψ ; (ii) the angle between the π + direction and the opposite direction of the B 0 candidate momentum in the π + π − rest frame, θ hh ; and (iii) the angle between the J/ψ and π + π − decay planes in the B 0 rest frame, χ. The angular variables are illustrated in Fig. 2.
In our previous study [3], we used the "Dalitz-plot" variables: the invariant mass squared of J/ψ π + , s 12 = m 2 (J/ψ π + ), and the invariant mass squared of the π + π − pair, s 23 = m 2 (π + π − ). Due to the J/ψ spin, the event distributions in the s 12 and s 23 plane do not directly show the effect of the matrix-element squared. Since the probability density functions (PDFs) expressed as functions of m hh and θ hh are easier to normalize, we use them instead. In this paper, the notation hh is equivalent to π + π − . The Dalitz-plot variables can be translated into (m hh , θ hh ), and vice versa. The formalism described below is for the decay sequence B 0 → J/ψ R, R → π + π − .

Amplitude formalism
The decay rate of ( -) B 0 → J/ψ π + π − has been described in detail in Ref. [4]. The differential decay width can be written in terms of the decay time t and the four other variables m hh , θ J/ψ , θ hh , and χ as [17] where N is a constant; ( -) A is the amplitude of ( -) B 0 → J/ψ π + π − at the decay time t = 0, which is itself a function of m hh , θ J/ψ , θ hh , and χ, summed over all resonant (and possibly non-resonant) components; ∆m is the mass difference between the heavy and light B 0 mass eigenstates, and ∆Γ the width difference; 3 q and p are complex parameters that describe the relation between mass and flavor eigenstates. In this analysis we take |p/q| to be equal to unity.
Forming the sum of B 0 and B 0 decay widths and integrating over decay time, yields the time-integrated and flavor-averaged decay width where we drop the term arising from quantum interference of the amplitudes in the last line. This results from the fact that the D factor is negligibly small for B 0 meson decays. Specifically, where α(t) is the decay time dependent detection efficiency. 4 Since ∆Γ/Γ is of the order of 1% for B 0 meson decays [18], the D term is about the same size.
We define A R (m hh ) to be the mass line shape of the resonance R, which in most cases is a Breit-Wigner function. It is combined with the decay properties of the B 0 and resonance to form the expression for the decay amplitude. For each resonance R: Here P R (P B ) is the scalar momentum of one of the two daughters of the resonance R (or the B 0 meson) in the R (or B 0 ) rest frame, J R is the spin of R, 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 h + h − resonance. F and F (L R ) R are the centrifugal barrier factors for the B 0 and the R resonance, respectively [19].
The factor √ P R P B results from converting the phase space of the Dalitz-plot variables m 2 hh and m 2 J/ψ h + to that of m hh and cos θ hh . The function defined in Eq. (5) is based on previous amplitude analyses [19,20].
We must sum over all final states, R, so for each J/ψ helicity, denoted by λ, equal to 0, +1 and −1 we have the overall decay amplitudes: where the Wigner-d functions are defined in Ref. [18] and ( -) h R λ are complex helicity coefficients. We note that the λ value of the J/ψ is equal to that of the R resonance. Finally, the total decay rate of ( -) B 0 → J/ψ π + π − at t = 0 is given by In order to determine the CP components, it is convenient to replace the complex helicity coefficients ( -) h R λ by the complex transversity coefficients Here ( -) a R 0 corresponds to the 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 the absence of direct CP violation, the relation between the B 0 and B 0 transversity coefficients isā R τ = η R τ a R τ , where η R τ is the CP eigenvalue of the τ transversity component for the intermediate state R, and τ denotes the 0, , or ⊥ components. Note that for the h + h − system both C and P are given by (−1) L R , so the CP of the h + h − system is always even. The total CP of the final state is (−1) L B , since the CP of the J/ψ is also even. The final state CP parities, for S, P, and D-waves, are listed in Table 1.
In this analysis a fit determines the amplitude modulus a R τ and the phase φ R τ of the amplitude for each resonance R, and each transversity component τ . For the τ =⊥ amplitude, the L B value of spin-1 (or spin-2) resonances is 1 (or 2). While the other transversity components, τ = 0 or , have two possible L B values of 0 and 2 (or 1 and 3) for spin-1 (or -2) resonances. We use only the smaller values for each. Studies show that our results for fractions of different interfering components are not sensitive to these L B choices.

Dalitz fit fractions
A complete description of the decay is given in terms of the fitted complex amplitudes.
Knowledge of the contribution of each component can be summarized by defining a fit To determine F R τ one needs to integrate over all the four variables: m hh , θ hh , θ J/ψ , χ. The interference terms between different helicity components vanish after integrating Eq. (7) over the two variables of cos θ J/ψ and χ, i.e.
The decay rate is the sum of the contributions from the three helicity terms. To define the transversity fractions, we need to write Eq. (10) in terms of transversity amplitudes. Since d J R −1,0 = −d J R 1,0 , the sum of the three helicity terms is equal to the sum of three transversities, given as Thus, we define the transversity fit fractions as where λ = 0 for τ = 0, and λ = 1 for τ =⊥ or . 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.

Selection requirements
In this analysis we adopt a two step selection. The first step, preselection, is followed by a multivariate selection based on a boosted decision tree (BDT) [21]. Preselection criteria are implemented to preserve a large fraction of the signal events, yet reject easily eliminated backgrounds, and are identical to those used in Ref. [3]. A B 0 → 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 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 [18] for subsequent use. Each pion candidate is required to have p T greater than 250 MeV, and that the scalar sum, p T (π + ) + p T (π − ) is required to be 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 the pion selection we require DLL(π − K) > −10 and DLL(π − µ) > −10. The B 0 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 • . 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 properties of vertex χ 2 , pointing angle, flight distance, p T and χ 2 IP . The BDT is trained on a simulated sample of two million B 0 → J/ψ π + π − signal events generated uniformly in phase space with unpolarized J/ψ → µ + µ − decays, and a background data sample from the sideband 5566 < m(J/ψπ + π − ) < 5616 MeV. Then the BDT is tested on independent samples from the same sources. The BDT can take any value from -1 to 1. The distributions of signal and background are approximately Gaussian shaped with r.m.s. of about 0.13. Signal peaks at BDT of 0.27 and background at -0.22. To minimize possible bias on the signal acceptance due to the BDT, we choose a loose requirement of BDT> 0, which has about a 95% signal efficiency and a 90% background rejection rate.

Fit model
We first select events based on their J/ψ π + π − invariant mass and then perform a full fit to the decay variables. The invariant mass of the selected J/ψπ + π − combinations is shown in Fig. 3. There is a large peak at the B 0 s mass and a smaller one at the B 0 mass on top of the background. A double Crystal Ball function with common means models the radiative tails and is used to fit each of the signals [22]. The known B 0 s − B 0 mass difference [18] is used to constrain the difference in mean values. Other components in the fit model take into account background contributions from B − → J/ψK − and , and combinatorial backgrounds. The exponential combinatorial background shape is taken from like-sign combinations, that are the sum of π + π + and π − π − candidates, and found to be a good description in previous studies [19,23].
The shapes of the other components are taken from the Monte Carlo simulation with their normalizations allowed to vary. The mass fit gives 18 841±204 signal and 10 207±178 background candidates within ±20 MeV of the B 0 mass peak. Only candidates within ±20 MeV of the B 0 mass peak are retained for further analysis. To improve the resolution of the mass and angular variables used in the amplitude analysis, we perform a kinematic fit constraining the B 0 and J/ψ masses to their PDG mass values [18], and recompute the final state momenta [24]. One of the main challenges in performing a mass and angular analysis is to construct a realistic probability density function, where both the kinematic and dynamical properties are modeled accurately. The PDF is given by the sum of signal, S, and background, B, functions. The B 0 signal includes events from the reaction B 0 → J/ψ K 0 S . These are described by a separate term in the PDF. The total PDF is where f sig is the fraction of the signal in the fitted region (f sig = (64.9 ± 1.2)% obtained from the mass fit in Fig. 3), ε is the detection efficiency described in Sec. 5.1, and B is the background PDF described in Sec. 5.2. The K 0 S component is modeled by a Gaussian function, G, with mean m K 0 S and width σ K 0 S . The Gaussian parameters together with the K 0 S fraction in the B 0 peak, f K 0 S , are determined in the fit. The normalization factors N sig for the signal and N K 0 S for the K 0 S candidates are efficiency-multiplied theoretical functions integrated over the four analysis variables, m hh , θ hh , θ J/ψ , and χ, given by Examination of the event distribution for m 2 (π + π − ) versus m 2 (J/ψ π + ) in Fig. 4 shows obvious structures in m 2 (π + π − ). To investigate if there are visible exotic structures in m 2 (J/ψ π + ), we examine the J/ψ π + invariant mass distribution as shown in Fig. 5(a) where we fit the m(J/ψ π + π − ) distribution to extract the background levels in bins of m(J/ψ π + ) (red points). Similarly, Fig. 5(b) shows the π + π − mass distribution. Apart reflection and the (blue) solid curve is the total. The points at the bottom show the difference between the data points and the total fit divided by the statistical uncertainty on the data. Figure 4: Distribution of m 2 (π + π − ) versus m 2 (J/ψ π + ) for all events within ±20 MeV of the B 0 mass.

Detection efficiency
The detection efficiency is determined from a sample of about four million simulated B 0 → J/ψπ + π − events that are generated uniformly in phase space with unpolarized The efficiency model can be expressed 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 dependence on χ is modeled by where The free parameters are determined by fitting the simulated χ distributions using Eq. (18) The fitted values of a are modeled by a second order polynomial function These variables are related to s 23 by Thus, ε 1 (s 12 , s 13 ) can be modeled by a two-dimensional fifth order polynomial function as ε 1 (s 12 , s 13 ) = 1 + 1 (x + y) + 2 (x + y) 2 + 3 xy + 4 (x + y) 3 + 5 xy(x + y) + 6 (x + y) 4 + 7 xy(x + y) 2 + 8 x 2 y 2 + 9 (x + y) 5 + 10 xy(x + y) 3 + 11 x 2 y 2 (x + y) where all the i are the fit parameters. The χ 2 /ndf is 313/299. The values of the parameters are given in Table 2.
The projections of the fit used to measure the efficiency parameters are shown in Fig. 6. The efficiency shapes are well described by the parametrization. The parameterized efficiency as a function of m(π + π − ) versus cos θ π + π − is shown in Fig. 7.

Background composition
The main background source in the B 0 signal region is combinatorial and can be taken from the like-sign combinations within ±20 MeV of the B 0 mass peak. In addition, there is background arising from partially reconstructed B π + π − γ and B 0 s → J/ψ φ with φ → π + π − π 0 ), reflections from misidentified Λ 0 b → J/ψ K − p and B 0 → J/ψK − π + decays, which cannot be present in the like-sign combinations. We use simulated samples of these decays to model their contributions. The Λ 0 b normalizations are determined from a previous analysis [25]. The background level in the opposite-sign combination (B 0 → J/ψ π + π − ) is studied by fitting the m(J/ψ π + π − ) distributions in bins of m(π + π − ). The resulting background distribution in the ±20 MeV When this data-simulation hybrid sample is used to extract the background parameters, a further re-weighting procedure is applied based on comparison of m(π + π − ) distributions between the overall fit and the background data points in Fig. 5(b).  Simulation LHCb Figure 7: The variation of ε 1 is shown as a function of m(π + π − ) and cos θ π + π − .
To better model the angular distributions in the ρ(770) mass region, the background is separated into the K * 0 reflection from B 0 , the ρ, and other backgrounds. The total background PDF is sum of these three components: where the N 's are normalizations, the contributing fractions having values of f K * 0 = (5.3 ± 0.2)% and f ρ = (9.5 ± 0.6)%; the other background is normalized as 1 The K * 0 background is modeled by the function where m 0 , Γ 0 , a, b, α 0 are free parameters determined by fitting to the B 0 → J/ψ K * 0 simulation. The last part (1 + p b1 cos χ + p b2 cos 2χ) is a function of the χ angle. We have verified that the three backgrounds have consistent χ distributions, thus the parameters p b1 and p b2 are determined by fitting all backgrounds simultaneously.

) [GeV]
π + π m( The ρ background is described by the function where m ρ , Γ ρ are free parameters. The parameters are obtained by fitting to simulated B 0 s → J/ψ η (→ ργ) events. The model for the remaining backgrounds is with the function Here the variable ζ = 2(m 2 hh − m 2 min )/(m 2 max − m 2 min ) − 1, where m min and m max are the fit boundaries, B 2 (ζ) is a fifth order Chebychev polynomial with coefficients b i (i = 1-5), and q(ζ) and p(ζ) are two first order Chebychev polynomials with parameters c j (j = 1-4). Figure 9 shows the projections of cos θ ππ and m(π + π − ) from the like-sign data combinations added with all the additional simulated backgrounds. The other background includes the Λ 0 b background and the combinatorial background which is described by the like-sign combinations. The fitted background parameters are given in Table 3. The cos θ J/ψ background distribution is shown in Fig. 10. Lastly, the χ background distribution, shown in Fig. 11 fit with the function 1 + p b1 cos χ + p b2 cos 2χ, determines the parameters p b1 = −0.004 ± 0.013 and p b2 = 0.070 ± 0.013.

Resonance models
To study the resonant structures of the decay B 0 → J/ψπ + π − we use 29 047 event candidates with invariant mass within ±20 MeV of the B 0 mass peak which include 10 207±178 background candidates. The background yield is fixed in the fit. Apart from non-resonant (NR) decays, the possible resonance candidates in the decay B 0 → J/ψπ + π − are listed in Table 4. We use Breit-Wigner (BW) functions for most of the resonances except f 0 (980). The masses and widths of the BW resonances are listed in Table 4. When used in the fit, they are fixed to these values except for the parameters of f 0 (500) which are allowed to vary by their uncertainties. For the f 0 (980) we use a Flatté shape [28]. Besides the mass, this shape has two additional parameters g ππ and g KK , which are fixed in the fit to the ones obtained from an amplitude analysis of B 0 s → J/ψπ + π − [27], where a large signal is evident. These parameters are m 0 = 945.4 ± 2.2 MeV, g ππ = 167 ± 7 MeV and g KK /g ππ = 3.47 ± 0.12. All background and efficiency parameters are fixed in the fit.
To determine the complex amplitudes in a specific model, the data are fitted maxi-  mizing the unbinned likelihood given as where N is the total number of candidates, and F is the total PDF defined in Eq. (15).  6 Fit results

Final state composition
In order to compare the different models quantitatively, an estimate of the goodness of fit is calculated from 4D partitions of the fitting variables. To distinguish between models, we use the Poisson likelihood χ 2 [29] 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. The χ 2 /ndf and the negative of the logarithm of the likelihood, −lnL, of the fits are given in Table 5 for various fitting models, where ndf, the number of degrees of freedom, is equal to N bin minus the number of fit parameters minus one. Here the five-resonance model (5R) contains the resonances: ρ(770), f 0 (500), f 2 (1270), ρ(1450) and ω(782), the "Best Model" adds a ρ(1700) resonance to the 5R model, the 7R model adds a f 0 (980) resonance to the Best Model, and the 7R+NR model adds a non-resonant component. We also give the change of lnL for various fits with respect to the 5R model in Table 5. The 7R model gives a slightly better likelihood compared to the Best Model, however, the decrease of the −lnL due to adding f 0 (980) is less than the expected ∆lnL at 3σ significance. Thus, we use the Best Model, which maintains a significance larger than 3σ for each resonance component, as our baseline fit, while the 7R model is only used to establish an upper limit on the presence of the f 0 (980). The Dalitz fit projections on the four observables: m(π + π − ), cos(θ π + π − ), cos θ J/ψ and χ are shown in Fig. 12 for the Best Model. Table 6 shows the summary of fit fractions of different components for various models. The fit fractions of the interference terms in the Best Model are computed using Eq. (13) and listed in Table 7. Table 8 shows the resonant phases from the Best Model. In the Best Model the CP -even components sum to (56.0±1.4)%, including the interference terms, so  that the CP -odd fraction is (44.0±1.4)%. The structure near the peak of the ρ(770) is due to ρ − ω interference. The fit fraction ratio is found to be where the uncertainties are statistical and systematic, respectively; wherever two uncertainties are quoted in this paper, they will be of this form. The systematic uncertainties will be discussed in detail in Sec. 6.3. The 7R model fit gives the ratio of observed decays into π + π − for f 0 (980)/f 0 (500) equal to (0.6 +0.7+3.3 −0.4−2.6 ) × 10 −2 . To determine the statistical uncertainty, the full error matrix and parameter values from the fit are used to generate 500 data-size sample parameter sets. For each set, the fit fractions are calculated. The distributions of the obtained fit fractions are described by bifurcated Gaussian functions. The widths of the Gaussians are taken as the statistical errors on the corresponding parameters. We will discuss the implications of this measurement in Sec. 7.
In Fig. 13 we show the fit fractions of the different resonant components in the Best Model. Table 9 lists the fit fractions and the transversity fractions of each contributing resonance. For a P -or D-wave resonance, we report its total fit fraction by summing all the three components. Table 10 shows the branching fractions of the resonant modes calculated by multiplying the fit fraction listed in Table 9 with B(B 0 → J/ψ π + π − ) = (3.97±0.09±0.11±0.16)×10 −5 , obtained from our previous study [3], where the uncertainties are statistical, systematic, and due to normalization, respectively. These branching fractions are proportional to the squares of the individual resonant amplitudes.

Angular moments
Angular moments are defined as an average of the spherical harmonics, Y 0 l (cos θ ππ ) , in each efficiency-corrected and background-subtracted π + π − invariant mass interval. The moment distributions provide an additional way of visualizing the effects of different resonances and their interferences, similar to a partial wave analysis. Figure 14 shows the distributions of the angular moments for the Best Model. In general the interpretation of these moments is that Y 0 0 is the efficiency corrected and background subtracted event distribution, Y 0 1 the sum of the interference between S-wave and P-wave and between Pwave and D-wave amplitudes,      3.4 ± 24.5 S-wave and D-wave amplitudes, Y 0 3 the interference between P-wave and D-wave, Y 0 4 the D-wave, and Y 0 5 results from an F-wave [19,30]. For the moments with odd-l, one will always find that B 0 and B 0 have opposite sign; thus the sum of their contributions is expected to be small.

Systematic uncertainties
The sources of systematic uncertainties on the results of the amplitude analysis are summarized in Table 11. Uncertainties due to particle identification and tracking are taken from Ref. [3] and are taken into account in the branching fraction results, but do not appear in the fit fractions as they are independent of pion kinematics. For the uncertainties due to the acceptance or background modeling, we repeat the data fit 100 times where the parameters of acceptance or background modeling are varied according to the corresponding error matrix. For the acceptance function, the error matrix is obtained by fitting the simulated acceptance as described in Sec. 5.1. For the background function, the error matrix is obtained by fitting the hybrid data-simulated sample as described in Sec. 5.2.
There is uncertainty on the fractions of sources in the hybrid MC-data sample for background modeling. Instead of using the fits to the π + π − mass distribution to determine the background fractions, we use the fractions found from the J/ψ π + π − mass fit shown in Fig. 3 that finds the Λ 0 b reflection is 9.6%, the B 0 reflection is 4.2%, the B 0 s background is 11.5% and the combinatorial part is 74.7%, instead of the ones found in Sec. 5.2. We then fit the new hybrid sample to get the background parameters. The data fit is repeated with the new background parameters; the changes on the fit results are added in quadrature with the uncertainties of the background modeling discussed above. The two background uncertainties have similar sizes.
We neglect the mass resolution in the fit where the typical resolution is 3 MeV. A previous study shows that the resolution effects are negligible except for the ω(782) resonance whose total fit fraction is underestimated by (0.09 ± 0.08)%. We take the quadrature of 0.09% and 0.08%, equal to 0.12%, as the systematic uncertainty of the total fit fraction of the ω(782). These uncertainties are included in the "Acceptance" category.
The uncertainties due to the fit model include adding each resonance that is listed in Table 4 but not used in the 7R model, varying the centrifugal barrier factors defined in Eq. (5) substantially, replacing f 0 (500) model by a Bugg function [31] and using the alternative Gounaris and Sakurai Model [32] for the various ρ mesons. The largest variation among those changes is assigned as the systematic uncertainty for modeling. We also find that increasing the default angular momentum L B for the P and D-wave cases gives negligible differences.
Finally, we repeat the amplitude fit by varying the mass and width of all the resonances except for the f 0 (980), in the 7R model within their errors one at a time, and add the changes in quadrature. For the f 0 (980) resonance, we change the resonance parameters m 0 , g ππ and g KK /g ππ to the values obtained from Solution II in [27] instead of using the ones obtained from Solution I. 7 Substructure of the f 0 (980) and f 0 (500) mesons The substructure of mesons belonging to the scalar nonet is controversial. Most mesons are thought to be formed from a combination of a q and a q. Some authors introduce the concept of qqqq states or superpositions of the tetraquark state with the qq state [33]. In either case, the I = 0 f 0 (500) and the f 0 (980) are thought to be mixtures of the underlying states whose mixing angle has been estimated previously. In the qq model, the mixing is parameterized by a normal 2×2 rotation matrix characterized by the angle ϕ m , so that the observed states are given in terms of the base states as |f 0 (980) = cos ϕ m |ss + sin ϕ m |nn |f 0 (500) = − sin ϕ m |ss + cos ϕ m |nn , In this case only the |dd part of the |nn wave function contributes (see Fig. 1). Thus we have where the Φ's are phase space factors [2,33,34]. 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. The 7R model fit gives the ratio of branching fractions B B 0 → J/ψ f 0 (980), f 0 (980) → π + π − B B 0 → J/ψ f 0 (500), f 0 (500) → π + π − = (0.6 +0.7+3.3 −0.4−2.6 ) × 10 −2 .
In order to set an upper limit on |ϕ m |, we simulate the final ϕ m measurement using as input the central value of the measured ratio, the full statistical error matrix obtained from the 7R model fit, and asymmetric Gaussian random variables different for the positive, +3.3%, and negative, −2.6%, systematic uncertainties (see Table 11). The resulting rate ratios of f 0 (980) to f 0 (500) are then multiplied by a factor of B (f 0 (500) → π + π − ) /B (f 0 (980) → π + π − ) × Φ(500) Φ(980) where a Gaussian random variable is used for B (f 0 (980) → π + π − ) to take into account the uncertainty in the measurement shown in Eq. (34). The upper limit at 90% confidence level is determined when 10% of the simulations exceed the limit value. We find tan 2 ϕ m ≡ r f σ = 1.1 +1.2+6.0 −0.7−0.7 × 10 −2 < 0.098 at 90% C.L which translates into a limit of |ϕ m | < 17 • at 90% CL , where we neglect the effect caused by the small systematic uncertainty on the ratio of phase space factors.
If the scalar meson substructure is tetraquark, the wave functions are: The ratio r f σ was predicted to be 1/2 for pure tetraquark states in Ref. [2]. The measured upper limit on r f σ of 0.098 at 90% CL deviates from the tetraquark prediction by 8 standard deviations.

Conclusions
We have studied the resonance structure of B 0 → J/ψπ + π − decays using a modified amplitude analysis. The decay distributions are formed by a series of final states described by individual π + π − interfering decay amplitudes. The data are best described by adding coherently the ρ(770), f 2 (1270), f 0 (500), ω(782), ρ(1450) and ρ(1700) resonances, with the largest component being the ρ(770). The final state is 56.0% CP -even, where we have taken into account both the fit fractions and the interference terms of the different components. Our understanding of the final state composition allows future measurements of CP violation in these resonant final states. These results supersede those obtained in Ref. [3].
There is no evidence for f 0 (980) resonance production. We limit the absolute value of the mixing angle between the lightest two scalar states, the f 0 (500) and the f 0 (980), in the qq model to be less than an absolute value of 17 • at 90% confidence level. We find that f 0 (980) production is much smaller than predicted for tetraquarks, which we rule out at the 8 standard deviation level using the model of Ref. [2]. Concern has been expressed [33] that if the f 0 (980) were a tetraquark state the measurement of the mixingdependent CP -violating phase in the decay B 0 s → J/ψ f 0 (980) could be affected due to additional decay mechanisms. Our result here alleviates this potential source of error. centres are supported by IN2P3 (France), KIT and BMBF (Germany), INFN (Italy), NWO and SURF (The Netherlands), PIC (Spain), GridPP (United Kingdom). We are indebted to the communities behind the multiple open source software packages on which we depend. We are also thankful for the computing resources and the access to software R&D tools provided by Yandex LLC (Russia).