Extraction of the Sivers function with deep neural networks
Deep Neural Networks (DNNs) are a powerful and flexible tool for information
extraction and modeling. In this study, we use DNNs to extract the Sivers
functions by globally fitting Semi- Inclusive Deep Inelastic Scattering (SIDIS)
and Drell-Yan (DY) data. To make predictions of this Transverse
Momentum-dependent Distribution (TMD), we construct a minimally biased model
using data from COMPASS and HERMES. The resulting Sivers function model,
constructed using SIDIS data, is also used to make predictions for DY
kinematics specific to the valence and sea quarks, with careful consideration
given to experimental errors, data sparsity, and complexity of phase space.
I. P. Fernando, D. Keller
# A Modern Global Extraction of the Sivers Function

###### Abstract

Deep Neural Networks (DNNs) are a powerful and flexible tool for information extraction and modeling. In this study, we use DNNs to extract the Sivers functions by globally fitting Semi-Inclusive Deep Inelastic Scattering (SIDIS) and Drell-Yan (DY) data. To make predictions of this Transverse Momentum-dependent Distribution (TMD), we construct a minimally biased model using data from COMPASS and HERMES. The resulting Sivers function model, constructed using SIDIS data, is also used to make predictions for DY kinematics specific to the valence and sea quarks, with careful consideration given to experimental errors, data sparsity, and complexity of phase space.

###### Contents

* I Introduction
* II Kinematics and Formalism
* II. SIDIS process
* II. DY process
* III Fitting \(\mathcal{N}_{q}(x)\)
* IV The Extraction Technique
	* IV.1 Data selection
	* IV.2 MINUIT fits for \(SU(3)_{\text{flavor}}\)
	* IV.3 DNN Method Testing
	* IV.4 DNN model from real data
* V Results
	* V.1 DNN fit to SIDIS data
	* V.2 Sivers in Momentum Space
	* V.3 Sivers First Transverse Moment
	* V.4 Projections
		* V.4.1 SIDIS Projections
		* V.4.2 DY Projections
		* V.4.3 The 3D Tomography of Proton
* VI Conclusion and Discussion

## I Introduction

The development of modern information extraction techniques has far outpaced experimental progress. Even with decades of experimental results, there is limited data on Transverse Momentum Dependent Parton Distribution Functions (TMDs) for global analysis that can be applied to a 3D phenomenological interpretation. Theoretical efforts are providing an ever-evolving framework and toolbox to interpret the data, leading to new areas of research that optimize information extraction unique to spin physics and the internal structure of hadrons. Artificial intelligence (AI) is accelerating data-driven research with the superior capacity of deep neural networks (DNNs) to be used as function approximators. DNNs can learn to approximate relationships contained within data, provided there are no limits to the number of neurons, layers, and computing power. With enough data, such approximations can be made with high accuracy.

The Sivers function is the most studied of the eight leading-twist TMD distributions that pertain to polarized nucleons. The Sivers distribution is naively time-reversal odd and is expected to be process-dependent, which leads to a distribution equal in magnitude but opposite in sign in Semi-Inclusive Deep Inelastic Scattering (SIDIS) compared to the Drell-Yan (DY) process [1; 2]. The Sivers distribution, as measured from the transverse single-spin asymmetries, provides information on the correlation between the nucleon's spin and the angular distribution of outgoing hadrons in SIDIS or the di-muons in DY. The quark Sivers function, \(\Delta^{N}f_{q/p^{\uparrow}}\), describes the number density in momentum space of unpolarized quarks inside the transversely polarized target nucleon with nuclear spin. A non-zero Sivers function indicates a contribution of quark orbital angular momentum to the target's spin.

TMD extraction and modeling with sensitivity to TMD evolution are critical as they provide predictive power in the collinear limit. Representing evolution also accounts for the momentum-dependent QCD interactions between partons inside the hadron, which affect their distribution and can significantly impact observables. Some theoretical approaches use the parton model approximation without TMD evolution and demonstrate good agreement with experimental results. These calculations assume that evolution effects in asymmetries are suppressed, as asymmetries are ratios of cross-sections where evolution and higher-order effects should cancel out [3; 4; 5]. This approach may not capture the full complexity of QCD evolution for TMD distributions. On the other hand, studies that incorporated TMD evolution,
as seen in [6; 7], encountered difficulties and did not achieve better agreement with the Drell-Yan data compared to earlier analyses [4; 5; 8]. This situation poses a challenge in establishing the status of the TMD factorization theorem since its principal components include scale dependence and the nonperturbative and universal Collins-Soper (CS) kernel. However, the work in [9; 10] demonstrated the conditional universality of the Sivers function using a simultaneous fit to SIDIS, Drell-Yan, and W/Z boson production data, including TMD evolution and the universal nonperturbative CS kernel extracted in SV19 [11] from unpolarized measurements. In contrast, the working principle of DNNs as a means of extraction can facilitate the necessary flexibility to capture not only the TMD and DGLAP evolution features but also the full complexity of QCD. The challenge, of course, is that once captured in the model, disentangling these feature components is nontrivial, and we expect it to be the focus of much future work.

The phenomenology used to interpret and analyze experimental data relies on TMD factorization [12; 13; 14; 15; 16; 17; 18; 19; 20; 21], proven for single-spin asymmetries (SSAs) described in terms of convolutions of TMDs. The Sivers function has been previously extracted from SIDIS data by several groups, with generally consistent results [4; 22; 23; 24; 25; 26]. However, all previous phenomenological fits of the Sivers function (and other TMDs) require a Gaussian ansatz characterizing the shape of the distribution combined with an assumed form of the Bjorken-x-dependent \(\mathcal{N}_{q}(x)\). This leads to ambiguity in determining both the quantitative results from the fit and the qualitative features of the momentum distributions and their associated dynamics. The function form of \(\mathcal{N}_{q}(x)\) is usually offered only as a placeholder and is assumed to at least contain the appropriate ingredients to facilitate the extraction. This is undoubtedly a considerable simplification but one that has permitted significant progress. In the following analysis, we perform a global fit with the goal of testing the extraction ability of DNNs to maximize information extraction and minimize both the fit error and the analytical ambiguity associated with the interpretation of \(\mathcal{N}_{q}(x)\).

The exceptional capacity of DNNs to be ideal for function approximation is rigorously provable through the Universal Approximation Theorem [27; 28]. This is the advantage of DNNs over other machine learning (ML) approaches. In this regard, even the mere existence of a function implies that DNNs can be used to represent it and work with it without actually knowing the function form. With such a high-level abstraction, one can make use of the available data and make assessments not otherwise possible, even given an arbitrary degree of complexity. The complexity can be contained in the data relationships as well as in the experimental uncertainties. In order to make optimal use of experimental uncertainties, a detailed analysis must be provided on their estimated scale and correlations.

DNNs are also Turing-complete, implying the potential to simulate the computational aspects of any real-world general-purpose computing process. The implications are that there is potential for a type of generalizable framework that can be utilized and further developed over time without the knowledge of the exact rigorous details of the underlying mechanisms. Provided appropriately detailed global models, higher-level symbolic regressions can then be performed to infer the strict mathematical form. However, such an approach requires access to a significant amount of experimental data that holds the necessary information. Without the necessary amount of quality data, no matter the number of nodes or sophistication of architecture, DNNs are limited in what useful information they can extract, along with any other technique. Even with this constraint, DNNs can make considerable advancements with the use of sparse data with large experimental uncertainties. Generally, fitting with large errors or performing computational tasks with inherent fuzzy logic are tasks that are difficult to make optimal use of modern computational resources. DNNs are uniquely suited for such challenges.

In the remainder of this article, a first-level DNN extraction of \(\mathcal{N}_{q}(x)\) is performed to deduce the Sivers function from a global analysis using HERMES and COMPASS data. This investigation is exploratory, with the intention of developing tools and techniques that minimize error and maximize utility, which we hope to expand upon in further work. In the next section, Section II, we present the formalism of the Sivers function and the kinematics for both SIDIS and DY. In Section III, a discussion of the fitting techniques of \(\mathcal{N}_{q}(x)\) is presented with a focus on the methodology of the DNN approach. Section IV explains in detail the extraction technique, starting with model testing. We perform a baseline fit using the classical Minuit \(\chi^{2}\) minimization algorithm and then perform the DNN fit, demonstrating with pseudo-data the fidelity of the procedure. We then walk through the final DNN fit to experimental data for the polarized proton and the deuteron separately. The results of the fits are presented in Section V, and finally, in Section VI, some concluding remarks are provided.

## II Kinematics and Formalism

With the spin of the proton perpendicular to the transverse plane, the Sivers function is expected to reflect an anisotropy of quark momentum distributions for the up and down quarks, indicating that their motion is in opposite directions [1; 2]. This is manifestly due to quark orbital angular momentum (OAM). The most interesting and relevant aspects of the OAM, such as magnitude and partonic distribution shape as a function of the proton's state, cannot be determined by the Sivers effect alone. However, systematic studies can be performed to investigate the full 3D momentum distribution of the quarks in a transversely polarized proton, which can be used in concert with other information to exploit multi-dimensional partonic degrees of freedom using a variety
of hard processes. Here, we focus specifically on SIDIS and DY, but it should be noted that there is significant potential for broader model development that can come from combining all available data from multiple processes with additional constraints using the simultaneous DNN fitting approach presented here.

The Sivers function describes a difference in probabilities, which implies that it may not be positive definite. Making a comparison between the Sivers function from the DY process and the SIDIS process is still the focus of much experimental and theoretical effort. Under time reversal, the future-pointing Wilson lines are replaced by past-pointing Wilson lines that are appropriate for factorization in the DY process. This implies that the Sivers function is not uniquely defined and cannot exhibit process universality, as it depends on the contour of the Wilson line. This feature of the Sivers function is directly tied to the QCD interactions between the quarks (or gluons) active in the process, resulting in a conditional universality, as shown in [29],

\[\Delta^{N}f_{q/p^{\uparrow}}\left(x,k_{\perp}\right)\bigr{|}_{\text{SIDIS}}=- \left.\Delta^{N}f_{q/p^{\uparrow}}\left(x,k_{\perp}\right)\right|_{\text{DY}}. \tag{1}\]

This fundamental prediction still needs to be tested. Direct sign tests [30; 4; 8] can be performed, but experimental proof would require an analysis over a broad phase space of both SIDIS and DY, with consideration given to flavor and kinematic sensitivity for both valence and sea quarks. Our analysis will, in part, rely on this relationship rather than making direct tests of the validity of the sign change.

### SIDIS process

The Semi Inclusive Deep Inelastic Scattering (SIDIS) process involves scattering a lepton off of a polarized nucleon and measuring the scattered lepton and a fragmented hadron. In the nucleon-photon center of mass frame, the nucleon three-momentum \(\vec{p}\) is along the \(z\)-axis and its spin-polarization \(\vec{S}_{T}\) is on the plane perpendicular (transverse) to the \(\hat{z}\)-axis. In Fig. 1 the structquark, virtual-photon (with four-momentum \(\vec{q}\)), and the lepton belong to a plane called "Lepton Plane" (represented in yellow). The fragmented-hadron with momentum \(\vec{p}_{h}\) and its projection onto the \(\hat{x}-\hat{y}\) (i.e. \(\vec{p}_{hT}\)) belong to so-called "Hadron Plane" (represented in transparent). Thus, the transverse momentum \(\vec{k}_{\perp}\) of the structquark and \(\vec{p}_{hT}\) are falling onto the transverse-plane (represented in transparent-blue) perpendicular to both lepton plane and hadron plane. The azimuthal angle \(\phi_{h}\) of the produced hadron \(h\), and is the angle between the lepton plane and the hadron plane [31]. The differential cross-section for the SIDIS process depends on both collinear parton distribution functions (PDFs) \(f_{q/p}(x;Q^{2})\) and fragmentation functions \(D_{h/q}(z;Q^{2})\), where \(q\) is the quark flavor, \(p\) represents the target proton, \(h\) is the hadron type produced by the process, and \(z\) is the momentum fraction of the final state hadron with respect to the virtual photon. A simplified version of the SIDIS differential cross-section can be written up to \(\mathcal{O}(k_{\perp}/Q)\) as [32; 25],

\[\frac{d^{5}\sigma^{lp\to lhX}}{dxdQ^{2}dzd^{2}p_{\perp}}=\sum_{q}e _{q}^{2}\int d^{2}\mathbf{k}_{\perp}\ \left(\frac{2\pi\alpha^{2}}{x^{2}s^{2}}\frac{\hat{s}^{2}+\hat{u}^{2}}{Q^{4}}\right)\] \[\times\hat{f}_{q/p^{\uparrow}}(x,k_{\perp})D_{h/q}(z,p_{\perp})+ \mathcal{O}(k_{\perp}/Q)\, \tag{2}\]

where \(\hat{s},\hat{u}\) are partonic Mandelstam invariants, and \(\hat{f}_{q/p^{\uparrow}}(x,k_{\perp})\) is the unpolarized quark distribution,

\[\hat{f}_{q/p^{\uparrow}}(x,k_{\perp}) =f_{q/p}(x,k_{\perp})+\frac{1}{2}\Delta^{N}f_{q/p^{\uparrow}}(x,k _{\perp})\vec{S}_{T}\cdot(\hat{p}\times\hat{k}_{\perp})\] \[=f_{q/p}(x,k_{\perp})-\frac{k_{\perp}}{m_{p}}f_{1T}^{\perp g}(x,k _{\perp})\vec{S}_{T}\cdot(\hat{p}\times\hat{k}_{\perp}) \tag{3}\]

with transverse momentum \(k_{\perp}\) inside a transversely polarized (with spin \(\vec{S}_{T}\)) proton with three-momentum \(\vec{p}\), where \(\Delta^{N}f_{q/p^{\uparrow}}(x,k_{\perp})\) denotes Sivers functions that carry the nucleon's spin-polarization effects on the quarks which can be considered as a modulation to the unpolarized quark PDFs [4],

\[\Delta^{N}f_{q/p^{\uparrow}}(x,k_{\perp})=2\mathcal{N}_{q}(x)h(k_{\perp})f_{q/ p}(x,k_{\perp}) \tag{4}\]

where,

\[f_{q/p}(x,k_{\perp}) =f_{q}(x)\frac{1}{\pi\langle k_{\perp}^{2}\rangle}e^{-k_{\perp}^ {2}/\langle k_{\perp}^{2}\rangle}, \tag{5}\] \[h(k_{\perp}) =\sqrt{2e}\frac{k_{\perp}}{m_{1}}e^{-k_{\perp}^{2}/m_{1}^{2}}. \tag{6}\]

Here \(\mathcal{N}_{q}(x)\) is considered as a factorized \(x\)-dependent function with a form that has yet to be formally established, and \(m_{1}\) is a parameter that allows the \(k_{\perp}\) Gaussian dependence of the Sivers function to be different

Figure 1: Kinematics of the SIDIS process in the nucleon-photon center-of-mass frame.
from that of the unpolarized TMDs [4]. \(f_{q}(x;Q^{2})\) is the co-linear PDF for flavor \(q\) that is obtained from CTEQ6l [33] grid through LHAPDF [34], whereas the fragmentation functions for \(\pi^{\pm,0}\) are from [35], and for \(K^{\pm}\) are from [36] (DSS _formalism_), from recent global analyses of fragmentation functions at next-to-leading-order (NLO) accuracy in QCD. In terms of the cross-section ratios, the Sivers asymmetry in the SIDIS process can be written as,

\[A_{UT}^{\sin(\phi_{h}-\phi_{S})}(x,y,z,p_{hT})=\frac{d\sigma^{l \uparrow p\to hlX}-d\sigma^{l\downarrow p\to lhX}}{d\sigma^{l\uparrow p \to hlX}+d\sigma^{l\downarrow p\to hlX}}, \tag{7}\]

and can be parameterized [4] and further re-arranged as,

\[A_{UT}^{\sin(\phi_{h}-\phi_{S})}(x,z,p_{hT})\] \[=\mathcal{A}_{0}(z,p_{hT},m_{1})\left(\frac{\sum_{q}\mathcal{N} _{q}(x)e_{q}^{2}f_{q}(x)D_{h/q}(z)}{\sum_{q}e_{q}^{2}f_{q}(x)D_{h/q}(z)}\right), \tag{8}\]

where,

\[\mathcal{A}_{0}(z,p_{hT},m_{1})\] \[\quad=\frac{\sqrt{2e}zp_{hT}}{m_{1}}\frac{[z^{2}\langle k_{\perp} ^{2}\rangle+\langle p_{\perp}^{2}\rangle]\langle k_{S}^{2}\rangle^{2}}{[z^{2} \langle k_{S}^{2}\rangle+\langle p_{\perp}^{2}\rangle]^{2}\langle k_{\perp}^{ 2}\rangle}\] \[\quad\times\exp\left[-\frac{p_{hT}^{2}z^{2}\left(\langle k_{S}^{ 2}\rangle-\langle k_{\perp}^{2}\rangle\right)}{(z^{2}\langle k_{S}^{2} \rangle+\langle p_{\perp}^{2}\rangle)\left(z^{2}\langle k_{\perp}^{2}\rangle+ \langle p_{\perp}^{2}\rangle\right)}\right], \tag{9}\]

\[\langle k_{S}^{2}\rangle=\frac{m_{1}\langle k_{\perp}^{2}\rangle}{m_{1}^{2}+ \langle k_{\perp}^{2}\rangle}, \tag{10}\]

and fragmentation functions \(D_{h/q}(z,p_{\perp})\) (before \(p_{\perp}\)-integration),

\[D_{h/q}(z,p_{\perp})=D_{h/q}(z)\frac{1}{\pi\langle p_{\perp}^{2}\rangle}\exp^ {-p_{\perp}^{2}/\langle p_{\perp}^{2}\rangle}, \tag{11}\]

with \(\langle k_{\perp}^{2}\rangle=0.57\pm 0.08\) GeV\({}^{2}\) and \(\langle p_{\perp}^{2}\rangle=0.12\pm 0.01\) GeV\({}^{2}\) from the fits [37; 38] to HERMES multiplicities [39]. Note that we use the shorthand notation for the PDFs, FFs as well as TMDs by omitting \(Q^{2}\) in the expressions for the sake of convenience as is done in the literature.

Through this azimuthal asymmetry, the SIDIS process provides information about the correlations between the transverse momentum of the partons leaving through the fragmented target and the spin of the target itself. In this regard, SIDIS allows one to study the structure of individual hadrons by selecting these decay fragments at the detection level. In general, SIDIS provides access to a wide range of TMDs, and allows for studying TMDs of hadrons carrying different flavors and polarizations.

For our present analysis, HERMES and COMPASS have the best-polarized proton target data for SIDIS, while COMPASS has the best-polarized neutron target data. In the COMPASS data, the neutron target is actually a polarized deuteron but the neutron carries over 90% of the deuteron polarization when polarized in solid-state form. The JLab data on polarized \({}^{3}\)He is of a different class of experiments and will not be combined with the polarized deuteron data from COMPASS. It is worth noting that the uncertainties in the experimental data can greatly differ depending on the choice of polarized target.

### DY process

Consider the Drell-Yan process \(A^{\uparrow}B\to l^{+}l^{-}X\), where \(A^{\uparrow}\) is a transversely polarized target, and \(B\) is the hadron beam. In the hadronic c.m frame, the 4-momentum \(q\) and the invariant mass squared (\(Q_{M}\)) of the final-state di-lepton pair, Feynman-\(x\) (\(x_{F}\)) and the Mandelstam variable \(s\) are related as,

\[q=(q_{0},q_{T},q_{L})\,, q^{2}=Q_{M},\qquad x_{F}=\frac{2q_{L}}{\sqrt{s}},\] \[s=(p_{A}+p_{B})^{2}\,. \tag{12}\]

In the kinematical region of,

\[q_{T}^{2}\ll Q_{M},\qquad k_{\perp}\simeq q_{T}, \tag{13}\]

at order \(\mathcal{O}(k_{\perp}/Q_{M})\), and in the hadronic c.m frame, the Sivers Single Spin Asymmetry can be given as [3; 40],

Figure 2: Kinematics of the DY process in the hadronic center-of-mass frame.
\[A_{N}^{\sin(\phi_{\gamma}-\phi_{S})} (x_{F},Q_{M},q_{T})=\frac{\int_{0}^{2\pi}d\phi_{\gamma}\left(d\sigma^ {A^{\dagger}B\to l^{+}l^{-}X}-d\sigma^{A^{\dagger}B\to l^{+}l^{-}X}\right)\sin( \phi_{\gamma}-\phi_{S})}{\frac{1}{2}\int_{0}^{2\pi}d\phi_{\gamma}\left(d\sigma ^{A^{\dagger}B\to l^{+}l^{-}X}+d\sigma^{A^{\dagger}B\to l^{+}l^{-}X}\right)}, \tag{14}\] \[=\frac{\int_{0}^{2\pi}d\phi_{\gamma}\left(\sum_{q}\int d^{2}k_{ \perp 2}d^{2}k_{\perp 1}\delta^{2}(k_{\perp 1}+k_{\perp 1}-q_{T})\Delta^{N}f_{q/A^{ \dagger}}(x_{1},k_{\perp 1})f_{\bar{q}/B}(x_{2},k_{\perp 2})\hat{\sigma}_{0}^{q \bar{q}}\right)\sin(\phi_{\gamma}-\phi_{S})}{\int_{0}^{2\pi}d\phi_{\gamma} \left(\sum_{q}\int d^{2}k_{\perp 2}d^{2}k_{\perp 1}\delta^{2}(k_{\perp 1}+k_{ \perp 1}-q_{T})f_{q/A}(x_{1},k_{\perp 1})f_{\bar{q}/B}(x_{2},k_{\perp 2})\hat{ \sigma}_{0}^{q\bar{q}}\right)} \tag{15}\]

where,

\[\hat{\sigma}_{0}^{q\bar{q}} =e_{q}^{2}\frac{4\pi\alpha^{2}}{9Q_{M}}, \tag{16}\] \[x_{1,2} =\frac{\pm x_{F}+\sqrt{x_{F}^{2}+4Q_{M}/s}}{2}. \tag{17}\]

Note that here we follow the same convention as in [41, 42, 3, 24] for the azimuthal angle in the \(A-B\) center-of-mass frame with the hadron \(A^{\dagger}\) moving along the positive \(z\)-axis and hadron \(B\) along negative \(z\)-axis. Thus the mixed product \(\vec{S}_{T}\cdot(\hat{p}\times\hat{k}_{\perp i})\) upon integration in \(k_{\perp}\), (where \(i=\{1,2\}\)) yields a \(\sin(\phi_{\gamma}-\phi_{S})=\cos\phi_{\gamma}\) (when \(\phi_{S}=\pi/2\)) dependence for the Sivers asymmetry, which implies an overall \([-\sin^{2}(\phi_{\gamma}-\phi_{S})]\) in Eq.(15). For the case in which polarized hadron \(A^{\dagger}\) moves along the \(-\hat{z}\) axis (i.e. for the processes \(BA^{\dagger}\to l^{+}l^{-}X\)), the corresponding overall factor is \([+\sin^{2}(\phi_{\gamma}-\phi_{S})]\). The analytical integration of the numerator and denominator of Eq. (15) can be written as,

\[A_{N}^{\sin(\phi_{\gamma}-\phi_{S})} (x_{F},Q_{M},q_{T})\] \[=\frac{\int d\phi_{\gamma}\:\mathcal{C}(x_{F},Q_{M},q_{T},\phi_{ \gamma})\sin(\phi_{\gamma}-\phi_{S})}{\int d\phi_{\gamma}\:\mathcal{D}(x_{F},Q _{M},q_{T})}, \tag{18}\]

where,

\[\mathcal{C} (x_{F},Q_{M},q_{T},\phi_{\gamma})\] \[\equiv\frac{d^{4}\sigma^{\uparrow}}{dx_{F}dQ_{M}d^{2}q_{T}}- \frac{d^{4}\sigma^{\downarrow}}{dx_{F}dQ_{M}d^{2}q_{T}} \tag{19}\] \[=\frac{4\pi\alpha^{2}}{9sQ_{M}}\left(\frac{q_{T}}{m_{1}}\sqrt{2e} \frac{\langle k_{S}^{2}\rangle^{2}\exp[-q_{T}^{2}/\left(\langle k_{S}^{2} \rangle+\langle k_{\perp 2}^{2}\rangle\right)]}{\pi\left(\langle k_{S}^{2}\rangle+ \langle k_{\perp 2}^{2}\rangle\right)^{2}\langle k_{\perp 2}^{2}\rangle}\right)\] \[\times\sin(\phi_{\gamma}-\phi_{S})\sum_{q}\frac{e_{q}^{2}}{x_{1}+x _{2}}\Delta^{N}f_{q/A^{\dagger}}(x_{1})f_{\bar{q}/B}(x_{2}), \tag{20}\]

and

\[\mathcal{D} (x_{F},Q_{M},q_{T})\] \[\equiv\frac{1}{2}\left[\frac{d^{4}\sigma^{\uparrow}}{dx_{F}dQ_{M }d^{2}q_{T}}+\frac{d^{4}\sigma^{\downarrow}}{dx_{F}dQ_{M}d^{2}q_{T}}\right]\] \[=\frac{4\pi\alpha^{2}}{9sQ_{M}}\left(\frac{\exp[-q_{T}^{2}/\left( \langle k_{\perp 1}^{2}\rangle+\langle k_{\perp 2}^{2}\rangle\right)]}{\pi\left( \langle k_{\perp 1}^{2}\rangle+\langle k_{\perp 2}^{2}\rangle\right)}\right)\] \[\times\sum_{q}\frac{e_{q}^{2}}{x_{1}+x_{2}}f_{q/A}(x_{1})f_{\bar{ q}/B}(x_{2}), \tag{21}\]

and it can be further simplified as,

\[A_{N}^{\sin(\phi_{\gamma}-\phi_{S})}(x_{F},M,q_{T})\] \[=\mathcal{B}_{0}(q_{T},m_{1})\frac{\sum_{q}\frac{e_{q}^{2}}{x_{1}+ x_{2}}\mathcal{N}_{q}(x_{1})f_{q/A}(x_{1})f_{\bar{q}/B}(x_{2})}{\sum_{q}\frac{e_{q}^{2}}{x_{1}+ x_{2}}f_{q/A}(x_{1})f_{\bar{q}/B}(x_{2})} \tag{22}\]

where,

\[\mathcal{B}_{0}(q_{T},m_{1})=\frac{q_{T}\sqrt{2e}}{m_{1}}\frac{Y_{1}(q_{T},k_{ S},k_{\perp 2})}{Y_{2}(q_{T},k_{\perp 1},k_{\perp 2})} \tag{23}\]

and,

\[Y_{1}(q_{T},k_{S},k_{\perp 2}) =\Biggl{(}\frac{\langle k_{S}^{2}\rangle^{2}}{\langle k_{\perp 2}^{2} \rangle\left(\langle k_{S}^{2}\rangle+\langle k_{\perp 2}^{2}\rangle\right)^{2}}\Biggr{)}\] \[\times\exp\bigg{(}\frac{-q_{T}^{2}}{\langle k_{S}^{2}\rangle+ \langle k_{\perp 2}^{2}\rangle}\Biggr{)}, \tag{24}\]

\[Y_{2}(q_{T},k_{\perp 1},k_{\perp 2}) =\Biggl{(}\frac{1}{\langle k_{\perp 1}^{2}\rangle+\langle k_{\perp 2}^{2} \rangle}\Biggr{)}\] \[\times\exp\bigg{(}\frac{-q_{T}^{2}}{\langle k_{\perp 1}^{2}\rangle+ \langle k_{\perp 2}^{2}\rangle}\Biggr{)}, \tag{25}\]

\[\frac{1}{\langle k_{S}^{2}\rangle}=\frac{1}{m_{1}^{2}}+\frac{1}{\langle k_{\perp 1 }^{2}\rangle}, \tag{26}\]

with the assumption \(\langle k_{\perp 1}^{2}\rangle=\langle k_{\perp 2}^{2}\rangle=\langle k_{\perp}^{2} \rangle=0.25\) GeV\({}^{2}\) as in [3].

Through this azimuthal asymmetry, the Drell-Yan process allows one to preferentially probe from the target and
beam hadrons to create the quark anti-quark annihilation process of interest resulting in a dimuon pair in the detector. SIDIS only permits the measurement of a convolution of the TMDs function with a fragmentation function, whereas Drell-Yan allows the direct measurement of the TMDs without the complications of fragmentation functions and final state interactions. Coupled with its innate sensitivity to sea quarks, Drell-Yan is a critical process for determining the TMDs of the sea quarks.

## III Fitting \(\mathcal{N}_{q}(x)\)

To obtain accurate three-dimensional tomographic information on quarks and gluons inside the nucleon, it is critical to extract TMDs with minimal model dependence and little to no unknown biases. Fitting and statistical analysis tools such as MINUIT rely on \(\chi^{2}\) minimization or log-likelihood functions to compute the best-fit parameter values and uncertainties, including correlations between parameters. This class of algorithms has well-established statistical methods that have been used for decades in various scientific fields, making them a reliable and trusted tool. In frequentist statistics, the reliability of a \(\chi^{2}\) minimization fitting method can be evaluated through the concept of hypothesis testing. The fitting method minimizes the difference between the observed data and the expected theoretical model, expressed through a \(\chi^{2}\) statistic. The \(\chi^{2}\) statistic follows a known distribution, and the probability of obtaining a value as extreme as the observation can be directly calculated. The reliability of the \(\chi^{2}\) minimization fitting method depends greatly on the ability to accurately estimate the theoretical uncertainties and the degree to which the model approximates the observed data. When these conditions are met, the method can provide a reliable estimate of the parameters that describe the model and its uncertainties. However, chi-square fits can be sensitive to the choice of initial parameter values and may not always converge to the correct solution. Fitting with DNNs can provide considerable advantages and does not inherently sacrifice the statistical framework provided by chi-square fits but it's worth touching on some key attributes needed in the method in order to best maintain quality statistical relevance and interpretation of resulting fits.

To preserve the statistical robustness and reliability of traditional \(\chi^{2}\) minimization fitting when using a DNN, it is important to carefully consider the data quality, model selection, validation, interpretation, and testing criteria. The quality of the data used to train the DNN should be as high as possible to ensure the DNN learns the correct relationships between inputs and outputs. This is crucial because the reliability of the DNN is only as good as the quality of the data it is trained on. Quantifying differences between the training data and the real data used in the fit can be challenging and can lead to unknown biases and systematic errors. In the method used here, Monte Carlo data, which has been tuned and matched to the experimental data, is utilized. This is done by successively extracting information from the experimental data to impose into the generated Monte Carlo data and then using the improved Monte Carlo data to further refine the extraction technique.

The choice of DNN architecture, activation functions, regularization techniques, and other hyperparameters should be carefully selected to minimize over-fitting and maximize generalization performance. Cross-validation techniques can be used to tune these hyperparameters and ensure the best possible fit to the data. The quality of the fit should be quantified with a metric that is well-defined and can be interpreted statistically. This could still be the \(\chi^{2}\) statistic but may also be a variety of possible loss functions. The trained DNN should be validated on an independent test dataset to ensure that it generalizes well to new data and that it does not over-fit the training data. This is critical because over-fitting can lead to an unreliable and unstable model.

When interpreting the results of the DNN fit, it's important to carefully examine the relationships between the inputs and outputs. To better understand how the DNN makes its predictions, techniques such as feature importance and attention mechanisms can be used. While DNNs can have a reliable statistical interpretation, it requires a more detailed analysis than traditional algorithms like MINUIT.

Directed testing of the model predictions and verification of reliability through multiple trials is crucial. Studies to test accuracy and precision are useful along with quantifying the robustness of the extraction method itself once a DNN architecture has been chosen. In this regard, it's important to prove that the method can be flexible as well as correct. It is not normally useful to have a model and method that yields highly accurate results but only for a fine slice of phase space under particular conditions.

The conventional chi-square minimization routines are limited in their flexibility and applicability to more complex problems because they assume a specific functional form of the relationship between inputs and outputs. In contrast, DNNs can learn complex and nonlinear relationships, making them suitable for tasks that require some degree of abstraction and where there is no known specific functional form of the relationship. DNNs are proven to be universal approximators and can handle large amounts of data while generalizing well to new data, which improves their accuracy and robustness. This is especially helpful in the present application, where DNNs can be used to build better models with new experimental data as it becomes available.

The ambiguity of \(\mathcal{N}_{q}(x)\) in the literature is one of the main motivations for applying the DNN technique to extract Sivers functions. Our goal in this paper is solely to lay out the procedure in order to strategically utilize the applications of DNNs. In most Sivers function extractions in the literature [24, 25, 26, 24, 25, 26, 43, 44, 45, 46, 47], the treatment of \(\mathcal{N}_{q}(x)\) differs either by its parameterization
or by the way \(q\) is treated in \(\mathcal{N}_{q}(x)\). For example, in [9; 10; 25; 26], all anti-quarks (\(\bar{u}\), \(\bar{d}\), and \(\bar{s}\)) were treated the same (or combined) and referred to as the "unbroken sea."

Our first step is a generalization of the MINUIT fit parameterization of \(\mathcal{N}_{q}(x)\) for all light quark flavors, and the section IV.2 summarizes the corresponding MINUIT fit results using \(iminuit\) (python version of MINUIT) [48]. In these fits we use the same dataset as in [4], and obtained the fit parameters for \(\mathcal{N}_{q}(x)\) defined as,

\[\mathcal{N}_{q}(x)=N_{q}x^{\alpha_{q}}(1-x)^{\beta_{q}}\frac{(\alpha_{q}+ \beta_{q})^{(\alpha_{q}+\beta_{q})}}{\alpha_{q}^{\alpha_{q}}\beta_{q}^{\beta_{q }}}. \tag{27}\]

This expression is generalized for all light quark flavors, where \(N_{q}\) is a scalar for quark flavor \(q\). After our consistency check, this parameterization is used as a pseudodata generator to train and test the DNN model. We emphasize that our original MINUIT fit parameterization is then used as a tool to demonstrate that the DNN model is capable of predicting (or confirming) the 19-parameter model used to generate pseudodata, as illustrated in section IV.3. After building confidence in the DNN model from these preliminary tests, we move towards extracting the Sivers functions from experimental data from the SIDIS process with a polarized-proton target (see section IV.4). Previous work on global fits to SIDIS data considered the data from polarized-proton targets, polarized-deuteron targets, and polarized \({}^{3}\)He gas targets as a combined dataset. In [9], isospin symmetry was assumed for \(f_{1T,u}^{\perp}\) and \(f_{1T,d}^{\perp}\) for COMPASS2009 [49] and COMPASS2017 [50] datasets. As the nuclear effects on the Sivers functions are not very well understood, separately extracting the same observable with different types of polarized targets that contain different nuclear effects could provide some insight. Significant data would be required for each target type. The DNN's unique capacity to manage abstraction allows for the absorption of additional complexity in a semi-model-independent way. To obtain the most information from the data and fit results we attempt to decompose some of this abstraction using the following conditions,

* DNN fits to SIDIS data from proton and deuteron targets are performed independently to obtain separate models.
* No kinematic cuts are applied to take full advantage of the available data and to allow the DNN to extract all possible features.
* A Sivers function for each light-quark flavor is obtained to ensure the \(SU(3)_{\text{flavor}}\) breaking effects in QCD are also contained.

The technique to achieve the extraction under the aforementioned conditions is somewhat novel, so explicit details are provided step by step for clarity in the next section.

## IV The extraction technique

\(\mathcal{N}_{q}(x)\) in Eq. (4), in the factorized form \(\mathcal{N}_{q}(x)h(k_{\perp})f_{q/p}(x,k_{\perp})\). This approach assumes the validity of Eq. (4), leading to model dependencies associated with this assertion and the Gaussian interpretation for the \(k_{\perp}\) distributions of the partons. This introduces a bias to the analysis, but it is a known bias, and something that can be revisited after the method has been thoroughly explored. As such, we do not consider this to be a fully model-independent extraction, but every attempt is made to make it minimally model-dependent. The DNN treatment of \(\mathcal{N}_{q}(x)\) enables the flexibility of handling all the light quark-flavors \(q=\{u,\bar{u},d,\bar{d},s,\bar{s}\}\) independently. The method presented is somewhat analogous to the approach of treating \(\mathcal{N}_{q}(x)\) as a function of \(x\) being parameterized with the added flexibility of the DNN. It is natural to extend this study towards a symbolic regression of \(\mathcal{N}_{q}(x)\) but that is not the scope of the present work. Additionally, we hope the DNN treatment of just the \(\mathcal{N}_{q}(x)\) term, among the other possible decompositions of the Sivers function, can shed light on how best to move forward and take full advantage of modern computational tools.

Figure 4: A generic representation of the DNN architecture for \(\mathcal{N}_{q}(x)\), where \(q=\{u,d,s,\bar{u},\bar{d},\bar{s}\}\), and \(a_{m}^{(n)}\) represent the node \(m\) in the hidden layer \(n\). The figure represents only up to \(n=3\) for demonstration purposes.
We use a relationship for \(\mathcal{N}_{q}(x)\), similar to the seminal work [4], as a tool to generate our pseudodata for testing accuracy and reproducibility only. Although this definition is not used in any aspect of the final DNN fit results, we acknowledge the bias in the approach, which we attempt to minimize and account for through studies of systematic uncertainty in the extraction methodology. The generic feedforward DNN structure for \(\mathcal{N}_{q}(x)\) that we use in this work is represented in Fig. (4). As we consider the \(SU(3)_{\text{flavor}}\) breaking in QCD, we have \(\mathcal{N}_{u}(x)\), \(\mathcal{N}_{\bar{u}}(x)\), \(\mathcal{N}_{d}(x)\), \(\mathcal{N}_{\bar{d}}(x)\), \(\mathcal{N}_{s}(x)\), and \(\mathcal{N}_{\bar{s}}(x)\) to handle the six light-quark-flavors independently. Bjorken-\(x\) is the only input as the initial layer of each \(\mathcal{N}_{q}(x)\), and the final layer is a single-node output. The \(m_{1}\) in \(\mathcal{A}_{0}(z,phT,m_{1})\) as defined in Eq. (9) is treated as a free parameter, with the initialization obtained from our first chi-square minimization fit (discussed in Section IV.2), and then allowed to vary throughout the DNN training process with SIDIS data, as shown in Fig. 3. The DNN model results are then used to infer the projections for both SIDIS kinematics and DY kinematics (see Fig. 13).

A deep feedforward architecture is used with the hidden layers expanded with multiple numbers of nodes which are initialized randomly with Gaussian sampling of weights around \(zero\) with a standard deviation of 0.1. The degree of potential Non-linearity is introduced into the network by the choice of the activation function. The selection of the activation function can have a substantial impact on DNN performance and training dynamics. We chose the \(Relu6\) activation function. This activation is a variant of the Rectified Linear Unit (\(ReLU\)) function. The \(ReLU6\) activation function has been shown to empirically perform better under low-precision conditions by encouraging the model to learn sparse features earlier, which is beneficial for learning complex patterns and relationships from the experimental data. We also use _Least Absolute Shrinkage and Selection Operator Regression_ which is a regularization technique used to prevent overfitting and improve the model's performance and generalization ability, while also encouraging sparsity and feature selection. We also use L1 regularization. L1 regularization encourages sparsity in the activation by adding a penalty term to the loss function that is proportional to the absolute value of the weights [51]. By adding this regularization term, the most important inputs are weighted the most, so that noisy or redundant information is discarded. The strength of the regularization is controlled by the magnitude of the regularization coefficient, which is set to \(10^{-12}\). Additionally, we use a dynamically decreasing learning rate. The learning rate is automatically reduced by 10% if the _training loss_ has not decreased within the last 200 epochs1 (i.e. _patience_ = 200). The optimizer used was _Adam_ while the loss function used was _Mean Squared Error_. During the hyperparameter optimization process, there were slight deviations in the number of layers, nodes per layer, initial learning rate, batch size, and the number of epochs, but the basics of the scheme just described remain consistent for all DNNs used.

Footnote 1: An epoch is a complete cycle of the passing of training data through the algorithm.

Our strategy is to first perform an exercise using only pseudodata to verify the extraction method that will ultimately be used on the real experimental data. First, we devise a _generating function_ for the SIDIS Sivers asymmetry data using a conventional \(\chi^{2}\)-minimization routine (MINUIT in this case), without following the popular assumption of "unbroken sea" [4] in order to generalize the treatment of quarks and antiquarks. We perform a series of conventional MINUIT fits step-wise to obtain the final 19 parameters for the case of broken \(SU(3)_{\text{flavor}}\) symmetry in QCD. Then, we produce pseudodata (or replicas) for the SIDIS asymmetry by sampling from the mock experimental errors using the _generating function_ with kinematics and binning in \(x\), \(z\) and \(p_{hT}\) as in the experimental data. Then a DNN model is constructed with all hyperparameters tuned in order to achieve the highest possible accuracy and precision. Here our nomenclature becomes quite specific and we refer to the resulting distribution of DNN fits as a _DNN model_. The first model obtained with the method for a particular set of data is referred to as the **First Iteration**. At this stage, we use the distribution of fit results to obtain the mean and the error band from the initial DNN model to re-parameterize the _generating function_ so that it produces more realistic pseudodata. The DNN fits are performed again improving the quality (both accuracy and precision) of the resulting fits to result in a **Second Iteration** DNN model. One can repeat the number of iterations until the resulting model is no longer improving within the experimental uncertainties. In this way, the DNN model approaches the best approximation of the Sivers functions in comparison to the _true_ values put into the _generating function_.

After confirming that the method works well, the extraction of Sivers function using the SIDIS experimental data is performed. We treat the data for a polarized proton-target and deuteron-target separately for two reasons. First, fitting these together would introduce another bias that would need to be managed directly. This is the case even assuming isospin symmetry in the \(u\) and \(d\)-quarks' Sivers functions. Second, our approach leaves open the possibility to explore the nuclear dependence of the Sivers functions. The construction of the DNN models for proton and deuteron is analogous. To perform the fit another **First Iteration**, as previously described, is performed by developing and tuning a DNN model using the data from the real experiment rather than pseudo-data from the generating function. In the proceeding iterations, we use the generating function in order to tune the hyperparameters to achieve the highest possible quality of fit in comparison with the results from the **First Iteration**. Once a tuned model is obtained, we perform an extended study for evaluating the algorithmic uncer
tainty2 as well as the systematic uncertainty of the DNN extraction method.

Footnote 2: Algorithmic uncertainty is the degree of increase in the distribution of the resulting fits that is not directly from propagated experimental error.

To elaborate on the pedagogy of this method, we organize the remainder of this section into the following subsections: (IV.1) Data selection, (IV.2) MINUIT fits for the case of \(SU(3)_{\text{flavor}}\), (IV.3) DNN model training with pseudodata, (IV.4) DNN model training with real experimental data.

### Data selection

No data points were left out of our dataset intentionally because they were suspect or classified as outliers. No kinematic cuts to exclude data points were applied. There is more world data available that could be included in our fit. Still, we chose to limit our data based on the similarity of process and experimental configuration to preserve consistency in this trial global fit with DNNs. In this regard, we focus our attention on the fixed target SIDIS and DY data. For the proton DNN fits, some data points will be left out of the training process for validation studies but they will be reincorporated after the appropriate number of epochs is determined for the optimal model performance. For the neutron DNN fits, the polarized \({}^{3}\)He data from Jefferson Lab [52] is used to test the new projections of the DNN model trained on the deuteron COMPASS data only.

Table 1 summarizes the kinematic coverage, the number of data points, and reaction types of the datasets that are considered in this work. In addition to the SIDIS datasets that are used in the fits, the polarized DY dataset from the COMPASS experiment is also listed in Table 1 as we demonstrate the predictive capability of the DNN model by comparing the projections with the real data points. The DY projections are made using the trained SIDIS DNN model assuming a sign change expected from _conditional_ universality. For the case of training the DNN model related to proton-target we use HERMES2009 [53], COMPASS2015 [54] and HERMES2020 [55] data points associated with 1D kinematic binning, leaving the HEREMES2020[55] data associated with the 3D kinematic binning to compare with the projections from the trained model. The COMPASS2009 [49] dataset with a polarized-deuteron target is used for the neutron Sivers extraction as a separate DNN model.

For the initial \(\chi^{2}\)-minimization fit with MINUIT the same datasets are used as in [4] for consistency which is HERMES2009 [53], COMPASS2009 [49], and COMPASS2015 [54]. This fit is described in the next subsection.

### MINUIT fits for \(Su(3)_{\text{flavor}}\)

The analysis begins with a \(\chi^{2}\)-minimization fit with MINUIT similar to the approach in [4] except we expand the number of parameters to treat each of the light-quark flavors separately. The results of the MINUIT fits are shown in Table 2. Fit 1 is from the original fit results from Anselmino _et al_ directly from [4]. Here the \(\mathcal{N}_{q}(x)\) for the \(u\) and \(d\) quark used is Eq. (27) but \(\mathcal{N}_{\bar{q}}(x)=N_{\bar{q}}\) for anti-quarks. In this fit, there are three parameters \(\alpha_{q}\) and \(\beta_{q}\) and \(N_{q}\) for each quark flavor, and for each anti-quark it's just \(N_{\bar{q}}\) plus \(m_{1}\). This results in a 9-parameter fit. Fit 2 is a test to reproduce the same parameterization as in Fit 1. We note that in Fit 2 non of the 9 parameters are fixed or have bounds imposed. Both of these first columns only consider \(u\) and \(d\) quarks and antiquarks. The Fit 1 parameters were used as the initial values to perform Fit 2. The difference in these two sets of fit parameters demonstrates the challenge of systematic consistency with this method though some pa

\begin{table}
\begin{tabular}{c c c c} Dataset & Kinematic coverage & Reaction points \\ \hline HERMES2009 & \(0.023<x<0.4\) & \(p^{\uparrow}+\gamma^{*}\rightarrow\pi^{+}\) & 21 \\ (SIDIS) & \(0.2<z<0.7\) & \(p^{\uparrow}+\gamma^{*}\rightarrow\pi^{-}\) & 21 \\
[53] & \(0.1<p_{nT}<0.9\) & \(p^{\uparrow}+\gamma^{*}\rightarrow\pi^{0}\) & 21 \\  & \(Q^{2}>1\) GeV\({}^{2}\) & \(p^{\uparrow}+\gamma^{*}\to K^{+}\) & 21 \\  & & \(p^{\uparrow}+\gamma^{*}\to K^{-}\) & 21 \\ \hline HERMES2020 & \(0.023<x<0.6\) & \(p^{\uparrow}+\gamma^{*}\rightarrow\pi^{+}\) & 27, **64** \\ (SIDIS) & \(0.2<z<0.7\) & \(p^{\uparrow}+\gamma^{*}\rightarrow\pi^{-}\) & 27, **64** \\
[55] & \(0.1<p_{nT}<0.9\) & \(p^{\uparrow}+\gamma^{*}\rightarrow\pi^{0}\) & 27 \\  & \(Q^{2}>1\) GeV\({}^{2}\) & \(p^{\uparrow}+\gamma^{*}\to K^{+}\) & 27, **64** \\  & \(p^{\uparrow}+\gamma^{*}\to K^{-}\) & 27, **64** \\ \hline COMPASS2015 & \(0.006<x<0.28\) & \(p^{\uparrow}+\gamma^{*}\rightarrow\pi^{+}\) & 26 \\ (SIDIS) & \(0.2<z<0.8\) & \(p^{\uparrow}+\gamma^{*}\rightarrow\pi^{-}\) & 26 \\
[54] & \(0.15<p_{nT}<1.5\) & \(p^{\uparrow}+\gamma^{*}\to K^{+}\) & 26 \\  & \(Q^{2}>1\) GeV\({}^{2}\) & \(p^{\uparrow}+\gamma^{*}\to K^{-}\) & 26 \\ \hline COMPASS2009 & \(0.006<x<0.28\) & \(d^{\uparrow}+\gamma^{*}\rightarrow\pi^{+}\) & 26 \\ (SIDIS) & \(0.2<z<0.8\) & \(d^{\uparrow}+\gamma^{*}\rightarrow\pi^{-}\) & 26 \\
[49] & \(0.15<p_{nT}<1.5\) & \(d^{\uparrow}+\gamma^{*}\to K^{+}\) & 26 \\  & \(Q^{2}>1\) GeV\({}^{2}\) & \(d^{\uparrow}+\gamma^{*}\to K^{-}\) & 26 \\ \hline \hline JLAB2011 & \(0.156<x<0.396\) & \({}^{3}He^{+}+\gamma^{*}\rightarrow\pi^{+}\) & 4 \\ (SIDIS) [52] & \(0.50<z<0.58\) & \({}^{3}He^{+}+\gamma^{*}\rightarrow\pi^{-}\) & 4 \\  & \(0.24<p_{nT}<0.43\) & & \\  & \(1.3<Q^{2}<2.7\) & & \\ COMPASS2017 & \(0.1<x_{N}<0.25\) & \(p^{\uparrow}+\pi^{-}\to l^{+}l^{-}X\) & 15 \\ (DY) [50] & \(0.3<x_{\pi}<0.7\) & & \\  & \(4.3<Q_{M}<8.5\) & & \\  & \(0.6<q_{T}<1.9\) & & \\ \hline \end{tabular}
\end{table}
Table 1: The SIDIS and DY datasets considered in the fits include the DY data, which is used to demonstrate the predictive capability of the DNN model. Specifically, we make projections using the trained SIDIS DNN model, assuming a sign change to predict the real experimental DY data points. For the HERMES2020 dataset, data is available with both 1D and 3D kinematic bins. The 3D bin numbers are indicated in bold font.
rameters match reasonably well. For fit 3 we use the same convention but add in the strange quark so there is an additional four parameters \(N_{s}\), \(\alpha_{s}\), \(\beta_{s}\), and \(N_{\bar{s}}\) which leads to a 13-parameter fit. In order to initialize the 13 parameters in Fit 3, we use the corresponding values for those parameters from Fit 2 and zeros for the rest. Fit 4 uses Eq. (27) for both quarks and antiquarks so that the treatment of all three light-quark-flavors is the same. In addition to the parameters from Fit 3, Fit 4 contains six more parameters for the antiquarks. The result of Fit 4 leads to a larger \(N_{\bar{q}}\) value to compensate for the fact that \(\alpha_{\bar{s}}\) and \(\beta_{\bar{s}}\) are now present in the fit. However, the motivation behind performing Fit 4 in this way is to generalize the \(\mathcal{N}_{q}(x)\) in a flavor-independent fashion for both quarks and anti-quarks. Fit 4 is the final fit that we will use to generate pseudodata for testing the DNN fits and for calculating the model's accuracy.

### DNN Method Testing

We develop a systematic method of constructing, optimizing, and testing the DNN fits by using pseudodata to ensure a quality extraction from the experimental data. Our approach uses a combination of Monte Carlo sampling and synthetic data generation. The pseudodata points are randomly generated by sampling within multi-Gaussian distributions centered around each experimental data point, with variance given by the experimental uncertainty. Many pseudodata DNN fits (instances) are performed together to obtain the uncertainty of the resulting DNN model (mean and distribution). The general approach is to use existing experimental data to parameterize a fit function and then use it to generate new synthetic data (replicas) with similar characteristics. The pseudodata is generated with a known Sivers function so that the extraction technique can be explicitly tested. An error bar is assigned to each new data point which is taken directly from the experimental uncertainties reported for the complete set of kinematic bins. This approach aims to produce pseudodata that simulates the experimental data as closely as possible with particular sensitivity to phase space. It does this so that the test metrics are also relevant for the real experimental extraction. To do this the pseudodata generator must be very well-tuned to the kinematic range of the experimental data. Hence, the _generating function_ contains as much feature space information as possible. It's important to emphasize here that the metrics that we use to quantify the improvement in the **Second Iteration** compared to the **First Iteration** are sensitive to phase space. The _accuracy_ (proximity of the mean of the DNN fits to the _true_ Sivers) is defined as,

\[\epsilon_{q}(x,k_{\perp})=\left(1-\frac{|\Delta^{N}f^{(\text{true})}_{q/p^{ \text{\tiny T}}}-\Delta^{N}f^{(\text{mean})}_{q/p^{\text{\tiny T}}}|}{\Delta ^{N}f^{(\text{true})}_{q/p^{\text{\tiny T}}}}\right)\times 100\%, \tag{28}\]

and _precision_ (the standard deviation of replicas), as

\[\sigma_{q}(x,k_{\perp})=\sqrt{\frac{\sum_{i}\left(\Delta^{N}f^{(i)}_{q/p^{ \text{\tiny T}}}-\Delta^{N}f^{(\text{mean})}_{q/p^{\text{\tiny T}}}\right)^{2 }}{N}}. \tag{29}\]

The _generating-function_ used to produce the _true_ value of the Sivers is improved in the process of optimizing the DNN hyperparameters. This approach improves the _generating-function_ and the DNN fit with each iteration. As a result, more realistic data can be generated in each iteration, which enables better hyperparameter optimization and testing for the experimental data in the subsequent iteration. Note that experimental data still refers to pseudodata replicas that are generated using the real experimental data rather than the _generating-function_.

In the pseudodata test, the same number of replicas are used in the **First Iteration** and in the **Second Iteration**. The number of replicas should be kept the same across consecutive iterations to control statistical error variation from the replicas. To most accurately propagate the experimental uncertainty using the replica approach requires a sufficient amount of replicas so that only negligible statistical error from the replicas is added. For the present study, our pseudodata uncertainty is simplified and is represented by a single error bar which contains the experimental statistical error and systematic error provided by the data publication. 3

\begin{table}
\begin{tabular}{c c c c c} Parameter & Fit 1 & Fit 2 & Fit 3 & Fit 4 \\ \hline \(m_{1}\) & 0.8\(\pm\)0.9 & 3.87\(\pm\)0.31 & 7.0\(\pm\)0.6 & 7.0\(\pm\)4.0 \\ \(N_{u}\) & 0.18\(\pm\)0.04 & 0.475\(\pm\)0.03 & 0.89\(\pm\)0.05 & 0.89\(\pm\)0.06 \\ \(\alpha_{u}\) & 1.0\(\pm\)0.6 & 2.41\(\pm\)0.16 & 2.78\(\pm\)0.17 & 2.75\(\pm\)0.11 \\ \(\beta_{u}\) & 6.6\(\pm\)5.2 & 15.0\(\pm\)1.4 & 19.4\(\pm\)1.6 & 20.0\(\pm\)2.0 \\ \(N_{\bar{u}}\) & -0.01\(\pm\)0.03 & -0.032\(\pm\)0.017 & -0.07\(\pm\)0.06 & -0.12\(\pm\)0.06 \\ \(\alpha_{a}\) & - & - & - & 0.4\(\pm\)0.5 \\ \(\beta_{\bar{s}}\) & - & - & - & 20.0\(\pm\)16.0 \\ \(N_{d}\) & -0.52\(\pm\)0.20 & -1.25\(\pm\)0.19 & -2.33\(\pm\)0.31 & -2.4\(\pm\)0.4 \\ \(\alpha_{d}\) & 1.9\(\pm\)1.5 & 1.5\(\pm\)0.4 & 2.5\(\pm\)0.4 & 2.7\(\pm\)0.6 \\ \(\beta_{d}\) & 10\(\pm\)11 & 7.0\(\pm\)2.6 & 15.8\(\pm\)3.2 & 17.0\(\pm\)4.0 \\ \(N_{\bar{d}}\) & -0.06\(\pm\)0.06 & -0.05\(\pm\)0.11 & -0.29\(\pm\)0.27 & -0.7\(\pm\)0.5 \\ \(\alpha_{\bar{d}}\) & - & - & - & 1.5\(\pm\)0.6 \\ \(\beta_{\bar{d}}\) & - & - & - & 20\(\pm\)17 \\ \(N_{s}\) & - & - & -14.0\(\pm\)10.0 & -20.0\(\pm\)40.0 \\ \(\alpha_{s}\) & - & - & 4.9\(\pm\)3.3 & 4.7\(\pm\)3.0 \\ \(\beta_{s}\) & - & - & 3.0\(\pm\)4.0 & 2.3\(\pm\)3.1 \\ \(N_{\bar{s}}\) & - & - & -0.1\(\pm\)0.2 & 20.0\(\pm\)5.0 \\ \(\alpha_{s}\) & - & - & - & 9.5\(\pm\)1.4 \\ \(\beta_{\bar{s}}\) & - & - & - & 20.0\(\pm\)14.0 \\ \hline \(\chi^{2}/N_{data}\) & 1.29 & 1.59 & 1.69 & 1.66 \\ \end{tabular}
\end{table}
Table 2: Collection of MINUIT fit results. Fit 1 is from Anselmino et al [4], Fit 2: Re-fit as similar to [4], Fit 3: fit results including strange-quarks, Fit 4: fit results with the same treatment for all three light-quark-flavors.
Footnote 3: The _Success_ is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as the _Success_, which is the same as _Success_, which is the same as the _Success
and \(\mathcal{C}_{d}^{f}\) for results from SIDIS data from experiments associated with polarized-deuterium target. Here \(i\) and \(f\) indicate the **First Iteration** and **Second Iteration** respectively. The listed learning rate (multiplied by \(10^{-4}\)) is the initial learning rate as a dynamically decreasing learning rate is used. (This is explained in Sec. IV). The accuracy \(\varepsilon_{q}(x,k_{\perp})\) is defined in Eq. (28), and the results in this table correspond to the maximum deviation of the mean of the replicas from the true values \(\varepsilon_{q}^{max}\); whereas the precision \(\sigma_{q}(x,k_{\perp})\) is defined in Eq. (29) and the results are the maximum standard deviations of the replicas \(\sigma_{q}^{max}\), and are in the units of \(\times 10^{-3}\).

It is worth noting that the improvement in both _accuracy_\(\varepsilon_{q}^{max}\) and _precision_\(\sigma_{q}^{max}\) can be observed from the closeness of the solid line (mean of the 1000 DNN replicas) to the dashed line (_generating function_) for quarks (upper plots) and antiquarks (lower plots). The improvement of the DNN model in each case is significant to the point where it is difficult to distinguish the solid line and the corresponding dashed line, indicating a high degree of accuracy and precision. Also, we observed that the training loss in the **Second Iteration** is about one order of magnitude less than the one from the **First Iteration**.

### DNN model from real data

In contrast to the testing of the DNN model with pseudodata from the _generating function_, we now describe the steps to apply the DNN fit method to real experimental data. The pseudodata test from Sec. IV.3 is using the combined proton and deuteron data as in all previous work on global fits of the Sivers function. In the following extraction with real experimental data, the proton and deuteron data are fitted separately. To take full advantage of the information provided by the model testing in the previous section, the steps from Sec. IV.3 are performed again separately for proton and deuteron data. The starting hyperparameters for the first DNN fit in this section including the architecture, initial learning rate, batch size, as well as optimal number of epochs, are all determined based on the best accuracy and precision using the well-tuned pseudodata from the generating function first in each case. The provide more information on the initial hyperparameter so that even the **First Iteration** is a quality fit. In the following steps, two distinct DNN models are developed, one for the _proton_ quarks Sivers asymmetry and one for _neutron_ quark Sivers asymmetry. The following procedure is common to both.

1. **First Iteration**: construct a DNN fit and tune its hyperparameters by training with the experimental data. Use 10% of the data for the validation in each epoch.5 Footnote 5: For this section we list this as **First Iteration** though the real experimental data has already been fit to improve the _generating function_ for the pseudodata test. With every subsequent iteration, the pseudodata test and the DNN fit to the experimental data improve.
2. Identify the optimum number of epochs when the validation loss exceeds the training loss (see Fig. (7) as an example).
3. Perform a DNN fit, with the optimized hyperparameters from Step 1 and the number of epochs determined from Step 2, using all the data without leaving any for validation.
4. Improve the _generating function_: 1. Use the tuned DNN model in step 3 to infer the asymmetry over the 3D kinematics \((x,z,p_{hT})\) in fine bins. 2. Perform a MINUIT fit to obtain the new _generating function_. 3. Produce pseudodata for the SIDIS asymmetry using the _generating function_ in the previous step.

\begin{table}
\begin{tabular}{c c c c c c c} Hyperparameter & \(\mathcal{C}_{0}^{i}\) & \(\mathcal{C}_{0}^{f}\) & \(\mathcal{C}_{p}^{i}\) & \(\mathcal{C}_{p}^{f}\) & \(\mathcal{C}_{d}^{i}\) & \(\mathcal{C}_{d}^{f}\) \\ \hline Hidden Layers & 5 & 7 & 5 & 7 & 5 & 8 \\ Nodes/Layer & 256 & 256 & 550 & 550 & 256 & 256 \\ Learning Rate & 1 & 0.125 & 5 & 1 & 10 & 1 \\ Batch Size & 200 & 256 & 300 & 300 & 100 & 100 \\ Number of epochs & 1000 & 1000 & 300 & 300 & 200 & 200 \\ Training Loss & 0.6 & 0.05 & 1.5 & 1 & 2 & 1 \\ \hline \(\varepsilon_{max}^{max}\) & 95.67 & 99.27 & 55.21 & 94.04 & 56.80 & 93.02 \\ \(\varepsilon_{a}^{max}\) & 42.62 & 98.09 & 52.57 & 96.70 & 34.83 & 91.40 \\ \(\varepsilon_{d}^{max}\) & 80.46 & 98.89 & 55.69 & 93.13 & 52.44 & 89.27 \\ \(\varepsilon_{d}^{max}\) & 74.59 & 97.08 & 55.37 & 95.44 & 46.60 & 92.58 \\ \(\varepsilon_{a}^{max}\) & 45.53 & 79.27 & 49.54 & 90.64 & 36.34 & 93.41 \\ \(\varepsilon_{s}^{max}\) & 59.27 & 91.13 & 33.89 & 82.51 & 65.57 & 91.45 \\ \hline \(\sigma_{uu}^{max}\) & 3 & 0.1 & 5 & 2 & 2 & 0.4 \\ \(\sigma_{uu}^{max}\) & 2 & 0.2 & 6 & 2 & 8 & 2 \\ \(\sigma_{d}^{max}\) & 10 & 1 & 20 & 6 & 2 & 1 \\ \(\sigma_{d}^{max}\) & 7 & 4 & 20 & 8 & 7 & 1 \\ \(\sigma_{u}^{max}\) & 2 & 0.2 & 4 & 1 & 6 & 2 \\ \(\sigma_{s}^{max}\) & 1 & 0.1 & 4 & 2 & 6 & 3 \\ \hline \end{tabular}
\end{table}
Table 3: The summary of the optimized sets of hyperparameters: The indications in the table are \(\mathcal{C}_{0}^{i}\) and \(\mathcal{C}_{0}^{f}\) for results from the pseudodata from the generating function, \(\mathcal{C}_{p}^{i}\), and \(\mathcal{C}_{p}^{f}\) for results from SIDIS data from experiments associated with polarized-proton target, \(\mathcal{C}_{d}^{i}\) and \(\mathcal{C}_{d}^{f}\) for results from SIDIS data from experiments associated with polarized-deuterium target, where \(i\) and \(f\) indicates the **First Iteration** and **Second Iteration** respectively. The initial learning rate is also listed (\(\times 10^{-4}\)) as is the final training loss (\(\times 10^{-3}\)). The accuracy and precision in each case are the maxima over the phase space.
5. Perform Step 1 and Step 2 with the pseudodata from Step 4 (from the improved _generating function_).
6. **Second Iteration**: perform a DNN fit with the optimized hyperparameters from Step 5 using experimental data without leaving any data for validation.
7. Perform a comparison of the Sivers functions extracted at the **First Iteration** vs **Second Iteration** in terms of the _accuracy_, _precision_ and the magnitude of the _loss-function_ at the final epoch.

Although the architectural specifics such as the number of hidden layers, nodes per layer, and learning rate may be modified in the hyperparameter optimization step, the number of epochs and the number of replicas remain the same. The overall feedforward architecture structure remains consistent as well for simplicity. Table 3 shows the optimized set of hyperparameter configurations are represented with \(\mathcal{C}_{\{0,p,d\}}^{\{i,f\}}\) where \(i,f\) represent the **First Iteration** and the **Second Iteration**, and \(\{0,p,d\}\) represent with _pseudodata_, with _proton_ data and with _deuteron_ data respectively.

Clearly, the pseudodata from the _generating function_ still plays a critical role in the tuning and testing of the fit of the real experimental data. Accuracy is necessarily determined using pseudodata so that the mean of the final DNN model can be compared directly to the _true_ Sivers asymmetry. This procedure for determining accuracy assumes that the _true_ Sivers asymmetry from the experiment can be approximated by the well-tuned _generating function_ after the final iteration. There is a systematic error associated with analyzing accuracy this way but this type of error can be estimated.

After completing the extraction method described above, we perform a systematic uncertainty assessment on the overall method of extraction. To test the reliability of the extraction, we adjust the parameters of the _generating function_ and repeat the extraction process by again following the full set of steps for several variations of pseudodata. By using the Sivers asymmetry generated from the optimized _generating function_, the absolute differences between the mean of the DNN model and the true value over \(k_{\perp}\) was used to estimate the systematic uncertainties of the final DNN model from this extraction technique.

## V Results

In this section, we present the results from two separate DNN models: the _proton_-DNN and the _deuteron_-DNN, along with their optimized hyperparameters. Only SIDIS data was used to train the DNN models in this exploratory AI-based extraction technique. The optimized hyperparameter configurations are provided for both models in columns \(\mathcal{C}_{p}^{i(f)}\) and \(\mathcal{C}_{d}^{i(f)}\), with the subscripts \(p\) and \(d\), respectively, in Table 3. To quantitatively represent the improvement made by performing the steps mentioned in the previous section, we present our accuracy and precision results in the lower part of Table

Figure 6: The qualitative improvement of the extracted Sivers functions for \(u\) (blue), \(d\) (red), and \(s\) (green) quarks at \(x=0.1\) and \(Q^{2}\)=2.4 GeV\({}^{2}\) using the optimized _proton_-DNN model at the **Second Iteration** (solid-lines with dark-colored error bands with 68% CL), compared to the **First Iteration** (dashed-lines with light-colored error bands with 68% CL).

Figure 7: An example for the determination of the optimum number of epochs (about 300 in this case: marked by the vertical red dashed-line) at the **First Iteration** based on cross-over coordinates between the train (blue curve) and the validation (orange curve) losses. Then, the train-loss behaviors until the optimum number of epochs with all data points at the **First Iteration** and at the **Second Iteration** are represented by the green-curve and the red-curve respectively.
### DNN fit to SIDIS data

We now explore the results and compare some of our final fits and projections with those of other global fits. To do this we calculate Pearson's reduced \(\chi^{2}\) statistic based on the results from the DNN model and experimental data point so that the \(\chi^{2}\) value provided from other work can be used to make a quantitative comparison. Note that the \(\chi^{2}\) values indicated in our plots are calculated after the analysis is complete, rather than as a part of the minimization process.

The plots of SIDIS Sivers asymmetry data and our resulting DNN models (for _proton_ and _deuteron_) are shown in Fig. 8. Each plot includes the partial \(\chi^{2}\) values for the particular \(x,z\), and \(p_{hT}\) bins for each hadron type. HERMES2009 [53] (_top-left_), HERMES2020 [55] (_top-right_), and COMPASS 2015 [54] (_bottom-left_) are described with the _proton_-DNN model with a total \(\chi^{2}/N_{pt}\) = 1.04, whereas the COMPASS2009 [49] (_bottom-right_) dataset is described with the _deuteron_-DNN with a total \(\chi^{2}/N_{pt}\) = 0.81. In comparison with [4], there are some improvements in describing the proton SIDIS data on \(\pi^{\pm}\) and \(K^{+}\) in HERMES 2009, which can be noticed quantitatively based on the partial \(\chi^{2}\) values from the DNN model. This indicates that the possible effects attributed to the TMD evolution [6; 44; 56] and were assumed to be the cause of the larger \(\chi^{2}\) values in [4] for proton data on \(\pi^{+}\) may have been somewhat integrated into the DNN model. Although HERMES 2020 [55] reported SIDIS data in 1D kinematic bins as well as with 3D kinematic bins, in our fits we use the data in the form of 1D kinematic bins to be consistent with the rest of the datasets in our fits.

The _deuteron_-DNN model's description of COMPASS2009 [49] data is shown in the bottom-left sub-figure of Fig. 8. Without applying any cuts on the data, the DNN model yields a total \(\chi^{2}/N_{pt}\) of 0.80, covering the full range in \(x\), \(z\), and \(p_{hT}\) kinematic projections from the COMPASS2009 dataset. This is in contrast to the limited kinematic coverage considered in [7; 5; 9; 47]; notably, the data points at \(p_{hT}>~{}1\) GeV are described somewhat better by the _deuteron_-DNN model compare to the fits in [4; 25]. This suggests that performing dedicated fits to data specific to polarized nucleon targets enables better information extraction, which is true for both DNN and other fitting approaches. The advantage of DNNs in this case is to perform well even with limited data.

We did not include JLab [52] data in our _deuteron_-DNN model fits to use it as a projection test for the neutron Sivers asymmetry. Our projection indicates good agreement with the \({}^{3}\)He data, but both the data and the projection are largely consistent with zero. It is also important to note that in this work we are _not_ imposing any isospin symmetry condition (\(f_{1T}^{\perp u}=f_{1T}^{\perp d}\) and/or \(f_{1T}^{\perp\bar{u}}=f_{1T}^{\perp\bar{d}}\)) for the SIDIS data with the deuteron target as was done in [9]. The successful construction of the two different proton and neutron Sivers functions may indicate that our DNN approach can be particularly useful for analyzing data from polarized nucleons in different nuclei, potentially opening up a new way of exploring the nuclear effects associated with TMDs.

### Sivers in Momentum Space

The extracted Sivers functions, including the systematic uncertainties from the DNN models at \(x=0.1\) and \(Q^{2}=2.4\) GeV, are shown in Fig. 9 represented by the _mean_ with 68% Confidence Level (CL) _error-bands_. The corresponding optimized hyperparameter configurations for the _proton_-DNN model and _deuteron_-DNN model are \(\mathcal{C}_{p}^{f}\) and \(\mathcal{C}_{d}^{f}\), respectively, as given in Table 3. The Sivers functions extracted using the _deuteron_-DNN model show consistency with zero, considering the accompanying systematic uncertainties. However, this is still a significant result, given the limitation in statistics from the SIDIS data with a deuteron target. The Sivers functions extracted using the _proton_-DNN model have small systematic uncertainties. Note that we use \(\Delta^{N}f_{q/p^{\uparrow}}(x,k_{\perp})\) notation, as in [4], to represent the Sivers functions in our plots, and one can use Eq. (3) to convert to \(f_{1T}^{\perp q}(x,k_{\perp})\) notation.

Comparing the extracted Sivers functions in \(k_{\perp}\)-space to other extractions in the literature can also be useful, although we have not included those curves in our plots. In summary, the _proton_-DNN model extractions are relatively precise, with narrower error bands compared to those in [3; 4; 6; 7; 9; 10; 25].

### Sivers First Transverse Moment

The first transverse moment of the Sivers functions can be obtained through \(d^{2}k_{\perp}\)-integration of the Sivers functions [4],

\[\Delta^{N}f_{q/p^{\uparrow}}^{(1)}(x)=\int d^{2}k_{\perp}\frac{k _{\perp}}{4m_{p}}\Delta^{N}f_{q/p^{\uparrow}}(x,k_{\perp})\] \[=-f_{1T}^{\perp(1)q}(x)=\frac{\sqrt{\frac{e}{2}}(k_{\perp}^{2})m_ {1}^{3}}{m_{p}((k_{\perp}^{2})+m_{1}^{2})^{2}}\mathcal{N}_{q}(x)f_{q}(x;Q^{2}) \tag{30}\]

The extracted first transverse moments of the Sivers functions including the systematic uncertainties from the DNN models are given in Fig. 11 with 68% CL _error-bands_ using the optimized hyperparameter configurations \(\mathcal{C}2\) and \(\mathcal{C}3\) in Table 3 respectively for _proton_-DNN model and _deuteron_-DNN model. The calculated moments using the _deuteron_-DNN model are consistent with zero, based on the systematic uncertainties.

Comparing the results in Fig. 1 of [47], we see that the \(xf_{1T}^{\perp(1)u}\) from the DNN model is more consistent with [5; 6] in the vicinity of \(x=0.1\), although it is consistent
Figure 8: The DNN fit results of the SIDIS Sivers asymmetries (red) accompanied by 68% CL error-bands in comparison with the actual data (blue). The _proton_-DNN model is trained with HERMES2009, HERMES2020, and COMPASS2015, whereas the _deuteron_-DNN model is trained with the COMPASS2009 data. The calculated partial \(\chi^{2}\) values are provided as quantitative assessments for all kinematic bins.
with [26] at \(x=0.01\). The \(xf_{1T}^{\perp(1)d}\), in general, is consistent with the extractions from [4; 5; 6; 25; 47; 57]. Additionally, the extracted behavior of \(xf_{1T}^{\perp(1)u}\) and \(xf_{1T}^{\perp(1)d}\) is consistent with the qualitative observation in [25],

\[\Delta^{N}f_{u/p^{\uparrow}}^{(1)}(x) = -\Delta^{N}f_{d/p^{\uparrow}}^{(1)}(x)\] \[\text{or}\ \ f_{1T}^{\perp(1)u}(x) = -f_{1T}^{\perp(1)d}(x) \tag{31}\]

which was originally a prediction from the large-\(N_{c}\) limit of QCD [58]. Most importantly, the DNN model is able to capture the feature of the \(u\) and \(d\) quarks orbiting in opposite directions without imposing this constraint directly as done in [42]. In terms of the quantitative assessment, Eq. (31) could be accurate at the large-\(N_{c}\) limit, if the isospin breaking effects are also included at the next to leading order in \(\mathcal{O}(1/N_{c})\).

In regards to the light sea-quarks, the _proton_-DNN model extracts the features such as \(\Delta^{N}f_{\bar{u}/p^{\uparrow}}^{(1)}(x)>0\) and \(\Delta^{N}f_{d/p^{\uparrow}}^{(1)}(x)<0\), even considering the scale of the uncertainties. Additionally, the _proton_-DNN model is consistent with

\[\Delta^{N}f_{\bar{u}/p^{\uparrow}}^{(1)}(x)=-\Delta^{N}f_{\bar{d}/p^{\uparrow} }^{(1)}(x) \tag{32}\]

which was also a similar observation from a theoretical calculation based on \(SU(2)\) chiral Lagrangian [59] and the predictions at large-\(N_{c}\) limit of QCD [58]. The central values extracted in [26] are qualitatively similar to the features seen in Fig. 11 which are small but non-zero within the uncertainties. Additionally, the corresponding central values extracted in [4] are both negative but consistent with zero.

The first transverse moments \(xf_{1T}^{\perp(1)q}(x)\), in the case of \(SU(3)_{\text{flavor}}\), from our DNN result, are more precise (narrow error bands) than those in [4; 22; 25]. However, the error bands are slightly larger than those in JAM20 [5], which includes more data from SIDIS, DY and SIA, \(pp\)-collisions, and parameterizations for Sivers, Collins, and Transversity TMDs together.

### Projections

#### iv.4.1 SIDIS Projections

In Fig. 10, we compare the SIDIS Sivers asymmetries (in red) projected onto the HERMES2020 3D kinematic bins with the experiment measurements (in blue). These results are obtained using our _proton_-DNN model and are accompanied by 68% CL error bands. We also provide calculated partial \(\chi^{2}\) values for each kinematic bin as a quantitative assessment. Unlike in [9], we have made projections for all the data points since we have not applied any data cuts. There are relatively larger partial \(\chi^{2}/N_{pt}\) (as was also observed in [9]), but only in a couple of \(K^{+}\) and \(K^{-}\) bins. However, the \(\chi^{2}/N_{pt}\) values for \(\pi^{+},\pi^{-},K^{+},K^{-}\) are, respectively, 1.14, 1.04, 1.03 and 0.88, which leads to the total \(\chi^{2}/N_{pt}\) of 1.02 for all the data points from HERMES2020 in 3D kinematic bins. In Fig. (12), we present the projected SIDIS Sivers asymmetries for the JLab kinematics [52], obtained using our _deuteron_-DNN model. The figure includes 68% CL error bands and a comparison with the JLab neutron Sivers asymmetry data. These results are consistent with those reported in [4; 6; 7; 57].

#### iv.4.2 DY Projections

The resulting DNN model based on the SIDIS Sivers asymmetries is capable of projecting the Sivers asymmetries in DY experiments which could be sensitive to either valence quarks or sea quarks depending on the relevant kinematic coverage. For example, the COMPASS2017-polarized DY Sivers asymmetry measurements [50] are dominated by the valence quarks, and the upcoming SpinQuest (E1039) experiment's polarized DY asymmetry measurements [60] will be dominated by the sea quarks. For these DY projections, we follow the block diagram represented in Fig. 13, which includes the assumption of the _sign-change_ of the Sivers function in DY relative to the SIDIS process mentioned in Eq. (1). There

Figure 9: The extracted Sivers functions from the _proton_-DNN model (upper) and _deuteron_-DNN model (lower) at \(x=0.1\) and \(Q^{2}=2.4\) GeV\({}^{2}\) with 1-\(\sigma\) (68%) CL _error-bands_, including systematic uncertainties.
fore, using the trained _proton_-DNN model, we make projections for the DY Sivers asymmetries for the COMPASS2017 experiment [50] with a proton target and a pion beam using CTEQ6l [33] and JAM21PionPDFnlo [61] for proton PDFs and pion PDFs respectively. Meanwhile, using both the _proton_-DNN model and _deuteron_-DNN model, we make predictions for the SpinQuest experiment6. The kinematic-inputs are \(x_{1}\) (beam), \(x_{2}\) (target), \(x_{F}(=x_{1}-x_{2})\), \(q_{T}\) (transverse component of the virtual photon), and \(Q_{M}\) (di-lepton invariant mass).

Footnote 6: We are using the kinematic bins from the SeaQuest experiment [62].

The projected DY Sivers asymmetries for the COMPASS2017 experiment [63] are shown in Fig. 10. The projections of the _proton_-DNN model including 68% CL error bands (in red) in comparison with the actual data points (in blue).

Figure 10: Projections of the of HERMES 2020 data for 3D kinematic bins, using the _proton_-DNN model including 68% CL error bands (in red) in comparison with the actual data points (in blue).
PASS2017 kinematics using the trained _proton_-DNN model in comparison with the data [50] points are represented in the Fig. 14. Although the projections are based on the assumption of conditional universality, it is worth noting that without this assumption, negative asymmetry projections were observed. However, for clarity, the projections without assuming conditional universality are not shown in Fig. 14. When comparing the projections of the _proton_-DNN model with the predictions from [4; 6; 45], it is evident that the mean of the _proton_-DNN model projection, in terms of \(x_{F}\), is more consistent with the measured mean Sivers asymmetry in the experiment (refer to Figure 6 in [50]). At the same time, the _proton_-DNN model projection has relatively smaller uncertainty. It is important to note that the predictions mentioned in the cited works were based on different \(Q^{2}-\) evolution schemes, while the _proton_-DNN model incorporates DGLAP evolution through LHAPDF [34]. In [9], only two data points from the COMPASS2017 data were included in their fits, resulting in larger uncertainties in the projected asymmetry values for the remaining data points when compared to the projections generated by the _proton_-DNN model. The increasing trend of the projected DY Sivers asymmetries with respect to the \(q_{T}\) kinematic variable in [9] is consistent with the projections presented in the middle-right plot of Fig. 14 generated by the _proton_-DNN model, while in [7], the corresponding trend exhibits a very small negative slope in relation to \(q_{T}\).

A non-zero sea-quark Sivers asymmetry is inferring that the sea quarks have non-zero orbital angular momentum. The _proton_-DNN model predictions exhibit consistency with the non-zero Sivers asymmetry from the _sea_-quarks, with higher precision compared to existing predictions [4; 6; 45], for the SpinQuest kinematics. Additionally, in this work, we report our projections for the polarized Drell-Yan Sivers asymmetries for a deuteron target at the SpinQuest experiment, as shown in Figure 15 by the orange-colored bands. The central lines are negative in all kinematic projections of \(x_{1},x_{2},x_{F}\), and \(q_{T}\), yet consistent with zero, which correlates to the extracted Sivers functions shown in Fig. 9 and Fig. 11. The _proton_-DNN model predicts a positive slope with respect to \(q_{T}\) for a _proton_ target and a relatively small negative slope for a _deuteron_ target, as shown in the lower-right plot of Fig. 15. To date, with the exception of this work, no predictions have been made for the polarized DY Sivers asymmetry using a _deuteron_ target, which will be measured during the SpinQuest experiment. A noteworthy aspect of the forthcoming SpinQuest experiment

Figure 11: The extracted first transverse moments of Sivers functions from the _proton_-DNN model (upper) and _deuteron_-DNN model (lower) at \(x=0.1\) and \(Q^{2}=2.4\) GeV\({}^{2}\) with 68% CL _error-bands_, including systematic uncertainties.

Figure 12: The _deuteron_-DNN model projections (red) with 68% CL error-bands, for JLab kinematics [52] in comparison with measured data (blue) without the systematic uncertainty.
is that, in addition to measuring the Sivers asymmetry from _proton_ and _deuteron_ targets, it will also ascertain the transversity distributions of both quarks and gluons, utilizing a tensor-polarized deuteron Spin 1 target, as proposed in [63].

### The 3D Tomography of Proton

The TMD density of unpolarized quarks inside a proton polarized in the \(\hat{y}\)-direction can be graphically represented using the relation [9; 47],

\[\rho^{a}_{p\uparrow}(x,k_{x},k_{y};Q^{2})=f^{a}_{1}(x,k_{\perp}^{2};Q^{2})- \frac{k_{x}}{m_{p}}f^{\perp a}_{1T}(x,k_{\perp}^{2};Q^{2}), \tag{33}\]

where \(k_{\perp}\) is a two-dimensional vector \((k_{x},k_{y})\), and the unpolarized TMD and the Sivers function for quark-flavor \(a\) are respectively represented as \(f^{a}_{1}(x,k_{\perp}^{2};Q^{2})\), and \(f^{\perp a}_{1T}(x,k_{\perp}^{2};Q^{2})\). The corresponding quark density distributions from our _proton_-DNN model for all light quark flavors in \(SU(3)_{\rm flavor}\) at \(x=0.1\) and \(Q^{2}=2.4\) GeV\({}^{2}\) are shown in Fig. 16. The observed shifts in each quark flavor are linked to the correlation between the OAM of

Figure 16: Quark density distributions \(\rho^{a}_{p\uparrow}\) from the _proton_-DNN model (average of 1000 replicas ) for the light quark flavor \(a=\{u,\bar{u},d,\bar{d},s,\bar{s}\}\) inside a proton polarized along the \(+y\) direction and moving towards the reader, as a function of \((k_{x},k_{y})\) at \(x=0.1\) and \(Q^{2}=2.4\) GeV\({}^{2}\).

Figure 14: The _proton_-DNN model’s predictions (red) including 68% CL error-bands, for Sivers asymmetries in \(x_{1},x_{2},x_{F},q_{T}\), and \(Q_{M}\) kinematic projections for COMPASS DY kinematics [50] in contrast with the measured data (blue).

Figure 15: The _proton_-DNN model (red) and the _deuteron_-DNN model (orange) predictions including 68% CL error-bands, for Sivers asymmetries in \(x_{1},x_{2},x_{F}\), and \(q_{T}\) kinematic projections for the SpinQuest DY kinematics [60; 62].
quarks and the spin of the proton. The results shown in Fig. 16 provide evidence of non-zero OAM in the wave function of the proton's valence and sea quarks. The _proton_-DNN model calculations for the \(u\) and \(d\) quarks are similar to those reported in [9; 47], where the distortion has a positive shift for the \(u\)-quark and a negative shift for the \(d\)-quark with respect to the \(+x\) direction. From the results in Fig. 16, the _proton_-DNN model demonstrates that a virtual photon traveling towards a polarized proton "sees" an enhancement of the quark distribution, in particular more \(u,\bar{u}\)-quarks to its right-hand side and more \(d,\bar{d}\)-quarks to its left-hand side in the momentum space. Moreover, the resultant shifts for \(\bar{u},s\) quarks from the _proton_-DNN model are also in agreement with [9]. In the low-x region, the momentum space quark density becomes almost symmetric [47], and it indicates that the Sivers effect becomes smaller and the corresponding experimentally observed asymmetry is small.

The forthcoming data from Jefferson Lab at 12 GeV, Fermilab SpinQuest experiment, and the anticipated future data from the Electron-Ion Collider [64; 65; 66], along with their extensive kinematic coverage, are expected to provide invaluable insights into the 3D structure of the nucleon. Obtaining a model-independent estimate of quark angular momentum requires parton distributions that simultaneously depend on both momentum and position [67; 68; 69; 70]. In addition to experimental observations, lattice QCD (LQCD) computations provide a valuable tool for QCD phenomenology from first principles. For instance, LQCD has been utilized to investigate the Sivers effect and other TMD observables at different pion masses [71] as well as the generalized parton distribution at the physical pion mass [72]. Additionally, LQCD results on the Collins-Soper kernel over a range of \(b_{T}\) (the Fourier transform of the transverse momentum) are useful for global fits of TMD observables from different processes [73]. In this way LQCD could complement the experimental data and open up an avenue to enhance the DNN method to explore the 3D structure of nucleons more directly.

## VI Conclusion and Discussion

In this paper, we propose a new method for performing global fits to extract the fundamental Sivers distributions of unpolarized quarks in both polarized protons and neutrons. Our approach integrates techniques from artificial intelligence (AI) and utilizes a generating function to ensure the quality of the extraction while maintaining a high level of accuracy and precision. The primary objective of our method is to ensure the quality of the extraction, while simultaneously maintaining an optimal level of accuracy and precision. By leveraging the power of AI, we are able to implement an algorithm that can effectively handle complex and sparse data sets, and generate outputs with higher accuracy and precision. Furthermore, the use of a generating function provides an additional layer of quality control by providing a direct check on the extraction value over the phase space of interest using pseudodata. The generating function can then be improved to better reflect the experimental information and applied again to iteratively improve the extraction.

We have analyzed the global data on the transverse single spin asymmetry measurements provided by multiple experimental collaborations, particularly for the SIDIS process. We have used this data to construct a DNN model that can effectively interpolate and extrapolate to some degree, enabling us to predict Sivers asymmetries in both SIDIS and DY processes. We have successfully tested the method of extraction and demonstrated through both pseudodata and real experimental data that the technique can generally outperform classical fitting procedures with accuracy and precision that can be directly quantified. The technique demonstrates good reproducibility as detailed by our systematic studies. Further avenues for advancement in this type of modeling are evident.

As the first attempt to extract the Sivers function using AI techniques, we chose the \(\mathcal{N}_{q}(x)\) parameterization of the deep neural net to incorporate all x-dependent features. Our results show promise, as the DNN method is capable of fitting the SIDIS Sivers asymmetry data and consistently provides extraction of the Sivers function not only for valence quarks but also for all light quark flavors in a flavor-independent fashion. The trained model is able to make projections on untrained data for both SIDIS and DY kinematics, and these predictions are consistent with experimental results demonstrating the predictive capability of the method with relatively higher precision.

Our method and techniques are intended to be simple and clear enough to be easily reproduced and expanded upon. This analysis effectively demonstrates the potential that AI holds for particular types of analysis and information extraction. The work also emphasizes the need for increased cooperation between experimentalists, theorists, and the growing computational efforts that will undoubtedly accelerate advancements in the years to come. A global effort to standardize data collection, organization, and storage in an unbinned format with detailed covariance information must be developed to take full advantage of the AI tools now available. AI and its emerging technologies will continue to accelerate the progress of data-driven physics, but the speed of that progress largely depends on the level of foresight and cooperation within the community.

###### Acknowledgements.

 This work was supported by the U. S. Department of Energy (DOE) contract DE-FG02-96ER40950. The authors would like to thank Micheal Wagman, Andrea Signori, Alexei Prokudin, and Filippo Delcarro for their
insightful discussions. Additionally, the authors would like to acknowledge Bakur Parsamyan for providing the COMPASS data, Luciano Pappalardo for providing the HERMES data, and Sassot Rodolfo for providing the Fragmentation Functions.
