First Measurement of Polarizations in the Decay D 0 → ωφ

Using a data sample corresponding to an integrated luminosity of 2 . 93 fb − 1 collected at a center-of-mass energy ﬃﬃﬃ s p ¼ 3 . 773 GeV by the BESIII detector, the decay D 0 → ωϕ is observed for the first time. The branching fraction is measured to be ð 6 . 48 (cid:2) 0 . 96 (cid:2) 0 . 40 Þ × 10 − 4 with a significance of 6 . 3 σ , where the first and second uncertainties are statistical and systematic, respectively. An angular analysis reveals that the ϕ and ω mesons from the D 0 → ωϕ decay are transversely polarized. The 95% confidence level upper limit on longitudinal polarization fraction is set to be less than 0.24, which is inconsistent with current theoretical expectations and challenges our understanding of the underlying dynamics in charm meson decays.

Using a data sample corresponding to an integrated luminosity of 2.93 fb −1 collected at a centerof-mass energy √ s = 3.773 GeV by the BESIII detector, the decay D 0 → ωφ is observed for the first time. The branching fraction is measured to be (6.48 ± 0.96 ± 0.40) × 10 −4 with a significance of 6.3σ, where the first and second uncertainties are statistical and systematic, respectively. An angular analysis reveals that the φ and ω mesons from the D 0 → ωφ decay are transversely polarized. The 95% confidence level upper limit on longitudinal polarization fraction is set to be less than 0.24, which is inconsistent with current theoretical expectations and challenges our understanding of the underlying dynamics in charm meson decays.
Understanding the long-distance contributions to D 0 − D 0 mixing is crucial in tests of the Standard Model (SM) [1]. These contributions arise in two-body hadronic decays of the D 0 meson, such as when the D 0 meson decays to two vector (V ) mesons in the process D 0 → V V , which is expected to account for 10% of the total D 0 decay width [1]. In constrast to scalar and pseudoscalar mesons, vector mesons can be produced in three polarization states. Therefore, the decay D 0 → V V produces a longitudinal partial-wave amplitude (H 0 ), which is CPeven, and two transverse partial-wave amplitudes (H ± ), which are superpositions of CP -even and CP -odd states. The polarization in the D 0 → V V decay is sensitive to the V -A structure of electroweak interactions in the SM, spin correlations, and final state interactions, among other effects [2][3][4]. Thus, in addition to the partial decay widths, the polarization is an interesting observable.
In the last two decades, a "polarization puzzle" has arisen in the decays of heavy mesons to two vectors. In the beauty sector, naive power counting predicts that B → V V decays are dominated by longitudinal polarization, since the transverse polarization amplitude suffers from a helicity flip suppression on the order of Λ QCD /m b . The B factories confirmed this prediction in the decay modes B 0 → ρ + ρ − , B + → ρ 0 ρ + , and B + → ρ 0 K * + , each of which were found to favor the longitudinal configuration [5][6][7]. However, in contrast, the longitudinal amplitude does not dominate in the decay B → φK * (892), where it provides only 50% of the rate [5,8]. The observed deviation in this case can be explained either by penguin annihilation [9,10], rescattering of the finalstate interactions [11,12], or new physics beyond the SM model [13][14][15][16][17]. In the charm sector, the situation is more complicated since the heavy quark expansion method is less reliable in this quark mass regime. The naive factorization model [18] and the Lorentz invariant-based symmetry model [19] predict the longitudinal polarization fraction, f L = H 2 0 /(H 2 0 + H 2 − + H 2 + ), to be ∼0.5 and 0.33 for D 0 → V V , respectively. However, a previous measurement of the decay D 0 →K * 0 ρ 0 actually shows transverse polarization dominates [20], though the measurement suffers from a large uncertainty. In addition, the most precise measurement of the decay D 0 → ρ 0 ρ 0 by the FOCUS collaboration shows large longitudinal polarization with f L = (71 ± 4 ± 2)% [21].
To date, all experimential measurements of the helicity in D 0 → V V decays have been contrary to theoretical predictions, and the puzzle needs to be confirmed with more precise measurements, and validated using more decay modes. The singly-Cabibbo-suppressed decay D 0 → ωφ, which is kinematically allowed, can occur via internal emission of a W + boson. The branching fraction (BF) of D 0 → ωφ is predicted to be at the level of 10 −3 -10 −4 by various phenomenological models [1,[22][23][24], and its polarization features are particularly attractive. No signal for D 0 → ωφ has yet been observed, and only an upper limit on the BF, B(D 0 → ωφ) < 2.1 × 10 −3 [25], is available. An observation of D 0 → ωφ is desired to shed light on the polarization puzzle, test different theoretical models [22,26,27], measure CP -violating parameters and strong phases [28,29], and explore the dynamics of D 0 −D 0 mixing [1,23,30]. It is worth noting that the narrow widths of the ω and φ signals allow for a straightforward signal extraction in the process D 0 → ωφ, which is unlike other D 0 → V V decay modes, such as D 0 → ρ 0 ρ 0 , which require complicated and model dependent analyses of multi-body decays.
In this Letter, we present the first measurement of D 0 → ωφ using a ψ(3770) data sample corresponding to an integrated luminosity of 2.93 fb −1 collected by the BESIII detector [31]. The measurement is performed using a single-tag technique, where only one D 0 meson in the ψ(3770) → D 0D0 decays is reconstructed. Thus, the BF of D 0 → ωφ is calculated by where N sig is the signal yield extracted from data, N D 0D0 = (10597 ± 28 ± 89) × 10 3 is the total number of ψ(3770) → D 0D0 decays quoted from Ref. [32], is the detection efficiency, and B sub is the product of BFs for the intermediate-state decays.
In the decay D 0 → ωφ, polarization amplitudes can be extracted from angular distributions. The differential decay width is given by where θ ω is the angle between p ω π + × p ω π − and −p ω D 0 in the ω rest frame, and θ K is the angle between p φ K − and −p φ D 0 in the φ rest frame. Here, p ω π + , p ω π − , p φ K − , and p ω/φ D 0 are the momenta of the π + , π − , K − and D 0 , in the rest frame of either the ω or φ meson, respectively. By integrating over cos θ ω or cos θ K from −1 to +1, Eq. (2) is simplified to be where θ can be either θ ω or θ K . A detailed illustration of the decay topology, which shows the definitions of the angles can be found in Fig. 1 in Ref. [33]. A detailed description of the design and performance of the BESIII detector can be found in Ref. [34]. A Monte Carlo (MC) simulation tool based on Geant4 [35] is implemented, in which the e + e − annihilation is simulated with the KKMC generator [36] incorporating the effects of beam-energy spread and initial-state-radiation (ISR). An inclusive MC sample, composed of the production of DD pairs, the non-DD decays of the ψ(3770), continuum processes of e + e − → qq (q = u, d, s), and the ISR production of the J/ψ and ψ(3686) states, is used to study the potential background. In the MC sample, the known decay modes are generated with EvtGen [37] using BFs from the Particle Data Group (PDG) [38], and the remaining unknown decays are generated with Lundcharm [39]. The signal sample of D 0 → φω decays is modeled by a pseudoscalar meson decaying into two vector mesons with transverse polarization using Evt-Gen [37].
The φ and ω candidates are reconstructed from their dominant decays φ → K + K − and ω → π + π − π 0 , respectively, where the π 0 is identified by a photon pair. The charged tracks must be within the main drift chamber (MDC) acceptance region by requiring the polar angle | cos θ| < 0.93, and must originate from the interaction point (IP) with a distance of closest approach within ±1 cm in the plane perpendicular to the beam and ±10 cm along the beam direction. Particle identification (PID) is performed by requiring L π > L K and L K > L π for the π ± and K ± candidates, respectively, where L π and L K are the likelihoods for the pion and kaon hypotheses calculated by combining the time-of-flight (TOF) information from the TOF detector and the dE/dx information from the MDC.
Photon candidates are selected from neutral showers deposited in the electromagnetic calorimeter (EMC) with energies larger than 25 MeV in the barrel region (| cos θ| < 0.80) and 50 MeV in the end-cap regions (0.86 < | cos θ| < 0.92). The EMC timing is required to be within 700 ns relative to the event start time to suppress electronic noise and deposited energy unrelated to the collision events. Furthermore, a photon candidate is required to be at least 10 • away from any charged tracks to avoid any overlap between them. A π 0 candidate is formed by a photon pair with invariant mass within (0.115, 0.150) GeV/c 2 . To improve the resolution, a kinematic fit is imposed on the selected photon pair by constraining their invariant mass at the nominal π 0 mass [38], and the resultant kinematic variables are used in the subsequent analysis.
To identify the D 0 signal, the energy difference ∆E = E D − E beam and the beam-constrained mass M BC = E 2 beam /c 4 − p 2 D /c 2 are calculated, where E beam is the beam energy, and E D (p D ) is the reconstructed energy (momentum) of the D 0 candidate in the e + e − center-ofmass system. The D 0 signal peaks around zero in the ∆E distribution and around the nominal D 0 mass (m D ) in the M BC distribution. The D 0 → ωφ signal is reconstructed from all possible π + π − π 0 K + K − combinations. If there is more than one combination, the one with a minimum value of |∆E| is selected. A D 0 candidate is required to satisfy M BC > 1.84 GeV/c 2 and −0.03 < ∆E < 0.02 GeV. The ∆E requirement corresponds to an interval of 4 standard deviations from the peak position. The asymmetric boundaries stem from the photon energy detection in the EMC. A prominent peak corresponding to the K 0 S in the M π + π − distribution, arising from the background process D 0 → K 0 S + anything, is rejected by removing the mass range (0.490, 0.503) GeV/c 2 . Figure 1 shows the M BC distribution for the surviving events of data and the background predictions from various MC samples with the K + K − invariant mass M K + K − < 1.05 GeV/c 2 and the π + π − π 0 invariant mass M π + π − π 0 > 0.65 GeV/c 2 , where the clear peak around m D in data refers to the signal of D 0 → π + π − π 0 K + K − . Fit to the MBC distribution of the candidate events for D 0 → π + π − π 0 K + K − . Black dots with error bars are data, dashed cyan curve for combinatorial background, long dashed-dotted pink curve for the D 0 signal, the solid blue curve for the total fit, and shadow histograms for the non-D 0 background predictions from various MC samples. The two black and two pink (red) arrows represent the MBC signal and low (high)-sideband regions, respectively.
The D 0 → ωφ signal is evident in Fig. 2, where the distribution of M π + π − π 0 versus M K + K − as well as their corresponding projection plots are shown for events in the  2. (Top) the distributions of M K + K − versus M π + π − π 0 in the MBC signal (left) and sideband (right) regions, and (middle and bottom) the corresponding 1-D projection plots of M π + π − π 0 (left) and M K + K − (right). In the projection plots, the black dots with error bars are data, the solid blue, dashed red, dotted green, dashed-dotted blue, and long dashed-dotted cyan curves represent total fit results, SIG-NAL, BKGI, BKGII, and BKGIII, respectively.
To extract the signal yield, a two-dimensional (2D) unbinned maximum likelihood fit is performed on the M π + π − π 0 versus M K + K − distributions. This fit is performed simultaneously in both the M BC signal and sideband regions, where the sideband events are used to constrain the background from non-D 0 decays. The fit includes a signal component, SIGNAL, which has both ω and φ intermediate states, and three backgrounds, BKGI, BKGII, and BKGIII. The BKGI (BKGII) contains only the ω (φ) intermediate state, and BKGIII includes neither the ω nor φ intermediate states. It is worth noting that the above four components may exist in both D 0 and non-D 0 decays. The yield of the signal D 0 → ωφ is extracted from the M BC signal region by subtracting the contribution from non-D 0 decays estimated from the M BC sideband region.
The SIGNAL is described by a distribution obtained from a 2D kernel estimation [40] of the unbinned signal MC samples. BKGI is parameterized with the product of a distribution obtained from a 1D kernel estimation [40] of the ω signal MC for the M π + π − π 0 distribution and a reversed ARGUS function [41] defined by the formula of Eq.(4) in Ref. [42] for the M K + K − distribution. Vice versa, BKGII is described with the product of an ARGUS function for the M π + π − π 0 distribution and a distribution obtained from a 1D kernel estimation of the φ signal MC for the M K + K − distribution. BKGIII is the product of an ARGUS function for the M π + π − π 0 distribution and a reversed ARGUS function for the M K + K − distribution. To compensate for the resolution difference between data and simulation, the shapes derived from simulation are convolved with (1D or 2D) Gaussian functions, which share the same parameters between different fit components and these parameters are floated during the fit. The endpoints of the ARGUS functions are fixed to the corresponding threshold values of (m D −m φ ) and 2m K ± , respectively, where m φ (m K ± ) is the nominal mass of the φ (K ± ) meson [38].
Detailed MC studies show that the non-peaking background shapes in the M K + K − distributions are identical in both the M BC signal and sideband regions, but slightly different for M π + π − π 0 distributions due to the threshold effect of kinematics. Thus, the reversed ARGUS parameterizations of the M K + K − distributions share the same parameters in both M BC signal and sideband regions, but no constraint is implemented for the ARGUS functions for the M π + π − π 0 distributions in different M BC regions. We float SIGNAL, BKGI, BKGII, and BKGIII components in both M BC signal and sideband regions during the fit. The final signal yield is also constrained to be N SG = N sig + f · N SB , where N SG and N SB are the numbers of the SIGNAL component in the M BC signal and sideband regions, respectively, as shown in Fig. 2. The factor f is the ratio of the corresponding yields from the non-D 0 decay in the M BC signal and sideband regions, and its value is determined to be (44.3 ± 0.9)% by fitting the M BC distribution, as shown in Fig. 1. In this fit, the D 0 signal is described by the simulated signal shape convolved with a Gaussian function while the non-D 0 background by an ARGUS function [41], where we fix signal shape and float rest of the other parameters during the fit. The 2D simultaneous fit yields N sig = 195.9±29.1, which includes the uncertainties from N SB and N SG . The detection efficiency is calculated to be (3.32 ± 0.04)% by the same 2D simultaneous fit approach with an inclusive MC sample, which is a mixture of the signal MC sample generated by considering the polarization of D 0 → ωφ as discussed below, and various backgrounds. The BF of D 0 → ωφ is determined to be (6.48 ± 0.96 ± 0.40) × 10 −4 according to Eq. (1), where the first and second uncertainties are statistical and sys-tematic, respectively. The corresponding significance is 6.3 σ calculated by −2 ln(L 0 /L max ) including both statistical and systematic unncertainties, where L max and L 0 are the likelihood values for the nominal fit and the alternative fit with zero signal assumption, respectively. Different contributions to the systematic uncertainty will be described later.
To study the polarization in the D 0 → ωφ decay, the efficiency-corrected signal yields are evaluated in five equal bins of | cos θ ω | and | cos θ K | as shown in Fig. 3. Here, we extract the signal yield in each bin using a procedure similar to the 2D simultaneous fit approach discussed above. The corresponding detection efficiency is obtained using a simulated signal sample generated uniformly over phase space (PHSP). A joint χ 2 fit on the | cos θ ω | and | cos θ K | distributions of data is performed with Eq. (3), where f L is floated between [−1, 1]. The fit yields f L = 0.00 ± 0.10 ± 0.08, which corresponds to f L < 0.24 at 95% confidence level computed by integrating the likelihood versus f L curve from zero to 95% of the total curve after including the systematic uncertainty as described below. This result indicates that the vector mesons are transversely polarized in the D 0 → ωφ decay. According to Eq. (1), the systematic uncertainties for the BF measurement include those from the reconstruction efficiency, MC modeling, signal yield, number of D 0D0 events, and the BFs of the intermediate-state decays. The systematic uncertainties associated with the reconstruction efficiency include tracking and PID of the charged tracks, π 0 reconstruction, ∆E requirement, and K 0 S veto. The systematic uncertainties from the MC modeling include those from the MC statistics, ω → π + π − π 0 modeling, quantum correlation (QC) effect, and the longitudinal polarization fraction f L , which varies from f L = 0 to 1σ upper bound uncertainty of 0.13. The systematic uncertainty due to the 2D simultaneous fit includes those from signal and background probabil-ity density functions (PDFs), the ratio of background between the M BC signal and sideband regions (f ), and the fit bias. All of the above systematic uncertainties are estimated with different approaches. Refer to the supplementary material [33] for details. The uncertainties of N D 0D0 and the BFs of the intermediate-state decays are from Ref. [32] and PDG [38], respectively. The total systematic uncertainty is 6.2% calculated by summing all individual uncertainties quadratically and assuming them to be independent.
The systematic uncertainty for the f L measurement includes those from MC modeling, M BC signal region, background fraction f , different bin size of cos θ ω,K , and signal and background PDFs. We replace the PHSP signal sample with a MC sample generated under the hypothesis of transverse polarization to evaluate the efficiency in each bin of the cos θ ω and cos θ K distributions. We also extract the signal yields with the alternative M BC signal region enlarged by 2 MeV/c 2 , background fraction f , different bin size of cos θ ω and cos θ K , and signal and background PDFs as done for the BF measurement. A joint χ 2 fit is performed to each set of the efficiency corrected signal yields versus cos θ ω and cos θ K data, and the resultant change in f L is considered as a systematic uncertainty. The total systematic uncertainty is 0.08 calculated as the quadratic sum of the individual ones.
In summary, the decay D 0 → ωφ is observed for the first time with a significance of 6.3σ by analyzing the ψ(3770) data taken by the BESIII experiment, corresponding to an integrated luminosity of 2.93 fb −1 . The measured BF is (6.48 ± 0.96 ± 0.40) × 10 −4 , which is consistent with the factorization model predictions [1,22], but inconsistent with predictions based on SU(3) symmetry with nonet symmetry [22], the factorization-assisted topological-amplitude method [23], and the heavy quark effective Lagrangian and chiral perturbation theory [24]. Our angular distribution studies reveal that the ω and φ in the decay D 0 → ωφ are transversely polarized, which is the same as that observed in the decay D 0 →K * 0 ρ 0 , but contradicts predictions from the naive factorization [1] and Lorentz invariant-based symmetry [19] models. The results challenge our understainding of the underlying dynamics in charmed meson decays, and also may help in searches for new physics.
The BESIII collaboration thanks the staff of BEPCII, the IHEP computing center and the supercomputing center of USTC for their strong support.  Figure 1 illustrates the decay topology of D 0 → ωφ as well as the definitions of θ ω and θ φ using in the polarization analysis.

II. THE DETAIL UNCERTAINTIES ASSOCIATED WITH THE RECONSTRUCTION EFFIENCY
The uncertainties associated with the reconstruction efficiency include tracking and PID of the charged tracks, π 0 reconstruction, ∆E requirement, and K 0 S veto. The uncertainty associated with the tracking efficiency is studied using a control sample of ψ(3770) → DD with hadronic D decays via a partial reconstruction method [1,2], where a small deviation between data and simulation is present for kaon tracks with momenta less than 0.35 GeV/c. The kaons from φ decay in the signal are of low momenta. Consequently, a correction factor of 1.06 for K + K − is applied in the detection efficiency, and an uncertainty of 0.5% is assigned for each kaon or pion. The correction factor is the ratio of the efficiencies of data and simulation weighted according to the kaon momentum distribution. We also utilize this control sample to compute the uncertainties associated with PID (0.5%) and π 0 reconstruction efficiency (2.0%) [3].
The uncertainty originating from the ∆E requirement is studied using a control sample of D 0 → 2(π + π − )π 0 decays, which has a similar final state as the signal except with a pion pair instead of a kaon pair. The control sample is selected by a relatively loose ∆E requirement, i.e., ∆E < 0.1 GeV, and the corresponding signal yield is extracted by fitting the M BC distribution. The nominal ∆E requirement is then implemented on the control sample, and the resultant ratio of signal yields is taken as the efficiency. The approach is implemented for both data and inclusive MC samples, and the resultant difference in the data and MC efficiencies, 1.4%, is taken as the uncertainty.
The uncertainty from the K 0 S veto is studied by varying the K 0 S mass window requirement within ±1σ, and the larger difference in the BF, 0.8%, is taken as the uncertainty.
The total uncertainties associated with the reconstruction efficiency is 3.8%, which is the quadratic sum of above individual ones.

III. THE DETAIL UNCERTAINTIES ASSOCIATED WITH THE MC MODELING
The uncertainties from the MC modeling includes those from the MC statistics (0.8%), ω → π + π − π 0 modeling, quantum correlation (QC) [4] effect, and the longitudinal polarization fraction f L . The uncertainty due to the ω → π + π − π 0 modeling is assigned to be 0.5% on the basis of two MC samples generated with two different models [5,6]. From the analysis, the decay D 0 → ωφ appears to be transversely polarized, thus it is a mixture of CP -even and CP -odd components. The uncertainties associated with the polarization is studied by an alternative signal MC sample generated with 1σ upper bound uncertainty, f L = 0.13, and the resultant change in the efficiency, 3.2%, is taken as the uncertainty.
The total uncertainties associated with the MC modeling is 3.3%, which is the quadratic sum of above individual ones.

IV. THE DETAIL UNCERTAINTIES ASSOCIATED WITH 2D SIMULTANEOUS FITS
The systematic uncertainty due to the 2D simultaneous fit includes those from signal and background probability density functions (PDFs), the ratio of background between the M BC signal and sideband regions (f ), and the fit bias. The uncertainty arising from the signal PDF, 1.2%, is evaluated with an alternative fit, in which the signal PDFs are described using a different non-parameterized modeling of the simulated shape, convolved with a Gaussian function. The uncertainty of the background PDF, 0.4%, is determined by replacing the ARGUS function [7] with a modified one as used in Ref. [8]. The uncertainty from f is 0.1%, evaluated by varying its value within 1 σ when calculating the signal yield. The uncertainty due to the choice of the M BC signal region is evaluated to be 2.7% by enlarging its region by 2 MeV/c 2 , which is the resolution of the M BC distribution. The fit bias, 1.0%, is estimated with a large number of pseudo-experiments. Each pseudo-experiment sample is a composition of the signal generated according to the signal PDF and background expectations from the inclusive MC sample. The resultant pull distribution for the BF is consistent with a normal distribution, and we consider the average fit bias as the uncertainty.
The total uncertainty associated with the 2D simultaneous fits is 3.2%, which is the quadratic sum of above individual ones.