††thanks<EMAIL_ADDRESS>- ORCID: 0000-0002-0762-0275

# Simultaneous Monte Carlo analysis of parton densities
and fragmentation functions

E. Moffat<EMAIL_ADDRESS>Department of Physics, Old Dominion University,
Norfolk, Virginia 23529, USA  W. Melnitchouk<EMAIL_ADDRESS>Jefferson Lab,
Newport News, Virginia 23606, USA
Jefferson Lab Angular Momentum (JAM) Collaboration  T. C. Rogers  Department
of Physics, Old Dominion University, Norfolk, Virginia 23529, USA  Jefferson
Lab, Newport News, Virginia 23606, USA
Jefferson Lab Angular Momentum (JAM) Collaboration  N. Sato<EMAIL_ADDRESS>Jefferson Lab, Newport News, Virginia 23606, USA
Jefferson Lab Angular Momentum (JAM) Collaboration

###### Abstract

We perform a comprehensive new Monte Carlo analysis of high-energy lepton-
lepton, lepton-hadron and hadron-hadron scattering data to simultaneously
determine parton distribution functions (PDFs) in the proton and parton to
hadron fragmentation functions (FFs). The analysis includes all available
semi-inclusive deep-inelastic scattering and single-inclusive $e^{+}e^{-}$
annihilation data for pions, kaons and unidentified charged hadrons, which
allows the flavor dependence of the fragmentation functions to be constrained.
Employing a new multi-step fitting strategy and more flexible parametrizations
for both PDFs and FFs, we assess the impact of different data sets on sea
quark densities, and confirm the previously observed suppression of the
strange quark distribution. The new fit, which we refer to as “JAM20-SIDIS”,
will allow for improved studies of universality of parton correlation
functions, including transverse momentum dependent (TMD) distributions, across
a wide variety of process, and the matching of collinear to TMD factorization
descriptions.

††preprint: JLAB-THY-21-3304

##  I Introduction

The standard parton correlation functions of QCD, such as collinear parton
distribution functions (PDFs) and fragmentation functions (FFs), are being
utilized in an increasingly diverse range of phenomenological applications.
Beyond their traditional role in predicting new high energy phenomena, they
also enter frequently into the study of more complex and extended objects like
transverse momentum dependent (TMD) PDFs and FFs and generalized parton
distributions (GPDs), where they are needed to understand the transition
between different factorization regions. Both TMDs and GPDs are central to the
study of the nonperturbative parton structure of hadrons, and understanding
how they encapsulate their longitudinal and transverse features will be
critical to current experimental programs at Jefferson Lab and elsewhere, as
well as to the future Electron-Ion Collider. These considerations provide one
of the main motivations for the study of collinear PDFs and FFs in this paper.

The great value of PDFs and FFs extracted from global QCD data analysis lies
with their predictive power, or “universality”. However, the translation from
experimental data to quark and gluon operator structures is a challenging
inverse problem. It is not possible to exactly constrain parton correlations
from data alone since this connection involves nontrivial convolution
integrals in a factorization formalism (whose accuracy itself is difficult to
quantify in any given instance), and because of the limited quantity of
available data. The complexity of the inverse problem is also magnified by the
number of flavor degrees of freedom involved.

Nevertheless, assessing and maximizing the universality of collinear PDFs and
FFs is crucial given the increasingly broad scenarios where they are used. A
major focus in the current effort by the Jefferson Lab Angular Momentum (JAM)
Collaboration is therefore to both test and broaden the predictive power of
parton correlation functions. This is achieved through a Bayesian inference
procedure in which PDFs and FFs are extracted simultaneously, and the
uncertainty quantification associated with particular parametrizations of
parton correlation functions is given in terms of a Bayesian posterior
distribution. To test universality, the system of equations relating
observables to parton correlation functions must of course exceed the total
number of correlation functions involved — a minimum requirement is that the
parton correlation functions be overconstrained by the data in the fit. Of
course, realizing this in practical analyses requires that all parton
correlation functions be truly fitted simultaneously. This is a major
numerical and technological challenge, and traditionally PDFs and FFs have
thus been extracted in separate procedures. However, simultaneous fits can be
achieved with the Bayesian Monte Carlo approach, and have been implemented
recently in the JAM17 Ethier  _et al._ (2017) analysis of helicity PDFs, and
in the JAM19 Sato  _et al._ (2020) analysis of unpolarized PDFs and FFs. The
same basic methodology was also applied in the three-dimensional JAM3D20
Cammarota  _et al._ (2020) study, in the first combined analysis of TMD
observables that satisfies the overconstraining criterion.

In this paper, we extend the previous work by performing the first
simultaneous and overconstrained fit of unpolarized PDFs and FFs that utilizes
both charged hadron production in semi-inclusive deep-inelastic scattering
(SIDIS) and single-inclusive $e^{+}e^{-}$ annihilation (SIA). This is partly
motivated by a number of recent observations associated with the study of TMD
PDFs. For example, significant tension has recently been found between fits
performed with standard sets of PDFs and FFs and fixed order perturbative QCD
calculations in processes including SIDIS Gonzalez-Hernandez  _et al._ (2018);
Wang  _et al._ (2019), Drell-Yan (DY) Bacchetta  _et al._ (2019), and SIA into
wide-angle hadron pairs Moffat  _et al._ (2019). A number of suggested
solutions and explanations have been proposed to account for this, including a
possible need for power suppressed corrections Liu and Qiu  at the moderate
scales of most SIDIS experiments. However, more tests of the limits of
applicability of standard collinear factorization are needed before it is
possible to draw firm conclusions. Given that the majority of data used to
constrain collinear correlation functions (both PDFs and FFs) are either
highly inclusive or exist are at very high scales, or both, it is perhaps not
surprising that tension arises when these are evolved downward and used to
make predictions at lower scales and for highly differential observables.
Indeed, there have been few tests that $Q^{2}$-scaling, a hallmark of the
collinear perturbative regime, actually holds to a reasonable approximation in
SIDIS measurements at moderate $Q^{2}$. Our hope is that the new combined fit,
which we refer to as “JAM20-SIDIS”, will help to shed light on this and
similar issues in the future.

In Sec. II we begin the discussion by summarizing the methodology used in our
simultaneous Monte Carlo analysis, including the parametrizations used for the
distributions and the multi-step Bayesian inference algorithm. Details of the
data sets included in the fit are summarized in Sec. III, while in Sec. IV we
discuss the criteria for universality and how these are met in this analysis.
A detailed discussion of the numerical results is given in Sec. V, where we
present the fitted PDFs and FFs, as well as detailed comparisons of data to
theory. Finally, in Sec. VI we summarize our conclusions and discuss the
implications of our analysis. Some formulas for SIDIS cross sections and
structure functions are collected in Appendix A.

##  II Theoretical framework

In this section we give an overview of the theoretical framework on which our
analysis is based, including the observables to be fitted, the
parametrizations used for the PDFs and FFs, details of the perturbative QCD
setup, and Bayesian inference strategy employed.

###  II.1 Observables and factorization

In this analysis we work in standard collinear factorization Collins  _et al._
(1989); Ellis  _et al._ (1996); Collins (2011), in which QCD cross sections
are separated into perturbatively calculable partonic hard factors convoluted
with nonperturbative PDFs and/or FFs. We perform calculations of all
observables consistently to order $\alpha_{s}$ in the QCD coupling. Details of
the basic theoretical setups for the inclusive DIS, inclusive Drell-Yan
lepton-pair production and SIA reactions are provided in the literature Ellis
_et al._ (1996); Devenish and Cooper-Sarkar (2004), and will not be repeated
here. However, since SIDIS is a comparatively novel addition to global QCD
analyses, we review it in more detail in Appendix A.

The processes considered in the present analysis can be summarized as follows:

$\displaystyle\ell+N$ $\displaystyle{}\to\ell+X,\qquad$ inclusive DIS
$\displaystyle\ell+N$ $\displaystyle{}\to\ell+h^{\pm}+X,\qquad$ semi-inclusive
DIS $\displaystyle N_{1}+N_{2}$ $\displaystyle{}\to\ell^{+}+\ell^{-}+X,\qquad$
Drell-Yan lepton-pair production $\displaystyle\ell^{+}+\ell^{-}$
$\displaystyle{}\to h^{\pm}+X,\qquad$ single-inclusive annihilation

where $h^{\pm}$ represent charged pions, kaons, or unidentified hadrons, and
the nucleon $N$ (or $N_{1,2}$) in the initial state can be either a proton or
a neutron (in practice, deuteron). Within the framework of collinear
factorization, the cross sections for each of these processes can be written
schematically as convolutions of hard functions and the nonperturbative parton
distribution and fragmentation functions,

$\displaystyle\frac{\mathrm{d}\sigma_{\text{\tiny
DIS}}}{\mathrm{d}Q^{2}{}\mathrm{d}x_{\mbox{\tiny\rm Bj}}{}}$
$\displaystyle{}=\sum_{i}\mathcal{H}_{i}^{\text{\tiny DIS}}\otimes
f_{i}\,,\qquad$ (1) $\displaystyle\frac{\mathrm{d}\sigma_{\text{\tiny
SIDIS}}{}}{\mathrm{d}Q^{2}{}\mathrm{d}x_{\mbox{\tiny\rm
Bj}}{}\mathrm{d}z_{h}{}}$
$\displaystyle{}=\sum_{ij}\mathcal{H}_{ij}^{\text{\tiny SIDIS}}\otimes
f_{i}\otimes D_{j}^{h}\,,\qquad$ (2)
$\displaystyle\frac{\mathrm{d}\sigma_{\text{\tiny
DY}}{}}{\mathrm{d}Q^{2}{}\mathrm{d}x_{\mbox{\tiny\rm F}}{}}$
$\displaystyle{}=\sum_{ij}\mathcal{H}_{ij}^{\text{\tiny DY}}\otimes
f_{i}\otimes f_{j}\,,\qquad$ (3)
$\displaystyle\frac{\mathrm{d}\sigma_{\text{\tiny
SIA}}{}}{\mathrm{d}Q^{2}{}\mathrm{d}z_{h}{}}$
$\displaystyle{}=\sum_{j}\mathcal{H}_{j}^{\text{\tiny SIA}}\otimes
D_{j}^{h}\,,\qquad$ (4)

where the symbols $\otimes$ represent the convolution integrals in
longitudinal momentum fractions of the hard scattering functions
$\mathcal{H}_{ij}$ and the PDFs $f_{i}$ and FFs $D_{j}^{h}$ for parton flavors
$i,j$ (see the Appendix). In each process, $Q$ represents the hard scale given
by the photon virtuality, $Q\gg$ hadron masses, which allows the observables
to be factorized into the short-distance perturbative and long-distance
nonperturbative parts.

For the inclusive DIS and SIDIS processes,

$x_{\mbox{\tiny\rm Bj}}=\frac{Q^{2}}{2p\cdot q}$ (5)

is the usual Bjorken scaling variable, while for the DY process the analogous
scaling variables are defined as

$x_{1}=\frac{Q^{2}}{2p_{1}\cdot q}\,,\qquad x_{2}=\frac{Q^{2}}{2p_{2}\cdot
q}\,,$ (6)

where $p_{1}$ and $p_{2}$ denote the incoming hadron momenta, with the Feynman
scaling variable given by

$x_{\mbox{\tiny\rm F}}=x_{1}-x_{2}.$ (7)

In the DY center of mass frame, and in the limit of negligible hadron masses
($\ll Q$), the virtual photon rapidity can be written in terms of $x_{1}$ and
$x_{2}$ as

$y=\frac{1}{2}\ln\frac{x_{1}}{x_{2}}.$ (8)

For the processes involving fragmentation to a hadron $h$ in the final state,
we have

$\hskip 113.81102ptz_{h}=\frac{p_{h}\cdot p}{q\cdot p}\hskip
113.81102pt\textrm{[SIDIS]}$ (9)

for SIDIS in Eq. (2), while

$\hskip 102.43008ptz_{h}=\frac{2p_{h}\cdot q}{Q^{2}}\hskip
108.12054pt\textrm{[SIA]}$ (10)

for SIA in Eq. (4).

###  II.2 Perturbative QCD and numerical setups

For our numerical analysis we make use of Mellin space techniques to enable
fast evaluations of observables needed for the Bayesian analysis. In
particular, we solve the DGLAP evolution equations analytically in Mellin
space Vogt (2005), which allows one to effectively render high-dimensional
momentum space convolutions from process-specific factorization theorems,
along with the integrals in the DGLAP equations, in the form of lower-
dimensional inverse Mellin transforms. For example, for the inclusive DIS
observables one can write schematically

$\displaystyle\frac{\mathrm{d}\sigma_{\text{\tiny
DIS}}}{\mathrm{d}Q^{2}\mathrm{d}x_{\mbox{\tiny\rm Bj}}}=\sum_{i}\frac{1}{2\pi
i}\int dN\,x_{\mbox{\tiny\rm Bj}}^{-N}\,\widetilde{\cal H}^{\text{\tiny
DIS}}_{i}(N,\mu)\,U^{\rm
S}_{ij}(N,\mu,\mu_{0})\,\widetilde{f}_{j}(N,\mu_{0}),$ (11)

where $N$ here is the conjugate variable to $x_{\mbox{\tiny\rm Bj}}$,
$\widetilde{f}_{j}(N,\mu_{0})$ is the Mellin moment of the PDF
$f_{j}(x,\mu_{0})$, defined by

$\widetilde{f}_{j}(N,\mu_{0})=\int_{0}^{1}dx\,x^{N-1}\,f_{j}(x,\mu_{0}),$ (12)

and $\widetilde{\cal H}^{\text{\tiny DIS}}_{i}(N,\mu)$ is the corresponding
moment of the partonic DIS cross section. The analytic solution for the DGLAP
evolution is entirely encoded in the evolution matrix $U^{\rm S}_{i,j}$ that
evolves the moments $\widetilde{f}_{j}(N,\mu_{0})$ of the PDFs from a given
input scale $\mu_{0}$ to the relevant DIS hard scale $\mu=Q$. A similar
expression can be written for the SIA cross section,

$\displaystyle\frac{\mathrm{d}\sigma_{\text{\tiny
SIA}}}{\mathrm{d}Q^{2}\mathrm{d}z_{h}}=\sum_{ij}\frac{1}{2\pi i}\int
dM\,z_{h}^{-M}\,\widetilde{\cal H}^{\text{\tiny SIA}}_{i}(M,\mu)\,U^{\rm
T}_{ij}(M,\mu,\mu_{0})\,\widetilde{D}_{j}^{h}(M,\mu_{0})$ (13)

where $M$ is the Mellin conjugate variable for $z_{h}$,
$\widetilde{D}_{j}^{h}(M,\mu_{0})$ is the moment of the FF, and
$\widetilde{\cal H}^{\text{\tiny SIA}}_{i}$ is the moment of the partonic SIA
cross section. The superscripts $\rm S$ and $\rm T$ in the evolution matrix
distinguish between the spacelike and timelike evolution for the PDFs and FFs,
respectively, which are encoded in the corresponding DGLAP splitting kernels.

The same procedure can be extended for the case of SIDIS, which gives

$\displaystyle\frac{\mathrm{d}\sigma_{\text{\tiny
SIDIS}}}{\mathrm{d}Q^{2}\mathrm{d}x_{\mbox{\tiny\rm Bj}}\mathrm{d}z_{h}}=$
$\displaystyle\sum_{ijkl}\frac{1}{(2\pi i)^{2}}\int dN\,x_{\mbox{\tiny\rm
Bj}}^{-N}\int dM\,z_{h}^{-M}\,\widetilde{\cal H}^{\text{\tiny
SIDIS}}_{ik}(N,M,\mu)$ $\displaystyle\times$ $\displaystyle\ U^{\rm
S}_{ij}(N,\mu,\mu_{0})\,\widetilde{f}_{j}(\mu_{0})\,U^{\rm
T}_{kl}(M,\mu,\mu_{0})\,\widetilde{D}_{j}^{h}(M,\mu_{0}).$ (14)

For the case of the Drell-Yan process, a special treatment is required since
the Mellin moments for the partonic cross sections are not known. For this we
employ the strategy developed by Stratmann and Vogelsang Stratmann and
Vogelsang (2001), where by the Mellin moments are numerically pre-calculated
and used as lookup tables during the analysis. The resulting expression can be
written schematically as

$\displaystyle\frac{\mathrm{d}\sigma_{\text{\tiny
DY}}}{\mathrm{d}Q^{2}\mathrm{d}x_{\mbox{\tiny\rm F}}}=$
$\displaystyle\sum_{ijkl}\frac{1}{(2\pi i)^{2}}\int dN\,\int
dM\,\widetilde{\cal H}^{\text{\tiny DY}}_{ik}(N,M,\mu)$ $\displaystyle\times$
$\displaystyle\ U^{\rm
S}_{ij}(N,\mu,\mu_{0})\,\widetilde{f}_{j}(\mu_{0})\,U^{\rm
S}_{kl}(M,\mu,\mu_{0})\,\widetilde{f}_{l}(\mu_{0}),$ (15)

where the relevant inverse Mellin factors $x_{1}^{-N}$ and $x_{2}^{-M}$
arising from the scaling variables $x_{1}$ and $x_{2}$ for the incident
nucleons $N_{1}$ and $N_{2}$, respectively, in Eq. (3) are integrated
numerically with the hard scattering cross section and contained inside
$\widetilde{\cal H}^{\text{\tiny DY}}_{ik}(N,M,\mu)$.

The analytic solutions for the evolution matrices are computed at next-to-
leading logarithmic accuracy using splitting kernels up to ${\cal
O}(\alpha_{s}^{2})$ and the truncated solution for the single evolution
operators (see Ref. Vogt (2005) for details). We employ the zero-mass variable
flavor scheme for solving the DGLAP evolution equations, setting the input
scale for the PDFs and FFs at $\mu_{0}=m_{c}$. The numerical values for the
mass thresholds are taken from the PDG values in the $\overline{\rm MS}$
scheme Zyla  _et al._ (2020): $m_{c}=1.28$ GeV and $m_{b}=4.18$ GeV. The
strong coupling is evolved numerically using the QCD beta function up to
${\cal O}(\alpha_{s}^{2})$, using the boundary condition
$\alpha_{s}(M_{Z})=0.118$ at the $Z$ boson mass, $M_{Z}=91.18$ GeV. Finally,
all the process specific hard coefficients are computed at fixed next-to-
leading order in pQCD, which are available in the literature Floratos  _et
al._ (1981); Kretzer (2000); Stratmann and Vogelsang (2001); Anastasiou  _et
al._ (2003).

###  II.3 Parametrization of nonperturbative functions

For the nonperturbative parton distribution and fragmentation functions we use
standard parametrizations that have been utilized in the literature. Namely,
for the dependence on the parton momentum fraction $x$ of the PDF $f(x)$ we
use the template function

$f(x,\mu_{0})\ \to\ T\left(x;\bm{a}\right)={\cal
M}\frac{x^{\alpha}\left(1-x\right)^{\beta}\left(1+\gamma\sqrt{x}+\delta{x}\right)}{\int_{0}^{1}\mathrm{d}x\,x^{\alpha+1}\left(1-x\right)^{\beta}\left(1+\gamma\sqrt{x}+\delta{x}\right)},$
(16)

where $\bm{a}=\\{{\cal M},\alpha,\beta,\gamma,\delta\\}$ is a vector
containing the shape parameters ($\alpha$, $\beta$, $\gamma$, and $\delta$)
and a normalization coefficient (${\cal M}$) to be fitted. The integral in the
denominator ensures that the value of the normalization coefficient ${\cal M}$
is equal to the second moment ($x$-weighted integral) of the function
$T(x;\bm{a})$. For fitting the PDFs, we assume isospin symmetry to relate the
PDFs in the neutron, $f_{i/n}(x)$, to those in the proton, $f_{i/p}(x)\equiv
f_{i}(x)$, switching the $u\leftrightarrow d$ and $p\leftrightarrow n$ labels
for the light quark flavors, and taking the PDFs for other flavors equal for
the proton and neutron.

In practice, we parametrize the valence $u$ and $d$ quark distributions,
$u_{v}\equiv f_{u}-f_{\bar{u}}$ and $d_{v}\equiv f_{d}-f_{\bar{d}}$, directly
using the template function (Eq. (16)). The gluon distribution, $g\equiv
f_{g}$, is also directly parametrized per Eq. (16). For the sea quark and
antiquark distributions, we use five functions parametrized as in Eq. (16).
These are a flavor symmetric sea function ($S$) that dominates at very low $x$
and flavor specific functions ($q_{0}(\bar{q}_{0})$) for the $s$, $\bar{u}$,
$\bar{d}$, and $\bar{s}$ that take into account the possible nonperturbative
origin of the sea. The distributions for $s$, $\bar{u}$, $\bar{d}$, and
$\bar{s}$ are constructed from these according to: $q(\bar{q})\equiv
f_{q(\bar{q})}=S+q_{0}(\bar{q}_{0})$. Note that $s$ and $\bar{s}$ are
parametrized separately because their contributions to the $K^{+}$ and $K^{-}$
SIDIS cross sections differ. We do not fit the charm and bottom PDFs, and
their contributions are generated purely from the DGLAP evolution. In total
there are 8 parametrized PDF functions being fitted. For the valence quark
PDFs $u_{v}$ and $d_{v}$ and the nonperturbative sea components $\bar{u}_{0}$
and $\bar{d}_{0}$, we use the four shape parameters as in Eq. (16); for all
other distributions we set the $\gamma$ and $\delta$ parameters to zero. This
gives 24 free shape parameters and 8 free normalization parameters. The number
of free parameters is further reduced by valence number sum rules, which
constrain the normalization parameters $\cal M$ for the $u_{v}$, $d_{v}$, and
$s-\bar{s}$ distributions, whose lowest moments are required to be 2, 1, and
zero, respectively. The normalization for the gluon PDF is determined using
the momentum sum rule. With these constraints, there is a total of 28 free
parameters for the PDFs.

For the $z$ dependence of FFs, the functional form follows a similar template,

$D(z,\mu_{0})\ \to\ T\left(z;\bm{a}\right)={\cal
M}\frac{z^{\alpha}\left(1-z\right)^{\beta}\left(1+\gamma\sqrt{z}+\delta{z}\right)}{\int_{0}^{1}\mathrm{d}z\,z^{\alpha+1}\left(1-z\right)^{\beta}\left(1+\gamma\sqrt{z}+\delta{z}\right)},$
(17)

where again the integral in the denominator ensures that the coefficient
${\cal M}$ corresponds to the second moment ($z$-weighted integral) of the
function. In addition to the fragmentation to pions and kaons studied in
earlier JAM analyses of SIA and SIDIS data Sato  _et al._ (2016, 2020), here
we consider also the inclusive production of unidentified charged hadrons,
$h^{\pm}$. Accounting for unidentified hadrons can be implemented in two ways.
First, the hadron FFs can be fit independently from those for pions and kaons,
as preferred by the NNPDF Collaboration Bertone  _et al._ (2018).
Alternatively, one can take advantage of existing knowledge of specified
hadron FFs and add a fitted residual correction to their sum. Such an approach
was adopted by de Florian, Sassot, and Stratmann (DSS) de Florian  _et al._
(2007), for example, in which a residual correction was fitted to the sum of
previously obtained pion, kaon, and proton fragmentation functions.

In our analysis we follow the latter approach, but include only the pion and
kaon FFs, so that the residual term $D_{i}^{\rm res^{+}}$ parametrizes the
difference between the total hadron FF $D_{i}^{h^{+}}$ and the
$D_{i}^{\pi^{+}}$ and $D_{i}^{K^{+}}$ functions,

$D_{i}^{h^{+}}=D_{i}^{\pi^{+}}+D_{i}^{K^{+}}+D_{i}^{\rm res^{+}}.$ (18)

To reduce the total number of residual FFs being fit, we assume SU(3) flavor
symmetry for light quarks and antiquarks,

$\displaystyle D_{q}^{\rm res^{+}}=D_{u}^{\rm res^{+}}$ $\displaystyle=$
$\displaystyle D_{d}^{\rm res^{+}}=D_{s}^{\rm res^{+}},\ \ \ \ \ \ \ $ (19a)
$\displaystyle D_{\bar{q}}^{\rm res^{+}}=D_{\bar{u}}^{\rm res^{+}}$
$\displaystyle=$ $\displaystyle D_{\bar{d}}^{\rm res^{+}}=D_{\bar{s}}^{\rm
res^{+}},$ (19b)

where $D_{q}^{\rm res^{+}}$ and $D_{\bar{q}}^{\rm res^{+}}$ are parametrized
per the template (Eq. (17)). To allow for differentiation between the residual
FFs for light quarks and antiquarks, we leave ${\cal M}$ and $\beta$ for
$D_{\bar{q}}^{\rm res^{+}}$ as free parameters, but fix $\alpha$, $\gamma$,
and $\delta$ to be the same as for $D_{q}^{\rm res^{+}}$. This achieves a
similar constraint on the parameters as the condition used by DSS de Florian
_et al._ (2007), $2D_{\bar{q}}^{\rm
res^{+}}=\left(1-z\right)^{\beta^{\prime}}D_{q+\bar{q}}^{\rm res^{+}}$.

For the pion FFs, $D_{i}^{\pi^{+}}$, we reduce the number of fitted functions
by grouping the light quarks into “favored” (valence) and “unfavored” (non-
valence) flavors,

$\displaystyle D_{\rm fav}^{\pi^{+}}$
$\displaystyle=D_{u}^{\pi^{+}}=D_{\bar{d}}^{\pi^{+}},$ (20a) $\displaystyle
D_{\rm unf}^{\pi^{+}}$
$\displaystyle=D_{d}^{\pi^{+}}=D_{s}^{\pi^{+}}=D_{\bar{u}}^{\pi^{+}}=D_{\bar{s}}^{\pi^{+}},$
(20b)

where $D_{\rm fav}^{\pi^{+}}$ and $D_{\rm unf}^{\pi^{+}}$ are parametrized as
in Eq. (17). For the parameters of the kaon FFs, $D_{i}^{K^{+}}$, we equate
the “unfavored” flavors,

$D_{\rm
unf}^{K^{+}}=D_{d}^{\pi^{+}}=D_{s}^{\pi^{+}}=D_{\bar{u}}^{\pi^{+}}=D_{\bar{d}}^{\pi^{+}},$
(21)

but leave the favored FFs $D_{u}^{K^{+}}$ and $D_{\bar{s}}^{K^{+}}$
independent. Here $D_{\rm unf}^{K^{+}}$, $D_{u}^{K^{+}}$, and
$D_{\bar{s}}^{K^{+}}$ are parametrized per Eq. (17). For the heavier flavors,
we assume the charm and bottom quark and antiquark FFs to be equivalent,
$D_{c}^{h^{+}}=D_{\bar{c}}^{h^{+}}$ and $D_{b}^{h^{+}}=D_{\bar{b}}^{h^{+}}$
for $h=\pi,K,{\rm res}$, with $D_{c}^{h^{+}}$ and $D_{b}^{h^{+}}$ parametrized
per Eq. (17). Finally, the gluon FFs $D_{g}^{h^{+}}$ for $h=\pi,K,{\rm res}$
are also parametrized according to Eq. (17). We use charge conjugation
symmetry to relate FFs for opposite charges by

$D_{q}^{h^{+}}=D_{\bar{q}}^{h^{-}},$ (22)

where $h=\pi,K,{\rm res}$. This results in 5 fitted functions for pions and
residual hadrons, and 6 for kaons.

At this point, there are 17 shape parameters and 5 normalization parameters
for residual hadrons, 20 shape parameters, and 5 normalization parameters for
pions, and 24 shape parameters and 6 normalization parameters for kaons. The
number of shape parameters is reduced further because throughout the fitting
procedure, the parameters $\gamma$ and $\delta$ for the gluon, charm, and
bottom FFs are fixed at zero. In the end there are 16 free parameters to be
fitted for residual charged hadron FFs, 19 free FF pion parameters, and 24
free parameters for the kaon FFs. Together with the 28 PDF parameters, we have
a total of 87 free parameters for the fitted functions. In addition, there are
also 42 free parameters associated with normalization of various data sets,
making for a total of 129 free parameters to be fitted in the analysis.

###  II.4 Bayesian inference

Our methodology for extracting nonperturbative PDFs and FFs is based on the
general premise of Bayesian inference. Namely, we use Bayes’ theorem to define
a multivariate probability distribution ${\cal P}$ for the shape parameters
characterizing the PDFs and FFs (the posterior) at a given input scale
$\mu_{0}$,

$\displaystyle{\cal P}(\bm{a}|{\rm data})\sim{\cal L}(\bm{a},{\rm
data})\,\pi(\bm{a}),$ (23)

where $\cal L$ is a standard Gaussian likelihood function,

$\displaystyle{\cal L}\left(\bm{a},{\rm
data}\right)=\exp\left(-\frac{1}{2}\chi^{2}\left(\bm{a},{\rm
data}\right)\right),$ (24)

with the $\chi^{2}$ function defined by

$\chi^{2}(\bm{a})=\sum_{i,e}\bigg{(}\frac{d_{i,e}-\sum_{k}r_{e}^{k}\beta^{k}_{i,e}-T_{i,e}(\bm{a})/N_{e}}{\alpha_{i,e}}\bigg{)}^{2}+\sum_{k}\big{(}r^{k}_{e}\big{)}^{2}+\left(\frac{1-N_{e}}{\delta
N_{e}}\right)^{2}.$ (25)

Here, $d_{i,e}$ is the value of the $i$-th data point for the experimental
dataset $e$, with $T_{i,e}$ the theoretical prediction for the data point;
$\alpha_{i,e}$ is the uncorrelated systematic and statistical uncertainty for
each data point added in quadrature; $\beta^{k}_{i,e}$ is the $k$-th source of
point-to-point correlated systematic uncertainties for the $i$-th bin of
dataset $e$, with $r^{k}_{e}$ the related weight; and $N_{e}$ and
$\delta{N}_{e}$ are the normalization and normalization uncertainty for each
data set, respectively. In Eq. (23), $\pi(\bm{a})$ is the prior distribution
for the set of parameters $\bm{a}$, which is used as input for a given fit to
the data.

In principle, given the Bayesian posterior distribution, one can estimate
confidence regions for a generic observable ${\cal O}$ (such as a PDF or a
function of PDFs or FFs) by integrating over an $d$-dimensional parameter
space,

$\displaystyle{\rm E}[\mathcal{O}]$
$\displaystyle=\int\mathrm{d}^{d}\bm{a}\,\mathcal{P}\left(\bm{a}|{\rm
data}\right)\mathcal{O}\left(\bm{a}\right),$ (26a) $\displaystyle{\rm
V}[\mathcal{O}]$
$\displaystyle=\int\mathrm{d}^{d}\bm{a}\,\mathcal{P}\left(\bm{a}|{\rm
data}\right)\left(\mathcal{O}\left(\bm{a}\right)-{\rm
E}[\mathcal{O}]\right)^{2},$ (26b)

where E and V are the expectation value and variance of the observable ${\cal
O}$, respectively. Due to the significant numerical expense of evaluating the
likelihood function, the explicit usage of Eqs. (26) is often not practical.
Instead, a more efficient option is to build Monte Carlo parameter samples
$\\{\bm{a}_{k};k=1,\ldots,n\\}$, which contain all parameters, including the
$N_{e}$ from Eq. (25), that are faithfully distributed according to the
posterior distribution. These can in turn be used to evaluate the integrals in
Eqs. (26) as Monte Carlo sums,

$\displaystyle{\rm E}[\mathcal{O}]$
$\displaystyle=\frac{1}{n}\sum_{k=1}^{n}\mathcal{O}\left(\bm{a}_{k}\right),$
(27a) $\displaystyle{\rm V}[\mathcal{O}]$
$\displaystyle=\frac{1}{n}\sum_{k=1}^{n}\left(\mathcal{O}\left(\bm{a}_{k}\right)-{\rm
E}[\mathcal{O}]\right)^{2}.$ (27b)

Figure 1: Schematic illustration of the multi-step workflow employed in our
simultaneous Monte Carlo analysis. Each box represents a collection of Monte
Carlo samples associated with a specific nonperturbative hadronic structure
(PDFs, FFs). The vertical arrows indicate the inclusion of additional datasets
from which new optimized Monte Carlo samples (posteriors) are generated as
input (priors) for the next step.

Our Monte Carlo sampling strategy is based on data resampling methodology,
whereby multiple maximum likelihood optimizations are carried out. Each
optimization consists of taking a random point in parameter space and fitting
the parameters to data that have been distorted away from the central values
by Gaussian shifts within the quoted uncertainties. To build the Monte Carlo
samples, we use the multi-step strategy developed in Ref. Sato  _et al._
(2020), where the PDF and FF parameters are pre-optimized to minimize
evaluating the likelihood in parameter regions that are strongly disfavored.
To that end we start by first considering PDF and FF parameters separately
using flat priors, with the resulting samples from each type of hadron
structure combined at a later stage to build new prior samples for the final
runs. The workflow is illustrated in Fig. 1, where each step is represented as
vertical arrows that accumulate additional experimental data from the previous
step, with the posterior samples at each step becoming the priors for the
subsequent step. This strategy allows the samples to become more optimized and
avoids unnecessary likelihood evaluations in regions of parameters space by
disfavoring those regions in earlier stages of the multi-step chain.

##  III Data Sets

The data sets used in the present analysis include the primary electromagnetic
processes that traditionally have been used in global QCD analyses, namely,
inclusive DIS, Drell-Yan lepton-pair production (which constrain PDFs), SIA
(which constrains FFs), and SIDIS (which constrains both PDFs and FFs). The
inclusive DIS data are measurements of the $F_{2}(x_{\mbox{\tiny\rm
Bj}},Q^{2})$ structure function performed by the BCDMS Benvenuti  _et al._
(1989, 1990) and New Muon Collaborations Arneodo  _et al._ (1997a, b) at CERN,
and from experiments at SLAC Whitlow  _et al._ (1992), as well as from reduced
electron and positron cross sections from the H1 and ZEUS Collaborations
Abramowicz  _et al._ (2015) at DESY. These include both proton Benvenuti  _et
al._ (1989); Whitlow  _et al._ (1992); Arneodo  _et al._ (1997a) and deuteron
Benvenuti  _et al._ (1990); Whitlow  _et al._ (1992); Arneodo  _et al._
(1997b) targets, and with both neutral and charged current probes Abramowicz
_et al._ (2015). For the kinematics we implement cuts of $W^{2}>10$ GeV2 and
$Q^{2}>m_{c}^{2}$, where $W^{2}=M^{2}+Q^{2}(1-x_{\mbox{\tiny\rm
Bj}})/x_{\mbox{\tiny\rm Bj}}$, in order to select DIS data that can be fitted
within leading power factorization.

For Drell-Yan lepton-pair production data we use differential cross section
measurements $\mathrm{d}^{2}\sigma_{\mbox{\tiny
DY}}/\mathrm{d}Q\mathrm{d}x_{\mbox{\tiny\rm F}}$ by the E866/NuSea
Collaboration Hawker  _et al._ (1998); Towell  _et al._ (2001); Webb (2003) at
Fermilab, which include proton scattering from proton and deuteron targets. We
include data in the range $Q^{2}>36$ GeV2. Excluding lower $Q^{2}$ data is
recommended by Ref. Alekhin  _et al._ (2006), which demonstrated that
inclusion of the lower $Q^{2}$ data results in deteriorated prediction quality
with no reduction in uncertainty when compared with fits to DIS data alone.

Figure 2: Kinematic coverage of data used in this analysis, with $Q^{2}$
versus the Bjorken scaling variable $x_{\mbox{\tiny\rm Bj}}$ for inclusive DIS
Benvenuti  _et al._ (1989, 1990); Whitlow  _et al._ (1992); Arneodo  _et al._
(1997a, b); Abramowicz  _et al._ (2015) and SIDIS data Adolph  _et al._
(2017a, b) (left panel), fragmentation variable $z$ for SIDIS and SIA data
Brandelik  _et al._ (1980); Althoff  _et al._ (1983); Braunschweig  _et al._
(1989); Lu (1986); Aihara  _et al._ (1984, 1988); Cowan (1988); Derrick  _et
al._ (1987); Itoh  _et al._ (1995); Abe  _et al._ (2004); Buskulic  _et al._
(1995); Akers  _et al._ (1994); Abreu  _et al._ (1998); Albrecht  _et al._
(1989); Leitgab  _et al._ (2013); Leitgab (2013); Lees  _et al._ (2013)
(central panel), and momentum fractions $x_{1},x_{2},x_{\mbox{\tiny\rm F}}$
for Drell-Yan data Hawker  _et al._ (1998); Towell  _et al._ (2001); Webb
(2003) (right panel).

All SIA measurements are of the normalized differential cross sections
$(\mathrm{d}\sigma_{\mbox{\tiny SIA}}/\mathrm{d}z_{h})/\sigma_{\rm tot}$ for
the reaction $e^{+}e^{-}\to(\pi^{\pm},K^{\pm},h^{\pm})\,X$. The data are from
experiments performed by the TASSO Brandelik  _et al._ (1980); Althoff  _et
al._ (1983); Braunschweig  _et al._ (1989) and ARGUS Albrecht  _et al._ (1989)
Collaborations at DESY, by the TPC Lu (1986); Aihara  _et al._ (1984, 1988);
Cowan (1988), HRS Derrick  _et al._ (1987), SLD Abe  _et al._ (2004) and BaBar
Lees  _et al._ (2013) Collaborations at SLAC, by the OPAL Akers  _et al._
(1994); Abbiendi  _et al._ (2000), ALEPH Buskulic  _et al._ (1995) and DELPHI
Abreu  _et al._ (1998) Collaborations at CERN, and by the TOPAZ Itoh  _et al._
(1995) and Belle Leitgab  _et al._ (2013); Leitgab (2013) Collaborations at
KEK. As shown in Fig. 2, the SIA data cover the large-$Q^{2}$ region where a
leading power description in terms of FFs should be accurate. Approximately
half of the SIA data points have $Q\approx M_{Z}$, while the Belle and BaBar
$B$ factories have lower $Q\approx 10.5$ GeV. To ensure applicability of the
leading power formalism, the SIA data in our fits are restricted to the range
$0.2<z_{h}<0.9$.

Identification of heavy quark flavors for some of the SIA datasets is achieved
through measurement of the total energy and momentum in secondary vertices.
The flavor tagged cross sections for a specific flavor $q=c$ or $b$ are
particularly sensitive to the $D_{q}^{h}$, $D_{\bar{q}}^{h}$ and $D_{g}^{h}$
fragmentation functions into the observed hadron $h$. In general, however,
care needs to be taken with the precise method for separating primary quark
flavors, and there are ongoing discussions regarding the optimal approach to
this. For more in-depth discussion see, for example, Ref. Sato  _et al._
(2016).

Finally, the critical addition in this work compared with the previous JAM19
analysis Sato  _et al._ (2020) is the inclusion of unidentified charged hadron
data, along with charged pions and kaons, in the SIDIS off deuterium targets
from the COMPASS Collaboration Adolph  _et al._ (2017a, b) at CERN. Since the
SIDIS data $\mathrm{d}\sigma_{\mbox{\tiny
SIDIS}}^{h^{\pm}}/\mathrm{d}Q^{2}\mathrm{d}x_{\mbox{\tiny\rm
Bj}}\mathrm{d}z_{h}$ are differential in $x_{\mbox{\tiny\rm Bj}}$ and $z_{h}$,
they combine information on both PDFs and FFs, which appear in the description
of SIA, Drell-Yan, and DIS data. Furthermore, as illustrated in Fig. 2, the
SIDIS data have significant overlap in $x_{\mbox{\tiny\rm Bj}}$ and $z_{h}$
with the $x_{\mbox{\tiny\rm Bj}}$ and $x_{\mbox{\tiny\rm F}}$ range of
inclusive DIS and Drell-Yan data, respectively, and the $z_{h}$ range of SIA
data, so that the combined analysis constitutes a genuine test of their
universality. For the COMPASS SIDIS data we use the same kinematic cuts on
$W^{2}$ and $Q^{2}$ as for inclusive DIS, and restrict the fragmentation
variable to $0.2<z_{h}<0.8$ in order to exclude data from the target
fragmentation region and avoid large-$z$ threshold corrections.

##  IV Assessing universality

Before proceeding to the results of our numerical analysis, we briefly discuss
the criteria for universality of the PDFs and FFs and how these are
implemented in our analysis. Extracting parton correlation functions, and
using the extractions to test models of parton structure, is a nontrivial
inverse problem, the detailed examination of which is beyond the scope of the
present paper. However, a claim that the success of a fit is a measure of the
predictive power of the PDFs and FFs requires a number of basic minimal
conditions to be met:

  1. 1.

The system of unknown correlation functions must be over-constrained, by which
we mean that the constraints on unknown correlations imposed by data (or other
theoretical constraints such as sum rules) must be greater than the total
number of functions involved.

  2. 2.

Each unknown correlation function must appear at least twice within the set of
factorization formulas relating the correlation functions to physical
observables.

  3. 3.

There must be reasonable kinematical overlap between the observables so that
correlation functions can be compared within similar ranges of parton momentum
fractions.

Using isospin invariance to relate the PDFs in the proton to those in the
neutron, we have seven independent PDFs: $f_{u}$, $f_{d}$, $f_{s}$,
$f_{\bar{u}}$, $f_{\bar{d}}$, $f_{\bar{s}}$ and $f_{g}$, with PDFs for heavy
flavors generated perturbatively. For the FFs, there are five functions for
$\pi^{+}$ production: $D^{\pi^{+}}_{u}$, $D^{\pi^{+}}_{\bar{u}}$,
$D^{\pi^{+}}_{c}$, $D^{\pi^{+}}_{b}$ and $D^{\pi^{+}}_{g}$, assuming that for
equal $u$ and $d$ quark masses we can equate
$D^{\pi^{+}}_{\bar{d}}=D^{\pi^{+}}_{u}$. Charge symmetry allows all the FFs
for $\pi^{-}$ production to be related to those for $\pi^{+}$ production. For
$K^{+}$ production, there are six independent FFs: $D^{K^{+}}_{u}$,
$D^{K^{+}}_{\bar{s}}$, $D^{K^{+}}_{\bar{u}}$, $D^{K^{+}}_{c}$, $D^{K^{+}}_{b}$
and $D^{K^{+}}_{g}$, where we differentiate between the $u$ and $\bar{s}$
functions. Again, using charge symmetry the FFs for $K^{-}$ can be obtained
from these six $K^{+}$ FFs. Finally, for the unidentified charged hadrons
$h^{\pm}$ or residual FFs, we identify five independent functions: $D^{\rm
res^{+}}_{u}$, $D^{\rm res^{+}}_{\bar{u}}$, $D^{\rm res^{+}}_{c}$, $D^{\rm
res^{+}}_{b}$ and $D^{\rm res^{+}}_{g}$. This makes then a total of 23
functions to be determined.

The quark and gluon PDFs are constrained by their appearance in several sum
rules; in particular, the number sum rules,

$\int_{0}^{1}\mathrm{d}x\,\big{(}f_{q}(x)-f_{\bar{q}}(x)\big{)}=n_{q},$ (28)

where $n_{u}=2$, $n_{d}=1$ and $n_{s}=0$, and the momentum sum rule,

$\sum_{i=q,\bar{q},g}\int_{0}^{1}\mathrm{d}x\,xf_{i}\left(x\right)=1.$ (29)

Note that in Section II.3 these constraints were specifically used to fix the
values of the normalization parameters for several fitted functions. However,
for the purpose of assessing universality, they are simply counted as
additional independent equations which include and thus constrain the PDFs.

The data sets discussed in Sec. III also constrain the light quark and gluon
PDFs since they appear in expressions for multiple independent observables.
Counting these and also the four sum rules (28) and (29),

$f_{i}(x)\stackrel{{\scriptstyle i\neq c,b}}{{\ \longrightarrow\ }}\
\begin{cases}6\qquad\text{DIS}\\\ 2\qquad\text{Drell-Yan}\\\
6\qquad\text{SIDIS}\\\ 4\qquad\text{sum rules}\end{cases}$ (30)

there is a total of 18 relations between the light quark PDFs. The heavy
quarks appear in an even greater number of observables. The light quark
fragmentation functions appear in at least one SIA observable and, because of
charge conjugation invariance, in 2 SIDIS observables,

$D_{i}^{\pi^{+}}(z)\stackrel{{\scriptstyle i\neq c,b,g}}{{\ \longrightarrow\
}}\begin{cases}1\qquad\text{SIA}\\\ 2\qquad\text{SIDIS}\end{cases}$ (31)

and similarly for the kaon and charged hadron fragmentation functions.

For a robust stress-test of universality, there should be reasonable overlap
of the ranges in parton momentum fraction for both the PDFs and the FFs. An
indication for how well this is achieved in the current fit can be be gleaned
from the kinematical coverage plots shown in Fig. 2. To lowest order in
$\alpha_{s}$, the kinematical variables $x_{\mbox{\tiny\rm Bj}}$, $x_{1}$,
$x_{2}$ and $z_{h}$ approximate the parton momentum fractions $x$ and $z$,
respectively, while QCD evolution relates all values of $Q^{2}$. Figure 2
confirms that PDFs and FFs are both constrained by multiple processes in
overlapping regions of momentum fractions.

In summary, our analysis does indeed fulfill the basic criterion for
qualifying as a test of universality, and retaining predictive power for the
PDFs and FFs more generally. Note, however, that the momentum sum rule for FFs
has not been imposed in the analysis. Instead, this will be used as a
consistency check for the final fit in Sec. V.

##  V Numerical analysis

In this section we present the results of our simultaneous Monte Carlo
analysis of PDFs and FFs. We begin with a survey of the fitted cross sections
for the various global datasets used in this study, focusing especially on the
quality of agreement with the SIDIS and SIA data on $\pi^{\pm}$ and $K^{\pm}$,
as well as unidentified $h^{\pm}$ production. We then present our final fitted
PDFs and FFs, and discuss the vital role played by the SIDIS and SIA datasets
in particular in constraining the strange quark distribution in the proton.

###  V.1 Data and theory agreement

Figure 3: Reduced $\chi^{2}_{\rm red}$ values for each DIS (red), DY (green),
SIDIS (orange) and SIA (blue) experiment considered in this analysis (left
column), along with the corresponding mean and standard deviation of the
residuals for each experiment, E [residual] (right column).

To assess the agreement of the fitted results with the various datasets, in
Fig. 3 we show the reduced $\chi^{2}$ for each individual experiment, which is
defined by

$\chi^{2}_{\rm
red}=\frac{1}{N}\sum_{i,e}\frac{1}{\alpha_{i,e}^{2}}\bigg{(}d_{i,e}-{\rm
E}\Big{[}\sum_{k}r_{e}^{k}\beta^{k}_{i,e}+T_{i,e}/N_{e}\Big{]}\bigg{)}^{2}.$
(32)

Here, the expectation value ${\rm E}[...]$, as defined in Eq. (27a),
represents the mean theory, including optimized multiplicative and additive
corrections to match the data, with $N$ the total number of data points. In
Fig. 3 we show the mean and standard deviation of the Monte Carlo residuals
for each experiment $e$, where the residual per data point is defined as

${\rm residual}\,(e,i)=\frac{1}{\alpha_{i,e}}\bigg{(}d_{i,e}-{\rm
E}\Big{[}\sum_{k}r_{e}^{k}\beta^{k}_{i,e}+T_{i,e}/N_{e}\Big{]}\bigg{)}.$ (33)

For the inclusive DIS, Drell-Yan and SIDIS datasets we find excellent overall
agreement between data and theory, with $\chi^{2}_{\rm red}$ values close to
1. The $\chi^{2}_{\rm red}$ for the SIA datasets are slightly higher, but
nonetheless the overall fit is very good, giving a total reduced
$\chi^{2}_{\rm red}$ = 1.26 for almost 5000 data points. The values of
$\chi^{2}_{\rm red}$ for each type of dataset and for each specific hadron in
the final state are summarized in Table 1, along with the number of data
points for each dataset.

Table 1: Reduced $\chi^{2}_{\rm red}$ values for each type of dataset (DIS,
Drell-Yan, SIDIS, SIA) considered in this analysis, together with the number
of data points $N_{\rm dat}$ for each dataset.
reaction |  $\chi^{2}_{\rm red}$ |  $N_{\rm dat}$
---|---|---
DIS |  | 1.29 | 2680
DY |  | 1.52 | 250
SIDIS | $\pi^{\pm}$ | 1.39 | 498
| $K^{\pm}$ | 1.38 | 494
| $h^{\pm}$ | 0.85 | 498
SIA | $\pi^{\pm}$ | 1.09 | 231
| $K^{\pm}$ | 1.37 | 213
| $h^{\pm}$ | 1.15 | 120
total |  | 1.26 | 4984
Figure 4: Comparison of the multiplicities $dM^{h}/dz_{h}$ for $h=\pi^{+}$
(dashed lines) and $\pi^{-}$ (dotted lines) production with the COMPASS data
Adolph  _et al._ (2017a, b) in various bins of $x_{\mbox{\tiny\rm Bj}}$ and
$y$ (offset by a factor $2^{i}$). Figure 5: As in Fig. 4, but for $K^{\pm}$
COMPASS SIDIS data Adolph  _et al._ (2017a, b). Figure 6: As in Fig. 4, but
for unidentified hadron $h^{\pm}$ COMPASS SIDIS data Adolph  _et al._ (2017a,
b).

The residuals profile for the DIS, Drell-Yan and SIDIS datasets is well
centered around zero, with variances $\sim 1$, indicating an average Gaussian
behavior of their associated likelihood function. The variance for the SIDIS
$h^{-}$ data from COMPASS, however, is found to be up to $\approx 50\%$ below
unity, suggesting a deviation from a Gaussian likelihood. This may be due to
the fact that these data are dominated by systematic uncertainties, which is
also reflected by the relatively small reduced $\chi^{2}_{\rm red}$ values,
especially for the COMPASS $h^{-}$ data relative to the rest of DIS and SIDIS
data sets.

A more detailed comparison with the COMPASS SIDIS is made in Figs. 4, 5 and 6,
where we show the $z_{h}$ dependence of the $\pi^{\pm}$, $K^{\pm}$ and
$h^{\pm}$ multiplicities, respectively, which are defined as ratios of SIDIS
to inclusive DIS cross sections at the same $x_{\mbox{\tiny\rm Bj}}$ and
$Q^{2}$,

$\frac{dM^{h}}{dz_{h}}=\frac{\mathrm{d}\sigma_{\mbox{\tiny
SIDIS}}^{h}/\mathrm{d}Q^{2}\mathrm{d}x_{\mbox{\tiny\rm
Bj}}\mathrm{d}z_{h}}{\mathrm{d}\sigma_{\mbox{\tiny
DIS}}/\mathrm{d}Q^{2}\mathrm{d}x_{\mbox{\tiny\rm Bj}}}.$ (34)

The agreement between theory and the experimental $z_{h}$ spectrum is quite
remarkable, given that it spans some 2 orders of magnitude, which suggests
that at these kinematics a leading power perturbative QCD factorization at
next-to-leading order provides sufficient accuracy to describe the data.
Interestingly, the differences between the multiplicities for positively and
negatively charged hadron species increase with $x_{\mbox{\tiny\rm Bj}}$,
especially for kaons, and in the valence region these can differ by an order
of magnitude for low values of $Q^{2}$. Such differences can enhance our
ability to extract flavor dependent effects in nonperturbative PDFs and parton
to kaon FFs from the data. The new data set included for the first time in the
present JAM analysis, namely the unidentified charged hadron data shown in
Fig. 6, are also well described by our nonperturbative ansatz for the
corresponding FFs. In contrast to the excellent agreement with the $z_{h}$
dependence of the data in Figs. 4–6, we note that analysis of the same data
differential in the hadron transverse momentum using existing PDFs and FFs
within TMD factorization results in poor agreement between predictions and
data Gonzalez-Hernandez  _et al._ (2018); Wang  _et al._ (2019), indicating
that further work is needed to understand the SIDIS transverse momentum
spectra.

Figure 7: Data to theory ratios for SIA $\pi^{\pm}$ production cross sections
versus $z_{h}$, with the bands indicating the uncertainty on the fitted
result. Figure 8: As in Fig. 7, but for SIA $K^{\pm}$ production. Figure 9: As
in Fig. 7, but for SIA unidentified charged hadron $h^{\pm}$ production.

For the SIA data sets, there is a somewhat wider spread in the data versus
theory comparisons, as seen in Figs. 7, 8 and 9 for the $\pi^{\pm}$, $K^{\pm}$
and unidentified charged hadron $h^{\pm}$ final states, respectively.
Generally, the $\pi^{\pm}$ data have the best agreement among the SIA
datasets, with a reduced $\chi^{2}_{\rm red}=1.09$, followed by the hadron
data with $\chi^{2}_{\rm red}=1.15$, and lastly the kaon data, which have an
overall reduced $\chi^{2}_{\rm red}=1.37$. For about 3/4 of the $\approx 40$
SIA datasets, we find very good agreement with the global fit, with
$\chi^{2}_{\rm red}\approx 1$ or below. For the remaining datasets that have
larger $\chi^{2}_{\rm red}$ values, to better understand the reasons for some
of the tensions between data and theory we discuss in the following some
individual cases ranked by the reduced $\chi^{2}_{\rm red}$ values.

Starting with the datasets that have the largest $\chi^{2}_{\rm red}$ values,
namely, $\chi^{2}_{\rm red}\gtrsim 3$, we identify the OPAL ($\pi^{\pm}$ and
$c\to K^{\pm}$), TPC ($K^{\pm}$), SLD ($\pi^{\pm}$ and $c\to K^{\pm}$), DELPHI
($K^{\pm}$), and TASSO ($\pi^{\pm}$ and $h^{\pm}$ at 35 GeV) datasets. For the
inclusive OPAL ($\pi^{\pm}$) data, we observe in Fig. 7 that for $z_{h}<0.5$
the data are indeed in tension with the corresponding inclusive ALEPH and SLD
results, and the overall trend of the data/theory ratio suggests a possible
normalization issue with this dataset. This can also be said for the DELPHI
($K^{\pm}$) which appears to have some tension with the corresponding
inclusive OPAL and ALEPH results. Similarly, from Fig. 8 we find that the TPC
($K^{\pm}$) spectrum lies below the theory, suggesting again a normalization
problem with these data. The situation for the TASSO ($\pi^{\pm}$) data is
less clear, as only the $Q=14$ GeV dataset seems to give a bad fit, while data
at other energies can be described fairly well. This again hints at a problem
with the overall normalization for this dataset. The same behavior appears
also in the TASSO ($h^{\pm}$) data in Fig. 9, where both the $Q=35$ and 45 GeV
datasets are above the theoretical cross sections. The case of SLD and OPAL
($c\to K^{\pm}$) data in Fig. 8 shows a clear overestimation of the $z_{h}$
spectra. While one can argue that this problem could be a reflection of the
need for a more sophisticated heavy quark treatment in our theory, the
description of $b$-tagged data from SLD, DELPHI and OPAL is relatively good,
so that an explanation in terms of a normalization uncertainty in the SLD and
OPAL ($c\to K^{\pm}$) data may be more relevant.

For SIA datasets that have smaller, but still large, $\chi^{2}_{\rm red}$
values, $2\lesssim\chi^{2}_{\rm red}\lesssim 3$, we identify the $b$-tagged
TPC ($b\to\pi^{\pm}$), OPAL ($b\to\pi^{\pm}$), and TASSO ($h^{\pm}$ at 44
GeV). For the case of the TPC ($b\to\pi^{\pm}$) data, we see from Fig. 7 that
for the largest $z_{h}$ bin the theory overestimates the data. On the other
hand, good agreement is found for the SLD ($b\to\pi^{\pm}$) data at the same
kinematics. It is possible that at the smaller $Q$ values of TPC relative to
SLD, the range in $z_{h}$ where leading power factorization is applicable is
narrower, in particular for the $b$-tagged data. The $z_{h}$ dependence of the
OPAL ($b\to\pi^{\pm}$) data appear to be clearly different from the theory,
even within the large uncertainties. We note here that the OPAL data are
presented as truncated moments as a function of the lower limit of the
integration, $z_{h}^{\rm min}$, and the inclusion of the very high $z_{h}$
bins may be problematic for the validity of factorization theorems at
$z_{h}\to 1$. Lastly, as with TASSO ($h^{\pm}$ at 35 GeV), the somewhat large
$\chi^{2}_{\rm red}$ values for the 44 GeV data is likely attributable to a
problem with overall normalization.

For datasets that have $\chi^{2}_{\rm red}\lesssim 2$, we consider the
agreement to be generally acceptable. Indeed, the vast majority of datasets in
this category have $\chi^{2}_{\rm red}\approx 1$ or below. These include all
of the recent high-statistics $B$-factory data from BaBar ($\pi^{\pm}$,
$K^{\pm}$) and Belle ($\pi^{\pm}$, $K^{\pm}$), most of the TASSO ($\pi^{\pm}$,
$K^{\pm}$), TPC ($\pi^{\pm}$, $c\to\pi^{\pm}$) and SLD ($h^{\pm}$,
$b\to\pi^{\pm}$, $b\to h^{\pm}$) datasets, all of the ALEPH ($\pi^{\pm}$,
$K^{\pm}$, $h^{\pm}$) and most of the DELPHI ($\pi^{\pm}$, $K^{\pm}$, $b\to
K^{\pm}$, $h^{\pm}$, $b\to h^{\pm}$) data, along with the older ARGUS
($K^{\pm}$), TOPAZ ($\pi^{\pm}$, $K^{\pm}$) and OPAL ($K^{\pm}$, $h^{\pm}$,
$c\to h^{\pm}$, $b\to h^{\pm}$) data. Slightly higher, but still reasonable,
$\chi^{2}_{\rm red}$ values are obtained for the ARGUS ($\pi^{\pm}$), TPC
($h^{\pm}$), DELPHI ($b\to\pi^{\pm}$), and SLD ($K^{\pm}$, $c\to\pi^{\pm}$,
$c\to h^{\pm}$) datasets.

Finally, we note that most of the large $\chi^{2}_{\rm red}$ values found in
this analysis were absent in the previous JAM Monte Carlo analysis of
fragmentation functions Sato  _et al._ (2016). The main reason is the
restriction of the SIA datasets here to the range $0.2<z_{h}<0.8$, chosen to
coincide with the range over which the SIDIS data in this work are able to be
described within collinear factorization. For the LEP data in particular there
are many data points at $z_{h}<0.2$ which can be well fitted within the
current framework, and which would reduce the overall $\chi^{2}_{\rm red}$. A
careful point by point comparison of the individual $\chi^{2}_{\rm red}$
values for the various datasets indeed confirms that similar discrepancies
also occurred in Ref. Sato  _et al._ (2016). However, for consistency in our
joint analysis of PDFs and FFs, we restrict the kinematic range to the region
where both SIA and SIDIS can be simultaneously described. The same choice for
the $z_{h}$ range was made in the recent JAM19 analysis, which required SIDIS
data to be restricted to $z_{h}\gtrsim 0.2$ to ensure separation of the target
and current fragmentation regions.

###  V.2 Parton distributions and fragmentation functions

The proton PDFs from our simultaneous fit are displayed in Fig. 10 at a scale
$\mu^{2}=10$ GeV2, where we focus on the kinematic region of parton momentum
fractions $x\gtrsim 0.01$ that is constrained by the SIDIS data. For
comparison, we contrast our results with other next-to-leading order PDF
parametrizations, namely, from the CJ15 Accardi  _et al._ (2016) and NNPDF3.1
Ball  _et al._ (2017) global analyses. Compared with the other fits, our
valence $u$ and $d$ quark distributions have slightly larger magnitude in the
intermediate-$x$ region, $x\sim 0.1$, with a compensating stronger suppression
at small $x$ needed to ensure that the valence number sum rules are respected.
The ratio $d/u$ is quite compatible with the results from the other groups, on
the other hand, but has a significantly larger uncertainty at large $x$
compared with the CJ15 result, reflecting the Monte Carlo nature of our
analysis.

Figure 10: Proton PDFs from the present JAM20-SIDIS analysis (red bands)
versus $x$ at a scale $\mu=10$ GeV2, compared with the CJ15 Accardi  _et al._
(2016) (blue bands) and NNPDF3.1 Ball  _et al._ (2017) (green bands)
parametrizations. The bands shown represent the mean $\pm 1\sigma$.

The intermediate-$x$ enhancement in the valence PDFs in our fit is correlated
with the slightly smaller $\bar{d}+\bar{u}$ light antiquark sea compared with
the CJ15 and NNPDF3.1 parametrizations. This in turn is correlated with the
behavior of the strange quark sea, as seen in the ratio

$R_{s}=\frac{s+\bar{s}}{\bar{d}+\bar{u}}$ (35)

of the strange to nonstrange sea quark PDFs. In Fig. 10 this ratio is
generally larger in our analysis than for the other parametrizations, with a
somewhat bigger uncertainty. This is understood from the fact that in the CJ15
fit $R_{s}$ is fixed to be 0.4 at the input scale, with deviations from the
constant value arising only from DGLAP evolution. For the NNPDF3.1 fit the
uncertainties are smaller because of their inclusion of the neutrino DIS data,
which we do not include in our analysis because of unknown nuclear corrections
in neutrino scattering Accardi  _et al._ (2010); Kovarik  _et al._ (2011);
Kalantarians  _et al._ (2017). Our light antiquark asymmetry $\bar{d}-\bar{u}$
is also compatible with the other groups, but again with a larger uncertainty,
which may be related to the absence of collider $W$ and lepton asymmetry data
in our fit. Finally, for the gluon distribution, the magnitude and
uncertainties are very similar across all the analyses, even though our fit
does not include jet production data from hadron colliders. This reflects the
fact that the HERA DIS data, which are included here, provide strong
constraints on the shape of the gluon PDF via scaling violations.

Figure 11: Parton to hadron FFs versus $z$ at $\mu^{2}=100$ GeV2 from the
JAM20-SIDIS analysis for various parton flavors fragmenting to $\pi^{+}$ (red
bands), $K^{+}$ (blue bands), unidentified hadrons $h^{+}$ (green bands), and
residual hadrons $\delta h^{+}$ (yellow bands), defined as the difference
between $h^{+}$ and the sum of $\pi^{+}$ and $K^{+}$. The bands shown are mean
$\pm 1\sigma$.

For the parton to hadron FFs, we show in Fig. 11 the $z$ dependence of the FFs
at a scale $\mu^{2}=100$ GeV2 for the positively charged $\pi^{+}$, $K^{+}$
and unidentified hadrons $h^{+}$, as well as for the residual hadrons $\delta
h^{+}$, defined as the difference between $h^{+}$ and the sum of $\pi^{+}$ and
$K^{+}$ (so that the total is given by $h^{+}=\pi^{+}+K^{+}+\delta h^{+}$).
For most of the flavors we find that the quark $\to$ $\pi^{+}$ fragmentation
dominates, as expected from the pion being the lightest hadron in the QCD
spectrum. Exceptions to this are for $\bar{s}\to K^{+}$ and $c\to K^{+}$ at
intermediate $z$ values, and for $b$ quark fragmentation into residual hadrons
$\delta h^{+}$.

For gluon fragmentation, pion production dominates for $z$ up to $\sim
0.5-0.6$, above which kaon fragmentation becomes as sizeable as the pion. This
is consistent with the findings of previous FF analyses Hirai  _et al._
(2007); Sato  _et al._ (2016), which observed that the production of heavier
particles such as kaons requires larger momentum fractions from the
fragmenting gluon compared to the production of lighter particles.

The production of hadrons heavier than kaons, as indicated in Fig. 11 by the
residual hadrons $\delta h^{+}$, can be sizable and comparable to that of
kaons, especially for the $d$ and $s$ quarks and at large values of $z$. The
relatively large $d\to\delta h^{+}$ FF can be understood in terms of the
fragmentation into protons. Note that we have imposed flavor symmetry for the
residual hadron fragmentation, so that $D_{d}^{\delta h^{+}}=D_{s}^{\delta
h^{+}}$. In principle, the presence of hyperons such as $\Sigma^{+}$ should
brake this relation, but we leave analysis of such effects for future work. As
the case for the $g\to K^{+}$, the fragmentation of gluons into heavier
particles peaks at large $z$, where larger momentum fractions from the
fragmenting gluons are need for the production of heavier particles.

For production of hadrons initiated by heavy quarks, we find similar
fragmentation of charm quarks into pions and kaons, but a rather different
pattern for the fragmentation of bottom quarks. Some of this difference can be
explained by the flavor-changing properties of $u$-type quarks decaying into
$d$-type quarks. While the charm quark can decay into strange quarks and hence
enhance $K^{+}$ production, the same does not occur for bottom quarks, which
suppresses kaon production relative to pion production due to the mass
difference. Interestingly, the production of other species of charged hadrons
is much larger for $b$ quarks than for $c$ quarks, which may be understood
from the greater phase space available for $b$ quarks to decay into heavier
hadrons to which charm quarks cannot transition.

Figure 12: Normalized yield of truncated moments $\langle z\rangle_{i}^{h}$ of
the $i\to h$ FFs $zD_{i}^{h}$, for the favored $\pi^{+}$ (red) and favored
$K^{+}$ (blue), unfavored $\pi^{+}$ (light red) and unfavored $K^{+}$ (light
blue), the total hadron $h^{+}$ (green) and residual hadron $\delta h^{+}$
(yellow) FFs, at a scale of $\mu^{2}=100$ GeV2.

In Fig. 12 we present truncated moments

$\langle z\rangle_{i}^{h}=\int_{z_{\rm min}}^{1}\,dz\,zD_{i}^{h}(z),$ (36)

for each flavor $i$ and final state hadron $h$, where we take the lower limits
on the $z$ integration $z_{\rm min}=0.2$ to restrict the moment to the region
of SIDIS kinematics. The truncated moment indicate how energetic is the
production different a hadron species $h$ relative to the parent parton $i$.
In general, we find that the production of hadron species heavier than pions
and kaons is typically produced with lower energies, which is consistent with
the physical picture whereby more energy is required to produce heavier
hadrons than lighter hadrons.

As expected, the favored fragmentation of $\bar{d}$ quarks is predominantly
into highly energetic pions, while for the antistrange $\bar{s}$ the
production rate of energetic kaons is slightly higher than that of pions. The
unfavored fragmentation of $d,s$ and $\bar{u}$ quarks follows a similar
pattern, with the lightest (pion) state produced at the highest energies
followed by kaons and other heavier charged hadrons. An exception to this
behavior is for charm and bottom quark fragmentation: for $c$ quarks kaons are
produced with energies comparable to those of pions, while for $b$ quarks kaon
production is suppressed with heavier mass hadrons produced at similar
energies as pions.

Interestingly, the production of hadrons from gluons follows the same pattern
as for $u$-quark fragmentation. While the latter can be explained in terms of
mass differences between the produced hadron species, the fact that $u$ quarks
and gluons give a similar average energy profile across hadron species is
intriguing. On perturbative grounds one can argue that gluon fragmentation is
enhanced because of the $C_{A}=3$ factor in the the gluon splitting function,
$P_{gg}$, relative to quark splitting functions, $P_{qq}$ and $P_{gq}$, which
are proportional to $C_{F}=4/3$. The absence of direct constraints on the
gluon FF beyond scaling violations, however, anything drawing more than
speculative conclusions at present.

Figure 13: Monte Carlo samples for the $R_{s}$ ratio (left) and
$zD_{\bar{s}}^{K^{+}}$ FF (right) at $\mu^{2}=10$ GeV2, color coded according
to the scaled $\chi^{2}_{\rm red}$ for the SIA ($K^{\pm}$) (top row) and SIDIS
($K^{+},K^{-}$) (bottom row) datasets.

We conclude the discussion of our numerical results by focusing on the
correlation between the strange to nonstrange PDF ratio $R_{s}$ and the
strange to kaon fragmentation function $D_{\bar{s}}^{K^{+}}$. In Fig. 13 we
show $R_{s}$ and the $\bar{s}\to K^{+}$ FF, with individual Monte Carlo
samples color coded by the scaled $\chi^{2}_{\rm red}$ intensity (with darker
replicas indicating higher likelihoods) computed for the specific cases of SIA
($K^{\pm}$) and SIDIS ($K^{+},K^{-}$) datasets. The SIA datasets have a clear
preference for a smaller $R_{s}$ and enhanced $D_{\bar{s}}^{K^{+}}$, as was
found in the previous JAM19 analysis Sato  _et al._ (2020). Interestingly, the
SIDIS ($K^{+},K^{-}$) data, which have smaller $\chi^{2}_{\rm red}$, have a
slight tendency to favor solutions with a larger $R_{s}$ and smaller
$D_{\bar{s}}^{K^{+}}$, however, this preference is much weaker than the
preference of the SIA data for smaller $R_{s}$ values.

We also note that in the current analysis we have extended the flexibility of
the PDF and FF parametrizations, which allowed us to obtain a more uniform
Monte Carlo distribution of $R_{s}$ compared JAM19, where a more restricted
parametrization gave rise to multiple solutions. Our new analysis confirms
that the most probable solutions found in JAM19 did not result from
parametrization bias, and corroborates the need for a suppressed strange quark
PDF in the proton in order to simultaneously describe both the SIA and SIDIS
datasets within leading power QCD factorization.

##  VI Conclusion

In this paper we have presented the results of a simultaneous Monte Carlo
analysis of PDFs and FFs constrained by a diverse array of data from inclusive
and semi-inclusive DIS, Drell-Yan lepton-pair production, and SIA in
$e^{+}e^{-}$ collisions. The analysis extends the previous JAM19 Sato  _et
al._ (2020) simultaneous fit by including in addition unidentified charged
hadrons in the final states of SIDIS and SIA, and increasing the flexibility
of the PDF and FF parametrizations.

The analysis — referred to as “JAM20-SIDIS” — represents the most
comprehensive determination of parton to hadron ($\pi^{\pm}$, $K^{\pm}$,
$h^{\pm}$) FFs fitted concurrently with spin-averaged parton distributions,
broadening the test of universality of parton correlation functions to more
observables. The more thorough exploration of the parameter space and reduced
$\chi^{2}_{\rm red}$ values for each of the $\approx 70$ datasets fitted in
this study confirmed the previous finding Sato  _et al._ (2020) that the
combination of SIA and SIDIS datasets have a strong preference for a smaller
strange to nonstrange PDF ratio, $R_{s}$, correlated with an enhanced
$D_{\bar{s}}^{K^{+}}$ FF. As further tests of this scenario, we plan in future
to extend the experimental datasets to include weak-boson and jet production
in hadronic collisions, from both Tevatron and LHC data, as well as to relax
the $W^{2}$ cuts for inclusive DIS to incorporate more fixed-target DIS data
at high $x_{\mbox{\tiny\rm Bj}}$ values Coc .

An important application of the current results will be in benchmark
calculations of transverse momentum dependent cross sections, and in
particular for the small transverse momentum region where the transition from
collinear factorization to TMD factorization is expected to set in. One
motivation for the present project was to assess the possible role of
limitations in collinear PDF and FF fits in explaining discrepancies between
theory and data in the range of intermediate and large transverse momentum
across a number of transversely differential processes Gonzalez-Hernandez  _et
al._ (2018); Wang  _et al._ (2019); Bacchetta  _et al._ (2019); Moffat  _et
al._ (2019). For this, a truly simultaneous analysis of parton distribution
and fragmentation functions across the standard set of electromagnetic
processes, integrated over all transverse momentum, is necessary. The general
success of the collinear fits for transverse momentum integrated cross
sections that we have examined here, and their evident predictive power,
suggests that factors unique to the transverse momentum differential treatment
are responsible for the tension with data. We plan to address this also in
future work.

###### Acknowledgements.

This work was supported by the US Department of Energy contract DE-
AC05-06OR23177, under which Jefferson Science Associates, LLC operates
Jefferson Lab. The work of N.S. was supported by the DOE, Office of Science,
Office of Nuclear Physics in the Early Career Program. E.M. and T.R. were
supported by the U.S. Department of Energy, Office of Science, Office of
Nuclear Physics, under Award Number DE-SC0018106.

##  Appendix A Semi-inclusive DIS

In this appendix we summarize the basic cross section and structure function
formulas relevant for our analysis of the semi-inclusive leptoproduction of a
hadron $h$ (with four-momentum $p_{h}$) in the deep-inelastic scattering of a
lepton $\ell$ (momentum $l$) from a nucleon $N$ (momentum $p$) via the
exchange of a virtual photon (momentum $q$),

$\ell+N\rightarrow\ell^{\prime}+h+X\,,$ (37)

where the final state hadron is integrated over all transverse momentum. The
formal setup follows standard methods described in Refs. Altarelli  _et al._
(1979); Nason and Webber (1994); Furmanski and Petronzio (1982); Graudenz
(1994), for example, and we closely follow the specific techniques in Ref. de
Florian  _et al._ (1998) utilizing the double Mellin moment method from Ref.
Stratmann and Vogelsang (2001).

The spin-averaged cross section is parametrized in terms of the semi-inclusive
structure functions $F_{1}^{h}$ and $F_{L}^{h}$,

$\frac{\mathrm{d}\sigma{}}{\mathrm{d}x_{\mbox{\tiny\rm
Bj}}\mathrm{d}y\,\mathrm{d}z_{h}}=\frac{4\pi\alpha^{2}}{Q^{2}}\bigg{[}\frac{1+(1-y)^{2}}{y}F_{1}^{h}(x_{\mbox{\tiny\rm
Bj}},z_{h},Q^{2})+\frac{(1-y)}{x_{\mbox{\tiny\rm
Bj}}\,y}\,F_{L}^{h}(x_{\mbox{\tiny\rm Bj}},z_{h},Q^{2})\bigg{]},$ (38)

which are functions of the Bjorken scaling $x_{\mbox{\tiny\rm Bj}}$, the
hadron fragmentation scaling variable $z_{h}$, and the four-momentum transfer
squared $Q^{2}$, defined in the standard way as

$x_{\mbox{\tiny\rm Bj}}=\frac{Q^{2}}{2p\cdot q}\,,\qquad
z_{h}=\frac{p_{h}\cdot p}{q\cdot p}\,,\qquad Q^{2}=-q^{2},$ (39)

with $y=q\cdot p/l\cdot p$ the inelasticity. Our conventions for the semi-
inclusive structure function definitions are directly related to the
conventions for the inclusive DIS structure functions, with $F_{T}=2xF_{1}$
and $F_{L}=F_{2}-2xF_{1}$ the transverse and longitudinal structure functions,
respectively. (Note that other conventions use instead $F_{1}\to 2F_{1}$ and
$F_{L}\to F_{L}/x$ Anderle  _et al._ (2013).)

In the current fragmentation region, the factorization formulas for the semi-
inclusive structure functions, expanded to order $\alpha_{s}$ in the hard
part, are given by

$\displaystyle 2F_{1}^{h}(x_{\mbox{\tiny\rm Bj}},z_{h},Q^{2})$
$\displaystyle=$
$\displaystyle\sum_{q,\bar{q}}e_{q}^{2}\,\bigg{\\{}f_{q}(x_{\mbox{\tiny\rm
Bj}},Q^{2})\,D_{q}^{h}(z_{h},Q^{2})$ (40)
$\displaystyle+\,\frac{\alpha_{s}(Q^{2})}{2\pi}\Big{[}f_{q}\otimes
C_{qq}^{1}\otimes D_{q}^{h}+f_{q}\otimes C_{gq}^{1}\otimes
D_{g}^{h}+f_{g}\otimes C_{qg}^{1}\otimes D_{q}^{h}\Big{]}(x_{\mbox{\tiny\rm
Bj}},z_{h},Q^{2})\bigg{\\}},$
$\displaystyle\frac{1}{x}F_{L}^{h}(x_{\mbox{\tiny\rm Bj}},z_{h},Q^{2})$
$\displaystyle=$ $\displaystyle\sum_{q,\bar{q}}e_{q}^{2}\,$ (41)
$\displaystyle\times\frac{\alpha_{s}(Q^{2})}{2\pi}\Big{[}f_{q}\otimes
C_{qq}^{L}\otimes D_{q}^{h}+f_{q}\otimes C_{gq}^{L}\otimes
D_{g}^{h}+f_{g}\otimes C_{qg}^{L}\otimes D_{q}^{h}\Big{]}(x_{\mbox{\tiny\rm
Bj}},z_{h},Q^{2}),$

where $f_{i}$ and $D_{j}^{h}$ label the PDF of flavor $i$ in the proton and
parton $j$ $\to$ hadron $h$ FF, respectively. The functions $C^{1}_{ij}$
($C^{L}_{ij}$) are the lowest-order hard scattering coefficient functions for
the $F_{1}$ ($F_{L}$) structure functions, and the symbol $\otimes$ denotes
the convolution integral over longitudinal momentum fractions,
$\left[A\otimes{B}\right](x)\equiv\int_{x}^{1}\,(\mathrm{d}z/z)\,A(z)B(x/z)$.
Explicit expressions for $C^{1}_{ij}$ ($C^{L}_{ij}$) are given, for example,
in the appendices of Ref. de Florian  _et al._ (1998).

In all of the expressions above, kinematical corrections from a nonzero target
and final state hadron mass are neglected. To focus attention on the current
fragmentation region, in our analysis we impose the kinematic cuts $z_{h}>0.2$
and $W^{2}>10$ GeV2, and for the hard scale we choose $Q^{2}>m_{c}^{2}$.
Therefore, the kinematical corrections may not be entirely negligible Boglione
_et al._ (2019) at the energies of some of the experiments, although for now
we set aside a fuller account of their effect to a future dedicated analysis
of mass corrections in SIDIS.

## References

  * Ethier  _et al._ (2017) J. J. Ethier, N. Sato, and W. Melnitchouk, Phys. Rev. Lett. 119, 132001 (2017), arXiv:1705.05889 [hep-ph] .
  * Sato  _et al._ (2020) N. Sato, C. Andres, J. Ethier, and W. Melnitchouk, Phys. Rev. D 101, 074020 (2020), arXiv:1905.03788 [hep-ph] .
  * Cammarota  _et al._ (2020) J. Cammarota, L. Gamberg, Z.-B. Kang, J. A. Miller, D. Pitonyak, A. Prokudin, T. C. Rogers, and N. Sato, Phys. Rev. D 102, 054002 (2020), arXiv:2002.08384 [hep-ph] .
  * Gonzalez-Hernandez  _et al._ (2018) J. O. Gonzalez-Hernandez, T. C. Rogers, N. Sato, and B. Wang, Phys. Rev. D 98, 114005 (2018), arXiv:1808.04396 [hep-ph] .
  * Wang  _et al._ (2019) B. Wang, J. Gonzalez-Hernandez, T. Rogers, and N. Sato, Phys. Rev. D 99, 094029 (2019), arXiv:1903.01529 [hep-ph] .
  * Bacchetta  _et al._ (2019) A. Bacchetta, G. Bozzi, M. Lambertsen, F. Piacenza, J. Steiglechner, and W. Vogelsang, Phys. Rev. D 100, 014018 (2019), arXiv:1901.06916 [hep-ph] .
  * Moffat  _et al._ (2019) E. Moffat, T. Rogers, N. Sato, and A. Signori, Phys. Rev. D 100, 094014 (2019), arXiv:1909.02951 [hep-ph] .
  * (8) T. Liu and J.-W. Qiu, arXiv:1907.06136 [hep-ph] .
  * Collins  _et al._ (1989) J. C. Collins, D. E. Soper, and G. F. Sterman, “Factorization of Hard Processes in QCD,” (1989) pp. 1–91, arXiv:hep-ph/0409313 .
  * Ellis  _et al._ (1996) R. K. Ellis, W. J. Stirling, and B. R. Webber, Camb. Monogr. Part. Phys. Nucl. Phys. Cosmol. 8, 1 (1996).
  * Collins (2011) J. C. Collins,  _Foundations of Perturbative QCD_ (Cambridge University Press, Cambridge, 2011).
  * Devenish and Cooper-Sarkar (2004) R. Devenish and A. Cooper-Sarkar,  _Deep inelastic scattering_ (Oxford University Press, Oxford, 2004).
  * Vogt (2005) A. Vogt, Comput. Phys. Commun. 170, 65 (2005), arXiv:hep-ph/0408244 .
  * Stratmann and Vogelsang (2001) M. Stratmann and W. Vogelsang, Phys. Rev. D 64, 114007 (2001), arXiv:hep-ph/0107064 .
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Floratos  _et al._ (1981) E. Floratos, C. Kounnas, and R. Lacaze, Nucl. Phys. B192, 417 (1981).
  * Kretzer (2000) S. Kretzer, Phys. Rev. D 62, 054001 (2000), arXiv:hep-ph/0003177 .
  * Anastasiou  _et al._ (2003) C. Anastasiou, L. J. Dixon, K. Melnikov, and F. Petriello, Phys. Rev. Lett. 91, 182002 (2003), arXiv:hep-ph/0306192 .
  * Sato  _et al._ (2016) N. Sato, J. Ethier, W. Melnitchouk, M. Hirai, S. Kumano, and A. Accardi, Phys. Rev. D 94, 114004 (2016), arXiv:1609.00899 [hep-ph] .
  * Bertone  _et al._ (2018) V. Bertone, N. P. Hartland, E. R. Nocera, J. Rojo, and L. Rottoli, Eur. Phys. J. C 78, 651 (2018), arXiv:1807.03310 [hep-ph] .
  * de Florian  _et al._ (2007) D. de Florian, R. Sassot, and M. Stratmann, Phys. Rev. D 76, 074033 (2007), arXiv:0707.1506 [hep-ph] .
  * Benvenuti  _et al._ (1989) A. Benvenuti _et al._ , Phys. Lett. B 223, 485 (1989).
  * Benvenuti  _et al._ (1990) A. Benvenuti _et al._ , Phys. Lett. B 237, 592 (1990).
  * Arneodo  _et al._ (1997a) M. Arneodo _et al._ , Nucl. Phys. B483, 3 (1997a), arXiv:hep-ph/9610231 .
  * Arneodo  _et al._ (1997b) M. Arneodo _et al._ , Nucl. Phys. B487, 3 (1997b), arXiv:hep-ex/9611022 .
  * Whitlow  _et al._ (1992) L. W. Whitlow, E. M. Riordan, S. Dasu, S. Rock, and A. Bodek, Phys. Lett. B 282, 475 (1992).
  * Abramowicz  _et al._ (2015) H. Abramowicz _et al._ , Eur. Phys. J. C 75, 580 (2015), arXiv:1506.06042 [hep-ex] .
  * Hawker  _et al._ (1998) E. Hawker _et al._ , Phys. Rev. Lett. 80, 3715 (1998), arXiv:hep-ex/9803011 .
  * Towell  _et al._ (2001) R. Towell _et al._ , Phys. Rev. D 64, 052002 (2001), arXiv:hep-ex/0103030 .
  * Webb (2003) J. C. Webb,  _Measurement of continuum dimuon production in 800-GeV/c proton nucleon collisions_ , Ph.D. thesis, New Mexico State U. (2003), arXiv:hep-ex/0301031 .
  * Alekhin  _et al._ (2006) S. Alekhin, K. Melnikov, and F. Petriello, Phys. Rev. D 74, 054033 (2006), arXiv:hep-ph/0606237 .
  * Adolph  _et al._ (2017a) C. Adolph _et al._ , Phys. Lett. B 764, 1 (2017a), arXiv:1604.02695 [hep-ex] .
  * Adolph  _et al._ (2017b) C. Adolph _et al._ , Phys. Lett. B 767, 133 (2017b), arXiv:1608.06760 [hep-ex] .
  * Brandelik  _et al._ (1980) R. Brandelik _et al._ , Phys. Lett. B 94, 444 (1980).
  * Althoff  _et al._ (1983) M. Althoff _et al._ , Z. Phys. C 17, 5 (1983).
  * Braunschweig  _et al._ (1989) W. Braunschweig _et al._ , Z. Phys. C 42, 189 (1989).
  * Lu (1986) X.-Q. Lu,  _Heavy quark jets from $e^{+}e^{-}$ annihilation at 29 GeV_, Ph.D. thesis, John Hopkins University (1986).
  * Aihara  _et al._ (1984) H. Aihara _et al._ , Phys. Rev. Lett. 52, 577 (1984).
  * Aihara  _et al._ (1988) H. Aihara _et al._ , Phys. Rev. Lett. 61, 1263 (1988).
  * Cowan (1988) G. D. Cowan,  _Inclusive $\pi^{\pm}$, $K^{\pm}$ and $p\bar{p}$ production in $e^{+}e^{-}$ annihilation at $\sqrt{s}=29$ GeV_, Ph.D. thesis, University of California, Berkeley (1988).
  * Derrick  _et al._ (1987) M. Derrick _et al._ , Phys. Rev. D 35, 2639 (1987).
  * Itoh  _et al._ (1995) R. Itoh _et al._ , Phys. Lett. B 345, 335 (1995), arXiv:hep-ex/9412015 .
  * Abe  _et al._ (2004) K. Abe _et al._ , Phys. Rev. D 69, 072003 (2004), arXiv:hep-ex/0310017 .
  * Buskulic  _et al._ (1995) D. Buskulic _et al._ , Z. Phys. C 66, 355 (1995).
  * Akers  _et al._ (1994) R. Akers _et al._ , Z. Phys. C 63, 181 (1994).
  * Abreu  _et al._ (1998) P. Abreu _et al._ , Eur. Phys. J. C 5, 585 (1998).
  * Albrecht  _et al._ (1989) H. Albrecht _et al._ , Z. Phys. C 44, 547 (1989).
  * Leitgab  _et al._ (2013) M. Leitgab _et al._ , Phys. Rev. Lett. 111, 062002 (2013), arXiv:1301.6183 [hep-ex] .
  * Leitgab (2013) M. Leitgab,  _Precision measurement of charged pion and kaon multiplicities in $e^{+}e^{-}$ annihilation at Q = 10.52 GeV_, Ph.D. thesis, UIUC (2013).
  * Lees  _et al._ (2013) J. Lees _et al._ , Phys. Rev. D 88, 032011 (2013), arXiv:1306.2895 [hep-ex] .
  * Abbiendi  _et al._ (2000) G. Abbiendi _et al._ , Eur. Phys. J. C 16, 407 (2000), arXiv:hep-ex/0001054 .
  * Accardi  _et al._ (2016) A. Accardi, L. T. Brady, W. Melnitchouk, J. F. Owens, and N. Sato, Phys. Rev. D 93, 114017 (2016), arXiv:1602.03154 [hep-ph] .
  * Ball  _et al._ (2017) R. D. Ball _et al._ , Eur. Phys. J. C 77, 663 (2017), arXiv:1706.00428 [hep-ph] .
  * Accardi  _et al._ (2010) A. Accardi, F. Arleo, W. Brooks, D. D’Enterria, and V. Muccifora, Riv. Nuovo Cim. 32, 439 (2010), arXiv:0907.3534 [nucl-th] .
  * Kovarik  _et al._ (2011) K. Kovarik, I. Schienbein, F. Olness, J. Yu, C. Keppel, J. Morfin, J. Owens, and T. Stavreva, Phys. Rev. Lett. 106, 122301 (2011), arXiv:1012.0286 [hep-ph] .
  * Kalantarians  _et al._ (2017) N. Kalantarians, C. E. Keppel, and M. E. Christy, Phys. Rev. C 96, 032201 (2017), arXiv:1706.02002 [hep-ph] .
  * Hirai  _et al._ (2007) M. Hirai, S. Kumano, T.-H. Nagai, and K. Sudoh, Phys. Rev. D 75, 094009 (2007), arXiv:hep-ph/0702250 .
  * (58) C. Cocuzza et al., in preparation (2021).
  * Altarelli  _et al._ (1979) G. Altarelli, R. K. Ellis, G. Martinelli, and S.-Y. Pi, Nucl. Phys. B160, 301 (1979).
  * Nason and Webber (1994) P. Nason and B. R. Webber, Nucl. Phys. B421, 473 (1994).
  * Furmanski and Petronzio (1982) W. Furmanski and R. Petronzio, Z. Phys. C 11, 293 (1982).
  * Graudenz (1994) D. Graudenz, Nucl. Phys. B432, 351 (1994), arXiv:hep-ph/9406274 .
  * de Florian  _et al._ (1998) D. de Florian, M. Stratmann, and W. Vogelsang, Phys. Rev. D 57, 5811 (1998), arXiv:hep-ph/9711387 .
  * Anderle  _et al._ (2013) D. P. Anderle, F. Ringer, and W. Vogelsang, Phys. Rev. D 87, 034014 (2013), arXiv:1212.2099 [hep-ph] .
  * Boglione  _et al._ (2019) M. Boglione, A. Dotson, L. Gamberg, S. Gordon, J. Gonzalez-Hernandez, A. Prokudin, T. Rogers, and N. Sato, JHEP 10, 122 (2019), arXiv:1904.12882 [hep-ph] .