Study of the $e^{+}e^{-}\to\pi^{+}\pi^{-}\omega$ process at center-of-mass energies between 4.0 and 4.6 GeV

The $\pi^{+}\pi^{-}\omega$ final sate is produced from the $e^{+}e^{-}$ annihilation into a virtual photon, followed by hadronization into the $\pi^{+}(p_{1})\pi^{-}(p_{2})\omega(p_{3})$ final sate, where $p_{i}$ ($i=$ $1$, $2$, $3$) denote particle momenta after the kinematic fit. The $\pi^{+}\pi^{-}\omega$ final state may be produced non-resonantly, or via an intermediate resonance and subsequent decay; the possible resonance diagrams are shown in Fig. 2.

The amplitudes for these diagrams are constructed using the helicity formalism. Taking the first diagram in Fig. 2 as an example, one may define the helicity rotation angles as in Fig. 3. For resonance $R_{1}$ the polar angle ($\theta_{[12]}^{[123]}$) is defined as the angle spanned between the $R_{1}$ momentum and the positron beam direction, the azimuthal angle ($\phi_{[12]}^{[123]}$) is the angle between the $R_{1}$ production plane formed by the $R_{1}$ momentum and the $z$ axis and the plane formed by the $x$ and $z$ axes. Here, $xyz$ denotes the laboratory coordinates. The helicity amplitude for $\gamma^{*}\to R_{1}(\lambda_{R})\omega(\lambda_{3})$ is denoted by $F^{\gamma^{*}}_{\lambda_{R},\lambda_{3}}$ with specified helicity $\lambda_{R}$ and $\lambda_{3}$. For the $R_{1}\to\pi^{+}\pi^{-}$ decay, the azimuthal angle ($\phi_{[1]}^{[12]}$) is defined as the angle between the $R_{1}$ production plane and its decay plane, formed by the momenta of $\pi^{+}\pi^{-}$ from $R_{1}$. After boosting the two pion momenta to the $R_{1}$ rest frame, they are still located in the same decay plane. The polar angle ($\theta_{[1]}^{[12]}$) for $\pi^{+}$ is defined as the angle between the $\pi^{+}$ and $R_{1}$ momenta in the $R_{1}$ rest frame. The helicity amplitude of this decay is denoted by $F^{R_{1}}_{0,0}$. Helicity angles for the processes (b) and (c) are defined analogously. Table 1 summarizes the helicity angles and amplitudes for the three processes.

The decay amplitude for the process (a) is

where $D^{J}_{m,\lambda}(\phi,\theta,0)$ is the Wigner $D$-function, $J$ is the spin quantum number of resonance $R_{1}$, and $BW$ denotes the Breit-Wigner function.

The decay amplitude for the process (b) is

where $J$ is the spin of $R_{2}^{+}$. Since the $\omega$ helicity defined in the $R_{2}^{+}$ helicity system is different from that defined in the process (a), one needs to perform a rotation by the angles ($\theta_{3}^{\prime},\phi_{3}^{\prime}$) to align the $\omega$ helicity to coincide with that in the process (a). This issue has been addressed in the analyses lhcb ; belle and derived in detail in Ref. pingrg .

The decay amplitude for the process (c) reads

where the Wigner $D^{1}_{\lambda_{3}^{\prime\prime},\lambda_{3}}(\phi_{3}^{\prime\prime},\theta_{3}^{\prime\prime},0)$ function is used to align the $\omega$ helicity to coincide with that defined in the process (a).

For the direct three-body process $e^{+}e^{-}\to\pi^{+}\pi^{-}\omega$, the helicity amplitude is written as chung2 :

where $\mu$ is the $z$ component of the spin $J$ of the virtual photon in the helicity system, and $m(\lambda_{3})$ is the helicity value for $\gamma^{*}(\omega)$. Here, $\alpha$, $\beta$, and $\gamma$ are the Euler angles as defined in chung2 (see Fig. 4). $F_{\mu,\lambda_{3}}$ is the helicity amplitude; parity conservation requires $F_{\pm,\lambda_{3}}=-F_{\pm,-\lambda_{3}}$ and $F_{0,\lambda_{3}}=F_{0,-\lambda_{3}}$. Parity conservation also requires $F_{\pm,\lambda_{3}}(E_{i})=-F_{\pm,-\lambda_{3}}(E_{i})$ and $F_{0,\lambda_{3}}(E_{i})=F_{0,-\lambda_{3}}(E_{i})$, where $E_{i}(i=1,2,3)$ corresponds to the energy of the final state $\pi^{+}\pi^{-}\omega$.

One usually expands the helicity amplitudes in terms of the partial waves for the two-body decay in the $LS$-coupling scheme chung2 . For a spin-$J$ particle decay $J\to s+\sigma$, it follows

where $\lambda$ and $\nu$ are the helicities of two final-state particles $s$ and $\sigma$ with $\delta=\lambda-\nu$, and $g_{lS}$ is a coupling constant, $S$ is the total spin ${\bf S=s+\sigma}$, $l$ is the orbital angular momentum, $r=|{\bf r}|$, where ${\bf r}$ is the relative momentum between the two daughter particles in their mother rest frame, ${\bf r}_{0}$ corresponds to the value at the resonance’s known mass. $B_{l}(r)$ is the Blatt-Weisskopf factor chung2 , which suppresses the contributions with higher angular momentum. The Blatt-Weisskopf factors up to $l$ = $4$ are

where $d$ is a constant fixed to $3$ GeV^-1 for the meson final states lhcb .

The differential cross section is given by

where $m=\pm 1$ due to the polarization of the virtual photon being produced from $e^{+}e^{-}$ annihilation, and $d\Phi$ is the element of standard three-body PHSP. The $\Omega(\lambda_{3})=|{\boldsymbol{\varepsilon}}(\lambda_{3})\cdot({\bf q_{1}}\times{\bf q_{2}})|^{2}$ is the $\omega$ decay matrix element into the $\pi^{+}\pi^{-}\pi^{0}$ final states, where ${\boldsymbol{\varepsilon}}$ is the $\omega$ polarization vector, and $\bf q_{1}(q_{2})$ is the momentum vector for $\pi^{+}(\pi^{-})$ from the $\omega$ decay. Here we factor out the $BW$ function describing the $\omega$ line shape into the MC integration when applying the amplitude analysis to the data events.

The relative magnitudes and phases of the coupling constants are determined by an unbinned maximum likelihood fit. The joint probability density function (PDF) for the events observed in the data sample is defined as

where $p_{i}$ ($i$ = $1$, $2$, …, $5$) denotes the four-vector momenta of the final state particles, and $P_{i}$ is a probability to produce the $i$-th event. The normalized $P_{i}$ is calculated from the differential cross section

where $\sigma_{\rm MC}$ is the normalization factor which is calculated with a large MC sample as

where $N_{\rm MC}$ is the number of events retained with the same selection criteria as for data sample.

For technical reasons, rather than maximizing $\mathcal{L}$, $S=-\ln\mathcal{L}$ is minimized using the package MINUIT minuit . To subtract the contribution of background, the $\ln\mathcal{L}$ function is replaced with

where $\mathcal{L}_{\textrm{data}}$ and $\mathcal{L}_{\textrm{bkg}}$ are the joint PDFs for data and background, respectively. The background events are obtained from the $\omega$ sideband regions mentioned in Section 3.

A simultaneous fit is performed to data sets collected at different c.m. energies. The common parameters for different data samples in this fit are the masses, widths, and $\rm Flatt\acute{e}$ parameters for the resonances. The total function is taken as the sum of individual ones, $i.e.$,

The signal yield for the $i$-th resonance, $N_{i}$, can be estimated by scaling its cross section ratio $R_{i}$ to the number of net events

where $\sigma_{i}$ is the cross section for the $i$-th resonance as defined in Eq.(7), $\sigma_{\textrm{tot}}$ is the total cross section, and $N_{\textrm{obs}}$ and $N_{\textrm{bkg}}$ are the numbers of observed events and background events, respectively. In the simultaneous fit, the background events are taken from the $\omega$ sideband regions, and the number $N_{\rm bkg}$ is estimated with the background PDF with the $\omega$ signal region (see Fig. 1).

The statistical uncertainty, $\Delta N_{i}$, associated with the signal yield $N_{i}$, is estimated according to the error propagation formula using the covariance matrix, $V$, obtained in the simultaneous fit, i.e.

where ${\bf X}$ is a vector containing parameters, and ${\bf\mu}$ contains the fitted values for all parameters. The sum runs over all $N_{\textrm{pars}}$ parameters.

In the $\pi^{+}\pi^{-}$ and $\omega\pi^{\pm}$ mass spectrum, the $f_{0}(500)$, $f_{0}(980)$, $f_{2}(1270)$, $f_{0}(1370)$, $b_{1}(1235)^{\pm}$, $\rho(1450)^{\pm}$, and $\rho(1570)^{\pm}$ resonances are included in the amplitude model. The $f_{0}(980)$ line shape is parameterized by the $\rm Flatt\acute{e}$ formula:

where $\rho(s)=2k/\sqrt{s}$ and $k$ is the momentum of the $\pi$ or $K$ in the resonance rest frame, $g_{1}$ and $g_{2}/g_{1}$ are fixed to the measured values $(0.138\pm 0.010)$ $\rm GeV^{2}$ and $4.45\pm 0.25$  besiia ; besiib , respectively. $M$ is the mass of $f_{0}(980)$ taken from the PDG pdg .

For the $BW_{2}$ function of a wide resonance, e.g., $f_{0}(500)$, there are many parametrizations for the energy-dependent width besiia ; besiib , and we take the one used by the E791 Collaboration in the nominal fit,

where $m_{0}$ is the nominal mass of the resonance, and $\Gamma_{0}$ is its width. For other resonances, such as $b_{1}(1235)^{\pm}$, $f_{0}(1370)$, $f_{2}(1270)$, $\rho(1450)^{\pm}$, $\rho(1570)^{\pm}$, their line shapes are described with the $BW_{3}$ function,

where the widths are fixed to the individual PDG values pdg .

Based on the signal events in the $\pi^{+}\pi^{-}\pi^{0}$ mass spectrum, we select twelve c.m. energy points with relatively large statistics. We divide these selected points into two groups. Group A includes the data sets taken at $\sqrt{s}$ = $4.0076$, $4.1780$, $4.1890$, $4.1990$, $4.2093$, and $4.2188$ GeV, and group B includes $\sqrt{s}$ = $4.2263$, $4.2358$, $4.2439$, $4.2580$, $4.2668$, and $4.4156$ GeV. To check the significance of each resonance and determine the nominal solution, a simultaneous fit is performed to the data from a given group. In each group, the cross sections of these intermediate states are regarded to be energy-dependent, so the parameters responsible for the virtual photon $\gamma^{*}$ coupling to a given state are allowed to vary in the fit for various energy points, while the coupling constant parameters for the subsequent decay are taken as the common parameters for all energies. The conjugate modes share the same coupling constants. The masses, widths or $\rm Flatt\acute{e}$ parameters for the resonances of $f_{0}(500)$, $f_{0}(980)$, $f_{2}(1270)$, $f_{0}(1370)$, $\rho(1450)^{\pm}$, and $\rho(1570)^{\pm}$ are fixed to the measured values from PDG pdg , as given in Table 2. The mass and width of $b_{1}(1235)^{\pm}$ are floated due to large uncertainties. Then its nominal solution is fixed as the fitted result.

The significance of each intermediate state is estimated by the changes of $-2\ln\mathcal{L}$ and the number of degrees of freedom (NDF) after removing it from the simultaneous fit. We take the intermediate states with statistical significances greater than 5$\sigma$ in two groups as the nominal solution, including $f_{0}(500)$, $f_{0}(980)$, $f_{2}(1270)$, $f_{0}(1370)$, $b_{1}(1235)^{\pm}$, and $\rho(1450)^{\pm}$, as shown in Fig. 5. It is found that the contributions from $f_{0}(500)$ and $b_{1}(1235)^{\pm}$ are the most significant, as shown in the $M_{\pi^{+}\pi^{-}}$ and $M_{\omega\pi^{\pm}}$ spectra, respectively. The statistical significances for various intermediate resonances are shown in Table 2.

For the simultaneous fit, the ratios and the signal yields of various intermediate states are obtained according to Eq. (13), as shown in Tables 4 and 5. And their statistical uncertainties are determined based on Eq. (14), in which the correlation among parameters is included. With the intermediate states in the nominal solution, we perform the simultaneous fit to the data samples for groups A and B. Taking the two data samples from $\sqrt{s}$ = $4.1780$ and $4.2263$ GeV with large integrated luminosity as examples, Figs. 5 and  6 show the fit results for groups A and B, respectively.

In $e^{+}e^{-}$ collision experiments, the observed cross section, $\sigma_{\rm obs}(s)$, at a c.m. energy point $\sqrt{s}$, is related to the corresponding Born cross section, $\sigma_{0}(s)$, by the ISR factor

with

where $\Pi(\sqrt{s})$ is the vacuum polarization (VP) function. The $M_{\rm th}$ corresponds to the $\pi^{+}\pi^{-}\omega$ mass threshold, and $x$ is the effective fraction of the beam energy carried by photons emitted from the initial state, $x=\frac{2E_{\gamma}}{\sqrt{s}}$, and $E_{\gamma}$ is the energy of the ISR photons. The initial state radiative function, $W(s,x)$, which uses the QED calculation up to next to leading order in Ref wx ; alpha ; wsx ,

where

and we use the calculated results including the leptonic and hadronic parts both in the space-like and time-like region isr11 ; isr12 ; isr13 ; isr14 ; isr15 .

We use the generator model ConExc besevtgen to produce signal MC events and then iterate the Born cross-section measurement, in which the radiative function takes the result of high-order QED calculation up to the $\alpha^{2}$ accuracy alpha . The Born cross sections from the $\pi^{+}\pi^{-}\omega$ mass threshold to 4.6 GeV are used to calculate the ISR factor. The Born cross sections $\sigma_{0}(s)$ in the c.m. energy ranges of below $3.0$ GeV and $(4.0,4.6)$ GeV are taken from the measurements in Ref. babarXS and this work, respectively. In the c.m. energy interval of $(3.0,4.0)$ GeV, however, the Born cross section of $e^{+}e^{-}\to$ continuum light hadrons is described by a polynomial, and the Born cross sections for $J/\psi$ and $\psi(3686)$ are described by the function

where $\sqrt{s}$ is the c.m. energy, $J=1$ is the spin of the resonance, and the numbers of polarization states of the two incident particles are $2S_{1}+1=2$ and $2S_{2}+1=2$, respectively. The maximum momentum of the final-state channel is denoted as $k$, $\sqrt{s}_{0}$ is the c.m. energy at the resonance, and $\Gamma$ is the width of the resonance. The branching fractions of the resonance decays into the initial-state and final-state channels are denoted as $B_{\rm in}$ and $B_{\rm out}$, respectively. The cross sections are smoothed by a fit to seven Gaussian functions in various energy intervals. Since the detection efficiency is affected by the radiative correction, an iteration over the cross section is done until the latest two results become stable; specifically, when the updated Born cross sections change by less than the statistical uncertainty. The ISR correction factor for each c.m. energy point is given in Table 3.

The Born cross section at each c.m. energy is calculated by

where $N^{\rm sig}$ is the number of observed signal events, $(1+\delta^{\gamma})$ and $\frac{1}{|1-\Pi|^{2}}$ are the ISR correction and VP corrections, respectively. The factors $Br(\omega\to\pi^{+}\pi^{-}\pi^{0})$ and $Br(\pi^{0}\to\gamma\gamma)$ are the branching fractions of $\omega\to\pi^{+}\pi^{-}\pi^{0}$ and $\pi^{0}\to\gamma\gamma$ from the PDG pdg . We use $\epsilon$ to denote the detection efficiency determined by the TOY MC sample with detector simulation of helicity amplitude model. The numerical results of Born cross sections are listed in Table 3.