$\pi^-p\to\eta^{(\prime)}\, \pi^- p$ in the double-Regge region

The production of $\eta^{(\prime)}\pi$ pairs constitutes one of the golden channels to search for hybrid exotics, with explicit gluonic degrees of freedom. Understanding the dynamics and backgrounds associated to $\eta^{(\prime)}\pi$ production above the resonance region is required to impose additional constraints to the resonance extraction. We consider the reaction $\pi^-p\to \eta^{(\prime)} \pi^- \,p$ measured by COMPASS. We show that the data in $2.4<m_{\eta^{(\prime)}\pi}<3.0$ GeV can be described by amplitudes based on double-Regge exchanges. The angular distribution of the meson pairs, in particular in the $\eta' \pi$ channel, can be attributed to flavor singlet exchanges, suggesting the presence of a large gluon content that couples strongly to the produced mesons.


I. INTRODUCTION
Since the early days of the quark model, hadron spectroscopy has remained central to our understanding of QCD. High precision data on various reactions that have recently been collected from experiments at CERN, JLab, B-and charm factories have produced tantalizing evidence for the existence of exotic states that do not naturally fit within the quark model classification [1][2][3], e.g. pentaquark and tetraquark candidates [4][5][6]. The quantum numbers of some exotics are manifestly incompatible with a simple qq assignment. For example, states with J P C = 1 −+ have long been speculated to be hybrids, i.e. mesons where gluons play the role of constituents [7,8]. The paucity of data and the need for a thorough partial wave analysis to disentangle resonance from nonresonant background can be a challenging endeavor. The COMPASS collaboration extracted the ηπ and η π partial waves as a function of the invariant mass, m η ( ) π < 3.0 GeV from the measurement of diffractive pion dissociation on a nucleon target at 191 GeV [9]. These odd waves carry exotic quantum numbers, J P C = 1 −+ , 3 −+ ,... The key observations are that even waves are similar in both reactions, while the Pwave is significantly larger in η π. This reflects in a larger forward-backward asymmetry of the η π. Both channels present peaking structures in the P -waves at seemingly different masses. For a a long time, the two structures * lukasz.bibrzycki@up.krakow.pl † cesar.fernandez@nucleares.unam.mx were interpreted as two different states, the lighter one coupling mainly to ηπ and the heavier one to η π. However, the coupled-channel analysis in [10] showed that the data is consistent with the existence of a single exotic resonance. These conclusions have been confirmed by a recent independent analysis [11] and are supported by the latest lattice QCD computations [12]. At higher invariant masses, the reaction is expected to be dominated by cross-channel Regge exchanges, which is consistent with the cross section peaking in the forward and backward directions, with the peaks shrinking with increasing η ( ) π mass cf. Fig. 2 of Ref. [9]. Since a forward-backward asymmetry arises from the interference between even and odd waves, the larger exotic Pwave in η π is consistent with the observed larger asymmetry. This connection between resonances and Regge exchanges can be formalized via dispersion relations, e.g. in the form of finite energy sum rules [13][14][15]. Such relations can be used to constrain fits in the resonance region which, in combination with forthcoming high precision data from GlueX [16,17] and COMPASS [18], could lead to a more accurate determination of the exotic meson resonance parameters. A necessary step in this procedure is to fit the high mass region with analytical amplitudes that respect Regge asymptotic behavior. This is the main purpose of this work.
The paper is organized as follows. In Section II we describe the COMPASS partial waves, the procedure to compute the intensity distribution from them, and the main features of said distribution. Section III describes the double-Regge model used to fit the data. In Section IV we discuss the consequences of truncating the partial wave expansion in the analysis of the data and how it impacts the comparison to a given model and the extraction of the dynamics. In Section V we discuss what are the relevant contributions to the amplitudes, needed to reproduce the features of the angular and mass dependencies. Section VI describes our fitting strategy, fit results, and comparison to data. Section VII provides the connection between the COMPASS partial waves and the partial waves obtained from the double-Regge model. Finally, in Section VIII we summarize our results. The kinematical description of the η ( ) π reactions, statistical analysis, error propagation from the COMPASS partial waves to the intensity distribution, and other details and complementary information are left to the Appendices.

II. COMPASS INTENSITIES
In this section we describe the data on reactions analyzed by COMPASS [9]. The unpolarized cross sections for both reactions depend on five kinematical variables. These are, for example, the total center of mass energy squared s = (q + p 1 ) 2 , the invariant mass of the produced meson pair m 2 = m 2 η ( ) π = (k π + k η ) 2 , the square of the momentum transfer between the target and the recoil nucleon t p = (p 1 − p 2 ) 2 , and the spherical angle Ω determining the direction of the relative momentum between the two mesons in the rest frame of the pair. The COMPASS experiment operated with a fixed beam momentum of 191 GeV; in the analysis of [9] t p was integrated in the region t p ∈ [−1.0, −0.1] GeV 2 . Furthermore, since there was no measurement of the initial flux, the normalization of the event distribution is unknown. In the partial wave analysis of [9] the angular dependence of the event distribution, aka intensity function, in bins of m, was expanded in terms of angular functions, given by Ψ =+ LM (Ω) = √ 2Y M L (θ, 0) sin M φ and Ψ =− LM (Ω) = √ 2Y M L (θ, 0) cos M φ, which are the real spherical harmonics with referred to as the reflectivity. The angular variables Ω ≡ (θ, φ) determine the direction of the η ( ) in the Gottfried-Jackson (GJ) frame (see Appendix A for the axes orientation). The complex functions f LM (m) are obtained by fitting to the angular distributions for each energy bin m. In the strict sense they are not partial waves, as they do not depend on the initial and final nucleon helicities. However, if a single helicity amplitude happens to dominate the reaction, the f 's can approach genuine partial waves. In general however, one should think of the f 's as defining an effective parametrization of the data at the amplitude level. Nevertheless, in the following we refer to the f 's as partial waves, as customary.
In practice, partial wave extraction requires the sum in Eq. (2) to be truncated. In the COMPASS analysis of ηπ, seven partial waves were used, (L = 1, . . . , 6; M = 1) and (L = 2; M = 2), while for the η π channel it was six partial waves, namely (L = 1, . . . , 6; M = 1). All the waves describing the η ( ) π system have positive reflectivity = +. In the Regge asymptotic limit, reflectivity coincides with naturality of the exchange; at the nucleon vertex, the natural IP and f 2 are the dominant exchanges [19]. A single negative reflectivity wave was included in the fit, (L, M, ) = (0, 0, −), that includes possible reducible backgrounds. It was found to contribute at the 0.5% (1.1%) level to the total ηπ (η π) intensity, and will be neglected here.
The partial waves in Eq. (2) are written as where I LM (m) are the partial wave intensities and the phases δ LM (m) are determined with respect to the phase of the L = 2, M = 1 wave, i.e. δ 21 (m) ≡ 0. In our analysis, we use the intensities and phases provided in the corrigendum to Ref. [9]. The simplest way to compare the COMPASS results with a theoretical model would be to compare the partial waves. However, for reasons that will be discussed later in Section IV, we instead fit our amplitude model using an integral form of extended negative log-likelihood (ENLL) method [20][21][22] to the intensity I(m, Ω) reconstructed from the partial waves. There are two ways to reconstruct the I(m, Ω) from the COMPASS partial waves. One approach is to use the mean values of the intensity and phase at a given m, and use Eq. (2) to obtain I(m, Ω). We call this the mean value reconstruction (MVR). However, this method ignores the experimental uncertainties. The second method, which we refer to as MCR, uses Monte Carlo reconstruction. This is done by associating a probability distribution to the intensity and phase at each m independently. In doing so, instead of a single intensity value for each (m, φ, cos θ) point, we obtain a distribution. We can then compute the expected (mean) value of the intensity and its associated uncertainty at a given confidence level. The statistical errors are thus propagated from the partial waves to the intensity. The details on the MCR can be found in Appendix B. What remains unknown, however, are the uncertainties associated to the systematics of the COMPASS fit and the correlations among partial waves. As a consequence, the intensities reconstructed using MVR and MCR differ. In Figs. 1 and 2 we show the density plots of I(m, Ω) at three fixed m as well as the φ-integrated distributions In Fig. 3 we plot I θ (m, cos θ) for m above 2.3 GeV, for a total of seventeen mass bins in each channel. This    Fig. 1 for η π. We note that the η π backward peak is broader than the ηπ one and that the forward-backward asymmetry is more pronounced.
can be compared to the plot of the experimental data shown in Fig. 2   2. the intensity peaks in the forward cos θ ∼ 1 and backward cos θ ∼ −1 regions. In the forward region, most of the beam momentum is carried by the η ( ) , and in the backward region by the π. We call these clusters the "fast-η" and the "fast-π" regions, respectively; 3. the backward (fast-π) peak is larger than the forward (fast-η) peak, resulting in a forward-backward asymmetry. This effect is more pronounced in the case of the η π channel; 4. the backward peak is broader in η π than in the ηπ; 5. both the forward and backward peaks become narrower as the invariant mass m increases; 6. the MVR intensities at backward peak are larger than those of the MCR, and in the small |cos θ| region the intensity profile becomes smeared out in the MCR, so more structures are visible in the MVR in the region where intensities are low. Appendix C provides more insight on the differences between MVR and MCR.
These features are typical of diffractive processes, indicating the dominance of double-Regge exchanges in the energy region m 2.3 GeV. In the SU (3) flavor symmetric limit the π and the octet η 8 are degenerate, and so are the a 2 and f 2 Regge trajectories. Furthermore, if the SU (3) singlet exchanges (e.g. the IP ) are neglected, the forward and backward intensities are identical [23] for the production of the octet, and only even (nonexotic) waves contribute. Since η is dominated by the SU (3) singlet We note that the slope of F (m) for both ηπ and η π is steepper than for B(m).
we expect the asymmetry to be larger for the production of η π. The broadness of the peaks is related to the relative strength of the different double-Regge contributions to the amplitudes and will be addressed in Sections V and VI.
To quantify the forward-backward asymmetry we define with F (m) and B(m) being the forward and backward intensities, respectively, and A(m) the forward-backward asymmetry. Figure 4 shows F (m), B(m), and their sum T (m) for both MVR and MCR for the two channels. We find that the slope of F (m) is steeper than that of B(m). These intensities show clearly the difference between the MVR and the MCR, even though the total intensity in the MVR and MCR are similar.

III. DOUBLE-REGGE MODEL
We present here a double-Regge exchange model for the reactions in Eq. (1). Multi-Regge exchange formalism has been extensively studied theoretically in the past [15,[24][25][26][27][28]. An application of such formalism was presented in [29] for a similar reaction, two-pseudoscalar mesons production in K ± and π ± beam diffraction. More recently the double-Regge exchange was used to describe the central meson production in the high energy protonproton collisions [30][31][32]. We will adopt the same model and quote in this Section its main features.
The fast-η and fast-π regions correspond to the fast-η and fast-π double-Regge exchange amplitudes depicted in Fig. 5. The model assumes the dominance of leading Regge trajectories. Although it is known that daughter poles and cuts also contribute, for example to polarization observables [33][34][35][36], present data does not seem to be sensitive to subleading exchanges.
The top exchange is saturated by the a 2 trajectory for the fast-η amplitude, and by the f 2 or IP trajectory for the fast-π amplitude. The bottom exchange is either the f 2 or IP for both types of amplitude. It is common lore that, at COMPASS energies, the IP is the only relevant bottom exchange; however, this hypothesis is incompatible with data, as we will show in Section V.
Consequently, the total amplitude A Th (m, Ω) is the sum of six possible double-Regge amplitudes A Th (m, Ω) = c a2I P A a2I P + c a2f2 A a2f2 + c f2I P A f2I P + c f2f2 A f2f2 + c I P I P A I P I P + c I P f2 A I P f2 , (6) where the {c} are unknown and will be fitted to data.
Regge amplitudes are expressed in terms of Lorentz invariants. In addition to s, t p and m, as depicted in Fig. 5, for the fast-η and π amplitudes, the GJ angles are related to the following Lorentz invariants fast-η : There are only five independent variables. The fast-π invariant t π and s ηp can be expressed as linear combinations of the five fast-η variables. Appendix A summarizes the relevant kinematical relations.
The analytic structure is the same for all double-Regge amplitudes. The dependence in the momentum transferred (t η , t p ) for fast-η and (t π , t p ) for fast-π enters only via the trajectories (α 1 , α 2 ), where α 1 corresponds to the top exchange and α 2 to the bottom one. Hence, for fastη amplitudes α 1 ≡ α a2 (t η ) and for fast-π amplitudes α 1 ≡ α f2 (t π ) or α 1 ≡ α I P (t π ). The bottom trajectory is α 2 ≡ α f2 (t p ) or α 2 ≡ α I P (t p ) for both types depending on the bottom exchange.
Regge theory predicts the dependence in the invariant masses squared (s 1 , s 2 ) with (s 1 , s 2 ) = (s ηπ , s πp ) for the fast-η amplitudes and (s 1 , s 2 ) = (s ηπ , s ηp ) for the fast-π amplitudes. Since the nucleons play a spectator role given the large total energy, their spins can be ignored. For five spinless particles with an odd number of pseudoscalars, the generic amplitude for a double-Regge exchange is [15,29] The double-Regge limit corresponds to s, s 1 , s 2 → ∞ with κ −1 ≡ s/(α s 1 s 2 ) fixed, which is related to the cosine of the Toller angle [15]. The variation of α induces a smooth exponential dependence on the momentum trans-fer variables (t π , t η and t p ), that effectively absorbs possible form factor contributions. In particular, since the COMPASS measurement is performed at integrated bottom exchange momentum, the data are not sensitive to t p . Therefore our model Eq. (9) does not require any additional momentum transferred dependence. We have found that fitting simultaneously {c} and α does not lead to stable solutions, as the coefficients and the scale parameter are strongly correlated. Moreover, α should be of the order of the hadronic scale, O(GeV −2 ). We let α vary in exploratory fits and found α = 0.8 GeV −2 to be the optimal choice. The kinematical factor K is detailed in Appendix A.
The presence of two symmetric terms in the bracket of Eq. (9) is imposed from general considerations of the analytic structure of double-Regge amplitudes. The interested readers will find the technical details in Section 3.3 of Ref. [15].
The double-Regge amplitude of Eq. (9) has poles for positive integer values of the trajectories α 1 and α 2 , which are related to the spins of the physical particles in the t channel. Since only poles with even signature (−1) J = +1 can couple to ηπ and ππ, odd signature poles are removed by the signature factors The vertex function V (α 1 , α 2 , κ) is an analytic function of its arguments. Its most general form involves an infinite number of Reggeon-Reggeon-particle couplings and reduces to a polynomial in κ −1 for integer α 1 and α 2 [24]. In a dual model, all Reggeon-Reggeon-particle couplings are equal and the vertex simplifies to [15,37] where 1 F 1 is the confluent hypergeometric function of the first kind.
The six contributions in Eq. (6) are obtained from the generic double-Regge amplitude in Eq. (9) with the following substitutions Since the momentum transferred between the initial and final nucleon has been integrated over in the COMPASS analysis, we do not have access to the t p distribution. This distribution would allow us to discriminate between the bottom exchanges. Since the amplitude decreases exponentially with t p , we fix t p close to the COMPASS lowest limit, t p = −0.2 GeV 2 . Results are stable against small variation of this value.
Finally, we need to specify the Regge trajectories where we adopted the standard parametrization for the IP [38] and the f 2 [39] trajectories. Phenomenologically, the a 2 trajectory is very similar to that of ρ, which is referred as exchange degeneracy (EXD) [28,40]. Our model is thus entirely specified by the six real parameters {c}.
Each {c} in Eq. (6) is a product of two particle-Reggeonparticle couplings (top and bottom vertices) and one Reggeon-particle-Reggeon coupling (middle vertex). The particle-Reggeon-particle couplings could be extracted from quasi-two-body reactions [41,42], but the Reggeonparticle-Reggeon couplings are largely unknown. In principle, all couplings have residual dependence on t's that cannot be disentangled. This prevents us from imposing further relations among the η and η amplitude parameters.

IV. PARTIAL WAVE TRUNCATION BEYOND THE RESONANCE REGION
COMPASS extracted partial waves under the assumption that only a limited number of them contribute. This is justifiable in the resonance region, but as the invariant mass of the η ( ) π system increases so does the number of relevant waves. Since the overall intensity decreases in the high energy region, the significance of higher waves (L > 6) could not be established and, hence, they were neglected.
The Regge model developed in the previous Section is not based on a partial wave expansion and therefore implicitly includes all partial waves. One can thus study whether the approximation to truncate to L 6 waves is appropriate for our model. In Fig. 6 we show how the truncation affects the total intensity in the ηπ channel. We expand each amplitude into partial waves and then sum back only the ones considered in the COMPASS analysis. For example, at m ηπ = 2.64 GeV, the seven partial waves considered by COMPASS account for ∼ 80% of the intensity at the peak for the f 2 /f 2 exchange. At this m ηπ , only for the IP /IP and IP /f 2 amplitudes this truncation adequately reproduces the intensity (> 99% and > 97%, respectively). If we include the partial waves up to L = 10 with M = 1 and M = 2, the intensity of the amplitudes is almost completely recovered (> 99% for all amplitudes except for a 2 /f 2 , which is > 93%). In Fig. 7 we show the same plots for the η π channel. In this case, the main disagreement happens in the forward peak (a 2 /IP and a 2 /f 2 amplitudes), where only between 60% and 80% of the peak strength is accounted for by the COMPASS partial waves.
Thus, as mentioned earlier, in this energy region COMPASS waves should be considered as an effective parametrization of the data, rather than being directly compared with genuine partial waves from a model that contains an infinite number of waves. However, given that they have been extracted under the constraint of summing up to the total intensity, we can reconstruct I(m, Ω) from the partial waves using Eq. (2) with the two methods (MVR and MCR) explained in Section II and fit them with our model. In Section VII we will discuss how the model partial waves compare to the f 's extracted from the data by COMPASS.

V. THE MINIMAL SET OF AMPLITUDES
Our model described in Section III is completely determined by the six coefficients of the double-Regge exchange amplitudes. As a first approach, we fitted the intensity with all the six parameters unconstrained. However, those fits did not lead to a unique solution, sometimes having coefficients compatible with zero. The error estimation was unreliable. In order to make the fits to reach stable solutions, we have to restrict the parameters to at most four. Consequently, to establish which amplitudes must be neglected or included in the fits, in this Section we compare the angular and mass dependen-cies of the individual exchanges to the experimental ones from MVR shown in Figs. 1-4. Conclusions are identical for MCR.
In the SU (3) limit, the event distribution of ηπ becomes symmetric in cos θ. At the amplitude level, this is manifested via EXD, meaning that the parameters of the a 2 and f 2 Regge poles are equal, including the couplings, c a2I P −c f2I P and c a2f2 −c f2f2 . 1 Deviations from the EXD relation are manifested in the nonvanishing forward-backward asymmetry. The cos θ dependence is correlated to the t π or t η dependence arising from the top exchange trajectory. We thus expect the amplitudes with the same top exchange to have similar cos θ behavior. Both a 2 /f 2 and a 2 /IP amplitudes will be collectively denoted as top-a 2 amplitudes, and similarly for the topf 2 and top-IP amplitudes. By EXD, we expect that both top-a 2 and -f 2 matter.
As shown in Figs. 6 and 7, the top-a 2 and top-f 2 amplitudes produce a narrow forward and a narrow backward peak, respectively. The top-IP amplitudes produce a wider backward peak, which is due to the smaller slope of the IP trajectory. Given that the widths of the experimental backward peaks in ηπ and η π are similar to what is expected from the top-f 2 exchange, we conclude that the f 2 /f 2 and/or f 2 /IP amplitudes should account for most of the backward intensity. The residual contribution from the IP /IP and IP /f 2 amplitudes may be needed to further widen the peak. In particular, the top-IP contributions might be necessary for the η π channel.
We next investigate the mass dependence of the top-a 2 and top-f 2 amplitudes. In Fig. 4  Another important feature is the φ dependence. In Fig. 9 we compare the φ dependence in the forward region of I(m, Ω) and the top-a 2 amplitudes. We see that a single amplitude cannot reproduce the experimental distributions at all m. Therefore, we will include both a 2 /f 2 and a 2 /IP amplitudes in our fits.
In Fig. 10 we make the same comparison for the backward region. The f 2 /IP and IP /f 2 amplitudes do not peak at the correct position at any m. On the other hand, the f 2 /f 2 and the IP /IP match better the data. As explained above, the f 2 /f 2 amplitude is already favored by the observed B(m) slope.
In conclusion, the minimal set of amplitudes (MIN) common to both channels, consists of a 2 /IP , a 2 /f 2 , and for ηπ (upper row) and η π (lower row) compared to the same distributions for the f2/IP (blue), f2/f2 (red), IP /IP (green), and IP /f2 (orange) amplitudes. Each distribution is normalized to its peak value.
Additionally, as discussed earlier, we may extend this set in order to take into account the width of the backward peak. In particular, the η π peak is broader than predicted by the f 2 /f 2 amplitude. Including the IP /f 2 would help. However, it may disrupt the φ distribution as shown in Fig. 10. An option would be to include both IP /f 2 and f 2 /IP . However, as stated earlier in this Section, including more than four amplitudes makes the fits unstable. For this reason we do not include the IP /f 2 amplitude in any fits.
Therefore, we are left with two options to broaden the backward peak: either IP /IP or f 2 /IP . The IP /IP amplitude allows to broaden the backward peak without af-fecting much the φ dependence. It also would make the backward peak broader than the f 2 /IP exchange. The f 2 /IP exchange shifts the φ distribution to peak below π/2, but by interfering with f 2 /f 2 this shift may be reduced. Hence, we explore adding either the f 2 /IP or IP /IP amplitudes to the MIN set for fitting the intensities.
To summarize, the sets of amplitudes we explore are: (i) MIN, that includes the a 2 /IP , a 2 /f 2 and f 2 /f 2 amplitudes, i.e. parameter set {c a2I P , c a2f2 , c f2f2 }; (ii) MIN+f /IP , with parameter set {c a2I P , c a2f2 , c f2I P , c f2f2 }; (iii) MIN+IP /IP , with parameter set {c a2I P , c a2f2 , c f2f2 , c I P I P }.

A. Extended negative log-likelihood fit
The contribution of each amplitude in a given set (MIN, MIN+f /IP , and MIN+IP /IP ) is determined by fitting the MVR and the MCR distributions for each η ( ) π channel independently.
We first discuss how to fit the MVR. In each mass bin, the intensity I(m, Ω) depends on two angles Ω = (φ, cos θ). The continuous variables prevent us from using a standard χ 2 fit. Besides, we need to take into account the fact that the total intensity is a fixed quantity. Hence, we use an integral form of the extended negative loglikelihood function (ENLL) [20,21]: where the experimental intensity I Exp (m, Ω) is the fitting objective function computed using Eq. (2) with MVR, and the theoretical intensity I Th (m, Ω|{c}) is computed from Eq. (7). The experimental distributions are fitted simultaneously in 15 bins of {m}, in the range 2.38 < m < 2.98 GeV. Varying slightly this interval leaves the results unchanged. We minimize L using MINUIT [43] to obtain the {c} parameters weighting each theoretical amplitude. We note that the ENLL makes the total intensity of the model as close as possible to the total intensity of data. We remind the normalization of data is unknown, thus the {c} cannot be directly compared to normalized couplings. The overall sign of the amplitude is also undetermined, so we fix c a2I P to be positive. As said above, we expect c f2I P to be negative. The best fits found are reported in Table I. Local minima that do not follow the sign expectations were found, although with worse L than the reported best fit values. We note that the absolute value of the parameters is not a measure of the importance of any given amplitude contribution, because the A's in Eq. (12) have largely different magnitudes. In particular, bottom-IP are much larger than bottom-f 2 .
Fitting the MCR is more challenging. Each pseudodataset j is fitted using Eq. (14), obtaining an independent set of parameters {c} j . We estimate the expectation value of the parameters by averaging over N = 10 4 fits and the uncertainties from the appropriate quantiles. This number of pseudodatasets allows us to obtain the probability distribution of each parameter and the correlations with ∼ 1% statistical uncertainty (more details Appendix B). Table I gives the value of the ENLL for the three models fitted to the MVR and MCR for both channels, as well as the resulting fit parameters. The distribution of the ENLL for the MCR fits is shown in Fig. 11. We see that all the models have similar ENLL, with a nonsignificant preference for MIN+IP /IP fit, in particular for the η π channel. In Appendix D we analyze this difference more systematically and conclude that, statistically, there is indeed a preference for the MIN+IP /IP model for the η π channel.

B. Fit results
In Appendices C and E we compare MVR and MCR observables and fits, respectively, finding that MCR fits are more reliable. Here we summarize the results of the MCR fits and leave the MVR fit results for Appendix E.

ηπ MCR fits
The three models give consistent values for c a2I P and c a2f2 , providing almost identical descriptions of the forward peak. This can be appreciated in Fig. 12, where the experimental I θ (m ηπ , cos θ) MCR for the fast-η region is compared to the three models for all the fitted m ηπ bins. The three models agree very well. Figure 13 shows the same results for the fast-π region. Here the differences among models can be appreciated. As expected from Fig. 6, the MIN+IP /IP provides a wider peak, and, since the normalization is fixed in a ENLL fit, the maximum intensity at the peak is smaller than the MIN and MIN+f /IP results. The latter two fits are similar, with their uncertainty bands overlapping, except in the highest m ηπ bin. We note that for some energies the MIN+IP /IP provides a better description of the experimental distribution, while for others the MIN and MIN+f /IP fits look better.
Further insight can be obtained by examining the three-dimensional I(m ηπ , Ω) distributions. We define whereĪ and ∆I are the mean and dispersion of the experimental and theoretical distributions as obtained from MCR. 2 This quantifies point-by-point how similar the 2 We do not take into account the fact that both the theoretical MCR and the theoretical distributions are. Figure 14 shows these distributions for the ηπ channel in two mass bins. Figure 15 provides the φ distribution at the backward and forward peaks for three energies. Comparing the model and data distributions, one concludes that the former has more structure in the forward region. The experimental peak has two symmetric blobs in φ while the theory is rather asymmetric. As shown in Fig. 9, this is due to the asymmetry in φ of the top-a 2 amplitudes. We remind that the symmetry of the experimental φ distribution is exact and stems from Eq. (2). This symmetry is not imposed in the model, and is approximately reached by having both top-a 2 amplitudes interfering. All the three models peak at roughly the correct φ = π/2 and 3π/2. The situation is different for the backward peak. The MIN fit peaks slightly below (above) the experimental value of φ = 3π/2 (φ = π/2). Hence, we do not favor the MIN model, i.e. the f 2 /f 2 amplitude is not enough to reproduce the φ dependence of the fast-π region. From the bootstrap fits, we can study the parameter distributions and their correlations, summarized in Appendix F. The correlations confirm that the fast-π and fast-η amplitudes are essentially independent. The parameters are generally well determined and exhibit Gaussian behavior, except the IP /IP coefficient c I P I P that has a bimodal distribution. and experimental distributions are evaluated out of the same pseudodatasets, and therefore correlated. However, this still gives a qualitative description of the discrepancy between theory and experiment. Finally, although including a bottom-IP amplitude is necessary to describe the backward region, data do not show a clear preference for either MIN+f /IP or MIN+IP /IP . Since the two models point to different values for the f 2 /f 2 coupling, the latter cannot be determined unambiguously either. Currently, we do not have enough precision in the data to determine the contributions from the individual exchanges in the backward region.

η π MCR fits
The fit parameters for the three models are presented in the MCR columns of Table I. As for ηπ, the three models give consistent values for c a2I P and c a2f2 , providing almost identical descriptions of the forward peak. However, c a2f2 is compatible with zero at a 2σ level. This suggests larger level of EXD breaking in the η π channel. Fig-ures 16 and 17, compare the experimental I θ (m η π , cos θ) with the three models in the forward and backward regions, respectively. The models completely agree in the forward region, while the MIN+IP /IP provides a wider backward peak, in better agreement with the data.
The three-dimensional distributions for η π are shown in Figure 18 for the three models and MCR at m η π = 2.60 and 2.80 GeV. Figure 19 provides the φ distribution at the backward and forward peaks for three energies. As for ηπ, the MIN does not peak at the correct value of φ in the backward region. Results and conclusions are qualitatively similar to the ηπ channel. In particular, the preference for a MIN+IP /IP is clear. This points to a large affinity of η to gluons as discussed in the literature [44].  Figures 20 and 21 show the forward, backward, and total intensities. We show that MIN+f /IP for ηπ and MIN+IP /IP for η π reproduce all the intensities rather well.

Forward and backward intensities and asymmetry
We also note that the integrated forward intensity is systematically larger for the MCR. The opposite is true for the the backward one, which is systematically larger for the MVR. Consequently, the forward-backward asymmetry A(m), defined in Eq. (5c), is, in absolute value, larger for the MVR than for the MCR, as shown in Fig. 22. The existence of the asymmetry is a consequence of the odd (exotic) partial waves contribution. Taking into account the uncertainties in the partial waves makes the asymmetry less acute, but it is still sizeable and negative for both channels. The asymmetry is larger for the η π reaction, making this channel appropriate to search for hybrid candidates [7,8,44].

VII. CONSTRAINED PARTIAL WAVE ANALYSIS
As discussed earlier in Section IV, the f + LM (m) Exp amplitudes extracted by COMPASS are not exactly genuine partial waves. It is rather a parametrization that minimizes the ENLL estimator used to fit the actual event distributions. The ENLL fit makes a finite set of amplitudes reproduce the total intensity. Hence, any contribution from higher partial waves gets redistributed into the set included in the fit. Our model contains an infinite number of partial waves, which leads to a mismatch between the model partial waves and the COMPASS ones. The comparison can still be done if we project the model onto partial waves applying the same constrained procedure implemented by COMPASS. For simplicity, we consider I Th (m, Ω) of Eq. (7), as obtained from MVR. We follow the conventions in Eq. (2). For each energy bin m i , we extract the constrained partial waves (cPW) by minimizing the ENLL estimator (16) where I cPW is given by The L, M employed are the same truncated set as COMPASS.
The truncated set of partial waves suffers from the problem of discrete ambiguities, the so-called Barrelet zeros [45].
The cPW, uPW and COMPASS waves are shown in Figs. 23 and 24. As anticipated, the cPW agree with the COMPASS data very well, while the uPW can be quite different. Indeed, the truncation of uPW to the COMPASS set would reduce the integrated intensity to 86% for ηπ and 95% for η π. This is in agreement with the expectations from Figs effects were shown to be critical for ηπ. The dominant (2, 1) intensity in the ηπ channel is noteworthy. The unconstrained wave is very small compared to the COMPASS one, but the cPW matches the data, showing how the truncation makes low-lying partial waves to absorb the intensity of higher waves.

VIII. SUMMARY AND CONCLUSIONS
We studied the COMPASS data on the π − p → η ( ) π − p reactions for m > 2.38 GeV where the dynamics is expected to be dominated by Regge phenomenology. We considered a double-Regge model composed of up to six amplitudes that account for the possible top/bottom Regge exchanges. In particular, we included a 2 /IP and a 2 /f 2 to describe the fast-η (forward) region, and f 2 /IP , f 2 /f 2 , IP /IP , IP /f 2 for the fast-π (backward) region.
The COMPASS collaboration reported partial waves extracted from data under the assumption that only seven (six) partial waves contributed to the ηπ (η π) channel. This is justifiable in the resonance region, i.e. m η ( ) π 2 GeV. For higher energies, the number of relevant partial waves increases. Our Regge model is not based on a partial wave expansion and therefore implicitly includes all partial waves. Truncating to the set of waves used by COMPASS is not appropriate for this energy region, as in our model the discarded higher partial waves amount to a nonnegligible contribution to the intensities. Nevertheless, we reconstructed the total intensities from the COMPASS partial waves and fitted with our double-Regge model. We found that the ηπ intensity can be well described with four amplitudes, a 2 /IP , a 2 /f 2 , f 2 /f 2 , and either f 2 /IP or IP /IP . The inclusion of either bottom-IP amplitude is necessary to describe the forward region, but the data do not show a clear preference for either f 2 /IP or IP /IP amplitudes. For this reason, we could not disentangle the contributions from the individual ex- changes in the backward direction. In the η π channel, we found that the best model to reproduce the data consists of a 2 /IP , a 2 /f 2 , f 2 /f 2 , and IP /IP amplitudes. The IP /IP contribution is necessary to describe the data and points to a large gluon affinity of the η π system, potentially related to the existence of hybrid mesons. This is also consistent with the observed breakdown of exchange degeneracy between a 2 and f 2 in η π production.
The importance of the bottom-f 2 exchange, as shown in the slope of the integrated backward intensity, contradicts the common lore that, at COMPASS energies, the η ( ) π pairs are produced via IP exchange only, at least for this range of m η ( ) π .
A consequence of having an amplitude model that contains an infinite number of partial waves is that these cannot match the truncated waves from COMPASS. To bridge this apparent contradiction, we performed a constrained partial wave analysis of the model, using the same procedure as COMPASS. We found that these waves indeed agree well with the COMPASS ones. This proves the importance of studying the full amplitude rather than a truncated partial wave decomposition once the double-Regge regime is reached. variants as ξ and are the angles between the beam and the target, and between the beam and the recoil, respectively. They are given by where t η and s πp are defined in Eq. (8a). The energies and the nucleon angles depend only on s, m 2 ηπ and t p . The polar and azimuthal angles of the η must then depend on the remaining independent invariants. Indeed we obtain s πp = m 2 π + m 2 p + 2E 2 E π − 2|p 2 ||k| (sin sin θ cos φ + cos cos θ) . (A4b) The invariants t π and s ηp needed in fast-π amplitudes and defined in Eq. (8b) are related to the other Mandelstam variables by The kinematic function K that appears in Eq. (9) is given by

Appendix B: MCR and bootstrap
We describe how the MCR is performed, and consequently the bootstrap fit to it. The first step is to associate a probability distribution to each value of intensity with µ equal to the mean value and σ the uncertainty reported by COMPASS. For phase shifts, we use the von Mises distribution, i.e. the periodic equivalent to the Gaussian distribution, where I 0 (κ) is the modified Bessel function of order 0. The µ parameter is set to the mean value of each phase shift. The concentration parameterκ is the reciprocal measurement of the dispersion. If the uncertainty is small, the Gaussian distribution with σ equal to the experimental uncertainty is almost equal to the von Mises distribution withκ = 1/σ 2 , as shown in Fig. 26. For larger values of the uncertainty, the Gaussian and the von Mises distribution withκ = 1/σ 2 are quite different, and we decided to refit the concentration parameter to the Gaussian distribution. We believe this gives a better description of the actual experimental uncertainties. Nevertheless, we computed the MCR with the three options for the phase shifts distributions and we did not find relevant differences among the results. We are assuming that the uncertainties are statistical only, that data at different energy bins are statistically independent, and that the partial waves intensities and phases are uncorrelated. Although this is clearly not true, the correlation information is not available from the COMPASS analysis. This assumption likely leads to an overestimation of the error bands. There is an additional caveat: some of the uncertainties are large enough to make the intensities negative, which is unphysical. Hence, if a resampling provides a negative intensity for a given m, we set that particular intensity to zero. Alternatively, we consider a negative intensity in Eq. (2) instead of setting the value to zero. We found the difference to be negligible in the resulting MCR distributions and fits.
Once the distributions are chosen, we resample them N = 10 4 times to achieve a precision of ∼ 1% in the extracted distributions. For each resampled pseudodataset we can compute the value of the intensity I(m, Ω), and we can fit it to obtain the corresponding parameters {c} (bootstrap fit) [46][47][48]. Then, the expected value (mean) of each observable and the associated uncertainty (16% and 84% quantiles to obtain the 68% error bands) can be computed.

Appendix C: MVR vs. MCR observables
In Section II we mentioned that there were meaningful differences between the MVR and MCR. In particular, for I θ (m, cos θ) the forward peak is smaller for the MVR than for the MCR. This can be noticed in both Fig. 1  and 2. It can also be noticed in the integrated forward intensity of Figs. 20 and 21, where F (m) is systematically larger for the MCR. The opposite is true for the backward peak and B(m), which are systematically larger for the MVR. Consequently, the forward-backward asymmetry A(m) is, in absolute value, larger for the MVR than for the MCR, as shown in Fig. 22. The total intensity T (m) is very similar for both MCR and MVR, as displayed in Fig. 4. The difference between MVR and MCR is also apparent if we inspect the small |cos θ| region. The inclusion of the uncertainties and the calculation of expected values in the MCR leads to a smearing that flattens the experimental distribution. This effect is shown in Fig. 27. Fits to the MVR would try to match these structures that are mostly killed by uncertainties. However, removing the small |cos θ| region from the fits is not feasible, given that the total intensity is an important experimental constraint to ENLL fits. Consequently, we consider the physics extracted from the MCR fits more reliable. η π (right, red) resamples. The colored line is the mean of the distribution and the band represents the 68% confidence level. While for η π the distribution sits safely above ∆Lj = 0 (gray line), for ηπ this happens only for 96.5% of the pseudodatasets, corresponding to a 2σ preference. It is worth noticing that the agreement with the naive difference of the mean values of the ENLL in Table I MIN+f/I P MIN fits as they do not give an appropriate description of the I(m, Ω) distributions, see Sections VI B 1 and VI B 2. From the correlation matrices the substantial independence of fast-π and fast-η amplitudes is apparent for both models and channels.
For the ηπ channel, the parameters for the MIN+f /IP model are well determined and exhibit a Gaussian behavior. However, this does not happen for the MIN+IP /IP model, where the IP /IP amplitude parameter is not well determined and presents a bimodal distribution. The fit to the MVR using the MIN+IP /IP model presents a single and isolated absolute minimum, so the appearance of a two peak structure in the fit to the MCR is entirely due to the inclusion of the uncertainties. The c I P I P parameter cannot be well determined. The c a2I P and c a2f2 distribu-tions for mostly overlap in both models. For c f2f2 , the distributions barely overlap, as a consequence of the differences between the f 2 /IP and IP /IP amplitudes.
For the η π channel, all the parameters for both models are well determined and exhibit a nice Gaussian behavior. The c a2f2 distribution is compatible with zero at a 2σ level for both models, indicating that it is possible that the associated amplitude vanishes. This would mean a large violation of the EXD between the a 2 and f 2 Regge poles. Given that both the statistical analysis in Appendix D and the comparison to the data in Section VI B 2 favor the MIN+IP /IP , we find that the contribution of all four amplitudes (a 2 /IP , a 2 /f 2 , f 2 /f 2 , and IP /IP ) to the η π process can be well established with reasonable uncertainties.