Impact of the heavyquark matching scales in PDF fits
Eur. Phys. J. C
Impact of the heavyquark matching scales in PDF fits
The xFitter Developers' Team:
V. Bertone 1 3
D. Britzger 0
S. Camarda 7
A. CooperSarkar 6
A. Geiser 0
F. Giuli 6
A. Glazov 0
E. Godat 5
A. Kusina 4 9
A. Luszczak 8
F. Lyonnet 5
F. Olness 5
R. Placˇakyte˙ 2
V. Radescu 0 7
I. Schienbein 4
O. Zenaiev 0
0 DESY Hamburg , Notkestraße 85, 22609 Hamburg , Germany
1 Nikhef Theory Group Science Park 105 , 1098 XG Amsterdam , The Netherlands
2 Institut für Theoretische Physik, Universität Hamburg , Luruper Chaussee 149, 22761 Hamburg , Germany
3 Department of Physics and Astronomy, VU University , 1081 HV Amsterdam , The Netherlands
4 Laboratoire de Physique Subatomique et de Cosmologie, Université Grenoble Alpes , CNRS/IN2P3, 53 avenue des Martyrs, 38026 Grenoble , France
5 SMU Physics , Box 0175, Dallas, TX 752750175 , USA
6 University of Oxford , 1 Keble Road, Oxford OX1 3NP , UK
7 CERN , 1211 Geneva 23 , Switzerland
8 T.Kosciuszko Cracow University of Technology , 30084 Kraków , Poland
9 Institute of Nuclear Physics, Polish Academy of Sciences , ul. Radzikowskiego 152, 31342 Kraków , Poland
We investigate the impact of displaced heavyquark matching scales in a global fit. The heavyquark matching scale μm determines at which energy scale μ the QCD theory transitions from NF to NF + 1 in the variable flavor number scheme (VFNS) for the evolution of the parton distribution functions (PDFs) and strong coupling αS(μ). We study the variation of the matching scales, and their impact on a global PDF fit of the combined HERA data. As the choice of the matching scale μm effectively is a choice of scheme, this represents a theoretical uncertainty; ideally, we would like to see minimal dependence on this parameter. For the transition across the charm quark (from NF = 3 to 4), we find a large μm = μc dependence of the global fit χ 2 at NLO, but this is significantly reduced at NNLO. For the transition across the bottom quark (from NF = 4 to 5), we have a reduced μm = μb dependence of the χ 2 at both NLO and NNLO as compared to the charm. This feature is now implemented in xFitter 2.0.0, an open source QCD fit framework. 1 Introduction . . . . . . . . . . . . . . . . . . . . . . 2 Variable flavor number scheme (VFNS) . . . . . . . . 2.1 The matching scale μm . . . . . . . . . . . . . . 2.2 Historical choice of μm = mc,b,t . . . . . . . . . 2.3 Smooth matching across flavor thresholds . . . .

Contents
1 Introduction
The global analysis of PDFs has progressed significantly in
recent years. On the experimental front, there is data
ranging from the fixedtarget regime at low energy, on to HERA
and the LHC at very high energies. On the theoretical front,
the analysis can be performed not only at NLO, but now
at NNLO. To capitalize on these advances, it is essential
to include a proper treatment of the heavy quarks to enable
highprecision phenomenological analysis of measurements.
The variable flavor number scheme (VFNS) allows us to
deal with the heavyquark mass scale across the full
kinematic range by varying the number of active flavors (NF ) in
the DGLAP QCD evolution [
1–11
]. At low energy scales, the
DGLAP evolution only involves NF light flavors, and there
is no PDF for the heavy quark. At high energy, the
heavyquark PDF is included in the DGLAP evolution so that there
are now NF + 1 active flavors. To combine the above NF and
NF + 1 subschemes into a single VFNS, we must define an
energy scale μm where we match these together; this will be
the scale where we introduce the heavyquark PDF.
Historically, the matching scale μm was taken to be the
heavyquark mass m H . At the matching scale, the PDFs and
αS(μ) for NF + 1 are defined in terms of the NF quantities
by the following boundary conditions:
(1)
(2)
fi(NF +1)(x , μm )
=
×
j
1 +
α(NF +1)(μm ) = αS(NF )(μm )
S
j
Mi ⊗ f j(NF )(x , μm )
∞
n
n=1 k=0
cn k α(NF )(μm )
S
n lnk μ2m
m2H
.
The matching matrix Mij and coefficients cn k can be
perturbatively computed.1
The new xFitter 2.0.0 program2 links to the APFEL
code [
18
] which has implemented generalized matching
conditions that enable the switch from NF to NF + 1 at an
arbitrary matching scale μm . This allows us to introduce the
heavyquark PDF at any scale—not just at μm = m H ; this
flexibility provides a number of advantages. For example,
as the matching scale moves to higher scales, the theory at
the lower scales effectively becomes a fixed flavor number
scheme (FFNS); yet we still retain a VFNS at the higher
scales.
The choice of the matching scale μm , like the choice of
VFNS or FFNS, amounts to a theoretical scheme choice.
As such, the variation of μm represents a source of
theoretical uncertainty. The variable matching scale implemented in
xFitter provides a new incisive tool to study the impact of
these choices across a broad kinematic region. Additionally,
as we move from NLO to NNLO calculations, new features
are encountered, and these compel us to reexamine some of
the foundational elements used to construct this theoretical
framework.
Reconsidering the historical choice μm = m H is of
particular relevance for heavyquark initiated processes at the
LHC. In this context, the benefits of the FFNS close to
the threshold region and of the VFNS at higher scales are
often simultaneously needed to describe the data. Therefore,
a careful choice of the matching scales could help formulate
1 The perturbative coefficients of Mij at NLO are available in Refs. [
12,
13
], and at NNLO in Ref. [14]. m H is the mass of the NF + 1 flavor
quark. For αS(μ), the cn k coefficients are available in the Particle Data
Group review of Quantum Chromodynamics [
15
].
2 Information on the xFitter program can be found at http://www.
xFitter.org, and in Refs. [
16,17
].
a matching prescription between FFNS and VFNS able to
achieve this goal in a very simple fashion [19].
This study will examine the combined HERA data set and
evaluate the impact of the matching scale on the features of
the fit of PDFs. In Sect. 2, we review the key elements of the
VFNS used in this study. Section 3, shows the impact of the
matching scale μm on the PDFs. In Sect. 4, we perform a fit
of the combined HERA data sets at both NLO and NNLO,
and investigate the effect of the matching scale μm . Section 5
presents an example of how the μm flexibility can be used
as a tool to evaluate a recent suggestion for a NF dependent
PDF. Section 6 summarizes the general characteristics and
conclusions of this study.
2 Variable flavor number scheme (VFNS)
Here we will outline the key concepts of the heavyquark
VFNS which are relevant for this investigation.
2.1 The matching scale μm
A generalized formulation of the VFNS factorization is based
on the Collins–Wilczek–Zee (CWZ) renormalization scheme
which involves a sequence of subschemes parameterized by
the number of active quark flavors (NF ) [
20,21
]. For each
subscheme, the NF (active) flavors are renormalized using
the MS scheme while the heavy (inactive) flavors are
renormalized using zeromomentum subtraction. This ensures that
to all orders in perturbation theory (i) the results are gauge
invariant, (ii) the results for the active NF flavors match the
standard MS results, and (iii) the heavy (inactive) flavors
manifestly decouple.3 Specifically, both the DGLAP
evolution kernels for the NF active PDFs and the renormalization
group equation for α(NF )(μ) are pure MS.
S
To connect the separate NF subschemes into a single
scheme that spans the full kinematic range, we must choose
a matching scale μm which will relate the subschemes. This
is where we define the PDFs and αS of the NF + 1 scheme
in terms of the NF scheme, cf. Eqs. (1) and (2). A schematic
representation of this is displayed in Fig. 1.
For example, at scales μc < μ < μb the scheme has
NF = 4 active flavors {u, d, s, c} with 4flavor PDFs and
αS(4)(μ); the bottom quark is not treated as a parton and
f (4)(x , μ) = 0.
b
At the scale μ = μb, we can compute the 5flavor PDFs
and αS(5)(μ) in terms of the 4flavor quantities; the boundary
conditions are nontrivial and the PDFs and αS(μ) are not
3 For the CWZ scheme with NF (active) flavors and an arbitrary number
of heavy (inactive) flavors, the evolution of the PDFs and α(NF )(μ) will
S
involve only the active NF flavors; the inactive heavy flavors can be
ignored.
necessarily continuous. This scheme has NF = 5 active
flavors {u, d, s, c, b}, and the bottom quark is included in the
DGLAP evolution.
2.2 Historical choice of μm = mc,b,t
Historically, the matching scale μm was commonly taken to
be exactly equal to the mass of the heavy quark, μm = mc,b,t ;
this was a convenient choice for a number of reasons.
For example, the generic NLO matching condition for the
PDFs at the NF = 4 to NF = 5 transition is [
22
]:
fi(5)(x , μb)
where c0ij and c1ij are perturbatively calculable coefficient
functions. Note that the righthand side uses 4flavor PDFs
and αS, while the lefthand side uses 5flavors.
The choice μb = mb will cause the logarithms to
vanish, and this greatly simplifies the matching relations.
Additionally, at NLO in the MS scheme the constant term c0ij in
the matching equation coincidentally vanishes [
14
]. The net
result is that, for μb = mb, the PDFs will be continuous (but
not differentiable) at NLO. This is historically why μm was
set to mc,b,t .
However, at NNLO and beyond the situation is more
complex; in particular, the higherorder terms corresponding to
c0ij will be nonzero, and the matching of both the PDFs and
αS(μ) will be discontinuous. Consequently, the freedom to
arbitrarily choose the matching scale μm (and decide where
to place the discontinuities) will have a number of
advantages, as the next subsection will demonstrate.
2.3 Smooth matching across flavor thresholds
To gauge the impact of the contributions of the heavyquark
PDFs in a process independent manner, we can compare the
DGLAP evolved heavyquark PDF fb(x , μ) with a
perturbatively computed quantity: fb(x , μ). At NLO, fb(x , μ) takes
a gluon PDF and convolutes it with a perturbative (DGLAP)
splitting g → bb¯ [
23,24
]; this can be thought of as a
“perturbatively” computed bottom PDF. The result at NLO is
αS P
fb(x , μ) = 2π
g→bb¯ ⊗ fg ln
μ2
m2
b
The difference between fb(x , μ) and fb(x , μ) is due to the
higherorder terms which are resummed by the heavyquark
DGLAP evolution.4
To better understand these quantities, we compute DIS
bottom production at NLO in a 5flavor VFNS, and find the
cross section to be [
3
]:
σVFNS = σb→b ⊗ fb(x , μ) − fb(x , μ) +σg→b ⊗ fg(x , μ) .
∼ σFFNS
(5)
Here, σb→b ⊗ fb is the LO term, and σb→b ⊗ fb is the
subtraction (SUB) term. The unsubtracted NLO term σg→b ⊗ fg
corresponds (approximately) to a FFNS calculation. Here,
σb→b is proportional to a delta function which makes the
convolution trivial.
Thus, the combination ( fb − fb) represents
(approximately) the difference between a VFNS and FFNS result.5
These quantities are displayed in Fig. 2. In the region μ ∼
mb, fb(x , μ) and fb(x , μ) match precisely; it is this
cancellation which (at NLO) ensures physical quantities will have
a smooth transition across the flavor threshold.
4 In Eq. (4), fb(x, μ) includes the single splitting (g → bb¯); in
contrast, the DGLAP evolution of fb(x, μ) sums an infinite tower of
splittings. Note, we have used the NNPDF30_lo_as_118_nf_6 PDFs to
precisely match the order of the splitting kernels in the NLO calculation.
5 The above correspondences are only approximate as the VFNS and
FFNS also differ in αS and the PDFs.
At larger μ scales, fb(x , μ) and fb(x , μ) begin to
diverge; this indicates that the resummed heavyquark
logarithms are becoming sizable. The details clearly depend
on the specific x values. For large x (x ∼ 0.1) we find
fb(x , μ) > fb(x , μ), while for small x (x ∼ 0.001) the
result is fb(x , μ) < fb(x , μ); finally, for intermediate x
(x ∼ 0.01) the two terms nearly balance even for sizable μ
scales.
While the QCD theory ensures proper matching, this is
not so easy to implement in a general numeric calculation for
all observables, especially for complex observables
involving multiple numeric integrations. In particular, the
cancellation of Fig. 2 requires that the quark masses mc,b,t , the
strong coupling αS, and the order of the PDF evolution are
exactly matched in (i) the DGLAP evolution that generates
the PDFs, (ii) the partonic cross sections that are
convoluted with the PDFs, and (iii) the fragmentation function
(if used).
In practice, there are almost always slight differences. A
typical analysis might use a variety of PDFs from different
PDF groups, together with a selection of fragmentation
functions; each of these will be generated with a specific set of
quark masses and αS values which are most likely different.
Thus, it is essentially inevitable that the cancellations
exhibited in Fig. 2 will be spoiled leading to spurious contributions
which can be substantive.
Instead of setting the matching scale at the heavyquark
mass μm = mc,b,t , xFitter provides the flexibility to
delay the matching scale μm to a few multiples of the
heavyquark mass; this will avoid the need for the delicate
cancellation in the μm ∼ mc,b,t region, and the results will be
numerically more stable.
As an extreme example, one could imagine delaying the
matching scale to infinity (μm → ∞) which would amount
to a FFNS; here, the disadvantage is that the FFNS does
not include the resummation of the higherorder heavyquark
logs which have been demonstrated to improve the fit to the
data [
25
]. Using the new flexibility of the xFitter program,
it is possible to investigate the tradeoffs between a large and
small value for the matching scale μm .
A separate example is present in the transverse momentum
( pT ) distributions for heavyquark production ( pp → bb¯)
using the (general mass) GMVFNS [
26,27
]. If we compute
this in an NF = 5 flavor scheme, the contribution from the
bb¯ → bb¯ subprocess with an exchanged t channel gluon
will be singular at pT = 0. For a scale choice of the transverse
mass μ = pT2 + m2b (a common choice), the singularity can
be cured by either a different scale choice, or by delaying the
switch to the 5flavor scheme to a higher scale, e.g., μb ∼
2mb.
2.4 Discontinuities
At NNLO both the PDFs and αS(μ) will necessarily have
discontinuities when matching between the NF to NF + 1
flavor schemes as specified by Eqs. (1) and (2). If we are
analyzing a highprecision experiment and arbitrarily impose a
matching at the quark masses μm = mc,b,t , this may well
introduce discontinuities within the kinematic range of some
precision data. While it is true that these discontinuities
simply reflect the theoretical uncertainties, it is disconcerting to
insert them in the middle of a precision data set.
The ability to vary the matching scale μm provides us
with the option to shift the location of these
discontinuities for a particular analysis. For example, to analyze the
highprecision charm production HERA data, we
necessarily are working in the region of the bottom mass scale
(∼ 4.5 GeV). Both the PDFs and αS(μ) will be
discontinuous at the matching scale which transitions between the
NF = 4 and NF = 5 schemes. If the matching scale is
chosen in the region μm ∼ mb, these discontinuity will
appear in the region of the data. Instead, we can shift the
matching μm to a higher scale (for example, set μm to 2mb
or 3mb) and thus analyze the charm production data in a
consistent NF = 4 flavor framework. Yet, we still retain
the transition to NF = 5 flavors so that processes such as
LHC data at high scales are computed including the bottom
PDF.
Fig. 3 We display the bquark
PDF x f (5)(x , μ) for different
b
choices of the matching scales
μm = {mb/2, mb, 2mb}
(indicated by the vertical lines)
computed at NLO (a) and
NNLO (b)
0.3
)
,
x
f(b 0.2
x
0.1
0
−0.1
−0.2
O( s) matching conditions
b = mb / 2
b = mb
b = 2 mb
(a)
0.3
)
,
x
f(b 0.2
x
3 The matching scale µm
Having sketched the characteristics of a flexible matching
scale μm , we will examine the specific boundary condition
and the impact on the global fit of the PDFs.
3.1 Impact of matching on the PDFs
where L = ln(μ2b/m2b). The superscripts {4, 5} identify the
number of active flavors NF . The gluon and the light quarks
also have matching conditions analogous to Eq. (6).
As already mentioned, if we choose to match at μb = mb
then L = 0 and f (5)(x , μm ) will start from zero at μb = mb.
This coincidentalb zero (c0ij = 0) is the historic reason why
most NLO analyses perform the matching at μb = mb; if
both the c0ij and the c1ij L terms can be ignored, then the PDFs
are continuous (but not differentiable) across the matching
scale.8
6 A first study of the impact of moving the bottom matching scale with
respect to the bottom mass was already done in Ref. [
28
] in the context
of bb H production at the LHC using a matched scheme. The approach
developed in this study was more recently applied to the 13 TeV LHC
in Ref. [
29
].
7 At NNLO, the bottomquark matching condition also receives
contributions from the light quarks as well as gluons; this has been included
in the calculation.
8 While the VFNS framework is compatible with an intrinsic charm
or bottom PDF, we do not introduce these into the current study. For
additional details, see Refs. [
30–33
].
At NNLO this is no longer the case; the NNLO
constant term at O(αS2 ) does not vanish and the PDFs will have
a discontinuity regardless of the choice of matching scale.
Although the difference is subtle, the (red) curve for μb = mb
does start exactly from zero for the NLO calculation (Fig. 3a),
while for the NNLO calculation (Fig. 3b) it starts from a small
nonzero value.
As we vary the matching μb in the vicinity of mb, the
i j
sign of fb(5)(x , μb) is controlled by the log term (c1 L ). For
μb < mb this combination will render fb(5)(x , μb) negative,
and this will be compensated (in the sum rule for example)
by a positive shift in the 5flavor gluon. Thus, QCD ensures
that both momentum and number sum rules are satisfied to
the appropriate order.
Comparing different fb(5)(x , μ) curves computed with the
NLO matching conditions (Fig. 3a) at large μ scales, there are
obvious differences in the curves. This reflects the difference
between the single log contribution (c1ij L ) computed by the
matching condition of Eq. (6) and the resummed
contributions computed by the DGLAP evolution equation.
Specifically, the NLO matching includes the αS L contribution, but
it is missing α2 L 2 and higher terms; this is what gives rise to
S
the differences of Fig. 3a. Obviously, the α2 L 2 contributions
S
can be important.
Comparing the different fb(5)(x , μ) curves computed with
the NNLO matching conditions (Fig. 3b) at large μ scales, the
differences in the curves are greatly reduced compared to the
NLO case. The NNLO result includes both the αS L and the
α2 L 2 contributions, but it is missing α3 L 3 and higher orders.
S S
Clearly the inclusion of the α2 L 2 contributions dramatically
S
reduces the effect of the different choices of the μm matching
scale.
Finally, we wish to emphasize that ultimately the choice
of μm amounts to a choice of scheme. In the limit that
perturbation theory is computed to all orders, the infinite tower
of logarithms resummed by the DGLAP evolution equations
(in the NF + 1flavor scheme) will be explicitly summed
(b)
)
,xQ0.04
b(2
F
in the matching conditions (in the NF flavor scheme). In a
practical sense, while the differences at NLO are substantive,
at NNLO the residual differences at large μ scale are much
smaller. This reduced sensitivity on the choice of μm
provides increased flexibility and precision in our fits, as will be
illustrated in the following sections.
3.2 Impact of matching on F2b(x , Q)
Having examined the PDFs in the previous section we now
turn to a physical observable, F2b(x , Q).
Figure 4a shows the NLO result for F2b(x , Q) which will
receive contributions from the LO process (γ b → b) as well
as the NLO (γ g → bb¯) process. For μ < μb, fb(5)(x , μ) = 0
and only the gluon initiated process contributes. For μ μb,
the bottom PDF turns on (cf. Fig. 3), and the heavyquark
initiated process now contributes. Because the PDFs, αS(μ),
and mb are all carefully matched in this calculation, the
cancellation outlined in Sect. 2.3 ensures that the prediction for
the physical observable is relatively smooth in this region.
Figure 4b shows the NNLO result for F2b(x , Q). As with
the PDF matching of Fig. 3b, the additional NNLO
contributions significantly reduce the impact of the different
matching scales so that the prediction for F2b(x , Q) is now very
insensitive to μb.
The above smooth transition of F2b(x , Q) from the NF =
4 to the NF = 5 scheme holds even though the PDFs and
αS(μ) have discontinuities. Because we have used
consistent choices for {mb, fi(NF ), αS}, the cancellation of Sect. 2.3
applies, and the effect of any discontinuities in the physical
observable will be of higher order. Conversely, a mismatch
in {mb, fi(NF ), αS} would spoil this cancellation and result
in unphysical large contributions when f (5)(x , μ) is
introb
duced. This is precisely the case where shifting the matching
scale μb to a higher scale such as 2mb or 3mb would help
avoid these problems.
It is interesting to note that, as we compute even higher
orders, the discontinuities in the PDFs and αS(μ) will
persist at lower order; but any discontinuities in the physical
observables will systematically decrease order by order.
4 The PDF fits
4.1 xFitter, APFEL, and data sets
To study the effects of varying the matching scales for the
charm and bottom quark we will perform a series of fits to
various data sets. Since we are varying the matching scales
in the vicinity of mc and mb, we want data that constrain the
PDFs in this region. For this purpose, we include the very
precise combined HERA data sets as these provide strong
constraints in the region μ ∼ mc,b, and they also extend up to
higher scales [
34–37
]. In particular, the HERA measurement
of the charm and bottom cross sections are included as they
are sensitive to the choice of μc and μb.
These fits are performed with the xFitter program
using the APFEL evolution code [
18,38,39
]. The DIS
calculations use the FONLLB scheme for the NLO
calculations, and the FONLLC scheme for the NNLO calculations;
these are both O(αS2) prescriptions, and the details are
specified in Ref. [
6
]. We use mc = 1.45 GeV, mb = 4.5 GeV,
αS(MZ ) = 0.118 for both the NLO and the NNLO
calculations. The fit is performed using pole masses, but the
formalism can be used equally well with the MS definition of
the heavyquark masses [
40
]. For the PDFs, we use a
HERAPDF 14parameter functional form with initial QCD
evolution scale Q2
0 = 1.0 GeV2 and strangeness fraction fs = 0.4;
the other QCD fit settings and constraints are similar to the
analysis of Ref. [
40
].
The minimization of the χ 2 is performed using MINUIT
[
41
]. The correlations between data points caused by
sysTable 1 The χ 2 values at NLO
for individual data sets for a
selection of the charm matching
scales μc. The contribution of
the charm data contained in the
“Correlated χ 2” and in the “Log
penalty χ 2” terms is indicated
separately in the parentheses
2 mc
61/47
2.8/12
12/17
44/39
47/42
228/159
70/70
433/377
217/204
224/254
91 (12.5)
− 0.7 (− 0.4)
1430/1207
2 mc
tematic uncertainties are taken into account in the
“Correlated χ 2” contribution. A “Log penalty χ 2” arises from the
likelihood transition to χ 2 when the scaling of the errors is
applied [
16, 42
].
The full sets of data are listed in Tables 1, 2, 3 and 4,
and the reference for each data set is cited in Table 1. The
combined inclusive HERA data (HERA1+2) from Ref. [
34
]
includes both neutral current (NC) and charged current (CC)
results for electrons (em) and positrons (ep) at a variety of
energies. The charm cross sections from Ref. [
36
] include
the combined H1ZEUS results. The bottom cross sections
from ZEUS are presented in Ref. [
37
] and those from H1 in
Ref. [
35
].
4.2 Impact of matching on the fits: charm
region, and this is reflected in the results of Figs. 5, 6 and
Tables 1, 2.
Figure 5 displays the results for varying the charm quark
matching scale μc both for the NLO and NNLO
calculations.9 Comparing the NLO and NNLO cases, the NLO result
ranges over ∼ 100 units in χ 2, while the NNLO varies over
∼ 25 units of χ 2. This difference in the χ 2 variation reflects
the effects of the higherorder terms; it is reassuring to see
that the μc dependence decreases at higher orders.
At NLO, the matching conditions pick up the contribution
of only the single log term L (Eq. (6)), while at NNLO we
pick up both the L and the L 2 terms. In contrast, the DGLAP
evolved charm PDF resums the above, as well as an infinite
tower of logs: n∞=1 kn=0 αSn L k .
Examining the NLO analysis of Fig. 5a, we find that at
low scales, the χ 2 increases with increasing μc scale. While
The charm cross section data are expected to be
sensitive to the treatment of the charm PDF in the threshold
9 For these scans we hold the bottom matching fixed at μb = mb and
keep μc < mb so the ordering of the mass thresholds is not inverted.
3 mc
46/47
3.2/12
12/17
44/39
52/42
219/159
65/70
410/377
221/204
216/254
86 (0.8)
+ 4.2 (− 0.1)
1379/1207
3 mb
Fig. 6 χ2 vs. the charm
matching scale μc at a NLO and
b NNLO for only the H1ZEUS
combined charm production
data; note, this includes the
correlated χ2 contribution from
Tables 1 and 2
NLO
HERA charm data only
NNLO
HERA charm data only
our plot extends slightly below the charm mass, it is not
obvious if there is actually a minimum in μc. It is problematic to
compute with μc values much lower than mc as αS becomes
large and the charm PDF negative. Thus, the optimal
computational range for μc appears to be in the region of mc.
Focusing on the charm data alone as shown in Fig. 6a,
the situation is not so clear; the χ 2 increases with
increasing μc, but again there does not appear to be a minimum
at low μc values. Moving to large μc, the χ 2 values
initially increase, but then decrease as μc approaches mb. As
we want to maintain the ordering μc < μb, we cannot go to
larger scales unless we increase μb. While this is allowed,
it is more complex to explore the twodimensional {μc, μb}
parameter space; hence, we limit the present study to
variation of a single scale.
The χ 2 results for each individual data set is summarized
in Table 1. The data sets with the largest effects are (i) the
H1ZEUS combined charm cross section data, and (ii) the
very precise “HERA1+2NCep 920” set. The sensitivity of
the “HERA1+2NCep 920” set is due to a large number of
data points with small uncertainties.
Turning to the NNLO analysis of Fig. 5b and the results
of Table 2, a number of points are evident. Again, the two
data sets with the largest impact are the H1ZEUS combined
charm cross section data, and the “HERA1+2NCep 920” set.
In Fig. 5 the vertical lines indicate the bin boundaries for the
“HERA1+2NCep 920” data set.
Scanning in χ 2, discrete jumps are evident. As we vary the
matching scale, certain data bins move between the NF = 3
and NF = 4 schemes, shifting the χ 2 by one or two units
which is visible in Fig. 5b. These jumps reflect the underlying
theoretical uncertainty arising from the choice of NF .
In Fig. 5b the total NNLO variation of χ 2 is reduced
compared to the NLO case, and the minimum global χ 2 is now
in the region μc ∼ 2mc. Focusing on the charm data alone in
Fig. 6b, again it is not obvious if there is actually a minimum
in μc. Given the limitations of computing with μc mc,
the optimal computational range again appears to be in the
general region of mc.
While it may be tempting to try and optimize the matching
scale for each data set, recall that μm represents a choice of
scheme, and thus reflects an inherent theoretical uncertainty;
a specific choice of μm will not reduce this uncertainty.
This situation can also be found in complex global fits
where the final result may be a compromise of data sets
which are in tension; this is why a tolerance factor is often
introduced. This complexity is evident when examining the
details of Tables 1 and 2 which demonstrate the minimum
χ 2 for individual data sets is not simply correlated; this will
be discussed further in Sect. 4.4. An additional challenge of
analyzing the charm case is that μc can only vary over the
limited dynamic range between ∼ mc and μb. This will not
be an issue for the bottom quark (because mt mb), which
is considered in the following section.
4.3 Impact of matching on the fits: bottom
Figure 7 presents the results for varying the bottomquark
matching scale μb both for the NLO and NNLO calculations.
This figure highlights the ranges of χ 2; the NLO result ranges
over approximately ∼ 10 units in χ 2, and the NNLO varies
by about the same amount.
The reduced χ 2 variation as compared to the charm
case reflects, in part, the decrease in the strong coupling
αS(mb) < αS(mc), which also diminishes the higherorder
contributions. Figures 5 with 7 there is a χ 2 range of ∼ 100
vs. ∼ 10 for NLO, and ∼ 15 vs. ∼ 10 for NNLO.
Examining the NLO analysis of Fig. 7a, there is a slight
minimum for χ 2 in the region μb ∼ 2mb with relatively
flat behavior at larger μb scales. Correspondingly, there is a
similar behavior when we focus on only the bottom data of
Fig. 8a. The χ 2 results for each individual data set is
summarized in Table 3.
The data sets with the largest effects are (i) the very
precise “HERA1+2NCep 920” set, and (ii) the separate H1 and
ZEUS bottom cross section data. The H1 and ZEUS
bottom cross sections display some minimal χ 2 variation in the
region μb ∼ mb, but then is relatively flat out to very high
scales (μb ∼ 14mb). It is primarily the “HERA1+2NCep
920” set which drives the shape of the χ 2 curve in the
μb ∼ mb region. Compared to the charm results, the
interpretation of the bottom cross section data requires some care
as the number of data points is smaller, and the relative
uncertainty larger.
Turning to the NNLO analysis of Fig. 7b, the variation
of the χ 2 curve is within ∼ 8 units across the range of
the plot. The resolution of the vertical χ 2 scale
accentuates the discrete jumps as the data bins move between the
NF = 4 and NF = 5 schemes. The bin boundaries for the
“HERA1+2NCep 920” data set are indicated with vertical
lines.
Focusing on the bottom data alone as shown in Fig. 8b,
the χ 2 profile is flat within one unit across the plot range.
For both Figs. 7b and 8b, the χ 2 variation is within a
reasonable “tolerance” factor for the global fit; thus, the
matching scale μb can vary within this range with minimal impact
on the resulting fit.
The scale μb can extend up to larger scales, and Tables 3
and 4 display the results for 10mb and 14mb. The pattern
across the various data sets is consistent, and the overall χ 2
values rise slowly.
4.4 Comparisons
To facilitate comparisons of the NLO and NNLO results,
Fig. 9 displays the ratio χ 2/χ02 for charm (on the left) and
bottom (on the right) where χ02 is the value of the χ 2 at
μm = m H . Similarly, Fig. 10 displays the same ratio for only
the heavyquark data sets. By plotting χ 2/χ02, we can better
compare the fractional variation of χ 2 across the matching
scale values.
The motivation for the scaled plot of Figs. 9 and 10 is that
the overall χ 2 values are different; specifically, those of the
NNLO are greater than the NLO. This counter intuitive result
has been observed in other analyses [
34,43
], and it has been
suggested that this may be improved by resumming the
singular ln[1/x ] terms in the higherorder splitting kernels [44].
Here, we first make some observations specific to Figs. 9
and 10.
– At NLO for the case of charm, the optimal
computational scale for μc is in the general range μc ∼ mc for
both the inclusive data set (Fig. 9a) and the charm data
set (Fig. 10a). For lower scales (μc mc), αS(μ) is
The triangles (blue filled triangle) are NLO and the diamonds (red filled
lozenge) are NNLO
Fig. 10 The ratio (χ 2/χ02) of partial χ 2 values (charm/bottom data
only) from Figs. 6 and 8 as a function of the a charm and b bottom
matching scale μc,b in GeV. χ02 is the χ 2 value for μm equal to the
large and the charm PDFs are negative. For higher scales
(μc mc), χ 2/χ02 increases.
– At NLO for the case of bottom, the optimal scale for μb is
in the general range μb ∼ 2mb. For the inclusive data set
(Fig. 9b) the χ 2/χ02 variation is very mild (∼ 1%), while
for the bottom data set (Fig. 10b) the χ 2/χ02 variation is
larger (∼ 10%).
– At NNLO for the case of charm, the χ 2/χ02 variation is
reduced. For the inclusive data set (Fig. 9a) the χ 2/χ02
variation is very mild (∼ 2%), while for the charm data
set (Fig. 10a) the χ 2/χ02 variation is larger (∼ 10%).
There is no obvious optimal choice for the μc scale.
– At NNLO for the case of bottom, the χ 2/χ02 variation
is reduced and a matching scale choice in the region
μb ∼ mb appears to be optimal. For the inclusive data set
(Fig. 9b) the χ 2/χ02 variation is very mild (∼ 1%), while
for the bottom data set (Fig. 10b) the χ 2/χ02 variation is
slightly larger (∼ 5%).
While the detailed characteristics of the above fits will
depend on specifics of the analysis, there are two general
patterns which emerge: (i) the χ 2 variation of the NNLO
results are generally reduced compared to the NLO results,
and (ii) the relative χ 2 variation across the bottom transition
is reduced compared to the charm transition. For example,
quark mass. The triangles (blue filled triangle) are NLO and the
diamonds (red filled lozenge) are NNLO
although the global χ 2 can be modified by different choices of
data sets and weight factors, these general properties persist
for each individual data set of Tables 1, 2, 3 and 4; in fact, we
see that the bulk of the data sets are quite insensitive to the
details of the heavyquark matching scale. Additionally, there
are a variety of prescriptions for computing the heavy flavor
contributions; these primarily differ in how the higherorder
contributions are organized. As a cross check, we performed
a NLO fit using the FONNLA scheme; while the absolute
value of χ 2 differed, the above general properties persisted.
The net result is that we can now quantify the
theoretical uncertainty associated with the transition between
different NF subschemes. In practical applications, if we choose
μc ∼ mc, the impact of the NF = 3 to NF = 4 transition is
reduced as this is often below the minimum kinematic cuts of
the analysis (e.g. Q2min and Wm2 in ). Conversely, the NF = 4 to
NF = 5 transition is more likely to fall in the region of fitted
data; hence, it is useful to quantify the uncertainty associated
with the μb choice.
5 An example: NF dependent PDFs
The variable matching scale μm can be used as an incisive
tool to explore various aspects of the PDFs and global fits.
As an example, Ref. [
22
] introduced an NF dependent PDF
fi (x , μ, NF ) where NF is the active number of flavors in the
VFNS. This extension provides additional flexibility in the
region of the heavyquark thresholds; however, the
implementation of Ref. [
22
] only used a fixed matching scale of
μm = m H . Using xFitter we can improve on this concept by
generating PDFs with a variable μm scale. We illustrate this
below and provide example grids at xFitter.org.
The usual PDF can be generalized to include an NF
dependence [
22
]: fi (x , μ) → fi (x , μ, NF ). In this approach,
the many NF = {3, 4, 5, . . .} flavor schemes coexist, and
they can be selected by specifying the number of active
flavors NF along with the other arguments of the PDF. This
concept is represented pictorially in Fig. 11. All the NF sets of
PDFs are available above the matching scale μm . For
example, with an NF dependent PDF, one could simultaneously
analyze selected data sets with NF = 4 and others with
NF = 5 even if they overlap kinematically; the user has the
flexibility (and responsibility) to select NF .
Note in Fig. 11 that the various NF grids are not
individual fits but are related analytically via the flavor threshold
matching conditions. Operationally, they are generated from
an initial PDF fi (x , μ0, NF = 3) and αS(μ0) at the starting
scale μ0. The NF = 3 grid is generated by evolving from
μ0 to μmax. The NF = 4 grid is then generated by matching
at μc (which may or may not equal mc), and evolving up
to scale μmax. The NF = 5 and NF = 6 grids are
generated in a similar manner.10 This process ensures that all the
PDFs fi (x , μ, NF ) are analytically related to the PDF and
αS boundary conditions at μ0.
To provide an explicit illustration of the above, we have
generated a set of PDF grids with a variety of matching scales
(μb) for the matching between the schemes with NF = 4 and
10 Note the NF = {3, 4, 5, 6} grids are stored in separate LHAPDF
data files; they can be combined into an effective NF dependent PDF
as illustrated in Refs. [
22,45
].
NF = 5 active flavors: μb = {1, 3, 5, 10, ∞}×mb. We focus
on μb as this is the flavor transition most likely to fall within
a particular data set. For the initial PDF we use the NNLO
bottom fit with μb = 1 mb of Table 4, and we evolve at
NNLO. The PDFs are fixed such that they all match at the
initial evolution scale μ0 = 1.0 GeV with the same value of
αS(μ0) = 0.467464.
This is illustrated in Fig. 12 where we display the
bottom quark and gluon PDFs as a function of μ in GeV. As
we evolve up in μ, we explicitly see the transition from
NF = 4 to NF = 5 flavors at each respective μb threshold.
For these particular kinematic values, the discontinuity of the
bottom PDF is positive while that of the gluon is negative;
this ensures the momentum sum rule is satisfied.
Furthermore, we observe the spread in the bottom PDF at large μ
is broader than that of Fig. 3. In Fig. 12, while the values of
αS all coincide at μ0, the evolution across the different μb
thresholds result in different αS values at large μ scales. This
is in contrast to Fig. 3 where the values of αS all coincide
at the large scale μ = MZ . Additionally, note that the
illustration in Fig. 3 is based on the NNPDF3.0 PDF set while
Fig. 12 is based on our fit from Table 4.
Because the NF = 4 and NF = 5 grids are available
concurrently, we can choose to analyze the HERA data in
an NF = 4 flavor scheme for arbitrarily large scales, but
simultaneously allow LHC data to be analyzed in a NF =
5 flavor scheme throughout the full kinematic region even
down to low scales.
In this illustration, the PDFs revert to NF = 4 below μb;
however, this is not required. For example the NF = 5 PDFs
could be evolved backwards from μb to provide values at
scales μ < μb. Both APFEL [
18
] and QCDNUM [
46,47
]
have this capability.11
For bottom at NNLO using the results from Table 4 for
the inclusive data set, we observe the μb variation is
minimal. Thus, a choice in the range μb ∼ [mb, 5mb] yields
a Δχ 2 ≤ (1457 − 1453) ∼ 4 units out of ∼ 1450. This
minimal χ 2 dependence means we can shift the μb matching
scale if, for example, we want to avoid a NF flavor
transition in a specific kinematic region. While these results should
be checked with additional data sets, the insensitivity to μb,
especially at NNLO, is an important result as the ability to
displace the NF = 4 and NF = 5 transition can be beneficial
when this threshold comes in the middle of a data set.
Combined with the variable heavyquark threshold, the
NF dependent PDFs provide additional flexibility to analyze
multiple data sets in the optimal theoretical context.
11 However, it is generally advisable not to backwards evolve too far
in μ as this can become unstable [
48,49
].
6 Conclusions
In this study we have examined the impact of the heavy flavor
matching scales μm on a PDF fit to the combined HERA data
set.
The choice of μm allows us to avoid delicate cancellations
in the region μm ∼ m H as illustrated in Fig. 2. Additionally,
the discontinuities associated with the NF = 4 to NF = 5
transition can be shifted so that these discontinuities do not
lie in the middle of a specific data set.
Using xFitter and APFEL to study the μm dependence
of a global PDF fit to the HERA data, we can extract the
following general features. For the charm matching scale, μc ,
there is a large variation of χ 2 at NLO, but this is significantly
reduced at NNLO. In contrast, for the bottom matching scale,
μb, there is a relatively small variation of χ 2 at both NLO
and NNLO.
These observations can be useful when performing fits.
While charm has a larger χ 2 variation (especially at NLO),
the charm quark mass mc ∼ 1.45 GeV lies in a region which
is generally excluded by cuts in Q2 and/or W 2.
On the contrary, the χ 2 variation for the bottom quark is
relatively small at both NLO and NNLO. Since the
bottomquark mass mb ∼ 4.5 GeV is in a region where there is
abundance of precision HERA data, this flexibility allows us
to shift the heavy flavor threshold (and the requisite
discontinuities) away from any particular data set. Functionally, this
means that we can analyze the HERA data using an NF = 4
flavor scheme up to relatively large μ scales, and then
perform the appropriate NNLO matching (with the associated
constants and log terms) so that we can analyze the highscale
LHC data in the NF = 5 or even NF = 6 scheme.
These variable heavy flavor matching scales μm allow us
to generalize the transition between a FFNS and a VFNS, and
as a function of μ in GeV. The vertical lines in the plots show the
transition from the NF = 4 to NF = 5 flavor scheme
provides a theoretical “laboratory” which can quantitatively
test proposed implementations. We demonstrated this with
the example of the NF dependent PDFs. Having the
quantitative results for the χ 2 variation of the μc,b scales, one
could systematically evaluate the impact of using different
matching scale choices for the fi (x , μ, NF ).
In conclusion, we find that the ability to vary the heavy
flavor matching scales μm , not only provides new insights
into the intricacies of QCD, but also has practical advantages
for PDF fits.
Acknowledgements The authors would like to thank M. Botje,
J. C. Collins, and J. Rojo for valuable discussions. V. B and A. G. are
particularly grateful to A. Mitov, A. Papanastasiou, and M. Ubiali for
many stimulating discussions on the role of the bottomquark threshold
for bottominitiated processes at the LHC. We acknowledge the
hospitality of CERN, DESY, and Fermilab where a portion of this work was
performed. We are grateful to the DESY IT department for their
support of the xFitter developers. This work was also partially supported
by the U.S. Department of Energy under Grant no. DESC0010129.
V. B. is supported by an European Research Council Starting Grant
“PDF4BSM”. A. L. is supported by the Polish Ministry under program
Mobility Plus, Grant no. 1320/MOB/IV/2015/0.
Open Access This article is distributed under the terms of the Creative
Commons Attribution 4.0 International License (http://creativecomm
ons.org/licenses/by/4.0/), which permits unrestricted use, distribution,
and reproduction in any medium, provided you give appropriate credit
to the original author(s) and the source, provide a link to the Creative
Commons license, and indicate if changes were made.
Funded by SCOAP3.
1. M.A.G. Aivazis , F.I. Olness , W.K. Tung , Phys. Rev. Lett . 65 , 2339 ( 1990 ). https://doi.org/10.1103/PhysRevLett.65.2339
2. M.A.G. Aivazis , F.I. Olness , W.K. Tung , Phys. Rev. D 50 , 3085 ( 1994 ). https://doi.org/10.1103/PhysRevD.50.3085
3. M.A.G. Aivazis , J.C. Collins , F.I. Olness , W.K. Tung , Phys. Rev. D 50 , 3102 ( 1994 ). https://doi.org/10.1103/PhysRevD.50.3102
4. R.S. Thorne , R.G. Roberts , Eur. Phys. J. C 19 , 339 ( 2001 ). https:// doi.org/10.1007/s100520100605
5. A.D. Martin , W.J. Stirling , R.S. Thorne , G. Watt, Eur. Phys. J. C 70 , 51 ( 2010 ). https://doi.org/10.1140/epjc/s1005201014628
6. S. Forte , E. Laenen, P. Nason , J. Rojo , Nucl. Phys. B 834 , 116 ( 2010 ). https://doi.org/10.1016/j.nuclphysb. 2010 . 03 .014
7. R.D. Ball , V. Bertone , F. Cerutti , L. Del Debbio , S. Forte , A. Guffanti , J.I. Latorre , J. Rojo , M. Ubiali , Nucl. Phys. B 849 , 296 ( 2011 ). https://doi.org/10.1016/j.nuclphysb. 2011 . 03 .021
8. S. Alekhin , J. Blümlein , S. Klein , S. Moch, Phys. Rev. D 81 , 014032 ( 2010 ). https://doi.org/10.1103/PhysRevD.81.014032
9. S. Alekhin , J. Blümlein , S. Moch, Phys. Rev. D 86 , 054009 ( 2012 ). https://doi.org/10.1103/PhysRevD.86.054009
10. S. Alekhin , J. Blümlein , S. Moch, Phys. Rev. D 89 ( 5 ), 054028 ( 2014 ). https://doi.org/10.1103/PhysRevD.89.054028
11. T. Stavreva , F.I. Olness , I. Schienbein, T. Ježo , A. Kusina , K. Kovarík , J.Y. Yu , Phys. Rev. D 85 , 114014 ( 2012 ). https://doi.org/ 10.1103/PhysRevD.85.114014
12. S. Qian, The CWZ Subtraction Scheme (A New Renormalization Prescription For QCD ) And Its Application . Ph.D. thesis , IIT, Chicago ( 1985 ). http://wwwlib.umi.com/dissertations/fullcit? p8517585
13. J.C. Collins, W.K. Tung , Nucl. Phys. B 278 , 934 ( 1986 ). https:// doi.org/10.1016/ 0550  3213 ( 86 ) 90425  6
14. M. Buza , Y. Matiounine , J. Smith , R. Migneron , W.L. van Neerven , Nucl. Phys. B 472 , 611 ( 1996 ). https://doi.org/10.1016/ 0550  3213 ( 96 ) 00228  3
15. C. Patrignani et al., Chin. Phys. C 40 ( 10 ), 100001 ( 2016 ). https:// doi.org/10.1088/ 1674 1137/40/10/100001
16. S. Alekhin et al., Eur. Phys. J. C 75 ( 7 ), 304 ( 2015 ). https://doi.org/ 10.1140/epjc/s100520153480z
17. O. Zenaiev, PoS DIS2016 , 033 ( 2016 )
18. V. Bertone , S. Carrazza , J. Rojo , Comput. Phys. Commun . 185 , 1647 ( 2014 ). https://doi.org/10.1016/j.cpc. 2014 . 03 .007
19. V. Bertone , A. Glazov , A. Mitov , A. Papanastasious , M. Ubiali , In preparation. https://indico.hep.anl.gov/indico/getFile.py/access? contribId =89&sessionId=28&resId=1&materialId=slides&confId =1161
20. J.C. Collins , F. Wilczek , A. Zee , Phys. Rev. D 18 , 242 ( 1978 ). https://doi.org/10.1103/PhysRevD.18.242
21. J.C. Collins, Phys. Rev. D 58 , 094002 ( 1998 ). https://doi.org/10. 1103/PhysRevD.58.094002
22. A. Kusina , F.I. Olness , I. Schienbein, T. Ježo , K. Kovarík , T. Stavreva , J.Y. Yu , Phys. Rev. D 88 ( 7 ), 074032 ( 2013 ). https://doi. org/10.1103/PhysRevD.88.074032
23. F. Maltoni , G. Ridolfi, M. Ubiali , JHEP 07 , 022 ( 2012 ). https:// doi.org/10.1007/JHEP04( 2013 ) 095 . https://doi.org/10.1007/ JHEP07( 2012 ) 022 . [Erratum: JHEP 04 , 095 ( 2013 )]
24. M. Lim , F. Maltoni , G. Ridolfi, M. Ubiali , JHEP 09 , 132 ( 2016 ). https://doi.org/10.1007/JHEP09( 2016 ) 132
25. R.D. Ball , V. Bertone , L. Del Debbio , S. Forte , A. Guffanti , J. Rojo , M. Ubiali , Phys. Lett. B 723 , 330 ( 2013 ). https://doi.org/10.1016/ j.physletb. 2013 . 05 .019
26. M. Cacciari , M. Greco , P. Nason, JHEP 05 , 007 ( 1998 ). https://doi. org/10.1088/ 1126  6708 / 1998 /05/007
27. B.A. Kniehl , G. Kramer, I. Schienbein, H. Spiesberger , Eur. Phys. J. C 75 ( 3 ), 140 ( 2015 ). https://doi.org/10.1140/epjc/ s1005201533606
28. M. Bonvini , A.S. Papanastasiou , F.J. Tackmann , JHEP 11 , 196 ( 2015 ). https://doi.org/10.1007/JHEP11( 2015 ) 196
29. M. Bonvini , A.S. Papanastasiou , F.J. Tackmann , JHEP 10 , 053 ( 2016 ). https://doi.org/10.1007/JHEP10( 2016 ) 053
30. R.D. Ball , V. Bertone , M. Bonvini , S. Carrazza , S. Forte , A. Guffanti , N.P. Hartland , J. Rojo , L. Rottoli , Eur. Phys. J. C 76 ( 11 ), 647 ( 2016 ). https://doi.org/10.1140/epjc/s100520164469y
31. R.D. Ball , V. Bertone , M. Bonvini , S. Forte , P. Groth Merrild , J. Rojo , L. Rottoli, Phys. Lett. B 754 , 49 ( 2016 ). https://doi.org/10. 1016/j.physletb. 2015 . 12 .077
32. R.D. Ball , M. Bonvini , L. Rottoli, JHEP 11 , 122 ( 2015 ). https:// doi.org/10.1007/JHEP11( 2015 ) 122
33. F. Lyonnet , A. Kusina , T. Ježo , K. Kovarík , F. Olness , I. Schienbein , J.Y. Yu , JHEP 07 , 141 ( 2015 ). https://doi.org/10.1007/ JHEP07( 2015 ) 141
34. H. Abramowicz et al., Eur. Phys. J. C 75 ( 12 ), 580 ( 2015 ). https:// doi.org/10.1140/epjc/s1005201537104
35. F.D. Aaron et al., Eur. Phys. J. C 65 , 89 ( 2010 ). https://doi.org/10. 1140/epjc/s1005200911900
36. H. Abramowicz et al., Eur. Phys. J. C 73 ( 2 ), 2311 ( 2013 ). https:// doi.org/10.1140/epjc/s1005201323113
37. H. Abramowicz et al., JHEP 09 , 127 ( 2014 ). https://doi.org/10. 1007/JHEP09( 2014 ) 127
38. S. Carrazza , A. Ferrara , D. Palazzo , J. Rojo , J. Phys . G 42 ( 5 ), 057001 ( 2015 ). https://doi.org/10.1088/ 0954 3899/42/5/057001
39. V. Bertone , S. Carrazza , N.P. Hartland , Comput. Phys. Commun . 212 , 205 ( 2017 ). https://doi.org/10.1016/j.cpc. 2016 . 10 .006
40. V. Bertone et al., JHEP 08 , 050 ( 2016 ). https://doi.org/10.1007/ JHEP08( 2016 ) 050
41. F. James , M. Roos , Comput. Phys. Commun . 10 , 343 ( 1975 ). https://doi.org/10.1016/ 0010  4655 ( 75 ) 90039  9
42. F.D. Aaron et al., JHEP 09 , 061 ( 2012 ). https://doi.org/10.1007/ JHEP09( 2012 ) 061
43. R.D. Ball et al., JHEP 04 , 040 ( 2015 ). https://doi.org/10.1007/ JHEP04( 2015 ) 040
44. M. Bonvini , S. Marzani , T. Peraro, Eur. Phys. J. C 76 ( 11 ), 597 ( 2016 ). https://doi.org/10.1140/epjc/s1005201644456
45. D.B. Clark , E. Godat , F.I. Olness , Comput. Phys. Commun . 216 , 126 ( 2017 ). https://doi.org/10.1016/j.cpc. 2017 . 03 .004
46. M. Botje, ( 2016 ). arXiv: 1602 .08383 [hepph]
47. M. Botje , Comput. Phys. Commun . 182 , 490 ( 2011 ). https://doi. org/10.1016/j.cpc. 2010 .020
48. S. Lomatch , F.I. Olness , J.C. Collins, Nucl. Phys. B 317 , 617 ( 1989 ). https://doi.org/10.1016/ 0550  3213 ( 89 ) 90535 X
49. F. Caola , S. Forte , J. Rojo , Nucl. Phys. A 854 , 32 ( 2011 ). https:// doi.org/10.1016/j.nuclphysa. 2010 . 08 .009