Measurement of the $\omega \to \pi^+ \pi^- \pi^0$ Dalitz plot distribution

Using the production reactions $pd\to {}^3\mbox{He}\,\omega$ and $pp\to pp\omega$, the Dalitz plot distribution for the $\omega \to \pi^+ \pi^- \pi^0$ decay is studied with the WASA detector at COSY, based on a combined data sample of $ (4.408\pm 0.042) \times 10^4$ events. The Dalitz plot density is parametrised by a product of the $P$-wave phase space and a polynomial expansion in the normalised polar Dalitz plot variables $Z$ and $\phi$. For the first time, a deviation from pure $P$-wave phase space is observed with a significance of $4.1\sigma$. The deviation is parametrised by a linear term $1+2\alpha Z$, with $\alpha$ determined to be $+0.147\pm0.036$, consistent with the expectations of $\rho$-meson-type final-state interactions of the $P$-wave pion pairs.


Introduction
The present work and foreseeable follow-ups are based on two motivations: 1. To check and improve on our understanding of the importance of hadronic final-state interactions for the structure and decays of hadrons. 2. To improve the Standard Model prediction for the gyromagnetic ratio of the muon [1]. The present work accomplishes the first task concerning the Dalitz decay of the ω meson into three pions; it also constitutes a significant step forward towards providing improved hadronic input for the second task. In the following we shall discuss the two tasks in more detail.
The ω-meson resonance was discovered in 1961 [2]. Its main decay branch is ω → π + π − π 0 , with a branching ratio of BR = (89.2 ± 0.7)%. By now it is well established that the ω meson has spin-parity J P = 1 − [3]. As a consequence, the combination of Bose, isospin, and parity symmetry of the strong interaction demands that for the decay ω → π + π − π 0 every pion pair is in a state of odd relative orbital angular momentum. Given the limited phase space of the decay one can safely assume the Pwave to be the dominant partial wave. 8 If a pion pair is in a Pwave state, then the third pion will be in P-wave state relative to the pair. This "P-wave phase space" distribution has been confirmed experimentally. Historically the P-wave dominance of the decay has actually been used to pin down the quantum numbers of the ω meson [5,6,7,8].
If the pions, once produced in the decay, did not interact further, then solely the P-wave phase space would shape the Dalitz plot of the decay ω → π + π − π 0 . However, a pion pair in a P-wave shows a very strong final-state interaction. The two-pion P-wave phase shift is dominated by the ρ meson, and is now known very accurately [9,10,11]; this is essential in particular for theoretical studies of these decays using dispersion theory, which use the phase shifts as input directly [4,12]. In the similar decay φ → π + π − π 0 one can see the ρ meson as a resonance in the corresponding Dalitz plot [13,14]. In the decay ω → π + π − π 0 there is not enough energy for a pion pair to reach the ρ resonance mass; yet already for the available invariant masses the two-pion P-wave phase shift is significantly different from zero. In fact, every theoretical approach * Corresponding authors Email addresses: lena.heijkenskjold@physics.uu.se (L. Heijkenskjöld), siddhesh.sawant@iitb.ac.in (S. Sawant) 1 8 Genuine F-wave corrections have been modelled theoretically, and found to be tiny [4]. that deals with the decay ω → π + π − π 0 includes this non-trivial phase shift and/or the ρ meson in one way or the other; see, e.g., Refs. [15,16,17,18,19,4,20,12] and references therein. In practice this leads to an increase of strength towards the boundaries of the Dalitz plot, superimposed with the pure P-wave phase space, which drops towards the boundaries. This increase of strength is on the level of about 20% [19,4].
The prediction of the presence of final-state interactions leading to an increase of strength towards the Dalitz-plot boundaries ought to be tested experimentally. Interestingly this has not been achieved so far. The highest statistics of a dedicated ω → π + π − π 0 Dalitz plot measurement from 1966 had 4208±75 signal events [21]. Due to the limited statistics, fits with a pure P-wave phase space could not be distinguished from a distortion by the final-state interactions, i.e. by intermediate ρπ states. Surprisingly there were no further dedicated Dalitz plot studies of the ω → π + π − π 0 decay. In the present work we will reveal that the universal final-state interactions of the pion pairs are indeed present in the ω Dalitz decay. In the analysis presented here we have produced an acceptance-corrected Dalitz plot and extracted experimental values for parameters describing the density distribution. This constitutes the first task spelled out in the beginning of this introduction.
The second motivation for a precise measurement of the ω → π + π − π 0 Dalitz plot consists in improving hadronic input for the theoretical assessment of the hadronic light-by-light scattering contribution to the anomalous magnetic moment of the muon. The largest individual contribution is given by the lightest hadronic intermediate state, the so-called π 0 pole term, whose strength is determined by the corresponding singly and doubly virtual transition form factors. One of the few possibilities to gain experimental access to the doubly virtual π 0 transition form factor with high precision consists in studying vector meson conversion decays, in particular ω → π 0 + −which is intimately linked, through dispersion relations, to the ω → π + π − π 0 decay amplitude [22,23,12]. However, these theoretical descriptions of the ω transition form factor (see also Ref. [24]) fail to describe the very precise data on ω → π 0 µ + µ − taken by the NA60 collaboration [25,26], which may violate very fundamental theoretical bounds [27,28]. In all these studies, the ω → π + π − π 0 decay amplitude is a potential loose end, as it is so far only theoretically modelled, not experimentally tested. In combination with the precisely measured φ → π + π − π 0 Dalitz plot information, it could be used to further constrain the amplitude analysis of e + e − → π + π − π 0 , and hence the π 0 transition form factor in a wider range [29].
Recently, an easy-to-use polynomial parametrisation of the ω → π + π − π 0 Dalitz plot distribution has been suggested [4], as a generalisation of the commonly used one for the decay η → 3π 0 (which has a similar crossing symmetry). In the present work we utilise the same parametrisation and compare to recent theoretical approaches [4,20,12] that have provided predictions for the corresponding Dalitz plot parameters. Refs. [4,12] are based on dispersion theory: both employ the pion-pion P-wave scattering phase shift as input and describe rescattering of all three final-state pions consistently to all orders. The Dalitz plot parameters cited from Ref. [12] in Table 1 α × 10 3 β × 10 3 Uppsala [20] 202 -Bonn [4] 84 . . . 102 -JPAC [12] 94 -  [20,4,12], where at most two parameters were used in the fit.
refer to a formalism with only the overall normalisation of the amplitude fixed from experiment, hence the energy dependence is fully predicted; the ranges as given in Ref. [4] include (in addition to variations in the phase shift input etc.) potential deviations from such a pure prediction as gleaned from the analysis of the analogous φ → π + π − π 0 Dalitz plot [13]. Ref. [20] is based on an effective Lagrangian for the lightest pseudoscalar and vector mesons. The strength of the initial ω-ρ-π interaction is fitted to the decay width of ω → π + π − π 0 and crosschecked with the decay width of ω → π 0 γ. The Lagrangian provides the kernel for a Bethe-Salpeter equation that generates the two-pion rescattering. In contrast to the dispersive approaches, crossed-channel rescattering of the three-pion system is not included. One way to describe a three-particle Dalitz decay distribution is to use invariant masses of particle pairs [3]. This would be particularly useful for reading off resonance masses if the decay was mediated by one or several resonances. However, there are no intermediate resonances in the kinematically accessible energy range of the decay ω → π + π − π 0 that would be compatible with the symmetries of the strong interaction. For our case of interest, we first split off the P-wave phase space (see, e.g., Refs. [19,4]) and parametrise the rest by a polynomial distribution following [4]. Denoting the polarisation vector of the ω meson by (P ω , λ ω ) and the momenta of the outgoing pions by P + , P − , and P 0 , we start with the most general matrix element compatible with the symmetries, The dynamics of the final-state interactions is encoded in the scalar function F [19,4]. After summation over the helicity λ ω of the ω meson one obtains a Dalitz plot distribution proportional to with the pure P-wave phase-space distribution Note that for the P term we can account for "kinematic" isospin violations due to the difference between the masses of the uncharged and charged pions, m 0 and m ± , respectively. For the remaining distribution F , which covers the dynamics of the final-state interaction, we ignore isospin breaking effects.
The quantity F and therefore also |F | 2 would be a constant if there were no final-state interactions between the produced pions. In reality |F | 2 is not a constant, but relatively flat. Instead of parametrising |F | 2 by invariant masses of pion pairs we follow Ref. [4] and utilise normalised variables X and Y, which have their origin at the centre of the Dalitz plot. They are defined by with Here T i are the kinetic energies of the pions in the ω rest frame (centre-of-mass frame of the three-pion system). Finally one introduces polar coordinates by The expansion for |F | 2 , valid in the isospin limit, reads where N is a normalisation constant and G contains the expansion in Z and φ [4]: The Dalitz plot distribution can then be fitted using this formula to extract the "Dalitz plot parameters" α, β, γ, . . . . The fit results to the theory predictions of Refs. [4,20,12], if Eq. (8) is truncated at order Z (one parameter fit) or at order Z 3/2 (two parameter fit), are shown in Table 1. The reproduction of the theoretical Dalitz plot distributions is improved significantly in all cases when including the term ∝ β.

The experiment
The experimental data was collected using the WASA setup, where the ω was produced in the pd → 3 He ω reaction and in the pp → ppω reaction. The WASA detector [30,31] is an internal target experiment at the Cooler Synchrotron (COSY) storage ring, Forschungszentrum Jülich, Germany. The COSY proton beam interacts with an internal target consisting of small pellets of frozen hydrogen or deuterium (diameter ∼ 35 µm).
The WASA detector consists of a Central Detector (CD) and a Forward Detector (FD), covering scattering angles of 20 • -169 • and 3 • -18 • , respectively. The CD is used to measure decay products of the mesons. A cylindrical straw chamber (MDC) is placed in a magnetic field of 1 T, provided by a superconducting solenoid. The electromagnetic calorimeter (SEC) consists of 1012 CsI(Na) crystals which are read out by photomultipliers. A plastic scintillator barrel (PSB) is placed between the MDC and the SEC, allowing particle identification and accurate timing for charged particles. The FD consists of thirteen layers of plastic scintillators for energy and time determination and a straw tube tracker providing a precise track direction. When the ω mesons were produced using the pd → 3 He ω reaction, two different proton kinetic energies were used: T A = 1.450 GeV and T B = 1.500 GeV. The cross section of the reaction is 84(10) nb at the lower energy [32] and was studied previously by the CELSIUS/WASA collaboration. Triggers select events with at least one track in the FD with a high energy deposit in the thin plastic scintillator layers. This condition allows for an efficient selection of 3 He ions and provides an unbiased data sample of ω meson decays. The proton beam energy was chosen so that the 3 He produced in the pd → 3 He ω reaction stops in the second thick scintillator layer of the FD. The correlation plot ∆E − ∆E from a thin layer and the first thick layer of the FD is shown in Fig. 1(top). The band corresponding to the 3 He ion is well separated from the bands for other particles and allows a clear identification of 3 He. The 3 He from the reaction of interest has kinetic energies up to 700 MeV and scattering angles ranging from 0 • to 10 • .
The pp → ppω experiment was performed at T C = 2.063 GeV beam kinetic energy, corresponding to 60 MeV centre-of-mass excess energy and cross section 5.7 µb [33]. In the pp collision experiment, the selected events were required, at trigger level, to contain at least two tracks reaching the second thick layer of the plastic scintillators in the FD, at least two hits in the PSB, and at least one cluster in the SEC. In the offline analysis, pairs of tracks corresponding to the ∆E − ∆E proton bands, shown in Fig. 1(bottom), in different thick layers of the FD are selected as proton pair candidates.
For the particles measured in the CD, a common analysis procedure is used for all three data sets. Events are selected if they contain at least one pair of opposite charge particle tracks in the MDC with scattering angles greater than 30 • and at least two neutral clusters with energy deposit above 20 MeV in the SEC. Relative time between the tracks is checked to minimise pile ups. The charged particle tracks are assigned the charged pion mass. Possible combinations of the charged and neutral particle tracks for the selected events are tested using a constrained kinematic fit assuming the conservation of energy and momentum with the pd → 3 He π + π − γγ or pp → ppπ + π − γγ hypothesis, respectively. The events with p-values less than 0.05 are rejected. For the case when more than one track combination fulfils the selection, the combination giving larger pvalue is selected. Finally, further background suppression is achieved by applying a kinematic fit with the contending hypothesis pd → 3 He π + π − or pp → ppπ + π − , respectively. If the resulting p-value is larger than for the first fit, the event is rejected.
The missing mass distributions, MM( 3 He) and MM(pp), for the three data sets are shown in Fig. 2. The missing masses, calculated from the variables corrected by the kinematic fit, are equivalent to the invariant mass of the π + π − γγ system. The observed ω peak position is shifted from the nominal value of the ω mass by +0.7 MeV for MM( 3 He) in the two pd data sets and by +1.1 MeV for MM(pp). The observed shifts correspond to deviations from the nominal beam energy by 0.55 MeV and 0.75 MeV, respectively, which is well within the uncertainty of the absolute energy scale of COSY. To reproduce the experimental ω peak position, the missing mass distributions were shifted accordingly. To also reach agreement between experiment and simulation for the width of the ω peak, the resolution from simulated detector responses were adjusted.
Both the background shape and the ω peak content are fitted simultaneously to the experimental distribution using the following fit function: (9) where µ = MM( 3 He) or MM(pp). H ω (µ) and H 3π (µ) represent reconstructed distributions of simulated signal and background and correspond to events that have passed through the same analysis steps as the experimental data. H ω (µ) is normalised such that the fit gives directly the number of signal events, N S , and the related error. The other parameters fitted are a 0 , a 1 , a 2 , and a 3 (in case of pd data a 3 is set to 0). The range in µ used for the fit is

Dalitz plot
The Dalitz plot density is represented using a twodimensional histogram in the Z and φ variables, defined in Eq. (6). The size of the selected bins is determined by the experimental resolution of Z and φ and the statistics of the collected data sample. The number of events in each bin should be sufficient for determining the signal yield and to carry out a χ 2 fit of the Dalitz plot density parametrisation. The φ variable range [−π, π] is divided into six bins to preserve the threefold isospin symmetry and to be sensitive to a possible sin 3φ dependence. The Z variable range [0, 1] is also divided into six bins. Only the 21 bins fully contained inside the kinematic limits of the decay are used. Figure 3 introduces the bin numbering used for the presentation of the results.
A small shift of the Dalitz plot along the Y axis is due to the mass difference between the neutral and charged pions. It is most visible in Fig. 3 when comparing the regions at φ = π/2 to the ones at φ = −π/6 and −2π/3. The picture shows also seven sectors I-VII that are used to test the consistency of the fit results.
For each Dalitz plot bin, the experimental missing mass distribution is constructed and the number of entries in the ω peak is extracted by fitting a simulated ω → π + π − π 0 signal along with background contributions using Eq. (9).
Since the P-wave distribution reproduces the general features of the ω → π + π − π 0 Dalitz plot very well and the deviations are expected to be small, the efficiency correction is obtained using signal simulation with the P-wave. The efficiency, i , is extracted using the ratio i = N i /N G i . N G i is the number of events with generated kinematic variables corresponding to bin i in the Dalitz plot. N i is the content of the bin i when the reconstructed values of the kinematic variables are used for events passing all analysis steps. The extracted efficiencies for the three data sets are shown in Fig. 4. For the pd data sets the overall efficiency is 11%, while for pp it is 20 times lower. This low efficiency for the pp data sample has the following well-understood causes. The colour plot shows the kinematically allowed region of the ω → π + π − π 0 reaction with ω nominal mass as well as the density distribution from P-wave dynamics.
In most cases, the two fast proton tracks deposit only a fraction of their kinetic energy in the detector, leading to a lower precision of the kinetic energy determination and an asymmetric resolution function. The events from the tails will likely be rejected by the kinematic fit procedure. On the other hand for the pd → 3 He ω reaction, there is only one doubly charged 3 He stopping in the detector. Another cause for the low efficiency is the larger centre-of-mass velocity in the pp reaction, which decreases the average emission angle for decay particles, in particular for the charged pions. The pions will be more often emitted at angles below 30 • and will therefore be rejected in the analysis procedure. The Dalitz plot parameters (α, β, . . .) and normalisation factors for the three data sets (N A , N B , N C ) are determined by minimising the following χ 2 = χ 2 A + χ 2 B + χ 2 C function, where N i andσ i are the efficiency corrected experimental Dalitz plot bin content and error, respectively. H i is given by an inte- The parametrisation procedure of the Dalitz plot is tested using 10 6 signal events simulated with P-wave phase space only (i.e. G = 1) and without detector smearing. The extracted parameters are found to be consistent with zero and therefore the procedure does not introduce any bias at the present statistical accuracy.
The three independent data sets and the Dalitz plot symmetries allow for detailed checks of the experimental efficiency and the background subtraction procedure since the background distributions and efficiencies are different in the corresponding bins.
The method of background subtraction for the missing mass µ distributions is tested by preparing simulated distributions after full detector reconstruction, consisting of a sum of π + π − π 0 production background events and the ω signal generated using a P-wave phase space distribution. The background is obtained from the a 0 + a 1 µ + a 2 µ 2 + a 3 µ 3 × H 3π (µ) distributions with a i determined from the fits using Eq. (9) and by setting the average signal-to-background ratio to be the same as in the experimental data. The generated µ distributions with the number of events similar as in the experiment are then subjected to the same background subtraction as the experimental data. The combined fit of the Dalitz plot parametrisation to the samples A and B with only the α parameter gives α = (10 ± 35) · 10 −3 and χ 2 = 36/39. For set C α = (25 ± 59) · 10 −3 and χ 2 = 24/19. Therefore the background subtraction procedure does not introduce any experimental bias.
One can also study the bias and accuracy of the efficiency determination by considering X-or Y-dependent corrections for the efficiency: where ξ A,... , ζ A,... are single parameters for each data set. Fits to separate data sets show that all ζ coefficients are consistent with zero and do not change the value of the χ 2 . On the other hand, ξ B and ξ C were found to significantly deviate from zero, although with opposite signs. Applying these two corrections to the efficiency bin number 5 10 15 20 Normalised bin content     plot bin contents are provided in Table 2.
The extracted Dalitz plot parameters and goodness of fit for each data set separately are reported in Table 3. There is a significant decrease of the χ 2 value when including the α parameter and the results from the three data sets are consistent. The results of the fits for all data sets combined are given in Table 4. The p-value significantly improves after including the α parameter, while inclusion of an additional parameter does not improve the p-value any further. The efficiency corrected Dalitz plots for the three data sets are shown in Fig. 5, where they are compared to the P-wave distribution as well as the fit with the α parameter. We consider the result with the α parameter our main finding. The difference between the results from the first and second fits in Table 4 indicates the onset of dynamics in the reaction on top of the P-wave phase space distribution. This follows the expected behaviour of an increase towards the edges of phase space due to the attractive ππ final-state interaction (approaching the ρ resonance), yielding a positive value for the α parameter. Figure 6 shows a test of data consistency for the same Dalitz plot sectors, marked with Roman numerals in Fig. 3. Arithmetical averages of the normalised residuals with respect to the α parameter fit from bins corresponding to the same sectors are calculated for the separate data sets and for all three data sets combined. The error bars correspond to the calculated root mean square values, which are expected to be 1 for a random data sample with correctly estimated uncertainties.
The conclusion of the checks for the systematic effects is that the accuracy is dominated by statistic uncertainty.

Summary and discussion
For the first time a deviation from a pure P-wave distribution in ω → π + π − π 0 is observed and quantified by the determination of the parameter α = (147 ± 36) · 10 −3 , i.e. a positive value with 4.1σ significance. Figure 7 compares the experimental result of the α parameter to theoretical predictions. The systematic effects were studied by comparing three data sets using two production reactions, which differ significantly in resolution and acceptance. The chosen Z, φ parametrisation together with isospin symmetry allows for more tests of systematic effects, and the precision of the result is found to be dominated by the statistical uncertainty.
(MPD), co-financed by the European Union within the European Regional Development Fund. We gratefully acknowledge the support given by the Swedish Research Council, the Knut and Alice Wallenberg Foundation, and the Forschungszentrum Jülich FFE Funding Program. This work is based on the PhD theses of Lena Heijkenskjöld and Siddhesh Sawant.