DESY-24-057

Efficient computation of Fourier-Bessel transforms
for transverse-momentum dependent parton distributions
and other functions

Markus Diehl and Oskar Grocholski

Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

We present a method for the numerical computation of Fourier-Bessel transforms
on a finite or infinite interval. The function to be transformed needs to be
evaluated on a grid of points that is independent of the argument of the
Bessel function. We demonstrate the accuracy of the algorithm for a wide range
of functions, including those that appear in the context of transverse-
momentum dependent parton distributions in Quantum Chromodynamics.

###### Contents

  1. 1 Introduction
  2. 2 Developing the method
    1. 2.1 The original Levin method
    2. 2.2 Rescaling: general form
    3. 2.3 Rescaling with $r(z)=z$
      1. 2.3.1 Solution if $f_{1}$ or $f_{2}$ are powers of $z$
      2. 2.3.2 Behaviour of solutions at large $z$
  3. 3 Final form of the method
    1. 3.1 Behaviour of solutions at small and large $z$
    2. 3.2 Discretisation on Chebyshev grids
      1. 3.2.1 Solving the linear equation system
      2. 3.2.2 Examples of variable transforms
    3. 3.3 Summary of the algorithm
    4. 3.4 Implementation in BestLime
  4. 4 Numerical studies
    1. 4.1 Control parameters for solving the linear equation system
    2. 4.2 Benchmarks: grids and precision
    3. 4.3 Computation time
  5. 5 Numerical studies for TMD-like functions
    1. 5.1 Integrands
    2. 5.2 Choosing interpolation grids
      1. 5.2.1 Variable transformations
      2. 5.2.2 Number of subgrids
      3. 5.2.3 Transformation parameters and number of grid points
    3. 5.3 Estimating the integration error
    4. 5.4 Comparison with the method of Kang et al.
  6. 6 Summary
  7. A Useful relations for Bessel and Lommel functions
    1. A.1 Bessel functions
    2. A.2 Lommel functions
  8. B Integrands with a Gaussian decrease in $z$
  9. C Chebyshev interpolation

##  1 Introduction

The Fourier-Bessel transform — also known as Hankel transform — appears in
many branches of physics; in particular it arises from the multi-dimensional
Fourier transform of functions with rotational symmetry. It has the form

$\displaystyle I(q)$
$\displaystyle=\int_{z_{a}}^{z_{b}}dz\,J_{\nu}(qz)\,\tilde{f}(z)\,,$ (1.1)

where $J_{\nu}$ is a Bessel function of the first kind. In general the
integration runs from $0$ to $\infty$, but it may be reduced to a finite
interval by the support properties of $\tilde{f}(z)$.

In Quantum Chromodynamics, the Fourier-Bessel transform prominently appears in
the context of transverse-momentum dependent parton distributions (TMDs),
where cross sections are proportional to integrals of the form (1.1) with
$\nu=0,1,2,3$ and $\tilde{f}(z)$ equal to $z^{\nu+1}$ times the product of two
TMDs,111We note that in the TMD literature, the spatial distance is usually
called $b$ rather than $z$. see e.g. section 2.11 in [1]. The integration
limits are $0$ and $\infty$ in that case. The method proposed in [2] involves
similar integrals over TMDs from $0$ to some finite value, and from that
finite value to $\infty$. A more involved application is given by two-parton
TMDs, which appear in the description of double parton scattering [3]. In this
case, two nested Fourier transforms occur, the radial parts of which can be
converted into two nested Fourier-Bessel integrals, $\int dz_{1}\int
dz_{2}\,J_{\nu_{1}}(q_{1}z_{1})\,J_{\nu_{2}}(q_{2}\mskip
1.5muz_{2})\,\tilde{f}(z_{1},z_{2})$.

In the cases just mentioned, one typically needs to compute the transform
(1.1) for many settings of external parameters, namely longitudinal parton
momentum fractions and hard scales. For uncertainty estimates on TMDs fitted
to data (using the Hessian method or replicas) one must repeat the calculation
for many TMD sets. The numerical computation of the integrands itself can be
rather costly; at small transverse distances $z$ it involves Mellin
convolutions in the longitudinal momentum fractions, with lengthy kernels at
higher perturbative orders (see section 2.8 in [1] and references therein).

In the TMD context, the variable $q$ in (1.1) is a measured transverse
momentum, for instance of the produced lepton pair in the Drell-Yan process.
The process to be described is characterised by a hard scale $Q$, which can
range from a few $\operatorname{GeV}$ for semi-inclusive deep inelastic
scattering or non-resonant Drell-Yan production to about
$100\operatorname{GeV}$ for the production of a $W$, $Z$, or Higgs boson. For
the production of a pair of such bosons, the hard scale $Q$ is given by the
invariant mass of the pair and can be significantly larger. The TMD formalism
describes cross sections in the region $q\ll Q$, whilst a description in terms
of conventional collinear parton distributions is adequate for $q\sim Q$. A
smooth transition between the two regimes typically requires evaluating the
TMD expression for $q$ up to order $Q$, see section 4.7 in [1] for an overview
and [4] for a recent example.

In summary, the calculation of cross sections with TMDs may involve a large
number of Fourier-Bessel transforms of complicated functions for a possibly
wide range of $q$ values. In such a scenario, integration methods with nodes
that depend on $q$ come at a rather high computational cost since the
integrand needs to be computed many times. This is the case for the method of
Ogata [5], which is used for TMDs for instance in the codes artemide [6] or
ResBos2 [7], the latter employing the algorithm of Kang et al. [8].222The
respective codes and documentation can be found at
https://github.com/VladimirovAlexey/artemide-public for [6] and at
https://github.com/UCLA-TMD/Ogata and https://ucla-tmd.github.io/Ogata for
[8]. A method using fixed nodes in $z$ for a whole set of $q$ values is the
discrete Hankel transform, see [9, 10] for a mathematical description and [11]
for an application in hadron physics. Accessing high $q$ requires a large
number of nodes in this method.

In the present work, we adapt a method for oscillatory integrands due to Levin
[12, 13] to the specific case of Fourier-Bessel transforms. We show that high
accuracy can be achieved by combining this method with Chebyshev interpolation
(which was already done in [14, 15], albeit not for Bessel weighted
integrals). Specifically, we present a method to compute integrals of the form

$\displaystyle I(q)$
$\displaystyle=\int_{z_{a}}^{z_{b}}dz\,J_{\nu}(qz)\,\biggl{(}\frac{1+z}{z}\biggr{)}^{\nu}\,f(z)\,,$
$\displaystyle\nu\geq 0,\ \,q>0,\ \,z_{a}<z_{b}$ (1.2)

with the following specifications:

  * •

$z_{a}$ may be zero or positive, and $z_{b}$ may be finite or $\infty$.

We find that it can be useful for both accuracy and computation time to
evaluate the integral over $z\in[0,\infty]$ as the sum of integrals over two
or more subintervals.

  * •

$\nu$ may be integer or non-integer. Notice that for half integer $\nu$, the
integral (1.2) can be rewritten in terms of the spherical Bessel functions
$j_{n}(x)=\sqrt{\pi/(2x)}\,J_{n+1/2}(x)$.

  * •

$f(z)$ must be finite on the full integration interval, including the end
points $z_{a}$ and $z_{b}$. To obtain accurate integrals, $f(z)$ should should
be sufficiently smooth and not have fast oscillations. We will see that a non-
analytic behaviour like $f(z)\sim z^{1-\delta}$ at $z=0$ is amenable to the
method if $\delta$ is not too large, although the accuracy tends to degrade
with increasing $\delta$.

  * •

The function $f(z)$ is evaluated on a discrete set of $n$ points from $z_{a}$
to $z_{b}$ (specified in section 3.2). We find that with a suitable choice of
discretisation grid, good integration accuracy can be obtained for a rather
wide range of $q$.

We have implemented this method in a C++ library named BestLime,333Bessel
transformation with Levin’s integration method. which can be downloaded from
[16] and is briefly discussed in sections 3.4 and 4.3 of the present paper.
Although our primary physics motivation is the application to TMDs, we think
that the method may be of interest to other domains. We show in section 4.2
that it can handle a large variety of functions $f(z)$.

Our paper is organised as follows. In section 2 we develop our method starting
from the original one in [13], and in section 3 we present its final form. We
discuss general numerical aspects in section 4 and integrals appearing in TMD
cross sections in section 5. Our findings are summarised in section 6. Various
mathematical details are given in appendix A, B, and C.

##  2 Developing the method

###  2.1 The original Levin method

To begin with, let us present the original method of Levin [13] for the case
of Fourier-Bessel integrals. We use a slightly different notation than the
cited work. The integral to be evaluated has the form

$\displaystyle I(q)$
$\displaystyle=\int_{z_{a}}^{z_{b}}dz\,\Bigl{[}\,\widetilde{\omega}_{1}(z,q)\,\tilde{f}_{1}(z)+\widetilde{\omega}_{2}(z,q)\,\tilde{f}_{2}(z)\,\Bigr{]}$
(2.1)

with weights

$\displaystyle\widetilde{\omega}_{1}(z,q)$ $\displaystyle=J_{\nu}(qz)\,,$
$\displaystyle\widetilde{\omega}_{2}$ $\displaystyle=J_{\nu+1}(qz)\,,$ (2.2)

where it is assumed that the functions $\tilde{f}_{1}(z)$ and
$\tilde{f}_{2}(z)$ are _not_ strongly oscillating. The Fourier-Bessel
transform (1.1) is obtained by setting $\tilde{f}_{1}=\tilde{f}$ and
$\tilde{f}_{2}=0$. Setting $\tilde{f}_{1}=0$ and $\tilde{f}_{2}=\tilde{f}$
instead, one obtains the analogue of (1.1) for $J_{\nu+1}$ instead of
$J_{\nu}$. This requires only a minor additional computational effort; a bonus
that is preserved by our adaptation of Levin’s method.

To proceed, we introduce the vector notation

$\displaystyle\boldsymbol{\widetilde{\omega}}$
$\displaystyle=\bigl{(}\widetilde{\omega}_{1}\,,\widetilde{\omega}_{2}\bigr{)}^{T}\,,$
$\displaystyle\boldsymbol{\tilde{f}}$
$\displaystyle=\bigl{(}\tilde{f}_{1}\,,\tilde{f}_{2}\bigr{)}^{T}\,,$ (2.3)

where a superscript $T$ denotes transposition. Using the relations (A.1), one
readily derives

$\displaystyle\frac{d}{dz}\,\boldsymbol{\widetilde{\omega}}(z,q)$
$\displaystyle=\widetilde{A}^{\,T}(z,q)\cdot\boldsymbol{\widetilde{\omega}}(z,q)$
(2.4)

with the matrix

$\displaystyle\widetilde{A}(z,q)$ $\displaystyle=\begin{pmatrix}\nu/z\,,&q\\\
-q\,,&-(\nu+1)/z\end{pmatrix}\,.$ (2.5)

If one finds a solution for the system

$\displaystyle\boldsymbol{\tilde{f}}(z)$
$\displaystyle=\biggl{[}\frac{d}{dz}+\widetilde{A}(z,q)\biggr{]}\,\boldsymbol{\tilde{h}}(z,q)$
(2.6)

of ordinary differential equations (ODEs), one can evaluate the integral as

$\displaystyle I(q)$
$\displaystyle=\int_{z_{a}}^{z_{b}}dz\,\boldsymbol{\widetilde{\omega}}^{T}(z,q)\cdot\boldsymbol{\tilde{f}}(z)=\Bigl{[}\mskip
1.5mu\boldsymbol{\widetilde{\omega}}^{T}(z,q)\cdot\boldsymbol{\tilde{h}}(z,q)\mskip
1.5mu\Bigr{]}_{z_{a}}^{z_{b}}\;\;,$ (2.7)

where we use the common notation

$\displaystyle\bigl{[}F(z,q)\bigr{]}_{z_{a}}^{z_{b}}$
$\displaystyle=F(z_{b},q)-F(z_{a},q)\,.$ (2.8)

This is readily seen by inserting (2.6) in the first expression of (2.7),
integrating by parts, and then using (2.4).

As argued in [12, 13], the system (2.6) has one particular solution with
functions $\tilde{h}_{i}(z,q)$ that are not rapidly oscillating. One can find
this solution by collocation, which consists in approximating
$\tilde{h}_{1}(z,q)$ and $\tilde{h}_{2}(z,q)$ by a linear combination of $n$
basis functions $b_{i}(z)$ and by determining the coefficients of these
functions by evaluating the system (2.6) at a set of $n$ suitably chosen
collocation points $z_{i}$. This yields a coupled system of linear equations
that can be solved by standard linear algebra methods (see section 3.2.1 for
details). For the method to give accurate results, it is essential that the
approximation of $\tilde{h}_{1}$ and $\tilde{h}_{2}$ is sufficiently precise,
in particular at the end points $z_{a}$ and $z_{b}$, where these functions are
evaluated in the final result (2.7).

As it stands, this method is not suitable for a lower integration boundary
$z_{a}=0$, where the diagonal elements of $\widetilde{A}$ diverge. We
circumvent this problem by rescaling, as described in the next subsection.

###  2.2 Rescaling: general form

The representation (2.7) of $I(q)$ remains valid if one replaces
$\boldsymbol{\tilde{f}}$, $\boldsymbol{\tilde{h}}$ and
$\boldsymbol{\widetilde{\omega}}$ by the rescaled functions

$\displaystyle\boldsymbol{f}(z)$
$\displaystyle=r^{\nu}(z)\,\boldsymbol{\tilde{f}}(z)\,,$
$\displaystyle\boldsymbol{h}(z)$
$\displaystyle=r^{\nu}(z)\,\boldsymbol{\tilde{h}}(z)\,,$
$\displaystyle\boldsymbol{\omega}(z,q)$
$\displaystyle=r^{-\nu}(z)\,\boldsymbol{\widetilde{\omega}}(z,q)\,,$ (2.9)

where $r(z)$ will be specified below. Using the explicit form (2.5) of
$\widetilde{A}$, one can rewrite the system of ODEs (2.6) as

$\displaystyle f_{1}$ $\displaystyle=\frac{d}{dz}\,h_{1}+qh_{2}+\nu\mskip
1.5mu\bigg{[}\mskip 1.5mu\frac{1}{z}-\frac{r^{\prime}}{r}\mskip
1.5mu\biggr{]}\,h_{1}\,,$ $\displaystyle f_{2}$
$\displaystyle=\frac{d}{dz}\,h_{2}-qh_{1}-\biggl{[}\mskip
1.5mu\frac{\nu+1}{z}+\nu\,\frac{r^{\prime}}{r}\mskip 1.5mu\biggr{]}\,h_{2}\,,$
(2.10)

where we have omitted the arguments $z$ or $(z,q)$ of the functions and
abbreviated $r^{\prime}=dr/dz$. The singularity at $z=0$ in the first equation
of (2.2) is removed if

$\displaystyle r(z)$ $\displaystyle\sim z$ $\displaystyle\text{ for }z\to
0\,,$ (2.11)

and the singularity in the second equation is avoided by setting

$\displaystyle h_{2}(z,q)$ $\displaystyle=r(z)\,h_{3}(z,q)\,.$ (2.12)

The final form of the ODE system then reads

$\displaystyle f_{1}$ $\displaystyle=\frac{d}{dz}\,h_{1}+q\mskip 1.5mur\mskip
1.5muh_{3}+\nu\mskip 1.5mu\bigg{[}\mskip
1.5mu\frac{1}{z}-\frac{r^{\prime}}{r}\mskip 1.5mu\biggr{]}\,h_{1}\,,$
$\displaystyle f_{2}$
$\displaystyle=r\,\frac{d}{dz}\,h_{3}-qh_{1}-\biggl{[}\mskip
1.5mu(\nu+1)\,\frac{r}{z}+(\nu-1)\,r^{\prime}\mskip 1.5mu\biggr{]}\,h_{3}\,,$
(2.13)

and the integral $I(q)$ is given by

$\displaystyle I(q)$
$\displaystyle=\int_{z_{a}}^{z_{b}}dz\;\,r^{-\nu}(z)\,\Bigl{[}J_{\nu}(qz)\,f_{1}(z)+J_{\nu+1}(qz)\,f_{2}(z)\Bigr{]}$
$\displaystyle=\biggl{[}\mskip
1.5muJ_{\nu}(qz)\,r^{-\nu}(z)\,h_{1}(z,q)+J_{\nu+1}(qz)\,r^{-\nu+1}(z)\,h_{3}(z,q)\mskip
1.5mu\biggr{]}_{z_{a}}^{z_{b}}\,.$ (2.14)

We will show that if $f_{1}(z)$ and $f_{2}(z)$ are sufficiently well-behaved
at $z\to 0$, the solutions $h_{1}(z,q)$ and $h_{3}(z,q)$ of (2.2) are well
behaved at $z\to 0$ as well. This enables one to solve the system by
collocation with basis functions that have a Taylor expansion in $z$ around
$z=0$. We anticipate that our method will use basis functions with this
property, postponing details to section 3.2.

###  2.3 Rescaling with $r(z)=z$

The main purpose of this subsection is to explicitly solve the system (2.2) of
ODEs for particular choices of $f_{1}$ and $f_{2}$. This is an intermediate
step, and a reader mainly interested in final version of our method may skip
forward to section 3.

In the present subsection, we set

$\displaystyle r(z)$ $\displaystyle=z\,,$ (2.15)

such that the system of ODE takes the rather simple form

$\displaystyle f_{1}$ $\displaystyle=\frac{d}{dz}\,h_{1}+qz\mskip
1.5muh_{3}\,,$ (2.16) $\displaystyle f_{2}$ $\displaystyle=z\mskip
1.5mu\frac{d}{dz}\,h_{3}-qh_{1}-2\nu\mskip 1.5muh_{3}\,.$ (2.17)

One can eliminate $h_{1}$ by taking the derivative of (2.17) and then
inserting $dh_{1}/dz$ from (2.16). The result is

$\displaystyle z^{2}\mskip 1.5mu\frac{d^{2}}{dz^{2}}\,(z^{-\nu}\mskip
1.5muh_{3})+z\mskip 1.5mu\frac{d}{dz}\,(z^{-\nu}\mskip
1.5muh_{3})+\bigl{(}q^{2}z^{2}-\nu^{2}\bigr{)}\,(z^{-\nu}\mskip 1.5muh_{3})$
$\displaystyle=z^{1-\nu}\mskip 1.5mu\bigl{(}qf_{1}+f^{\mskip
1.5mu\prime}_{2}\mskip 1.5mu\bigr{)}$ (2.18)

with $f^{\mskip 1.5mu\prime}_{2}=d\mskip 0.1muf_{2}/dz$.

If one sets $f_{1}=f_{2}=0$, this is the Bessel equation for the function
$z^{-\nu}h_{3}$. The _homogeneous_ solution of equations (2.16) and (2.17) is
therefore given by

$\displaystyle\begin{pmatrix}h_{1}\\\ h_{3}\end{pmatrix}\Bigg{|}_{\text{hom}}$
$\displaystyle=c_{J}\mskip
1.5mu\begin{pmatrix}{}-(qz)^{\nu+1}\,J_{\nu+1}(qz)\\\
q\,(qz)^{\nu}\,J_{\nu}(qz)\end{pmatrix}+c_{Y}\mskip
1.5mu\begin{pmatrix}{}-(qz)^{\nu+1}\,Y_{\nu+1}(qz)\\\
q\,(qz)^{\nu}\,Y_{\nu}(qz)\end{pmatrix}$ (2.19)

with arbitrary coefficients $c_{J}$ and $c_{Y}$, where $Y_{\nu}$ denotes the
Bessel functions of the second kind. We obtained the solution for $h_{1}$ by
inserting the one for $h_{3}$ into (2.17) and using the relations (A.1) and
(A.2).

At small $z$, the Bessel functions $J_{\nu}$ and $Y_{\nu}$ can be represented
as sums of different power series in $z^{2}$, some of which are multiplied by
functions of $z$. To represent this in a compact form, we write

$\displaystyle\sigma(x)$ $\displaystyle=\sum_{k=0}^{\infty}a_{k}\mskip
1.5mux^{k}$ (2.20)

to denote a power series in the variable $x$. It is understood that the series
may terminate, so that one has a polynomial in $x$. We use the notation (2.20)
in a similar way as the familiar symbol $\mathcal{O}(x)$, i.e., different
occurrences of $\sigma(x)$ (even within a single formula) may have different
coefficients $a_{k}$, and these coefficients may depend on variables other
than $x$ (such as $q$ and $\nu$ in our context).

With this notation, the Bessel functions can be written in the form (A.1),
which implies a small-$z$ behaviour

$\displaystyle\begin{pmatrix}h_{1}\\\ h_{3}\end{pmatrix}\Bigg{|}_{\text{hom}}$
$\displaystyle=\begin{pmatrix}\sigma(z^{2})+z^{2\nu+2}\,\sigma(z^{2})+z^{2\nu+2}\mskip
1.5mu\ln(z)\,\sigma(z^{2})\\\
\sigma(z^{2})+z^{2\nu}\,\sigma(z^{2})+z^{2\nu}\mskip
1.5mu\ln(z)\,\sigma(z^{2})\end{pmatrix}$ (2.21)

of the homogeneous solution (2.19), where we assumed $\nu\geq 0$. The terms
with $\ln(z)$ only appear for integer values of $\nu$. We see that the
homogeneous solutions are non-analytic at $z=0$, except for half-integer
$\nu$. Recalling that the basis functions we will use for collocation have a
Taylor expansion around $z=0$, we require

$\displaystyle\nu$ $\displaystyle\geq 1\,,$ (2.22)

when using (2.2), which ensures that the functions (2.21) and their first
derivatives in $z$ are finite at that point. Integrals with $J_{\nu}$ for
$\nu<1$ can be evaluated using integration by parts, as specified in equation
(3) below.

####  2.3.1 Solution if $f_{1}$ or $f_{2}$ are powers of $z$

Let us first consider the case where $f_{2}=0$ and

$\displaystyle f_{1}(z)$ $\displaystyle=c_{1}\mskip 1.5muz^{\mu+\nu}\,,$
(2.23)

which corresponds to the integral of $c_{1}\mskip 1.5muz^{\mu}J_{\nu}(qz)$ in
(2.2). Notice that a finite value of $f_{1}(z)$ at $z=0$ requires $\mu+\nu\geq
0$.

The r.h.s. of (2.18) now equals $c_{1}qz^{1+\mu}$, and we recognise the Lommel
differential equation (A.6) for the function $z^{-\nu}h_{3}$ (up to a constant
on the r.h.s.). The equation is solved by $c_{1}\mskip 1.5muq^{-\mu}$ times
the Lommel functions of the first or second kind, $s_{\mu,\nu}(qz)$ or
$S_{\mu,\nu}(qz)$, which differ by a linear combination $c_{J}\mskip
1.5muJ_{\nu}(qz)+c_{Y}\mskip 1.5muY_{\nu}(qz)$ of the homogeneous solutions,
see (A.2).

It turns out that it is the Lommel functions of the second kind that are free
of fast oscillations in $z$. This is illustrated for selected values of $\mu$
and $\nu$ in figure 1. In the plots, we have scaled the functions as

$\displaystyle x^{\nu}\,(1+x)^{1-\mu-\nu}\,S_{\mu,\nu}(x)$ (2.24)

for clarity. According to (A.15) the scaled functions tend to $1$ for
$x\to\infty$, and according to (A.20) they have a finite limit for $x\to 0$ if
$\mu+\nu\geq 0$ and $\nu\geq 1$.

(a)

(b)

(c)

(d)

Figure 1:  Lommel functions of the second kind, scaled as specified in (2.24).
The functions for $\mu+\nu=1$ are not shown, since according to (A.13) they
are identical to $1$.

The corresponding solution for $h_{1}$ is readily obtained by inserting the
one for $h_{3}$ into (2.17). Using the relations (A.2), we obtain

$\displaystyle\begin{pmatrix}h_{1}\\\
h_{3}\end{pmatrix}\Bigg{|}_{\text{$f_{1}$}}$
$\displaystyle=\tilde{c}_{1}\begin{pmatrix}(\mu-\nu-1)\,(qz)^{\nu+1}\,S_{\mu-1,\nu+1}(qz)\\\
q\,(qz)^{\nu}\,S_{\mu,\nu}(qz)\end{pmatrix}$ (2.25)

with $\tilde{c}_{1}=c_{1}\mskip 1.5muq^{-\mu-\nu-1}$.

We now turn to the case where $f_{1}=0$ and

$\displaystyle f_{2}(z)$ $\displaystyle=c_{2}\mskip 1.5muz^{\mu+\nu}\,,$
(2.26)

which corresponds to the integral of $c_{2}\mskip 1.5muz^{\mu}J_{\nu+1}(qz)$
in (2.2). The r.h.s. of (2.18) then equals $(\mu+\nu)\,c_{2}\mskip
1.5muz^{\mu}$, and the equation is again solved by Lommel functions, but with
the first index decreased by $1$ and with a different overall factor:

$\displaystyle\begin{pmatrix}h_{1}\\\
h_{3}\end{pmatrix}\Bigg{|}_{\text{$f_{2}$}}$
$\displaystyle=(\mu+\nu)\,\tilde{c}_{2}\begin{pmatrix}(\mu-\nu-2)\,(qz)^{\nu+1}\,S_{\mu-2,\nu+1}(qz)\\\
q\,(qz)^{\nu}\,S_{\mu-1,\nu}(qz)\end{pmatrix}-\tilde{c}_{2}\begin{pmatrix}(qz)^{\mu+\nu}\\\
0\end{pmatrix}$ (2.27)

where $\tilde{c}_{2}=c_{2}\,q^{-\mu-\nu-1}$. The solution for $\mu+\nu=0$ is
trivial to verify.

By virtue of (A.20), the functions in (2.25) and (2.27) have series expansions

$\displaystyle\begin{pmatrix}h_{1}\\\
h_{3}\end{pmatrix}\Bigg{|}_{\text{$f_{1}$}}$
$\displaystyle=\begin{pmatrix}\sigma(z^{2})+z^{\mu+\nu+1}\,\sigma(z^{2})+z^{2\nu+2}\,\sigma(z^{2})+z^{2\nu+2}\mskip
1.5mu\ln(z)\,\sigma(z^{2})\\\
\sigma(z^{2})+z^{\mu+\nu+1}\,\sigma(z^{2})+z^{2\nu}\,\sigma(z^{2})+z^{2\nu}\mskip
1.5mu\ln(z)\,\sigma(z^{2})\end{pmatrix}$ (2.28) and
$\displaystyle\begin{pmatrix}h_{1}\\\
h_{3}\end{pmatrix}\Bigg{|}_{\text{$f_{2}$}}$
$\displaystyle=\begin{pmatrix}\sigma(z^{2})+z^{\mu+\nu}\,\sigma(z^{2})+z^{2\nu+2}\,\sigma(z^{2})+z^{2\nu+2}\mskip
1.5mu\ln(z)\,\sigma(z^{2})\\\
\sigma(z^{2})+z^{\mu+\nu}\,\sigma(z^{2})+z^{2\nu}\,\sigma(z^{2})+z^{2\nu}\mskip
1.5mu\ln(z)\,\sigma(z^{2})\end{pmatrix}\,,$ (2.29)

respectively. In addition to terms already present in the homogeneous solution
(2.21), we thus have a further power series in each function $h_{1}$ and
$h_{3}$, with a prefactor $z^{\mu+\nu+1}$ for $f_{1}$ and a prefactor
$z^{\mu+\nu}$ for $f_{2}$.

This implies that the power $\mu+\nu$ of $z$ in (2.23) or (2.26) must be large
enough if $h_{1}$ and $h_{3}$ are to be sufficiently well approximated by
basis functions that have a Taylor expansion around $z=0$. In particular, a
finite first derivative of $h_{1}$ and $h_{3}$ at $z=0$ requires $\mu+\nu\geq
0$ for $f_{1}$ and $\mu+\nu\geq 1$ for $f_{2}$.

####  2.3.2 Behaviour of solutions at large $z$

We now turn to the behaviour of the functions $h_{1}$ and $h_{3}$ in the limit
$z\to\infty$. If $f_{1}(z)$ or $f_{2}(z)$ is proportional to $z^{\mu+\nu}$, we
can use the the solutions (2.25) or (2.27) derived in the previous subsection.
According to (A.15), they behave like

$\displaystyle\begin{pmatrix}h_{1}\\\
h_{3}\end{pmatrix}\Bigg{|}_{\text{$f_{1}$}}$
$\displaystyle=\begin{pmatrix}z^{\mu+\nu-1}\;\sigma(z^{-2})\\\
z^{\mu+\nu-1}\;\sigma(z^{-2})\end{pmatrix}\,,$
$\displaystyle\begin{pmatrix}h_{1}\\\
h_{3}\end{pmatrix}\Bigg{|}_{\text{$f_{2}$}}$
$\displaystyle=\begin{pmatrix}z^{\mu+\nu}\;\sigma(z^{-2})\\\
z^{\mu+\nu-2}\;\sigma(z^{-2})\end{pmatrix}$ (2.30)

for $z\gg 1$. We postpone the discussion of this result to section 3.1.

Let us now consider functions with an exponential falloff in $z$, namely

$\displaystyle f_{1}(z)$ $\displaystyle=c_{1}\mskip
1.5muz^{\mu+\nu}\,\exp(-\kappa z)\,,$ $\displaystyle f_{2}(z)$
$\displaystyle=c_{2}\mskip 1.5muz^{\mu+\nu}\,\exp(-\kappa z)\,.$ (2.31)

with $\kappa>0$. To solve the system of ODEs in (2.16) and (2.17), we make the
ansatz

$\displaystyle h_{i}(z,q)$ $\displaystyle=\hat{h}_{i}(z,q)\,\exp(-\kappa z)$
$\displaystyle(i=1,3).$ (2.32)

This leads to the new system

$\displaystyle\begin{pmatrix}c_{1}\mskip 1.5muz^{\mu+\nu}\\\ c_{2}\mskip
1.5muz^{\mu+\nu}\end{pmatrix}$ $\displaystyle=-\begin{pmatrix}\kappa&-q\\\
q&\kappa\end{pmatrix}\cdot\begin{pmatrix}\hat{h}_{1}\\\
z\hat{h}_{3}\end{pmatrix}+\begin{pmatrix}d\mskip 1.5mu\hat{h}_{1}/dz\\\
z\mskip 1.5mud\mskip 1.5mu\hat{h}_{3}/dz-2\nu\mskip
1.5mu\hat{h}_{3}\end{pmatrix}\,,$ (2.33)

which can be solved with the ansatz

$\displaystyle\hat{h}_{1}(z,q)$
$\displaystyle=z^{\mu+\nu}\,\sum_{k=0}^{\infty}b_{1,k}(q)\mskip
1.5muz^{-k}\,,$ $\displaystyle\hat{h}_{3}(z,q)$
$\displaystyle=z^{\mu+\nu-1}\,\sum_{k=0}^{\infty}b_{3,k}(q)\mskip
1.5muz^{-k}\,.$ (2.34)

For every single $k$ in the sums, the first term on the r.h.s. of (2.33) is
leading at large $z$, whereas the second term on the r.h.s. is down by one
power of $z$. For the leading coefficients, we readily get

$\displaystyle\begin{pmatrix}b_{1,0}\\\ b_{3,0}\end{pmatrix}$
$\displaystyle=-\frac{1}{\kappa^{2}+q^{2}}\begin{pmatrix}\kappa&q\\\
-q&\kappa\end{pmatrix}\cdot\begin{pmatrix}c_{1}\\\ c_{3}\end{pmatrix}\,,$
(2.35)

whereas the subleading coefficients can be obtained recursively as

$\displaystyle\begin{pmatrix}b_{1,k+1}\\\ b_{3,k+1}\end{pmatrix}$
$\displaystyle=\frac{1}{\kappa^{2}+q^{2}}\begin{pmatrix}\kappa&q\\\
-q&\kappa\end{pmatrix}\cdot\begin{pmatrix}(\mu+\nu-k)\,b_{1,k}\\\
(\mu-\nu-1-k)\,b_{3,k}\end{pmatrix}\,,$ (2.36)

starting with $k=0$.

The upshot of this calculation is that if the functions $f_{1}$ and $f_{2}$
are of the form “power times exponential” at large $z$, then the solutions at
large $z$ have the same form, with the same damping factor $\kappa$ in the
exponential.

In appendix B we show that a similar statement holds if $f_{1}$ and $f_{2}$
have a Gaussian decrease in $z$ (possibly modified by a power and an
exponential in $z$).

##  3 Final form of the method

In the previous section, we have shown that the rescaling function $r(z)=z$
leads to well behaved solutions $h_{1}$ and $h_{3}$ of the differential
equations (2.2) in the limit $z=0$. This choice is however not optimal for the
opposite limit $z\to\infty$, where the requirement that the functions $f_{1}$
and $f_{2}$ in (2.2) must remain finite at the end points of the integration
interval would reduce the class of integrals that can be evaluated. Our final
choice of rescaling functions is therefore

$\displaystyle r(z)$ $\displaystyle=\frac{z}{1+z}\,,$ (3.1)

which interpolates smoothly between the choice $r(z)=z$ at small $z$ and no
rescaling as $z\to\infty$. A more general version would be $r(z)=z/(1+\alpha
z)$ with $\alpha>0$, but we shall not explore this option here.

We can then compute the integral

$\displaystyle I(q)$
$\displaystyle=\int_{z_{a}}^{z_{b}}dz\;\biggl{(}\frac{1+z}{z}\biggr{)}^{\nu}\,\biggl{[}J_{\nu}(qz)\,f_{1}(z)+J_{\nu+1}(qz)\,f_{2}(z)\biggr{]}$
(3.2)

as

$\displaystyle I(q)$
$\displaystyle=\biggl{[}J_{\nu}(qz)\,\biggl{(}\frac{1+z}{z}\biggr{)}^{\nu}\,h_{1}(z,q)+J_{\nu+1}(qz)\,\biggl{(}\frac{1+z}{z}\biggr{)}^{\nu-1}\,h_{3}(z,q)\biggr{]}_{z_{a}}^{z_{b}}\;,$
(3.3)

where $h_{1}$ and $h_{3}$ solve the system

$\displaystyle f_{1}(z)$
$\displaystyle=\frac{d}{dz}\,h_{1}(z,q)+q\,\frac{z}{1+z}\,h_{3}(z,q)+\frac{\nu}{1+z}\,h_{1}(z,q)\,,$
$\displaystyle f_{2}(z)$
$\displaystyle=\frac{z}{1+z}\,\frac{d}{dz}\,h_{3}(z,q)-qh_{1}(z,q)-\biggl{[}\mskip
1.5mu\frac{\nu-1}{(1+z)^{2}}+\frac{\nu+1}{1+z}\mskip
1.5mu\biggr{]}\,h_{3}(z,q)$ (3.4)

of ODEs. For convenience, we will call (3) the _Levin equations_ in the
following.

For the reasons given in section 2.3, the relations (3.3) and (3) are used for
$\nu\geq 1$. To evaluate integrals with Bessel functions of order below $1$,
we instead use the relation (A.2) and subsequent integration by parts, which
gives

$\displaystyle\int_{z_{a}}^{z_{b}}dz\;J_{\nu-1}(qz)\,\biggl{(}\frac{1+z}{z}\biggr{)}^{\nu-1}\,f_{0}(z)$
$\displaystyle\qquad=\frac{1}{q}\int_{z_{a}}^{z_{b}}dz\;\frac{d}{dz}\mskip
1.5mu\biggl{[}z^{\nu}\mskip
1.5muJ_{\nu}(qz)\biggr{]}\;\frac{1}{z^{\nu}}\,\biggl{(}\frac{1+z}{z}\biggr{)}^{\nu-1}\,f_{0}(z)$
$\displaystyle\qquad=\frac{1}{q}\biggl{[}J_{\nu}(qz)\,\biggl{(}\frac{1+z}{z}\biggr{)}^{\nu-1}\,f_{0}(z)\biggr{]}_{z_{a}}^{z_{b}}-\frac{1}{q}\int_{z_{a}}^{z_{b}}dz\;J_{\nu}(qz)\,\biggl{(}\frac{1+z}{z}\biggr{)}^{\nu}\,f_{1}(z)$
(3.5)

with

$\displaystyle f_{1}(z)$ $\displaystyle=\frac{z}{1+z}\mskip
1.5mu\frac{d}{dz}\mskip
1.5muf_{0}(z)-\biggl{[}\frac{\nu-1}{(1+z)^{2}}+\frac{\nu}{1+z}\biggr{]}\mskip
1.5muf_{0}(z)\,.$ (3.6)

###  3.1 Behaviour of solutions at small and large $z$

From equation (2.2) it is clear that one obtains the same integral with
different rescaling factors by a suitable transformation of the functions
$f_{i}$ and $h_{i}$ in the system of ODEs. For the choices made in the present
and in the previous section, we have

$\displaystyle f_{i}(z)\mskip 1.5mu\big{|}_{r=z/(1+z)}$
$\displaystyle=(1+z)^{-\nu}\;f_{i}(z)\mskip 1.5mu\big{|}_{r=z}$
$\displaystyle(i=1,2),$ $\displaystyle h_{1}(z,q)\mskip
1.5mu\big{|}_{r=z/(1+z)}$ $\displaystyle=(1+z)^{-\nu}\;h_{1}(z,q)\mskip
1.5mu\big{|}_{r=z}\,,$ $\displaystyle h_{3}(z,q)\mskip
1.5mu\big{|}_{r=z/(1+z)}$ $\displaystyle=(1+z)^{-\nu+1}\;h_{3}(z,q)\mskip
1.5mu\big{|}_{r=z}\,.$ (3.7)

We thus readily obtain the small and large $z$ behaviour of $h_{i}$ with our
final rescaling choice from our results in section 2.3.

For the limit $z\to 0$, we deduce from (2.28) and (2.29) that

$\displaystyle f_{1}(z)$ $\displaystyle\sim z^{\mu+\nu}\,\sigma(z)$
$\displaystyle\Rightarrow\;h_{3}(z,q)$
$\displaystyle\sim\sigma(z)+z^{\mu+\nu+1}\,\sigma(z)+z^{2\nu}\,\sigma(z)+z^{2\nu}\mskip
1.5mu\ln(z)\,\sigma(z)$ (3.8) and $\displaystyle f_{2}(z)$ $\displaystyle\sim
z^{\mu+\nu}\,\sigma(z)$ $\displaystyle\Rightarrow\;h_{3}(z,q)$
$\displaystyle\sim\sigma(z)+z^{\mu+\nu}\,\sigma(z)+z^{2\nu}\,\sigma(z)+z^{2\nu}\mskip
1.5mu\ln(z)\,\sigma(z)\,,$ (3.9)

where we recall that $\sigma(z)$ denotes a power series in $z$ in the generic
sense explained after equation (2.20). The behaviour of $h_{1}(z,q)$ is as the
one of $h_{3}(z,q)$ in both (3.1) and (3.1), except that the power $z^{2\nu}$
is replaced with $z^{2\nu+2}$. It is understood that $\mu+\nu\geq 0$, which
ensures that $f_{1}$ and $f_{2}$ remain finite at $z=0$ and that the integral
(3.2) with $z_{a}=0$ actually exists.

We recall from our discussion in section 2.3.1 that small non-integer values
of $\mu+\nu$ lead to a non-analytic behaviour of $h_{1}$ and $h_{3}$ at $z=0$
that is difficult to approximate by functions having a Taylor expansion around
that point. We therefore expect our method to become less accurate in this
case. The comparison between (3.1) and (3.1) reveals that for an integrand
behaving like $z^{\mu}J_{\rho}(qz)$ around $z=0$, one should use (3.2) with
$\nu=\rho$ and $f_{2}=0$ if $1<\mu+\rho<2$, rather than with $\nu=\rho-1$ and
$f_{1}=0$. The same holds if $0\leq\mu+\rho<1$. For $\mu+\rho=1$ or
$\mu+\rho\geq 2$, both options are possible.

For the large-$z$ limit, we can discuss several cases. From (2.30) we deduce
that for a power-law behaviour of $f_{i}$ we have

$\displaystyle f_{1}(z)$ $\displaystyle\sim
z^{\mu}\,\sigma(z^{-1})\quad\Rightarrow\;\begin{pmatrix}h_{1}(z,q)\\\
h_{3}(z,q)\end{pmatrix}\sim\begin{pmatrix}z^{\mu-1}\mskip
1.5mu\sigma(z^{-1})\\\ z^{\mu}\mskip 1.5mu\sigma(z^{-1})\end{pmatrix}\,,$
$\displaystyle f_{2}(z)$ $\displaystyle\sim
z^{\mu}\,\sigma(z^{-1})\quad\Rightarrow\;\begin{pmatrix}h_{1}(z,q)\\\
h_{3}(z,q)\end{pmatrix}\sim\begin{pmatrix}z^{\mu}\mskip 1.5mu\sigma(z^{-1})\\\
z^{\mu-1}\mskip 1.5mu\sigma(z^{-1})\end{pmatrix}\,.$ (3.10)

Both the functions $f_{i}$ and the solutions $h_{i}$ remain finite for
$z\to\infty$ provided that $\mu\leq 0$. Our method thus allows the computation
of Fourier-Bessel integrals not only for functions $f_{i}$ with a power-law
decrease but also for functions that tend to a finite value at $z\to\infty$.

For functions with an exponential falloff at large $z$, we get

$\displaystyle f_{1}(z)\text{ or }f_{2}(z)$ $\displaystyle\sim
z^{\mu}\,\sigma(z^{-1})\mskip 1.5mu\exp(-\kappa
z)\quad\Rightarrow\;\begin{pmatrix}h_{1}(z,q)\\\
h_{3}(z,q)\end{pmatrix}\sim\begin{pmatrix}z^{\mu}\mskip 1.5mu\sigma(z^{-1})\\\
z^{\mu}\mskip 1.5mu\sigma(z^{-1})\end{pmatrix}\exp(-\kappa z)$ (3.11)

for $\kappa>0$ according to section 2.3.2, whereas for functions with a
Gaussian falloff, the results of appendix B imply that

$\displaystyle f_{1}(z)$ $\displaystyle\sim z^{\mu}\,\sigma(z^{-1})\mskip
1.5mu\exp(-\lambda^{2}z^{2}-\kappa z\mskip
1.5mu\bigr{)}\;\quad\Rightarrow\;\begin{pmatrix}h_{1}(z,q)\\\
h_{3}(z,q)\end{pmatrix}\sim\begin{pmatrix}z^{\mu-1}\mskip
1.5mu\sigma(z^{-1})\\\ z^{\mu-2}\mskip
1.5mu\sigma(z^{-1})\end{pmatrix}\exp(-\lambda^{2}z^{2}-\kappa z\mskip
1.5mu\bigr{)}\;\,,$ $\displaystyle f_{2}(z)$ $\displaystyle\sim
z^{\mu}\,\sigma(z^{-1})\mskip 1.5mu\exp(-\lambda^{2}z^{2}-\kappa z\mskip
1.5mu\bigr{)}\;\quad\Rightarrow\;\begin{pmatrix}h_{1}(z,q)\\\
h_{3}(z,q)\end{pmatrix}\sim\begin{pmatrix}z^{\mu-2}\mskip
1.5mu\sigma(z^{-1})\\\ z^{\mu-1}\mskip
1.5mu\sigma(z^{-1})\end{pmatrix}\exp(-\lambda^{2}z^{2}-\kappa z\mskip
1.5mu\bigr{)}\;$ (3.12)

for $\lambda>0$. Both in the exponential and in the Gaussian case, the
solutions $h_{1}$ and $h_{3}$ thus have the same leading behaviour for $z\to
0$ as the functions $f_{1}$ or $f_{2}$, except for a possible downward shift
of the power $z^{\mu}$ by one or two units.

##### Using integration by parts.

It remains to discuss the case where one uses integration by parts (3) and
then evaluates the integral on its r.h.s. by Levin’s method.

With (3.6) we find

$\displaystyle f_{0}(z)$ $\displaystyle\sim
z^{\mu+\nu-1}\,\sigma(z)\quad\Rightarrow\;f_{1}(z)\sim
z^{\mu+\nu-1}\,\sigma(z)$ (3.13)

for $z\to 0$, so that our method can be applied to integrands behaving like
$z^{\mu}\mskip 1.5muJ_{\rho}(qz)$ with $\mu+\rho\geq 0$ for all $\rho\geq 0$.

At large $z$ we have

$\displaystyle f_{0}(z)$ $\displaystyle\sim z^{\mu}\,\sigma(z^{-1})\mskip
1.5mu\exp(-\lambda^{2}z^{2}-\kappa z\mskip
1.5mu\bigr{)}\;\quad\Rightarrow\;f_{1}(z)\sim
z^{\rho}\,\sigma(z^{-1})\,\exp(-\lambda^{2}z^{2}-\kappa z\mskip
1.5mu\bigr{)}\;$ (3.14)

with

$\displaystyle\rho$ $\displaystyle=\begin{cases}\mu-1&\text{ if
}\lambda=\kappa=0\,,\\\ \mu&\text{ if }\lambda=0\text{ and }\kappa>0\,,\\\
\mu+1&\text{ if }\lambda>0\,.\end{cases}$ (3.15)

One thus gets the same leading exponential or Gaussian behaviour for $f_{0}$
and $f_{1}$. A pure power law $f_{0}\sim z^{\mu}$ with $\mu\leq 0$ implies a
power-law decrease $f_{1}\sim z^{\mu-1}$.

##### Grand summary.

As an upshot of our detailed discussion, we retain that if a set of basis
functions is good for approximating the functions $f_{1}$ or $f_{2}$ in (3.2),
or $f_{0}$ in (3), it can be expected to be good for approximating the
solutions $h_{1}$ and $h_{3}$ of the Levin equations (3) and hence allow for
an accurate computation of the integrals. A small-$z$ behaviour like
$z^{1-\delta}$ with $\delta<1$ is tolerable for $f_{0}$ or $f_{1}$, because
the corresponding power of $z$ for $h_{1}$ and $h_{3}$ is larger by one unit.

###  3.2 Discretisation on Chebyshev grids

We now discuss in detail how we implement the method of collocation for
solving the Levin equations (3).

In most examples of Levin’s original work [12, 13], the chosen basis functions
are polynomials $b_{i}(z)=(z-z_{0})^{i}$, and the system of ODEs is evaluated
at _equispaced_ points between $z_{a}$ and $z_{b}$. It was found that the
accuracy of the method could be improved by dividing the integration domain
into several subintervals. Other options for the basis functions are mentioned
in [12].

It is well known that polynomial approximation on equispaced grids suffers
from Runge’s phenomenon: the approximation becomes poor in the vicinity of the
interval end points when one increases the number of grid points (i.e. the
polynomial order). This problem is avoided by using suitable grids that are
not equispaced. A prominent example is _Chebyshev interpolation_ , which is
described in detail in [17]. This was indeed used in the context of Levin’s
method in [14, 15].

Motivated by our experience with interpolating (single and double) parton
distributions on Chebyshev grids [18, 19], we adopt Chebyshev interpolation to
our present case. Our method involves three steps.

##### Variable transform.

We first define a transformation from $z$ to the variable $u$ that will be
used for interpolation with Chebyshev polynomials. The values $u_{a}=u(z_{a})$
and $u_{b}=u(z_{b})$ must both be finite. We require $u(z)$ to be an analytic
function with $du/dz>0$ for all $z\in[z_{a},z_{b}]$, except for the point
$z_{b}=\infty$, where necessarily $du/dz=0$.

A primary purpose of this transformation is to map an infinite integration
interval $[z_{a},\infty]$ onto a finite interval in $u$. On a finite interval
in $z$, one may take the trivial transformation $u=z$. However, even in this
case a non-trivial variable transformation can be beneficial, for instance
when $f_{i}(z)$ rises or decreases very steeply.

Depending on the behaviour of $f_{i}(z)$, different variable transformations
$u(z)$ lead to good results of our method. A selection of such transformations
is presented in section 3.2.2.

##### Subintervals.

We find that it is often useful to split $[z_{a},z_{b}]$ into $k$ subintervals
$[z_{0},z_{1}]$, $[z_{1},z_{2}]$, …, $[z_{k-1},z_{k}]$ with $z_{0}=z_{a}$ and
$z_{k}=z_{b}$ and to evaluate the integral on each of them separately. We will
see that this benefits both the integration accuracy and the computation time.
Following the notation introduced in [18, 19], we write

$\displaystyle[z_{0},z_{1},\ldots,z_{k}]_{(n_{1},\ldots,n_{k})}$ (3.16)

to characterise a grid and its subgrids, where $z_{j}$ are the interval
boundaries, and $n_{j}$ is the number of points on each subgrid. The total
number of grid points is then $\sum_{j=1}^{k}n_{j}-(k-1)$, because adjacent
subgrids share their end points. In the notation (3.16) it is understood that
one and the same variable transform is used on the full interval from $z_{0}$
to $z_{k}$.

For simplicity, we will continue to present our method for a single
integration interval $[z_{a},z_{b}]$, bearing in mind that this may be a
subinterval in the sense just described.

##### Chebyshev interpolation.

On each subinterval, we approximate the functions $h_{1}(z,q)$ and
$h_{3}(z,q)$ using Chebyshev interpolation. A compendium of this method with
the formulae needed for our method is given in appendix C. In brief, we
approximate

$\displaystyle h_{i}(z,q)$ $\displaystyle\approx
p_{i}\bigl{(}u(z),q\bigr{)}\,,$ (3.17)

where $p_{i}(u,q)$ is the unique polynomial of order $N$ that reproduces the
function

$\displaystyle h_{i}(z_{j},q)$ $\displaystyle=p_{i}\bigl{(}u_{j},q\bigr{)}\,,$
$\displaystyle u_{j}$ $\displaystyle=u(z_{j})$ (3.18)

at the $n=N+1$ points $z_{j}$ that correspond to the Chebyshev points $u_{j}$
on the interval $[u_{a},u_{b}]$, with $u_{j}$ given in (C.11). Note that these
points include the interval boundaries, $z_{a}=z_{0}$ and $z_{b}=z_{N}$, where
the functions $h_{i}$ are required for evaluating the Fourier-Bessel integral
with the master formula (3.3) of our method.

The derivatives of $dh_{i}/dz$ appearing in the Levin equations (3) can be
approximated by the derivatives of the interpolants $dp_{i}/dz$. When
evaluating the equations at the points $z_{j}$, these derivatives can be
evaluated via a matrix multiplication. In summary, we discretise (3) by making
the following replacements:

$\displaystyle f_{i}(z)$ $\displaystyle\to f_{i}(z_{j})$
$\displaystyle(i=1,2)\,,$ $\displaystyle h_{i}(z,q)$ $\displaystyle\to
p_{i}(u_{j},q)$ $\displaystyle(i=1,3)\,,$ $\displaystyle\frac{d}{dz}\mskip
1.5muh_{i}(z,q)$
$\displaystyle\to\frac{du}{dz}(z_{j})\sum_{k}D^{u}_{jk}\,p_{i}(u_{k},q)$
$\displaystyle(i=1,3)\,,$ (3.19)

where $D^{u}_{jk}$ is given in (C.12) and $du/dz$ is the Jacobian of the
variable transform. The Levin equations are thus replaced by a system of
linear equations, which can be written as

$\displaystyle B(q)\cdot\boldsymbol{P}(q)$ $\displaystyle=\boldsymbol{F}$
(3.20)

with the $2n$-dimensional vectors

$\displaystyle\boldsymbol{F}$
$\displaystyle=\bigl{(}f_{1}(z_{0}),f_{1}(z_{1}),\ldots,f_{1}(z_{N}),\;f_{2}(z_{0}),f_{2}(z_{1}),\ldots,f_{2}(z_{N})\bigr{)}^{T}\,,$
$\displaystyle\boldsymbol{P}(q)$
$\displaystyle=\bigl{(}p_{1}(z_{0},q),p_{1}(z_{1},q),\ldots,p_{1}(z_{N},q),\;p_{3}(z_{0},q),p_{3}(z_{1},q),\ldots,p_{3}(z_{N},q)\bigr{)}^{T}$
(3.21)

and the $2n$-dimensional square matrix

$\displaystyle B(q)$ $\displaystyle=\begin{pmatrix}B^{aa}&B^{ab}(q)\\\
B^{ba}(q)&B^{bb}\end{pmatrix}$ (3.22)

with diagonal blocks

$\displaystyle B^{aa}_{jk}$
$\displaystyle=\frac{du}{dz}(z_{j})\,D^{u}_{jk}+\frac{\nu}{1+z_{j}}\,\delta_{jk}\,,$
$\displaystyle B^{bb}_{jk}$
$\displaystyle=\frac{z_{j}}{1+z_{j}}\,\frac{du}{dz}(z_{j})\,D^{u}_{jk}-\biggl{[}\mskip
1.5mu\frac{\nu-1}{(1+z_{j})^{2}}+\frac{\nu+1}{1+z_{j}}\mskip
1.5mu\biggr{]}\mskip 1.5mu\delta_{jk}$ (3.23)

and off-diagonal blocks

$\displaystyle B^{ab}_{jk}(q)$
$\displaystyle=q\,\frac{z_{j}}{1+z_{j}}\,\delta_{jk}\,,$ $\displaystyle
B^{ba}_{jk}(q)$ $\displaystyle=-q\mskip 1.5mu\delta_{jk}\,,$ (3.24)

where $j$ and $k$ run from $0$ to $N$ in all cases.

Note that up to this point, the only approximation we have made was to replace
the derivatives of $h_{i}$ at the points $z_{j}$ by the derivatives of their
approximants $p_{i}$. If we use the integration-by-parts identity (3), then we
also take Chebyshev interpolation to approximately compute the function
$f_{1}$ as

$\displaystyle f_{1}(z_{j})$
$\displaystyle\approx\sum_{k}C_{jk}\,f_{0}(z_{k})\,,$ $\displaystyle C_{jk}$
$\displaystyle=\frac{z_{j}}{1+z_{j}}\mskip
1.5mu\frac{du}{dz}(z_{j})D^{u}_{jk}-\biggl{[}\frac{\nu-1}{(1+z_{i})^{2}}+\frac{\nu}{1+z_{i}}\biggr{]}\mskip
1.5mu\delta_{jk}\,,$ (3.25)

which is the discretised version of (3.6).

In section 3.1 we investigated in detail how the small-$z$ behaviour of
$h_{1}$ and $h_{3}$ is related to the small-$z$ behaviour of $f_{1}$ or
$f_{2}$. Since we require the variable transformation between $u$ and $z$ to
be analytic and to have a finite derivative at $z=0$, the interpolation
polynomials $p_{i}\bigl{(}u(z),q\bigr{)}$ in $u$ have a Taylor expansion
around $z=0$. To obtain a good approximation when discretising the Levin
equations, we hence need that $h_{1}$ and $h_{3}$ can be sufficiently well
approximated by a Taylor series around $z=0$, with the first $n$ Taylor
coefficients being free and the remaining ones fixed by the variable
transform.

A corresponding statement holds for the small-$z$ behaviour of $f_{0}(z)$ when
using integration by parts. Notice that if one has $f_{0}\sim z^{1-\delta}$
with $0<\delta<1$, one gets a divergent first derivative of $f_{0}$ at $z=0$,
which is however multiplied with $z$ in (3.6) and thus can be accommodated
within our setup.

##### Using quadrature.

We will see that the matrix $B$ in (3.20) becomes ill-conditioned when $q$
becomes small. As described in the next subsection, this can in part be dealt
with by linear algebra methods. On the other hand, there is no compelling need
to use the Levin method when the integrand in (3.2) is not rapidly oscillating
in the integration interval $[z_{a},z_{b}]$. In our algorithm, we therefore
replace the Levin method with Clenshaw-Curtis quadrature if for a given $q$
one has

$\displaystyle qz_{b}$ $\displaystyle\leq j_{\nu}\,,$ (3.26)

where $j_{\nu}$ is the first zero of the Bessel function $J_{\nu}(x)$. We find
that this yields a high integration accuracy and also saves computing time.
The quadrature result is obtained as the sum

$\displaystyle I(q)$
$\displaystyle\approx\sum_{j}w^{z}_{j}\;\biggl{(}\frac{1+z_{j}}{z_{j}}\biggr{)}^{\nu}\,\biggl{[}J_{\nu}(qz_{j})\,f_{1}(z_{j})+J_{\nu+1}(qz_{j})\,f_{2}(z_{j})\biggr{]}\,,$
(3.27)

where $I(q)$ is defined in (3.2) and the weights are given by

$\displaystyle w^{z}_{j}$
$\displaystyle=\biggl{(}\frac{du}{dz}(z_{j})\biggr{)}^{-1}\;w^{u}_{j}$ (3.28)

with $w^{u}_{j}$ from (C.13).

####  3.2.1 Solving the linear equation system

The system (3.20) has a unique solution if the matrix $B$ is nonsingular. By
taking increasingly fine interpolation grids, i.e. increasingly high order of
the interpolation polynomials $p_{i}(u,q)$, one should eventually approach the
limit of the original Levin equations, which has an infinity of solutions. One
can hence expect that for sufficiently fine interpolation grids, the matrix
$B$ is close to singular. This poses a challenge to the numerical solution of
(3.20).

The work in [14, 15] contains a detailed analysis of the linear equation
system for a different set of integrals, along with a strategy for its
solution. We follow a rather similar strategy in the present work, selecting
one out of two algorithms to solve (3.20).

Our default algorithm uses the LU decomposition (see e.g. [20, Chapter 2.3]),
which represents $B$ in the form

$\displaystyle B$ $\displaystyle=\Pi\cdot L\cdot U\,,$ (3.29)

where $\Pi$ is a permutation matrix, $L$ is a lower diagonal matrix with all
diagonal elements equal to $1$, and $U$ an upper diagonal matrix. The product
of the diagonal elements $U_{ii}$ is equal to $\det(B)$, so that at least one
of these elements is zero when $B$ is singular. If this is not the case, the
solution of $B\cdot\boldsymbol{P}=\boldsymbol{F}$ is obtained by successively
solving

$\displaystyle L\cdot\boldsymbol{G}$
$\displaystyle=\Pi^{-1}\cdot\boldsymbol{F}\,,$ $\displaystyle
U\cdot\boldsymbol{P}$ $\displaystyle=\boldsymbol{G}\,,$ (3.30)

which can be done by forward and backward substitution, respectively. The
second step involves the inverse diagonal elements $1/U_{ii}$ of $U$ and
becomes numerically unstable when $B$ is close to singular. As criterion for
this case, we use

$\displaystyle\min_{i}\bigl{(}|U_{ii}|\bigr{)}$
$\displaystyle\;\leq\;r_{LU}^{\text{max}}\,\max_{i}\bigl{(}|U_{ii}|\bigr{)}\,,$
(3.31)

where $r_{LU}^{\text{max}}$ can be chosen by the user.

If (3.31) is satisfied for a given matrix $B$, we use the singular value (SV)
decomposition. In this case, $B$ is represented as

$\displaystyle B$ $\displaystyle=U\cdot S\cdot V^{T}\,,$ (3.32)

where $U$ and $V$ are orthogonal444Note that $U$ denotes different matrices in
(3.32) and (3.29). and $S=\operatorname{diag}(S_{ii})$ is diagonal and
positive semidefinite.

For nonsingular $B$, the solution of $B\cdot\boldsymbol{P}=\boldsymbol{F}$ is
then readily given by

$\displaystyle\boldsymbol{P}$
$\displaystyle=V\cdot\operatorname{diag}(1/S_{ii})\cdot
U^{T}\cdot\boldsymbol{F}\,.$ (3.33)

If $B$ is singular, then some of the $S_{ii}$ vanish, and one replaces the
corresponding entries $1/S_{ii}$ in (3.33) with zero. As shown e.g. in [20,
Chapter 2.6], this yields an approximate solution of (3.20) in the sense that
$\boldsymbol{P}$ obtained in this way minimises
$|B\cdot\boldsymbol{P}-\boldsymbol{F}|$, i.e. the distance between the vectors
on the two sides of the equation. The numerical stability of this procedure is
greatly enhanced by also replacing $1/S_{ii}$ with zero when $S_{ii}$ is very
small; this is called a _truncated_ SV decomposition. We choose to replace an
entry $1/S_{ii}$ with zero if it satisfies

$\displaystyle S_{ii}$
$\displaystyle\;<\;r_{SV}^{\text{max}}\,\max_{j}(S_{jj})\,,$ (3.34)

with $r_{SV}^{\text{max}}$ to be chosen by the user.

A numerical computation is typically much faster for the LU decomposition
(3.29) than for the SV decomposition (3.32). One should hence take
$r_{LU}^{\text{max}}$ as small as can be done without compromising the
accuracy of the result. Our choice of the control parameters
$r_{LU}^{\text{max}}$ and $r_{SV}^{\text{max}}$ is discussed in section 4.1.

The computation time for both the LU and the SV decomposition grows very
quickly with the dimension of the matrix $B$. From this point of view, it is
beneficial to split the overall integration region $[z_{a},z_{b}]$ into
subintervals, since this allows one to work with smaller matrices $B_{k}$ on
each subinterval $k$.

####  3.2.2 Examples of variable transforms

A selection of variable transformations that are suitable for integration over
$z$ from $0$ to $\infty$ is presented in table 1. The table has been adapted
from table 1 in [19], where the same transformations were used to describe the
dependence of double parton distributions on the transverse distance between
the two partons.

The following brief characterisation of the different transforms is the result
of detailed numerical investigations, part of which will be presented in the
next sections.

  * •

The inverse power law transformation in the first row of the table is suitable
for integrals with a large-$z$ behaviour like $f(z)\sim z^{\mu}$ with $\mu\leq
0$.

  * •

The exponential transform in the second row works well for functions that have
a Gaussian behaviour $f(z)\sim\exp(-\lambda^{2}z^{2})$ at large $z$, possibly
modified by an exponential $\exp(-\kappa z)$ times a power of $z$.

  * •

The exponential transform with a square root in the third row works well for
$f(z)\sim\exp(-\kappa z)$ at large $z$, possibly modified by a power of $z$.

In our study [19] of double parton distributions, we have shown that the
Gaussian transformation in the last row yields good _interpolation_ accuracy
for functions with a Gaussian behaviour at large $z$. Likewise, we found that
the exponential transform in the second row allows for an accurate
interpolation of functions with an exponential falloff. As we will see in
section 5, the same combinations of variable transforms and integrands work
with Levin’s method but are less powerful than the combinations specified
above.

shorthand | ${}-u(z)$ | $z(u)$ | $du/dz$
---|---|---|---
inv pow | $(z+z_{0})^{-\alpha}$ | $|u|^{-1/\alpha}-z_{0}$ | $\alpha\,(z+z_{0})^{-1-\alpha}$
exp | $\exp\Bigl{[}-\frac{1}{4}\mskip 1.5mumz\Bigr{]}$ | $\frac{4}{m}\mskip 1.5muL(u)$ | $\frac{1}{4}\mskip 1.5mum|u|$
exp sqrt | $\exp\Bigl{[}\mskip 1.5mu1-\sqrt{1+mz/2}\,\Bigr{]}$ | $\frac{2}{m}\mskip 1.5mu\Bigl{[}L^{2}(u)+2L(u)\mskip 1.5mu\Bigr{]}$ | $\frac{1}{4}\mskip 1.5mum|u|\,\Bigl{[}L(u)+1\mskip 1.5mu\Bigr{]}^{-1}$
Gauss | $\exp\Bigl{[}-\frac{1}{4}\mskip 1.5mu(m^{2}z^{2}+mz)\Bigr{]}$ | $\frac{1}{2m}\,\bigl{[}\sqrt{16L(u)+1}-1\mskip 1.5mu\bigr{]}$ | $\frac{1}{4}\mskip 1.5mum|u|\,\sqrt{16L(u)+1}$
Table 1: Variable transformations suitable for integration over $z$ from $0$
to $\infty$. It is understood that $\alpha>0$, $z_{0}>0$, and $m>0$. All
transformations satisfy $u\leq 0$ and $du/dz>0$ for $u<0$. At $u=0$ one has
$z=\infty$ and $du/dz=0$. We abbreviate $L(u)=\ln\bigl{(}1/|u|\bigr{)}$.

###  3.3 Summary of the algorithm

Our complete algorithm naturally proceeds in several steps.

  1. 1.

The user must specify an interpolation grid, i.e. a variable transform $u(z)$,
a subinterval division in $z$, and the number of grid points per subinterval.
Then

    1. (a)

compute the grid points $z_{j}$,

    2. (b)

evaluate the values at the grid points for the function(s) to be integrated
over, i.e. for $f_{1}(z),f_{2}(z)$ or $f_{0}(z)$ in (3.2) or (3),

    3. (c)

perform steps 2 to 4 for each subinterval and add up the integrals from these
subintervals.

  2. 2.

For a given value of $\nu$

    1. (a)

compute the values of $\bigl{[}z_{j}/(1+z_{j})\bigr{]}^{\rho}$ with $\rho=1$,
$\nu$, and $\nu-1$ for each grid point,

    2. (b)

compute the diagonal blocks (3.2) of the matrix $B$, as well as the matrix $C$
in (3.2),

    3. (c)

compute the Clenshaw-Curtis weights $w^{z}_{j}$ in (3.28).

We recall that for a given $\nu$, the method allows the computation of
integrals with $J_{\nu-1}(qz)$, $J_{\nu}(qz)$, and $J_{\nu+1}(qz)$.

  3. 3.

For a given value of $q$

    1. (a)

determine from (3.26) whether to use quadrature or the Levin method,

    2. (b)

compute the values of $\bigl{[}(1+z)/z\mskip 1.5mu\bigr{]}^{\nu}\,J_{\nu}(qz)$
and $\bigl{[}(1+z)/z\mskip 1.5mu\bigr{]}^{\nu}\,J_{\nu+1}(qz)$, for all grid
points $z_{j}$ if quadrature is used, and for the interval boundaries
otherwise,

    3. (c)

if the Levin method is used: compute the off-diagonal blocks (3.24) of the
matrix $B$,

    4. (d)

if the Levin method is used: set up the linear equation solver by computing
the LU and — if needed — the SV decomposition of $B$.

  4. 4.

For given values $f_{i}(z_{j})$ on the interpolation grid ($i=0,1$, or $2$)

    1. (a)

for integrals with $J_{\nu-1}(qz)\mskip 1.5mu$: evaluate the boundary terms at
$z_{a}$ and $z_{b}$ in (3) and compute $f_{1}(z_{j})$ from $f_{0}(z_{j})$
using (3.2),

    2. (b)

if the Levin method is used: solve the matrix equation (3.20),

    3. (c)

evaluate the integral from the master formula (3.3) or (3).

The time required for step 1b crucially depends on the functions to be
integrated; an important feature of our method is that the functions need not
be evaluated separately for different values of $q$. The evaluation of Bessel
functions is limited to step 3b and required only once for a given $q$.

The operations in step 2 need to be done only once for a given choice of $\nu$
and the interpolation grid – they should be irrelevant for the overall
computing time in a typical workflow. We find step 3d to be by far the most
time consuming one. By comparison, all operations in step 4 are fast.

###  3.4 Implementation in BestLime

We have implemented this algorithm in the C++ library BestLime, which is
available on [16]. It provides a number of variable transforms $u(z)$,
including those in table 1, as well as the possibility to work with user-
defined transformations. In addition to routines to compute the integrals in
(3.2) or (3) for given discretised functions $f_{1}$, $f_{2}$ or $f_{0}$, it
contains analogous calls for integrands $J_{\rho}(qz)\,\tilde{f}(z)$ with
$\rho=\nu$, $\nu+1$, or $\nu-1$. These calls simply multiply $\tilde{f}(z)$
with $[z/(1+z)]^{\nu}$ or $[z/(1+z)]^{\nu-1}$ and then use the master formulae
(3.2) or (3).

For linear algebra operations, BestLime uses the Eigen library [21],
specifically the part for dense matrices in Eigen version 3.4.0. In
particular, the linear equations are solved using the methods PartialPivLU (LU
decomposition with partial pivoting) or BDCSVD (Bidiagonal Divide and Conquer
SVD).

For the computation of Bessel functions (and of the hypergeometric function
${}_{2}F_{1}$ in the benchmark comparison of section 4.2), the GNU Scientific
Library (GSL) [22] is used.

The methods of BestLime that deal with Chebyshev interpolation have been taken
form the ChiliPDF library [18, 19], which is under development.

##  4 Numerical studies

In this section, we investigate several numerical and computational aspects of
our integration algorithm. To characterise the interpolation grids, we use the
shorthand names for variable transforms in the first column of table 1, along
with the notation for grid boundaries and grid points in equation (3.16).

All results in this and the next sections have been obtained with the library
BestLime. We begin with a brief discussion of the parameters that control the
algorithm used for solving the discretised Levin equations (3.20). After this,
we present benchmarks for the integration accuracy. Finally, we give an
indication of the computation time required by BestLime for representative
integration setups.

###  4.1 Control parameters for solving the linear equation system

As explained in section 3.2.1, we use the LU decomposition (3.29) to solve the
discretised Levin equations, provided that the matrix $B$ is sufficiently far
from being singular. A measure for this is the ratio

$\displaystyle r_{LU}$
$\displaystyle=\frac{\min_{i}\bigl{(}|U_{ii}|\bigr{)}}{\max_{i}\bigl{(}|U_{ii}|\bigr{)}}\,,$
(4.1)

where $U_{ii}$ denotes the diagonal elements of the matrix $U$ in (3.29). Let
us see when this ratio becomes very small.

In the following exercise, we consider grids of the form
$[0,z_{1},\infty]_{(n_{1},n_{2})}$ with $(n_{1},n_{2})$ equal to $(16,32)$ or
$(32,64)$ . We use either the exp sqrt transform specified in table 1 with the
parameter $m$ ranging from $0.32$ to $3.85$, or the exp transform with $m$
ranging from $0.19$ to $2.24$. (These settings correspond to the parameter
scans presented in section 5.2.3.) Apart from some details, our results are
very similar for different values of $\nu$, and in figures 2 and 3 we show
them for $\nu=1$ and $\nu=2$, respectively.

We compute $B$ for a wide range of $z_{1}$ and $q$ values, with linearly
spaced $z_{1}$ values between $0.03$ and $1$ and logarithmically spaced $q$
values between $0.01$ and $100$. We find $r_{LU}$ to be rather strongly
correlated with the product $qz_{1}$ of these variables, as shown for the
lower subinterval $[0,z_{1}]$ in the top panels of figures 2 and 3 and for the
upper subinterval $[z_{1},\infty]$ in the bottom panels of these figures.

(a)

(b)

(c)

(d)

Figure 2:  The ratio $r_{LU}$ defined in (4.1) for the matrix $B$ on a given
subgrid. The plots in this figure are for $\nu=1$. See the text for details of
the variable transform $u(z)$.

(a)

(b)

(c)

(d)

Figure 3:  As figure 2, but for $\nu=2$.

On the upper subinterval, $r_{LU}$ is always relatively large, decreasing only
slightly when the number of grid points is increased. One can expect the LU
decomposition to give reliable numerical results in this range, and our
integration benchmarks confirm that this is indeed the case.

The situation is very different in the lower subinterval. Here $r_{LU}$ drops
sharply as $qz_{1}$ decreases to values of order $1$, and it remains very
small when $qz_{1}$ decreases even further. This confirms our expectation
given in section 3.2.1: in a $z$ interval where the Bessel functions
$J_{\nu}(qz)$ and $J_{\nu+1}(qz)$ have only few or no oscillations, the
discretised version of the Levin equations approaches the “continuum limit”,
in which these equations have an infinite number of solutions.

We checked separately that with a single subinterval $[0,\infty]$ the value of
$r_{LU}$ stays above $10^{-6}$ for $n\leq 96$ points with the variable
transformations and $q$ values specified above. For a single interval from $0$
to $\infty$, one can hence expect the LU method to work well throughout.

When $r_{LU}\leq r_{LU}^{\text{max}}$, we no longer use the LU decomposition
to solve the linear equations and resort to the computationally more stable
(and more expensive) SV decomposition. When the latter is used, we also need
to choose the threshold $r_{SV}^{\text{max}}$ for truncating small singular
values as specified in (3.34).

To set these two control parameters of our method, we disabled the use of
quadrature in intervals with $qz_{1}\leq j_{\nu}$ (see the discussion around
(3.26)) and monitored how the integration accuracy depends on
$r_{LU}^{\text{max}}$ and $r_{SV}^{\text{max}}$. We performed this monitoring
for the integrals in sections 4.2 and 5 and found

$\displaystyle r_{LU}^{\text{max}}$
$\displaystyle=r_{SV}^{\text{max}}=10^{-12}$ (4.2)

to be a satisfactory setting. Our choice for $r_{LU}^{\text{max}}$ is
conservative in the sense that decreasing it by two orders of magnitude does
not degrade the accuracy in the test cases just mentioned.

When the use of quadrature in intervals with $qz_{1}\leq j_{\nu}$ is enabled,
the integration accuracy is even less dependent on the precise settings of
$r_{LU}^{\text{max}}$ and $r_{SV}^{\text{max}}$. This is because small values
of $r_{LU}$ or $r_{SV}$ tend to appear only for small values of $qz$, when
quadrature can be used.

We find that the only settings where the SV decomposition is needed involve
very fine subgrids and intermediate values of $q$ (small enough to give a
small $r_{LU}$ but large enough to fail our criterion for using quadrature).
An example is given in the rightmost column of table 4.

###  4.2 Benchmarks: grids and precision

To assess the performance of our algorithm, we compute a number of Fourier-
Bessel integrals whose exact form is known analytically. The integrals
considered are

$\displaystyle I(q)$
$\displaystyle=\int_{0}^{z_{b}}dz\,J_{\nu}(qz)\,\tilde{f}(z)$
$\displaystyle\text{ with }z_{b}=\begin{cases}\infty&\text{ for cases 1 to
6a,}\\\ 10&\text{ for cases 6b and 7,}\end{cases}$ (4.3)

where $\tilde{f}(z)$ and $I(q)$ for the different cases are given in table 3.
Unless specified otherwise in table 3, the integrals are evaluated for

$\displaystyle\nu$ $\displaystyle=\begin{cases}0,\,0.5,\,1,\,1.5,\,2&\text{
with $J_{\nu-1}$ calls,}\\\ 1,\,1.5,\,2,\,2.5,\,3&\text{ with $J_{\nu}$
calls,}\\\ 2,\,2.5,\,3&\text{ with $J_{\nu+1}$ calls,}\end{cases}$ (4.4)
$\displaystyle q$
$\displaystyle=0.01,\,0.1,\,1,\,2,\,3,\,5,\,10,\,15,\,20,\,25,\,30$ (4.5)

The different “calls” specified in (4.4) refer to the master equations (3.2)
and (3) of our method, i.e. the $J_{\nu-1}$, $J_{\nu}$, and $J_{\nu+1}$ calls
respectively correspond to a nonzero $f_{0}$, $f_{1}$, and $f_{2}$ in these
equations.

case | $\tilde{f}(z)$ | $I(q)$
---|---|---
1 | $z^{\mu}\,e^{-\kappa z}$ | $\frac{\Gamma(\mu+\nu+1)}{\Gamma(\nu+1)}\,\bigl{(}\frac{q}{2}\bigr{)}^{\nu}\,\kappa^{-\mu-\nu-1}\,{}_{2}F_{1}\Bigl{(}\frac{\mu+\nu+1}{2},\frac{\mu+\nu+2}{2};\nu+1;-\frac{q^{2}}{\kappa^{2}}\Bigr{)}$
2 | $e^{-\lambda^{2}z^{2}}$ | $\frac{\sqrt{\pi}}{2\lambda}\,\exp\Bigl{(}-\frac{q^{2}}{8\lambda^{2}}\Bigr{)}\,I_{\nu/2}\Bigl{(}\frac{q^{2}}{8\lambda^{2}}\Bigr{)}$
3 | $z^{\nu+1}\,e^{-\lambda^{2}z^{2}}$ | $\frac{1}{q}\,\bigl{(}\frac{q}{2\lambda^{2}}\bigr{)}^{\nu+1}\,\exp\Bigl{(}-\frac{q^{2}}{4\lambda^{2}}\Bigr{)}$
4a | $\bigl{(}\frac{z}{z^{2}+b^{2}}\bigr{)}^{\nu+1}$ | $\frac{1}{\Gamma(\nu+1)}\,\bigl{(}\frac{q}{2}\bigr{)}^{\nu}\,K_{0}(qb)$
4b | $z^{\nu+1}\big{/}(z^{2}+b^{2})^{\mu+\nu+1}$ | $\frac{1}{\Gamma(\mu+\nu+1)}\,\bigl{(}\frac{q}{2}\bigr{)}^{\nu}\,\bigl{(}\frac{q}{2b}\bigr{)}^{\mu}\,K_{\mu}(qb)$
5a | $1$ | $\frac{1}{q}$
5b | $z^{-\nu}$ | $\frac{\Gamma(1/2)}{2\mskip 1.5mu\Gamma(\nu+1/2)}\,\bigl{(}\frac{q}{2}\bigr{)}^{\nu-1}$
6a | $z^{-\nu+1}$ | $\frac{1}{2\mskip 1.5mu\Gamma(\nu)}\,\bigl{(}\frac{q}{2}\bigr{)}^{\nu-2}$
6b | $z^{-\nu+1}$ | $\frac{1}{2\mskip 1.5mu\Gamma(\nu)}\,\bigl{(}\frac{q}{2}\bigr{)}^{\nu-2}\,-\frac{1}{q}\,z_{b}^{-\nu+1}J_{\nu-1}(qz_{b})$
7 | $z^{\nu+1}$ | $\frac{1}{q}\,z_{b}^{\nu+1}\,J_{\nu+1}(qz_{b})$

Table 2:  The functions $\tilde{f}(z)$ used for our benchmark exercise and
their Fourier-Bessel transforms $I(q)$ specified in (4.3). Cases 4a, 5a, and
6a are special cases of 4b, 5b, and 6b, respectively, but are handled with
different integration settings (see table 3). Here $I_{\nu}$ and $K_{\nu}$ are
the modified Bessel functions of the first and second kind, and ${}_{2}F_{1}$
is the Gaussian hypergeometric function.

case | parameters | variable transform | remarks
---|---|---|---
1 | $\kappa=1.5,\,\mu=\kappa^{2}$ |  exp sqrt, $m=2\kappa$ |
2 | $\lambda=2$ |  exp, $m=4\lambda$ | not $\nu=0.5$
3 | $\lambda=2$ | ———— ” ———— | $q\leq 10$
4a | $b=1.2$ |  inv pow, $z_{0}=1,\alpha=0.5$ | $q\leq 10$
4b | $b=1.2,\,\mu=2.5$ |  inv pow, $z_{0}=1,\alpha=1$ | $q\leq 10$
5a |  |  inv pow, $z_{0}=1,\alpha=0.5$ | not $\nu=0.5$
5b |  | ———— ” ———— | only $\nu\geq 0.5$, not using $J_{\nu+1}$ call
6a |  | ———— ” ———— | only $\nu\geq 1$
6b |  | ———— ” ———— | also for $q=0.001,\,50,\,100,\,200,\,300$
7 |  | trivial transform $z=u$ | also for $q=0.001,\,50,\,100,\,200,\,300$

Table 3:  Details for computing the integrals in table 3. The variable
transforms are specified in table 1 and discussed in section 3.2.2.

We find grid settings both for intermediate or for high precision.
Specifically, we obtain a relative integration error $\varepsilon$

$\displaystyle\varepsilon$ $\displaystyle\leq 5\times 10^{-4}$
$\displaystyle\text{ with \ }[0,1,\infty]_{(20,25)}\text{ \ and \
}[0,\infty]_{(42)}\,,$ $\displaystyle\varepsilon$ $\displaystyle\leq 10^{-6}$
$\displaystyle\text{ with \ }[0,1,\infty]_{(34,44)}$ (4.6) for cases 1 to 6a,
whereas for cases 6b and 7 we obtain $\displaystyle\varepsilon$
$\displaystyle\leq 5\times 10^{-4}$ $\displaystyle\text{ with \
}[0,10]_{(24)}\,,$ $\displaystyle\varepsilon$ $\displaystyle\leq 10^{-6}$
$\displaystyle\text{ with \ }[0,10]_{(34)}\,.$ (4.7)

To find these settings, we did _not_ perform a dedicated fine tuning of the
interval boundary $z_{1}$ or of the grid transformation parameters given in
table 3. Instead, we decreased the number of grid points in each setting until
the accuracy goal of $5\times 10^{-4}$ or $10^{-6}$ was reached for all
integrals.

We note that in cases 4a and 4b the integral $I(q)$ decreases exponentially
with $q$, whilst in case 3 it decreases like a Gaussian in $q$. We limit the
benchmark to $q\leq 10$ in these cases. For larger value the accuracy degrades
rather quickly with our grid settings, which is not too surprising because the
integrals become extremely small as a result of numerical cancellations.

Notice also that we omit the value $\nu=0.5$ in our benchmark for cases 2 and
5a. In these cases, one has an integrand behaving like $z^{\mu}\mskip
1.5muJ_{\nu}(qz)$ with $\mu+\nu=0.5$ at small $z$, for which we expect our
method to be less accurate as discussed in section 3.1. This is borne out by
our numerical study: the relative integration error for $\nu=0.5$ in cases 2
and 5a is $\varepsilon\leq 3\times 10^{-3}$ for the grids
$[0,1,\infty]_{(20,25)}$ and $[0,\infty]_{(42)}$, and $\varepsilon\leq 2\times
10^{-4}$ for the grid $[0,1,\infty]_{(34,44)}$. Although this is worse than
the benchmark accuracy in (4.2), it shows that our method can handle these
cases as well.

###  4.3 Computation time

Indicative computation times needed for integration calls under different
conditions are given in table 4. The first two grids are those in our
benchmark exercise of the previous subsection; and the grid with
$(n_{1},n_{2})=(40,50)$ may be used for estimating the integration error of
the grid with $(n_{1},n_{2})=(20,25)$, as will be discussed in section 5.3.
The times in the table do _not_ include the time required for computing the
integrand on the grid points (step 1b in section 3.3), since that time
crucially depends on the function for the integrand and not on our algorithm.

As discussed in section 4.1, the integration method used in the first of two
subintervals depends on $q$. At sufficiently low $q$, quadrature can be used,
whilst for large $q$ the matrix $B$ becomes well conditioned and the LU
decomposition can be employed. For intermediate $q$ values, the LU method can
be used for less fine grids, but it will become unstable for finer grids, when
the SV decomposition must be taken instead.

grid | $[0,\infty]_{(42)}$ | $[0,1,\infty]_{(20,25)}$ | $[0,1,\infty]_{(40,50)}$
---|---|---|---
method | LU | (CC, LU) | (LU, LU) | (CC, LU) | (LU, LU) | (SV, LU)
$q$ range | low | low | intermediate | low | high | intermediate
new $q$ [$\operatorname{\mu s}$] | 36 | 12 | 17 | 60 | 92 | 1500
same $q$ [$\operatorname{\mu s}$] | 1.6 | 0.8 | 1.3 | 2.4 | 3.6 | 4.4
Table 4:  Indicative computation times of the integration call for $J_{\nu-1}$
(which includes integration by parts). The second row specifies the
integration method used in each subinterval. The $q$ ranges are $0.095\ldots
0.1$ (low), $14.25\ldots 15$ (intermediate), and $38\ldots 40$ (high). Timings
have been obtained on a recent laptop with an Intel® Core™ i7-1355U processor,
32 GiB memory, and code compiled with g++ version 11.4.0.

The _absolute_ times in the table should not be over-interpreted, since they
refer to a particular computer (where they fluctuate at the level of 20%
between different runs). Running the same code on a range of other laptops and
desktop PCs, we find that some times are larger than those in the table by a
factor up to 5, whilst in a few cases calls are faster by a factor up to 1.4.

Independently of the used computer, we observe a clear hierarchy for the
timing of calls under different conditions, and thus for the different
elements of our algorithm:

  * •

integration is significantly faster if the $q$ value is the same as in the
previous call,

  * •

times for integration at a new value of $q$ increase faster with the number of
grid points than times for integration at the same value of $q$,

  * •

quadrature is faster than the Levin method, and within the Levin method, an SV
decomposition takes more time than an LU decomposition (especially for the $q$
update),

  * •

for a comparable total number of grid points, using subgrids reduces the
computation time.

These findings are in line with our expectations formulated in sections 3.2.1
and 3.3. We also note that the time for initialising the integration routine
for a specific grid and $\nu$ value (step 2 in section 3.3) is well below
$100\operatorname{\mu s}$ in all cases. As expected, this is negligible for
the overall computing time budget.

The times in the table are for $J_{\nu-1}$ calls, i.e. they require
integration by parts and hence include the matrix multiplication (3.2).
Correspondingly, the integration calls for $J_{\nu}$ or $J_{\nu+1}$ are
slightly faster, namely between 10% to 20% for integration without update of
$q$.

##  5 Numerical studies for TMD-like functions

In this section, we explore in some detail how our method performs for
integrals that appear in TMD cross sections. Up to a global factor, such cross
sections are proportional to

$\displaystyle I(q)$ $\displaystyle=\int_{0}^{\infty}\\!dz\,J_{0}(qz)\,z\mskip
1.5muW(z,Q)$ (5.1)

where $W(z,Q)$ is product of TMDs evolved to the scale $Q$. For brevity, we do
not display the additional dependence on longitudinal momentum fractions. We
will also consider the cumulative $q$ spectrum

$\displaystyle K(q)$
$\displaystyle=2\pi\int_{0}^{q}d\tilde{q}\;\tilde{q}\,I(\tilde{q})=2\pi
q\int_{0}^{\infty}\\!dz\,J_{1}(qz)\,W(z,Q)\,.$ (5.2)

Notice that the integrands in (5.1) and (5.2) have the same behaviour at $z\to
0$.

We recall that $q$ is a measured transverse momentum and $z$ is a distance
between parton fields in the plane transverse to the collision axis in the
process. We will therefore give $q$ in units of $\operatorname{GeV}$ and $z$
in units in $\operatorname{GeV}^{-1}$ (using a system in which the speed of
light is unity).

###  5.1 Integrands

In this subsection, we specify different forms of $W(z,Q)$ used in our
numerical study. A reader not familiar with TMDs (but interested in Fourier-
Bessel transforms) may skip over the details, but should note two important
features of the integrand: $(i)$ it is non-analytic at $z=0$, and $(ii)$ its
leading behaviour at large $z$ is either given by an exponential $\exp(-\kappa
z)$ or by a Gaussian $\exp(-\lambda^{2}z^{2})$.

Our purpose is not a detailed phenomenological investigation of TMDs, but to
understand the interplay between their qualitative behaviour and our
integration algorithm. We therefore use leading-order evolution and tree-level
short-distance matching for the TMDs. For details and references, we refer to
[1, chapter 2]. The product of two evolved TMDs has the form

$\displaystyle W(z,Q)$
$\displaystyle=\bigl{[}f_{\text{np}}(z)\bigr{]}^{2}\,\exp\bigl{[}-2S(z,Q)\bigr{]}\,,$
(5.3)

where

$\displaystyle S(z,Q)$
$\displaystyle=\int_{\mu_{z}}^{Q}\frac{d\mu}{\mu}\,\biggl{[}4C_{F}\,\frac{\alpha_{s}(\mu)}{2\pi}\,\ln\frac{Q}{\mu}-3C_{F}\,\frac{\alpha_{s}(\mu)}{2\pi}\biggr{]}$
$\displaystyle=\biggl{[}\,\frac{c_{1}}{\alpha_{s}(Q)}-c_{2}\,\biggr{]}\,\ln\frac{\alpha_{s}(\mu_{z})}{\alpha_{s}(Q)}-c_{1}\biggl{[}\,\frac{1}{\alpha_{s}(Q)}-\frac{1}{\alpha_{s}(\mu_{z})}\,\biggr{]}$
(5.4)

with coefficients

$\displaystyle c_{1}$ $\displaystyle=8\pi C_{F}/\beta_{0}^{2}\,,$
$\displaystyle c_{2}$ $\displaystyle=3C_{F}/\beta_{0}\,,$
$\displaystyle\beta_{0}$ $\displaystyle=11-2n_{f}/3\,,$ $\displaystyle
C_{F}=4/3$ (5.5)

and with one-loop running of the strong coupling:

$\displaystyle\alpha_{s}(\mu)$
$\displaystyle=\frac{2\pi}{\beta_{0}\ln(\mu/\Lambda)}\,.$ (5.6)

The scale parameter $\Lambda$ is chosen such that $\alpha_{s}(M_{Z})=0.13$ at
the mass of the $Z$ boson. For simplicity, we take $n_{f}=5$ flavours
everywhere; a more realistic treatment would not affect the qualitative
behaviour of our integrands. The initial conditions for evolution are
formulated at the scale

$\displaystyle\mu_{z}$ $\displaystyle=b_{0}\mskip
1.5mu\sqrt{1/z^{2}+1/z^{2}_{\text{max}}\rule{0.0pt}{7.66385pt}}\,,$
$\displaystyle b_{0}$ $\displaystyle=2e^{-\gamma}\approx 1.12\,,$
$\displaystyle z_{\text{max}}$ $\displaystyle=0.5\operatorname{GeV}^{-1}\,.$
(5.7)

For $z\to 0$, we have a non-analytic behaviour in $z$, namely

$\displaystyle\exp\bigl{[}-2S(z,Q)\bigr{]}$
$\displaystyle\;\sim\;\biggl{(}\frac{z\Lambda}{b_{0}}\biggr{)}^{8C_{F}/\beta_{0}}\biggl{(}\ln\frac{b_{0}}{z\Lambda}\biggr{)}^{2c_{1}/\alpha_{s}(Q)}\approx\biggl{(}\frac{z\Lambda}{b_{0}}\biggr{)}^{1.4}\;\biggl{(}\ln\frac{b_{0}}{z\Lambda}\biggr{)}^{1.14/\alpha_{s}(Q)}\,.$
(5.8)

Note that if $f_{\text{np}}(z)$ is sufficiently regular at $z=0$, the function
$z\mskip 1.5muW(z)$ appearing in the integral (5.1) has a finite first
derivative at that point.

We performed all numerical studies for three values of $Q$, namely

$\displaystyle Q$
$\displaystyle=2\operatorname{GeV},\,20\operatorname{GeV},\,100\operatorname{GeV}\,,$
(5.9)

with corresponding values $\alpha_{s}(2\operatorname{GeV})\approx 0.33$,
$\alpha_{s}(20\operatorname{GeV})\approx 0.17$, and
$\alpha_{s}(100\operatorname{GeV})\approx 0.13$. The resulting power of
$\ln(1/z)$ in (5.8) is about $3.5$, $6.7$, and $8.9$, respectively and thus
rather high for the two larger $Q$ values.

For the TMDs at scale $\mu_{z}$, we consider two different forms that have
been used in phenomenological analyses. The first one is

$\displaystyle f_{\text{np, Yukawa}}(z)$
$\displaystyle=\cosh\biggl{[}\biggl{(}\frac{2c^{2}}{\kappa}-\frac{\kappa}{4}\biggr{)}z\,\biggr{]}\,\bigg{/}\cosh\biggl{[}\biggl{(}\frac{2c^{2}}{\kappa}+\frac{\kappa}{4}\biggr{)}z\,\biggr{]}\,,$
$\displaystyle\kappa$ $\displaystyle=0.642\operatorname{GeV}\,,$
$\displaystyle c$ $\displaystyle=0.521\operatorname{GeV}$ (5.10)

and corresponds to equation (3.7) and the first line of table 12 in [6]. At
large $z$ this has a Yukawa-type behaviour like $\exp(-\kappa z/2)$. The
second form has a Gaussian behaviour and reads

$\displaystyle f_{\text{np, Gauss}}(z)$
$\displaystyle=\bigl{(}1-c^{2}z^{2}\bigr{)}\exp\bigl{(}-\lambda^{2}z^{2}/2\bigr{)}\,,$
$\displaystyle\lambda$ $\displaystyle=0.374\operatorname{GeV}\,,$
$\displaystyle c$ $\displaystyle=0.117\operatorname{GeV}\,.$ (5.11)

This was used in the study [23], see equation (34) and table XI in that paper.
Notice that we define the parameters $\kappa$ and $\lambda$ such that the
_squared_ TMDs and hence the integrands in (5.1) and (5.2) decrease like
$\exp(-\kappa z)$ or $\exp(-\lambda^{2}z^{2})$.

In addition to these “Yukawa” and “Gauss” TMDs, we also consider the “toy TMD”
investigated in [8], which corresponds to

$\displaystyle z\mskip 1.5muW_{\text{toy}}(z)$
$\displaystyle=\frac{1}{\Gamma\bigl{(}\beta^{2}/\sigma^{2}\bigr{)}}\,\biggl{(}\frac{\beta
z}{\sigma^{2}}\biggr{)}^{\beta^{2}/\sigma^{2}}\exp\biggl{[}-\frac{\beta
z}{\sigma^{2}}\biggr{]}\,.$ (5.12)

This has the same form as the first entry in our table 3, with an exponential
decay parameter $\kappa=\beta/\sigma^{2}$. The function $W_{\text{toy}}(z)$
takes its maximum when $1/z$ is equal to

$\displaystyle Q$ $\displaystyle=\frac{\beta}{\beta^{2}-\sigma^{2}}\,,$ (5.13)

which in [8] was taken as a representative of the physical scale $Q$ of a hard
process in this simplified TMD setup. We choose the parameters in (5.12) as

$\displaystyle\beta/\sigma^{2}$
$\displaystyle=\kappa=0.642\operatorname{GeV}\,,$ $\displaystyle\beta$
$\displaystyle=1/\kappa+1/Q\,,$ (5.14)

thus taking the same value of $\kappa$ as in (5.1) for all $Q$.

Unless specified otherwise, the plots in the present section are obtained with
the following interpolation grids:

exp sqrt grid with $\displaystyle[0,z_{1},\infty]_{(16,32)}$ for the toy and
Yukawa TMDs, exp grid with $\displaystyle[0,z_{1},\infty]_{(16,32)}$ for the
Gauss TMD. (5.15)

It is natural to specify the mass parameters $m$ in the grid transformations
(see table 1) as multiples of the characteristic exponential decay parameter
of the integrands. We therefore define

$\displaystyle r_{\kappa}$ $\displaystyle=m/\kappa$ for the toy and Yukawa
TMDs, $\displaystyle r_{\lambda}$ $\displaystyle=m/\lambda$ for the Gauss TMD.
(5.16)

Our default grid parameters are

$\displaystyle z_{1}$ $\displaystyle=0.05\,,$ $\displaystyle r_{\kappa}$
$\displaystyle=3$ for the toy and Yukawa TMDs, $\displaystyle z_{1}$
$\displaystyle=0.05\,,$ $\displaystyle r_{\lambda}$ $\displaystyle=5$ for the
Gauss TMD (5.17)

and have been chosen such that one obtains good integration accuracy for all
three $Q$ values in (5.9), see table 5 and the plots in section 5.2.3.

The functions $z\mskip 1.5muW(z)$, which multiply $J_{0}(qz)$ in the integral
$I(q)$, are plotted on the l.h.s. of figure 4. When computing $I(q)$,
integration by parts is used as specified in (3), and the resulting function
$f_{1}(z)$ that appears in the master equation of the Levin method (3.2) is
shown on the r.h.s. of the figure. We see non-trivial structures of $f_{1}(z)$
appearing as a consequence of taking the derivative and weighting with powers
of $z/(1+z)$. The points on the curves for $f_{1}(z)$ correspond to the
interpolation grids just specified.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4:  Left: the transformed function $z\mskip 1.5muW(z)$ in (5.1) for
different values of $Q$ and the three choices of TMDs specified in the text.
Right: the corresponding function $f_{1}(z)$ in our method, see equations
(3.2) and (3). The points on the curves correspond to the interpolation grids
specified by (5.1) and (5.1).

The Fourier-Bessel transform $I(q)$, which gives the transverse-momentum
spectrum of the TMD process, is shown in figure 5. We observe that the
transform becomes negative at some value $q\sim Q$. This is a clear indication
that TMD factorisation does not describe the cross section in such a region.
We recall from the introduction that the evaluation of $I(q)$ may however be
necessary even in that region, in order to achieve a smooth transition between
the TMD regime and the regime in which the cross section can be computed in
terms of collinear parton distributions and fixed-order perturbative cross
sections.

(a)

(b)

(c)

Figure 5:  Absolute values of the Fourier-Bessel transform $I(q)$ for the
functions $z\mskip 1.5muW(z)$ shown on the l.h.s. of figure 4. Dashed lines
indicate that $I(q)$ is negative. The curve for $Q=100\operatorname{GeV}$ in
panel (a) has a zero crossing at $q$ between $100$ and
$101\operatorname{GeV}$.

###  5.2 Choosing interpolation grids

We now compare the performance of different settings for the interpolation
grids in our method. We generally find that high accuracy is easier to reach
for $K(q)$ than for $I(q)$, which is why we mostly focus on $I(q)$ in the
following.

To quantify the performance of a particular grid setting for a given choice of
$W(z)$, we evaluate the Fourier-Bessel integrals for a set of logarithmically
spaced values $q_{i}$ between $0.01\operatorname{GeV}$ and
$100\operatorname{GeV}$. We then determine the root mean square of the
relative integration error $\varepsilon$, averaging over all $q_{i}$ that are
not too close to the point $q_{0}$ where $I(q)$ has a zero crossing.
Specifically, we define

$\displaystyle\varepsilon^{2}_{\text{rms}}$
$\displaystyle=\langle\varepsilon^{2}\rangle_{q_{i}\notin R}$
$\displaystyle\text{ with }R$ $\displaystyle=(q_{0}-\Delta,q_{0}+\Delta)\,,$
(5.18)

where

$\displaystyle\Delta$ $\displaystyle=\begin{cases}1\operatorname{GeV}&\text{
for }Q=2\operatorname{GeV},\\\ 5\operatorname{GeV}&\text{ for
}Q=20\operatorname{GeV},\\\ 10\operatorname{GeV}&\text{ for
}Q=100\operatorname{GeV}.\end{cases}$ (5.19)

To compute the integration error, we take the analytic form of the Fourier-
Bessel transform for the toy TMD (see table 3) and the result of an adaptive
Gauss-Kronrod integrator for the Yukawa and Gauss TMDs.555We use the routine
scipy.integrate.quad of SciPy version 1.7.3, which employs the QAGS algorithm
of the QUADPACK library [24]. The relative target accuracy is set to
$10^{-12}$.

####  5.2.1 Variable transformations

We begin with a comparison between different variable transformations, fixing
the grid structure to the form $[0,z_{1},\infty]_{16,32}$. For each
transformation we need to choose parameters $(z_{1},r_{\kappa})$ or
$(z_{1},r_{\lambda})$. We do this by performing a parameter scan in these
values, with $z_{1}\in[0.01,0.3]$ in steps of $0.01$ and with $r_{\kappa}$ or
$r_{\lambda}$ varying from $0.25$ to $8$ in steps of $0.25$. The parameters
that give the smallest value of $\varepsilon_{\text{rms}}$ are listed in table
5.

toy TMD
---
| exp | exp sqrt
$Q[\operatorname{GeV}]$ | $(z_{1},r_{\kappa})$ | $r_{\kappa}$ | $(z_{1},r_{\kappa})$ | $(z_{1},z_{2},r_{\kappa})$
2 | (0.15, 1.0) | 7.0 |  1(0.1, 2.0) |  1(0.1, 1.0, 1.5)
20 | (0.15, 1.0) | 4.0 |  1(0.1, 2.0) |  1(0.1, 1.5, 1.5)
100 | (0.05, 1.5) | 6.0 | (0.05, 3.0) | 0.05, 0.5, 1.0)

Yukawa TMD
---
| exp | exp sqrt
$Q[\operatorname{GeV}]$ | $(z_{1},r_{\kappa})$ | $r_{\kappa}$ | $(z_{1},r_{\kappa})$ | $(z_{1},z_{2},r_{\kappa})$
2 |  1(0.1, 1.5) |  18.0 |  1(0.1, 2.0) |  1(0.1, 1.0, 1.5)
20 | (0.15, 1.5) |  18.0 |  1(0.1, 2.0) |  1(0.1, 0.8, 1.5)
100 | (0.05, 1.0) | 10.0 | (0.05, 4.0) | (0.05, 0.8, 6.0)

Gauss TMD
---
| Gauss | exp
$Q[\operatorname{GeV}]$ | $(z_{1},r_{\lambda})$ | $r_{\lambda}$ | $(z_{1},r_{\lambda})$ | $(z_{1},z_{2},r_{\lambda})$
2 | (0.15, 1.5) | 8.0 |  1(0.1, 3.0) |  1(0.2, 3.0, 5.0)
20 | (0.15, 1.5) | 7.0 |  1(0.1, 4.0) |  1(0.1, 1.0, 4.0)
100 | (0.05, 2.0) | 4.0 | (0.05, 5.0) | (0.15, 1.2, 7.0)

Table 5:  Grid parameter values used for the plots in figures 6 to 8. These
values give the best average integration error $\varepsilon_{\text{rms}}$ in a
parameter scan. The labels exp, exp sqrt, and Gauss refer to the variable
transformations in table 1. The corresponding grids have the structure
$[0,\infty]_{(48)}$ for one parameter, $[0,z_{1},\infty]_{(16,32)}$ for two
parameters, and $[0,z_{1},z_{2},\infty]_{(16,16,16)}$ for three parameters.
$z_{1}$ and $z_{2}$ are given in units of $\operatorname{GeV}^{-1}$.

The relative integration error for the best parameter setting is shown in
figure 6 as a function of $q$. We see that for the TMDs with an exponential
falloff (toy and Yukawa TMD) the exp sqrt transform performs somewhat better
than the exp transform. In turn, for the Gauss TMD , the exp transform
performs much better than the Gauss transform. We note that in both cases, the
preferred transformation is the one that leads to a stronger decrease of the
integrand for $u\to 0$ (which corresponds to $z\to\infty$). We use this
preferred transformation in all subsequent comparisons.

(a)

(b)

(c)

Figure 6:  Relative integration error achieved with different variable
transforms on grids of the form $[0,z_{1},\infty]_{(16,32)}$. Details are
given in the text.

####  5.2.2 Number of subgrids

Let us now compare the accuracy of our algorithm with one, two, or three
subgrids for a comparable total number of points. As in the previous
subsection, we determine the best parameter settings for each grid by a
parameter scan; the resulting values are given in table 5. We note in passing
that for a single subgrid, the preferred values of $r_{\kappa}$ or
$r_{\lambda}$ are much larger than with two or more subgrids.

More importantly, we find that for our integrands and $q$ values, two subgrids
are strongly preferred, as seen in figures 7 and 8. We hence adopt this
setting for the subsequent investigations. We also tried settings with two
subgrids and a different partition $(n_{1},n_{2})$ for the number of grid
points, finding no significant improvement compared to the choice $(16,32)$.

(a)

(b)

(c)

Figure 7:  Relative integration error for grids of the form
$[0,\infty]_{(48)}$ (left) and $[0,z_{1},\infty]_{(16,32)}$ (right). The
variable transform for a given TMD is as in the right panels of figure 6.

(a)

(b)

(c)

Figure 8:  The same as figure 7, but for grids of the form
$[0,z_{1},z_{2},\infty]_{(16,16,16)}$ (left) and $[0,z_{1},\infty]_{(16,32)}$
(right).

####  5.2.3 Transformation parameters and number of grid points

Our method would be of limited practical use if its accuracy changed rapidly
with the choice of grid parameters, since this would require a detailed
parameter scan for each new integrand. Fortunately, this is not the case, as
we will show now.

Figures 9 to 11 show the dependence of the relative integration error on the
parameters $z_{1}$ and $r_{\kappa}$ or $r_{\lambda}$. They give the _worst_
accuracy for all $q_{i}$ values away from the zero crossing of $I(q)$, i.e.

$\displaystyle\varepsilon_{\text{max}}$ $\displaystyle=\max_{q_{i}\notin
R}(\varepsilon)$ $\displaystyle\text{ with }R$
$\displaystyle=(q_{0}-\Delta,q_{0}+\Delta)$ (5.20)

with $\Delta$ given in (5.19). We obtain very similar plots for
$\varepsilon_{\text{rms}}$ defined in (5.18). Note that
$\varepsilon_{\text{max}}$ is a more stringent measure for the performance of
a parameter setting.

We see in figures 9 and 10 that the best accuracy for $I(q)$ at $Q=2$ or
$20\operatorname{GeV}$ is obtained with $z_{1}$ between $0.05$ and
$0.15\operatorname{GeV}^{-1}$ and with $r_{\kappa}$ or $r_{\lambda}$ between
$2$ and $4$. For $Q=100\operatorname{GeV}$, the region of favourable
parameters is somewhat shifted, with $z_{1}$ between $0.05$ and
$0.10\operatorname{GeV}^{-1}$ and with $r_{\kappa}$ or $r_{\lambda}$ between
$4$ and $6$. The parameter region of good accuracy for $K(q)$ is much wider in
$r_{\kappa}$ or $r_{\lambda}$ and includes the good region for $I(q)$.
Corresponding plots for the toy TMD are similar to those for the Yukawa TMD
and not shown here.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9:  The relative integration error $\varepsilon_{\text{max}}$ defined
in (5.20), evaluated for the Yukawa TMD as a function of the grid parameters
with our default grid structure $[0,z_{1},\infty]_{(16,32)}$. The accuracy for
$I(q)$ is shown on the left and the one for $K(q)$ on the right.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10:  As figure 9, but for the Gauss TMD.

With fewer grid points, namely $(12,24)$ instead of $(16,32)$, we find that
the favourable parameter regions become smaller, as is shown in the left
panels of figure 11. Still, a relative accuracy of $10^{-3}$ or better can be
achieved in the same parameter region that is preferred for $(16,32)$ points.
In the right panels of figure 11 we see that with even more points, namely
$(24,48)$, the favourable parameter region becomes much wider, although for
$Q=100$ there is still a clear preference for low $z_{1}$. The corresponding
plots for the Yukawa and the toy TMD look quite similar and are not shown
here.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11:  The relative integration error $\varepsilon_{\text{max}}$ of
$I(q)$ for the Gauss TMD, evaluated with fewer or more grid points that in the
left panels of figure 10.

###  5.3 Estimating the integration error

To determine whether a particular grid setup in our method is satisfactory,
one needs a practical method to estimate the integration accuracy. A
straightforward way to achieve this is to re-compute the integral with twice
as many points on each subgrid, keeping all other settings the same. We have
evaluated the corresponding error estimates for the grids
$[0,z_{1},\infty]_{(16,32)}$ and $[0,\infty]_{(48)}$, both for the the best
and for the worst parameter settings identified in our parameter scans. A
selection of our results in shown in figure 12.

(a)

(b)

Figure 12:  Actual and estimated relative integration accuracy, as described
in the text. The left panels are for the best setting of the corresponding
parameter scan and the right panels for the worst one. The worst setting in
the upper row is $(z_{1},r_{\kappa})=(0.03,1.0)$, the worst setting in the
lower row is $r_{\kappa}=0.5$. The best settings are given in table 5.

In general we find that the error estimates so obtained are extremely close to
the actual integration errors. The only notable exception to this statement is
shown in the top left panel, where we see that for most $q$ values the error
is _underestimated_ by doubling the number of points. Since this problem
appears in a region where the actual relative error is well below $10^{-6}$,
we consider it to be of minor importance in practice. We therefore do not
pursue this issue further here.

We also observe that sometimes the error estimate fails only at specific $q$
values, as seen in the bottom left panel. This happens when the difference
between approximate and exact integrals has a zero crossing as a function of
$q$. As long as one monitors the error estimate also in the vicinity of such
points, this should not pose any problem.

###  5.4 Comparison with the method of Kang et al.

Let us finally compare our method with the algorithm of Kang et al. [8], which
implements the Ogata method [5] with a particular optimisation of the
parameter $h$ that controls the density of integration nodes. We use the C++
implementation of the algorithm that links to the GSL library, see
https://github.com/UCLA-TMD/Ogata for details. The optimisation algorithm
requires an estimate for the value of $x$ at which the function $x^{2}\mskip
1.5muW(x/q)$ takes its maximum; we verified that our results are stable with
respect to variations of this parameter. The other input parameter to the
algorithm is the number $N$ of function calls to $W(z)$.

For our method, we take the preferred parameter transform for each TMD,
together with a grid $[0,z_{1},\infty]_{(n_{1},n_{2})}$ with our default grid
parameters (5.1). For a total number $N$ of grid points, and thus of calls to
the function $W(z)$, we take $n_{1}=\operatorname{floor}(N/3)+1$ and
$n_{2}=N-\operatorname{floor}(N/3)$ points on the respective subgrids. We
recall that the grid parameters in (5.1) were determined for
$(n_{1},n_{2})=(16,32)$, i.e. for $N=47$.

When comparing the accuracy of the two algorithms at the same value of $N$,
one should bear in mind that with our method the computation of $W(z)$ at $N$
values of $z$ allows one to evaluate the Fourier-Bessel transform for a broad
range of $q$ values, which is not the case for the Ogata method. We
nevertheless use $N$ as a measure of the computational cost in both cases.

The following plots show the accuracy achieved for the toy and Gauss TMDs. The
corresponding accuracy for the Yukawa TMD is typically between these two
cases, both for our method and the one of [8]. We omit this in our plots for
the sake of clarity.

We see in figure 13 that for low $q$ the algorithm of [8] works very well. It
often performs better than our method for small $N$. The situation is very
different for $q=20\operatorname{GeV}$ and $q=100\operatorname{GeV}$, as shown
in figure 14. With $N=60$ points our method achieves relative errors below
$10^{-4}$, except for the toy TMD at $q=Q=100\operatorname{GeV}$, where $I(q)$
is close to its zero crossing. By contrast, the algorithm of [8] requires a
much higher number of nodes for $q=20\operatorname{GeV}$ (well above $100$ for
the toy TMD), and at $q=100\operatorname{GeV}$, no satisfactory result is
obtained for the toy TMD even with $320$ nodes.

We observe that the intermediate parameter $h_{u}$ in the optimisation
procedure of Kang et al. grows linearly with $q$ until it reaches its maximum
allowed value $h_{\text{max}}=2$. As discussed in section 2 of [8], the value
of $h_{u}$ used in the algorithm is increasingly far away from its optimum as
$q$ increases, and a reliable result can then only be achieved by increasing
the number $N$ of function calls.

In summary, we find that both methods require more grid points as $q$
increases, but that these requirements are much less severe for our method.

(a)

(b)

(c)

Figure 13:  Comparison of our method with the algorithm of Kang et al. at
$q=0.01\operatorname{GeV}$ (left) and $q=2\operatorname{GeV}$ (right) for
different values of $Q$. The relative integration error is evaluated for $N$
calls to the function $W(z)$, with $N$ varying in steps of $5$ or $10$.
Continuous curves are shown as a guide for the eye.

(a)

(b)

(c)

Figure 14:  As figure 13 but for $q=20\operatorname{GeV}$ (left) and
$q=100\operatorname{GeV}$ (right). $N$ varies in steps of $20$ for $N\geq
100$.

##  6 Summary

We have presented a method to compute Fourier-Bessel transforms, $\int
dz\,J_{\nu}(qz)\,\tilde{f}(z)$, with $z$ ranging from $0$ to $\infty$ or being
restricted to a smaller interval. The method works for general index $\nu\geq
0$ of the Bessel function. The function $\tilde{f}(z)$ may be non-analytic at
$z=0$ (details are given in section 3.1), and for $z\to\infty$ it should
decrease or tend to a constant.

Our method is a modification of an algorithm by Levin [13]. For $\nu\geq 1$
the integration problem is transformed to the system of ODEs given in (3), and
integration by parts provides access to Bessel functions of order $\nu<1$ as
specified in (3). The system of ODEs is solved by collocation, with Chebyshev
polynomials in a transformed variable $u(z)$ as basis functions. The resulting
matrix equation is handled with standard linear algebra methods. As a
complement, we use Clenshaw-Curtis quadrature in $z$ intervals where $qz$ does
not exceed the first zero of $J_{\nu}(x)$. This avoids most of the cases in
which this matrix equation is poorly conditioned, so that its solution would
require a numerically expensive singular value decomposition. The details of
our algorithm are laid out in section 3.3 and implemented in the C++ library
BestLime, which can be downloaded from [16].

The method requires a number of presets to be decided by the user. The first
one is a variable transform $u(z)$ that maps the integration interval in $z$
onto a finite interval in $u$ and should be adapted to the large-$z$ behaviour
of $\tilde{f}(z)$. The $z$-grid on which $\tilde{f}(z)$ must be evaluated is
then specified by the number of points, with the possibility of dividing the
total integration range into several subintervals. As we have shown, such a
subdivision is beneficial for both the accuracy and the computational cost of
integration. This set of choices is admittedly more involved than for other
methods (the Ogata method [5] for instance requires a step parameter $h$ and
the total number $N$ of function evaluations). On the other hand, we found
that with about $50$ to $80$ grid points, a high integration accuracy can be
reached for a wide range of grid parameters and over a wide range of $q$, so
that a detailed tuning of the grid settings should in general not be
necessary. Settings that we found to be satisfactory are given in section 4.2
for a wide class of functions, and in equations (5.1) to (5.1) for integrands
appearing in TMD cross sections.

Our method is most efficient for computing Bessel-weighted integrals of
several functions $\tilde{f}(z)$ at the same value of $q$. Changing the value
of $q$ requires somewhat more computation time for linear algebra
manipulations, but does _not_ necessitate a new evaluation of the function(s)
$\tilde{f}(z)$, since the grid in $z$ does not depend on $q$. This presents an
important advantage if the function $\tilde{f}(z)$ is expensive to compute. In
an ongoing separate project, we found that computing two-fold Fourier-Bessel
transforms $\int dz_{1}\int
dz_{2}\,J_{\nu_{1}}(q_{1}z_{1})\,J_{\nu_{2}}(q_{2}\mskip
1.5muz_{2})\,\tilde{f}(z_{1},z_{2})$ can be done accurately and efficiently
with the implementation of our method in BestLime.

## Acknowledgements

It is a pleasure to thank Frank Tackmann and Alexey Vladimirov for useful
discussions. The implementation of our method in BestLime is heavily drawing
on the ChiliPDF library [18, 19], which is under development. We gratefully
acknowledge the contributions of our collaborators Florian Fabry, Peter Plößl,
Riccardo Nagar, and Frank Tackmann to that project. Special thanks go to
Riccardo Nagar and Frank Tackmann for their consent to releasing the present
code before a full release of ChiliPDF.

This work is in part supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) – grant number 409651613 (Research Unit FOR 2926)
and grant number 491245950. O.G. is supported by the German Academic
Scholarship Foundation.

##  Appendix A Useful relations for Bessel and Lommel functions

In this appendix, we collect results about Bessel and Lommel functions that
are used in the main body of our work.

###  A.1 Bessel functions

The properties of Bessel functions are well documented in the literature. The
following relations can for instance be found in [25, Chapter 10].

##### Relations.

Functions of different order are related by

$\displaystyle Z_{\nu}(x)$
$\displaystyle=\frac{x}{2\nu}\,\Bigl{[}Z_{\nu-1}(x)+Z_{\nu+1}(x)\Bigr{]}\,,$
$\displaystyle\frac{d}{dx}\,Z_{\nu}(x)$
$\displaystyle=\frac{1}{2}\,\Bigl{[}Z_{\nu-1}(x)-Z_{\nu+1}(x)\Bigr{]}\,,$
(A.1)

which implies

$\displaystyle\frac{d}{dx}\mskip 1.5mu\Bigl{[}x^{\nu}Z_{\nu}(x)\Bigr{]}$
$\displaystyle=x^{\nu}\,Z_{\nu-1}(x)\,,$ (A.2)

where $Z_{\nu}$ stands for $J_{\nu}$ or $Y_{\nu}$, i.e. for the Bessel
functions of the first or the second kind. The two sets of functions are
related by

$\displaystyle Y_{\nu}(x)$
$\displaystyle=\frac{\cos(\nu\pi)\,J_{\nu}(x)-J_{-\nu}(x)}{\sin(\nu\pi)}$
(A.3)

for non-integer $\nu$. For $n\in\mathbb{Z}$, one can take the limit $\nu\to n$
of the r.h.s. to define $Y_{n}$.

##### Behaviour at large arguments.

The leading asymptotic behaviour for $x\to\infty$ of $J_{\nu}(x)$ and
$Y_{\nu}(x)$ consists of oscillations and a global decrease like $x^{-1/2}$,
namely

$\displaystyle J_{\nu}(x)$ $\displaystyle\sim\sqrt{\frac{2}{\pi
x}}\,\cos\Bigl{(}x-\nu\mskip 1.5mu\frac{\pi}{2}-\frac{\pi}{4}\Bigr{)}\,,$
$\displaystyle Y_{\nu}(x)$ $\displaystyle\sim\sqrt{\frac{2}{\pi
x}}\,\sin\Bigl{(}x-\nu\mskip 1.5mu\frac{\pi}{2}-\frac{\pi}{4}\Bigr{)}\,.$
(A.4)

##### Series representations.

For $\nu\geq 0$ one has

$\displaystyle x^{\nu}\mskip 1.5muJ_{\nu}(x)$
$\displaystyle=x^{2\nu}\,\sigma(x^{2})\,,$ $\displaystyle x^{\nu}\,Y_{\nu}(x)$
$\displaystyle=\begin{cases}\sigma(x^{2})+x^{2\nu}\mskip
1.5mu\ln(x)\,\sigma(x^{2})&\text{ for }\nu\in\mathbb{Z}\,,\\\
\sigma(x^{2})+x^{2\nu}\,\sigma(x^{2})&\text{ for
}\nu\notin\mathbb{Z}\,,\end{cases}$ (A.5)

where we use our generic notation $\sigma(x^{2})$ for a power series in
$x^{2}$, explained after equation (2.20).

###  A.2 Lommel functions

The literature about Lommel functions is much more sparse than for Bessel
functions. This holds especially for the Lommel functions of the second kind,
which are relevant for the discussion in section 2.3. We therefore provide
some more detail here and in particular derive a representation in terms of
power series. Unless explicitly mentioned, the following relations can be
found in [25, Chapter 11.9] or in [26, Chapter 8.57].

The Lommel functions of the first and second kind, $s_{\mu,\nu}(x)$ and
$S_{\mu,\nu}(x)$, solve the second-order differential equation

$\displaystyle
x^{2}\,w^{\prime\prime}(x)+x\,w^{\prime}(x)+\bigl{[}x^{2}-\nu^{2}\bigr{]}\,w(x)$
$\displaystyle=x^{\mu+1}\,.$ (A.6)

The Lommel functions of first kind are given in terms or a power series as

$\displaystyle s_{\mu,\nu}(x)$ $\displaystyle=x^{\mu+1}\,\sigma(x^{2})$ (A.7)

and oscillate like the Bessel functions for large $x$. For the Lommel
functions of the second kind one has

$\displaystyle S_{\mu,\nu}(x)$
$\displaystyle=s_{\mu,\nu}(x)+2^{\mu-1}\,\Gamma\biggl{(}\frac{\mu+\nu+1}{2}\biggr{)}\,\Gamma\biggl{(}\frac{\mu-\nu+1}{2}\biggr{)}$
$\displaystyle\quad\times\biggl{[}\sin\Bigl{(}\mskip-1.5mu(\mu-\nu)\frac{\pi}{2}\,\Bigr{)}\,J_{\nu}(x)-\cos\Bigl{(}\mskip-1.5mu(\mu-\nu)\frac{\pi}{2}\,\Bigr{)}\,Y_{\nu}(x)\biggr{]}\,.$
(A.8)

If $\mu+\nu$ or $\mu-\nu$ is a negative odd integer, several coefficients of
the power series in (A.7) become infinite, as do the $\Gamma$ functions in
(A.2). The function $S_{\mu,\nu}$ can be defined by a limiting procedure in
these cases, as shown in [27, Chapter 3.4] or [28, Chapter 10.73].

##### Relations.

The Lommel functions of the second kind satisfy

$\displaystyle S_{\mu,\nu}(x)$ $\displaystyle=S_{\mu,-\nu}(x)$ (A.9)

and

$\displaystyle S_{\mu,\nu}(x)$
$\displaystyle=\frac{x}{2\nu}\,\Bigl{[}(\mu+\nu-1)\,S_{\mu-1,\nu-1}(x)-(\mu-\nu-1)\,S_{\mu-1,\nu+1}(x)\Bigr{]}\,,$
$\displaystyle\frac{d}{dx}\,S_{\mu,\nu}(x)$
$\displaystyle=\frac{1}{2}\,\Bigl{[}(\mu+\nu-1)\,S_{\mu-1,\nu-1}(x)+(\mu-\nu-1)\,S_{\mu-1,\nu+1}(x)\Bigr{]}\,,$
(A.10)

which implies

$\displaystyle\frac{d}{dx}\mskip 1.5mu\Bigl{[}x^{\nu}S_{\mu,\nu}(x)\Bigr{]}$
$\displaystyle=(\mu+\nu-1)\,x^{\nu}\,S_{\mu-1,\nu-1}(x)\,.$ (A.11)

Notice the close analogy of these relations with (A.1) and (A.2). The identity

$\displaystyle S_{\mu+2,\nu}(x)$
$\displaystyle=x^{\mu+1}-(\mu+\nu+1)(\mu-\nu+1)\,S_{\mu,\nu}(x)$ (A.12)

implies that

$\displaystyle S_{\nu+1,\nu}(x)$ $\displaystyle=x^{\nu}\,,$ $\displaystyle
S_{-\nu+1,\nu}(x)$ $\displaystyle=x^{-\nu}\,,$ (A.13)

and by recursion

$\displaystyle x^{-\nu}\,S_{\nu+1+2k,\nu}(x)$
$\displaystyle=\sum_{j=0}^{k}c_{\nu,\mskip 1.5muk,\mskip 1.5muj}\,x^{2j}\,,$
$\displaystyle x^{\nu}\,S_{-\nu+1+2k,\nu}(x)$
$\displaystyle=\sum_{j=0}^{k}c_{-\nu,\mskip 1.5muk,\mskip 1.5muj}\,x^{2j}$
(A.14)

for $k\in\mathbb{N}$, with coefficients $c_{\nu,\mskip 1.5muk,\mskip 1.5muj}$
whose explicit form we do not need.

##### Behaviour at large arguments.

For $x\to\infty$, the Lommel functions of the second kind have an asymptotic
expansion

$\displaystyle x^{\nu}\mskip 1.5muS_{\mu,\nu}(x)$ $\displaystyle\sim
x^{\mu+\nu-1}\,\sigma(x^{-2})\,,$ (A.15)

where the first term in the series $\sigma$ is equal to $1$. This expansion is
also valid for the special cases mentioned below equation (A.2).

##### Series representation.

As in the case of Bessel functions, we assume $\nu\geq 0$. Combining (A.1),
(A.7), and (A.2), one obtains a series representation

$\displaystyle x^{\nu}\mskip 1.5muS_{\mu,\nu}(x)$
$\displaystyle=\begin{cases}\sigma(x^{2})+x^{\mu+\nu+1}\,\sigma(x^{2})+x^{2\nu}\mskip
1.5mu\ln(x)\,\sigma(x^{2})&\text{ for }\nu\in\mathbb{Z}\,,\\\
\sigma(x^{2})+x^{\mu+\nu+1}\,\sigma(x^{2})+x^{2\nu}\,\sigma(x^{2})&\text{ for
}\nu\notin\mathbb{Z}\,,\end{cases}$ (A.16)

provided that neither $\mu+\nu$ nor $\mu-\nu$ is a negative odd integer. For
these exceptional values of $\mu$ and $\nu$, we need separate derivations. We
limit ourselves to values satisfying

$\displaystyle\mu+\nu>-1\,,$ (A.17)

which is sufficient for the discussion in section 2.3. This leaves us with the
possibility that $\mu-\nu$ is a negative odd integer. We then can distinguish
two cases:

  * •

$\nu\in\mathbb{Z}$. Then $\mu+\nu$ is a positive odd integer, because we
impose (A.17). We can therefore use the second representation in (A.14) and
find that $x^{\nu}\,S_{\mu,\nu}(x)$ is a polynomial in $x^{2}$.

  * •

$\nu\notin\mathbb{Z}$. In this case we can use a representation derived in
[27, Chapter 3.4] and also in [28, Chapter 10.73]. It has the form

$\displaystyle x^{\nu}\mskip 1.5muS_{\nu-1,\nu}(x)$
$\displaystyle=x^{2\nu}\mskip 1.5mu\sigma(x^{2})+2^{\nu-1}\mskip
1.5mu\Gamma(\nu)\biggl{[}x^{\nu}J_{\nu}(x)\,\ln\mskip-1.5mu\frac{x}{2}-\frac{\pi}{2}\,x^{\nu}\mskip
1.5muY_{\nu}(x)\biggr{]}$
$\displaystyle=\sigma(x^{2})+x^{2\nu}\,\sigma(x^{2})+x^{2\nu}\mskip
1.5mu\ln(x)\,\sigma(x^{2})$ (A.18)

and is valid for all $\nu>0$. Starting from this, we can repeatedly apply the
recursion relation (A.12) to obtain

$\displaystyle x^{\nu}\mskip 1.5muS_{\mu=\nu-1-2k,\mskip 1.5mu\nu}(x)$
$\displaystyle=d_{\mu,\nu}\,x^{\nu}\,S_{\nu-1,\nu}(x)+x^{\mu+\nu+1}\mskip
1.5mu\sigma(x^{2})\,,$ $\displaystyle k=1,2,3,\ldots,$ (A.19)

where $d_{\mu,\nu}$ is some number and the series multiplying $x^{\mu+\nu+1}$
is in fact a polynomial.

Putting everything together, we find that one can always write

$\displaystyle x^{\nu}\mskip 1.5muS_{\mu,\nu}(x)$
$\displaystyle=\sigma(x^{2})+x^{\mu+\nu+1}\,\sigma(x^{2})+x^{2\nu}\,\sigma(x^{2})+x^{2\nu}\mskip
1.5mu\ln(x)\,\sigma(x^{2})$ (A.20)

for $\mu+\nu>-1$ and $\nu\geq 0$. If $\mu+\nu$ is a positive odd integer, then
the r.h.s. reduces to a polynomial in $x^{2}$. Otherwise, the term with
$\ln(x)$ is present if $\nu$ is an integer or if $\nu-\mu$ is a negative odd
integer.

##  Appendix B Integrands with a Gaussian decrease in $z$

In section 2.3.2 we solved the system (2.16), (2.17) of differential equation
for functions $f_{1}$ or $f_{2}$ with an exponential decrease at large
arguments. Let us do the same for the case of functions with a Gaussian
decrease. We treat the functions $f_{1}$ and $f_{2}$ in turn, since the
corresponding solutions are slightly different.

We start with $f_{2}=0$ and

$\displaystyle f_{1}(z)$ $\displaystyle=c_{1}\mskip
1.5muz^{\mu+\nu}\,\exp\bigl{(}-\lambda^{2}z^{2}-\kappa z\mskip
1.5mu\bigr{)}\,,$ (B.1)

where $\lambda>0$, whilst $\kappa$ may be positive, negative, or zero. With
the ansatz

$\displaystyle h_{i}(z,q)$
$\displaystyle=\hat{h}_{i}(z,q)\,\exp(-\lambda^{2}z^{2}-\kappa z)$
$\displaystyle(i=1,3)$ (B.2)

the system (2.16), (2.17) of ODEs becomes

$\displaystyle\begin{pmatrix}c_{1}\mskip 1.5muz^{\mu+\nu}\\\ 0\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}2\lambda^{2}z\hat{h}_{1}\\\ q\mskip
1.5mu\hat{h}_{1}+2\lambda^{2}z^{2}\mskip
1.5mu\hat{h}_{3}\end{pmatrix}-\kappa\begin{pmatrix}\hat{h}_{1}\\\
z\hat{h}_{3}\end{pmatrix}+\begin{pmatrix}d\mskip 1.5mu\hat{h}_{1}/dz+q\mskip
1.5muz\hat{h}_{3}\\\ z\mskip 1.5mud\mskip 1.5mu\hat{h}_{3}/dz-2\nu\mskip
1.5mu\hat{h}_{3}\end{pmatrix}\,,$ (B.3)

which can in turn be solved with the ansatz

$\displaystyle\hat{h}_{1}(z,q)$
$\displaystyle=z^{\mu+\nu-1}\,\sum_{k=0}^{\infty}b_{1,k}(q)\mskip
1.5muz^{-k}\,,$ $\displaystyle\hat{h}_{3}(z,q)$
$\displaystyle=z^{\mu+\nu-3}\,\sum_{k=0}^{\infty}b_{3,k}(q)\mskip
1.5muz^{-k}\,.$ (B.4)

We find that for each $k$ in the sums, the first term on r.h.s. of (B.3) is
leading at large $z$, whilst the second term is down by $z^{-1}$ and the third
term is down by $z^{-2}$. One can therefore solve for the coefficients
recursively, starting with

$\displaystyle b_{1,0}$ $\displaystyle=-c_{1}\mskip 1.5mu/(2\lambda^{2})\,,$
$\displaystyle b_{3,0}$ $\displaystyle=-q\mskip 1.5mub_{1,0}\mskip
1.5mu/(2\lambda^{2})\,.$ (B.5)

We now turn to the case where $f_{1}=0$ and

$\displaystyle f_{2}(z)$ $\displaystyle=c_{2}\mskip
1.5muz^{\mu+\nu}\,\exp\bigl{(}-\lambda^{2}z^{2}-\kappa z\mskip
1.5mu\bigr{)}\,.$ (B.6)

If we make again the ansatz (B.2), the system (2.16), (2.17) of ODEs takes the
form

$\displaystyle\begin{pmatrix}0\\\ c_{2}\mskip 1.5muz^{\mu+\nu}\end{pmatrix}$
$\displaystyle=-\begin{pmatrix}2\lambda^{2}z\hat{h}_{1}-q\mskip
1.5muz\hat{h}_{3}\\\ 2\lambda^{2}z^{2}\mskip
1.5mu\hat{h}_{3}\end{pmatrix}-\kappa\begin{pmatrix}\hat{h}_{1}\\\
z\hat{h}_{3}\end{pmatrix}+\begin{pmatrix}d\mskip 1.5mu\hat{h}_{1}/dz\\\
z\mskip 1.5mud\mskip 1.5mu\hat{h}_{3}/dz-2\nu\mskip 1.5mu\hat{h}_{3}-q\mskip
1.5mu\hat{h}_{1}\end{pmatrix}\,.$ (B.7)

We now set

$\displaystyle\hat{h}_{1}(z,q)$
$\displaystyle=z^{\mu+\nu-2}\,\sum_{k=0}^{\infty}b_{1,k}(q)\mskip
1.5muz^{-k}\,,$ $\displaystyle\hat{h}_{3}(z,q)$
$\displaystyle=z^{\mu+\nu-2}\,\sum_{k=0}^{\infty}b_{3,k}(q)\mskip
1.5muz^{-k}\,,$ (B.8)

such that the three terms on the r.h.s. of (B.7) are leading, subleading in
$z^{-1}$, and subleading in $z^{-2}$, respectively. One can then solve for the
coefficients iteratively, starting with

$\displaystyle b_{3,0}$ $\displaystyle=-c_{2}\mskip 1.5mu/(2\lambda^{2})\,,$
$\displaystyle b_{1,0}$ $\displaystyle=q\mskip 1.5mub_{3,0}\mskip
1.5mu/(2\lambda^{2})\,.$ (B.9)

##  Appendix C Chebyshev interpolation

In this appendix, we collect formulae for Chebyshev interpolation that are
relevant in our method. We use the same notation as in [18, 19], where
references or derivations for the following results can be found.

We start with functions defined on the interval $[-1,1]$, as is customary in
the discussion of Chebyshev polynomials. The generalisation to any other
finite interval is trivial and given at the end of this appendix.

##### Functions with support $t\in[-1,1]$.

The Chebyshev polynomials of the first kind, $T_{k}(t)$, can be defined by

$\displaystyle T_{k}(\cos\theta)$ $\displaystyle=\cos(k\mskip 1.5mu\theta)\,,$
(C.1)

where $k\geq 0$ is an integer. Using the multiple-angle formula for the
cosine, one readily verifies that $T_{k}(t)$ is a polynomial in $t$ of order
$k$. The set of all $T_{k}(t)$ is orthonormal with respect to the integration
measure $1/\sqrt{1-t^{2}}$.

For a given integer $N$, the _Chebyshev points_ are given by

$\displaystyle t_{j}$
$\displaystyle=\cos\theta_{j}\,,\qquad\theta_{j}=\frac{j\pi}{N}\qquad\text{with}\quad
j=0,\ldots,N\,.$ (C.2)

They form a descending series from $t_{0}=1$ to $t_{N}=-1$ and obey the
reflection property $t_{N-j}=-t_{j}$. The polynomial $T_{N}(t)$ assumes its
maxima $+1$ and minima $-1$ at the Chebyshev points, as is easily seen from
(C.1). We note that the density of Chebyshev points increases from the centre
toward the end points of the interval $[-1,1]$. This feature is crucial to
avoid Runge’s phenomenon for equispaced interpolation grids.

The _Chebyshev interpolant_ $p_{N}(t)$ of a function $f(t)$ is the unique
polynomial that satisfies

$\displaystyle p_{N}(t_{j})$ $\displaystyle=f(t_{j})\qquad\text{for}\quad
j=0,\ldots,N$ (C.3)

at the Chebyshev points $t_{j}$. It can be written in the form

$\displaystyle p_{N}(t)$ $\displaystyle=\sum_{k=0}^{N}\,\beta_{k}\mskip
1.5muc_{k}\mskip 1.5muT_{k}(t)$ (C.4)

with coefficients

$\displaystyle c_{k}$
$\displaystyle=\frac{2}{N}\sum_{j=0}^{N}\,\beta_{j}\mskip
1.5muT_{k}(t_{j})\,f(t_{j})\,,$ (C.5)

where $\beta_{0}=\beta_{N}=1/2$ and $\beta_{j}=1$ otherwise.

For functions $f(t)$ that are sufficiently regular in the interval $[-1,1]$
(including the end points), one can show that the difference $|f(t)-p_{N}(t)|$
uniformly tends to zero on the full interval in the limit $N\to\infty$. A
precise statement of this convergence theorem can be found in [18, section 2].
We note that a function is _not_ sufficiently regular for the purpose of this
theorem if its first derivative diverges at an endpoint of the interval. This
explains our insistence on a finite first derivative of $h_{1}$ and $h_{3}$ in
the discussion of section 2.3.

##### Differentiation.

Given the Chebyshev interpolant $p_{N}(t)$ for a function $f(t)$, one can
approximate the derivative $f^{\prime}(t)=df/dt$ of the function by the
derivative $dp_{N}/dt$ of its interpolant. Note that in general $df/dt$ is not
equal to $dp_{N}/dt$ at the Chebyshev points; this is obvious because one
cannot compute the exact values of $f^{\prime}(t_{j})$ from the discrete set
of function values $f(t_{j})$. At the Chebyshev points, the derivative of the
interpolating polynomial can be computed as

$\displaystyle\frac{d}{dt}\mskip 1.5mup_{N}(t_{j})$
$\displaystyle=\sum_{k=0}^{N}\,D_{jk}\,f(t_{k})$ (C.6)

with a differentiation matrix given by

$\displaystyle D_{00}$ $\displaystyle=-D_{NN}=\frac{2N^{2}+1}{6}\,,$
$\displaystyle D_{jj}$
$\displaystyle=-\dfrac{\cos\theta_{j}}{2\sin^{2}\theta_{j}}$
$\displaystyle\text{for}\quad j\neq 0,N,$ $\displaystyle D_{jk}$
$\displaystyle=\dfrac{\beta_{k}}{\beta_{j}}\,\dfrac{(-1)^{j+k}}{t_{j}-t_{k}}$
$\displaystyle\text{for}\quad j\neq k\,.$ (C.7)

Note that the matrix multiplication (C.6) maps a vector $f(t_{k})=\text{\it
const}$ onto the zero vector, as it must be because the derivative of a
constant function is zero.

##### Integration.

Approximating a function by its Chebyshev interpolant, one can derive an
integration rule

$\displaystyle\int_{-1}^{1}dt\,f(t)$
$\displaystyle\approx\sum_{\genfrac{}{}{0.0pt}{}{k=0}{\text{even}}}^{N}\,\frac{2\mskip
1.5mu\beta_{k}\mskip 1.5muc_{k}}{1-k^{2}}=\sum_{j=0}^{N}\,w_{j}\mskip
1.5muf(t_{j})$ (C.8)

with weights

$\displaystyle w_{j}$
$\displaystyle=\frac{4\beta_{j}}{N}\,\sum_{\genfrac{}{}{0.0pt}{}{k=0}{\text{even}}}^{N}\,\beta_{k}\,\frac{\cos(k\mskip
1.5mu\theta_{j})}{1-k^{2}}\,.$ (C.9)

This is known as _Clenshaw-Curtis quadrature_. A detailed discussion of its
accuracy, as well as a comparison with Gauss quadrature, is given in [17,
chapter 19] and in [29].

##### Functions with support $u\in[u_{a},u_{b}]$.

It is easy to generalise the preceding results to an arbitrary finite interval
in the variable $u$ by introducing the linear transformation

$\displaystyle u(t)$
$\displaystyle=\frac{u_{a}-u_{b}}{2}\,t+\frac{u_{a}+u_{b}}{2}\,.$ (C.10)

With the Chebyshev points

$\displaystyle u_{j}$
$\displaystyle=u(t_{j})=\frac{u_{a}-u_{b}}{2}\,\cos\frac{j\pi}{N}+\frac{u_{a}+u_{b}}{2}$
(C.11)

in the new variable, we obtain the following rules for differentiation and for
integration:

$\displaystyle\frac{d}{du}\mskip 1.5muf(u_{j})$
$\displaystyle\approx\sum_{k=0}^{N}\,D^{u}_{jk}\,f(u_{k})\,,$ $\displaystyle
D^{u}_{jk}$ $\displaystyle=\frac{2}{u_{a}-u_{b}}\,D_{jk}$ (C.12)

with $D_{jk}$ from (C) and

$\displaystyle\int_{u_{a}}^{u_{b}}du\,f(u)$
$\displaystyle\approx\sum_{j=0}^{N}\,w^{u}_{j}\mskip 1.5muf(u_{j})\,,$
$\displaystyle w^{u}_{j}$ $\displaystyle=\frac{u_{b}-u_{a}}{2}\,w_{j}$ (C.13)

with $w_{j}$ from (C.9).

## References

  * [1] R. Boussarie et al., _TMD Handbook_ , arXiv:2304.03302.
  * [2] M. A. Ebert, J. K. L. Michel, I. W. Stewart and Z. Sun, _Disentangling long and short distances in momentum-space TMDs_ , _JHEP_ 07 (2022) 129 [arXiv:2201.07237].
  * [3] M. G. A. Buffing, M. Diehl and T. Kasemets, _Transverse momentum in double parton scattering: factorisation, evolution and matching_ , _JHEP_ 01 (2018) 044 [arXiv:1708.03528].
  * [4] P. Cal, R. von Kuk, M. A. Lim and F. J. Tackmann, _The $q_{T}$ spectrum for Higgs production via heavy quark annihilation at N3LL’ + aN3LO_, arXiv:2306.16458.
  * [5] H. Ogata, _A numerical integration formula based on the Bessel functions_ , _Publ. Res. Inst. Math. Sci._ 41 (2005) 949–970.
  * [6] I. Scimemi and A. Vladimirov, _Analysis of vector boson production within TMD factorization_ , _Eur. Phys. J. C_ 78 (2018) 89 [arXiv:1706.01473].
  * [7] J. Isaacson, Y. Fu and C. P. Yuan, _Improving ResBos for the precision needs of the LHC_ , arXiv:2311.09916.
  * [8] Z.-B. Kang, A. Prokudin, N. Sato and J. Terry, _Efficient Fourier Transforms for Transverse Momentum Dependent Distributions_ , _Comput. Phys. Commun._ 258 (2021) 107611 [arXiv:1906.05949].
  * [9] N. Baddour and U. Chouinard, _Theory and operational rules for the discrete Hankel transform_ , _J. Opt. Soc. Am. A_ 32 (2015) 611.
  * [10] N. Baddour, _The Discrete Hankel Transform_ , in _Fourier Transforms – Century of Digitalization and Increasing Expectations_.  IntechOpen, 2019.  https://doi.org/10.5772/intechopen.84399.
  * [11] S. Venkat, J. Arrington, G. A. Miller and X. Zhan, _Realistic Transverse Images of the Proton Charge and Magnetic Densities_ , _Phys. Rev. C_ 83 (2011) 015203 [arXiv:1010.3629].
  * [12] D. Levin, _Procedures for computing one- and two-dimensional integrals of functions with rapid irregular oscillations_ , _Math. Comp._ 38 (1982) 531.
  * [13] D. Levin, _Fast integration of rapidly oscillatory functions_ , _Journal of Computational and Applied Mathematics_ 67 (1996) 95.
  * [14] J. Li, X. Wang, T. Wang and S. Xiao, _An improved Levin quadrature method for highly oscillatory integrals_ , _Applied Numerical Mathematics_ 60 (2010) 833.
  * [15] J. Li, X. Wang, T. Wang, S. Xiao and M. Zhu, _On an improved-Levin oscillatory quadrature method_ , _Journal of Mathematical Analysis and Applications_ 380 (2011) 467.
  * [16] M. Diehl and O. Grocholski, _BestLime: a C++ library for computing Fourier- Bessel transforms with Levin’s integration method_ , 2024.  10.5281/zenodo.11113673.
  * [17] L. N. Trefethen, _Approximation Theory and Approximation Practice_. Society for Industrial and Applied Mathematics, 2012.  https://people.maths.ox.ac.uk/trefethen/ATAP.
  * [18] M. Diehl, R. Nagar and F. J. Tackmann, _ChiliPDF: Chebyshev interpolation for parton distributions_ , _Eur. Phys. J. C_ 82 (2022) 257 [arXiv:2112.09703].
  * [19] M. Diehl, R. Nagar, P. Plößl and F. J. Tackmann, _Evolution and interpolation of double parton distributions using Chebyshev grids_ , _Eur. Phys. J. C_ 83 (2023) 536 [arXiv:2305.04845].
  * [20] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, _Numerical Recipes_. Cambridge University Press, Cambridge, MA, 2007.  https://numerical.recipes/book.html.
  * [21] Gaël Guennebaud and Benoît Jacob and others, _Eigen v3_ , 2010.  https://eigen.tuxfamily.org.
  * [22] _GNU Scientific Library Documentation_.  https://www.gnu.org/software/gsl/doc/html/index.html.
  * [23] A. Bacchetta, F. Delcarro, C. Pisano, M. Radici and A. Signori, _Extraction of partonic transverse momentum distributions from semi-inclusive deep-inelastic scattering, Drell-Yan and Z-boson production_ , _JHEP_ 06 (2017) 081 [arXiv:1703.10157].  [Erratum: JHEP 06, 051 (2019)].
  * [24] R. Piessens, E. de Doncker-Kapenga, C. Überhuber and D. Kahaner, _QUADPACK. A subroutine package for automatic integration_. Springer, 1983, 10.1007/978-3-642-61786-7.
  * [25] F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier et al., _NIST Digital Library of Mathematical Functions_.  Release 1.1.12 of 2023-12-15. https://dlmf.nist.gov.
  * [26] I. S. Gradshteyn, I. M. Ryzhik, D. Zwillinger and V. Moll, _Table of integrals, series, and products_. Academic Press, Amsterdam, 2014, 10.1016/C2010-0-64839-5.
  * [27] Y. L. Luke, _Integrals of Bessel Functions_. McGraw-Hill, New York, 1962.  Reprinted in 2014.
  * [28] G. N. Watson, _A Treatise on the Theory of Bessel Functions_. Cambridge University Press, Cambridge, UK, 2nd ed., 1944.  Reprinted in 1995.
  * [29] L. N. Trefethen, _Is Gauss Quadrature Better than Clenshaw-Curtis?_ , _SIAM Review_ 50 (2008) 67.