Constrained fit of the valence transversity distributions from dihadron production

We present a constrained analysis of the valence transversity Parton Distribution Functions from dihadron production in semi-inclusive DIS. While usual extractions of the transversity distributions rely explicitly on the fulfilment of the Soffer bounds, the present analysis releases that implicit restriction to implement further explicit constraints through the Lagrange multipliers method. The results are quantitatively comparable to previous analyses in the kinematical range of data ; the qualitative impact of the chosen fitting strategy translates into an increased flexibility in the functional form.


Introduction
Parton Distribution Functions (PDF) describe the structure of hadrons in terms of their constituents. At leading-twist, this structure is defined by the distributions corresponding to the vector, axial-vector and tensor Dirac bilinears. The first two, the unpolarized and the helicity Parton Distribution Functions, while under constant study and improvement, are relatively well-determined. In the last decade, the third leading-twist PDF, that corresponds to the tensor bilinear [1,2], has been explored in phenomenology. The transversity PDF describes the distribution of the transverse polarization of quarks as related to that of their parent hadron. As such, it is a chiral-odd object who is paired, in physical processes, to another chiral-odd object.
The best proposals for an optimizing partner list: the single-hadron Fragmentation Function (FF) depending on the fraction of longitudinal momentum of the fragmenting quark carried by the outgoing hadron as well as on the trace of the intrinsic transverse momentum of that active quark ; the second candidate is the dihadron FF which involves a fragmentation process into a pair of hadrons defining a production plane whose orientation w.r.t. the orientation of the transverse polarization of the active quark generates modulations at the cross-section level. Thanks to the data for semi-inclusive pion production in electron-positron annihilation from the B factories, both proposals have been pursued down to the level of either global fits and single-process fits.
Dihadron Fragmentation Functions (DiFFs) were proposed by Jaffe et al. [3] to access the transversity PDF. A history of successes in the most recent years have led to a bright phenomenology, as demonstrate the recent determinations of the transversity PDF [4,5] -that is, the collinear PDF as opposed to the Transverse Momentum Distribution [6][7][8]-, effort to which this manuscript is summing.
Yet another possible way to the transverse polarization of quarks is found through Generalized Parton Distributions relevant in Deeply Virtual Meson Production [9].
Even though the experimental statistics for processes involving transversities do not reach that of either the unpolarized or the longitudinally polarized PDFs, it would be timely to adopt new developments in methodology, as has been done by the PDF community, for the determination of the transversity PDF. With this paper we will take a step in that direction. For that, it is important to review the fundamental properties of the transversity PDF.
. PDFs are subject to first principle constraints. In the usual helicity basis representation for Dirac particles, the transversity distributions, h 1 (x), do not correspond to probability densities -unlike both the unpolarized, f 1 (x), and longitudinally polarized distributions, g 1 (x). As a consequence, the following inequality was proposed by Soffer [10] |h q for each flavor of quark and anti-quarks separately. The Soffer bound is conserved under QCD evolution towards higher scales [11,12]. The bound in Eq. (1) has been treated as first principle in all extractions of transversities until now. Together with the support -x ∈ [0, 1]it constitutes the main set of first principles to which the distributions must obey. As a result of the fundamental characteristic of the transversity PDF to describe the internal structure of hadrons in a given spin configuration, its first Mellin moment -the moment that, for conserved currents, embodies the fundamental charges-provides interesting information for processes involving hadrons that are also relevant for searches of New Physics [13][14][15]. In view of the phenomenological results [5,7] and conflicting agreements with the latest lattice QCD determinations of the tensor charge [16][17][18], we analyze the statistical probability of the -expression of the-Soffer bound given the semi-inclusive DIS data. Subsequently this statistical significance will serve as a weight for the fit of the transversities for which no first principle constraint is assumed other than the support in Bjorken x, implemented through the parameterization. The methodology is complemented by the choice of the mentioned functional form such to optimize the assessment of an uncertainty that reflects as objectively as possible the error coming from the preceeding steps. The behavior at large values of the momentum fraction will be guided through a further iteration of the minimization imposing inequality constraints that simulate a smooth fall-of for x → 1 through the Lagrange multipliers method. This paper is organized as follows. In Section 2, we overview the formalism related to dihadron production in semi-inclusive DIS. Our methodology is described in Section 3. In Section 4 we discuss our results and then conclusions and possible extensions of this work are drawn in Section 5.

Formalism
In this work we consider the access to the transversity PDF that is made through a Single Spin Asymmetry in the dihadron semi-inclusive production off a proton/deuteron target process, l(k) + N (P ) → l(k ) + (H 1 (P 1 ) + H 2 (P 2 )) + X .
The (unpolarized) lepton beam is denoted by l with 4-momentum k and the outgoing lepton k . The transversely polarized target N will be a proton or deuteron that is accessed through nuclear targets, i.e. hydrogen at HERMES [19] and 6 LiD (deuteron) and NH 3 (proton) at COMPASS [20]. The produced hadrons are labelled H 1 and H 2 , their invariant mass squared is P 2 h = M 2 h , with P h total momentum of the pair. Their relative momentum isdefined, R = (P 1 − P 2 )/2. The momentum transferred to the nucleon target is q = k − k . The usual DIS invariants are defined, x = Q 2 /(2P · q), y = P · q/P · k and z = z 1 + z 2 = P · P h /(P · q). We refer to Ref. [4] and references therein for further details.
In this process, the chiral-odd PDF couples to a Fragmentation Function whose chiral behavior is understood as the transverse polarization of the fragmenting quark related to the orbital angular motion of the transverse component of the pair relative transverse momentum R T via the azimuthal angle φ R of the hadron pair plane [21]. The transverse momentum R T can be expressed in terms of the invariant mass of the pair and the polar angle θ of the one pion relative to P h in their center-of-mass.
In the limit M 2 h << Q 2 , to leading-order and after having selected the dominant contribution through Partial Wave Expansion in the polar angle θ [22], the relevant -Single-Spin-asymmetry reads where e q is the fractional charge of a parton with flavor q, A(y) = 1 − y + y 2 /2, B(y) = 1 − y.
The D q 1 is the DiFF describing the hadronization of an unpolarized parton with flavor q into an unpolarized hadron pair. The H q 1 ≡ H q 1,sp is a chiral-odd DiFF describing the correlation between the transverse polarization of the fragmenting parton with flavor q and the azimuthal orientation of the plane containing the momenta of the detected hadron pair. Both have been determined from Belle's data and corresponding PYTHIA-generated multiplicities [4,23].
The 3D asymmetry can be projected onto each of the relevant variable by integrating the two others over their corresponding kinematical intervals. As the PDF only depends on x, we here only use the x-projected part of the triptic. The relevant combinations for the proton target are [24] x h p and for the deuteron target are where h qv 1 is the valence combination, i.e. h q 1 − hq 1 , and f q+q 1 ≡ f q 1 + fq 1 . Unless explicitly stated, we will use the MSTW08LO PDF set [25] to evaluate the unpolarized PDFs. The quantities n q (Q 2 ) and n ↑ u (Q 2 ) are, respectively, the integrals over the respective range of (z, M h ) at a given value Q 2 of the unpolarized DiFF for a flavor q and the chiral-odd DiFF for u. We will use the bootstrapped version of the DiFFs [4]. In the present manuscript, we will not take into account the Q 2 dependence of the transversity PDF but rather consider its determination as given at an average scale Q 2 = 5GeV 2 -that corresponds to the average Q 2 value of the data. The dependence on the hard scale is kept in all the quantities on the r.h.s of Eqs. (4-5) as we do not expect cancellations between numerator and denominator, in that sense our approximation is slightly different from that adopted in Ref. [8]. Consequences of this approximation will be discussed in Section 4.
In this paper, we use the relevant Single Spin Asymmetry detected at HERMES [19] and COMPASS for identified pairs [26,27]. Those data sets coincide with those used in Ref. [4]. The kinematical ranges are 1.2 < Q 2 < 31.5 GeV 2 , 0.0064 < x < 0.2871 -with increasing values of Q 2 for increasing values of x-and the DiFF variables are 2m π < M h 1.3GeV and 0.2 z < 0.9.

Methodology
The determination of the valence transversities is obtained in three main steps. Following the hypothesis that first principles constraints can be implemented through the fitting procedure, we first assess their goodness given the data using a Bayesian approach. This first step allows a reweighting of each data point, taking into account implicitly, e.g., the Soffer bound. Then the fulfilment of other theoretically justified behaviors is achieved through a constrained fit using the Lagrange multipliers method. Finally, the constrained fit is repeated N times using the bootstrap technique. In particular, we will use the exact same bootstrapped dihadron fragmentation functions as in Ref. [4].

Probability of the Soffer Bound
The practical implementation of the Soffer bound relies on the other two leading-twist PDFs, to be precise on the fits of the unpolarized and helicity PDFs. An uncertainty is attached to each PDF value at the corresponding (x, Q 2 )-bin that should combine the proper uncertainty -output of the fitting procedure-when provided and an error coming from the various choices for physical parameters and hypothesis, e.g. the value of α s (M 2 Z ) and the order in perturbation theory.
In the standard parametrizations of h 1 (x), the conservation of the bound is ensured by construction at the level of the functional form at Q 0 and is conservatively fulfilled under QCD evolution [12]. Hence, a theoretical error could analogously be attributed to the Soffer bound (SB) in the first place, introducing an extra flexibility on the parametrization.
On the other hand, the confrontation of the SB with the combination of transversities, Eqs. (4)(5), also depends on the approximations used in the extraction of the DiFFs. The mentioned theoretical error should ideally comprise both sources of uncertainty. Thus, working within the hypothesis that all the above mentioned is properly assessed, we would like to estimate the goodness of the Soffer bound -expressed in terms of PDF through fit outputs-given the transversity combinations extracted from HERMES and COMPASS data.
For that purpose we use a Bayesian approach to evaluate the probability of the PDF lying outside the SB. We first map each experimental point, which has a continuous -Gaussiandistribution, into an in-out case. Suppose that we know the limits of the bound exactly, i.e., we can identified unequivocally whether the data are inside or outside the bound: θ j represents that probability, where A j is the true value of the transversity combinations. The region is the area comprised by the inequality of Eq. (1) evaluated in each corresponding kinematical bin using MSTW08 at LO for f 1 (x, Q 2 ) and JAM15 [28], at NLO, for g 1 (x, Q 2 ). We assume an error from going from NLO to LO for the helicity PDF based on the NLO-LO difference for unpolarized PDFs. We define an hyperparameter t corresponding to the probability that the data value lies outside Figure 1: Distribution of the hyperparameter t as given by Eq. (9). The vertical blue line corresponds to t bay , the blue shaded area to the 68% C.L. evaluated from the upper/lower bound of the integral of F (t). The orange shaded area corresponds to the 95% C.L.. The y-axis has arbitrary units.
namely a prior distribution for θ j given t. The true values for the transversity combinations can only be inferred through the actual data, such that The desired probability is expressed as with N = 1 0 F (t) dt, the norm, and π(t) is the prior for t that is chosen to be flat and the other probabilities are defined above. Evaluating F (t) assuming the A i are gaussianly distributed, we find a distribution with a mode at t = 0 and a central value for t, and the 68% confidence interval is [0.009, 0.09]. It is illustrated in Fig. 1. How should we interpret this result? It is highly unlikely that the bound is incorrect by more than 10%. There exists a window through which the agreement of the implemented SB and the transversity combinations is poorer. Since the SB is considered as first principle, this result is translated here into a relaxation of the expression of the bound through a statistical reweighting bin-by-bin in the objective function -i.e. here the χ 2 . The Bayesian weight is obtained as follows, In Fig. 2, the Bayesian weight is represented for each bin in the bottom plots. It can be appreciated that the proton combination is almost not affected by this procedure while the Figure 2: Combinations of transversities for the proton (left) and the deuteron (right) as extracted in Ref. [4]. The bottom plots show the statistical weight for each bin as evaluated through the t-distribution.
lowest and highest x-bins of the deuteron combination happen to statistically lie inside that window of poor agreement, as expected from previous publications [4].
This result is included in the fitting procedure through the bootstrap method. As explained in, e.g., Ref. [4], N = 200 replicas of the extracted combinations x h p/D 1 j , with j = 1, 22 the data label, of Eqs. (4)(5) are generated randomly within the data errorbar. As reweighting the overall chi-square function can be understood as modifying the error by σ 2 j → σ 2 j /w j , the replicas are generated within the corresponding extended gaussian. The 200 thus generated data sets are shared among four functional forms, minimizing each one for 50 replicas. The objective function can be written, for each replica r, as where {p I } is the set of free parameters to be determined. In the next Section, the function h 1 theo will be extensively described.
The same exercise has been repeated with CT09MC1 [29] in both Eqs. (4)(5) and the expression of the SB. The lightest weights for the deuteron combination slightly change, repercussions will be discussed in Section 4. The overall conclusions of this section are not affected by the choice of LO unpolarized PDF.

Parameterization and Constraints
Now that the main first principle based constraint has been implicitly included through a Bayesian reweighting -that is effectively implemented at the stage of bootstrapping the data-, the parameterization can take a more adaptable form: the data alone lead the determination of the free parameters of the chosen form. The trend of the data, confined in a small kinematical range of x, consists in a growth with increasing values of Bjorken x for x h u V 1 . Such a behavior is not appealing knowing that the support of the parton distributions should be restricted to x ≤ 1. The desired support can be ensured through a judicious choice of the functional form and/or a constrained fit, see e.g. [30]. The behavior of x h d V 1 is affected by bigger errors w.r.t. the up. In the same fashion, a flexible parameterization can be guided to decrease smoothly with increasing x, as expected by the underlying QCD evolution ocurring at higher x values, see e.g. [31], and power-law fall-off ruling the two other leading-twist PDFs. This observation indirectly relates to the Soffer bound as expressed in terms of PDF fits.
The minimization itself contains two steps. First, the objective function -the χ 2 function expressed in terms of Eqs. (4)(5) and the respective propagated error-is minimized via a nonlinear least-square implemented through scipy.optimize.curve_fit in Python, determining the set of q best fit parameters, {p}. The definition of the χ 2 ({p}) requires an expression for a functional form. The first considerations while determining the parameterization is to ensure integrability and support at x = 0, where the exponent has been chosen based on previous outputs and the expected small-x behavior [32] to be slightly modified by the polynomial in x. The latter, P q n (g(x)), of order n, can be as flexible as the data allow for. Namely a balance between flexibility and loss of degrees of freedom has to be found. We choose to express P n in terms of Bernstein polynomials as has been done in Ref. [33]. They are defined as with n k the binomial coefficients. Those polynomials have the advantage of selecting particular regions in g(x) and as such can be employed to emphasize particularly relevant kinematical regions. That is, for semi-inclusive DIS, the low-x region becomes relevant w.r.t. the valence and large-x region. A rescaling of the variable will help the parameterization adjusting the data, we choose g(x) = x 0.3 . A desirable feature of the polynomials in Eq. (13) is a statistically representative error outside the data range, yet in agreement with the first principle constraints at hand. In order to span the Bjorken variable range as significantly as possible, we use 4 different degrees for the polynomial P n -thus distinguishing the 4 different functional forms by their order n-and optimize the number of Bernstein polynomials by trial and error. In Fig. 3, we show the polynomials that we have adopted, respectively, for the up and the down parameterization, expressed as with i = 1, · · · , 4, n = {10, 20, 30, 40} and where the p q i,k are the free parameters which number varies for each i. The set of values {κ q,i } has length n κ u,i = {2, 3, 3, 3} and n κ d,i = {4, 3, 4, 3} and their values have been chosen to flexibly span the relevant region consistently with the data, as can be seen in Fig. 3.
In the following step, we will clarify our choice for a limited coverage of the functional form for values of x 0.6 for the up and x 0.5 for the down. As mentioned above, we need to insure a smooth fall-off of the transversity in the limit x → 1 that cannot be achieved exclusively from the choice of functional form. In most cases, a second step will be required to constrain the functional form in an allowed region. When the objective function is subject to m constraints of the form C j ({p }) = 0, the later are imposed through the Lagrangian to which a stationary point of L is found minimizing with respect to the parameters {p } and the Lagrange parameters {λ}.
Once the convergence of the first step guaranteed and given the linearity of our functional form in terms of the parameters, it is sufficient to define our new objective function as follows The vector of parameters p q I i,k , of length n κ u,i + n κ d,i , corresponds to {p I }, the set of best fit parameters obtained through the main minimization, and the covariance matrix V also comes from step I. The chi-square function depends on the new set of best fit parameters, {p II }, which consists in the set made of p q II i,k . In previous -unpolarized and longitudinally polarized-PDF determinations, the method of the Lagrange multipliers has been made popular for error estimation [34]. In the present approach, this method is used to impose limits on the fit parameters. In particular, we use the more general inequality constraints through scipy.optimize.minimize in Python, which is based on the Lagrange multipliers method above described [35]. We guide the large-x behavior of the down parameterization only using the following N c = 4 constraints

Results
Both steps of the minimization procedure are carried out in Python. The convergence for both the main χ 2 minimization Eq. (12) and the optional Lagrange multipliers method are fast. Twenty-five percents of the replicas converge inside the established bound at the first minimization. The 75% left is hence constrained by the Lagrange multipliers method. For the higher-order polynomials, less than 15 additional iterations are necessary while the two lowerorder polynomials require about 30-40 extra iterations. Although the constraints are imposed for the down distribution only, their fulfilment affects directly the parameterization of the up. The balance between u V and d V is clearly passed on through the objective function. Notice that not all 4 constraints are active for each replica r.
In Fig. 4, we show the envelopes for the obtained valence transversities at 68% CL, respectively, for the up and the down distributions. Following the color code introduced in Fig. 3, it can be observed that the higher order polynomials -n = 30, 40 in yellow and green-allow a larger error bar at smaller values of x. On the other hand, the lower order polynomialsn = 10, 20 in red and purple-are confined in the mid-to large-x region. All 4 equally contribute  Figure 6: Combinations of transversities for the proton (left) and the deuteron (right) compared to the global 68% in blue (with full blue line contour) and 95% in orange (with dashed contours) for the present fit compared to the extracted data [4]. The dark blue square come from HERMES data while the dark orange dots are extracted from COMPASS data.
at valence values of the Bjorken variable. The final envelope of xh q V 1 (x, 5GeV 2 ) is built using the four versions of the functional form together. It is shown at 68% and 95% CL in Fig. 5. For consistency, we also show the Soffer bound at NLO using the helicity PDF of JAM15 [28] and CT18 [36]. Were there central values for the transversity PDFs, they would be enclosed inside the Soffer bound. In that sense, our result is similar to the first Hessian approach of Ref. [24].
The small-x behavior predicted in Ref. [32] is remarkably fulfilled for the up transversity PDF up to, at least, x = 0.1 ; no clear conclusion can be drawn for the down above x ∼ 0.03.
In Fig. 6, the 68% and 95% CL proton and deuteron combinations are depicted and compared to the point-by-point extraction of Eqs. (4)(5). The final χ 2 /d.o.f. is evaluated against those extracted data. The number of degrees of freedom here is ( j w j )×22+N c −(n κu +n κ d ). The corresponding stacked histogram is shown in Fig. 7.
We have repeated the fitting procedure using the CT09MC1 set for the unpolarized PDF [29], it is a LO set fitted on real data and NLO pseudo data with α s at 1-loop. There is no qualitative difference in the final transversity PDF. We notice a slightly smaller average value for the isovector tensor charge, quantity that will be defined here below.
Compared to previous extractions, we might comment that our result for the down transversity presents a smaller errorband in the data region -as we build our functional form on that purpose. It is achieved at the expense of a slightly wider band for the up distribution in that region. The errorband increases for small-and large-x values. In this respect, our down distribution differs from that of Refs. [4,5] and the less flexible versions of the parameterization of Ref. [24], all of them being Dihadron-based extractions. The same comment is in order for comparisons with single-hadron semi-inclusive DIS [6,7].
We also compare our results to the lattice QCD evaluation of the isovector transversity PDF [37]. The order of magnitude as well as mid-to large-x behavior are in a reasonable agreement. The lower-x region of our result is more structured, as dictated by the data. However, our analysis has been carried out without QCD evolution for the transversity PDF. It is to be understood as given at an average value of Q 2 = 5GeV 2 . On the other hand, the lattice result corresponds to a renormalization scale of 4GeV 2 . The appearent inadequacy at x ∼ 0.001 − 0.002 could be a shortcoming of our approximation. To that regard we have further performed our analysis considering QCD evolution only for the Fragmentation Functions. In  Figure 8: Comparison with the isovector combination of the transversity PDF from lattice QCD [37], in gray. The cyan band represent our 95% CL transversity combination.
fixing the unpolarized PDF's scale to 5GeV 2 in Eqs. (4)(5), we observe that the down distribution is narrower in the valence region and the up errorband is wider at small values of x. The chisquare increases of about 15%. A wider distribution for g T with an higher mean value is found. However, the comparison with the isovector transversity from lattice QCD is not substantially improved. Finally, we notice that the PDF+Lattice result for the isovector combination [8] is systematically higher than ours for the whole support in x, except for x → 0.
We next evaluate the first Mellin moment of the transversity PDF to get the tensor charge The distributions correspond best to skewed distributions and the obtained values are In particular, we are interested in the isovector combination, g T (Q 2 ), as it is of particular interest for, e.g., beta decay observables [13,14]. We show the corresponding stacked histogram in Fig. 9. The Gaussian distribution showed on the r.h.s. of Fig. 9 is given by at 1σ ; the 2σ error is 0.42. The truncated tensor charge is found to be g T ( 5GeV 2 ) = 0.48 +0.13 −0.11 , corresponding to a skewed distribution.
The trend observed above leads to discernible consequences here: the more flexible the parameterization, the wider the distribution of the tensor charge. Our result is compatible with lattice QCD evaluations -the latest results are g T = 0.926(32) [16], g T = 0.972(41) [17] and g T = 0.965(38)(+13 − 41) [18]. It is at the same time wide enough to encompass all of the other determinations. In a region where scarce to no data are avaible, the choice of functional form generates a crucial uncertainty. The role of the small-x region, where no strong expressions of first principle bounds nor data are available, matters. The addition of proton-proton collision data [38], used in the first global fit of the transversity [5], does unfortunately not increase the coverage at lower x values -but it certainly does in Q 2 . In that sense, data from an Electron Ion Collider would be extremly helpful, be it for a qualitative improvement over a quantitative reduction of the uncertainty. Would the tensor charge become relevant in view of search for New Physics, the appropriate observables would not necessarily linearly depend on the uncertainty of the tensor charge, as shown in the case of beta decay in Ref. [14].

Conclusions
The transversity Parton Distribution Function is the least known of the three leading-twist PDFs. Its phenomenological determination has been made possible thanks to independent data for their partners in semi-inclusive DIS, the chiral-odd fragmentation functions. In particular, the formalism developed around dihadron Semi-Inclusive DIS allows for a collinear extraction of combinations of valence transversity PDFs [21]. It requires the knowledge on the Dihadron Fragmentation Functions that have been studied in Ref. [4,23]. Since its first point-by-point extraction from Dihadron semi-inclusive DIS [39], enormous efforts towards the extraction of the transversity PDF accessible in processes involving fragmentation to a pion-pair [5] have been made.
In this paper we have presented a fit of the valence transversity PDFs from dihadron semiinclusive DIS data without explicitly involving the Soffer bound. That important constraint on the transversity distribution is implemented through a Bayesian reweighting that reflects its goodness in outlying kinematical bins. It is similar to assigning a theoretical uncertainty or a prior to the expression of the bound. Then, the desired behavior at x → 0 and x → 1 is found by interplay between functional form for both the up and down PDFs -setting the support in x-and constraints for the down distribution -dictating the large-x behavior-, imposed through the Lagrange multipliers method. While the Soffer bound was only implicitly considered, a guided fall-off of the PDFs for larger x values indirectly guarantees that the down distribution stays inside an envelope that was designed based on a plausible power-law fall-off.
The obtained 68% and 95% CL envelopes for the valence transversities fulfill the expected small-x behavior, the up distribution fully respects the Soffer bound and the edges of the down envelope reflect the relaxation of the bound combined with the lack of data in the large-x region. Our results globally show wider errorbands outside the data range w.r.t the errorband inside that range than found in previous extractions. That trend, allowed by the -absence of-data, translates into a wide distribution of the tensor charge yet with a more comfortable extrapolation in the whole support. We believe that this excercise supports the idea that the choice of functional form contributes to the global error of the PDF determination, especially in kinematical regions with limited data coverage.
The effect of QCD evolution will become important when including data from protonproton collision from RHIC. A natural extension of this work could be the inclusion of the corresponding DGLAP routine and the aforementioned data. More semi-inclusive DIS data in the valence region are expected from CLAS12 and SoLID at JLab ; data from the future Electron Ion Collider would improve the uncertainty in the low-x region.