Efficient Fourier Transforms for Transverse Momentum Dependent Distributions

TL;DR

This paper introduces a fast, accurate numerical Hankel transform algorithm for TMD calculations in particle physics, improving computational efficiency and precision across various processes.

Contribution

The paper develops and provides open-source implementations of a novel Hankel transform algorithm tailored for TMD phenomenology, enhancing numerical accuracy and speed.

Findings

Improved numerical accuracy in TMD calculations.

Open-source libraries in multiple programming languages.

Enhanced computational efficiency for Fourier transforms.

Abstract

Hadron production at low transverse momenta in semi-inclusive deep inelastic scattering can be described by transverse momentum dependent (TMD) factorization. This formalism has also been widely used to study the Drell-Yan process and back-to-back hadron pair production in $e^+e^-$ collisions. These processes are the main ones for extractions of TMD parton distribution functions and TMD fragmentation functions, which encode important information about nucleon structure and hadronization. One of the most widely used TMD factorization formalism in phenomenology formulates TMD observables in coordinate $b_\perp$-space, the conjugate space of the transverse momentum. The Fourier transform from $b_\perp$-space back into transverse momentum space is sufficiently complicated due to oscillatory integrands that it requires a careful and computationally intensive numerical treatment in order to avoid potentially large numerical errors. Within the TMD formalism, the azimuthal angular dependence is analytically integrated and the two-dimensional $b_\perp$ integration reduces to a one-dimensional integration over the magnitude $b_\perp$. In this paper we develop a fast numerical Hankel transform algorithm for such a $b_\perp$-integration that improves the numerical accuracy of TMD calculations in all standard processes. Libraries for this algorithm are implemented in Python 2.7 and 3, C++, as well as FORTRAN77. All packages are made available open source.

Equations

$$\frac{d{\sigma }^h}{d{x}_{\mathrm{bj}}dydz{d}^2{q}_{\perp }}=\frac{2\pi {\alpha }_{\mathrm{EM}}^2}{Q^2}\frac{1+(1-y{)}^2}{y}W(q_{\perp },x_{\mathrm{bj}},z,Q),$$

$$q=\ell -{\ell }^{\prime },Q^2=-q^2,x_{\mathrm{bj}}=\frac{Q^2}{2P\cdot q},y=\frac{P\cdot q}{P\cdot \ell },z=\frac{P\cdot P_h}{P\cdot q}.$$

$$W(q_{\perp },x_{\mathrm{bj}},z,Q)=$$

$$\tilde{W}(b_{\perp },x_{\mathrm{bj}},z,Q)=\int d^2{\mathbf{q}}_{\perp }{e}^{i{\mathbf{q}}_{\perp }\cdot {\mathbf{b}}_{\perp }}W(q_{\perp },x_{\mathrm{bj}},z,Q),$$

$$W(q_{\perp },x_{\mathrm{bj}},z,Q)=$$

$$=$$

$${\int }_0^{2\pi }\frac{d\varphi }{2\pi }{e}^{-i{q}_{\perp }{b}_{\perp }\cos (\varphi )}=J_0(q_{\perp }{b}_{\perp }),$$

$${\mathbf{q}}_{\perp }^{\alpha }W(q_{\perp },\cdots ),{\mathbf{q}}_{\perp }^{\alpha }{\mathbf{q}}_{\perp }^{\beta }W(q_{\perp },\cdots )$$

$${\int }_0^{\infty }\frac{d{b}_{\perp }{b}_{\perp }^{\nu +1}}{2\pi }{J}_{\nu }(q_{\perp }{b}_{\perp })\tilde{W}(b_{\perp }),$$

$${\int }_0^{\infty }dxf(x)J_n(x),$$

$${\int }_{-\infty }^{\infty }dx|x{|}^{2n+1}f(x)=h\sum _{j=-\infty ,j\ne 0}^{\infty }{w}_{nj}|x_{nj}{|}^{2n+1}f(x_{nj})+\mathcal{O}\left(e^{-c/h}\right),$$

$$x_{nj}=h{\xi }_{nj},w_{nj}=\frac{2}{{\pi }^2{\xi }_{n|j|}{J}_{n+1}(\pi {\xi }_{n|j|})},$$

$${\int }_0^{\infty }dx{x}^{2n+1}f(x)=h\sum _{j=1}^{\infty }{w}_{nj}{x}_{nj}^{2n+1}f(x_{nj})+\mathcal{O}\left(e^{-c/h}\right).$$

$${\mathcal{I}}_{nN+1}=h\sum _{j=N+1}^{\infty }{w}_{nj}{x}_{nj}^{2n+1}f(x_{nj})$$

$${\int }_0^{\infty }dx{x}^{2n+1}f(x)=h\sum _{j=1}^N{w}_{nj}{x}_{nj}^{2n+1}f(x_{nj})+\left[{\mathcal{I}}_{nN+1}+\mathcal{O}\left(e^{-c/h}\right)\right].$$

$$\displaystyle x=\frac{\pi}{h}\psi(t)\qquad\mbox{with \; $\psi(t)=t\tanh\left(\frac{\pi}{2}\sinh t\right)$}\;\;,$$

$${\int }_0^{\infty }dxf(x)J_n(x)=$$

$$=$$

$${\int }_0^{\infty }dxf(x)J_n(x)=\pi \sum _{j=1}^N{w}_{nj}f\left(\frac{\pi }{h}\psi (x_{nj})\right)J_n\left(\frac{\pi }{h}\psi (x_{nj})\right){\psi }^{\prime }(x_{nj})+\left[{\mathcal{I}}_{nN+1}^{\prime }+\mathcal{O}\left(e^{-c/h}\right)\right],$$

$${\mathcal{I}}_{nN+1}^{\prime }=\pi \sum _{j=N+1}^{\infty }{w}_{nj}f\left(\frac{\pi }{h}\psi (x_{nj})\right)J_n\left(\frac{\pi }{h}\psi (x_{nj})\right){\psi }^{\prime }(x_{nj}).$$

$$\frac{\pi }{h}\psi (x_{nj})\approx \pi {\xi }_{nj}\left[1-2\exp \left(-\frac{\pi }{2}{e}^{x_{nj}}\right)\right],$$

$$J_n\left(\frac{\pi }{h}\psi (x_{nj})\right)\approx 2\pi {\xi }_{nj}{J}_{n+1}(\pi {\xi }_{nj})\exp \left(-\frac{\pi }{2}{e}^{x_{nj}}\right).$$

$$\frac{\partial }{\partial h}\left(h(h{\xi }_{n1}{)}^{2n+1}f(h{\xi }_{n1})\right)=0.$$

$$h_{\mathrm{u}}{\xi }_{nN}=\frac{\pi }{h}\psi (h{\xi }_{nN}).$$

$$h=\frac{1}{{\xi }_{nN}}{\sinh }^{-1}\left(\frac{2}{\pi }{\tanh }^{-1}\left(\frac{h_{\mathrm{u}}}{\pi }\right)\right),$$

$$W(q_{\perp })={\int }_0^{\infty }\frac{d{b}_{\perp }{b}_{\perp }}{2\pi }\tilde{W}(b_{\perp })J_0(b_{\perp }{q}_{\perp }),$$

$$\tilde{W}(b_{\perp },\beta ,\sigma )=\frac{1}{b_{\perp }}{\left(\frac{\beta b_{\perp }}{{\sigma }^2}\right)}^{{\beta }^2/{\sigma }^2}\frac{e^{-\frac{\beta b_{\perp }}{{\sigma }^2}}}{\Gamma \left(\frac{{\beta }^2}{{\sigma }^2}\right)}.$$

$$W^{\mathrm{exact}}(q_{\perp },\beta ,\sigma )=\frac{1}{2\pi }\left(\frac{{\sigma }^2}{\beta }\right)\frac{\Gamma \left(\frac{{\beta }^2}{{\sigma }^2}+1\right)}{\Gamma \left(\frac{{\beta }^2}{{\sigma }^2}\right)}{}_2{\tilde{F}}_1\left(\frac{1}{2}\left(\frac{{\beta }^2}{{\sigma }^2}+1\right),\frac{1}{2}\left(\frac{{\beta }^2}{{\sigma }^2}+2\right);1;-\frac{q_{\perp }^2{\sigma }^4}{{\beta }^2}\right),$$

$$b_{\perp }^{\mathrm{peak}}=\frac{{\beta }^2-{\sigma }^2}{\beta }.$$

$$Q=\frac{1}{b_{\perp }^{\mathrm{peak}}}=\frac{\beta }{{\beta }^2-{\sigma }^2}.$$

$$W^{\mathrm{Ogata}}(q_{\perp },\beta ,\sigma )=\pi \sum _{j=1}^N{w}_{nj}f\left(\frac{\pi }{h}\psi (h{\xi }_{nj})\right)J_n\left(\frac{\pi }{h}\psi (h{\xi }_{nj})\right){\psi }^{\prime }(h{\xi }_{nj})$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: Fast Bessel Transform (FBT)

CPC Library link to program files: (to be added by Technical Editor)

Developer’s respository link: https://github.com/UCLA-TMD/Ogata

Code Ocean capsule: (to be added by Technical Editor)

Licensing provisions: MIT

Programming language: C++, FORTRAN77, Python 2/3

Nature of Problem: In order to perform extractions of transverse momentum distribution functions, numerical Hankel transforms must be performed from $b_{\perp}$-space to momentum space. However, these numerical Hankel transforms are a huge bottleneck in these extractions, making these extractions extremely computationally intensive.

Solution method: We develop a numerical Hankel transform algorithm by optimizing Ogata quadrature formula. This algorithm improves the performance of these numerical Hankel transforms by nearly an order of magnitude.

1 Introduction

The transverse momentum dependent (TMD) parton distribution functions (PDFs) and fragmentation functions (FFs) have received great attention from both theoretical and experimental communities in recent years. These TMD PDFs and FFs, or in short called TMDs, provide new information on hadron structure: the three-dimensional imaging of hadrons in both longitudinal and transverse momentum space [1, 2, 3, 4]. Significant progress has been made in the last few years in terms of measuring transverse momentum dependent unpolarized and polarized cross sections in experiments, as well as extracting the associated spin-independent and spin-dependent TMDs in phenomenology, see Refs. [5, 6] and references therein.

TMDs are non-perturbative objects in Quantum Chromodynamics (QCD) and thus they have to be either computed on the lattice, or extracted from experimental data. For recent developments on lattice computation of TMDs, see Ref. [7]. On the other hand, in order to extract TMDs from the experimental data, one relies on proper QCD factorization theorems [8]. TMD factorization [9, 10, 11, 12] describes cross sections in scattering events where the relevant transverse momentum $q_{\perp}$ of the observed final state is much smaller than the hard scale $Q$: $q_{\perp}\ll Q$. In such a region, the cross section can be factorized in terms of TMD PDFs and/or FFs and perturbatively calculable short distance hard coefficients. The relevant processes that have been extensively studied include semi-inclusive deep inelastic scattering (SIDIS) [13, 14], Drell-Yan type process in proton-proton collisions [15, 16, 17, 18, 19, 20, 21, 22, 23], and back-to-back dihadron production in $e^{+}e^{-}$ collisions [24, 25]. There are also other new opportunities in studying TMDs which are proposed recently in e.g. Refs. [26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43], and usually involve jet measurements.

Within the TMD factorization formalism, the cross section is written as a convolution of the relevant transverse momentum dependent functions. To motivate our discussion and thus make the case more concrete, let us take SIDIS as an example. The differential cross section for the unpolarized scattering process of $e(\ell)+p(P)\to e(\ell^{\prime})+h(P_{h})+X$ can be written as [44]

[TABLE]

where the standard SIDIS variables are defined as

[TABLE]

The unpolarized structure functions $W$ in Eq. (1) can be factorized as follows

[TABLE]

where ${\bf q}_{\perp}=-{\bf P}_{h\perp}/z$, $e_{q}$ is the fractional electric charge for the quarks, and $H(Q,\mu)$ is the hard function to be given by Eq. (35) in Sec. 4. On the other hand, the vectors ${\bf k}_{\perp}$ and ${\bf p}_{\perp}$ are the momentum of the produced quark relative to the parent proton and the momentum of the produced hadron with respect to the fragmenting quark, respectively. The function $f_{q/p}(x_{\rm bj},k_{\perp}^{2})$ is the TMD PDF while $D_{h/q}(z,p_{\perp}^{2})$ is the TMD FF. Here we have suppressed the additional scale parameters in the TMDs, which are associated with QCD evolution of the TMDs [45, 46, 47, 48, 49]. In general, the convolution and integration of TMDs over the momenta ${\bf k}_{\perp}$ and ${\bf p}_{\perp}$ are quite involved. Thus in the original Collins-Soper-Sterman (CSS) approach [50], one takes a Fourier transformation from the momentum space to the coordinate ${\bf b}_{\perp}$ space,  111There are also other approaches in the literature that do not work in the ${\bf b}_{\perp}$ space, see e.g. Refs. [51, 52]. Notice that we drop explicit dependence on $x_{\rm bj}$ and $z$ for the rest of this paper.

[TABLE]

and thus one can write

[TABLE]

where $b_{\perp}=|{\bf b}_{\perp}|$ is the magnitude of the vector ${\bf b}_{\perp}$, $J_{0}$ is the Bessel function of the first kind of order [math], and $f_{q/p}(x_{\rm bj},b_{\perp}),~{}D_{h/q}(z,b_{\perp})$ are the Fourier transform of the TMD PDF and FF, respectively. In going from the first to second line in Eq. (5), we perform analytically the integration over the azimuthal angle $\phi$ ,

[TABLE]

and thus the two-dimensional Fourier transform reduces to a one-dimensional Hankel transform.

For the polarized scattering, the generic structure of the cross sections can be written as [10, 44]

[TABLE]

with $W(q_{\perp},\cdots)$ representing a generic function of $q_{\perp}=|{\bm{q}}_{\perp}|$ and “$\cdots$” denoting the other kinematic variables. The Fourier transform of such functions will lead to Bessel functions of order 1 and 2. In fact, as shown in [53], all the spin-dependent structure functions at leading-power can be expressed in terms of an integration over ${\bm{b}}_{\perp}$ multiplied by the Bessel functions of $J_{0}$, $J_{1}$, or $J_{2}$. Generically the integrals become

[TABLE]

where the function $\widetilde{W}(b_{\perp})$ contains the $b_{\perp}$ space TMD physics.

Without loss of generality the integration form stemming from the Fourier transform encountered in TMD observables can be written as

[TABLE]

where $x=b_{\perp}q_{\perp}$, $J_{n}(x)$ is the Bessel function of order $n$, and $f(x)$ is usually a smooth function of $x$ that slowly decays as $x\to\infty$. Such an integration can be extremely computationally intensive and time consuming with standard integration routines based on adaptive Gaussian quadratures or Monte Carlo integration methods due to the oscillatory nature of the Bessel functions. In the context of TMD global analysis, one has to compute the above integration many times, and for different kinematic regions, in order to find the best fit for the non-perturbative TMDs. This has become a huge hurdle for TMD phenomenology in the past for carrying out the global QCD analysis on TMDs using the data from HERMES, COMPASS, JLAB 6 GeV, Relativistic Heavy Ion Collider (RHIC) and BELLE experiments and it will become even more challenging for the large amount data that is expected from the JLab 12 GeV program and the future Electron Ion Collider (EIC). Because of this, and because of the complexity of the TMD evolution improving the efficiency and the speed of the numerical integration of Eq. (9) is extremely desirable.

Ogata has introduced a quadrature formula in Ref. [54] that is optimized for integrands with Bessel functions for exactly the same integral as in Eq. (9). As will be discussed in Sec. 2, this quadrature method has two parameters, $h$ and $N$, which control the node spacing and the truncation of the quadrature sum, respectively. While this formalism has been previously used in TMD analysis in Ref. [55]; the scheme used in this reference suffers from large numerical inefficiencies. Namely Ref. [55] method takes the parameter $h$ to be held fixed for the entire global analysis, while increasing $N$ to reach convergence of the integration. We will therefore refer to this method as a ‘fixed $h$ Ogata scheme’. In this paper, we propose a new numerical algorithm by which the parameter $h$ is optimized for an input number of nodes $N$, referred to as ‘optimized $h$ Ogata scheme’. Fig. 1 illustrates the difference between these two schemes. We note that while the fixed $h$ Ogata method used in Ref. [55] is sufficient for the extraction of TMDs from several hundred data points, future high precision global extraction of TMDs will utilize tens of thousands data points. Moreover, the reliable estimation of errors of parameters may warrant Monte Carlo fitting methods, such as ones used in JAM15 extraction from Ref. [56].

As the Fourier transforms are a huge bottleneck in these extractions, the numerical inefficiencies of fixed $h$ Ogata methods will be massively amplified to the point where precision extractions may no longer be feasible. We demonstrate in this paper that the optimized Ogata quadrature is more efficient than the traditional adaptive Gaussian quadrature method, Monte Carlo method, as well as fixed $h$ Ogata quadrature method. The rest of this paper is organized as follows. In Sec. 2, we summarize the relevant formalism for the Ogata quadrature method and describe our optimized numerical algorithm in detail. In Sec. 3, we give a demonstration of the optimized Ogata quadrature and benchmark the algorithm against adaptive Gaussian quadrature using an exponential function which has an analytic Fourier transform. In Sec. 4, we apply our numerical method to an example of a phenomenological form of TMDs. We conclude our paper in Sec. 5.

2 Optimized Ogata Quadrature Formalism

In this section we first review the original Ogata quadrature formalism and then discuss our optimization scheme for performing high efficiency numerical integrals relevant to TMD analysis. The Ogata method is based on a quadrature formalism that was first introduced in Ref. [57] by Frappier and Olivier. The quadrature formula for the integrand of the form $|x|^{2n+1}f(x)$ reads:

[TABLE]

where the function $f(x)$ must be an integrable function for the sum to be finite. The nodes $x_{nj}$ and the weights $w_{nj}$ of the quadrature are given by

[TABLE]

with $j=\pm 1,\pm 2,\cdots$, and $\xi_{nj}$ the zeros of the Bessel function $J_{n}(\pi x)$ of order $n$, i.e. $J_{n}(\pi\xi_{nj})=0$, and the parameter $1/h$ represents the node density. The term ${\mathcal{O}}\left(e^{-c/h}\right)$ accounts for the error of the quadrature sum approximation at a finite $h$ as described in equation (2.2) of Ref. [54], and $c$ is a positive constant, whose precise value depends on the functional form of $f(x)$. For the time being, we will be interested in the case of $f(x)$ being an even function of $x$ which results in the following quadrature formula

[TABLE]

In practice the sum in Eq. (12) is truncated at a given $j=N$ which introduces an error of

[TABLE]

and the quadrature formula becomes

[TABLE]

The following change of variables, see Ref. [54], optimizes the convergence of integrals of the typical TMD functional form from Eq. (9):

[TABLE]

so that Eq. (9) becomes

[TABLE]

where $\psi^{\prime}(t)=d\psi(t)/dt$. At this point, it is important to realize that the part of the integrand beside the factor $|t|^{2n+1}$ is an even function of $t$, and thus we can apply Eq. (12) and obtain the following quadrature formula

[TABLE]

where $\mathcal{I}^{\prime}_{n\,N+1}$ are the same truncation errors defined in Eq. (13) but with the transformed integrand,

[TABLE]

Eq. (17) is the aforementioned Ogata quadrature formula, which we advocate in our current paper. The variable substitution has the useful asymptotic behavior

[TABLE]

such that the asymptotic behavior for the Bessel function becomes

[TABLE]

This variable substitution then enforces the double exponential convergence of the quadrature sum in $j$.

The quadrature sum has two parameters, $h$ and $N$, which control the efficiency and the magnitude of the error terms. To generate a high efficiency algorithm, the numerical integration must be performed with a small $N$ while at the same time the error terms must also be small, to ensure reliable results. By inspecting Eqs. (18) and (20) one notes that for a small number of function calls the truncation errors will be large if $h$ is too small. To compensate for this, the numerical algorithm would need to sum a large number of nodes $N$ to minimize truncation errors. Small values of $h$ will then tend to generate numerical inefficiencies. This has been the leading cause of numerical inefficiencies in previous implementations of the Ogata quadrature method. At the same time for larger values of $h$ the quadrature error grows as $\sim e^{-c/h}$, see Eq. (10). These observations indicate the need to find optimal values for $h$ and $N$ that keep the error term in Eq. (17) as small as possible. We found that such optimal values can be obtained by enforcing the largest contribution to the quadrature to be the first term in the truncated sum of Eq. (12) which can be achieved by maximizing the contribution of the first node, i.e.

[TABLE]

By solving numerically Eq. (21) for $h$ one finds the optimal value of $h$ for the quadrature method in Eq. (12). We will refer to this optimal value as $h_{\rm u}$.

It is now worth noting that $h_{\rm u}$ will tend to be a large value. This makes intuitive sense since minimizing truncation errors can be achieved by using a large spacing parameter. However, taking a large value of $h$ introduces quadrature errors which behave like $e^{-c/h}$ and tend to be large for $h=h_{\rm u}$. This issues can be mitigated by using the following scheme. We first use the condition in Eq. (21) to minimize truncation errors in Eq. (12). We then impose the condition that the final nodes of Eqs. (12) and (17) are placed at the same location by enforcing that

[TABLE]

This ensures that the quadrature in Eq. (17) has the same truncation errors as Eq. (12) with $h_{\rm u}$. The solution for $h$ in the above equality is given by

[TABLE]

when $h_{u}<\pi$. This value, labeled as $h_{t}$ is the optimal value for $h$ to be used in Eq. (17). Note that $h_{\rm t}$ is suppressed by a large factor of $\pi\xi_{nN}$ so that $h_{\rm t}\ll h_{\rm u}$. In Fig. 2 we plot the ratio $h_{\rm t}/h_{\rm u}$ as a function of $h_{\rm u}$ for $N=10,20$, and $40$. We find that in all cases $\frac{h_{\rm t}}{h_{\rm u}}\ll 1$ which avoids large errors in Eq. (21). We note that Eq. (22) only has a real solution when $h_{u}<\pi$. It turns out that when the variance of the input function is very large, the value of the parameter $h_{u}$ which was determined from Eq. (21) may be larger than $\pi$. To ensure that there is a real solution to Eq. (22) as well as to ensure that $h_{t}<h_{u}$ for all values of $N$, we must set an upper boundary on the parameter $h_{u}$, which we set to be $h_{max}=2$. When this occurs, the value of $h_{u}$ that we use is smaller than the optimal value. This issue could lead to truncation errors if the variance of the input function is sufficiently large. Beyond the range, the user needs to rely on the number of sampling points $N$ in order to decrease the errors stemming from the truncation errors.

In summary, for a given choice of number of integrand evaluations $N$, our procedure minimizes the error contribution in $h$ as well as truncation errors by applying the conditions Eqs. (21) and (22). The application of these conditions determines an optimal choice for $h$ in implementing the quadrature formula of Eq. (17). We shall refer to this procedure as “the optimized Ogata quadrature formula.” We will demonstrate below the efficiency of our formalism, first through the use of toy TMDs, and then through QCD based TMDs.

3 Benchmarking the Numerical Precision

In this section, we demonstrate the efficiency of the optimized Ogata quadrature method using toy TMDs for which the exact Fourier-Bessel transform is known. We will compare the numerical efficiency of the optimized Ogata quadrature against adaptive Gaussian quadrature, which is available in QUADPACK integration routine in Ref. [58]. It is important to emphasize that even though we mainly demonstrate the method for the integration involving Bessel function $J_{0}(x)$, we have checked that it works equally well for integration involving either $J_{1}(x)$ or $J_{2}(x)$, relevant for TMD studies in polarized scattering.

To assess the efficacy of our quadrature method we will map the error of the integration relative to the exact known result as a function of number of integrand calls. As discussed before we are interested in performing integrals of the form

[TABLE]

where the function $\widetilde{W}(b_{\perp})$ contains the $b_{\perp}$ space TMD physics. Therefore we will use a toy $\widetilde{W}(b_{\perp})$ which mimics the ${\bm{b}}_{\perp}$ space behavior of realistic unpolarized TMDs that has an exact analytic Fourier-Bessel transform. Specifically we choose the gamma distributions which are given in terms of the distribution’s mean, $\beta$, and variance, $\sigma$, as

[TABLE]

This function has an exponential $b_{\perp}$-dependence, and has been used in the literature for TMD studies [59]. Its exact Fourier-Bessel transform is given by

[TABLE]

where ${}_{2}\tilde{F}_{1}\left(a,b;c;d\right)$ is the regularized Gaussian hyper-geometric function.

The function $\widetilde{W}(b_{\perp},\beta,\sigma)$ has a single peak in $b_{\perp}$ space, which is given in terms of $\beta$ and $\sigma$ as

[TABLE]

We further introduce a parameter $Q$ to our toy TMD model, which is the inverse of the $b_{\perp}^{\rm peak}$, i.e.

[TABLE]

Such a $Q$-dependence mimics the hard scale $Q$ encountered in QCD based TMDs such as the photon virtuality in SIDIS reactions  [50, 60]. Notice that it is the quantity $q_{\perp}/Q$ that controls how oscillating the toy TMD is. The larger $q_{\perp}/Q$ is, the more oscillating the integrand is in $b_{\perp}$ space and the more numerically intensive the integration becomes.

For our demonstration, we take $Q=2$ (GeV) and $\sigma=1$ (GeV*-1*) similar to the usual JLab kinematics. We choose $q_{\perp}=0.2,2$, and $4$ (GeV), and plot the integrands on the left hand side of Fig. 3. As one can see clearly, the integrands do become more oscillating as $q_{\perp}/Q$ increases. To test the precision of our formalism, we take as an example in our optimized Ogata formula in Eq. (17).

[TABLE]

The relative percent error is defined as

[TABLE]

where the exact result $W^{\rm exact}(q_{\perp},\beta,\sigma)$ is given in Eq. (26). For this analysis, we increase the number of nodes at $q_{\perp}=0.2,2,4$ GeV until the best relative error for the numerical inversion is of the order of one percent. This requires $4,7$ and $10$ nodes at $q_{\perp}=0.2,2,4$ GeV, respectively. On the right panel of Fig. 3, we plot the relative percent error of the numerical integration as a function of the parameter $h$ for $q_{\perp}=0.2$, $2$, and $4$ GeV. One can see that in each case, there is an optimal value of the parameter $h$, which minimizes the measured error.

Intuitively having a small node spacing $h$ should result in a small error, since the error in $h$ is of the order $\mathcal{O}(e^{-c/h})$ in Eq. (17). However, the truncation errors $\mathcal{I}^{\prime}_{n\,N+1}$ in Eq. (17) will generate large errors in the numerical integration, due to the factors of $\,f\left(\frac{\pi}{h}\psi(h\xi_{nj})\right)\,J_{n}\left(\frac{\pi}{h}\psi(h\xi_{nj})\right)$, unless one increases $N$ significantly. Therefore, for a small and fixed number of nodes $N$, there is an optimized $h$ that minimizes the errors as argued in Sec. 2. On the right panel of Fig. 3, we indicate with stars the values of $h$ which are determined by the optimization conditions Eqs. (21) and (22). We find that within this range of kinematics, our optimization conditions indeed determine suitable values of $h$ for our quadrature method which is key to achieve high efficiency in the numerical integration.

We now compare the efficiency of the optimized Ogata quadrature against adaptive Gaussian quadrature. Note that the optimized Ogata quadrature first samples the integrand to determine the optimal value of $h$ using Eq. (21). The integration routine then samples the integrand at $N$ nodes to perform the sum. This results in a total number of integrand calls $N_{tot}$. To control $N_{tot}$ for adaptive Gaussian quadrature, we change the relative error tolerance of the integrator. Likewise, we measure the total number of function calls of adaptive Gaussian quadrature requested by the numerical routine. In Fig. 4 we plot the relative error as a function of $N_{tot}$ for small, intermediate and large values of $q_{\perp}/Q$. As is evident, the optimized Ogata quadrature method is more efficient than adaptive Gaussian quadrature, for all three probed regions of $q_{\perp}/Q$ with relative errors that go below $0.1\%$.

4 Application to TMDs

In this section we use the optimized Ogata quadrature to calculate the SIDIS $q_{\perp}$-differential cross sections in QCD TMD factorization framework. We then use these calculations to describe COMPASS charged hadron multiplicity data [14]. In addition, we use adaptive Gaussian quadrature and Vegas Monte Carlo algorithm for the same calculations to benchmark the performance.

Let’s first summarize the basic ingredients for the implementation of the unpolarized SIDIS structure function $W$ in Eq. (5) in the CSS TMD framework [9, 47]. In such context the TMD PDFs and TMD FFs can be expressed as

[TABLE]

where $\mu$ is the renormalization scale, $\zeta$ is the rapidity scale, $C_{q/j}^{pdf}$ and $C_{j/q}^{ff}$ are perturbatively calculable coefficient functions (see Ref. [47]), and $f_{j/p}\left(\hat{x},\mu_{b*}\right)$ and $d_{h/j}\left(\hat{z},\mu_{b*}\right)$ are the standard collinear PDFs and FFs, respectively. We will use the initial scale $Q_{0}^{2}=2.4$ GeV2. We follow the usual $b_{*}$-prescription [9] to avoid the Landau pole of $\alpha_{s}$, with

[TABLE]

The perturbative Sudakov factor $S_{\rm pert}$ is given by

[TABLE]

i.e., the evolution is done from the auxiliary scale $\mu_{b_{*}}=2e^{-\gamma_{E}}/b_{*}$ to the scale $\mu$. In the actual phenomenology, we set the rapidity scale $\zeta=Q^{2}$ and set the renormalization scale $\mu=Q$. We will implement the TMD evolution at next-to-leading-logarithmic (NLL) accuracy, and use the coefficient functions $C$ at next-to-leading order (NLO). All the relevant NLO coefficients and NLL anomalous dimensions can be found in Refs. [9, 47]. In addition, we use NLO expression for hard function $H(Q,\mu)$ in Eq. (5) from Ref. [47] reads

[TABLE]

and we set $\mu=Q$ so that the logarithmic terms vanish. Finally, we choose the parametrizations for the non-perturbative factors used in Refs.  [61, 62] which are given by

[TABLE]

with $g_{q}=0.106$ GeV2, $g_{2}=0.84$, and $g_{h}=0.042$ GeV2.

The expression for the $W$ term in Eq. (5) is given by

[TABLE]

where $\otimes$ is the convolution operator given in Eqs. (38) and (39) in [61].

Having established the QCD based TMD setups, let’s examine the behavior of the SIDIS cross section in $b_{\perp}$ space. In Fig. 5 the $b_{\perp}$ space integrand given in Eq. (5) for the SIDIS differential cross section is plotted, for three different values of $q_{\perp}/Q=0.1,~{}1.0$ and $2.0$, respectively. We take the lepton-proton center-of-mass energy square $S_{\ell p}=52.7$ (GeV2), $x_{\rm bj}=0.25$, $z=0.5$, and $Q^{2}=2.5$ (GeV2). These kinematics are within the coverage of the pion production in unpolarized lepton-proton SIDIS data at the HERMES experiment  [13]. Just like in the case of the toy TMDs in Sec. 3, the integrand becomes more oscillating as $q_{\perp}/Q$ increases. As a consequence, the numerical estimation of the Fourier-Bessel transform from $b_{\perp}$-space to $q_{\perp}$-space becomes increasingly more challenging for larger values of $q_{\perp}/Q$.

We next analyze the performance of our optimized Ogata quadrature to get the $q_{\perp}$-space cross sections against the fixed $h$ Ogata, adaptive Gaussian quadrature, and the Vegas Monte Carlo integration routines. For that we consider a related experimental observable, the hadron multiplicity which has been measured by both HERMES [13] and COMPASS [14] experiments. The COMPASS hadron multiplicity is defined as [14]

[TABLE]

where the numerator is the SIDIS cross section for the production of a hadron $h$ that we have been discussing so far, and the denominator is the inclusive DIS cross section. In Fig. 6, we plot the absolute value of the theoretical prediction for the hadron multiplicity, $|M^{h^{+}}(q_{\perp},x_{\rm bj},z,Q)|$, as a function of $q_{\perp}/Q$, using the above mentioned integration algorithms. For illustration, we choose the kinematics to be consistent with hadron multiplicity data from COMPASS experiment: $0.02<x_{\rm bj}<0.032$, $z=0.2$, and 1.7 GeV2 $<Q^{2}<$ 3 GeV2. The four panels correspond to optimized Ogata in the top left (labeled as ‘Opt. Ogata’), standard Ogata in the top right (labeled as ‘Fixed $h$ Ogata’), adaptive Gaussian quadrature in the bottom left (labeled as Quad), and Vegas Monte Carlo in the bottom right (labeled as Vegas), respectively. We also plot the COMPASS experimental data in Fig. 6 for comparison. 222In order to describe the data, the normalization issue with the COMPASS data must be resolved. We follow the work done in [63] to normalize the COMPASS multiplicities such that the data and theory are equal at the lowest values of the transverse momentum in each $z$ bin.

Note that at relatively large hadron transverse momentum $q_{\perp}/Q\gtrsim 2$, the theoretical calculations in TMD factorization formalism would become negative. There, one has to include the so-called $Y$-term [9], or switch/match onto the usual collinear factorization formalism [50, 64]. For very large values of $Q$, such an switching/matching from TMD factorization to simply a collinear factorization is straightforward and usually happens when $q_{\perp}$ becomes very large $q_{\perp}\sim Q$. It turns out that for small values of $Q$ (order of several GeV) which are most relevant for data at HERMES, COMPASS, and JLab, the matching is very tricky and does not occur as a sharp transition in the usual $W+Y$ prescription. Due to this, it is very important that in phenomenological applications the matching can be implemented without large numerical errors stemming from the Bessel transform using an efficient algorithm such as the one we are proposing. In the vicinity of $q_{\perp}\sim Q$ the calculation should be specifically precise in order to allow for the transition to happen. Moreover since the precise location of where the transition occurs can only be confronted phenomenologically, it is important that the ability to transform the $W$ term into $q_{\perp}$ space is as precise as possible to avoid any biases for the TMD extraction. See, e.g. Refs. [65, 66] for more details. It is because of this reason that our demonstration in Fig. 6 covers the broad region of $q_{\perp}$, from the small $q_{\perp}\ll Q$ to much larger $q_{\perp}\gtrsim Q$.

For each integration method in this Fig. 6, the number of nodes is increased until the relative error of the inversion associated with doubling the number of nodes is smaller than the total experimental uncertainty for all data points in the set 333In other applications of our optimized Ogata quadrature algorithm, the same method can be used to estimate the values for $N$ for a given integrand and desired precision.. This curve is plotted in orange, while in green and blue we plot the inversions with roughly half and double the number of nodes as the orange curve. For the optimized and fixed $h$ Ogata methods, we provide the total number of calls to the integrand $N_{tot}$. We note that the total number of calls to the integrand for the optimized Ogata contains two contributions. The algorithm first uses the gradient of the $h$-space function $h(h\xi_{n1})^{2n+1}\widetilde{W}(x_{\rm bj},z,h\xi_{n1}/q_{\perp},Q)$ to determine $h_{opt}$ in Eqs. (21) and (22). The numerical algorithm then samples the nodes to perform the quadrature sum. The first contribution tends to require more sampling points than the second contribution. However, the perturbative factors as well as the collinear distribution functions in this expression vary much more slowly in $h$-space than the non-perturbative factors. We find that this tends to be true in general for processes with small $Q$, where non-perturbative TMD structure is more sensitive. One can then sample only the expression $h(h\xi_{n1})^{2n+1}(h\xi_{n1}/q_{\perp}e^{-S_{NP}})$ to determine the value of $h_{opt}$. Since sampling this expression requires little computational power, we do not count this contribution to $N_{tot}$. We also note that in principle in an actual fit where a $\chi^{2}$ minimization is performed, one does not need to estimate the $h$ values for every step of $\chi^{2}$ evaluation. Instead these values can be updated every certain number of steps during the minimization. This will further enhance the efficiency of our method compared to the other approaches. In the adaptive quadrature and Vegas Monte Carlo methods, we indicate the average number of calls to the integrand $\left<N_{tot}\right>$.

As one can note, in the limit of large sampling, all the numerical integrators converge to the same result. However, the optimized Ogata quadrature converges to this result nearly an order of magnitude faster than the fixed $h$ method, nearly two orders of magnitude faster than adaptive Gaussian quadrature, and nearly three orders of magnitude faster than Vegas Monte Carlo integration. This result demonstrates that our optimized Ogata algorithm can improve significantly the numerical efficiency of the Fourier-Bessel integration encountered in the TMD analysis.

Finally in Fig. 7, we plot four multiplicity distributions at different values of $\langle z\rangle=0.2,0.3,0.4,0.6$, respectively for the bins $0.02<x_{\rm bj}<0.032$ and 1.7 GeV2 $<Q^{2}<$3 GeV2 using the optimized Ogata algorithm with the number of nodes $N=6$. $N=6$ was chosen as it gave reliable results over the entire region in Fig. 6. It is worthwhile to emphasize again that the theory predictions 444Note that the $q_{\perp}$-distribution in the TMD formalism gets narrower as $\langle z\rangle$ increases. This is expected from theoretical consideration, see e.g. Ref. [67]. become extremely efficient, thanks to the optimized Ogata quadrature. This gives us a great confidence that the optimized Ogata method would be ideal in the future for performing efficient numerical calculations and/or for the global analysis of TMDs. 555Note that we are not presenting a new fit here. Rather we just display, using fixed parameters from Refs. [61, 62], that this numerical method can be used to perform efficient numerical calculations for describing TMD data.

5 Conclusions

In this paper we have developed a high performance numerical algorithm for Hankel transforms for TMD factorization formalism from position $b_{\perp}$ space to transverse momentum $q_{\perp}$ space using the optimized Ogata quadrature method, which uses the zeros of Bessel functions as nodes. For a relatively small and fixed number $N$ of functional calls to the integrand, we derived conditions to find the optimal parameter $h$, which controls the node density. Such an optimized Ogata quadrature ensures a small number of calls while achieving a high accuracy at the same time, and thus becomes extremely efficient in TMD studies. We use both toy TMDs, and parametrizations of QCD based TMDs to demonstrate the efficiency of our integration algorithm. We found that the optimized Ogata quadrature performs nearly an order of magnitude faster than standard Ogata methods, nearly two orders of magnitude faster than adaptive Gaussian quadrature, and nearly three orders of magnitude faster than Vegas Monte Carlo integration for all regions of transverse momentum in semi-inclusive deep inelastic scattering. Our algorithm thus can have wide application in the future TMD computations and/or TMD global analysis. The code which illustrates the optimized Ogata quadrature is available for download in Python2, Python3, C++ with Boost dependency, C++ with GSL dependency, and a standalone Fortran 77 with an open source license at https://github.com/UCLA-TMD/Ogata. For information on installation and usage visit https://ucla-tmd.github.io/Ogata/.

Acknowledgements

Z.K. is supported by the National Science Foundation under Grant No. PHY-1720486 and No. PHY-1945471. A.P. is supported by the National Science Foundation under Grant No. PHY-1623454 and the DOE Contract No. DE-AC05-06OR23177, under which Jefferson Science Associates, LLC operates Jefferson Lab. N.S. was supported by the DOE contract DE-SC0018106. J.T. is supported by the National Science Foundation under Grant No. DGE-1650604. This work is also supported within the framework of the TMD Topical Collaboration.