New analysis of $\eta\pi$ tensor resonances measured at the COMPASS experiment

We present a new amplitude analysis of the $\eta\pi$ $D$-wave in $\pi^- p\to \eta\pi^- p$ measured by COMPASS. Employing an analytical model based on the principles of the relativistic $S$-matrix, we find two resonances that can be identified with the $a_2(1320)$ and the excited $a_2^\prime(1700)$, and perform a comprehensive analysis of their pole positions. For the mass and width of the $a_2$ we find $M=(1307 \pm 1 \pm 6)$~MeV and $\Gamma=(112 \pm 1 \pm 8)$~MeV, and for the excited state $a_2^\prime$ we obtain $M=(1720 \pm 10 \pm 60)$~MeV and $\Gamma=(280\pm 10 \pm 70)$~MeV, respectively.

i Supported by the DFG cluster of excellence 'Origin and Structure of the Universe' (www.universecluster.de) (Germany) j Supported by the Laboratoire d'excellence P2IO (France) k Also at Chubu University, Kasugai, Aichi 487-8501, Japan ab l Also at Dept. of Physics, National Central University, 300 Jhongda Road, Jhongli 32001, Taiwan m Present address: LP-Research Inc., Tokyo, Japan n Also at KEK, 1-1 Oho, Tsukuba, Ibaraki 305-0801, Japan o Also at Moscow Institute of Physics and Technology, Moscow Region, 141700, Russia p Present address: RWTH Aachen University, III. Physikalisches Institut, 52056 Aachen, Germany q Supported by CERN-RFBR Grant 12-02-91500 r Present address: LP-Research Inc., Tokyo, Japan s Also at Dept. of Physics, National Kaohsiung Normal University, Kaohsiung County 824, Taiwan t Also at University of Eastern Piedmont, 15100 Alessandria, Italy u Present address: Uppsala University, Box 516, 75120 Uppsala

Introduction
The spectrum of hadrons contains a number of poorly determined or missing resonances, the better knowledge of which is of key importance for improving our understanding of Quantum Chromodynamics (QCD). Active research programs in this direction are being pursued at various experimental facilities, including the COMPASS and LHCb experiments at CERN [1,2,3,4], CLAS/CLAS12 and GlueX at JLab [5,6,7], BESIII at BEPCII [8], BaBar, and Belle [9]. To connect the experimental observables with the QCD predictions, an amplitude analysis is required. Fundamental principles of S-matrix theory, such as unitarity and analyticity (which originate from probability conservation and causality), should be applied in order to construct reliable reaction models. When resonances dominate the spectrum, which is the case studied here, unitarity is especially important since it constrains resonance widths and allows us to determine the location of resonance poles in the complex energy plane of the multivalued partial wave amplitudes.
In 2014, COMPASS published high-statistics partial wave analyses of the π − p → η ( ) π − p reaction, at p beam = 191 GeV [2]. The odd angular-momentum waves have exotic quantum numbers and exhibit structures that may be compatible with a hybrid meson [10]. The even waves show strong signals of non-exotic resonances. In particular, the D-wave of ηπ, with I G (J PC ) = 1 − (2 ++ ), is dominated by the peak of the a 2 (1320) and its Breit-Wigner parameters were extracted and presented in Ref. [2]. The D-wave also exhibits a hint of the first radial excitation, the a 2 (1700) [11].
In this letter we present a new analysis of the D-wave based on an analytical model constrained by unitarity, which extends beyond a simple Breit-Wigner parameterization. Our model builds on a more general framework for a systematic analysis of peripheral meson production, which is currently under development [12,13,14]. Using the 2014 COMPASS measurement as input, the model is fitted to the results of the mass-independent analysis that was performed in 40 MeV wide bins of the ηπ mass. The a 2 and a 2 resonance parameters are extracted in the single-channel approximation and the coupled-channel effects are estimated by including the ρπ final state. We determine the statistical uncertainties by means of the bootstrap method [15,16,17,18,19], and assess the systematic uncertainties in the pole positions by varying model-dependent parameters in the reaction amplitude.
To the best of our knowledge, this is the first precision determination of pole parameters of these resonances that includes the recent, most precise, COMPASS data.

Reaction Model
We consider the peripheral production process π p → ηπ p ( Fig. 1(a)), which is dominated by Pomeron (P) exchange. High-energy diffractive production allows us to assume factorization of the "top" vertex, so that the πP → ηπ amplitude resembles an ordinary helicity amplitude [20]. It is a function of s and t 1 , the ηπ invariant mass squared and the invariant momentum transfer squared between the incoming pion and the η, respectively. It also depends on t, the momentum transfer between the nucleon target and recoil. In the Gottfried-Jackson (GJ) frame [21], the Pomeron helicity in πP → ηπ equals the ηπ total angular momentum projection M, and the helicity amplitudes a M (s,t,t 1 ) can be expanded in partial waves a JM (s,t) with total angular momentum J = L. The allowed quantum numbers of the ηπ partial waves are J P = 1 − , 2 + , 3 − , . . .. The Pomeron exchange has natural parity. Parity relates the amplitudes with opposite spin projections a JM = −a J−M [22]. That is, the M = 0 amplitude is forbidden and the two M = ±1 amplitudes are given, up to a sign, by a single scalar function.
The assumption about the Pomeron dominance can be quantified by the magnitude of unnatural partial waves. In the analysis of ref. [2], the magnitude of the L = M = 0 wave was estimated to be < 1%, and it also absorbs other possible reducible backgrounds. The patterns of azimuthal dependence in the central production of mesons [23,24,25,26,27] indicate that at low momentum transfer, t ∼ 0, the Pomeron behaves as a vector [28,29], which is in agreement with the strong dominance of the |M| = 1 component in the COMPASS data. 1 We are unable to further address the nature of the exchange from the data of ref. [2] since they are integrated over the momentum transfer t 2 . However, from this single energy analysis we cannot be sure the exchange is purely Pomeron. Analyses such as Ref. [32] suggest there could be an f exchange, but for our analysis the natural parity exchanges will be similar, so we consider an effective Pomeron which may be a mixture of pure Pomeron and f . We note here that COMPASS has published data in the 3π channel, which are binned both in 3π invariant mass and momentum transfer t [3].
The COMPASS mass-independent analysis [2] is restricted to partial waves with L = 1 to 6 and |M| = 1 (except for the L = 2 where also the |M| = 2 wave is taken into account). The lowest mass exchanges in the crossed channels of πP → ηπ correspond to the a (in the t 1 channel) and the f (in the u 1 channel) trajectories, thus higher partial waves are not expected to be significant in the ηπ mass region of interest, √ s < 2 GeV. However, the systematic uncertainties associated with an analysis based on a truncated set of partial waves is hard to estimate.
To compare with the partial wave intensities measured in Ref. [2], which are integrated over t from t min = −1.0 GeV 2 to t max = −0.1 GeV 2 , we use an effective value for the momentum transfer t eff = −0.1 GeV 2 and a JM (s) ≡ a JM (s,t eff ). The effect of a possible t eff dependence is taken into account in the estimate of the systematic uncertainties. The natural parity exchange partial wave amplitudes a JM (s) can be identified with the amplitudes A ε=1 LM (s) as defined in Eq. (1) of Ref. [2], where ε = +1 is the reflectivity eigenvalue that selects the natural parity exchange.
In the following we consider the single, J = 2, |M| = 1 natural parity partial wave, which we denote by a(s), and fit its modulus squared to the measured (acceptance corrected) number of events [2]: Here, I(s) is the intensity distribution of the D wave, p = λ 1/2 (s, m 2 η , m 2 π )/(2 √ s) the ηπ breakup momentum, and q = λ 1/2 (s, m 2 π ,t eff )/(2 √ s), which will be used later, is the π beam momentum in the ηπ rest frame with λ (x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz being the Källén triangle function. Since the physical normalization of the cross section is not determined in Ref. [2], the constant N on the right hand side of Eq. (1) is a free parameter.
In principle, one should consider the coupled-channel problem involving all the kinematically allowed intermediate states (see Fig. 1(b)). The PDG reports the important final states for the 2 ++ system are the 3π (ρπ, f 2 π) and ηπ systems [11]. Far from thresholds, a narrow peak in the data is generated by a pole in the closest unphysical sheet, regardless of the number of open channels. The residues (related to the branching ratios) depend on the individual couplings of each channel to the resonance, and therefore their extraction requires the inclusion of all the relevant channels. However, the pole position is expected to be essentially insensitive to the inclusion of multiple channels. This is easily understood in the Breit-Wigner approximation, where the total width extracted for a given state is independent of the branchings to individual channels. Thus, when investigating the pole position, we restrict the analysis to the elastic approximation, where only ηπ can appear in the intermediate state. We will elaborate on the effects of introducing the ρπ channel, which is known to be the dominant one of the decay of a 2 (1320) [11], as part of the systematic checks.
In the resonance region, unitarity gives constraints for both the ηπ interaction and production. Denoting the ηπ → ηπ scattering D-wave by f (s), unitarity and analyticity determine the imaginary part of both amplitudes above the ηπ threshold s th = (m η + m π ) 2 : From the analysis of kinematical singularities [33,34,35] it follows that the amplitude a(s) appearing in Eq.
(1) has kinematical singularities proportional to K(s) = p 2 q, and f (s) has singularities proportional to p 4 . The reduced partial waves in Eqs. (2) and (3) are free from kinematical singularities, and defined by e.g.â(s) = a(s)/K(s),f (s) = f (s)/p 4 , with ρ(s) = 2p 5 / √ s being the two-body phase space factor that absorbs the barrier factors of the D-wave. Note that Eq. (2) is the elastic approximation of Fig. 1(b).
We writef in the standard N-over-D form,f (s) = N(s)/D(s), with N(s) absorbing singularities from exchange interactions, i.e. "forces" acting between ηπ also known as left hand cuts, and D(s) containing the right hand cuts that are associated with direct channel thresholds. Unitarity leads to a relation between D and N, Im D(s) = −ρ(s)N(s), with the general once-subtracted integral solution Here, the function D 0 (s) is real for s > s th and can be parameterized as Note that the subtraction constant has been absorbed into c 0 of D 0 (s). The rational function in Eq. (5) is a sum over two so-called Castillejo-Dalitz-Dyson (CDD) poles [36], with the first pole located at s = ∞ (CDD ∞ ) and the second one at s = c 3 . The CDD poles produce real zeros of the amplitudef and they also lead to poles off on the complex plane (second sheet). Since these poles are introduced via parameters like c 1 , c 2 , rather than being generated through N (cf. Eq. (4)), they are commonly attributed to genuine QCD states, i.e. states that do not originate from effective, long-range interactions such as pion exchange [37]. In order to fix the arbitrary normalization of N(s) and D(s), we set c 0 = (1.23) 2 since it is expected to be numerically close to the a 2 mass squared expressed in GeV. One also expects c 1 to be approximately equal to the slope of the leading Regge trajectory [38]. The quark model [39] and lattice QCD [40] predict two states in the energy region of interest, so we use only two CDD poles. It follows from Eq. (4) that the singularities of N(s) (which originate from the finite range of the interaction) will also appear on the second sheet in D(s), together with the resonance poles generated by the CDD terms. We use a simple model for N(s), where the left hand cut is approximated by a higher order pole, Here, g and s R effectively parameterize the strength and inverse range of the exchange forces in the Dwave, whereas the power n = 7 makes the integral in Eq. (4) account for the finite range of interactions with the appropriate powers to regulate the threshold singularities, and additional powers as the model for the left hand singularities. The parameterization of N(s) removes the kinematical 1/s singularity in ρ(s). Therefore, dynamical singularities on the second sheet are either associated with the particles represented by the CDD poles, or the exchange forces parameterized by the higher order pole in N(s).
The general parameterization forâ(s), which is constrained by unitarity in Eq. (2), is obtained following similar arguments and is given by a ratio of two functionŝ where D(s) is given by Eq. (4) and brings in the effects of ηπ final state interactions, while n(s) describes the exchange interactions in the production process πP → ηπ and contains the associated left hand singularities. In both the production process and the elastic scattering no important contributions from light meson exchanges are expected since the lightest resonances in the t 1 and u 1 channels are the a 2 and f 2 mesons, respectively. Therefore, the numerator function in Eq. (7) is expected to be a smooth function of s in the complex plane near the physical region, with one exception: the CDD pole at s = c 3 produces a zero inâ(s). Since a zero in the elastic scattering amplitude does not in general imply a zero in the production amplitude, we write n(s) as where the function to the right of the pole is expected to be analytical in s near the physical region. We parameterize it using the Chebyshev polynomials T j , with ω(s) = s/(s + Λ) approximating the left hand singularities in the production process, πP → ηπ. The real coefficients a j are determined from the fit to the data. In the analysis, we fix Λ = 1 GeV 2 . We choose an expansion in Chebyshev polynomials as opposed to a simple power series in ω to reduce the correlations between the a j parameters. Since we examine the partial wave intensities integrated over the momentum transfer t, we assume that the expansion coefficients are independent of t. The only t-dependence comes from the residual kinematical dependence on the breakup momentum q.
A comment on the relation between the N-over-D method and the K-matrix parameterization is worth noting. If one assumes that there are no left hand singularities, i.e. let N(s) be a constant, then Eq. (4) is identical to that of the standard K-matrix formalism [41]. Hence we can relate both approaches through K −1 (s) = D 0 (s). It is also worth noting that the parameterization in Eq. (5) automatically satisfies causality, i.e. there are no poles on the physical energy sheet.    The absolute value is taken as there is a phase ambiguity because we fit only the intensity ∼ |a(s)| 2 . Note that each curve is an independent fit for a specific number of terms n p . The curves for n p = 4, 5, and 6 all coincide in the resonance region, as shown in the inset. Table 1: Best fit denominator and production parameters for the fit with two CDD poles, s R = 1.5 GeV 2 , N = 10 6 , c 0 = (1.23) 2 , and the number of expansion parameters n p = 6, leading to χ 2 /d.o.f. = 1.91. Denominator uncertainties are determined from a bootstrap analysis using 10 5 random fits. We report no uncertainties on the production parameters as they are highly correlated.

Methodology
We fit our model to the intensity distribution for π − p → ηπ − p in the D-wave (56 data points) [2], as defined in Eq. (1), by minimizing χ 2 . We fix the overall scale, N = 10 6 (see Eq. (1)), and fit the coefficients a j (see Eq. (8)), which are then expected to be O(1), and also the parameters in the D(s) function. In the first step we obtain the best fit for a given total number of parameters, and in the second step we estimate the statistical uncertainties using the bootstrap technique [15,16,17,18,19]. That is to say, we generate 10 5 pseudodata sets, each data point being resampled according to a Gaussian distribution having as mean and standard deviation the original value and error, and we repeat the fit for each set. In this way, we obtain 10 5 different values for the fit parameters, and we take the means and standard deviations as expected values and statistical uncertainties, respectively.
In order to assess the systematic uncertainties we study the dependence of the pole parameters on variations of the model, specifically we change i) the number of CDD poles from 1 to 3, ii) the total number of terms n p in the expansion of the numerator function n(s) in Eq. (8), iii) the value of s R in the left hand cut model, iv) the value of t eff of the total momentum transfered, and v) the addition of the ρπ channel to study coupled-channel effects.
The parameters corresponding to the best fit with two CDD poles are given in Table 1. The addition of another CDD pole does not improve the fit, as the limited precision in the data is incapable of indicating any further resonances. Specifically the residue of the additional pole turns out to be compatible with zero, leaving the other fit parameters unchanged. We associate no systematic uncertainty to that.
As discussed earlier, an acceptable numerator function n(s) should be "smooth" in the resonance region, i.e. without significant peaks or dips on the scale of the resonance widths. The parameters c i and g of the denominator function are related to resonance parameters, while s R controls the distant second sheet singularities due to exchange forces. The expansion in n(s), shown in Fig. 3 for s R = 1.5 GeV 2 and two CDD poles, has a singularity occurring at s = −1.0 GeV 2 because of the definition of ω(s) and our choice of Λ 3 . For variations in n(s) between n p = 3 and n p = 7, we find the pole positions are relatively stable, which we discuss later in our systematic estimates.
The dependence on t eff is expected to affect mostly the overall normalization. Indeed, the variation from t eff = −1.0 GeV 2 to −0.1 GeV 2 gives less than 2% difference for the a 2 (1700) parameters, and < 1 for the a 2 (1320), and can be neglected compared to the other uncertainties.

Results
This analysis allows us to extract the ηπ → ηπ elastic amplitude in the D-wave. By construction, the amplitude has a zero at s = c 3 . Figure 4 shows the real and imaginary parts off (s), with the 3σ error bands estimated by the bootstrap analysis. Resonance poles are extracted by analytically continuing the denominator of the ηπ elastic amplitude to the second Riemann sheet (II) across the unitarity cut using D II (s) = D(s) + 2iρ(s)N(s). By construction, no first-sheet poles are present. We find three second-sheet poles in the energy range of (m π + m η ) ≤ √ s ≤ 3 GeV, two of which can be identified as resonances, as shown in Fig. 5 for n p = 6 and s R = {1.0, 1.5, 2.0, 2.5} GeV 2 .
The mass and width are defined as m = Re √ s p and Γ = −2 Im √ s p , respectively, where s p is the pole position in the s plane. Two of the poles found can be identified as the a 2 (1320) and a 2 (1700) resonances in the PDG [11]. The lighter of the two corresponds to the a 2 (1320). For s R = 1.5 GeV 2 , the pole has mass and width m = (1307 ± 1) MeV and Γ = (112 ± 1) MeV, respectively. The nominal value is the best fit pole position, and the uncertainty is the statistical deviation determined from the bootstrap. Values of s R between 1.0 and 2.5 GeV 2 lead to pole deviations of at most ∆ m = 2 MeV and ∆ Γ = 3 MeV. The heavier pole corresponds to the excited a 2 (1700). For s R = 1.5 GeV 2 , the resonance has mass and width The a 2 (1320) and a 2 (1700) poles (see Fig. 5) are found to be stable under variations of s R , which modulates the left hand cut. As expected, there is a third pole that depends strongly on s R and reflects the singularity in N(s) modeled as a pole. Its mass ranges from 1.4 to 3.3 GeV, and its width varies between 1.3 and 1.8 GeV as s R changes from 1 GeV 2 to 2.5 GeV 2 .
In the limit g → 0, this pole moves to −s R as expected, while the other two migrate to the real axis above threshold [42].
Changing the number of expansion terms between n p = 3 and n p = 7 does not in any significant way affect the a 2 (1320) or a 2 (1700) pole positions. The maximal deviations are ∆ m(a 2 ) = 5 MeV, ∆ Γ(a 2 ) = 7 MeV and ∆ m(a 2 ) = 40 MeV, ∆ Γ(a 2 ) = 30 MeV between three and seven terms in the n(s) expansion.
In order to demonstrate that coupled-channel effects do not influence the pole positions, we consider an extension of the model to include a second channel also measured by COMPASS, ρπ [3], and simultaneously fit the ηπ [2] and the ρπ [3] final states. The branching ratio of the a 2 (1320) is saturated at the level of ∼85% by the ηπ and 3π channels [11], with the ρπ S-wave having the dominant contribution. For simplicity we consider the ρ to be a stable particle with mass 775 MeV, the finite width of the ρ being relevant only for √ s < 1 GeV. The amplitude is thenâ j (s) = ∑ k [D(s)] −1 jk (s) n k (s). The denominator is now a 2 × 2 matrix, whose diagonal elements are of the form given by Eq. (4), with the appropriate phase space for each channel. The off-diagonal term is parameterized as a single real constant. The production elements n k (s) are as in Eq. (8), with independent coefficients for each channel. We also performed a K matrix coupled-channel fit and obtained very similar results that are shown in Figure 6. The coupledchannel effects produce a competition between the parameters in the numerators to fit the bump at 1.6 GeV in ηπ and the dip at 1.8 GeV in ρπ at the same time. The ρπ data prefers not to have any excited a 2 (1700), which conversely is evident in the ηπ data. Therefore, the uncertainty in the a 2 (1700) pole position increases, as it is practically unconstrained by the ρπ data. Note, however, that in Ref. [3] the dip at √ s ∼ 1.8 GeV in the ρπ data is t-dependent, while we use the t-integrated intensity, so it may be expected that the effects of the a 2 are suppressed in our combined fit.
We find the following deviations in the pole positions relative to the single-channel fit: ∆m(a 2 ) = 2 MeV, ∆Γ(a 2 ) = 3 MeV, ∆m(a 2 ) = 20 MeV and ∆Γ(a 2 ) = 10 MeV. These deviations are rather small and we quote them within our systematic uncertainties.

Summary and Outlook
We describe the 2 ++ wave of π p → ηπ p reaction in a single-channel analysis emphasizing unitarity and analyticity of the amplitude. These fundamental S-matrix principles significantly constrain the possible form of the amplitude making the analysis more stable than standard ones that use sums of Breit-Wigner resonances with phenomenological background terms.
The robustness of the model allows us to reliably reproduce the data, and to extract pole positions by analytical continuation to the complex s-plane. We use the single-energy partial waves in Ref. [2] to extract the pole positions. We find two poles that can be identified as the a 2 (1320) and the a 2 (1700) resonances, with pole parameters m(a 2 ) = (1307 ± 1 ± 6) MeV, m(a 2 ) = (1720 ± 10 ± 60) MeV, Γ(a 2 ) = (112 ± 1 ± 8) MeV, Γ(a 2 ) = (280 ± 10 ± 70) MeV, where the first uncertainty is statistical (from the bootstrap analysis) and the second one systematic. The systematic uncertainty is obtained adding in quadrature the different systematic effects, i.e. the dependence on the number of terms in the expansion of the numerator function n(s), on s R , on t eff (negligible), and on the coupled-channel effects. The a 2 results are consistent with the previous a 2 (1320) results found in Ref. [2]. We note that a new mass-dependent COMPASS analysis of the 3π final state using Breit-Wigner forms in 14 waves is in progress.
The third pole found tends to −s R in the limit of vanishing coupling, indicating that this pole arises from the treatment of the exchange forces, and not from the CDD poles that account for the resonances.
In the future this analysis will be extended to also include the η π channel [43], where a large exotic P-wave is observed [2].
Additional material is available online through an interactive website [44,45].