Azimuthal asymmetries of back-to-back $\pi^\pm-(\pi^0,\eta,\pi^\pm)$ pairs in $e^+e^-$ annihilation

This work reports the first observation of azimuthal asymmetries around the thrust axis in $e^+e^-$ annihilation of pairs of back-to-back charged pions in one hemisphere, and $\pi^0$ and $\eta$ mesons in the opposite hemisphere. These results are complemented by a new analysis of pairs of back-to-back charged pions. The $\pi^0$ and $\eta$ asymmetries rise with the relative momentum $z$ of the detected hadrons as well as with the transverse momentum with respect to the thrust axis. These asymmetries are sensitive to the Collins fragmentation function $H_1^{\perp}$ and provide complementary information to previous measurements with charged pions and kaons in the final state. In particular, the $\eta$ final states will provide additional information on the flavor structure of $H_1^{\perp}$. This is the first measurement of the explicit transverse-momentum dependence of the Collins fragmentation function from Belle data. It uses a dataset of 980.4~fb$^{-1}$ collected by the Belle experiment at or near a center-of-mass energy of 10.58 GeV.

This work reports the first observation of azimuthal asymmetries around the thrust axis in e + e − annihilation of pairs of back-to-back charged pions in one hemisphere, and π 0 and η mesons in the opposite hemisphere. These results are complemented by a new analysis of pairs of back-to-back charged pions. The π 0 and η asymmetries rise with the relative momentum z of the detected hadrons as well as with the transverse momentum with respect to the thrust axis. These asymmetries are sensitive to the Collins fragmentation function H ⊥ 1 and provide complementary information to previous measurements with charged pions and kaons in the final state. In particular, the η final states will provide additional information on the flavor structure of H ⊥ 1 . This is the first measurement of the explicit transverse-momentum dependence of the Collins fragmentation function from Belle data. It uses a dataset of 980.4 fb −1 collected by the Belle experiment at or near a centerof-mass energy of 10.58 GeV.

I. Introduction
A description of the three-dimensional partonic structure of the nucleon is an essential test for our understanding of quantum chromodynamics (QCD). Successful tools for the study of the nucleon have been semi-inclusive hard reactions, particularly the use of leptonic probes such as electrons and muons. At high enough momentum transfers, QCD factorization theorems can be applied, and the process can be described using a convolution over parton distribution functions (PDFs), fragmentation func-tions (FFs), and the matrix element describing the elementary hard scattering of the probe off the parton inside the nucleon. PDFs [1] can be interpreted as the leading coefficients of the wave function of the nucleon on the light-cone in a Q 2 expansion, where Q 2 is the squared 4-momentum transfer, and have an probabilistic interpretation in the parton model as the probability of finding a parton q in the nucleon carrying a momentum fraction x of the parent nucleon. So-called unintegrated PDFs also carry a dependence on the transverse momentum of the struck quark. Fragmentation functions [2], on the other hand, describe the hadronization of a quark into a final-state hadrons containing at least one detected hadron. Fragmentation functions depend on the dimensionless variable z, which, in a partonic picture, can be interpreted as the momentum fraction of the struck quark carried by the detected hadron. In addition, unintegrated FFs depend on the transverse momentum P h⊥ of the hadron with respect to the initial quark direction. Since FFs encode the dependence of the properties of the detected hadron with the quantum numbers of the struck quark, knowledge of them is essential for the extraction of information on the partonic structure of the nucleon from semi-inclusive hard scattering experiments. This is in particular true for the transverse spin structure of the nucleon. The large single transverse spin asymmetries of π 0 and η mesons observed in pp collisions were at odds with the expectation that they would vanish due to the suppression of spin-flip amplitudes in the hard scattering [3]. However, Collins showed [4] that spin-flip amplitudes for soft components of the cross section, the PDFs and FFs, are not necessarily suppressed. In the collinear picture, in which the dependence of the PDFs and FFs on intrinsic transverse momenta is integrated over, the PDF that corresponds to the spin-flip amplitude is the so-called transversity PDF h 1 [5][6][7][8]. This can be interpreted as the probability of finding a transversely polarized quark in a transversely polarized nucleon with its polarization direction along the polarization of the parent nucleon and is one of the three leading-twist PDFs needed to describe the nucleon in a collinear picture. It is a chiral-odd function, and since chiral-odd amplitudes are strongly suppressed in perturbative QCD [3], h 1 has to be coupled to another chiral-odd function to construct a chiral-even observable such as a cross section. Experimentally, the most relevant channels to access transversity are transverse single spin asymmetries in semi-inclusive deep-inelastic scattering (SIDIS) or pp scattering. Here transversity couples, for instance, to the transverse polarization dependent chiral-odd Collins FF H ⊥ 1 [4] or the di-hadron interference FF H 1 [9,10]. Since both the transversity PDF as well as the transverse polarization dependent FFs are a priori unknown, an independent measurement of the FF is needed. Such a measurement can be performed in e + e − annihilation, where a back-to-back qq pair is created and hadronizes. The azimuthal dependence of the cross section of backto-back production of hadrons can be described by the product of the quark and anti-quark H ⊥ 1 together with the polarization averaged FFs. This allows access to the Collins FF without the complication of other, potentially unknown, functions that cannot be calculated in perturbative QCD. A disadvantage of e + e − annihilation at the energies relevant for FF measurements is the small sensitivity to gluon fragmentation as well as to the flavor of the fragmenting quark. This is because the production probability of all light quarks solely depend on e 2 q , where e q is the electric charge of the quark, and it is assumed that e + e − annihilation into virtual photons dominates, as in selected Belle data.
The first unambiguous observation of the Collins effect came from SIDIS off transversely polarized protons [11]. The behavior of the observed π + and π − asymmetries indicated that the Collins FF had opposite signs for favored versus disfavored fragmentation [cf. Eq. (11)], motivated also by the Schäfer-Teryaev sum rule for the Collins FFs [12]. These results spurred a wide range of both theoretical and experimental activities. The first measurement sensitive to the Collins FF for charged pions in e + e − annihilation was performed at Belle [13,14]. It was subsequently used, together with SIDIS data, for the first extraction of transversity in a global fit [15]. The Belle results were confirmed by BaBar [16]. Later, BaBar also reported the transverse momentum dependence as well as the observation of a significant signal for asymmetries involving kaons [17]. At lower energies, Collins asymmetries in e + e − annihilation have been measured by the BESIII collaboration [18]. The Q 2 dependence of the Collins function might provide interesting insight into the non-trivial evolution of transverse momentum dependent functions (cf. Ref. [2] and references therein).
Here, we report the first measurement of azimuthal asymmetries in back-to-back production of hadron pairs, where one hadron is a charged pion and the other hadron a π 0 or an η. We report the fractional-energy and the transverse-momentum dependence of these asymmetries as well as of asymmetries for charged pions. These results provide additional constraints on the Collins function in global fits. The final states including η mesons will provide sensitivity to the fragmentation of strange quarks and are also of interest since there are hints that the transverse spin asymmetries of π 0 and η mesons in pp collisions are different [19,20]. This paper is structured as follows: In Sec. II the observables are introduced, Sec. III briefly describes the Belle detector. Section IV details the analysis steps, Sec. V reports the result, and Sec. VI provides the summary and conclusion. Data tables are provided in two Appendices. In the following we set c = 1.

II. Formalism
The probability of a transversely polarized quark q ↑ to fragment into an unpolarized hadron h is given by [21] where S ⊥ is the transverse polarization of the quark,k a unit vector with the direction of the quark momentum k, M h is the hadron mass, and D q/h 1 is the polarizationaveraged fragmentation function. Here, the fragmenting quark of flavor q, as well as the identified hadron h in the final state, has been added to the notation of the FFs in order to indicate the dependence of FFs on the final hadron to describe the cross section of back-to-back production discussed below. Equation (1) describes an azimuthal modulation of the hadron momenta around the quark axis, with the strength of the modulation given by the Collins FF H ⊥ 1 . As described in the introduction, a measurement of the effect given by Eq. (1) in single inclusive hadron production in e + e − annihilation, i.e., in the process e + e − → h + X, is not possible due to the chiral oddness of H ⊥ 1 . Instead, the process e + e − → h 1 h 2 | back-to-back +X is considered, where two back-toback hadrons are detected. In this case, the Collins effect can be probed because it appears in a product of two chiral-odd quantities: the quark and antiquark Collins FF. The specific azimuthal modulation is in turn sensitive to the correlation of the transverse polarizations of the produced quark and anti-quark.
The corresponding cross section for inclusive back-toback production of two hadrons can be expressed as with ⊗ signifying convolutions over transverse momenta. The invariant y = (P 1 · l)/(P 1 · (l + l )) can be calculated from the 4-momenta of h 1 , the electron, and the positron, P 1 , l, and l , respectively. The dependence on the quark polarization appearing in Eq. (1) is now contained in the dependence on the azimuthal angles φ 1 and φ 2 , which are measured between the hadron planes and the event plane as shown in Fig. 1. The observable transverse momenta of the hadrons with respect to the thrust axis, which is defined below in Eq. (4), are denoted P ti and serve as a proxy for the parton level P i⊥ .
Equation (2) can be written more compactly as In the e + e − center-of-mass (c.m.) system, used in the following for all calculations, the kinematic factors A and B can be expressed as A = 1 4 (1 + cos 2 θ) and B = 1 4 (sin 2 θ). The angle θ is the angle between the qq axis and the beam axis [22]. Since the transverse projection of the polarization can be calculated in QED as (sin 2 θ)/(1 + cos 2 θ), the appearance of these factors is a reflection of the transverse-polarization dependence of H ⊥ 1 . In a leadingorder partonic picture, the angles φ i would be measured around the qq axis. As this quantity is not accessible, it is approximated by using the thrust axis. The thrust axis is defined as the unit vectorn that maximizes the thrust T : The sum runs over all charged tracks and photons in the event.
Using the thrust axis, it can be determined whether or not the hadrons h 1 and h 2 in a given pair are in different hemispheres ("back-to-back") by requiring for their respective three-momenta P i : The azimuthal angles φ i are calculated as Here,ẑ is the unit vector along the e + beam direction. In the following, the Collins angle of a hadron pair is defined as φ 12 ≡ φ 1 + φ 2 . In terms of φ 12 , the hadron pair yield over all events for a given kinematic bin is given by N 12 ≡ N 12 (φ 12 ). The normalized yield is computed from N 12 by dividing by the average yield: R 12 (φ 12 ) = (N 12 (φ 12 ))/( N 12 ). Considering only a cos(φ 12 ) modulation, R 12 can be parameterized Coordinate system used for this measurement. The thrust axis is denotedn and forms the angle θ in the c.m. system with the beam axis (blue, color online). The thrust axis and beam axis span the event plane. The back-to-back hadrons with momenta P i (i = 1, 2) form the azimuthal angles φi with the event plane. The transverse momenta of the hadrons with respect to the thrust axis are denoted P ti.
as R 12 = 1 + a 12 (θ, z 1 , z 2 , P 2 t1 , P 2 t2 ) cos(φ 12 ), with the azimuthal asymmetry 1 Note that in the expression for a 12 above, the full dependence of the asymmetry a 12 on θ, z i , and P 2 ti is kept. In the measurements presented in this work, at most two variables are kept differential, the other ones are integrated over their accepted ranges.
Measured azimuthal distributions can be strongly distorted due to acceptance and radiation effects. To remedy those effects the double ratio (DR) method can be used. A DR is the ratio of normalized distributions from different kinds of hadron pairs. Under the assumption that the effects are quark-/hadron-flavor independent, they largely cancel in double ratios [14,23,24]. In the previous charged-pion analysis [13,14,16], one double ratio was defined as the ratio of the normalized yield of unlike-sign (π + π − ) to that of like-sign pairs (π + π + and π − π − ). In the current analysis this is extended to include 1 The parameters of the functional forms of the single ratios are denoted with small letters while capital letters are used for the parametrization of the later-introduced double ratios.
neutral mesons: Here, R 0± 12 (R η± 12 , R L 12 ) denote the normalized yields of π 0 π + +π 0 π − (ηπ + +ηπ − , π + π + +π − π − ) pairs and the '+' sign between different combinations means that both pair combinations are considered for the yields. For charged pions, asymmetries of like-sign pairs (L), unlike-sign pairs (U), or pairs that are summed over both charges (C) can be considered. From these combinations the following two double ratios have traditionally been constructed: Analogue to the definition of R L 12 for like-sign pairs, R U 12 and R C 12 denote the normalized yields of the unlike-sign and charge-summed pairs. From R C 12 and R L 12 the double ratio
The double ratios (8)-(10) contain the fragmentation functions of interest in various combinations. To simplify expressions, fragmentation functions are often categorized into favored and disfavored, depending on whether or not the fragmenting-quark flavor is part of the valence structure of the hadron formed. For pions, employing charge and isospin symmetry, the non-strange FFs are [25,26] (11) Besides up and down quarks, the contribution of strange quarks is considered here. 2 Employing the same symmetry arguments as before, the probability for strange-quark fragmentation is the same for all pion states, thus In a similar way the number of FFs for η production can be reduced to Since strange quarks are part of the η valence structure, the respective fragmentation function is not disfavored as is the case of the π 0 fragmentation functions.
The various double ratios can then be expressed in terms of these FFs [25]. Using only the first term of a Taylor expansion in cos(φ 12 ) one obtains R U C 12 ≈ 1 + cos(φ 12 ) sin 2 (θ) 1 + cos 2 (θ) and in particular 2 Charm is qualitatively different due to its mass and the dominance of weak decay channels in pion production. In particular the Collins effect for charm quarks is expected to be small and found so in charm enhanced data samples at Belle and BaBar [13,14,16].
Using Eq. (13) results in the following expression for the η double ratio: In the measurement presented here, a parametrization of the form 1 + A 12 cos(φ 12 ) is fitted to the double ratios. The amplitude A 12 of the cos(φ 12 ) modulation is the azimuthal asymmetry that is presented for various meson combinations and binnings in z and P t .

III. Experiment
The Belle experiment [27] at the KEKB storage ring [28] recorded about 1 ab −1 of e + e − annihilation data. The data were taken mainly at the Υ(4S) resonance at √ s = 10.58 GeV, but also at other Υ(1S) to Υ(5S) resonances and at a continuum setting of √ s = 10.52 GeV. This analysis used data from all these sources for a total integrated luminosity of 980.4 fb −1 . The Belle instrumentation used in this analysis includes a central drift chamber (CDC) and a silicon vertex detector, which provide precision tracking for tracks in 0.30 rad< θ Lab < 2.62 rad, and electromagnetic calorimeters (ECL) [29] covering the same region. The complete ECL consists of 8736 CsI(Tl) counters, which are subdivided into the barrel region (0.56 rad < θ Lab < 2.25 rad) and the endcaps. This analysis uses the barrel ECL for the reconstruction of π 0 and η mesons. Particle identification is performed using information on dE/dx in the CDC, a time-of-flight system in the barrel, aerogel Cherenkov counters in the barrel and the forward endcap, as well as a muon and K L identification system embedded in the flux return steel outside the superconducting solenoid coils. The magnet provides a 1.5 T magnetic field. Using these systems, the selection of charged pions in the barrel, which is used in this analysis, achieves a purity of 97% over all kinematic bins.

IV. Analysis
As in previous similar Belle extractions of azimuthal asymmetries of hadrons and di-hadron pairs [13,14,30], hadronic events are selected by requiring a minimum visible energy of 7 GeV and a thrust T > 0.8. These constraints reduce the contribution of τ leptons and B mesons to below 1% and allow the inclusion of all on-and off-resonance data in the analysis. A number of fiducial constraints are applied in the c.m. system with the goal to minimize effects from variations of the acceptance of the detector on the extracted asymmetries. For this reason only mesons reconstructed from tracks and photons in the barrel region of the detector are considered. Table I lists the fiducial as well as the other constraints applied. This work expands the previous charged-pion analysis [13,14] to π 0 and η mesons, which requires adaptation of several differing or additional selection requirements. They are highlighted in Table I. No correction of the asymmetries for these kinematic restrictions are applied, i.e., the asymmetries extracted are averages in the so-defined phase space.
To minimize the impact of the fiducial constraints on the extracted asymmetry, a hierarchical set of openingangle constraints on photons, hadron momenta, and the thrust axis is applied. This ensures that the detector acceptance of all mesons is radially symmetric around the thrust axis and the acceptance in z and P t of charged and neutral mesons is approximately equal. All photons used for the reconstruction of π 0 and η mesons have a maximal opening angle of 0.5 rad from the thrust axis. All charged and reconstructed neutral mesons used in the asymmetry computation are required to have a maximal opening angle of 0.3 rad from the thrust axis in the c.m. system. Finally, dictated by the geometric acceptance of the ECL, the thrust-axis polar angle is restricted to 1.34 rad < θ < 2.03 rad to ensure the radial symmetry of the acceptance for photons inside the barrel around the thrust axis. To reconstruct π 0 and η mesons, pairs of photons are used for which a minimum energy of 50 MeV and 150 MeV, respectively, is required to reduce background due to combinatorics.
The yields of π 0 and η mesons in each kinematic bin are extracted from a fit to the two-photon invariantmass distribution, with a Crystal-Ball [31] function for the signal and a fifth-order polynomial for the background. The signal to background ratio determined in this way is then used to correct the measured raw asymmetry for the background contribution in the respective kinematic bin in the way described below. Some exemplary fits for π 0 and η mesons are shown in Fig. 2. The measured invariant-mass distributions from experimental data were compared with those from simulations. The simulations used in this analysis employ Pythia [32] TABLE I. Constraints applied in the analysis. The ones that are different in this analysis compared to previous Belle Collins analyses [13,14] are set in bold. (See text for description.)

Description
Constraint Minimum visible energy Evis Evis > 7 GeV Thrust T T > 0.8 Opening angle αO of reconstructed meson w.r.t.n αO < 0.3 rad Thrust axis polar angle θ 1.34 rad < θ < 2.03 rad Minimum photon energy E γ,π 0 for π 0 E γ,π 0 > 50 MeV Minimum photon energy Eγ,η for η Eγ,η > 150 MeV Opening angle αO,γ for photons w.r.t.n αO,γ < 0.5 rad and EvtGen [33] for various physics processes not including the polarization-dependent Collins effect, and GEANT3 [34] for the detector effects. For low-z bins some disagreement between the shape of the invariantmass distributions of reconstructed π 0 s in experimental data and simulation was observed. Therefore an almost non-parametric method, which does not rely on the fit of the signal, was evaluated as well. The method is based on the observation that the background, defined as any pair of electromagnetic clusters in the ECL that do not come from the same π 0 , is well described by the simulation in the sideband region both in magnitude and shape. Hence, instead of fitting the entire invariant-mass spectrum with a background and a signal component, a background description using a quadratic function fitted to 20 points in the upper and lower sidebands obtained from MC, respectively, is used. Once determined in this way, the background is subtracted from the measured invariant-mass spectrum leaving the remaining yield as the signal. The difference between the two extraction methods for the final asymmetry is small, typically less than one per mille in absolute asymmetry value, and is added to the systematic uncertainties. Using the reconstructed π 0 and η mesons, as well as charged pions that are reconstructed using the Belle tracking and particle identification subsystems described in Sec. III, pairs of "back-to-back" hadrons are constructed. This is done by assigning a hemisphere to each meson in the event based on the projection on the thrust axisn and then considering all combinations of hadrons in the first hemisphere with those in the second. Utilizing the thrust axis, the azimuthal angles φ 1 and φ 2 for these "back-to-back" pairs of mesons are computed using Eq. (6).
Double ratios of φ 12 -dependent yields are constructed for the various meson pairs. A cosine function is fitted to the data in order to extract raw asymmetries binned in various combinations of z i and P ti . Here, i = 1 always refers to the neutral meson in the pair when applicable. For pairs of charged pions, the assignment of the first and second pion in a pair is random. Since smearing effects are largest and the Collins effect is smallest at low z, a constraint of z 1 > 0.2 is used, with the exception of the results that are binned in both z 1 and z 2 , where z i > 0.1 is used. The bin boundaries for the P t binning are 0, 0.15, 0.3, 0.5, and 3 GeV. For the binning in z i , bin boundaries differ between results only binned in z 1 and those binned in both z 1 and z 2 . In the for- For the η, due to its higher mass, an additional constraint of z > 0.3 is added for all mesons in the respective pairs.
To arrive at the final asymmetries, several corrections are applied to the raw asymmetries as explained below.
First, the raw asymmetries for π 0 and η mesons are corrected for the contribution from the combinatorial background. The background contribution is determined by calculating asymmetries using γ pairs with a reconstructed mass in the sideband region of the π 0 (η) invariant-mass distribution. Given the limited statistics in this region, four values of the asymmetry are calculated, two in the lower sideband and two in the upper sideband. The observed background asymmetries on both sides of the π 0 (η) signal are consistent with each other and we use a linear fit to extract the contribution of the background to the asymmetry in the signal region using the signal-to-background ratio extracted from the fits to the invariant-mass spectra described earlier.
Second, false asymmetries, determined from simulations, are subtracted. Since the simulation does not contain the Collins effect, any residual asymmetry is a systematic error. These residual asymmetries are consistent with zero within their statistical uncertainties, which are added to our final systematic uncertainties. The relative contribution of these uncertainties ranges from the sub-percent level at low z to a few percent at high z.
Finally, the asymmetries are corrected for thrustsmearing and bin-migration effects. The smearing of the reconstructed z values is negligible due to the excellent momentum reconstruction of the Belle apparatus. In contrast, bin migration is significant for the reconstructed P t . The reason for this is that P t is defined with respect to the thrust axis, the latter suffering from sizable misreconstruction due to particles missed outside of the detector acceptance.
To estimate and correct for the effect of the smearing in Typical two-photon invariant-mass distributions, fit using a Crystal-Ball function for the signal and a polynomial background function, for π 0 (top plots) and η (bottom plots) mesons. In each plot, the green dash-dotted line represents the fitted background using a polynomial of fifth order, the red dashed line the fitted signal, and the blue dotted line is the combined background and signal fit. The combined fit agrees well with the experimental data in black. The vertical dashed lines indicate the boundaries used in the analysis for signal events. P t , a reweighted simulation sample was used. Reweighting the existing simulation is necessary, as the original simulation does not contain the Collins effect. The procedure used weights for each reconstructed hadron pair by assigning a weight w i = 1 + A cos(φ i 12 ), where A is the amplitude of the injected Collins effect and φ i 12 the Collins angle of the i th pair.
The goal of the reweighting of the simulation is the reproduction of the shape of the double-ratio asymmetries observed in the data. The P t dependence of the extracted asymmetries, discussed in more detail in Sec. V, is well described by a linear function in each z bin. Therefore, a (P t1 , P t2 )-dependent amplitude of the form A(P t1 , P t2 ) = 1+a N,D P t1 P t2 was chosen for the reweighting in each z bin. The observed double ratios determine the amplitudes of modulation in the numerator (a N ) and denominator (a D ) only up to a common scaling factor. The dependence of the smearing factor on this scaling factor and on reasonable variations of the ratio a N /a D was observed to be negligible. Using this reweighted simulation, a correction factor f S for each bin is calculated as the ratio of the input double-ratio asymmetries and the reconstructed double-ratio asymmetries. For the former, the generated kinematics of the detected hadrons are used and the thrust axis is computed taking all generated particles in the event into account, including those that are outside of the acceptance of the spectrometer.
The statistical uncertainties in f S contribute to the final systematic uncertainty. Values for f S are between f S = 1.2 and f S = 1.3, with the exception of the kinematic boundaries in the lowest P t bin or when both particles in the pair are in the highest z bin. Here, the hadrons are close to the thrust axis, enhancing smearing effects, and the correction factor takes values between f S = 1.4 and f S = 1.5, depending on the particle species. The relative uncertainty on f S is again driven by the Monte Carlo statistics and is below 2% in the single-z binning, while for the binning in the z values of both hadrons it is below 3% for most bins, but reaches 10% for the highest (z 1 ,z 2 ) bin.
The applied corrections for smearing effects, background contributions, and false asymmetries can be summarized by Here, A raw,bg-corrected is the raw asymmetry after background correction. A MC is the false asymmetry measured in simulation. Finally, the asymmetry is corrected for smearing using the smearing correction f S . Similarly, systematic uncertainties that arise from the statistical uncertainties on the smearing effects, the background contribution, and the false asymmetries can be summarized by Here, δF is the systematic uncertainty stemming from the differences in extracted raw asymmetries using the two different fit procedures.

V. Results and Discussion
Azimuthal asymmetries are measured for double ratios involving charged pions, neutral pions, and eta mesons. Their cosine amplitudes are extracted in various kinematic binnings including z, P t , and a mixed z-P t binning. Significantly non-zero cosine amplitudes are found for all double ratios examined, with magnitudes of mainly a few percent but reaching up to 20% in certain kinematic corners, as pointed out further below.
One novelty of the measurements presented here compared to previous Belle analyses [13,14] is the inclusion of explicit transverse-momentum dependence of the asymmetries. This should help significantly to better constrain the transverse-momentum dependence of the Collins fragmentation function. Figure 3 shows the dependence of both A U C 12 and A U L 12 on the transverse momentum of each of the two pions, where the superscripts U C and U L denote the charge sign combination as defined in (9). In general, A U L 12 is found to be about double the size of A U C 12 , consistent with previous analyses of these asymmetries [13,14,16]. Both asymmetries exhibit a clear rise with increasing, P t1 and P t2 without showing any indication of leveling out at larger values of P t1 and P t2 . In contrast, the largest asymmetry (in this projection) of around 10% for A U L 12 is found in the last (P t1 , P t2 ) bin. This behavior is similar to what was found by BaBar [16], which can be explained perhaps by the limited reach in P t . A direct quantitative comparison of these results with those by BaBar is hampered by the significantly different binning used here. Only in the case of the (z 1 , z 2 ) binning, a few bins at large z 1 and z 2 can be made out that have similar average z and P t . Still, the polar angular range of the thrust axis covered by the two measurements is quite different leading to a sin 2 θ/(1 + cos 2 θ) scaling of the cosine modulations [cf. Eqs. (14)- (16)] that are in variance with each other. However, those are simple scale factors that can be divided out, leaving asymmetries that can be directly compared. In the end, a discrepancy between Belle and BaBar is apparent that cannot be explained easily by charm contributions included here but corrected for at BaBar. Such discrepancy between Belle and BaBar is not new and was observed already before for the large-z region [35]. It is thought to be caused by differences in the applied constraints, e.g., differences in the methodology for removing τ contributions.
Since there are already published results from Belle for charged-pion pairs for the (z 1 , z 2 ) binning, which cover roughly the same kinematic region, a comparison between the results presented here and those from the previous publications [13,14] is provided. The previous results use a smearing correction to correct back to the qq axis extracted from simulation. Since this is not an observable and can be defined cleanly only at leading order, this correction is replaced with a correction back to the thrust axis in the present analysis. Therefore the comparison is performed for asymmetries for which the smearing corrections are removed. This corresponds to a division by the mean smearing correction factor 1.66 for the previous analysis whereas the available bin-by-bin correction is used for this analysis. Further, the compared asymmetry values have been corrected for the kinematic factor sin 2 (θ)/(1 + cos 2 (θ)) bin-by-bin, which differs between the two analyses as a result of the different fiducial constraints. The analysis in Ref. [14] uses a constraint on the z projection of the thrust axis of |T z | < 0.75, which corresponds to 0.72 rad < θ < 2.42 rad. Hence, for the previous analysis the mean kinematic factor is 0.77 whereas it is 0.91 for the presented analysis. The results after adjustments for both the smearing and kinematic factors for the asymmetry values and their uncertainties is the comparison shown in Fig. 4.
There are two further noteworthy differences between the two analyses: (i) The previous analysis does not apply opening-angle constraints. One effect of this difference is that the sampled P t range is different, since high-z hadrons tend to be closer to the thrust axis.
(ii) The previous Belle analysis corrects for the charm contribution using a D * sample. In this analysis, the charm contribution was not corrected for, since using the D * sample can introduce a bias in phase space and introduces larger uncertainties. Instead, the fractional contribution from charm to the event sample is given for each bin in Appendix A, so it can be used for a global extraction.
For the comparison in Fig. 4, it is assumed that the Collins signal coming from charm fragmentation vanishes. In that case, the charm contribution reduces to a simple dilution of the asymmetry of size (1 − f c ), where f c is the ratio of the number of events coming from cc production compared to the sum from cc and light quarks (uds), which in this analysis is extracted from Monte Carlo simulations (see Appendix A for more details). As such the dilution factor can be divided out.
Before discussing the comparison with the previous Belle results, one word of caution on such a charm correction is in place here: The observable of interest in this analysis is the cosine moment of a double ratio, the ]. Clearly, the charm correction sketched above works when both hadron pairs suffer the same amount of dilution. However, it does not work in general when the charm contribution is different for the two hadron pairs, as in that case the dilution factors do not factor out. While this is of a lesser problem for the π 0 asymmetries presented here, as the charm fractions are similar for charged-pion pairs and those involving a π 0 (cf . Tables II-VI), it is certainly more difficult to make this argument for the η asymmetries. It is also for that reason that both the π 0 and η asymmetries discussed further below are not corrected for charm contributions. Figure 5 shows an example comparison of the model-dependent charm fractions in the (z 1 , z 2 ) used for the A π 0 12 and A η 12 asymmetries extracted from the Belle Monte Carlo. Here the superscripts refer to the charge combinations as defined in (8). The charm fractions become small and similar at large z, but deviate from each other for π 0 and η at lower values of z, where the charm fraction gets as large as 20% in the case of π ± η pairs.
Coming back to the comparison presented in Fig. 4, in general a good agreement is visible with the exception of one point in the third z 1 and z 2 bin, which seems to be an outlier. However, a quantification of the agreement is difficult, since the uncertainties of the measurements are correlated. Disregarding this correlation and excluding the outlier, one arrives at a χ 2 per degree of freedom of 1.2. The consistency between the results indicates that the assumption of a vanishing asymmetry for charm quarks is justified.
A second novelty of this measurement is the inclusion of double ratios involving neutral mesons, more specifically π 0 and η. The fragmentation functions for neutral pions are related to those of charged pions through isospin symmetry. Similarly, the η fragmentation functions can be related to those of pions through SU(3) flavor symmetry, which, however, is known to be violated due to the substantially larger mass of strange quarks. Figure 6 displays the dependence of A π 0 12 on z 1 and z 2 . As expected from the charged-pion results, significant asymmetries that rise with z are observed. In the highest  z2) after undoing the different smearing corrections, and after correcting for the different average transverse polarization of the qq pairs in the two measurements due to differences in the θ ranges probed. To make the comparison, the contribution of charm quarks was corrected for by assuming a vanishing charm asymmetry. To avoid confusion with the corrected results, the symbolÂ U L 12 has been used to denote the asymmetry. The lowest z bin was omitted, since the previous analysis used a constraint of z > 0.2. In the figure, data points of the previous analysis are offset horizontally by 0.02 for better visibility.
(z 1 , z 2 ) bin, for which one expects the largest correlation between the fragmenting quark, including its polarization and the final-state hadron, they are reaching 20%. In the lowest z bin, where a large amount of disfavored fragmentation contributes, the asymmetries are consistent with zero within statistical and systematic precision on the sub-percent level.
For the double ratios involving neutral mesons, the asymmetries do not have to be symmetric under interchange of the hadron subscript on z and P t as the neutral meson in the numerator of the double ratios is identified as hadron 1 and the charged pion in the opposite hemisphere as hadron 2. As a result the z 1 and P t1 dependences provide the most sensitivity to the π 0 and η fragmentation functions.
The transverse-momentum dependence is explored in both a mixed z 1 -P t1 binning and a P t1 -P t2 binning. Figure 7 shows the results for A π 0 12 versus z 1 and P t1 , and Fig. 8 the results versus P t1 and P t2 . For P t1 approaching zero, the asymmetry vanishes. The continuous rise with P t1 is consistent with a linear behavior. Higher values of z 1 are again associated with larger values of A π 0 12 , following the same behavior encountered for the charged-pion case.
The results for the η asymmetries have significantly larger uncertainties than those from π 0 . They are extracted from the Belle data imposing a minimum z of 0.3 for both the η and the charged pions involved in the construction of the double ratios. Figure 9 shows the results of A η 12 binned in (z 1 , z 2 ). The rise with z is much less pronounced than the one for charged and neutral pions. Indeed, for the sole z 1 dependence, integrating over P t1 as well as the kinematics of the hadrons in the opposite hemisphere, the asymmetry appears almost constant as shown in Fig. 10. Figure 12 shows the results of A η 12 binned in (P t1 , P t2 ). A clear rise of the asymmetry with transverse momentum can be identified that reaches up to 0.05 for the largest values of P ti . Within large uncertainties, these results for A η 12 are mostly consistent with those of A π 0 12 . In the case of the mixed (z 1 , P t1 ) binning, displayed in Fig. 11, no definite behavior is visible. While clearly rising with P t1 for the last z 1 bin (z 1 > 0.7), the asymmetry is otherwise nearly consistent with a constant, especially as one approaches the lowest z 1 bin. Nevertheless, within the much larger uncertainties the η asymmetries are consistent with the A π 0 12 results, which is shown explicitly in Fig. 13 for the (z 1 , P t1 ) binning, and for which the z > 0.3 requirement was also applied to the π 0 asymmetries. One caveat of this direct comparison is the difference in charm contributions to the π 0 and η, which are about 20-30% larger for the η sample and cannot be eliminated easily as discussed above. On the other hand, for bins with similar enough charm contributions, a comparison is better motivated. Considering Tables II-V, the best candidates appear to be the first few bins in the (P t1 , P t2 ) binning, for which the η and π 0 asymmetries are fully consistent.
Direct extraction of the fragmentation functions for π 0 and η from the double ratio results for comparison with those for charged pions requires further assumptions on the charged-pion fragmentation functions, and is hampered by the complexity of the double ratios. This becomes apparent when recalling the rather involved parton-model expressions (14)- (17) for the various meson combinations. The expression for A π 0 12 is equal to that of A U L 12 −A U C 12 as a result of the isospin relations (11) and (12). Figure 14 displays both A π 0 12 and the difference between A U L 12 and A U C 12 , and indeed good agreement is found. The comparison is to be taken with caution as not all potential correlations between the three asymmetries are taken into account.
The non-vanishing asymmetries for double ratios involving π 0 and η mesons do not necessarily point to nonvanishing Collins fragmentation functions for these two. It is plausible for non-vanishing asymmetries to arise in the case of vanishing Collins functions for π 0 and η due to the presence of the second ratio term in Eqs. (16)  (17), which involves only the charged pions. 3 The first ratio term can be rewritten in terms of products of only π 0 fragmentation functions (in the case of A π 0 12 ) or of π 0 and 3 As a reminder, the second term enters because of using chargedpion pairs in the denominator of the double ratios.
η fragmentation functions (in the case of A η 12 ), i.e., the first ratio is governed by neutral-meson fragmentation functions only, while the second term by charged-pion fragmentation functions. Taking into account that the favored and disfavored pion Collins fragmentation functions are on average of similar magnitude but opposite in sign, thus leading to cancellation effects in the com- bination relevant for the π 0 , a scenerio is plausible in which the π 0 Collins fragmentation is small and the observed signal is due to the term containing the chargedpion fragmentation functions. This is also consistent with the vanishing π 0 Collins asymmetries observed in semiinclusive DIS [36]. The non-vanishing results for A π 0 12 and A η 12 would then mainly be a reflection of the nonvanishing azimuthal modulation in the denominator of those double ratios.

VI. Summary and Conclusion
An analysis of azimuthal asymmetries related to the Collins mechanism has been presented for pairs of backto-back neutral and charged pions as well as η mesons and charged pions. The analysis substantially differs from previous Belle analyses in that results are only presented in the thrust-axis frame without correcting to the qq axis, the opening angle of the hadrons to the thrust axis was limited to 0.3 (which effectively corresponds to a z-dependent upper limit on P t ), and asymmetries were not corrected for charm contributions. Instead, the charm fraction is included and its impact can more properly be treated in future analyses when relevant results on charm azimuthal asymmetries become available, e.g., from Belle II [37]. More importantly, this measurement significantly expands the scope of previous Belle measurements by a) including π 0 and η mesons; and b) exploring the transverse-momentum dependence of the azimuthal asymmetries. Significant asymmetries for all channels are observed. Asymmetries mostly rise, within the given kinematic coverage, with z and P t . The signal for η and π 0 mesons agrees within uncertainties. We show the results for charged-pion pairs agree well with previous Belle measurements [13,14].   The fraction of events originating from charm production is given for the various meson combinations and kinematic binning listed in Tables II-VI. Here, the charm fraction is defined as the ratio of meson pairs that come out from cc production over those coming out of qq (q = u, d, s, c) production as determined from Pythia and EvtGen Monte Carlo simulations employing the Belle default tune. The charm fractions generally are largest at low values of z, reaching fractions as large as 40%, and decrease rapidly with increasing z to a negligible level in the very last z bins. A much milder dependence on P t is observed for all hadron pairs. The fractions are in average larger for pairs involving η mesons compared to those involving only pions. π ± π ± π 0 π ± ηπ ± π 0 π ± (z > 0. TABLE II. Charm fraction in z1 bins. All numbers are in percent. The minimum zi for pions is raised to z1,2 > 0.3 in the last two columns to align with the zi constraint for pairs involving η mesons.    π ± π ± π 0 π ± ηπ ± π 0 π ± (z > 0.  TABLE VI. Charm fraction in (z1, Pt1) bins. All numbers are in percent.

Appendix B. Tables of results
In this section, all asymmetry results are tabulated together with the averages in the kinematic variables z 1 , z 2 , P t1 , and P t2 , as well as of the quantity sin 2 θ/(1 + cos 2 θ), which corresponds to a measure of the size of transverse polarization of the quark-anti-quark pair produced. The tabulated values are obtained from the hadron pairs with the same kinematics that are used to bin the data. Then the average of hadron pairs that appear in the double ratio is taken.