Studies of Transverse Momentum Dependent Parton Distributions and Bessel Weighting

In this paper we present a new technique for analysis of transverse momentum dependent parton distribution functions, based on the Bessel weighting formalism. The procedure is applied to studies of the double longitudinal spin asymmetry in semi-inclusive deep inelastic scattering using a new dedicated Monte Carlo generator which includes quark intrinsic transverse momentum within the generalized parton model. Using a fully differential cross section for the process, the effect of four momentum conservation is analyzed using various input models for transverse momentum distributions and fragmentation functions. We observe a few percent systematic offset of the Bessel-weighted asymmetry obtained from Monte Carlo extraction compared to input model calculations, which is due to the limitations imposed by the energy and momentum conservation at the given energy/Q2. We find that the Bessel weighting technique provides a powerful and reliable tool to study the Fourier transform of TMDs with controlled systematics due to experimental acceptances and resolutions with different TMD model inputs.


Introduction
The study of the spin structure of protons and neutrons is one of the central issues in hadron physics, with many dedicated experiments, recent (HERMES at DESY, CLAS and Hall-A at JLAB), running (COMPASS at CERN, STAR and PHENIX at RHIC), approved (JLab 12 GeV upgrade [1], COMPASS-II [2]) or planned (Electron Ion Collider [3][4][5]). The Transverse Momentum Dependent (TMD) parton distribution functions and fragmentation functions play a crucial role in gathering and interpreting information of a true "3-dimensional" imaging of the nucleon. These Transverse Momentum Dependent distribution and fragmentation functions (collectively here called "TMDs") can be accessed in several types of processes, one of the most important is single particle hadron production in Semi-Inclusive Deep Inelastic Scattering (SIDIS) of leptons on nucleons. A significant amount of data on spin-azimuthal distributions of hadrons in SIDIS, providing access to TMDs has been accumulated in recent years by several collaborations, including HERMES, COMPASS and Halls A,B and C at JLab [6][7][8][9][10][11][12][13][14][15]. At least an order of magnitude more data is expected in coming years of running of JLab 12 [1].
A rigorous basis for studies of TMDs in SIDIS is provided by TMD factorization in QCD, which has been established in Refs. [16][17][18][19][20][21][22][23] for leading twist single hadron production with transverse momentum of the produced hadron being much smaller than the hard scattering scale, and the order of Λ QCD , that is Λ 2 QCD < P 2 h⊥ Q 2 . In this kinematic domain the SIDIS cross section can be expressed in terms of structure functions encoding the strong-interaction dynamics of the hadronic sub-process γ * + p → h + X [24][25][26][27], which are given by convolutions of a hard scattering cross section and TMDs. However the extraction of TMDs as a function of the light-cone fraction x and transverse momentum k ⊥ from single and double spin azimuthal asymmetries is hindered by the fact that observables are complicated convolutions in momentum space making the flavor decomposition of the underlying TMDs a model dependent procedure.
Based on TMD factorization theorems, experimentally measured cross sections are expressed as convolutions of TMDs where k ⊥ dependence is integrated over and related to measured value of P h⊥ . A reliable method to directly access the k ⊥ dependence of TMDs is very desirable. However, various assumptions involved in modern extractions of TMDs from available data rely on conjectures of the transverse momentum dependence of distribution and fragmentation functions [28][29][30][31][32][33][34][35][36][37][38] making estimates of systematic errors due to those assumptions extremely challenging.
In a paper by Boer, Gamberg, Musch, and Prokudin [39], a new technique has been proposed called Bessel weighting, which relies on a model-independent deconvolution of structure functions in terms of Fourier transforms of TMDs from observed azimuthal moments in SIDIS with polarized and unpolarized targets. In this paper, we apply the Bessel weighting procedure to present an extraction of Fourier transforms of TMDs from a Monte Carlo event generator. As an application of this procedure we consider the ratio of helicity, g 1L and unpolarized, f 1 TMDs from the double longitudinally polarization asymmetry. This paper is organized as follows: We begin our discussion in Section 2 with a brief review of the formalism of the SIDIS cross section and its representation in both momentum and Fourier conjugate b T space [39], while also presenting a description of the experimental procedure to study the TMDs based on the Bessel weighting. In Section 3 we introduce a fully differential Monte Carlo generator which has been developed to test the procedure for extraction of TMDs from SIDIS. As a test of the quality of our constructed Monte Carlo, in Section 3.2 we present a study of the Cahn effect [40,41] contribution to the average cos φ moment in SIDIS. In Section 4 we present the extraction of the double spin asymmetry A LL (b T ), defined as the ratio of the difference and the sum of electro-production cross sections for anti-parallel and parallel configurations of lepton and nucleon spins using the Bessel weighting procedure. The effects of different model inputs and experimental resolutions and acceptances on extracted TMDs are investigated. Finally in Section 5 we draw some conclusions of the present analysis and outline steps for future work.

Extraction of TMDs using Bessel Weighting
The SIDIS cross section can be expressed in a model independent way in terms of a set of 18 structure functions [24,25,27,[42][43][44], where the first two subscripts of the structure functions F XY indicate the polarization of the beam and target, and in certain cases, a third sub-script in F XY,Z indicates the polarization of the virtual photon. The structure functions depend on the the scaling variables x, z, the four momentum Q 2 = −q 2 , where q = l − l is the momentum of the virtual photon, and l and l are the 4-momenta of the incoming and outgoing leptons, respectively. P h⊥ is the transverse momentum component of the produced hadron with respect to the virtual photon direction. The scaling variables have the standard definitions, x = Q 2 /2(P · q), y = (P · q)/(P · l), and z = (P · P h )/(P · q). Further, in Eq. (2.1) α is the fine structure constant; the angle ψ is the azimuthal angle of around the lepton beam axis with respect to an arbitrary fixed direction [44], and φ h is the azimuthal angle between the scattering plane formed by the initial and final momenta of the electron and the production plane formed by the transverse momentum of the observed hadron and the virtual photon, whereas φ S is the azimuthal angle of the transverse spin in the scattering plane [45]. Finally, ε is the ratio of longitudinal and transverse photon fluxes [27].
At tree-level, in a parton model factorization framework [25,27,43], the various structure functions in the cross section are written as convolutions of the TMDs which relate transverse momenta of the active partons and produced hadron. For our purposes, the unpolarized and double longitudinal polarized structure functions are where k ⊥ is the intrinsic transverse momentum of the struck quark, and p ⊥ is the transverse momentum of the final state hadron relative to the fragmenting quark k (see Fig. 1). f a 1 (x, k 2 ⊥ ), g 2 1L and D a (z, p 2 ⊥ ) represent TMD PDFs and fragmentation functions respectively of flavor a, e a is the fractional charge of the struck quark or anti-quark and the summation runs over quarks and anti-quark flavors a.
Measurements of transverse momentum P h⊥ of final state hadrons in SIDIS with polarized leptons and nucleons provide access to transverse momentum dependence of leading twist TMDs. Recent measurements of multiplicities and double spin asymmetries as a function of the final transverse momentum of pions in SIDIS at COMPASS [46], HERMES [47], and JLab [13][14][15] suggest that transverse momentum distributions depend on the polarization of quarks and possibly also on their flavor [38] (see also discussion in Ref. [48]). Calculations of transverse momentum dependence of TMDs in different models [49][50][51][52] and on the lattice [53,54] also indicate that the dependence of transverse momentum distributions on the quark polarization and flavor maybe significant. Larger intrinsic transverse momenta of sea-quarks compared to valence quarks have been discussed in an effective model of the low energy dynamics resulting from chiral symmetry breaking in QCD [55].
As stated above, the various assumptions on transverse momentum dependence of distributions on spin and flavor of quarks however make phenomenological fits very challenging. To minimize these model assumptions, Kotzinian and Mulders [56] suggested using so called P h⊥ -weighted asymmetries, where the unknown k ⊥ -dependencies of TMDs are integrated out, thus providing access to moments of TMDs. However, the P h⊥ -weighted asymmetries introduce a significant challenge to both theory and experiment. For example, the weighting with P h⊥ emphasizes the kinematical region with higher P h⊥ , where the statistics are poor and systematics from detector acceptances are difficult to control and at the same time theoretical description in terms of TMDs breaks down.
The method of Bessel weighting [39] addresses these experimental and theoretical issues. First, Bessel weighted asymmetries are given in terms of simple products of Fourier transformed TMDs without imposing any model assumptions of the their transverse momentum dependence. Secondly, Bessel weighting regularizes the ultraviolet divergences resulting from unbound momentum integration that arises from conventional weighting. Further, in this paper we will demonstrate that they provide a new experimental tool to study the TMD content to the SIDIS cross section that minimize the transverse momentum model dependencies inherent in conventional extractions of TMDs. Also they suppress the kinematical regions where cross sections are small and statistics are poor [39].
We begin the review by re-writing the SIDIS cross section as a Bessel weighted integral in b T space [39]: where the structure functions F XY,Z are now given as simple products of Fourier Transforms of TMDs.
Here we consider the unpolarized and double longitudinal structure functions, where the Fourier transform of the TMDs are defined as

Bessel Weighting of Experimental Observables
It is now straightforward to express Bessel weighting of experimental observables. They are quantities which can be presented as simple products of Fourier transforms of distribution and fragmentation functions, allowing the application of standard flavor decomposition procedures. Noting that one can project out the unpolarized and double longitudinally polarized structure functions F LL , and F U U,T , by integrating Eq. (2.4) with the zeroth order Bessel function J 0 (|b T ||P h⊥ |) over the transverse momentum of the produced hadron P h⊥ , we arrive at an expression for the longitudinally polarized cross sectionσ where dΦ ≡ dx dy dψ dz dP h⊥ P h⊥ represents shorthand notation for the phase space differential and |b T | ≡ b T , and |P h⊥ | ≡ P h⊥ . Now we form the double longitudinal spin asymmetry Note that in our definition b T is Fourier conjugate variable to P h⊥ [39].
The experimental procedure to study the structure functions in b T -space amounts to discretizing Eq. (2.8). Now Eq. (2.9) results in an expression of sums and differences of Bessel functions for a given set of experimental events. The details on this formulation are given in Appendix A. The resulting expression for the spin asymmetry is where j ± indicates a sum on events and where N ± is the number of events with positive/negative products of lepton and nucleon helicities.
The cross sectionsσ ± (b T ) can be extracted for any given b T using sums over the same set of data. These cross sections contain the same information as the cross sections, dσ/dΦ in Eq. (2.8) differential with respect to the outgoing hadron momentum. The momentum dependent and the b T -dependent representations of the cross section are related by a 2-D Fourier-transform in cylinder coordinates.
In the next Section we describe a new dedicated Monte Carlo generator which includes quark intrinsic transverse momentum within the generalized parton model.

Fully Differential Monte Carlo for SIDIS
A Monte Carlo generator is a crucial component in testing different experimental procedures. The Monte Carlo generator we use was developed to study partonic intrinsic motion using the framework of the so-called generalized parton model described in detail in Ref. [29]. We consider the SIDIS process where is the incident lepton, N is the target nucleon, and h represents the observed hadron, and the four-momenta are given in parenthesis. Following the Trento conventions [45], the spatial component of the virtual photon momentum q is along the positive z direction and the proton momentum P is in the opposite direction, as depicted in Fig. 1. In the parton model, the virtual photon scatters off an on-shell quark where the initial quark momentum k and scattered quark momentum k have the same intrinsic transverse momentum component k ⊥ with respect to the z axis, and where the initial quark has the fraction x of the proton momentum. The produced hadron momentum, P h has the fraction z of scattered quark momentum k in the (x,ỹ,z) frame and p ⊥ is the transverse momentum component with respect to the scattered quark k (see also, Appendix C). A great deal of phenomenological effort (see for example [29,34,57]) has been devoted to using the generalized parton model, with intrinsic quark transverse momentum, to account for experimentally observed spin and azimuthal asymmetries as a function of the produced hadron's transverse momentum P h⊥ in SIDIS processes. In order to take into account non-trivial kinematic effects arising from the standard approximations [25,27] we develop a Monte Carlo based on the fully differential SIDIS cross section [29] which is given by, where the summation runs over quarks flavors, and the kinematic factors K(x, y) and ε, and the Jacobian J(x, Q 2 , k ⊥ ) are defined in Appendix C. λ is the product of target polarization and beam helicity (λ = ±1), f 1,q (x, k 2 ⊥ ) and g 1L,q (x, k 2 ⊥ ) are the unpolarized and helicity TMDs , and D 1,q (z, p 2 ⊥ ) is the unpolarized fragmentation function, φ l is the scattered lepton azimuthal angle 1 . We adopt the parton kinematics in [29,58] with the additional requirements, that the kinematics of the initial and final parton momenta are kept exact [59], and the nucleon mass is not set to zero. Also the hard scattering matrix elements are calculated for on-shell scattered partons.
In the Monte Carlo generator software, we used the general-purpose, self-adapting event generator, Foam [60], for drawing random points according to an arbitrary, userdefined distribution in n-dimensional space. Figure 1. Kinematics of the process. q is the virtual photon, k and k are the initial and struck quarks, k ⊥ is the quark transverse component. P h is the final hadron with a p ⊥ component, transverse with respect to the fragmenting quark k direction.

Kinematical Distributions
Implementing the Monte Carlo we generate kinematical distributions in x, z, k ⊥ , and p ⊥ of SIDIS events for several model inputs of TMDs. These distributions are then used to check the consistency of dependence of extracted quantities under different model assumptions, including, for example Gaussian and non-Gaussian distributions in transverse momentum.
In case the dependence is assumed to be a Gaussian, x and z dependent widths are assumed, so that TMDs take the following form, where f (x) and D(z) are corresponding collinear parton distribution and fragmentation functions and the widths are x and z dependent functions. In our studies we adopt the modified Gaussian distribution functions and fragmentation functions from Eq. (3.3)-(3.5), in which x and k ⊥ dependencies are inspired by AdS/QCD results [61,62], with k 2 where the constants C and D may be different for different flavors and polarization states (see for example [38]). Similarly such non-factorized x,k ⊥ distribution functions are also suggested by the diquark spectator model [63] and the NJL-jet model [36,64].
For the x and z dependence in Eqs. (3.3) and (3.5) we use the parametrizations, using input values C = 0.54 GeV 2 and D = 0.5 GeV 2 . We also assume that k 2 ⊥ g 1L = 0.8 k 2 ⊥ f 1 ; this assumption is consistent with lattice studies [54] and experimental measurements [14].
As an example of a non-Gaussian k ⊥ distribution we implement the following one inspired by the shape of the resulting distribution in the light-cone quark model [65,66] where the coefficients for g 1L (x, k 2 ⊥ ) are chosen in such a way that effectively We then generate events using the cross section from Eq. (3.2) for both Gaussian and non-Gaussian initial distributions respectively, and we show the resulting transverse momentum distributions in Figs. 2 and 3. Note that the generator we construct is implemented with on mass-shell partons and with four momentum conservation imposed. While this choice is not compulsory we adopt it as it allows us to fully reconstruct kinematics for a given event. At the same time limitations due to available phase space integration will modify the reconstructed distributions with respect to the input distributions. We analyze the effect of the available phase space in the Monte Carlo on the average k 2 ⊥ for finite beam energies as a function of x by calculating the effective k 2 ⊥ from the following formula, where the index j runs over the N Monte Carlo generated events. Note, dσ M C is the cross section of the Monte Carlo simulation, that is Eq.  Fig. 3 we study the effect of the non-Gaussian distribution Eq. (3.6). Integrating Eq. (3.7) over k ⊥ gives a value of k 2 ⊥ = 0.084 GeV 2 , and the dashed curve represents the fit to the Monte Carlo distribution with a value of k 2 ⊥ = 0.064 GeV 2 for the 6 GeV initial lepton beam energy.
In Fig. 4. the average k 2 ⊥ versus x from the Monte Carlo for different incoming beam energies, for 0.5 < z < 0.52, is presented. For the modified Gaussian distribution function with the input value k 2 ⊥ (x) = 0.54 x(1−x) GeV 2 , the suppression of the generated k 2 ⊥ (x) compared to input distributions (solid line) is greater for the lower beam energy. In Fig. 5 the constraints of four momentum conservation also affect the p 2 ⊥ distributions, which in turn also affect the observed P h⊥ distribution.
The systematics of the extraction of the TMDs in momentum space due to the kinematic constraints has been studied in detail using our fully differential Monte Carlo. We conclude this section with the general observation that imposing four momentum conservation in the event generator effectively modifies the initial distributions due to the limitations of the available phase space in the generator. This deformation is more pronounced at lower energies or Q 2 . A shift of a few percent is visible for 160 GeV cm energy, while for the lower 6 GeV cm energy the effective k 2 ⊥ is lower than the input value by approximately ∼ 20%.

The Cahn effect in the Monte Carlo Generator
As an example of an application of our constructed Monte Carlo we present a study of the Cahn effect [40,41] contribution to the average cos φ moment in SIDIS. We generate Monte Carlo events using the following expression for the cross section [29], Fig. 1). As stated above, in the Monte Carlo we impose four momentum conservation with target mass corrections. In Fig. 6 we present output from the Monte Carlo using the non-factorized Gaussian distribution function and fragmentation function (Eqs. (3.3) and (3.5)). We also compare our results to the HERMES data, and Ref. [58]. It is clear that the results of our Monte Carlo are comparable to that of [58] and close to HERMES data. For the red triangles we used k 2 ⊥ = 0.54 x(1 − x) and p 2 ⊥ = 0.5 z(1 − z) GeV 2 . As one can see for HER-MES kinematics the modified Gaussian TMDs reduces the contribution of the Cahn effect contribution to the cos(φ h ) moment. In Ref. [58] this effect is achieved by imposing a so-called direction cut (that the quark moves in the forward direction with respect of the proton). In this Monte Carlo there are two main factors that modify the distribution; the four-momentum conservation and x(z) dependent values of k 2 ⊥ ( p 2 ⊥ ). In the next Section we apply the Bessel weighting formalism for the double longitudinal spin asymmetry in semi-inclusive deep inelastic scattering to data from our Monte Carlo generator.

Bessel Weighted Double Spin Asymmetry
In this Section, we present an extraction of the Bessel weighted double longitudinal spin asymmetry in b T space. We also carry out a study of the accuracy of such an extraction. We use the dedicated fully differential SIDIS single hadron Monte Carlo to generate events based on the input TMDs. For simplicity we perform this comparison in a one flavor approximation.

Results from the Monte Carlo
The Monte Carlo generated events are used like experimental events to extract both the Bessel weighted asymmetry, A J 0 (b T P h⊥ ) LL , and the ratio of the Fourier transform of g 1L and f 1 using the Bessel weighting method described in [39]. The results are then compared to the Monte Carlo input. The Bessel moments are extracted from the Monte Carlo with 6 GeV beam energy using both the modified Gaussian type of functions (see Eqs. (3.3)-(3.5)) and power law-tail like function (see Eq. (3.6)).
The numerical results of our studies are summarized and displayed in Figs. 7 and 8 for the modified Gaussian distribution function and for the power law-tail like distribution function inputs respectively. In the left panel of Fig. 7 we show the Bessel-weighted asymmetry versus b T . The blue curve labeled "BW Input", is the asymmetry calculated analytically using the right hand side of Eq. (2.9) and the Fourier transformed input distribution functions (one can compare this with the model calculation in Ref. [67]).
We now compare various distributions generated from the Monte Carlo. We plot the generated distribution using Eq. (2.10) (full red points) labeled "BW(P h⊥ ) Generated", and the black triangles labeled "BW(P h⊥ ) Sm + Acc", which represents the same extraction after experimental smearing and acceptance (using the CLAS detector [68]). Next we consider the Fourier transform ratiog 1L tof 1 , the (green) curve with triangles up labeled "BW(k ⊥ )" obtained from numerically Fourier transforming the k ⊥ distributions from the  Monte Carlo generator on an event by event basis (see Eq. (2.6)), .

(4.1)
This quantity corresponds to the right hand side of Eq. (2.9) in the one flavor approximation, where the fragmentation functions are expected to cancel out. For the Monte Carlo generated events, this cancellation is only an approximation, leading to the deviation between the red (or black) points and the green curve at large b T . The reasons for the imperfect cancellation are discussed in section 4.2.
In order to quantitatively assess the deviation between the curves in the left panel of Fig. 7, we plot ratios of these values (see right panel). The red points represent the deviation from unity that is due to the imperfect cancellation of the fragmentation function. The black triangles represent the same after experimental smearing and acceptance are taken into account. Finally the open blue squares represent the deviation between the analytic result from the input distributions and the Monte Carlo generated events, Eq. (2.10).
The error bars in b T space for each point give the statistical standard deviation (see Appendix B for more details). If we use the same data set to integrate over P h⊥ for all b T points, the errors in b T space are correlated and one should calculate the correlation/covariance matrix before performing any global fit to the data points. To circumvent this problem, we used different Monte Carlo samples for each b T -point to avoid correlations in the error calculations.

Interpretation of the Results
One primary question addressed in this study is how robust the Bessel-weighting technique behaves under simulating real experimental conditions. Comparing the round (red) data points with the triangular (black) ones in Figs. 7 and 8, we see that switching on experimental smearing and acceptance in our simulation does not change the results significantly.
An altogether different question concerns the validity of the generalized parton model at the relatively low beam energies available in experiments today. The generalized parton model is an approximation that assumes certain components of the intrinsic parton momenta are suppressed for large beam energies and can thus be integrated out from the distributions. This becomes apparent in Eqns. (2.2) and (2.3), where a delta function is present only in the two transverse dimensions. An explicit four-momentum conservation law embedded in the formula of the cross section is thus lost. As an example of a particularly striking consequence of this, one observes that there is no explicit mechanism that prevents events at values of P h⊥ larger than allowed by the finite beam energy. Naturally, the lower the beam energy becomes, the more serious these inaccuracies of the generalized parton model have to be taken.
Our Monte Carlo simulation allows us to take the factorized form of the generalized parton model cross section Eq. (3.2) as a basis and then to impose full four-momentum conservation for the partons, assuming the initial quarks are on-shell with non-zero mass. We also take a non-zero target mass into account. This procedure does not necessarily lead to a more accurate description of the underlying physics, because it still rests on the simplified picture of the generalized parton model and involves the approximation of an on-shell quark. Nonetheless, implementing these modifications can serve us to get an indication for the magnitude of uncertainties resulting from the aforementioned kinematic approximations in the generalized parton model. Analyzing our MC results with four momenta conservation and target mass correction, we are able to distinguish two effects in the left panels of Figs. 7 and 8: 1. Solid (blue) curve versus triangular (green) data points: The distributions realized in the MC simulation differ from the input distributions. In the MC, the four-momentum conservation does not allow the variables k ⊥ and p ⊥ in Eq. (3.2) to be sampled independently over the whole integration range, as it would have to be done to reproduce the unmodified generalized parton model Eqns. 2. Triangular (green) data points vs. circular (red) data points: inadequacy of the generalized parton model to describe the data. In a single flavor scenario, the distribution functionsD a 1 cancel exactly on the right hand side of Eq. (2.9). Therefore, there should not be any difference between the full asymmetry A J 0 (b T P h⊥ ) LL (b T ) of Eqs. (2.9), (2.10) and the ratio of TMD PDFs Eq. (4.1). However, we do observe a difference between the circular (red) data points and the triangular (green) data points in the left panels of Figs. 7 and 8. Again, the four-momentum conservation we have implemented is the reason for the observed difference. Since k ⊥ and p ⊥ are no longer sampled in accordance with Eqns. (2.2) and (2.3), the right hand side of Eq. (2.9) looses the prerequisites for its derivation and is violated to some degree. Therefore, we see only an incomplete cancellation of FFs for the Monte Carlo events.
To an experimentalist who is concerned about systematic errors attributed to the observables he or she extracts, the first of the two effects above is not an issue. The purpose of the generalized parton model is to provide a parametrization of the data one observes. Any effect of the underlying scattering mechanism that can be absorbed into the distributions does not contradict the validity of the model. The only concern one might have is that the distributions become beam energy/Q 2 dependent, an issue that should be addressed using TMD evolution equations.
On the other hand, the second effect presented above can be taken as an indication for systematic uncertainties. If, indeed, the physical reality does not generate events in accordance with the functional shape of the generalized parton model, then using the model for the extraction of distributions necessarily involves systematic errors. Again, we point out that it is unclear whether the modifications we have implemented in our MC bring us closer to the physical reality. Nonetheless, the modifications are reasonable and so we believe they can give us a hint about the order of magnitude of systematic errors from the corresponding approximations in the model. One can then estimate that for calculations such as those performed in Ref. [67], systematic errors in the comparison with experimental data for b T < 6 GeV −1 are of the order of a few percent. For the data with b T > 6 GeV −1 , the effects of four-momentum conservation (difference between red and green points) becomes more pronounced, and a fit of data using the generalized parton model without manifest four-momentum conservation therefore becomes less accurate.

Conclusions
We have presented the first studies of Bessel-weighted asymmetries using a multi-dimensional Monte Carlo generator based on the fully differential cross section for TMD studies using the tree level parton model [29]. Two models have been used in the simulation; a modified Gaussian and a power law tail, for the distribution and fragmentation functions. The Bessel-weighted sums of double polarization observables, in particular, provide access to transverse momentum dependencies of partonic distributions f 1 and g 1L . Bessel-weighted asymmetries (described in [39]) have been extracted from the generated Monte Carlo events and studies of systematic uncertainties have been performed. We observe a few percent systematic offset of the Bessel-weighted asymmetry obtained from Monte Carlo extraction compared to input model calculations, which is due to the limitations imposed by the energy and momentum conservation at the given energy/Q 2 .
We find that the Bessel weighting technique provides a powerful and reliable tool to study the Fourier transform of TMDs with controlled systematics due to experimental acceptances and resolutions with different TMD models inputs. We plan to expand our studies with more advanced parton shower and fragmentation mechanisms, as well as to include nuclear modifications in our Monte Carlo and extraction procedure.
A Monte Carlo generator including spin-orbit correlations, quark-gluon interactions and correlations between the current and target fragmentation region, which is applicable in a wide range of kinematics, will be crucial for both experimental techniques and phenomenology of Fourier transformed TMDs. Moreover, evolution equations for the distributions are typically formulated directly in coordinate (Fourier) space [16][17][18][19][20][21][22][23]. Phenomenological studies can then performed in this space, see for example [69,70]. Thus, the study of the scale dependence of Bessel weighted asymmetries should prove important in studies of evolution of TMDs. For the above stated reasons we propose Bessel weighted asymmetries as clean observables to study the scale dependence of TMD PDFs and FFs at existing (HERMES, COMPASS, JLab) and future facilities (Electron Ion Collider, JLab 12 GeV).

A Bessel Weighting
In this Appendix we review the Bessel weighting framework, and the procedure to calculate the Bessel-weighted asymmetry for the longitudinally polarized beam and target, for a given set of experimental events which is expressed in Eq. (2.10).
From Eq. (2.4) the SIDIS cross section written in terms of the Fourier transformed TMD PDFs and FFs [39] for the leading twist unpolarized and doubly longitudinal polarized structure functions is given by where K(x, y) is given below in Eq. (C.2) and where |P h⊥ | ≡ P h⊥ and |b T | ≡ b T .
Using the Bessel weighting procedure, which in this case amounts to weighting with J 0 , we write the cross sectionσ(B T ) in B T space, in terms of the structure functions 2 F U U,T and F LL σ(B T ) = 2π dP h⊥ P h⊥ J 0 (B T P h⊥ ) dσ dx dy dψ dz dφ h dP h⊥ P h⊥ where the structure functions in b T space are given by the products of Fourier transformed TMDs [39], Labeling the cross section with ± for the double longitudinal spin combinations S || λ e = ±1 we haveσ The Bessel weighted double spin asymmetry is b T space is, .

(A.5)
Now we derive the formula to extract Bessel-weighted asymmetries by means of an event by event weighting in P h⊥ , while binning in x, y, and z. First we express the unpolarized and doubly polarized helicity structure functions in B T space as using the shorthand notation for the differential phase space factor dΦ ≡ dx dy dψ dz dP h⊥ P h⊥ . Re-expressing the cross sections in terms of the number of events in the differential phase space "volume", Eq. (A.6) is given by, and where dn ± are the number of events in a differential phase space volume, dΦ, and N ± 0 is the standard normalization factor, that is the product of the number of beam and target particles with ± polarization per unit target area. Now we discretize the momentum integration in Eq. (A.7) and (A.8) for a fixed phase space cell in x, y, z such that the corresponding differential dx dy dz becomes the bin volume ∆x∆y∆z. Eqs. (A.7) and (A.8) thus become (A.10) where we sum over the discrete momentum index i, and ∆n ± i are the number of events for polarization ± as a function of P h⊥i .
Substituting Eqs. (A.9) and (A.10) into Eq.(A.4), the experimental procedure to calculate the Bessel weighted asymmetry, A J 0 (b T P h⊥ ) LL (b T ), becomes, where j ± indicates a sum on events and where N ± is the number of events with positive/negative products of lepton and nucleon helicities for a given x, y and z, and wherẽ S ± is the sum over events for ± helicities.

B Error calculations
To get an error estimate for the binned Bessel-weighted sums it is straightforward to show that the error for the Bessel-weighted sums are, is given by Eq. (A.11).

C Equations in the Monte Carlo
Eq. 3.2 in the Monte Carlo is implemented in the form, dσ dxdydz LC dp 2 ⊥ dk 2 ⊥ dφ l dφ k dφ where the summation runs over quarks flavors and, K(x, y) = α 2 xyQ 2 y 2 2(1 − ε) 1 + γ 2 2x , ε = 1 − y − 1 4 γ 2 y 2 1 − y + 1 2 y 2 + 1 4 γ 2 y 2 , (C.2) and the Jacobian J is given by In above equations x -is the Bjorken variable, while x LC = k − /P − is the light-cone (LC) fraction of the proton momentum carried by the quark k [29]. The quark four momentum is given by, where k 0 is the quark energy, and k {x,y,z} are the x, y and z components of the quark momentum in the CM frame of virtual photon and proton. P = 0.5(E p + |P pz |), where proton energy in the CM is E p = P 2 pz + M 2 . Taking into account the nucleon mass, the on-shell condition for the final quark implies where k ⊥ is the parton transverse momentum. The scattered quark momentum k is constructed using k = k + q (see Fig 1) and p ⊥ is the transverse momentum of the hadron with respect to the scattered quark k . φ k is the initial quark azimuthal angle Fig 1. z LC = P + h /k + is the light-cone fraction of the quark momentum carried by the resulting hadron in the (x,ỹ,z)-system [29], wherez is aligned along the scattered quark k Fig 1. The final hadron momentum is constructed using, wereφ is the angle between quark and hadron planes. φ h is the angle between leptonic and hadronic planes according to the Trento convention and P h⊥ is the final hadron transverse momentum [45]. The final hadron SIDIS variables φ h , P h⊥ and z are calculated after event generation. Here we should note, that theoretical or phenomenological DFs and fragmentation functions are expected to be in the light cone coordinate system (see Eq. C.1). Motivated by the fact that x LC x and z LC z is a widely used approximation we implement f 1 (x) and D 1 (z) instead of f 1 (x LC ) and D 1 (z LC ) in Monte Carlo.