1 Introduction

Parton Distribution Functions (PDFs) describe the structure of hadrons in terms of their constituents. At leading twist, the hadronic 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, the one corresponding to the tensor bilinear [1,2,3,4], 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. Transversity is a chiral-odd object that is paired, in physical processes, to another chiral-odd object.

One such candidate object is the single-hadron Transverse Momentum Dependent Fragmentation Function (TMD 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 (DiFF), measured in a fragmentation process into a pair of hadrons whose relative momenta tag the transverse polarization of the active quark via modulations at the cross-section level. Thanks to the data on semi-inclusive pion production in electron–positron annihilation from the B factories, both types of FFs have been measured and used to assess transversity PDFs in either global fits and single-process fits.

Dihadron Fragmentation Functions were proposed by Jaffe et al. [5] to access the transversity PDF and successfully implemented in phenomenological studies [6], as demonstrated by the recent determinations of the transversity PDF [7,8,9] – that is, the collinear PDF discussed in this article, as opposed to the Transverse Momentum Dependent (TMD) transversity distributions that have been recently studied in Refs. [10,11,12].

Yet another possible way to access the transverse polarization of quarks is found through Generalized Parton Distributions relevant for Deeply Virtual Meson Production [13, 14].

The analyses for the transversity are based on reduced experimental inputs with respect to the data sets used for the other two leading-twist PDFs. All three are subject to analogous first principles, e.g. positivity constraints and support requirement. In this paper, we propose to broaden the methodology to incorporate the assessment of first principles directly in the fitting procedure. This methodology palliates the low experimental accuracy of processes involving transversity PDFs.

To leading-order, the unpolarized, \(f_1(x)\), and longitudinally polarized distributions, \(g_1(x)\), can be interpreted as probability densities in terms of the helicity representation for Dirac particles – basis of the positivity bound \(|g_1^q(x)|<f_1^q(x)\). The interpretation of transversity distribution, \(h_1(x)\), must be done in the transversity basis – leading to the bound \(|h_1^q(x)|<f_1^q(x)\). Another positivity constraint for \(h_1(x)\) has been stated by Soffer [15]

$$\begin{aligned} |h_1^q(x)|\le \frac{1}{2}\left( f_1^q(x)+g_1^q(x)\right) , \end{aligned}$$
(1)

for each flavor of quark and anti-quarks separately. If the Soffer bound is realized at a low factorization scale Q, it is preserved under QCD evolution towards higher scales [16, 17]. The bound in Eq. (1) has been treated as a first-principle constraint in all extractions of transversities until now. Together with the support \(x\in [0,1]\) and the corresponding end-point behavior, they constitute the basic properties that any parametrization of the transversity PDF must satisfy. The implementation of positivity constraints generates a hierarchical dependence on the unpolarized PDF for the helicity PDF and both the unpolarized and helicity PDFs for the transversity.

As a result of the fundamental characteristic of the PDFs to describe the internal structure of hadrons in a given spin configuration, their first Mellin moments embody the fundamental charges when originated from conserved currents. The first moment of the transversity distribution is known as the tensor charge, in analogy to the vector and axial charge. The tensor charge provides interesting information for processes involving hadrons that are also relevant for searches of New Physics [18,19,20,21].

On the one hand, recent analyses of the transversity highlight a tension between the data and the expression of the Soffer bound in some kinematical regions [7, 9, 11]. On the other hand, the tensor charge for the valence and isovector flavor combinations as determined through phenomenology and lattice [22,23,24,25] differ, with the largest discrepancy coming from the u-valence contribution. In light of the phenomenological conflicts with the Soffer bound and the latest lattice QCD determinations of the tensor charge, we analyze the level of probability of the Soffer bound given the semi-inclusive DIS data. Subsequently the derived significance of the bound will be reflected in the objective function. The Soffer bound only appears in that step, hence leveraging the choice of the parametrization in the minimization procedure. The parametrization will be chosen, first, to fulfill the support in Bjorken x and, second, so as to optimize the assessment of the uncertainty that reflects as objectively as possible the error coming from both the experimental data and the theory, i.e. the implementation of the positivity bounds here. The methodology includes a further iteration of the minimization imposing inequality constraints to ensure a smooth behavior for \(x\rightarrow 1\) through the method of Lagrange multipliers. We carry out the new analysis based on the point-by-point extraction for proton and deuteron SIDIS data of Ref. [8]. There have been no recent improvements from the theory side or new determinations of DiFFs that would justify a complete reanalysis of the two-hadron SIDIS asymmetry.

This paper is organized as follows. Our methodology is described in Sect. 2. In Sect. 3 we discuss our results –with a particular emphasis on the tensor charge, and then conclusions and possible extensions of this work are drawn in Sect. 4. The state-of-the-art formalism related to dihadron production in semi-inclusive DIS is overviewed in Appendix A.

2 Methodology

The determination of the valence transversities is obtained in three main steps. We work within the hypothesis that first principles constraints can be implemented through the fitting procedure to optimize the collective information of the \(N_d=22\) data points and the first principles. We first assess the probability that the Soffer bound is valid using a Bayesian approach. This first step allows a reweighting of each data point, taking into account the Soffer bound implicitly. Then the fulfillment of other theoretically justified behaviors is achieved through a constrained fit using the method of Lagrange multipliers. Finally, the constrained fit is repeated N times using the bootstrap technique as an integral part of the error treatment [8].

2.1 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. A PDF parametrization should come with an uncertainty that combines the statistical uncertainty obtained in the fit and an error coming from the various choices for physical parameters and hypothesis, e.g. the value of \(\alpha _s(M_Z^2)\) or the order in perturbation theory.

In the standard parameterizations of \(h_1(x)\), the validity of the bound is ensured by construction at the level of the functional form at \(Q_0\) and is preserved under QCD evolution [17]. 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 extracted from proton and deuteron targets, Eqs. (A.3)–(A.4),Footnote 1 also depends on the approximations used in the extraction of the DiFFs. The mentioned theoretical error should ideally comprise both sources of uncertainty. Thus, we would like to estimate the goodness of the Soffer bound – expressed in terms of PDF parameterizations – given the transversity combinations extracted from HERMES [26] and COMPASS [27, 28] data and implicitly include the bound as a source of uncertainty in the fitting procedure.

For that purpose we evaluate the probability of the transversity PDF lying outside the SB. We first map each experimental point, which has a Gaussian distribution, 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: \(\theta _{j}\) represents that probability,

$$\begin{aligned} p\left( \theta _{j}|A_{j}\right) ={\left\{ \begin{array}{ll} \delta \left( \theta _{j}-1\right) &{} \text{ for }\quad A_{j}\in \text{ region }\\ \delta \left( \theta _{j}\right) &{} \text{ for }\quad A_{j}\notin \text{ region } \end{array}\right. }, \end{aligned}$$
(2)

where \(A_{j}\) is the true value of the transversity combinations of Eqs. (A.3A.4) evaluated for \({j=1,\dots ,N_d}\). The region is the area comprised by the inequality of Eq. (1) evaluated in each corresponding kinematical bin using MSTW08LO at LO [29] for \(f_1(x_j, Q_j^2)\) and JAM15 [30], at NLO, for \(g_1(x_j, Q_j^2)\).

We define a hyperparameter t corresponding to the probability that the data value lies outside the bound,

$$\begin{aligned} p\left( \theta _{j}|t\right) =\left( 1-t\right) \delta \left( \theta _{j}-1\right) +t\delta \left( \theta _{j}\right) , \end{aligned}$$
(3)

namely a prior distribution for \(\theta _{j}\) given t. The true values for the transversity combinations can only be inferred through the actual data, such that

$$\begin{aligned} p\left( \theta _{j}|\text{ data }\right) =\int p\left( \theta _{j}|A_{j}\right) p\left( A_{j}| \text{ data }\right) dA_{j}, \end{aligned}$$
(4)

where the probability of \(A_{j}\) given the data evaluates the distance of the data point j to the bound, considering both the uncertainty on the data and on the bound. The distributions are assumed to be Gaussians centered on the experimental value, for the data, or as explained after Eq. (2), for the bound, with the corresponding standard deviation. We assume an error from going from NLO to LO for the helicity PDF based on the NLO-LO difference for unpolarized PDFs.

The final desired probability is expressed, using Bayes theorem, as

$$\begin{aligned} p\left( t|\text{ data }\right)= & {} N\, \pi (t)\int \, \displaystyle \prod _{j=1}^{N_d}p\left( \theta _{j}|t\right) \, p\left( \theta _{j}|\text{ data }\right) \, d\theta _{j}\; ,\nonumber \\\equiv & {} F(t), \end{aligned}$$
(5)

\(N=\int _0^1\, F(t) \, dt\) is the norm, \(\pi (t)\) is the prior for t that is chosen to be flat, and the probabilities are defined above. Evaluating F(t), we find a distribution with a mode at \(t=0\) and a central value for t,

$$\begin{aligned} {\bar{t}}= & {} \int _0^1\, t\,F(t) \, dt=0.049\pm 0.040, \end{aligned}$$
(6)

for the \(68\%\) confidence interval. It is illustrated in Fig. 1.

Fig. 1
figure 1

Distribution of the hyperparameter t as given by Eq. (5). The vertical blue line corresponds to the central value \({\bar{t}}\), the blue shaded area to the \(68\%\) CL evaluated from the upper/lower bound of the integral of F(t). The orange shaded area corresponds to the \(95\%\) CL. The y-axis has arbitrary units

Fig. 2
figure 2

Combinations of transversities for the proton (upper) and the deuteron (lower) as extracted in Ref. [8]. The bottom plots show the statistical weight for each bin as evaluated through F(t)

How should we interpret the result shown in Fig. 1? It is unlikely that the bound is incorrect with a probability higher than \(t=10\%\) as most of the contribution to the integral of F(t) comes from the range \(0<t<0.1\). Still there exists a window through which the agreement of the implemented SB and the transversity combinations is poorer. Since the SB is introduced here as first principle, we translate this result into a relaxation of the expression of the bound through a Bayesian reweighting bin-by-bin in the objective function – i.e. here the \(\chi ^2\). The weight is obtained as follows,

$$\begin{aligned} w_j = \int _{1-p\left( A_{j}|\text{ data }\right) }^1\, F(t) \, dt\quad . \end{aligned}$$
(7)

In Fig. 2, the statistical 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 lowest and highest x-bins of the deuteron combination happen to statistically lie inside that window of poor agreement, as expected from previous analyses [7, 8].

In the bootstrap technique, the minimization of the objective function is performed \(N=200\) times. As explained in, e.g., Ref. [8], N replicas of the extracted combinations \(\left( x\, h_1^{p/D}\right) _j\), with \({j=1,...,N_d}\), of Eqs. (A.3A.4) are generated randomly within the data \(1\sigma \) error bars. As reweighting the overall chi-square function can be understood as a scaling of the error,

$$\begin{aligned} \sigma _j^2\rightarrow \sigma _j^2/w_j,\quad j=1,\ldots ,N_d\, \end{aligned}$$
(8)

the N replicas are generated within the corresponding extended Gaussian error bars. The objective function can be written, for each replica r, as

$$\begin{aligned}&\chi ^{2}_r\left( \{p^{I}\}\right) \nonumber \\&\quad =\sum _j \,w_j \frac{\left[ x_jh_{1\,\tiny {\text{ theo }}}^{p/D}\left( x_j ; \{p^{I}\}\right) - \left( xh_{1}^{p/D}(x)\right) _{j,\,r}\right] ^2}{\sigma _j^2}\;, \end{aligned}$$
(9)

where \(\{p^{I}\}\) is the set of free parameters to be determined for each replica r. In the next Section, the function \(h_{1\,\tiny {\text{ theo }}}^{p/D}\) will be extensively described.

The same exercise has been repeated with CT09MC1 [31] in both Eqs. (A.3A.4) and the expression of the SB. The lightest weights for the deuteron combination slightly change, repercussions will be discussed in Sect. 3. The overall conclusions of this section are not affected by the choice of LO unpolarized PDF.

2.2 Parameterization and constraints

Now that the main first principle based constraint has been included through a statistical reweighting – which 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. An unbiased parametrization is particularly welcome in kinematical ranges with little to no data or, a fortiori, kinematical ranges that exhibit apparent conflict with the positivity bounds. The latter are best accommodated through a flexible form that could realistically be contrasted with future data. We wish to impose the required support through the functional form. It is indeed possible for the up distribution. However, the behavior of \(x\,h_1^{d_v}\) is affected by bigger errors as it is dominated by the deuteron data, and a judicious choice of parametrization is not sufficient to ensure the expected behavior towards the end-point \(x\rightarrow 1\). A smooth fall-off in x is expected by the underlying QCD evolution occurring at higher x-values – see e.g. [32]. This observation directly relates to the power-law fall-off ruling the two other leading-twist PDFs. We tame undesired behaviors in the large-x region for the down distribution through a constrained fit using the method of Lagrange multipliers – see e.g. [33]. The minimization itself contains two steps. First, the objective function is minimized via a non-linear least-square, determining the set of best fit parameters, \(\{p^I\}\). The definition of the \( \chi ^2(\{p^I\})\), Eq. (9), requires a functional form for \(h_{1\,\tiny {\text{ theo }}}^{p/D}(x)\). The first considerations while determining the parameterization is to guarantee integrability and support at \(x=0\),

$$\begin{aligned} x\,h_1^{q_v}(x)= & {} x^{1.25}\times P_n^q\left( g(x)\right) , \end{aligned}$$
(10)

where the exponent has been chosen based on previous outputs and the expected small-x behavior [34] to be slightly modified by the polynomial in x. The latter, \(P_n^q\left( g(x)\right) \), of order n, can be as flexible as the data allow for. We choose to express \(P_n\) in terms of Bernstein polynomials as has been done in Ref. [35]. They are defined as

$$\begin{aligned} B_{k, n}(x)= & {} \left( \begin{array}{l}n\\ k \end{array}\right) x^k (1-x)^{n-k} , \end{aligned}$$
(11)

with \(\left( \begin{array}{l}n\\ k \end{array}\right) \) 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 with respect to 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. (10) 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 four different degrees for the polynomial \(P_n\) – thus using four different functional forms distinguishable 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. The functional form is now expressed as

$$\begin{aligned} x\,h_{1,i}^{q_v}\left( x ; p_{i,k}^q\right)= & {} x^{1.25}\,\sum _{k=\{\kappa _{q,i}\}} \,p_{i,k}^q\, B_{k, n_i}\left( g(x)\right) , \end{aligned}$$
(12)

with \(i=1,\dots ,4\), \(n=\{10, 20, 30, 40\}\) and where the \(p_{i,k}^q\) are the free parameters which number varies for each i and are listed in Table 1. The set of values \(\{\kappa _{q,i}\}\) has length \(n_{\kappa _{u,i}}=\{2,3,3,3\}\) and \(n_{\kappa _{d,i}}=\{4,3,4,4\}\) and their values have been chosen to flexibly cover the relevant region consistently with the data, as can be seen in Fig. 3. The four functional forms are shared among the N replicas, each \(x\,h_{1,i}^{q_v}\) will be evaluated N/4 times.

Table 1 Table with free parameters for each degree of the polynomials as described by Eq. (12)

In the following step, we will clarify our choice for a limited coverage of the functional form for values of \(x\gtrsim 0.6\) for the up and \(x\gtrsim 0.5\) for the down. As mentioned above, we need to ensure a smooth fall-off of the transversity in the limit \(x\rightarrow 1\) that, for the d contribution, 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_l(\{p'\})=0\) with \(l=1,\dots m\), the later are imposed through the Lagrangian

$$\begin{aligned} {\mathcal {L}}(\{p'\},\{\lambda \})= & {} \chi ^2(\{p'\})+\sum _{l}^m \lambda _l C_l(\{p'\}), \end{aligned}$$
(13)

to which a stationary point of \({\mathcal {L}}\) is found minimizing with respect to the parameters \(\{p'\}\) and the Lagrange parameters \(\{\lambda \}\).

Fig. 3
figure 3

Bernstein polynomials \(B_{k, n_i}\left( g(x)\right) \) used in the functional form for the valence up transversity (upper) and the valence down transversity (lower). The degree of the polynomials is, respectively, \(n=~\{10, 20, 30, 40\}\) in red with dot-dashed contours, purple/dashed, yellow/full and green/dotted. See text

We define a new objective function that depends on the new set of best fit parameters, \(\{p^{II}\}\), which consists in the set made of \(p^{q\,II}_{i,k}\), and replaces the set obtained through the first minimization, \(\{p^{I}\}\). We guide the large-x behavior of the down parameterization only, using the following \(N_c=4\) constraints

$$\begin{aligned} C_{i,l}^{d_v}\left( p^{d\,II}_{i,k}\right)= & {} x_l\,h_{1,i}^{d_v}\left( x_l ; p^{d\,II}_{i,k}\right) < \epsilon _l\nonumber \\&\qquad \text{ for } \quad l=1,\dots ,N_c/2, \nonumber \\ C_{i,l}^{d_v}\left( p^{d\,II}_{i,k}\right)= & {} x_l\,h_{1,i}^{d_v}\left( x_l ; p^{d\,II}_{i,k}\right) >- \epsilon _l \nonumber \\&\qquad \text{ for } \quad l=1,\dots ,N_c/2\;, \end{aligned}$$
(14)

with \(x_l=\{0.3, 0.55\}\) and \(\epsilon _l=\{0.2, 0.1\}\). In other words, we add 4 degrees of freedom to our problem. The values for \(\epsilon _l\) have been set considering the steepness of the functional forms and the trend of \(f_1(x,Q^2)\) and \(g_1(x,Q^2)\) through which is emulated the shift to small values of x induced by DGLAP.

In previous – unpolarized and longitudinally polarized – PDF determinations, the method of the Lagrange multipliers has been made popular for error estimation [36]. In the present approach, this method is used to impose limits on the fit parameters. This last step completes a methodology that focuses on the adaptability of the parametrization to constraints, in opposition to a parametrization that is constrained a priori.

Fig. 4
figure 4

Valence transversities: \(68\%\) envelope for the present fit for the four different degrees of the Bernstein polynomial of the functional form, respectively, \(n=10, 20, 30, 40\) in red with dot-dashed contours, purple/dashed, yellow/full and green/dotted, for the up (upper) and the down (lower)

Fig. 5
figure 5

Valence transversities for the up (upper) and the down (lower) at \(68\%\) in blue (with full blue line contour) and \(95\%\) in orange (with dashed contours) for the present fit. The gray bands represent the Soffer bound at NLO, here using JAM and MSTW08 at NLO

3 Results

3.1 Transversity PDF

Both steps of the minimization procedure are carried out in Python. The convergence for both the main \(\chi ^2\) minimization Eq. (9) 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 method of Lagrange multipliers. For the higher-order polynomials, less than 15 additional iterations are necessary while the two lower-order polynomials require about 30–40 extra iterations. Although the constraints are imposed for the down distribution only, their fulfillment 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 polynomials, \(n=10,20\) in red and purple, are confined in the mid- to large-x region. All 4 equally contribute at valence values of the Bjorken variable. The final envelope of \(xh_1^{q_v}(x,\langle 5 \text{ GeV }^2\rangle )\) 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 [30] and the unpolarized PDF of MSTW08 at NLOFootnote 2 [29]. Were there central values for the obtained transversities, they would be enclosed inside the Soffer bound. In that sense, our result is similar to the first Hessian approach of Ref. [7].

Fig. 6
figure 6

Combinations of transversities for the proton (upper) and the deuteron (lower) 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 [8]. The dark blue square come from HERMES data while the dark orange dots are extracted from COMPASS data

The small-x behavior predicted in Ref. [34] is fulfilled for the up transversity PDF up to, at least, \(x=0.1\) ; no clear conclusion can be drawn for the down above \(x\sim 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. (A.3A.4). The final \(\chi ^2/\)d.o.f. is evaluated against those extracted data. The number of degrees of freedom here is \((\sum _j\, w_j)\times N_d+N_c-(n_{\kappa _{u}}+n_{\kappa _{d}})\). The average value for the total \(\chi ^2/\)d.o.f. is 1.35. Each functional form contributes to average values of \(\chi _{n=10}^2/\)d.o.f.\(=1.32\), \(\chi _{n=20}^2/\)d.o.f.\(=1.44\), \(\chi _{n=30}^2/\)d.o.f.\(=1.37\), \(\chi _{n=40}^2/\)d.o.f.\(=1.26\). The four corresponding histograms are shown in Fig. 7.

The analysis has been carried out using the MSTW08 LO parametrization. We have repeated the fitting procedure using the CT09MC1 set for the unpolarized PDF [31], it is a LO set fitted on real data and NLO pseudo data with \(\alpha _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 error band in the data region – as our functional form was built on that purpose. It is achieved at the expense of a slightly wider band for the up distribution in that region. The error band increases for small- and large-x values. In this respect, our down distribution differs from that of Refs. [8, 9] and the less flexible versions of the parameterization of Ref. [7], all of them being dihadron-based extractions. The same comment is in order for comparisons with single-hadron semi-inclusive DIS [10, 11]. The combination \(h_1^{u_v}-h_1^{d_v}\) of the two available collinear extractions is compared in Fig. 8. The small-x behavior differs and the error band associated to the present analysis is wider in the valence region.

We also compare, in Fig. 8, our results to the lattice QCD evaluation of the isovector transversity PDF [38, 39]. There is an acceptable agreement in the region \(0.1<x<0.3\). The lower-x region of our result is more structured, as dictated by the data. Both lattice evaluations exceed the phenomenological determinations at large-x – region in which the parameterization of the fits are strongly bound by the positivity limits.

Fig. 7
figure 7

Histograms of chi-square per degree of freedom for each parametrization. The color code represent, as above, the four different degree of the Bernstein polynomial of the functional form, respectively, \(n=10, 20, 30, 40\) in red, purple, yellow and green

As explained in Appendix A, our analysis of the transversity PDF through Eqs. (A.3, A.4) has been carried out at an average value of \(Q^2=\langle 5\, \)GeV\(^2\rangle \). To that regard we have further performed our analysis considering QCD evolution only for the Fragmentation Functions. We find a 15% increase in the \(\chi ^2\) when fixing the scale of the unpolarized PDFs in the asymmetry to 5 GeV\(^2\). 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 [12] is systematically higher than the present results for the whole support in x, except for \(x\rightarrow 0\).

Fig. 8
figure 8

Upper plot: Comparison with the isovector combination of the transversity PDF of this work, in cyan meshed crosses, with the result at 4 GeV\(^2\) obtained through the global fit of Ref. [9], in meshed pink, both at 90 % CL. These collinear extractions are compared to lattice QCD evaluations from the ETM Collaboration [38], in gray, and the LP\(^3\) Collaboration [39], in dashed blue. Lattice results correspond to a scale of 4 GeV\(^2\). Lower plot: same as upper plot with logarithmic x-scale

3.2 Tensor charge

We next evaluate the first Mellin moment of the transversity PDF to get the tensor charge

$$\begin{aligned} {\delta {q}_{v}} \left( 5\, \text{ GeV }^2\right)= & {} \int _0^1 dx \, h_1^{q_v} \left( x; 5 \,\text{ GeV }^2\right) . \end{aligned}$$
(15)

The distributions correspond best to skewed Gaussian distributions and the obtained values are

$$\begin{aligned} {\delta u_{v}} \left( 5 \,\text{ GeV }^2\right)= & {} 0.28\begin{array}{c} +0.17 \\ -0.20 \end{array},\nonumber \\ {\delta d_{v}} \left( 5\, \text{ GeV }^2\right)= & {} -0.40\begin{array}{c} +0.41 \\ -0.31 \end{array}, \end{aligned}$$
(16)

to \(1\sigma \). The corresponding \(2\sigma \) uncertainties are, for the up, \((-0.42,0.31)\) and, for the down, \((-0.57,0.87)\). 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 [18, 19]. We show the corresponding stacked histogram in Fig. 9. The Gaussian distribution showed on the r.h.s. of Fig. 9 is given by

$$\begin{aligned} g_T \left( 5 \,\text{ GeV }^2\right)= & {} 0.57\pm 0.21, \end{aligned}$$
(17)

at \(1\sigma \) ; the \(2\sigma \) error is 0.42. The truncated tensor charge is found to be \(g_T^{\urcorner }( 5\, \)GeV\(^2)~=~0.48\begin{array}{c} +0.13 \\ -0.11 \end{array}\), corresponding to a skewed Gaussian distribution.

Fig. 9
figure 9

Stacked histogram of isovector tensor charge. The color code represent, as above, the four different degrees of the Bernstein polynomial of the functional form, respectively, \(n=10, 20, 30, 40\) in red, purple, yellow and green. Notice that the various contributions are represented on a single bar, i.e. stacked. On the lower panel, we can see that \(g_T\) is gaussianly distributed, as shown in red with the corresponding \(1\sigma \) error in gray

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 of the isovector tensor charge – the latest results are \(g_T=0.926(32)\) for the ETM Collaboration [23], \(g_T = 0.972(41)\) [24] and \(g_T = 0.965 (38) (+13 -41)\) [25]. It is at the same time wide enough to encompass all of the other phenomenological determinations. As for the separate up and down contributions, e.g., the ETM Collaboration reports \(\delta u=0.716(28)\) and \(\delta d=-0.210(11)\). The tensor charge for the up quark differs from our result as well as most phenomenologically determined tensor charges. The relaxation of the Soffer bound in the range of the data did not ease that discrepancy, as also discussed in Ref. [40].

In a region where scarce to no data are available, 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 [41], used in the first global fit of the transversity [9], 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 extremely 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. [19].

4 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 its 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 [6]. It requires the knowledge on the Dihadron Fragmentation Functions that have been studied in Ref. [8, 42]. Since its first point-by-point extraction from dihadron semi-inclusive DIS [43], enormous efforts towards the extraction of the transversity PDF, accessible in processes involving fragmentation to a pion-pair [9], have been made.

In this paper we have presented a constrained fit of the valence transversity PDFs from dihadron semi-inclusive DIS data. The adopted methodology is characterized by the flexibility to accommodate constraints from the fitting procedure. It consists in three steps. First we have examined the positivity constraints on the transversity distribution: a priori from the expression of the bound combined with the data is introduced as a theoretical uncertainty on the PDF that is then implemented through a reweighting of the data. Second we have chosen flexible functional forms – free from ad hoc expressions of the Soffer bound – to enhance the information obtained from the data and first principles. It results in an improved treatment of the PDF uncertainties. Finally, we have guided the valence down transversity PDF to follow a fall-off behavior at \(x\rightarrow 1\) through the method of Lagrange multipliers, forcing the set of parameters to observe the theoretical constraints.

The obtained \(68\%\) and \(95\%\) CL envelopes for the valence transversity distribution functions fulfill the expected small-x behavior, the up distribution fully realizes 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 error bands outside the data range with respect to the error band 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. The resulting isovector tensor charge is in agreement with the determinations of lattice QCD. We believe that this exercise 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.

At SIDIS scale the relative effect of DGLAP is small compared to the accurracy of the obtained PDFs. The effect of QCD evolution will become important when including data from proton-proton collision from RHIC for which the role of NLO corrections might become important, as discussed in Ref. [9]. Our methodology could be applied, as a natural extension of this work, to the extended set of data including the aforementioned proton-proton data and the inclusion of the corresponding DGLAP routine. 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.