First simultaneous global QCD analysis of dihadron fragmentation functions and transversity parton distribution functions

We perform a comprehensive study within quantum chromodynamics (QCD) of dihadron observables in electron-positron annihilation, semi-inclusive deep-inelastic scattering, and proton-proton collisions, including recent cross section data from Belle and azimuthal asymmetries from STAR. We extract simultaneously for the first time $\pi^+\pi^-$ dihadron fragmentation functions (DiFFs) and the nucleon transversity distributions for up and down quarks as well as antiquarks. For the transversity distributions we impose their small-$x$ asymptotic behavior and the Soffer bound. In addition, we utilize a new definition of DiFFs that has a number density interpretation to then calculate expectation values for the dihadron invariant mass and momentum fraction. Furthermore, we investigate the compatibility of our transversity results with those from single-hadron fragmentation (from a transverse momentum dependent/collinear twist-3 framework) and the nucleon tensor charges computed in lattice QCD. We find a universal nature to all of this available information. Future measurements of dihadron production can significantly further this research, especially, as we show, those that are sensitive to the region of large parton momentum fractions.


I. INTRODUCTION
Phenomenological analyses of hadrons, at the level of quarks and gluons (partons), based on high-energy collision measurements rely on two key ingredients: parton distribution functions (PDFs) and fragmentation functions (FFs).The former provides insight into the momentum-space partonic structure of the incoming nucleon while the latter encodes how a parton forms multiple (colorless) particles in the final state.With regard to FFs, there are two common situations: experiments tagging on one hadron or on a hadron pair (dihadron) within a parton-initiated jet.The single-hadron unpolarized collinear FF D h/i 1 (z) (h denoting the hadron type, i the parton, and z the momentum fraction carried by the hadron) has been studied in observables described by leading-twist collinear factorization and extracted by numerous groups over the last nearly four decades .When transverse momentum dependent (TMD) observables are considered, then, for unpolarized hadrons, not only is D h/i 1 (z, z 2⃗ k 2 T ) relevant ( ⃗ k T giving the transverse momentum of the parton relative to the hadron), but also the chiral-odd Collins TMD FF H T ) [38][39][40][41][42][43][44][45] and [46][47][48][49][50][51][52][53][54] have been extracted using processes described by leading-power TMD factorization.We refer to Ref. [55] for a review on FFs.
On the dihadron side, the FFs depend on more variables [56,57], which can be chosen to be the dihadron total momentum fraction z, the relative momentum fraction ζ, the relative transverse momentum ⃗ R T (whose magnitude can be related to the invariant mass M h of the dihadron [56,57]), along with the parton transverse momentum ⃗ k T .Consequently, even after integrating over ⃗ k T , two dihadron FFs (DiFFs) remain that can be studied within leading-twist collinear factorization [56][57][58][59][60][61]: D h1h2/i 1 (z, M h ) and H ∢ h1h2/i 1 (z, M h ), with the latter being chiral-odd and sometimes referred to as the interference FF (IFF).(Although both D 1 and H ∢ 1 are DiFFs, we will often use the term IFF to distinguish H ∢ 1 from D 1 .)Unlike the single-hadron unpolarized (collinear and TMD) FFs, only one group has extracted these DiFFs [62,63], and only one measurement exists (dσ/dz dM h cross section from Belle for e + e − → (h 1 h 2 ) X [64], which postdates Refs.[62,63]) that is directly sensitive to D h1h2/i 1 (z, M h ).We note that a new definition of DiFFs has recently been introduced [65] that, in particular, allows D h1h2/i 1 (z, M h ) to retain a number density interpretation, from which one can meaningfully calculate expectation values.
With regard to PDFs, at leading twist the longitudinal momentum structure of the nucleon is characterized by the spin-averaged PDF f 1 (x), the helicity PDF g 1 (x), and the transversity PDF h 1 (x), with x being the momentum fraction carried by the parton.The transversity PDF quantifies the degree of transverse polarization of quarks within a transversely polarized nucleon.It is also used to calculate the tensor charges of the nucleon: where h qv 1 ≡ h q 1 − h q 1 are the valence distributions, and µ is the renormalization scale.In addition to phenomenological analyses of experimental data [46-54, 63, 66-69], the tensor charges can be found through ab initio computations in lattice QCD (LQCD) [70? -80] and model calculations [81][82][83][84][85][86][87][88][89][90][91][92].
Given the above discussion, the motivation for this work is the following: 1. Perform the first simultaneous global QCD analysis (JAMDiFF) of the π + π − DiFFs and transversity PDFs (for up and down quarks and antiquarks) from e + e − annihilation, SIDIS, and pp data.
2. Include, for the first time, the Belle cross section data [64], the latest pp dihadron measurements from STAR [96], and all kinematic variable binnings for the relevant processes under consideration, making this the most comprehensive study of dihadron observables to date.
3. Examine the compatibility between phenomenological results for h 1 (x) and/or the tensor charges based on the dihadron approach, those in the TMD/collinear twist-3 framework, and LQCD computations by implementing theoretical constraints at small x [97] and large x [98] (where experimental data is absent) in order to meaningfully calculate the integrals in Eq. ( 1) and include LQCD data into the analysis.
Connected to point (3), the ability to probe certain low-energy beyond the Standard Model (BSM) physics relies on knowledge of the isovector tensor charge g T = δu − δd as well as the individual quark tensor charges δu and δd (see, e.g., Refs.[99][100][101][102][103][104][105][106][107]).Therefore, it is crucial to have precise values for δu, δd, and g T and to test how compatible they are between different approaches used for their determination.To highlight the importance of this issue, we have assembled essential results and discussion regarding our extraction of transversity and calculation of the tensor charges, as well as a comparison to other phenomenological work and lattice QCD, in Ref. [108], with some supplemental details on that aspect of the analysis provided here alongside points (1) and (2).
We organize the paper as follows.In Sec.II we summarize the factorization formulas for all relevant observables, while in Sec.III we discuss the parameterization choices for the DiFFs, IFFs, and transversity PDFs as well as review the Bayesian methodology employed for the analysis.In Sec.IV we compare our theoretical calculations to the data and discuss the quality of the fit.In Sec.V we show and discuss the extracted DiFFs and IFFs as well as calculate expectation values for z and M h .In Sec.VI we present the extracted transversity PDFs and tensor charges, comparing to findings from other phenomenological analyses and LQCD, and providing additional information beyond Ref.[108].We also elaborate on the importance of high-x measurements in further testing the compatibility between phenomenology and LQCD results for the tensor charges.Finally, in Sec.VII we summarize our results and discuss future directions for this analysis, including the role of new experimental data.

II. EXPERIMENTAL OBSERVABLES USED IN THE GLOBAL ANALYSIS
In this analysis we study π + π − dihadron production in e + e − annihilation, SIDIS, and proton-proton collisions with the theoretical formulas given at leading order (LO) in the strong coupling α s using leading-twist collinear factorization.Since we only consider π + π − pairs, we will drop the superscript from our notation for the DiFFs.In the subsections below we summarize each of the three processes and the relevant equations.We note that the formulas in this section use a new definition of the DiFFs that has a number density interpretation [65].

A. Dihadron Production in Electron-Positron Annihilation
In the e + e − annihilation process, e + e − → (π + π − )X, an electron and positron annihilate to form an inclusive spectrum of π + π − pairs with invariant mass M h and fractional energy z.The center of mass (COM) energy is denoted by √ s.The cross section for this process is given by [65] where the sum is over quarks and antiquarks, α em is the fine structure constant, and the hard scale of the process is set to be µ = √ s.In addition, one has where e q is the charge of the quark of flavor q in units of the elementary charge, G F is the Fermi constant, and M Z and Γ Z are the mass and decay width of the Z boson, respectively.Furthermore, a Z = −1 + 4 sin 2 θ W , and a q = 1 − 8/3 sin 2 θ W for q = u, c, ū, c or a q = −1 + 4/3 sin 2 θ W for q = d, s, b, d, s, b, where θ W is the weak mixing angle.As we will argue below (see Eq. ( 13) and Appendix B), there are five independent DiFFs (D u 1 , D s 1 , D c 1 , D b 1 , and D g 1 ), and this observable is only capable of constraining one of them (which we choose to be D u 1 ).Thus, we supplement the experimental data from Belle with information from PYTHIA [109], which will be discussed in more detail in Sec.IV A.
We also consider the process e + e − → (π + π − )(π + π − )X where two dihadron pairs are detected, with the second having invariant mass M h and fractional energy z.By measuring an azimuthal correlation of two hadron pairs detected in opposite hemispheres, one obtains an observable that is sensitive to the IFF H ∢ 1 (z, M h ).This modulation (known as the Artru-Collins asymmetry [110] and denoted as a 12R by Belle [111]) is given by [62,63,110,[112][113][114]] where θ is defined as the polar angle between the beam axis and the reference axis in the COM system [112].As we will argue later (see Eq. ( 18) and Appendix B), there is only one IFF that is needed for phenomenology, which we choose to be H ∢,u 1 .Thus, only one observable is needed to constrain the IFFs.However, one can see from Eq. ( 4) that A e + e − is proportional to [H ∢,u 1 ] 2 .Consequently, the asymmetry cannot uniquely determine the sign of H ∢,u 1 , and it must be fixed by hand.Given the sign of the SIDIS asymmetries from experimental measurements (see the following section) and assuming that the transversity up quark PDF must be positive (as is found in all phenomenological, model, and lattice QCD studies) leads to the conclusion that H ∢,u 1 must be negative (which is also supported by model calculations for H ∢ 1 [115]).In the analysis we therefore choose H ∢,u 1 to be negative.

B. Dihadron Production in SIDIS
The SIDIS process is given by ℓN ↑ → ℓ ′ (π + π − )X, where a lepton with 4-momentum k scatters off a transversely polarized proton (N = p) or deuteron (N = D) with 4-momentum P and is detected with 4-momentum k ′ along with a π + π − pair with 4-momentum P h .The 4-momentum transfer squared is given by Q 2 ≡ −(k − k ′ ) 2 , and we define the Bjorken scaling variable x bj ≡ Q 2 /(2P • q) and the inelasticity y ≡ P • q/P • k.We denote the virtualphoton 4-momentum and relative 4-momentum between the two hadrons (divided by 2) by q and R, respectively.The transverse component ⃗ R T of R makes an azimuthal angle ϕ R about the virtual-photon direction, and the transverse nucleon spin ⃗ S T makes an azimuthal angle ϕ S .COMPASS [116] defines the angle ϕ RS as ϕ RS ≡ ϕ R + ϕ S − π, while HERMES [117] defines it as ϕ RS ≡ ϕ R + ϕ S .
The asymmetry (usually denoted as A sin (ϕ R +ϕ S ) sin θ U T [117] or A sin (ϕ R +ϕ S ) U T [116]) can be written as [57,63,66,118,119] The factor c(y) is given by The opposite sign between the asymmetries for COMPASS and HERMES comes from the different definitions for the angle ϕ RS mentioned above.The hard scale is set to be µ = Q, and the unpolarized PDF f 1 is taken from Ref. [36] (we utilize only the mean value).
Both COMPASS and HERMES measured the asymmetry for proton targets, while COMPASS also used a deuteron target.For the deuteron target we neglect nuclear corrections and take the (unpolarized and transversity) PDFs in Eq. ( 5) to be the average of proton and neutron PDFs, using isospin symmetry to relate the neutron PDFs to those of the proton.Due to the symmetry of the π + π − IFFs (see Eq. ( 18)), one has in the numerator for a proton target q e 2 q h q 1 H ∢,q while for the deuteron one has q e 2 q h q/D 1 where h q/D 1 are the deuteron transversity PDFs.From this one immediately sees that the SIDIS data is sensitive only to the valence transversity PDFs [63,66,69], with the proton asymmetry primarily sensitive to u v and the deuteron asymmetry equally sensitive to u v and d v .

C. Dihadron Production in Proton-Proton Collisions
In the process p ↑ p → (π + π − )X, a transversely polarized proton with 4-momentum P A collides with an unpolarized proton with 4-momentum P B producing an outgoing dihadron pair with 4-momentum P h (with P hT the magnitude of the transverse component) and pseudo-rapidity η.The measured asymmetry (denoted as A U T by STAR [120]) is given by [121] A pp U T = H(M h , P hT , η) where the numerator and denominator read [122] H(M h , P hT , η) = 2P hT i a,b,c,d with x a (x b ) the momentum fraction of the parton coming from the proton with momentum P A (P B ).We note that A pp U T is sensitive to both the quark and antiquark transversity PDFs.The sum i is over all partonic sub-processes, with a,b,c,d summing over all possible parton flavor combinations for a given channel.The relevant Mandelstam variables are s = (P A + P B ) 2 , t = (P A − P h ) 2 , and u = (P B − P h ) 2 .The total momentum fraction of the dihadron is fixed at LO to be The hard scale is set to be µ = P hT .The hard (perturbative) partonic cross sections dσ and d∆σ [122] depend on ŝ = x a x b s, t = x a s/z, and û = x b u/z.The explicit expressions for d∆σ can be found in Appendix A. The limits on the integration are III. PHENOMENOLOGICAL METHODOLOGY In this section we discuss the parameterization of the DiFFs, IFFs, and transversity PDFs.All functions are evolved using the DGLAP equations [123][124][125] appropriate for the transversity PDFs [126][127][128][129] and DiFFs/IFFs [65] with LO splitting functions.We mention that evolution equations were previously derived for "extended" DiFFs (at LO) in Ref. [130] (and commented on in Ref. [113]) and for collinear DiFFs (functions of (z 1 , z 2 )) in Refs.[131][132][133][134] at LO and recently in Refs.[135,136] at next-to-leading order (NLO).The strong coupling α s is evolved at leading-logarithmic accuracy with the boundary condition α s (M Z ) = 0.118.We use the zero mass variable flavor scheme for DGLAP evolution.The values of the heavy quark mass thresholds for the evolution are taken from the PDG as m c = 1.28 GeV and m b = 4.18 GeV in the MS scheme [137].

A. DiFF Parameterization
The DiFF D 1 for π + π − production satisfies (see Appendix B and Ref. [62]) which means, including the gluon, there are five independent D 1 functions to be fitted.D u 1 , D s 1 , and D g 1 are parameterized at the input scale µ 0 = 1 GeV, while D c 1 and D b 1 are parameterized at µ = m c and µ = m b , respectively.The primary challenge with extracting DiFFs is their dependence on both z and M h .Here we choose to discretize the M h dependence while using an explicit functional form for the z dependence.Such a choice has the computational advantage of being convertible to Mellin space.For D u 1 , we choose the following M h grid: We note that the grid is not uniform and is instead chosen in a way to best describe the detailed resonance structure of the e + e − cross section [64].For s, c, b and g we choose grids that are less dense: We find that these grids, and the parameters associated with them (see Eq. ( 15)), are necessary and sufficient to describe the Belle and PYTHIA data.At each value of M h on the grid, denoted by M i,j h , the z dependence is parameterized as where i = u, s, c, b, g, a ≡ {N i jk , α i jk , β i jk } is the set of parameters to be inferred, and , where B(a, b) is the beta function, normalizes the function to the first moment to largely decorrelate the normalization and shape parameters.For the up quark, it is necessary to include terms up to k = 3 in order to describe the Belle cross section data.The other functions will be constrained by PYTHIA-generated data (see Sec. IV A) and we find that only the k = 1 term is necessary to describe them.This leads to 14 × 9 = 126 parameters for D u 1 , 7 × 3 = 21 parameters each for D s 1 and D c 1 , 5 × 3 = 15 parameters for D b 1 , and 4 × 3 = 12 parameters for D g 1 , for a total of 195 free parameters.The parameters are generally restricted within the ranges 0 < N i jk < 15, −2 < α i jk < 10, and 0 < β i jk < 10, but in practice these ranges may be restricted further depending on i, j, and k.These ranges are chosen to give the fit sufficient but not unnecessary flexibility.
The functional form in Eq. ( 15) is easily converted into Mellin space at each M i,j h .The inverse Mellin transform requires a complex contour integration (see Eq. ( 33)), which is performed by discretizing the contour.To obtain the function at any value of M h , we interpolate both the real and imaginary components of the function in Mellin space using the values calculated at M i,j h .Once interpolated, the real and imaginary components are recombined and the inverse Mellin transform is performed to obtain the function in momentum space at any value of M h .
We contrast our choice of parameterization with that of Ref. [62], which used a considerably more complicated functional form for u, s, and c with different functions for the continuum and three resonant channels.The bottom quark DiFF was not included in the analysis of Ref. [62], and the gluon DiFF was only generated perturbatively.In Ref. [68] three choices for the gluon DiFF at the input scale were used: Our results in Fig. 11 suggest that none of these scenarios seem to be correct across all z and M h .Overall, our parameterization has the advantages of requiring fewer assumptions about the functional forms, reducing complicated correlations between parameters, and being convertible to Mellin space.The latter two features make the analysis far more computationally efficient.
We also enforce the positivity bound [57] approximately on each Monte Carlo replica by using a Bayesian prior (see Sec. III D) that in effect imposes a penalty on the χ 2 function when the bounds are violated [138].For each replica and at each step of the χ 2 minimization we first calculate D 1 at the input scale µ 0 = 1 GeV at 300 points spaced linearly in the (z, M h ) plane, with 0.2 < z < 1 and 2m π < M h < 2.0 GeV, noting that if positivity is enforced at the input scale it will automatically hold at larger scales for any function that evolves through the DGLAP equation [139,140].We repeat this process for u, s, c, b, and g.Any DiFFs that are negative contribute to the overall where the sum z,M h represents summing over the 300 (z, M h ) points, and Θ(X) is the step function.If the contribution is non-zero, it is proportional to the size of the violation.The weight w is chosen to be 3 so that the initial contribution to the χ 2 is generally O(1000), although the size of the initial contribution can vary significantly depending on the starting parameters.This ensures that the contribution is non-negligible but that it also does not dominate the entire χ 2 function which otherwise is also O(1000).This choice ensures that the violation of the bound is small.

B. IFF Parameterization
The IFF H ∢ 1 for π + π − production satisfies (see Appendix B and Ref. [62]) Consequently, there is just a single independent IFF (since the transversity cannot couple to the IFF for gluons at LO), which we choose to be H ∢,u 1 .The IFF H ∢,u 1 is parameterized similarly to D i 1 .Since the relevant data for H ∢,u 1 is comparatively sparse and has larger errors, far fewer parameters are needed.Thus, we are able to choose a less dense M h grid M u h = [2m π , 0.50, 0.70, 0.85, 1.00, 1.20, 1.60, 2.00] GeV, and at each M u,j h , the z dependence is given by at the input scale µ 0 = 1 GeV.This leads to a total of 8 × 6 = 48 free parameters for the IFF.The parameters are restricted within the ranges −1 < N u jk < 0 (see discussion below Eq. ( 4)), −1 < α u jk < 5, and 0 < β u jk < 10 for all j and k.As with D i 1 , this form is easily converted into Mellin space.We also enforce the positivity bound [57] which, due to the symmetry relations Eqs. ( 13) and (18), only needs to be applied to i = u.Since we have fixed D u 1 to be positive and H ∢,u 1 to be negative, we only need to check the condition D u 1 > −H ∢,u 1 .This is implemented similarly to the positivity constraint on D 1 in Eq. ( 16), with a Bayesian prior (see Sec. III D) that in effect leads to a χ 2 penalty equal to with w again chosen to be 3 (see discussion below Eq. ( 17)).

C. Transversity PDF Parameterization
The transversity PDFs are parameterized at the input scale µ 0 = 1 GeV using the form normalized to the first moment . We choose to parameterize the valence distributions h uv 1 and h dv 1 , as well as the antiquark distributions h ū 1 = −h d 1 .Since we only have three unique observables to constrain the transversity PDFs (proton SIDIS, deuteron SIDIS, and pp collisions), we choose the relation between the antiquarks based on predictions from the large-N c limit [141].We note that only one previous phenomenological analysis has included antiquark transversity PDFs, which was an exploratory study in Ref. [54] that found h ū 1 and h d 1 to be small and consistent with zero (although had set h ū 1 = h d 1 , at variance with the large-N c limit).The N parameter is restricted between 0 < N i < 1 for i = u (see discussion below Eq. ( 4)) and −1 < N i < 1 for i = d, ū, while the other parameters are always restricted between 0.09 < α i < 0.26 (see discussion below surrounding Eq. ( 23)), 0 < β i < 20, −20 < γ i < 20, and −20 < δ i < 20.We have tested the model dependence of our input transversity functions by adding a second shape of the form in Eq. ( 22) to all three quark flavors.We find that the results only change marginally, indicating that the experimental data do not call for a more expressive form for these PDFs.We caution, however, that this test is not exhaustive of all possible functional forms one could choose for the transversity functions.We have a total of 3 × 5 = 15 parameters for the transversity PDFs.Combined with the 195 parameters for D 1 , the 48 parameters for H ∢ 1 , and 7 normalization parameters (see Eq. ( 28) below), we end up with a total of 265 fitted parameters.
We also place a constraint on the small-x behavior of the transversity PDFs (governed by the α i parameter in Eq. ( 22)), where no experimental data is available.Theoretical calculations have placed limits on this parameter as x → 0 (ignoring saturation effects) [97]: We apply this limit to both the valence quarks and antiquarks, and there is a roughly 50% uncertainty on this value from 1/N c and NLO corrections [142].At the input scale, we calculate α i → 0.17, and so we limit the α i parameter to the range 0.085 ≤ α i ≤ 0.255 for all quark flavors.Strictly speaking, this limit only applies as x → 0, while our approach places a limit on the entire range of x.We find, however, that limiting α i as such has no impact on the resulting transversity PDFs in the measured region or on our ability to describe the experimental data.Thus, this simplified approach is sufficient to capture the x → 0 behavior while not affecting results at moderate or high x.
We also enforce the Soffer bound [98] on the transversity PDFs, given by where f i 1 and g i 1 are taken from Ref. [36].This is again enforced through a Bayesian prior (see Sec. III D) that in effect leads to a χ 2 penalty, similar to Eqs. ( 17), ( 21): where i = u, d, ū, d and The sum over ± is due to the fact that the positivity bound applies to the absolute value of h i 1 .The sum x represents summing over 100 points in x linearly spaced in the range 0.001 < x < 0.99.The weight w is chosen to be 10 (see discussion below Eq. ( 17)).In order to account for the errors on f 1 and g 1 , we add their 1σ errors in quadrature, add this error to the mean of the Soffer bound, and use these values in the fit.

D. Bayesian Analysis
Our Bayesian analysis consists of sampling the posterior distribution given by with a likelihood function of Gaussian form, and a prior function π(a).The χ 2 function in Eq. ( 27) is defined for each replica as where d e,i is the data point i from experimental dataset e, and T e,i is the corresponding theoretical value.Note that the PYTHIA and LQCD data are not included in this sum.All uncorrelated uncertainties are added in quadrature and labeled by α e,i , while β k e,i represents the k th source of point-to-point correlated systematic uncertainties for the i th data point weighted by r e,k .The latter are optimized per values of the parameters a via ∂χ 2 /∂r e,k = −2r e,k , which can be solved in a closed form.We include normalization parameters N e for each dataset e as part of the posterior distribution per data set.
The prior π(a) for each replica is given by where l is a product over the parameters, and f is a product over the LQCD and PYTHIA datasets.The step function Θ(X) forces the parameters to be within the chosen range between a min l and a max l , which also automatically enforces the small-x constraint Eq. (23).A Gaussian penalty controlled by the experimentally quoted normalization uncertainties δN e is applied when fitting N e .The final three exponentials utilize Eqs. ( 17), (21), and (25).
The posterior distribution is sampled via data re-sampling, whereby multiple maximum likelihood optimizations of P(a|data) are carried out by adding Gaussian noise with width α e,i to each data point across all data sets (including the PYTHIA and LQCD priors).The resulting ensemble of replicas {a k ; k = 1, . . ., n} is then used to obtain statistical estimators for a given observable (which are functions of DiFFs/PDFs) O(a), such as the mean and variance, The agreement between data and theory is assessed by using the "reduced" χ 2 which is defined as where N dat is the number of data points under consideration and E[• • • ] represents the mean theory as defined in Eq. (30).We note that in Eq. ( 31) we also consider the LQCD and PYTHIA-generated data in order to quantify their compatibility with the experimental data.

E. Mellin Space Techniques
We have chosen parameterizations for all the functions in a way that can be converted to Mellin space, which has two advantages.First, as is well known, it allows us to solve the DGLAP evolution equations analytically rather than using an iterative solution as is required in momentum space.We now demonstrate the second advantage of this approach when calculating the pp asymmetry [143].One can relate the DiFFs, IFFs, and transversity PDFs in momentum space and Mellin space through where F (N ; µ 2 ) is the N th moment of F (x; µ 2 ) and F (M, M h ; µ 2 ) is the M th moment of F (z, M h ; µ 2 ).The integration is over a contour in the complex N plane that intersects the real axis to the right of the rightmost poles of F (N ; µ 2 ).Then, the numerator Eq. ( 9) and denominator Eq. ( 10) of the pp asymmetry can be written as Since we do not fit f 1 and instead take it from Ref. [36], the quantities in brackets are independent of any fitting parameters and so only need to be calculated once prior to the actual fit.Thus, the integrals over x a and x b need not be calculated at every step of the fitting procedure, massively reducing the computation time, especially for the double convolution in the numerator.We take advantage of this computational efficiency, allowing us to perform a simultaneous fit of the DiFFs, IFFs, and transversity PDFs.

IV. GLOBAL ANALYSIS DATA AND QUALITY OF FIT
In this section we discuss the experimental and LQCD results that enter the analysis as well as the PYTHIAgenerated data.

A. PYTHIA-Generated Data
This analysis must be supplemented by PYTHIA [144,145] data due to the fact that there are five independent D 1 functions (see Sec. III A), but only one experimental observable (e + e − cross section dσ/dz dM h ) currently available to constrain them.(Although D 1 appears in the denominators of the SIDIS and pp asymmetries, those observables are being used to constrain the transversity PDFs and offer no significant sensitivity to D 1 .)Thus, we will discuss how we generate the PYTHIA data and how it is used to help constrain D i 1 (z, M h ).Our goal with the PYTHIA-generated data is to provide reasonable constraints on the D 1 functions for the strange, charm, bottom, and gluon in the absence of experimental data beyond that from Belle.To do this, we first generate data at the Belle energy √ s = 10.58GeV for the ratio σ q /σ tot with q = s, c, b, where σ tot ≡ q σ q with the sum over all quark flavors and σ q is the flavor-dependent cross section.This data provides priors that constrain the strange, charm, and bottom DiFFs.In order to provide some constraint on D This strategy allows us to gain sensitivity to D g 1 scaling violations.We do not include the bottom quark contribution at the lowest Belle energy √ s = 10.58GeV but include it at all energies above this, due to the fact that our massless quark formalism would be highly inadequate at such low energies comparable to 2m b .
This leads to a total of 14 generated datasets from PYTHIA, with each one covering the same (z, M h ) region as the Belle data.The same cuts used on the real Belle data are used on the PYTHIA data.The cut in Eq. ( 35) below is used with √ s = 10.58GeV regardless of the value of √ s, as there is no benefit within our analysis of going beyond the kinematic region of the Belle data.
The PYTHIA [144,145] data can be generated with arbitrarily high statistics, so we neglect the tiny statistical error.In order to propagate model uncertainties in PYTHIA, we use a collection of tunes ("PYTHIA 6 default", "PYTHIA 6 ALEPH", "PYTHIA 6 LEP/Tevatron", "PYTHIA 8") to estimate the systematic uncertainties on the flavor-dependent cross sections.This list of tunes is similar to that used in Ref. [64].The tunes are summarized in Table IV in Appendix C. The resulting spread of the flavor-dependent cross sections is used to estimate the errors due to the tunes, with the mean taken as the central value and the minimum and maximum defining the uncorrelated Collaboration Ref.
Observable Process Non-perturbative function(s) Belle [64] dσ/dz dM h e + e − → (π systematic error.We choose to fit the ratio σ q /σ tot as it leads to smaller variance between the different tunes compared to taking the absolute cross sections σ q .Taking such a ratio should also help to cancel NLO and thrust-cut effects.The χ 2 red for the PYTHIA data is shown in Table V in Appendix C, along with comparisons between the PYTHIA-generated data and theory calculations.

B. Experimental Data and Lattice QCD Results
A summary of the experimental data and lattice results included in our analysis is given in Table I, where we also indicate which fitted non-perturbative functions enter the observables.The kinematic coverage of the measurements is shown in Fig. 1.We use the π + π − production in e + e − annihilation cross section data from Belle [64] at √ s = 10.58GeV.(We note that there is a typo in the Belle publication regarding the units of the cross section.The units are labeled as µb/GeV but the cross section is actually shown in units of nb/GeV.)Since we can only extract information on the IFF up to M h = 2.0 GeV (see below), we place a cut of M h < 2.0 GeV on the cross-section data.The data are provided in bins of M h , with width ∆M h = 0.02 GeV, and bins of z, with width ∆z = 0.05.We evaluate the cross section at the mean of each M h bin ⟨M h ⟩ and z bin ⟨z⟩.In order to avoid the sharp kaon and D 0 resonances, we also cut out all data points with ⟨M h ⟩ = 0.49 GeV and ⟨M h ⟩ = 1.87 GeV.The kinematic limit of the Belle data is given by M h < √ s 2 z, which is reached in the bins of z below z = 0.45.We place a cut of 2⟨M h ⟩ ⟨z⟩ √ s < 0.7 (35) in order to avoid the region near the kinematic limit.We also cut out the lowest z bin 0.20 < z < 0.25 as we find difficulty in describing the data in this bin, which may be an indication of the need to resum small-z logarithms.In total, with these cuts we include 1,094 of the 1,468 data points provided by Belle, and the kinematic ranges are 0.25 < z < 1.0 and 0.3 < M h < 2.0 GeV.The data has an overall normalization uncertainty of 1.6%, and the systematic errors are uncorrelated point-to-point [64].
We include the e + e − annihilation Artru-Collins asymmetry data provided by Belle [111].The same cuts on M h and z discussed above are applied to both hadron pairs.The data is available in three binnings: (z, M h ), (M h , M h ), and (z, z).We include all three binnings in our analysis and evaluate the theory formulas using the mean values ⟨M h ⟩, ⟨z⟩, ⟨M h ⟩, ⟨z⟩.The systematic errors on the data are uncorrelated point-to-point and there is no overall normalization error [111].We note that in all Belle measurements in this analysis, a cut is placed on the thrust T [146] of T > 0.8.In our LO formalism, it is not possible to take this thrust cut into account because T < 1 can only be generated via NLO effects.We expect that the errors from this omission largely cancel out in the ratio of H ∢ 1 /D 1 which appears in all of the asymmetries, thus reducing the impact on the extraction.However, since the extraction of D 1 relies on cross-section data rather than asymmetries, its extraction could be noticeably influenced by next-to-leading order (NLO) corrections.
The transversity PDFs are constrained by SIDIS asymmetry and pp asymmetry data.SIDIS data is available from HERMES [117] and COMPASS [116].We use all three available binnings in x bj , z, and M h , and, as for the e + e − annihilation data, apply the cuts ⟨z⟩ > 0.25 and ⟨M h ⟩ < 2.0 GeV.We evaluate the theory formulas using the mean values ⟨x bj ⟩, ⟨z⟩, ⟨M h ⟩, ⟨y⟩.The HERMES data has an overall normalization uncertainty of 8.1%, while the COMPASS data has 2.2% for proton and 5.4% for deuterium.The systematic errors on the HERMES data are correlated point-to-point [117], while those on the COMPASS data are uncorrelated [116].
The pp collision data is provided by STAR, at both √ s = 200 GeV [120] and √ s = 500 GeV [96].We apply the cut ⟨M h ⟩ < 2.0 GeV to all of the data.The 0.3, 0.4) on the opening angle R of the pion pair, with 0.3 treated as the default.This cut is used to filter out pion pairs that do not originate from a single parton.We use the data corresponding to R < 0.3 and have verified that the changes to the resulting functions by instead using data from the other cuts are negligible compared to the statistical uncertainties.The √ s = 500 GeV data is provided with an opening angle of R < 0.7.A larger opening angle cut is acceptable here, as the increased energy means that gluon radiation is occurring at wider angles, allowing the dihadron pair to still be considered as originating from a single parton even with a larger R-cut value.The data is provided binned in P hT , M h , and η, with the results (often) provided for both η > 0 and η < 0 when binned in P hT or M h .We include all binnings and take the central values ⟨P hT ⟩, ⟨M h ⟩, and ⟨η⟩ when evaluating our theory formulas.All systematic errors are uncorrelated, and the 200 GeV and 500 GeV data have normalization uncertainties of 4.8% and 4.5%, respectively [96,120].
We also consider the inclusion of LQCD data on tensor charges in the fit and treat them as Bayesian priors (see Sec. III D).We restrict ourselves to results at the physical pion mass with 2 + 1 + 1 flavors, where calculations are available from ETMC [77] and PNDME [71] on δu, δd, and g T .We choose to fit δu and δd rather than g T in order to provide flavor separation.The reported uncertainties are treated as uncorrelated.

C. Data vs. Theory
We now discuss the results of our global analysis, where we consider three different scenarios.The first includes all experimental and LQCD data discussed above and will be referred to as "JAMDIFF (w/ LQCD)."The "JAMDIFF (no LQCD)" fit removes the LQCD data.The "JAMDiFF (SIDIS only)" fit excludes all of the STAR pp and LQCD data, and also sets the antiquark transversity PDFs to zero as they cannot be constrained in this scenario.The χ 2 red for the three scenarios is summarized in Table II.In the plots that follow comparing our theory results with the experimental measurements, we show the JAMDiFF (w/ LQCD) fit.We will reserve a discussion about our calculated tensor charges and comparison to those from LQCD computations (and other values from phenomenology) to Sec.VI.We re-emphasize that we have performed a simultaneous global analysis of DiFFs/IFFs and transversity PDFs, where, unlike previous work [63,66,68,69], the parameters for the DiFFs are not fixed (from a fit of only e + e − annihilation) but allowed to be free along with the transversity PDF parameters.We have also studied an exhaustive set of available data on dihadron observables, which includes, for the first time, the Belle cross section [64], the √ s = 500 GeV measurements from STAR [96], and all kinematic variable binnings for the relevant processes under consideration, amassing 1471 experimental data points.We are able to describe all of the experimental data very well.
The results for the Belle cross section are shown in Fig. 2. We find that our parameterization discussed in Sec.III A is able to account for the resonance structure in the data and the general behavior of the measurement across the kinematic range in z and M h , with χ 2 red = 1.01.There are some slight discrepancies in the lowest z bin at small M h .The fitted normalization N e (see Eq. ( 28)) for this dataset is 1.00 (1).In Figs.3-5 we display the Belle Artru-Collins asymmetry results.Again our theory curves are able to describe the data very well across all kinematics, with χ 2 red = 1.27, 0.60, 0.42 for the (z, M h ), (M h , Mh ), (z, z) binnings, respectively.The only discrepancy is in the 0.77 < M h < 0.90 GeV bin at high z, which is what causes the χ 2 red for the (z, M h ) binning to be much larger than the other two.There is no fitted normalization for this dataset.As seen in Table II, the χ 2 red for the Belle data (both the cross section and Artru-Collins asymmetry) is nearly identical across the three fit configurations.
The results for the SIDIS asymmetry are shown in Fig. 6.The fitted normalizations N e (see Eq. ( 28)) are 0.98(4) for HERMES, 1.014 (6) for COMPASS proton, and 1.007 (8) for COMPASS deuteron.The theory calculations are generally in agreement with the data for all three projections (x bj , M h , z).However, we focus the reader's attention on the highest x bj COMPASS proton point (data/curve in red in the left panel of Fig. 6).The trend is for JAMDiFF (w/ LQCD) to increase at larger x bj more rapidly than the COMPASS proton data.The JAMDiFF (no LQCD) fit follows the trend of the data much better (see Fig. 19 in Sec.VI B), which is apparent in the reduction of χ 2 red for the COMPASS proton x bj binning from 1.98 to 0.65 when the lattice data is removed.In Sec.VI B we will elaborate on the importance of high-x experimental measurements in testing the compatibility between phenomenology and LQCD results for the tensor charges.We lastly remark that the SIDIS only fit has χ 2 red values similar to the analyses that included pp data, demonstrating clear compatibility between the two reactions.
The results for the pp asymmetry are shown in Figs.28)) are 1.00(3) for the 200 GeV STAR data and 1.14 (6) for the 500 GeV STAR data.As with the e + e − and SIDIS observables, we again find generally good agreement with the pp reaction across all kinematic binnings.The one discrepancy of note is with the η binning for the 500 GeV STAR data at the highest η value (see Fig. 10) that has more sensitivity to larger x.The JAMDiFF (w/ LQCD) fit overshoots that point and has a trend of continuing to increase with increasing η.Data at larger η are needed in order to test this behavior.The JAMDiFF (no LQCD) fit still increases with η but falls below the JAMDiFF (w/ LQCD) curve (see Fig. 19).This accounts for the improvement in the χ 2 red for the 500 GeV STAR data in the η binning, going from 2.97 to 1.83, once the LQCD data is removed.Interestingly, the 200 GeV STAR data for the η binning improves when LQCD data is included, going from 1.46 to 0.52.So it seems that, within our analysis, there are competing trends at larger η in pp that measurements at more forward rapidities may be able to resolve.
FIG. 2. dσ/dz dM h cross section data from Belle [64] at √ s = 10.58GeV (blue circles) plotted as a function of M h against the mean JAMDiFF (w/ LQCD) result (using Eq. ( 2)) with 1σ uncertainty (red lines with bands).Each panel is a different bin of z. [111] binned in (z, M h ) at √ s = 10.58GeV (blue circles) plotted as a function of z against the mean JAMDiFF (w/ LQCD) result (using Eq. ( 4)) with 1σ uncertainty (red lines with bands).The different panels show different bins of M h .FIG. 4. Artru-Collins asymmetry data from Belle [111] binned in (M h , M h ) at √ s = 10.58GeV (blue circles) plotted as a function of M h against the mean JAMDiFF (w/ LQCD) result (using Eq. ( 4)) with 1σ uncertainty (red lines with bands).The different panels show different bins of M h .FIG. 5. Artru-Collins asymmetry data from Belle [111] binned in (z, z) at √ s = 10.58GeV (blue circles) plotted as a function of z against the mean JAMDiFF (w/ LQCD) result (using Eq. ( 4)) with 1σ uncertainty (red lines with bands).The different panels show different bins of z.FIG. 6. SIDIS asymmetry data from HERMES [117] (green triangles) and COMPASS [116] (red circles for proton, blue squares for deuteron) plotted against the mean JAMDiFF (w/ LQCD) result (using Eq. ( 5)) with 1σ uncertainty (colored lines with bands).The data is binned in x bj (top left), M h (bottom left), and z (bottom right).We note that the asymmetries from HERMES and COMPASS are defined such that they have opposite signs (see Eq. ( 5) and the surrounding discussion).FIG. 7. Proton-proton asymmetry data from STAR at √ s = 200 GeV [120] with opening angle cut R < 0.3 plotted against the mean JAMDiFF (w/ LQCD) result (using Eq. ( 8)) with 1σ uncertainty (colored lines with bands).The left panel shows the results binned in M h while the right panel shows them binned in P hT , for both η > 0 (red circles) and η < 0 (blue squares).

A pp U T
√ s = 500 GeV (R < 0.7) Proton-proton asymmetry data from STAR at √ s = 500 GeV [96] with opening angle cut R < 0.7 (red circles for η > 0, blue squares for η < 0) plotted as a function of M h against the mean JAMDiFF (w/ LQCD) result (using Eq. ( 8)) with 1σ uncertainty (colored lines with bands).The different panels show different bins of P hT .

A pp U T
√ s = 500 GeV (R < 0.7) STAR η > 0 JAMDiFF (w/ LQCD) FIG. 9. Proton-proton asymmetry data from STAR at √ s = 500 GeV [96] with opening angle cut R < 0.7 and η > 0 (red circles) plotted as a function of P hT against the mean JAMDiFF (w/ LQCD) result (using Eq. ( 8)) with 1σ uncertainty (red lines with bands).The different panels show different bins of M h .FIG. 10.Proton-proton asymmetry data from STAR (black circles) plotted as a function of η against the mean JAMDiFF (w/ LQCD) result (using Eq. ( 8 36)), with M h integrated from 2mπ to 2 GeV.Right panel: Similar to the left panel but for ⟨z|M h ⟩ i (see Eq. ( 37)), with z integrated from 0.25 to 1.
Similarly, the average value of z for a single π + π − pair with a given M h formed from the fragmentation of a parton i is These quantities are shown in Fig. 13.The general feature is that the average M h (z) increases as z (M h ) becomes larger.This can perhaps be understood by the fact that the number of dihadrons produced is greater at smaller z and decreases with increasing z (see Fig. 11).We notice from the left panel of Fig. 13 that the u, d = u, s, and c quarks typically fragment into π + π − pairs with similar invariant masses when the dihadron carries a small to moderate momentum fraction z.As z → 1, there is a tendency for the c quark to produce a heavier π + π − pair than the u, d, and s quarks.When more π + π − pairs are produced (smaller z), we expect on average they will have a smaller mass, in which case the mass of the quark flavor becomes less relevant.As less π + π − pairs are produced as one nears threshold (z → 1), then the mass of the quark more directly correlates to the mass of the dihadron.The mass of the b quark being much heavier than the others, it separates itself to yield heavier dihadrons even at small z.The large uncertainties for the gluon expectation values are a reflection of the large uncertainties on D g 1 .The right panel of Fig. 13 displays a clear hierarchy of the average dihadron momentum fractions, with z decreasing from the lightest (u, d, s quarks) to the heaviest (c, b quark).This pattern continues until around m c , where π + π − pairs produced from the c quark start to carry slightly larger momentum fractions that are comparable to dihadrons produced by the s quark.The trend in the plot aligns with the idea that for a given M h , π + π − pairs need to carry smaller momentum fractions if they arise from the fragmentation of a heavier quark.The gluon produces dihadron pairs with the same z (≈ 0.35 − 0.4) independent of the invariant mass.
In Figs. 14 and 15 we show the IFF at the scale µ 2 = 100 GeV 2 .As mentioned in Sec.II A, the sign of H ∢,u 1 cannot be fixed by the experimental data alone, and we have thus forced H ∢,u 1 to be negative to be consistent with the expected sign of h uv 1 .We generally find that the IFF is not constrained by the positivity bound |H ∢,u 1 | < D u 1 , except at large M h and z where its magnitude begins to be limited by the bound.The IFF decreases with increasing z and peaks at M h ≈ 0.8 GeV.

A. Main Results
As discussed in Sec.IV B, we perform three different analyses for the transversity PDFs to assess the impact of experimental datasets and LQCD on the results.In Fig. 16 we show the JAMDiFF (no LQCD) and SIDIS only results for the transversity PDFs and the tensor charges.The two fits agree within errors, with the JAMDiFF (no LQCD) fit preferring a slightly larger u v distribution and the SIDIS only fit preferring a more negative d v .These results translate directly into the tensor charge values, with the SIDIS only result being slightly smaller for δu and more negative for δd.
From the SIDIS only fit one can easily understand the signs of the extracted transversity PDFs.Referring to Fig. 6, one sees that the proton asymmetry from COMPASS is negative.In Eq. ( 6), we showed that the proton asymmetry is dominated by h uv 1 .Since H ∢,u 1 is (chosen to be) negative (see Fig. 14), h uv 1 must be positive to lead to an overall negative asymmetry.A similar argument applies to the HERMES proton asymmetry, except that the negative sign in Eq. ( 5) makes the asymmetry positive.For the deuteron asymmetry from COMPASS one sees that it is largely consistent with zero.From Eq. ( 7) the asymmetry is proportional to the sum of h uv 1 and h dv 1 .Since the proton asymmetry fixes h uv 1 to be positive, the deuteron asymmetry then fixes h dv 1 to be negative so that the sum largely cancels.This qualitative finding nicely agrees with a model-independent QCD analysis for large N c [141].Taking the SIDIS results for h uv 1 and h dv 1 and neglecting antiquarks, one can also understand the sign of the pp asymmetry.The measurements are at mid-rapidity where ŝ ≈ − t ≈ −û.The q ↑ q(q ′ ) → q ↑ q(q ′ ) and q ↑ g → q ↑ g channels give the main contributions and their hard factors in this region are negative (see Appendix A).The numerator of the asymmetry for these channels involves the combinations (cf.Eq. ( 9)) < 0, respectively.Thus, when multiplied by the negative hard factors, one ultimately obtains a positive asymmetry.
In Fig. 17 we compare our results with and without LQCD for the transversity PDFs to those from Radici, Bacchetta [68] (whose analysis did not consider the inclusion of lattice data).We also compare to a version of JAM3D that has been slightly updated from Ref. [54] (see the footnote [148]) that we will refer to as JAM3D * .For the no LQCD results we agree with RB18 within errors, but with a larger h uv 1 in the region 0.04 ≲ x ≲ 0.3.The overall smaller errors on our analysis can be partially attributed to the inclusion of all three SIDIS binnings (x bj , z, M h ), while the RB18 analysis only included the x bj binning.We note that the inclusion of the small-x constraint (Eq.( 23)) and antiquarks in this analysis (neither of which were included in RB18) have no significant impact on the valence transversity distributions in the measured region.For details on the comparison to JAM3D * , see Ref. [108].
Within our analysis, the increase in h uv 1 in the x ≳ 0.3 region when LQCD is included is a consequence of the LQCD results for δu being larger than the result of the fit without LQCD (see Fig. 2 of Ref. [108]).The experimental data provides strong constraints in the 0.005 ≲ x ≲ 0.2 region (see Fig. 1), while the small-x constraint of Eq. ( 23) prevents h uv 1 from becoming significantly larger at x ≲ 0.005.Similarly, the Soffer bound does not allow the transversity PDF to get large at very high x.Thus, in order to increase δu, the best option for the fit is to increase h uv 1 in the x ≈ 0.3 region, although this leads to a slight deterioration in the description of the SIDIS data (see Table II), especially the COMPASS x bj binning.The lack of overlap within errors between the no LQCD and w/ LQCD results will be discussed in detail in Sec.VI B. For h dv 1 , the largest change occurs below x ≈ 0.05, where, after the inclusion of the LQCD data, it now tends to be negative.This is a consequence of the LQCD result for δd being more negative than the result of the fit without LQCD.
In Fig. 17 we also show in the bottom row the JAMDiFF (no LQCD) result for the antiquark transversity distributions and compare to JAM3D * , where both analyses assume h ū 1 = −h d 1 (see Sec. III C).The two results are in agreement.The Soffer bound forces the antiquarks to be very small above x ≳ 0.3.Below that region, they still remain small and consistent with zero.The smaller uncertainties on the JAM3D * result are partly due to the less flexible parameterization used in that analysis.The inclusion of the LQCD measurements causes h ū 1 to become slightly negative at small x, though the result is still compatible with zero.For JAM3D * , h ū 1 is negative up to x ≈ 0.2.While the observation of nonzero antiquark transversity PDFs is interesting, it would be premature to assign any significance to this result.
We now move on to the tensor charges.We refer the reader to Ref. [108] for a detailed discussion of the resulting tensor charges from the no LQCD and with LQCD fits.In Fig. 18 we present a full comparison of the JAMDiFF (w/ LQCD) fit to recent LQCD results [70-73, 75? -80].The phenomenological results of JAMDiFF (w/ LQCD) and JAM3D * (w/ LQCD) are generally on the lower end for g T when considering all LQCD values, but nevertheless are in good agreement with ETMC, NME, Mainz, and LHPC.They are in reasonable agreement with PNDME, RQCD, and JLQCD, and the worst agreement is found with QCDSF/UKQCD/CSSM, PACS, and χQCD.We thus find that our extracted g T is compatible with most LQCD calculations.Our tensor charge results for the JAMDiFF (w/ LQCD), JAMDiFF (no LQCD), and JAMDiFF (SIDIS only) fits are summarized in Table III.definition that now has a number density interpretation [65], we made the first calculations of expectation values for the π + π − dihadron invariant mass (as a function of z) and momentum fraction (as a function of M h ), studying also how these vary with parton flavor.
Furthermore, we have tested the compatibility of different techniques used for determining the tensor charges of the nucleon.By incorporating theory constraints on the transversity PDFs at small x and large x we were able to meaningfully include LQCD results for δu and δd from ETMC [74] and PNDME [71] into our analysis.We found compatibility with this data while maintaining a very good description of the experimental measurements.In addition, our results for δu, δd and the x dependence of transversity match closely to those from the single-hadron TMD/collinear twist-3 analysis of JAM3D [54,148].We have thus demonstrated, for the first time, the universal nature of all available information on the transversity PDFs and tensor charges of the nucleon.We also note that we extracted information about the antiquark transversity PDFs, which had not been done previously for dihadron production.
There are several future directions for extending this analysis.More data is expected from STAR on not only the dihadron azimuthal asymmetry but also the unpolarized cross section.The latter will be important for providing the first experimental constraint on the gluon D 1 (z, M h ).We also anticipate measurements of dihadron multiplicities in SIDIS to offer further experimental quark flavor separation of D 1 (z, M h ).Furthermore, it is mandatory to extend the theoretical framework to next-to-leading order (NLO).Generally, NLO corrections to cross sections can be large, and they will certainly modify our quantitative results.However, we expect that our main qualitative findings about the transversity PDFs and the tensor charges will not be affected, as the relevant observables are asymmetries for which typically significant cancellations of higher-order corrections occur.Nevertheless, a NLO analysis is needed to provide a definitive answer to that question.Overall, future measurements at STAR, Jefferson Lab, and the Electron-Ion Collider can provide crucial information to help reduce uncertainties on the transversity PDFs and tensor charges [153,154].In particular, we identified the large-x (in the case of SIDIS) and forward rapidity (in the case of pp) regions as key to further test the compatibility between LQCD results for the tensor charges and experiment.Measurements in the small-x region will also help test the small-x asymptotic behavior of transversity imposed in our analysis and give further insight on the antiquark transversity PDFs.Moreover, dihadron, single-hadron TMD/collinear twist-3 observables, and LQCD (tensor charge and x-dependent transversity computations [77,155]) should eventually be fit simultaneously in a "universal" analysis.û ↔ t relative to Ref. [122].This leads to the following for the partonic cross sections involving quark flavors q, q ′ : d∆σ q ↑ q→q ↑ q d t = 8πα The hard factors involving antiquark fragmentation are identical to the hard factors of the corresponding chargeconjugated partonic processes.We note these are also consistent with the partonic cross sections associated with the "derivative term" of the fragmentation part of A N in p ↑ p → h X [93].
Appendix B: Symmetry of π + π − DiFFs We provide details on deriving the symmetry relations for π + π − pairs that allow us to reduce the number of DiFFs to be fitted, which is an important aspect of the analysis.For this discussion, it is necessary to introduce the notation D h1h2 , where h 1 h 2 indicates the dihadron pair.First, under the interchange of h 1 and h 2 one has the relations for all quarks q.The sign change for H ∢ 1 occurs due to the fact that interchanging the hadrons switches the sign of ⃗ R T , which appears in the prefactor of H ∢ 1 in its correlator definition, but the correlator itself is independent of the ordering of h 1 and h 2 [56,65].From charge-conjugation symmetry, one also has the relations for all quarks q.Combining Eqs.(B1) and (B2) to eliminate the π − π + DiFFs, one has From isospin symmetry (which holds to sufficient accuracy for this discussion), one also has the relations for the light quarks

FIG. 1 .
FIG. 1. Kinematic coverage of the datasets included in this analysis.The top panel shows the e + e − annihilation and SIDIS data as a function of z and M h .The bottom panel shows the SIDIS and pp data as a function of x and µ 2 .The variable x represents x bj for SIDIS and xa for pp collisions, while the scale µ 2 represents Q 2 for SIDIS and P 2 hT for pp collisions.For STAR, the solid points are x min a and the light shaded region is to indicate that the xa-integration extends up to xa = 1.
7 and 10 for the STAR √ s = 200 GeV data and in Figs.8-10 for the √ s = 500 GeV data.The fitted normalizations N e (see Eq. (

FIG. 16 .
FIG. 16.Left panel: Transversity PDFs xh uv 1 and xh dv 1 plotted as a function of x at the scale µ 2 = 4 GeV 2 .We show results for the JAMDiFF (no LQCD) fit (red) and SIDIS only fit (violet) with 1σ uncertainty.The Soffer bound is indicated by the dashed black lines.Right panel: The tensor charges δu and δd at the scale µ 2 = 4 GeV 2 for the JAMDiFF (no LQCD) and SIDIS only fits with 1σ uncertainty.

FIG. 21 .
FIG. 21.PYTHIA-generated data at √ s = 30.73GeV.The strange (green points), charm (orange points), and bottom (pink points) cross section ratios are plotted as a function of M h and compared to the mean JAMDiFF (w/ LQCD) results with 1σ uncertainty bands (colored bands).The different panels show different bins of z.

TABLE I .
Data from experiment and LQCD used in this analysis for the DiFFs D1 and H ∢ 1 and the transversity PDF h1.

TABLE II .
Summary of χ 2 red values for the different fit configurations defined in Sec.IV C.

TABLE IV .
Summary of the different PYTHIA 6 tunes used to generate the PYTHIA data.The "PYTHIA 6 default" tune uses the corresponding default parameters shown above.The "PYTHIA 6 ALEPH" and "PYTHIA 6 LEP/Tevatron" tunes modify those parameters as shown in their respective columns (if the row is blank then the default value is used).For the "PYTHIA 8" tune (not shown in the table) we use all default parameters.

TABLE V .
Summary of χ 2 red values from the JAMDiFF (w/ LQCD) fit for the PYTHIA datasets where the observable is σ q /σ tot .