Dalitz plot analysis of $B^0 \to \overline{D}^0 \pi^+\pi^-$ decays

The resonant substructures of $B^0 \to \overline{D}^0 \pi^+\pi^-$ decays are studied with the Dalitz plot technique. In this study a data sample corresponding to an integrated luminosity of 3.0 fb$^{-1}$ of $pp$ collisions collected by the LHCb detector is used. The branching fraction of the $B^0 \to \overline{D}^0 \pi^+\pi^-$ decay in the region $m(\overline{D}^0\pi^{\pm})>2.1$ GeV$/c^2$ is measured to be $(8.46 \pm 0.14 \pm 0.29 \pm 0.40) \times 10^{-4}$, where the first uncertainty is statistical, the second is systematic and the last arises from the normalisation channel $B^0 \to D^*(2010)^-\pi^+$. The $\pi^+\pi^-$ S-wave components are modelled with the Isobar and K-matrix formalisms. Results of the Dalitz plot analyses using both models are presented. A resonant structure at $m(\overline{D}^0\pi^-) \approx 2.8$ GeV$/c^{2}$ is confirmed and its spin-parity is determined for the first time as $J^P = 3^-$. The branching fraction, mass and width of this structure are determined together with those of the $D^*_0(2400)^-$ and $D^*_2(2460)^-$ resonances. The branching fractions of other $B^0 \to \overline{D}^0 h^0$ decay components with $h^0 \to \pi^+\pi^-$ are also reported. Many of these branching fraction measurements are the most precise to date. The first observation of the decays $B^0 \to \overline{D}^0 f_0(500)$, $B^0 \to \overline{D}^0 f_0(980)$, $B^0 \to \overline{D}^0 \rho(1450)$, $B^0 \to D_3^*(2760)^- \pi^+$ and the first evidence of $B^0 \to \overline{D}^0 f_0(2020)$ are presented.


Introduction
The study of the Cabibbo-Kobayashi-Maskawa (CKM) mechanism [1,2] is a central topic in flavour physics. Accurate measurements of the various CKM matrix parameters through different processes provide sensitivity to new physics effects, by testing the global consistency of the Standard Model. Among them, the CKM angle β is expressed in terms of the CKM matrix elements as arg(−V cd V * cb /V td V * tb ). The most precise measurements have been obtained with the B 0 → (cc)K ( * )0 decays by BaBar [3], Belle [4] and more recently by LHCb [5]. The decay 1 B 0 → D 0 π + π − through the b → cūd transition has sensitivity to the CKM angle β [6][7][8][9][10] and to new physics effects [11][12][13][14].
The Dalitz plot analysis [15] of B 0 → D 0 π + π − decays, with the D 0 → K + π − mode, is presented as the first step towards an alternative method to measure the CKM angle β. Two sets of results are given, where the π + π − S-wave components are modelled with the Isobar [16][17][18] and K-matrix [19] formalisms. Dalitz plot analyses of the decay B 0 → D 0 π + π − have already been performed by Belle [20,21] and BaBar [22]. Similar studies for the charged B decays B − → D ( * )+ π − π − have been published by the Bfactories [23,24]. The LHCb dataset offers a larger and almost pure signal sample. Feynman diagrams of the dominant tree level amplitudes contributing to the decay B 0 → D 0 π + π − are shown in Fig. 1. In addition to the interest for the CKM parameter measurements, the analysis of the Dalitz plot of the B 0 → D 0 π + π − decay is motivated by its rich resonant structure. The decay B 0 → D 0 π + π − contains information about excited D mesons decaying to Dπ, with natural spin and parity J P = 0 + , 1 − , 2 + , ... A complementary Dalitz plot analysis of the decay B 0 s → D 0 K − π + was recently published by LHCb [25,26], and constrains the 1 The inclusion of charge conjugate states is implied throughout the paper.
phenomenology of the D 0 K − (D − sJ ) and K − π + states. The spectrum of excited D mesons is predicted by theory [27,28] and contains the known states D * (2010), D * 0 (2400), D * 2 (2460), as well as other unknown states not yet fully explored. An extensive discussion on theory predictions for the cū, cd and cs mass spectra is provided in Refs. [26,29]. More recent measurements performed in inclusive decays by BaBar [30] and LHCb [29], have led to the observation of several new states: D * (2650), D * (2760), and D * (3000). However, their spin and parity are difficult to determine from inclusive studies. Orbitally excited D mesons have also been studied in semi-leptonic B decays (see a review in Ref. [31]) with limited precision. These are of prime interest both in the extraction of the CKM parameter |V cb |, where longstanding differences remain between exclusive and inclusive methods (see review in Ref. [32]), and in recent studies of B → D ( * ) τν τ [33] which have generated much theoretical discussion (see, e.g., Refs. [34,35]).
A measurement of the branching fraction of the decay B 0 → D 0 ρ 0 is also presented. This study helps in understanding the effects of colour-suppression in B decays, which is due to the requirement that the colour quantum numbers of the quarks produced from the virtual W boson must match those of the spectator quark to form a ρ 0 meson [36][37][38][39][40]. Moreover, using isospin symmetry to relate the decay amplitudes of B 0 → D 0 ρ 0 , B 0 → D − ρ + and B + → D 0 ρ + , effects of final state interactions (FSI) can be studied in those decays (see a review in Refs. [37,41]). The previous measurement for the branching fraction of B 0 → D 0 ρ 0 has limited precision, (3.2 ± 0.5) × 10 −4 [21], and is in agreement with theoretical predictions that range from 1.7 to 3.4 × 10 −4 [38,42].
The analysis of the decay B 0 → D 0 π + π − presented in this paper is based on a data sample corresponding to an integrated luminosity of 3.0 fb −1 of pp collision data collected with the LHCb detector. Approximately one third of the data was obtained during 2011 when the collision centre-of-mass energy was √ s = 7 TeV and the rest during 2012 with √ s = 8 TeV.
The paper is organised as follows. A brief description of the LHCb detector as well as the reconstruction and simulation software is given in Sec. 2. The selection of signal candidates and the fit to the B 0 candidate invariant mass distribution used to separate and to measure signal and background yields are described in Sec. 3. An overview of the Dalitz plot analysis formalism is given in Sec. 4. Details and results of the amplitude analysis fits are presented in Sec. 5. In Sec. 6 the measurement of the B 0 → D 0 π + π − branching fraction is documented. The evaluation of systematic uncertainties is described in Sec. 7. The results are given in Sec. 8, and a summary concludes the paper in Sec. 9.

The LHCb detector
The LHCb detector [59] 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 [60], 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 [61] placed downstream of the magnet. The tracking system provides a measurement of momentum, p, with a relative uncertainty that varies from 0.4% at low momentum to 0.6% at 100 GeV/c. The minimum distance of a track to a primary vertex, the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of p transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [62]. 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 [63].
The online event selection is performed by a trigger which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. At the hardware trigger stage, events are required to have a muon with high p T or a hadron, photon or electron with high transverse energy in the calorimeters. For hadrons, the transverse energy threshold is 3.5 GeV. The software trigger requires a two-, three-or four-track secondary vertex with a significant displacement from the primary pp interaction vertices (PVs). At least one charged particle must have a transverse momentum p T > 1.7 GeV/c and be inconsistent with originating from a PV. A multivariate algorithm [64] is used for the identification of secondary vertices consistent with the decay of a b hadron. The p T of the photon from D * − s decay is too low to contribute to the trigger decision.
Simulated events are used to characterise the detector response to signal and certain types of background events. In the simulation, pp collisions are generated using Pythia [65] with a specific LHCb configuration [66]. Decays of hadronic particles are described by EvtGen [67], in which final-state radiation is generated using Photos [68]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [69] as described in Ref. [70].

Event selection
Signal B 0 candidates are formed by combining D 0 candidates, reconstructed in the decay channel K + π − , with two additional pion candidates of opposite charge. Reconstructed tracks are required to be of good quality and to be inconsistent with originating from a PV. They are also required to have sufficiently high p and p T and to be within kinematic regions where reasonable particle identification (PID) performance is achieved, as determined by calibration samples of D * + → D 0 (K − π + )π + decays. The four final state tracks are required to be positively identified by the PID system. The D 0 daughters are required to form a good quality vertex and to have an invariant mass within 100 MeV/c 2 of the known D 0 mass [32]. The D 0 candidates and the two charged pion candidates are required to form a good vertex. The reconstructed D 0 and B 0 vertices are required to be significantly displaced from the PV. To improve the B 0 candidate invariant mass resolution, a kinematic fit [71] is used, constraining the D 0 candidate to its known mass [32].
By requiring the reconstructed D 0 vertex to be displaced downstream from the reconstructed B 0 vertex, backgrounds from both charmless B decays and direct prompt charm production coming from the PV are reduced to a negligible level. Background from D * (2010) − decays is removed by requiring m(D 0 π ± ) > 2.1 GeV/c 2 . Backgrounds from doubly mis-identified D 0 → K + π − or doubly Cabibbo-suppressed D 0 → K − π + decays are also removed by this requirement.
To further distinguish signal from combinatorial background, a multivariate analysis based on a Fisher discriminant [72] is applied. The sPlot technique [73] is used to statistically separate signal and background events with the B 0 candidate mass used as the discriminating variable. Weights obtained from this procedure are applied to the candidates to obtain signal and background distributions that are used to train the discriminant. The Fisher discriminant uses information about the event kinematic properties, vertex quality, IP and p T of the tracks and flight distance from the PV. It is optimised by maximising the purity of the signal events.
Signal candidates are retained for the Dalitz plot analysis if the invariant mass of the B 0 meson lies in the range [5250, 5310] MeV/c 2 and that of the D 0 meson in the range [1840, 1890] MeV/c 2 (called the signal region). Once all selection requirements are applied, less than 1 % of the events contain multiple candidates, and in those cases one candidate is chosen randomly.
Background contributions from decays with the same topology, but having one or two mis-identified particles, are estimated to be less than 1 % and are not considered in the Dalitz analysis. These background contributions include decays like Partially reconstructed decays of the type B 0 → D 0 π + π − X, where one or more particles are not reconstructed, have similar efficiencies to the signal channel decays. They are distributed in the region below the B 0 mass. By requiring the invariant mass of B 0 candidates to be larger than 5250 MeV/c 2 , these backgrounds are reduced to a negligible level, as determined by simulated samples of B 0 → D * 0 π + π − and B 0 → D * 0 ρ(770) with D * 0 decaying into D 0 γ or D 0 π 0 under different hypotheses for the D * 0 helicity. The signal and combinatorial background yields are determined using an unbinned extended maximum likelihood fit to the invariant mass distribution of B 0 candidates. The invariant mass distribution is shown in Fig. 2, with the fit result superimposed. The fit uses a Crystal Ball (CB) function [76] convoluted with a Gaussian function for the signal distribution and a linear function for the combinatorial background distribution in the mass range of [5250, 5500] MeV/c 2 . Simulated studies validate this choice of signal shape and the tail parameters of the CB function are fixed to those determined from simulation. Table 1 summarises the fit results on the free parameters, where µ B 0 is the mean peak position and σ G is the width of the Gaussian function. The parameter σ CB is the width of the Gaussian core of the CB function. The parameters f CB and p 1 give the fit fraction of the CB function and the slope of the linear function that describes the background distribution. The yields of signal (ν 0 s ) and background (ν 0 b ) events given in Table 1 are calculated within the signal region. The purity is (97.8 ± 0.2) %.

Dalitz plot analysis formalism
The analysis of the distribution of decays across the Dalitz plot [15] allows a determination of the amplitudes contributing to the three-body B 0 → D 0 π + π − decay. Two of the three Table 1: Results of the fit to the invariant mass distribution of B 0 → D 0 π + π − candidates. Uncertainties are statistical only.

Parameter
Value possible two-body invariant mass-squared combinations, which are connected by are sufficient to describe the kinematics of the system. The two observables m 2 (D 0 π − ) and m 2 (π + π − ), where resonances are expected to appear, are chosen in this paper. These observables are calculated with the masses of the B 0 and D 0 mesons constrained to their known values [32]. The invariant mass resolution has negligible effect and therefore it is not modeled in the Dalitz plot analysis. The total decay amplitude is described by a coherent sum of amplitudes from resonant or nonresonant intermediate processes as The complex coefficient c i and amplitude A i ( x) describe the relative contribution and dynamics of the i-th intermediate state, where x represents the (m 2 (D 0 π − ), m 2 (π + π − )) coordinates in the Dalitz plot. The Dalitz plot analysis determines the coefficients c i . In addition, fit fractions and interference fit fractions are also calculated to give a conventionindependent representation of the population of the Dalitz plot. The fit fractions are defined as and the interference fit fractions between the resonances i and j (i < j) are defined as where the integration is performed over the full Dalitz plot with m(D 0 π ± ) > 2.1 GeV/c 2 . Due to these interferences between different contributions, the sum of the fit fractions is not necessarily equal to unity.
The amplitude A i ( x) for a specific resonance r with spin L is written as The functions F (L) B (q, q 0 ) and F (L) r (p, p 0 ) are the Blatt-Weisskopf barrier factors [77] for the production, B 0 → rh 3 , and the decay, r → h 1 h 2 , of the resonance, respectively. The parameters p and q are the momenta of one of the resonance daughters (h 1 or h 2 ) and of the bachelor particle (h 3 ), respectively, both evaluated in the rest frame of the resonance. The value p 0 (q 0 ) represents the value of p (q) when the invariant mass of the resonance is equal to its pole mass. The spin-dependent F B and F r functions are defined as The value for the radius of the resonance, r BW , is taken to be 1.6 GeV −1 × c (= 0.3 fm) [78]. The function T L ( x) represents the angular distribution for the decay of a spin L resonance. It is defined as L = 0 : T 0 = 1, L = 1 : T 1 = 1 + y 2 cos θ × qp, L = 2 : T 2 = (y 2 + 3/2)( L = 3 : T 3 = 1 + y 2 (1 + 2y 2 /5)(cos 3 θ − 3 cos(θ)/5) × q 3 p 3 , L = 4 : T 4 = (8y 4 /35 + 40y 2 /35 + 1)(cos 4 θ − 30 cos 2 (θ)/35 + 3/35) × q 4 p 4 .
The helicity angle, θ, of the resonance is defined as the angle between the direction of the momenta p and q. The y dependence accounts for relativistic transformations between the B 0 and the resonance rest frames [79,80], where Finally, R( x) is the resonant lineshape and is described by the relativistic Breit-Wigner (RBW) function unless specified otherwise, where s = m 2 (h 1 h 2 ) and m r is the pole mass of the resonance; Γ (L) (s), the mass-dependent width, is defined as where Γ 0 is the partial width of the resonance, i.e., the width at the peak mass s = m r . The lineshapes of ρ(770), ρ(1450) and ρ(1700) are described by the Gounaris-Sakurai (GS) function [81], where The ρ − ω interference is taken into account by where Γ 0 is used, instead of the mass-dependent width Γ (L) (s), for ω(782) [82]. The D * (2010) − contribution is vetoed as described in Sec. 3. Possible remaining contributions from the D * (2010) − RBW tail or general D 0 π − P-waves are modelled as where β 1 and β 2 are free parameters. The π + π − S-wave contribution is modelled using two alternative approaches, the Isobar model [16][17][18] or the K-matrix model [19]. Contributions from the f 0 (500), f 0 (980), f 0 (2020) resonances and a nonresonant component are parametrised separately in the Isobar model and globally by one amplitude in the K-matrix model. In the Isobar model, the f 0 (2020) resonance is modelled by a RBW function and the modelling of the f 0 (500), f 0 (980) resonances and the nonresonant contribution are described as follows. The Bugg resonant lineshape [83] is employed for the f 0 (500) contribution, where The parameters are fixed to m r = 0.953 GeV/c 2 , s A = 0.41 m 2 π , b 1 = 1.302 GeV 2 /c 4 , b 2 = 0.340, A = 2.426 GeV 2 /c 4 and g 4π = 0.011 GeV/c 2 [83]. The phase-space factors of the decay channels ππ, KK and ηη correspond to ρ 1,2,3 (s), respectively and are defined as s , 1, 2, and 3 = π, K and η.
The Flatté formula [84] is used to describe the f 0 (980) lineshape, where The parameters g 1(2) and m r [46] are m r = 939.9 ± 6.3 MeV/c 2 , g 1 = 199 ± 30 MeV and g 2 /g 1 = 3.0 ± 0.3. The nonresonant contribution is described by Its modulus equals unity, and a slowly varying phase over m 2 (π + π − ) accounts for rescattering effects of the π + π − final state and α is a free parameter of the model.
The K-matrix formalism [19] describes the production, rescattering and decay of the π + π − S-wave in a coherent way. The scattering matrix S, from an initial state to a final state, is where I is the identity matrix, ρ is a diagonal phase-space matrix and T is the transition matrix. The unitarity requirement SS † = I gives The K-matrix is a Lorentz-invariant Hermitian matrix, defined as K −1 = T −1 + iρ. The amplitude for a decay process, is computed by combining the K-matrix obtained from scattering experiment with a production vector to describe process-dependent contributions. The K-matrix is modelled as a five-pole structure, where the indexes i, j = 1, 2, 3, 4, 5 correspond to five decay channels: ππ, KK, ηη, ηη and multi-meson (mainly 4π states) respectively. The coupling constant of the bare state α to the decay channel i, g α i , is obtained from a global fit of scattering data and is listed in Table 2 , regulates the singularities near the π + π − threshold, the so-called "Adler zero" [85,86]. The Hermitian property of the K-matrix imposes the relation f scatt ij = f scatt ji , and since only π + π − decays are considered, if i = 1 and j = 1, f scatt ij is set to 0. The production vector is modelled with where f prod 1j and β α are free parameters. The singularities in the K-matrix and the production vector cancel when calculating the amplitude matrix element.

Dalitz plot fit
An unbinned extended maximum likelihood fit is performed to the Dalitz plot distribution. The likelihood function is defined by The background probability density function (PDF) is given by f b ( x; θ b ) and is described in Sec. 5.1. The signal PDF, f s ( x i ; θ s ), is described by where the decay amplitude, M ( x; θ s ), is described in Sec. 4 and the efficiency variation over the Dalitz plot, ε( x), is described in Sec. 5.2. The fit parameters, θ s and θ b , include complex coefficients and resonant parameters like masses and widths. The value N is the total number of reconstructed candidates in the signal region. The number of signal and background events, ν s and ν b , are floated and constrained by the yields, ν 0 s and ν 0 b , determined by the D 0 π + π − mass fit and shown with their statistical uncertainties in Table 1.

Background modelling
The only significant source of candidates in the signal region, other than B 0 → D 0 π + π − decays, is from combinatorial background. It is modelled using candidates in the upper m(D 0 π + π − ) sideband ([5350, 5450] MeV/c 2 ) with a looser requirement on the Fisher discriminant, and is shown in Fig. 3. The looser requirement gives a similar distribution in the Dalitz plane but with lower statistical fluctuations. The Dalitz plot distribution of the combinatorial background events lying in the upper-mass sideband is considered to provide a reliable description of that in the signal region, as no dependence on m(D 0 π + π − ) is found by studying the Dalitz distribution in a different upper-mass sideband region. The combinatorial background is modelled with an interpolated non-parametric PDF [87,88] using an adaptive kernel-estimation algorithm [89].

Efficiency modelling
The efficiency function ε( x) accounts for effects of reconstruction, triggering and selection of the B 0 → D 0 π + π − signal events, and varies across the Dalitz plane. Two simulated samples are generated to describe its variation with several data-driven corrections. One is uniformly distributed over the phase space of the Dalitz plot and the other is uniformly distributed over the square Dalitz plot, which models efficiencies more precisely at the kinematic boundaries. The square Dalitz plot is parametrised by two variables m and θ that each varies between 0 and 1 and are defined as where m(π + π − ) max = m B 0 − m D 0 , m(π + π − ) min = 2m π and θ(π + π − ) is the helicity angle of the π + π − system. The two samples are fitted simultaneously with common fit parameters. A 4th-order polynomial function is used to describe the efficiency variation over the Dalitz plot. As the efficiency in the simulation is approximately symmetric over m 2 (D 0 π + ) and m 2 (D 0 π − ), the polynomial function is defined as with  The efficiency is corrected using dedicated control samples with data-driven methods. The corrections applied to the simulated samples include known differences between simulation and data that originate from the trigger, PID and tracking.

Results of the Dalitz plot fit
The Dalitz plot distribution from data in the signal region is shown in Fig. 5. The analysis is performed using the Isobar model and the K-matrix model. The nominal fit model in each case is defined by considering many possible resonances and removing those that do not significantly contribute to the Dalitz plot analysis. The resulting resonant contributions are given in Table 3 while the projections of the fit results are shown in Fig. 6 (Fig. 7) for the Isobar (K-matrix) model.
The comparisons of the S-wave results for the Isobar model and the K-matrix model are shown in Fig. 8. The results from the two models agree reasonably well for the amplitudes and phases. In the π + π − mass-squared region of [1.5, 4.0] GeV 2 /c 4 , small structures are seen in the K-matrix model, indicating possible contributions from f 0 (1370) and f 0 (1500) states. These contributions are not significant in the Isobar model and are thus not included in the nominal fit: adding them results in marginal changes and shows similar qualitative behaviour to the K-matrix model as displayed on Fig. 8. The measured S-waves from both models qualitatively agree with predictions given in Ref. [91].

Resonance
Spin  Fig. 10. There is a significant contribution from the D * J (2760) − resonance observed in Ref. [29] and a spin-3 assignment gives the best description. A detailed discussion on the determination of the spin of D * J (2760) is provided in Sec. 8.2. The fit quality is evaluated by determining a χ 2 value by comparing the data and the fit model in N bins = 256 bins that are defined adaptively to ensure approximately equal  6 Measurement of the B 0 → D 0 π + π − branching fraction Measuring the branching fractions of the different resonant contributions requires knowledge of the B 0 → D 0 π + π − branching fraction. This branching fraction is normalised relative to the B 0 → D * (2010) − π + decay that has the same final state, so systematic uncertainties are reduced. Identical selections are applied to select The kinematic constraints remove backgrounds from doubly mis-identified D 0 → K + π − or doubly Cabibbo-suppressed D 0 → K − π + decays and no requirement is applied on m(D 0 π + ). The invariant mass distributions of m(D 0 π − ) and m(D 0 π + π − ) for the B 0 → D * (2010) − π + candidates are shown in Fig. 11 and are fitted simultaneously to determine the signal and background contributions. The D * (2010) − signal distribution is  modelled by three Gaussian functions to account for resolution effects while its background is modelled by a phase-space factor. The modelling of the signal and background shapes in the m( The efficiencies for selecting B 0 → D * (2010) − π + and B 0 → D 0 π + π − decays are obtained from simulated samples. To take into account the resonant distributions in the Dalitz plot, the B 0 → D 0 π + π − simulated sample is weighted using the model described in the previous sections.
where the first uncertainty is statistical and the second uncertainty comes from the branching fraction of the normalisation channel.

Common systematic uncertainties and checks
Two categories of systematic uncertainties are considered, each of which is quoted separately. They originate from the imperfect knowledge of the experimental conditions and from the assumptions made in the Dalitz plot fit model. The Dalitz model-dependent uncertainties also account for the precision on the external parameters. The various sources are assumed to be independent and summed in quadrature to give the total. Experimental systematic uncertainties arise from the efficiency and background modelling and from the veto on the D * (2010) − resonance. Those corresponding to the signal efficiency are due to imperfect estimations of PID, trigger, tracking reconstruction effects, and to the finite size of the simulated samples. Each of these effects is evaluated by the differences between the results using efficiencies computed from the simulation and from the data-driven methods. The systematic uncertainties corresponding to the modelling of the small residual background are estimated by using different sub-samples of backgrounds. The systematic uncertainty due to the veto on the D * (2010) − resonance is assigned by changing the selection requirement from m 2 (D 0 π ± ) > 2.10 GeV/c 2 to 2.05 GeV/c 2 .
The systematic uncertainties related to the Dalitz models considered (see Sec. 4) include effects from other possible resonant contributions that are not included in the nominal fit, from the modelling of resonant lineshapes and from imperfect knowledge of the parameters of the modelling, i.e., the masses and widths of the π + π − resonances considered, and the resonant radius.
The radius of the resonances (r BW ) is set to a unique value of 1.6 GeV −1 × c in the nominal fit. In the systematic studies, it is floated as a free parameter and its best fit value is 1.84 ± 0.05 GeV −1 × c (1.92 ± 0.31 GeV −1 × c) for the Isobar (K-matrix) model. The value 1.85 GeV −1 × c is chosen to estimate the systematic uncertainties due to the imperfect knowledge of this parameter.
The masses and widths of the π + π − resonances considered are treated as free parameters with Gaussian constraints according to the inputs listed in Table 3. The differences between the results from those fits and those of the nominal fits are assigned as systematic uncertainties.
For the Isobar model, additional systematic uncertainties due to the modelling of the f 0 (500) and f 0 (980) resonances are considered. The Bugg model [83] for the f 0 (500) resonance and the Flatté model [84] for the f 0 (980) resonance, used in the nominal fit, are replaced by more conventional RBW functions. The masses and widths, left as free parameters, give 553 ± 15 MeV/c 2 and 562 ± 39 MeV, for the f 0 (500) meson and 981 ± 13 MeV/c 2 and 191 ± 39 MeV, for the f 0 (980) meson. The resulting differences to the nominal fit are assigned as systematic uncertainties.
The kinematic variables are calculated with the masses of the D 0 and B 0 mesons constrained to their known values [32]. These kinematic constraints affect the extraction of the masses and widths of the D 0 π − resonances. The current world average value for the B 0 meson mass is 5279.58 ± 0.17 MeV/c 2 and for the D 0 meson is 1864.84 ± 0.07 MeV/c 2 [32]. A conservative and direct estimation of the systematic uncertainties on the masses and  Several cross-checks have been performed to study the stability of the results. The analysis was repeated for different Fisher discriminant selection criteria, different trigger requirements and different sub-samples, corresponding to the two data-taking periods and to the two half-parts of the D 0 π + π − invariant mass signal region, above and below the B 0 mass [32]. Results from those checks demonstrate good consistency with respect to the nominal fit results. No bias is seen, therefore no correction is applied, nor is any related uncertainty assigned.

Systematic uncertainties on the
The systematic uncertainties related to the measurement of the B 0 → D 0 π + π − branching fraction are listed in Table 4. The systematic uncertainties on the PID, trigger, reconstruction and statistics of the simulated samples are calculated in a similar way to those of the Dalitz plot analysis. Other systematic uncertainties are discussed below. The systematic uncertainty on the modelling of the D 0 π − and D * (2010) − π + invariant mass distributions is estimated by counting the number of signal events in the B 0 signal region assuming a flat background contribution. The D * (2010) − mass region is restricted to the range [2007,2013] MeV/c 2 for this estimate. The calculated branching fraction is nearly identical to that from the mass fit and thus has a negligible contribution to the systematic uncertainty. The signal purity of B 0 → D * (2010) − π + is more than 99%. To account for the effect of resonant structures on the signal efficiency, the data sample is divided using an adaptive binning scheme. The average efficiency is calculated in a model independent way as where N i is the number of events in bin i and ε i is the average efficiency in bin i calculated from the efficiency model. The difference between this model-independent method and the nominal is assigned as a systematic uncertainty.

Significance of resonances
The Isobar and K-matrix models employed to describe the Dalitz plot of the B 0 → D 0 π + π − decay include all of the resonances listed in Table 3. The statistical significances of wellestablished π + π − resonances are calculated directly with their masses and widths fixed to the world averages. They are computed as the relative change of the minimum of the negative logarithm of the likelihood (NLL) function with and without a given resonance. Besides the π + π − resonances listed in Table 3, the significances of the f 0 (1370), f 0 (1500) and f 2 (1525) are also given. The results, expressed as multiples of Gaussian standard deviations (σ), are summarised in Table 5. All of the other π + π − resonances not listed in this Table have large   . The current analysis does not preclude a charged spin-1 D * state at around the same mass, but it is not sensitive to it with the current data sample size.
Studies have also been performed to validate the spin-0 hypothesis of the D * 0 (2400) − resonance, as the spin of this state has never previously been confirmed in experiment [32]. When moving to other spin hypotheses, the minimum of the NLL increases by more than 250 units in all cases, which confirms the expectation of spin 0 unambiguously.

Results of the Dalitz plot analysis
The shape parameters of the π + π − resonances are fixed from previous measurements except for the nonresonant contribution in the Isobar model. The fitted value of the parameter α defined in Eq. (20) is −0.363 ± 0.027, which corresponds to a 10 σ statistical significance compared to the case where there is no varying phase. An expansion of the model by including a varying phase in the D 0 π − axis is also investigated but no significantly varying  Table 6. The present precision on the mass and width of the D * 0 (2400) − resonance is improved with respect to Refs. [29,32]. The result for the width of the D * 2 (2460) − meson is consistent with previous measurements, whereas the result for the mass is above the world average which is dominated by the measurement using inclusive production by LHCb [29]. In the previous LHCb inclusive analysis, the broad D * 0 (2400) − component was excluded from the fit model due to a high correlation with the background lineshape parameters, while here it is included. The present result supersedes the former measurement. The Dalitz plot analysis used in this paper ensures that the background under the D * 2 (2460) − peak and the effect on the efficiency are under control, resulting in much lower systematic uncertainties compared to the inclusive approach.
The moduli and the phases of the complex coefficients of the resonant contributions, defined in Eq. (2), are displayed in Tables 7 and 8. Compatible results are obtained using both the Isobar and K-matrix models. The results for the fit fractions are given in Table 9, while results for the interference fit fractions are given in Appendix C. Pseudo experiments are used to validate the fitting procedure and no biases are found in the determination of parameter values. Table 6: Measured masses (m in MeV/c 2 ) and widths (Γ in MeV) of the D * 0 (2400) − , D * 2 (2460) − and D * 3 (2760) − resonances, where the first uncertainty is statistical, the second and the third are experimental and model-dependent systematic uncertainties, respectively.

Branching fractions
The measured branching fraction of the B 0 → D 0 π + π − decay in the phase-space region m(D 0 π ± ) > 2.1 GeV/c 2 is taking into account the systematic uncertainties reported in Table 4. The first uncertainty is statistical, the second systematic, and the third the uncertainty from the branching fraction of the B 0 → D * (2010) − π + normalisation decay channel. The result agrees with the previous Belle measurement (8.4 ± 0.4 ± 0.8) × 10 −4 [21] and the BaBar measurement (8.81 ± 0.18 ± 0.76 ± 0.78 ± 0.11)×10 −4 [22], obtained in a slightly larger phase-space region. A multiplicative factor of 94.5% (96.2%) is required to scale the Belle (BaBar) results to the same phase-space region as in this analysis.  The branching fraction of each quasi-two-body decay, B 0 → r i h 3 , with r i → h 1 h 2 , is given by where the resonant states (h 1 h 2 ) = (D 0 π − ), (π + π − ). The fit fractions F i , defined in Eq. (3), are obtained from the Dalitz plot analysis and are listed in Table 9. The correction factors, ε corr i , account for the cut-off due to the D * (2010) − veto. They are obtained by generating pseudo experiment samples for each resonance over the Dalitz plot and applying the same requirement (m(D 0 π ± ) > 2.1 GeV/c 2 ). They are summarised in Table 10. The correction factors are the same for the Isobar model and the K-matrix model. The effects due to the uncertainties of the masses and widths of the resonances are included in the uncertainties given in the table.
Using the overall B 0 → D 0 π + π − decay branching fraction, the fit fractions (F i ) and the correction factors (ε corr i ), the branching fractions of quasi two-body decays are calculated in Table 11 When accounting for the branching fractions of the ω(782) and f 2 (1270) to π + π − , one obtains the following results for the Isobar model and For the K-matrix model, one obtains and In both models, the fifth uncertainty is due to knowledge of the π + π − decay rates [32].
The parameter r f is related to the mixing angle by the equation in the qq model and by in the [qq ][qq ] tetraquark model [57,58]. The form factors F (B 0 → f 0 (980)) and F (B 0 → f 0 (500)) are evaluated at the four-momentum transfer squared equal to the square of the D 0 mass. Finally, values of the mixing angles as a function of form factor ratio are obtained in Fig. 13 for the qq model and the [qq ][qq ] tetraquark model. Such angles have also been computed by LHCb for the decays B 0 (s) → J/ψπ + π − [47][48][49]. The expectation is that the ratio of form factors should be close to unity. However, LHCb has recently performed a search for the decay B 0 s → D 0 f 0 (980) [95]. The limit set on this decay is below the value expected in a simple model based on our measured value of B(B 0 → D 0 f 0 (500)) and assuming equal form factors. More complicated models may be needed in order to explain all results.
The above discussion is one possible interpretation of the results. Another possible mechanism [91,96] involves the generation of pseudo-scalar resonances through the interactions of π + π − mesons.

Isospin analysis of the B → Dρ system
The measured branching fraction of the B 0 → D 0 ρ(770) 0 decay, presented in Table 11, can be used to perform an isospin analysis of the B → Dρ system. Isospin symmetry relates the amplitudes of the decays B + → D 0 ρ(770) + , B 0 → D − ρ(770) + , and B 0 → D 0 ρ(770) 0 , which can be written as linear combinations of the isospin eigenstates A I with I = 1/2 and 3/2 [37,41] The strong phase difference between the amplitudes A 1/2 and A 3/2 is denoted by δ Dρ . Final-state interactions between the states D 0 ρ 0 and D − ρ + may lead to a value of δ Dρ different from zero and through constructive interference, to a larger value of B(B 0 → D 0 ρ 0 ) than the prediction obtained within the factorisation approximation. In the heavy-quark limit, the factorisation model predicts [97,98] δ Dρ = O(Λ QCD /m b ) and the amplitude ratio and With a frequentist statistical approach [99], R Dρ and cos δ Dρ are calculated for the Isobar and K-matrix models in Table 13. These results are not significantly different from the predictions of factorisation models. As opposed to the theoretical expectations [37,41] and in contrast to the B → D ( * ) π system [40], non-factorisable final-state interaction effects do not introduce a sizeable phase difference between the isospin amplitudes in the B → Dρ system . The precision on R Dρ and cos δ Dρ is dominated by that of the branching fractions of the decays B + → D 0 ρ(770) + (14%) and B 0 → D − ρ(770) + (17%) [32]. The precision of the branching fraction of the B 0 → D 0 ρ(770) 0 decay is 7.3% (9.2%) for the Isobar (K-matrix) model (see Table 11).

Conclusion
A Dalitz plot analysis of the B 0 → D 0 π + π − decay is presented. The decay model contains four components from D 0 π − resonances, four P-wave π + π − resonances and one D-wave π + π − resonance. Two models are used to describe the S-wave π + π − resonances. The Isobar model uses four components, including the f 0 (500), f 0 (980), f 0 (2020) resonances and a nonresonant contribution. The K-matrix approach describes the π + π − S-wave using a 5 × 5 scattering matrix with a production vector. The overall branching fraction of B 0 → D 0 π + π − and quasi-two-body decays are measured. Significant contributions from the f 0 (500), f 0 (980), ρ(1450) and D * 3 (2760) − mesons are observed for the first time. For the latter, this is a confirmation of the observation from previous inclusive measurements, and the spin-parity of this resonance is determined for the first time to be J P = 3 − . This suggests a spectroscopic assignment of 3 D 3 , and shows that the 1D family of charm resonances can be explored in Dalitz plot analysis of B-meson decays in the same way as recently seen for the charm-strange resonances [25,26]. Evidence for the f 0 (2020) meson is also seen for the first time. The measured branching fractions of two-body decays are more precise than the existing world averages and there is good agreement between values from the Isobar and K-matrix models.
The masses and widths of the D 0 π − resonances are also determined. The measured masses and widths of the D * 0 (2400) − and D * 3 (2760) − states are consistent with the previous measurements. The precision on the D * 0 (2400) − meson is much improved. For the measurement on the mass and width of the D * 2 (2460) − meson, the broad D * 0 (2400) − component was excluded from the fit model in the former LHCb inclusive analysis [29], due to a high correlation with the background lineshape parameters, while here it is included. The present result therefore supersedes the former measurement.
The significant contributions found for both the f 0 (500) and f 0 (980) allow us to constrain on the mixing angle between the f 0 (500) and f 0 (980) resonances. An isospin analysis in the B → Dρ decays using our improved measurement of the branching fraction of the decay B 0 → D 0 ρ 0 is performed, indicating that non-factorisable effects from final-state interactions are limited in the Dρ system.      Figure 17: The first eight unnormalised Legendre polynomial weighted moments (0 to 7 correspond to (a) to (h)) for background-subtracted and efficiency-corrected B 0 → D 0 π + π − data and the Dalitz plot fit results as a function of m 2 (π + π − ). Only results in the region m 2 (π + π − ) < 3 GeV 2 /c 4 are shown.
B Systematic uncertainties on the parameters in the Dalitz plot analysis B.1 Systematic uncertainties for the Isobar model Table 14: Systematic uncertainties on the D 0 π − resonant masses (MeV/c 2 ) and widths (MeV) for the Isobar model.     Table 18: Systematic uncertainties on the D 0 π − resonant masses (MeV/c 2 ) and widths (MeV) for the K-matrix model.

C Results for the interference fit fractions
The central values of the interference fit fractions for the Isobar (K-matrix) model are given in Table 22 (Table 23). The statistical, experimental systematic and model-dependent uncertainties on these quantities are given in Tables 24, 25 and 26 (Tables 27, 28 and 29).   Table 9.    Table 9.    Table 9.    Table 9.  Table 9.  Table 9.  Table 9.  Table 9.

D Results of the K-matrix parameters
The moduli and phases of the K-matrix parameters in Eq. (25) are listed in Table 30. The break-down of systematic uncertainties are shown in Tables 31 and 32.