Amplitude analysis and branching fraction measurement of the decay $D_{s}^{+} \to \pi^{+}\pi^{0}\pi^{0}$

Using a data set corresponding to an integrated luminosity of 6.32~$\rm fb^{-1}$ recorded by the BESIII detector at center-of-mass energies between 4.178 and 4.226~GeV, an amplitude analysis of the decay $D_{s}^{+} \to \pi^{+}\pi^{0}\pi^{0}$ is performed, and the relative fractions and phases of different intermediate processes are determined. The absolute branching fraction of the decay $D_{s}^{+} \to \pi^{+}\pi^{0}\pi^{0}$ is measured to be $(0.50\pm 0.04_{\text{stat}}\pm 0.02_{\text{syst}})\%$. The absolute branching fraction of the intermediate process $D_{s}^{+} \to f_0(980)\pi^{+}, f_0(980)\to\pi^{0}\pi^{0}$ is determined to be $(0.21\pm 0.03_{\text{stat}}\pm 0.03_{\text{syst}})\%$.


Introduction
The constituent quark model has been very successful in explaining the composition of hadrons in the past few decades. In this model, the observed meson spectrum is described as bound qq states grouped into SU(n) flavor multiplets. The nonets of pseudo-scalar, vector and tensor mesons have been well identified. Nevertheless, the identification of the scalar-meson nonet is still ambiguous. Distinguishing scalar mesons from non-resonant background is rather difficult due to their broad widths and non-distinctive angular distribution. There are copious candidates for the J P C = 0 ++ nonets [1]. The case with isospin zero states, e.g. f 0 (500), f 0 (980), f 0 (1370), f 0 (1500), and f 0 (1710), is the most complicated from both experimental and theoretical points of view. Among them, the f 0 (980) meson, as a possible tetraquark candidate [2][3][4], is particularly interesting and can be studied via the hadronic decays D + s → π + π 0 π 0 , D + s → π + π + π − and D + s → K + K − π + . Charge conjugation is implied throughout in this paper. The current published branching fraction (BF) of D + s → f 0(2) π + from the D + s → π + π + π − decays has large discrepancies [1,5,6] with that measured from the D + s → K + K − π + decays. The f 0(2) contributions may suffer from the contaminations of a 0 (980) → K + K − or ρ → π + π − and the D + s → π + π 0 π 0 decays offer a cleaner environment due to absence of these contributions.
Furthermore, hadronic D + s decays can be used to probe the interplay of short-distance weak-decay matrix elements and long-distance QCD interactions, and the measured BFs provide valuable information concerning the amplitudes and phases that induce in decay processes [7][8][9][10][11].

Detector and data sets
The BESIII detector [13] records symmetric e + e − collisions provided by the BEPCII storage ring in the range from √ s = 2.00 GeV to √ s = 4.95 GeV [14,15]. The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a magnetic field of 1.0 T. The solenoid is supported by an octagonal fluxreturn yoke with resistive plate counter muon identification modules interleaved with steel. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region. The time resolution in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps. The end cap TOF system was upgraded in 2015 using multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [16][17][18].
The data samples used in this analysis are listed in Table 1 [19,20]. Since the cross section of D * ± s D ∓ s production in e + e − annihilation is about a factor of twenty larger than that of D + s D − s [21], and the D * + s meson decays to γD + s with a dominant BF of (93.5 ± 0.7)% [1], the signal events discussed in this paper are selected from the process e + e − → D * ± s D ∓ s → γD + s D − s . Simulated data samples produced with a geant4-based [22] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector and the detector response, are used to determine detection efficiencies and to estimate backgrounds. The simulation models the beam energy spread and initial state radiation (ISR) in the e + e − annihilations with the generator kkmc [23,24]. The inclusive MC sample includes the production of open charm processes, the ISR production of vector charmonium(-like) states, and the continuum processes incorporated in kkmc [23,24]. The known decay modes are modelled with evtgen [25,26] using BFs taken from the Particle Data Group [1], and the remaining unknown charmonium decays are modelled with lundcharm [27,28]. Final state radiation (FSR) from charged final state particles is incorporated using photos [29].

Event selection
The data samples were collected just above the D * ± s D ∓ s threshold, which allows to extract relatively pure samples for amplitude analysis and measurements of absolute BFs of the hadronic D + s meson decays with a tag method. The tag method has single-tag (ST) and double-tag (DT) candidates. The ST candidates are those D ± s mesons without further requirements on the remaining tracks and EMC showers. The DT candidates are identified by fully reconstructing the D + s D − s mesons, where one of the D ± s mesons decays into the signal mode D + s → π + π 0 π 0 and the other to a tag mode. The D ± s mesons are reconstructed through the final state particles, i.e. π ± , K ± , η, η , K 0 S and π 0 , whose selection criteria is discussed below.
For charged tracks not originating from K 0 S decays, the distance of closest approach to the interaction point is required to be less than 10 cm along the beam direction and less than 1 cm in the plane perpendicular to the beam. Particle identification (PID) for charged tracks combines measurements of the specific ionization energy losses in the MDC (dE/dx) and the flight time in the TOF to form a likelihood L(h) (h = K, π) for the hypothesis of being a hadron h. A charged hadron is identified as a kaon if L(K) is larger than L(π), otherwise it is identified as a pion.
The K 0 S candidates are reconstructed from two oppositely charged tracks satisfying |cos θ| < 0.93 and the distance of closest approach along the beam direction must be less than 20 cm. The two charged tracks coming form the K 0 S are assigned as π + π − without imposing further PID criteria. They are constrained to originate from a common vertex and are required to have an invariant mass within |M π + π − − m K 0 S mass taken from PDG [1]. Photon candidates are identified using showers in the EMC. The deposited energy of each shower must be more than 25 MeV in the barrel region (|cos θ| < 0.80) and more than 50 MeV in the end cap region (0.86 < |cos θ| < 0.92). The angle between the position of each shower in the EMC and any charged track must be greater than 10 degrees to exclude showers originating from charged tracks. The difference between the EMC time and the event start time is required to be within [0, 700] ns to suppress electronic noise and showers unrelated to the event.
The π 0 (η) candidates are reconstructed through π 0 → γγ (η → γγ) decays, with at least one photon in the barrel. The invariant masses of the photon pairs for π 0 and η candidates must be in the ranges [0.115, 0.150] GeV/c 2 and [0.490, 0.580] GeV/c 2 , respectively, which are about three times the resolution of the detector. A kinematic fit that constrains the γγ invariant mass to the π 0 or η nominal mass [1] is performed to improve the mass resolution. The χ 2 of the kinematic fit is required to be less than 30. The η candidates are formed from the π + π − η combinations with an invariant mass within a range of [0.946, 0.970] GeV/c 2 .
Seven tag modes are used to reconstruct the tag D − s candidate and its mass (M tag ) is required to fall within the mass window listed in Table 2. The recoiling mass of the tag D − is calculated in the e + e − center-of-mass system, where p Ds is the momentum of the D − s candidate in the e + e − center-of-mass frame and m Ds is the known D − s mass [1]. The value of M rec is required to be within the region listed in Table 1.

Further selections
The following selection criteria are further applied in order to obtain data samples with high purity for the amplitude analysis. The selection criteria discussed in this section are not used in the BF measurement. An eight-constraint kinematic fit is performed to select photon from D * ± s decays and the best DT candidates assuming D − s candidates decaying to one of the tag modes and D + s decaying to the signal mode with two hypotheses: the signal D + s comes from a D * + s or the tag D − s comes from a D * − s . In this kinematic fit, the total four-momentum is constrained to the initial four-momentum of the e + e − system, and the invariant masses of (γγ) π 0 , (π + π − ) K 0 S , tag D − s , and D * +(−) s candidates are constrained to the corresponding known masses [1]. The best combination is chosen with the minimum χ 2 . After the selection, an additional constraint of the signal π + π 0 π 0 invariant mass to the known D + s mass is added and the updated four-momenta of final-state particles from the kinematic fit are used for the amplitude analysis in order to ensure that all candidates fall within the phase-space boundary.
The energy of the transition photon from D * + s → γD + s is required to be smaller than 0.18 GeV. The recoiling mass against this photon and the signal D + s is required to fall in the range [1.952, 1.995] GeV/c 2 . The D + s → π + π 0 η decay contributes to the background when π 0 η is misreconstructed as π 0 π 0 . This background is reduced via an "η" veto to reject events which simultaneously satisfy where M γ 1 γ 3 and M γ 2 γ 4 are the invariant masses of any combinations of the photons used to reconstruct the two π 0 s in the signal decay. There is also background originating from decays by exchanging K − and π 0 2 (K + and π 0 1 ). This background is excluded by rejecting events which simultaneously satisfy Figure 1 shows the fits to the invariant-mass distributions of the D + s candidates reconstructed in the signal mode, M sig , for the two data samples. The signal is described by a MC-simulated shape convolved with a Gaussian resolution function and the background is described by a simulated shape based on inclusive MC samples. Finally, a mass window

Fit method
The intermediate-resonant composition is determined by an unbinned maximum-likelihood fit to data. The likelihood function L is constructed with a signal-background combined probability density function (PDF), which depends on the momenta of the three final state particles: where i and j indicate the data sample groups and the final-state particles, respectively, N D,i is the number of candidate events in the data i, f S (f B ) is the signal (background) PDF and w is the purity of signal.
The signal PDF is written as where (p j ) is the detection efficiency modeled by a RooNDKeysPdf derived from phase space MC sample, A(p j ) represents the total amplitude, and R 3 is the standard element of three-body phase space. The isobar formulism is used to model the total amplitude. The total amplitude is the coherent sum of individual amplitudes of intermediate processes, A = ρ n e iφn A n where magnitude ρ n and phase φ n are the free parameters to be determined by data. The amplitude of the n th intermediate process (A n ) is given by where S n is the spin factor (Sec The background PDF is given by is the efficiency-corrected background shape. The background events in the signal region from the generic MC sample are used to derive the background shape B(p j ) with RooNDKeysPdf [30]. RooNDKeysPdf is a kernel estimation method [31] implemented in RooFit [30] which models the distribution of an input dataset as a superposition of Gaussian kernels. The M π + π 0 and M π 0 π 0 distributions of events outside the M sig signal region between the data and the generic MC samples are compared to check validity of the background from the generic MC samples. The distributions of background events from the generic MC samples within and outside the M sig signal region are also examined. They are found to be compatible within statistical uncertainties. Note that the (p j ) term in Eq. (4.4) is explicitly written out as it is independent of the fitted variables and is dropped during the log-likelihood fit. The normalization integral terms in the signal and background PDF are handled by MC integration, where X(p j ) is |A(p j )| 2 or B (p j ), k is the index of the k th event, N G is the number of the generated MC events and N M is the number of the selected MC events. The D + s meson in the MC samples used here decays to π + π 0 π 0 according to the PDF M g (p j ), while the D − s meson decays into one of the tag modes. These MC samples are generated with different √ s according to the luminosities and cross sections, and satisfy all selection criteria as those of the data samples. At the beginning, a preliminary PDF is used, and then a recursive process is performed until the result converges. To account for any bias caused by differences in PID or tracking efficiencies between data and MC simulation, each signal MC event is weighted with a ratio, γ (p), of the efficiency of data to that of MC simulation and the MC integration then becomes (4.6)

Spin factors
The spin-projection operators are defined as [32] P  The quantities p a , p b , and p c are the momenta of particles a, b, and c, respectively, and r a = p b − p c . The covariant tensors are given bỹ The spin factors for S, P , and D wave decays are where theT (l) factors have the same definition ast (l) . The tensor describing the D + s decay is denoted byT and that of the a decay is denoted byt.

Blatt-Weisskopf barrier factors
For the process a → bc, the Blatt-Weisskopf barrier F L (p j ) [33] is parameterized as a function of the angular momentum L and the momentum q of the final-state particle b or c in the rest system of a,

Propagator
The intermediate resonances f 2 (1270) and f 0 (1370) are parameterized as relativistic Breit-Wigner functions, where s a denotes the invariant-mass squared of the two final-state particles considered; m 0 and Γ 0 are the mass and the width of the intermediate resonance, respectively, and are fixed to the PDG values [1]. The f 0 (980) resonance is represented by the Flatté formula [34], where s π 0 π 0 is the π 0 π 0 invariant-mass squared and g 1,2 are coupling constants to the corresponding final states. The parameters are fixed to g 1 = 0.165 GeV/c 2 , g 2 /g 1 = 4.21 and m f 0 (980) = 955 MeV/c 2 , as reported in Ref.

Fit results
The Dalitz plot of M 2 π + π 0 versus M 2 π + π 0 for the data samples is shown in Fig. 2(a) and that for the signal MC samples generated based on the results of the amplitude analysis is shown in Fig. 2(b). In the fit, the magnitude and phase of the reference amplitude D + s → f 0 (980)π + are fixed to 1.0 and 0.0, respectively, while those of other amplitudes are left floating. The masses and widths of all resonances are fixed to the corresponding PDG averages [1], and w i are fixed to the purities discussed in Sec. 4.1. The systematic uncertainties associated with these fixed parameters are considered by repeating the fit after variation of the fixed parameters according to their uncertainties.
Besides the dominant amplitudes D + s → f 0 (980)π + , D + s → f 0 (1370)π + , and D + s → f 2 (1270)π + , we have tested all possible intermediate resonances including ρ(1450), f 0 (1500), ρ(1700), (ππ) S , (ππ) P , (ππ) D etc., where the subscript denotes a relative S (P or D) wave between final-state particles. We have also examined all possible combinations of these intermediate resonances to check their significances, correlations, and interferences. By requiring a significance larger than 3σ, eventually, D + s → f 0 (980)π + , D + s → f 0 (1370)π + , D + s → f 2 (1270)π + , D + s → π + (π 0 π 0 ) D , and D + s → (π + π 0 ) D π 0 are chosen for the nominal set. Note that D + s → f 0 (500)π + is tested but it has a significance less than 2σ.  In the calculation of fit fractions (FFs) for individual amplitudes, the phase-space MC truth information is involved with neither detector acceptance nor resolution. The FF for the n th amplitude is defined as Ngen |A| 2 , (4.14) where N gen is the number of the phase-space MC events at generator level. Interference between the n th and the n th amplitudes (IN) is defined as (for n < n only) The statistical fluctuations of FFs are obtained by randomly sampling the fit variables according to their fitted values and covariance matrix. The distribution of each FF is fitted with a Gaussian function and the width of the Gaussian function is defined as the statistical uncertainty of the FF. The phases, FFs, and statistical significances for the amplitudes are listed in Table 3. The interferences between amplitudes are listed in Table 4. The Dalitz plot projections are shown in Fig. 3. The sum of the FFs is not unity due to interferences between amplitudes. Other tested amplitudes, but not included in the nominal fit, and their significances are listed in Table 5. Table 3. The phases, FFs, and statistical significances for the amplitudes. The first and second uncertainties in the phases and FFs are statistical and systematic, respectively. The total FF is 111.4%.

Systematic uncertainties for the amplitude analysis
The systematic uncertainties for the amplitude analysis are summarized in Table 6, with their definitions described below: i Resonant parameters. The masses and the widths of f 0 (1370) and f 2 (1270) are varied by their corresponding uncertainties [1]. The mass and coupling constants of the f 0 (980) Flatté formula are varied according to Ref. [34]. The changes of the phases φ and FFs are assigned as the associated systematic uncertainties. Amplitude  ii R values. The associated systematic uncertainties are estimated by repeating the fit procedure by varying the radii of the intermediate state and D + s mesons within 1 GeV −1 .
iii Background estimation. First, the purities of signals for the two sample groups, i.e. w in Eq. (4.1) are varied by their corresponding statistical uncertainties to study uncertainties associated with backgrounds. The differences caused by the variation are assigned as the uncertainties. Second, an alternative MC-simulated shape is used to examine the uncertainty arising from the background shape modeling. Alternative background shapes are extracted with the relative fractions of the dominant backgrounds from e + e − → qq and non-D * ± s D ∓ s open-charm processes varied by the statistical uncertainties of their cross sections.
iv Resonances with significances less than 3σ. The corresponding uncertainties are taken to be the differences of the phases φ and FFs with and without the intermediate resonances with statistical significances less than 3σ.  v Experimental effects. To estimate the systematic uncertainty related to the difference between MC simulation and data associated with the PID and tracking efficiencies, γ in Eq. (4.6), the amplitude fit is performed varying the PID and tracking efficiencies according to their uncertainties. The differences from the nominal results are so tiny that this source of systematic uncertainty is negligible.

Branching fraction measurement
In addition to the selection criteria for final-state particles described in Sec. 3, it is required that π + must have momentum greater than 100 MeV/c to remove soft π + from D * + decays. The best tag candidate with M rec closest to the D * + s known mass [1] is chosen if there are multiple ST candidates. The data sets are organized into three sample groups, 4.178 GeV, 4.189-4.219 GeV, and 4.226 GeV, that were acquired during the same year under consistent running conditions.
The yields for various tag modes are obtained by fitting the corresponding M tag distributions and listed in Table 7. As an example, the fits to the M tag spectra of the ST candidates in the data sample at √ s = 4.178 GeV are shown in Fig. 4. In the fits, the signal is modeled by a MC-simulated shape convolved with a Gaussian function to take into account the data-MC resolution difference. The background is described by a second-order Chebyshev function. MC studies show that there is no significant peaking background in any tag mode, except for D − → K 0 S π − and D − s → ηπ + π − π − faking the D − s → K 0 S K − and D − s → π − η tags, respectively. Therefore, the MC-simulated shapes of these two peaking background sources are added to the background models.  Once a tag mode is identified, the signal decay D + s → π + π 0 π 0 is searched for at the recoiling side. In the case of multiple candidates, the DT candidate with the average mass, (M sig + M tag )/2, closest to the D + s nominal mass is retained. A K 0 S → π 0 π 0 mass veto, M π 0 π 0 / ∈ (0.458, 0.520) GeV/c 2 , is applied on the signal D + s to remove the peaking background D + s → K 0 S π + . To measure the BF, we start from the following equations for each tag mode: is the ST yield for the tag mode; N DT tag,sig is the DT yield; B tag and B sig are the BFs of the tag and signal modes, respectively; ST tag is the ST efficiency to reconstruct the tag mode; and DT tag,sig is the DT efficiency to reconstruct both the tag and the signal decay modes. In the case of more than one tag modes and sample groups,

Tag mode
where α represents tag modes in the i th sample group. By isolating B sig , we find where N ST α,i and ST α,i are obtained from the data and inclusive MC samples, respectively.
DT α,sig,i is determined with signal MC samples with D + s → π + π 0 π 0 events are generated according to the results of the amplitude analysis. The BF for π 0 → γγ is introduced to account for the fact that the signal is reconstructed through this decay.
The DT yield N DT total is found to be 587 ± 44 from the fit to the M sig distribution of the selected D + s → π + π 0 π 0 candidates. The fit result is shown in Fig. 5, where the signal shape is described by a MC-simulated shape convolved with a Gaussian function to take into account the data-MC resolution difference. The background is described by a simulated shape from the inclusive MC sample. A small peaking background originating from D 0 → K − π + π 0 is considered in the inclusive MC sample. Taking the difference in π 0 reconstruction efficiencies for each signal mode between data and MC simulation into account by multiplying the efficiencies by a factor of 99.5% for each π 0 , we determine the BF of D + s → π + π 0 π 0 to be (0.50 ± 0.04 stat ± 0.02 syst )%. The relative systematic uncertainty for the total yield of the ST D − s mesons is assigned to be 0.4% by examining the changes of the fit yields when varying the signal shape, background shape, and taking into account the background fluctuation in the fit. The systematic uncertainty due to the signal shape is studied by repeating the fit without the convolved Gaussian. The MC-simulated background shape is altered by varying the relative fractions of the dominant backgrounds from e + e − → qq or non-D * + s D − s open-charm processes by their statistical uncertainties of their related cross sections. The largest change is taken as the corresponding systematic uncertainty. The π + tracking (PID) efficiency is studied with the processes e + e − → K + K − π + π − (e + e − → K + K − π + π − (π 0 ) and π + π − π + π − (π 0 )). The systematic uncertainty due to tracking (PID) efficiency is estimated to be 1%(1%). The systematic uncertainty of the π 0 reconstruction efficiency is investigated by using a control sample of the process e + e − → K + K − π + π − π 0 . The selection criteria listed in Sec. 3 are used to reconstruct the two kaons and the two pions. The recoiling mass distribution of K + K − π + π − is fitted to obtain the total number of π 0 's and the π 0 selection is applied to determine the number of reconstructed π 0 's. The average ratio between data and MC efficiencies of π 0 reconstruction, weighted by the corresponding momentum spectra, is estimated to be 0.995 ± 0.008. After correcting the simulated efficiencies to data by this ratio, the residual uncertainty 0.8% is assigned as the systematic uncertainty arising from each π 0 reconstruction. The uncertainty due to the limited MC statistics is obtained by where f i is the tag yield fraction, and i and δ i are the signal efficiency and the corresponding uncertainty of tag mode i, respectively. The uncertainty from the amplitude analysis model is estimated by varying the model parameters based on their error matrix. The distribution of 600 efficiencies resulting from this variation is fitted by a Gaussian function and the fitted width divided by the mean value is taken as a relative uncertainty. All of the systematic uncertainties are summarized in Table 8. Adding them in quadrature gives a total systematic uncertainty in the BF measurement of 4.0%.