Asymmetric heavy-quark hadroproduction at LHCb: Predictions and applications

We present a phenomenological analysis of asymmetric bottom- and charm-quark production within the LHCb acceptance relevant for $pp$ collisions at $\sqrt{s} = 13 \, {\rm TeV}$. Predictions are provided for both anti-$k_t$ bottom- and charm-jet pairs, which are kept differentially with respect to the invariant mass of the jet pair. It is quantified how data in this region can provide sensitivity to the couplings of the $Z$ boson to heavy quarks, and we investigate what precision is needed to compete with LEP. We also discuss how asymmetry and rate measurements can provide constraints on a particular class of new-physics models, which contains gauge bosons with small/moderate couplings to light/heavy quarks and masses of the order of $100 \, {\rm GeV}$. Predictions are obtained including all relevant QCD and QED/weak contributions up to next-to-leading order, which have been implemented in a Fortran code which allows to directly compute the asymmetric cross sections. We provide all relevant analytic formulas for our computations.


Introduction
The production of bottom and charm quarks at high-energy colliders is a topic of considerable interest. While not directly observed, these quarks fragment into unstable bottom and charm hadrons with a typical mean lifetime of 10 −12 s. As a consequence of the short but finite lifetime, bottom and charm hadrons decay within the detector at a location which is displaced from the primary collision point. This distinct experimental signature can be used to associate the production of a particle jet in the collision with that originating from a heavy quark, or to improve the efficiency for exclusively reconstructing the heavy-flavour hadron, which in turn has allowed detailed studies of heavy-quark production. A relevant example is the pair-production of bottom-and charm-quarks in e + e − collisions in the vicinity of the Z pole, as studied at both LEP and SLC. Precision measurements JHEP03(2019)166 of both the production rates and the asymmetries in angular distributions of the produced heavy-quarks has allowed to perform precision tests of the Standard Model (SM), and has led to the most stringent constraints on the coupling structure of the Z boson to all quarks but the top quark [1]. Similar studies of the angular asymmetries in heavy-quark production have also been carried out at hadron colliders. In pp collisions at the Tevatron, a measurement of the asymmetry in b-quark pair production has been performed for B-hadrons by the DØ collaboration [2], and also for bottom-quark jet (b-jet) pairs by the CDF collaboration [3]. A measurement of the b-jet pair asymmetry has also been achieved by the LHCb collaboration in pp collisions at the LHC [4].
The asymmetric hadroproduction of heavy-quarks provides important information as compared to what is accessible in e + e − collisions. First, the production mechanisms are entirely different in these collisions, and therefore unique information is provided in hadron collisions. In addition, a measurement of the asymmetry can be performed differentially in the invariant mass of the bb system across a large range of values. This information allows to test a number of new-physics scenarios which are not accessible in e + e − collisions, and there have been a number of relevant phenomenological studies both in the SM and beyond (cf. [5][6][7][8][9][10][11][12][13][14] for instance). It is, however, important to note that the prediction and measurement of heavy-quark asymmetries at hadron colliders also come with a number of challenges. Experimentally it is necessary to account for the effects of pile-up, and to suppress the extremely large background contributions from light-flavour jet production. In addition, the absolute value of the predicted asymmetry is typically quite small. This is mainly a consequence of the large suppression introduced by the symmetric gluon-fusion subprocess for heavy-quark pair production. On the theoretical side, the evaluation of QCD corrections (which are dominant) to heavy-quark production are more complicated at hadron colliders because all external particles are coloured. Obtaining predictions are furthermore computationally more intensive, as the partonic cross sections have to be convoluted with parton distribution functions (PDFs).
The purpose of this work is to provide robust predictions for both bottom-and charmquark jet-pair production in pp collisions at the 13 TeV LHC in the forward direction. There are at least two motivations for focussing on this specific kinematic region. First, the forward regime provides unique opportunities to measure heavy-quark asymmetries at the LHC, because of the increased asymmetry between q andq PDFs present when these partons carry large energy fractions, and the reduced dilution of the symmetric gluon-fusion contribution. Second, the LHCb experiment is a forward detector [15] and able to perform both charge-and flavour-tagging of heavy-quark jets [4,16]. In fact, the recent LHCb measurement of the Z → bb production cross section [17] indicates that finely binned heavy-quark asymmetry measurements in the Z-pole region should be possible as well. As we will show in this article, the latter point is of relevance as there is a long-standing tension between the measured and the SM value of the bb forward-backward asymmetry in e + e − collisions [1]. It is also discussed how measurements of bottom-and charm-quark pair production can provide constraints on new-physics models, which contain gauge bosons with masses of around 100 GeV and small/moderate couplings to light/heavy quarks.

JHEP03(2019)166
The remainder of the paper is laid out as follows. In section 2 we provide details of the theoretical set-up that we use to obtain our numerical results. The SM predictions for the cross sections and the asymmetries are given in section 3 and section 4, respectively. Two applications of our results are presented in section 5 and conclude our article. The technical details of our calculations and their numerical implementation can be found in appendix A.

Theoretical framework
As discussed in the introduction, the goal of the experimental analysis [4] is to measure an asymmetry in the rapidity distributions of b-andb-quarks produced in pp collisions. Experimentally, this has been achieved by requiring the presence of two anti-k t jets [18] which are both charge-(in the presence of a semi-leptonic B decay) and flavour-tagged. This procedure allows to differentiate between b-andb-quark jets and to construct asymmetric observables. Practically, the asymmetry is measured differentially with respect to the invariant mass of the b-andb-jet pair system.
The corresponding theoretical predictions for the inclusive process pp → QQX with Q referring to either a bottom or charm quark in this work are obtained assuming a standard factorisation theorem [19], whereby the hadron-level cross section can be computed by convoluting the individual partonic cross sections with the relevant PDFs. Theoretical predictions for heavy-quark production can be characterised in terms of the perturbative accuracy of the partonic cross sections according to where dσ (n, m) denotes the coupling-stripped differential partonic cross section and α (α s ) is the QED (QCD) coupling. The leading order (LO) contributions to (2.1) correspond to n + m = 2, while the next-to-leading order (NLO) contributions have n + m = 3, and so forth. In this work, we include all numerically relevant NLO corrections to the distributions. 1 The technical details of the calculation and implementation of the various contributions to the partonic cross sections are discussed in appendix A. The techniques to obtain NLO corrections to 2 → 2 processes are by now standard, and in the case of pp → QQX all relevant NLO contributions are known since some time [20][21][22][23][24][25][26][27][28][29] (see also references therein for partial results). We therefore refrain from giving NLO expressions for (2.1) in the main text. Instead we provide an overview of the numerical implementation of our calculations in the following, and discuss the details of the various inputs and scheme choices, which are relevant to the numerical predictions provided in this article.

General set-up
The numerical predictions in this paper are obtained by means of a private Fortran code, which is linked to a number of external libraries: the evaluation of the input PDFs is performed with LHAPDF [30], the numerical integration algorithms of CUBA [31] are used, and

JHEP03(2019)166
one-loop scalar integrals are evaluated with OneLOop [32,33]. An important aspect of the implementation of our analytic calculations is that we separate the contributions to the partonic cross sections into parts that are symmetric and asymmetric under interchange of the final-state heavy quarks. The numerical integration of the symmetric and asymmetric contributions can therefore be performed independently. This approach significantly improves the stability of the numerical integration of the asymmetric cross sections, as only the asymmetric contributions of the partonic cross sections are integrated and the numerical adaption of the integration is specifically optimised for these contributions. The mass of the considered heavy quark is included in our calculations, and we therefore work in a scheme with N F = 4 (3) massless quarks for the bottom-quark (charm-quark) predictions. The O(α 3 s ) corrections to the symmetric cross sections are obtained with the matrix elements [20] implemented in POWHEG BOX [34]. The calculation of the weak boxdiagram corrections of O(αα 2 s ) have instead been obtained using MadLoop [35] as part of the loop-induced module [36] available in MadGraph5 aMC@NLO [37]. 2 All other contributions have been computed with the aid of FeynArts [39] and FormCalc [40], and the relevant analytic formulas for the asymmetric contributions to the partonic cross section are collected in appendix A. In these cases, we have used the technique of phase-space slicing [41] or dipole subtraction [42] to regulate the explicit (implicit) divergences present in the virtual (real) phase-space.
We add that differential O(α 4 s ) results have been first presented for top-quark production in [43], as well as for massless partons to leading colour in [44]. At present, a calculation of b-tagged jets (either massive or massless) is instead not available. However, it can be expected that such predictions will become available in the future when issues related to numerical stability or flavour-tagging of subtraction terms in the next-to-next-to-leading order (NNLO) QCD calculations have been resolved.

Observables
To match the experimental definition of jet observables, we construct anti-k t jets with a radius parameter R = 0.5 and tag them as a Q-jet (Q-jet) if they contain a Q (Q), with Q being the heavy quark. Throughout this work, all observables are computed in terms of these flavour-tagged jets. The label "jet" will however be suppressed, meaning for example that the invariant mass of a b-andb-jet pair will be simply called m bb . If not stated otherwise, we will always place the following kinematic cuts on the flavour-tagged heavy-quark jets referring to these selections as "LHCb kinematic cuts". The requirements on the transverse momentum (p T ) and the pseudorapidity (η) ensure that the jets are reconstructed according to the flavour tagging algorithm in use at LHCb [16]. The cut on the angular separation (φ) between the two flavour-tagged jets in the azimuthal plane ensures that the two heavyquark jets are well separated. This cut therefore avoids configuration which can appear for

JHEP03(2019)166
instance at O(α 3 s ) where both the Q andQ are contained within a single jet (at LO this cannot occur as the heavy quarks are produced back-to-back). The impact of the choice of the angular cut on the predicted asymmetries is discussed in section 4.
The two primary observables of interest are the heavy-quark production cross sections and the corresponding asymmetries. The cross sections are computed differentially in m QQ within the LHCb fiducial region (2.2). The asymmetries are also computed differentially in m QQ , and defined according to Here dσ S (A) refers to the convolution integral of the differential (a)symmetric partonic cross sections with the relevant PDFs, and ∆y = y Q − yQ is the difference between the rapidities of the Q and theQ.

Heavy-quark mass effects
As mentioned above, we retain the effects of the heavy-quark mass throughout our calculations. The following choices for the heavy-quark masses in the on-shell scheme are adopted These values are broadly consistent with the recommendations of the LHC Higgs Cross Section Working Group [45]. When we provide predictions for either cross sections or asymmetries, we do not consider the uncertainties associated to (2.4). The motivation for this is that the mass corrections within the considered fiducial region are typically small (although not negligible), and the resulting ambiguities are small compared to the scale uncertainties. This statement is corroborated in figure 1 (left), which shows LO differential bb cross sections within the LHCb fiducial region (2.2) for different choices of m b . These distributions are obtained with the LUXqed15 [46] central PDF set member with factorisation (µ F ) and renormalisation (µ R ) scales set dynamically to m bb , and the distributions have been normalised to the result obtained with m b = 4.75 GeV. As can be seen from the distribution obtained with m b = 0 the mass corrections amount to 3% to 10% within m bb ∈ [40, 100] GeV. On the other hand, a variation of m b in the range m b ∈ [4.5, 5.0] GeV results in cross-section changes below the percent level. While this study focuses on the case of symmetric bb production within the LHCb fiducial region, similar corrections (although with opposite sign) are observed for asymmetric bb production within this region. In the cc case, the mass corrections remain always below 2%. We also note that the inclusion of mass effects leads to a positive correction to the symmetric cross section within the LHCb fiducial region, while the inclusive cross section within the same invariant mass region receives negative corrections. While to achieve precision predictions in the region of m bb ∈ [40, 100] GeV including mass corrections is clearly important, at larger values of m bb one could alternatively perform the calculation taking the heavy quarks to be massless. Employing a massless scheme JHEP03(2019)166   would have the advantage that logarithmic mass corrections could be resummed, but also has some weaknesses. In this context it is important to recall that the measurement is performed by requiring the presence of two well separated flavour-tagged anti-k t jets according to (2.2). With such kinematic requirements, the phase-space regions where the NLO fixed-order calculation receives large logarithmic corrections (for example, due to the presence of g → bb collinear enhancements) are avoided, and the effects of resumming these types of contributions is therefore typically small. 3 Another consideration is that the prediction of flavour-tagged anti-k t jets is only infrared (IR)-safe for the massive calculation. Due to the presence of wide-angle g → QQ splittings of soft gluons, the massless calculation is IR-unsafe [47]. We therefore provide our predictions including the full mass effects up to NLO, such that it is possible that numerical predictions computed with the massive NNLO calculation [48] can be added consistently at a later date. Alternatively one could consider a flavour-tagging algorithm which is IR-safe and achievable experimentally.

PDFs, input parameters and scale variation
As a baseline PDF set in this work, we use the variant of NNPDF31 nlo as 0118 [49] where the charm-quark PDF is generated purely perturbatively. This is a variable flavour number scheme set with N F = 5, which is therefore evolved (both PDFs and α s ) with five active flavours above the b-quark mass threshold. As discussed throughout this section, we deliver predictions for bottom-and charm-quark pair production at NLO including both QCD and QED/weak corrections. For consistency, these predictions should be obtained by convoluting the partonic cross sections with PDFs which have been extracted from a PDF fit including both QCD and QED effects. There are two important points related to the choice of PDFs which we describe below.

JHEP03(2019)166
First, when calculating bb (cc) production using N F = 4 (3) active flavours, there is a mis-match at O(α 3 s ) between the perturbative cross-section calculation and our input PDF set which uses N F = 5. To account for this we include the relevant compensation terms following [50]. Second, our baseline PDF does not include a photon PDF or the effects of a joint QCD-QED evolution [51]. There has recently been quite some activity in precisely determining the photon PDF [46,[52][53][54][55] and also a number of studies of electroweak (EW) effects in top-quark pair production have been presented [56][57][58]. We have studied the impact of the latter two types of contributions for bottom-and charm-quark pair production, and found these effects of very limited importance. This is demonstrated in figure 1 (right) where the contributions of gluon-fusion, quark-annihilation, and gluon-photon scattering to the bb cross section are shown. The given predictions are obtained at LO with the LUXqed15 PDF set using µ F = µ R = m bb and employing the reference cuts (2.2) at the 13 TeV LHC. In the considered invariant mass range, we find that the photon-induced contributions lead to effects at the permille level. Compared to the uncertainty of the total cross section (which is around 10%), these effects are thus entirely negligible. 4 We have therefore chosen to use a PDF set based on NLO QCD which does not include a photon PDF. A consequence of ignoring the mixed QCD-QED evolution effects is to slightly overestimate the uncertainty due to µ F variation used to assess the theoretical uncertainty of our predictions.
In this work, we use the following input parameters: m h = 125 GeV, m t = 173 GeV, M W = 80.385 GeV, Γ W = 2.085 GeV, M Z = 91.1876 GeV, Γ Z = 2.4952 GeV as well as G F = 1.16638 · 10 −5 GeV −2 . Employing this input and including the dominant one-and two-loop universal corrections to the ρ-parameter [59], we have derived the following values for the square of the sine of the weak mixing angle sin 2 θ w = 0.2293 and the electromagnetic coupling α = 1/128.55. For the values of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, we take |V us | = |V cd | = 1 − |V ud | 2 = 1 − |V cs | 2 = 0.23, |V tb | = 1, while all other elements are set to zero. For the evaluation of the O(αα 2 s ) corrections we use a complex massscheme [60], accounting for the width effects of the Z boson. In the latter case, the (complex) value of the weak mixing angle is derived from the complex W -and Z-boson masses.
To assess the uncertainty due to missing higher-order corrections, scale variations are performed by changing both µ F and µ R independently by a factor of two around a reference scale µ 0 , with the constraint that 1/2 < µ F /µ R < 2. Predictions are obtained for the two following choices of the reference scale corresponding to the invariant mass of the heavy-quark jet pair and the mean transverse energy of the heavy-quark jets, respectively. When observables such as the asymmetry defined in (2.3) or a cross-section ratio between heavy quarks are considered, the scale variations are computed by correlating the scales between numerator and denominator.
To conclude this section, we note that it is straightforward to also provide predictions for stable top quarks with our numerical set-up. The production of top quarks at forward rapidities is of considerable interest [61][62][63], and there has been significant experimental

JHEP03(2019)166
progress towards performing precise forward measurements of this process [64,65]. However, as only partial event reconstruction is possible, studying the phenomenology of stable top quarks within the LHCb acceptance is not very practical. For a detailed discussion of the leptonic top-quark asymmetries at LHCb, we refer the interested reader to [63].

Cross-section predictions
The main goal of this work is to provide reliable predictions for asymmetric heavy-quark production in the fiducial region (2.2). Having a clear understanding of the associated cross sections is, however, an important ingredient of this analysis as well. From the theoretical point of view, it is important to validate the absolute heavy-quark jet rates as well as the shape of the invariant mass distributions, in particular in the region around the Z pole. Experimentally, measurements of the cross sections may give handles on the (charged) flavour-tagging efficiency and the mis-tag rates, as well as providing an important validation of the jet-energy scale and resolution corrections. The differential heavy-quark cross sections may also lead to constraints on new-physics models which contain light gauge bosons, a point we will return to in section 5. The remainder of the current section is dedicated to the study of the symmetric distributions. In the figure we have focussed on the region of m QQ ∈ [60, 300] GeV, where the differential cross sections span several orders of magnitude. The scale uncertainties of the NLO distributions are about 10%, which represents a marked improvement with respect to the LO results. We also find that the NLO distributions corresponding to the two scale choices (2.5) lead to consistent results, and tend to lie within the uncertainty bands of the LO distributions. In the region of m QQ ≥ 100 GeV the cross sections are entirely dominated by the QCD contributions, and there is a 5% to 10% difference between the central values of the NLO results obtained with µ 0 = m QQ or µ 0 = E T,Q . An improvement in the perturbative stability of the predictions in this region would require the inclusion of O(α 4 s ) corrections, either through a fixed-order calculation or by performing resummation (see for example [66]). The fiducial cross sections within the invariant mass bin m QQ ∈ [250, 300] GeV are approximately 30 pb. Assuming an integrated luminosity of 5 fb −1 and a signal efficiency QQ = 0.6%, these numbers imply that a relative statistical uncertainty of about 1% may be achievable with future LHCb data.

Cross sections
LHCb has recently performed a measurement of the process pp → Z → bb + jets at √ s = 8 TeV [17]. This measurement is performed differentially with respect to m bb in bins of width 4 GeV in the Z-pole region, suggesting that future measurements of inclusive heavy-quark pair production will also be possible with similarly fine binning. In figure 3 we provide predictions for both bottom (left) and charm (right) jet-pair production, focussing on the invariant-mass region of m QQ ∈ [60, 120] GeV. Besides the total rates also spectra of the various subprocesses are shown. At LO we display the purely QCD (α 2 s ) and EW (α 2 ) contributions, while at NLO we have chosen to depict various combinations of mixed QCD-EW corrections. When considering the EW corrections to the LO QCD processes (αα 2 s ) only the values of the NLO coefficient is displayed, where we have separated the impact of the QED and weak corrections. The QED corrections in this case are negative, and thus the absolute values of the NLO coefficient are shown. In the case of the QCD corrections to the LO EW processes (α 2 α s ) the sums of the LO and NLO coefficients are given, where we have also displayed the result when including only initial-state radiation (ISR) from QCD (labelled as ISR only).
The LO QCD contributions are by far dominant, while the LO EW contribution only becomes relevant (reaching roughly 10%) in the region of m QQ ∈ [85, 95] GeV. The QED corrections to the LO QCD process are negative and more important in the case of charmquark production, 5 where these effects amount to half a percent of the total cross section. For both bottom-and charm-quark production, the weak corrections are negligibly small. The QCD corrections to the LO EW process have a more important impact on the obtained results. This is primarily due to the contribution of hard QCD corrections to the heavyquark lines, where an emitted gluon is not reconstructed as part of the heavy-flavour jet. A consequence of such resonantly enhanced events which lose energy via the emission of such  As a final comment to figure 3 (right), we note that the scale variation present for the cc prediction of O(α 2 α s ) in the bin m cc ∈ [60,65] GeV is a genuine effect in our calculation. It originates from the NLO contributions associated to qg → QQq (and the corresponding collinear counterterm) with purely photon exchange -this term can be considered as part of the QED correction to the LO transition gγ → QQ. The scale dependence of the NLO corrections would normally compensate that of the LO contribution, which is however absent in our computation due to the missing photon PDF. This increased scale dependence has a permille effect on the total cross section.
To conclude the discussion of the differential cross sections, we present predictions for b-jet pair production within the LHCb acceptance as a function of the minimum angular separation φ min bb between the b-jets. In the LHCb analysis [4], a cut φ bb > 2.6 is imposed, and said to provide improved sensitivity to the asymmetry of the signal process by enhancing "non-gg production mechanisms". In section 4 we will show that the value of this cut is not so important for the signal process, however it is still likely to be relevant for reducing the background contribution from light-flavoured jets. An experimental study of this distribution may therefore be useful when studying/reducing the contamination of background events. The relevant predictions for pp collisions at

Cross-section ratios
In addition to the measurements of the bb and cc cross sections discussed in the last section, it is also of interest to perform measurements of the cross-section ratio between the different heavy-quark types. As the theoretical predictions for the cross-section ratios are very precise, these measurements will provide an important experimental benchmark for testing and validating the (charged) flavour-tagging efficiency and mis-tag rates.   Our predictions for the cross-section ratio of c-and b-jet pairs at the 13 TeV LHC in the phase-space region (2.2) are shown in figure 5. These distributions are obtained with µ 0 = m QQ , and the uncertainty has been evaluated by correlating the scale variations between the charm-and bottom-quark predictions. In the considered m QQ range, the ratio between the cc and bb cross sections is below 1. The observed deviation of the ratio from 1 can be attributed to the mass dependence of the LO cross section -see figure 1 (left) -and also to the different EW charges of up-and down-type quarks which affects the ratio both close to and away from the Z peak. In figure 5 (left), the NLO ratio is also displayed for the case that the O(αα 2 s ) corrections have been removed. These corrections arise dominantly in the form of QED corrections to the gg → QQ subprocess. They are negative and amount to effects of the order of e 2 Q · 1% on the spectra, where e Q denotes the electric charge of the heavy quark Q. While the O(αα 2 s ) contributions thus have a negligible impact at the level of the cross sections, they have a visible effect on the ratio of the symmetric rates.
As discussed in section 2, we have chosen to use PDFs that do not include a photon PDF, and as a result photon-initiated contributions are not included in our computations. To assess the potential uncertainty due to these missing contributions, we have recomputed the ratio of the cc and bb spectra at LO with the LUXqed15 PDF set. An uncertainty is then calculated according to where the second ratio is computed excluding all photon-initiated contributions. The modulus of this uncertainty is then added linearly to the scale uncertainty, both in the positive and negative directions. The corresponding results are computed at NLO accuracy with µ 0 = m QQ and shown in figure 5   an intrinsic charm-quark PDF. Using either of these PDF sets for the current predictions would introduce some level of inconsistency. At present, we therefore recommend to use the conservative uncertainty including δR γ when comparing our results to future data.

Asymmetry predictions
In this section we provide differential predictions for the asymmetries as defined in (2.3). These predictions are obtained by separately computing both the numerator and denominator of this expression at NLO, i.e. including terms up to n + m = 3 in the expansion defined in (2.1). The corresponding LO results are obtained by including terms up to n + m = 2. The exception is that we also take into account the O(α 3 s ) contribution to the numerator when evaluating the asymmetry at LO. This procedure is motivated by the well-known fact that the O(α 2 s ) corrections do not generate an asymmetry in the SM, so that the O(α 3 s ) terms should be considered the leading contribution to the asymmetry. In fact, these terms are numerically dominant except for m QQ values close to the Z pole.
In figure 6 our results for the bb (left) and cc (right) asymmetries are presented. The NLO predictions corresponding to the two different scales choices (2.5) lead to consistent results across the considered mass range. Close to the Z peak it is found that the NLO corrections have an important impact on the absolute value of the asymmetries as well as the uncertainty estimates. In both the low-and high-mass regions, the uncertainty estimate from the LO prediction is artificially small and should not be considered robust.
Compared to the results presented in [14], our current predictions include the following improvements. First, a jet definition consistent with the LHCb flavour-tagging algorithm [16] is used throughout. Second, based on the recent measurement of bb production in the vicinity of the Z peak [17], numerical predictions for both the cross section and asymmetry are provided in fine bins in this invariant-mass region. Third, more precise theoretical predictions are obtained by including a number of subleading NLO corrections, which were  previously absent. Fourth, the numerical predictions are computed with updated PDFs which include LHC data, and have been obtained with two reference scale settings which lead to a more reliable uncertainty estimate. These improvements should allow for a more precise comparison to data, which can in turn be used to perform more stringent tests of the SM as well as new physics. Two applications along these lines are discussed in section 5. Before discussing these applications, it is important to estimate the potential sensitivity of future experimental measurements. The original measurement of the b-jet pair asymmetry at LHCb was performed with an integrated luminosity of 1 fb −1 collected at 7 TeV [4]. Results were presented in the m bb bins of [40,75] GeV, [75,105] GeV and [105, 300] GeV, and the measurement was statistically limited in each bin. To estimate the statistical sensitivity expected at 13 TeV, we compute the corresponding uncertainties via where A denotes the central value of the NLO asymmetry obtained with µ 0 = m QQ , and N is the expected number of events within the data set. To calculate N , we use our cross-section predictions, assume a data set of 5 fb −1 , and further apply experimental efficiencies for the reconstruction of a pair of charged-and flavour-tagged jets of bb = 0.6% and cc = 0.3%. The values of these efficiencies are obtained by inverting (4.1) for the 7 TeV measurement [4] using the corresponding central NLO prediction at 7 TeV. We note that the value of bb = 0.6% corresponds to a factor of two improvement compared to what has been achieved in the original measurement.
The results of our sensitivity study are shown in figure 7, where the projections for the statistical uncertainties (4.1) are overlaid on the predictions for the b-(left) and c-jet (right) pair asymmetries. This study indicates that a significant improvement in statistical precision will be achievable with future data sets, and that finely binned measurements close to the Z pole should be possible. This is a consequence of the higher cross sections, the increased data sample size, and the assumption about the improved signal efficiency. In the event that a data sample of 50 fb −1 is collected at LHCb [67] (such as in the high-luminosity phase of the LHC), it is likely that measurements of the heavy-quark asymmetries will be systematically and not statistically limited. We conclude this section by returning to the choice of the angular cut φ QQ used in defining the fiducial region (2.2). As mentioned in section 3, the motivation for introducing this cut is to increase the sensitivity to the asymmetry by enhancing "non-gg production mechanisms". To assess this statement, we study the impact of the choice of φ min QQ on the observable σ A /σ 1/2 S , where σ S (A) is the (a)symmetric production cross section. The motivation behind this definition is that the significance of a statistically limited measurement of the asymmetry is approximately A/δA stat . Our definition is therefore useful as it estimates the overall statistical sensitivity to the asymmetry measurement itself, rather than just the asymmetry. This is relevant because, while the value of the asymmetry may increase as the value of the cut φ min QQ is increased, the number of analysed events simultaneously decreases. Our predictions for σ A /σ 1/2 S as a function of φ min bb are shown in figure 8. The two different sets of predictions correspond to the results restricted to the invariant mass bins m bb ∈ [75, 105] GeV and m bb ∈ [105, 300] GeV. The obtained distributions are close to flat as φ min bb increases, indicating that from a statistical point of view the sensitivity to the asymmetry is not improved by requiring larger φ min bb values. We have therefore chosen to provide predictions for φ min QQ = 2.6, which matches the original value advocated in [4]. It is worth noting that the choice of this cut may also be important for background rejection (i.e. from light jets). In the far future, if the asymmetry measurements becomes systematically limited, it may be useful to perform a dedicated study of this issue.

Applications
In this section we present two applications of our calculations of heavy-quark production. We will first discuss the model-independent constraints that future LHCb measurements of the ratio of the bb and cc asymmetry may allow to set on the couplings of the Z boson to bottom-and charm-quark pairs. We will compare the results of our sensitivity studies to the constraints on the Zbb and Zcc couplings that arise from the Z-pole measurements performed at LEP [1]. Our second application consists in using the recent LHCb measurement of Z → bb production [17] to put constraints on the new-physics model proposed in [68] which aims at explaining the long-standing LEP anomaly of the forward-backward asymmetry of the bottom quark.

Constraints on Zbb and Zcc couplings
Similarly to the top-quark asymmetry, the dominant contribution to the asymmetry of bottom-and charm-quark arises from QCD for most values of the invariant mass of the heavy-quark pair. An important exception is the mass region close to the Z-pole, where the double-resonant contribution from Z-Z interference becomes dominant, and accounts for the bulk of the total asymmetry [12][13][14]. Measurements of the bottom-and charm-quark asymmetry can therefore be used to set limits on the Zbb and Zcc couplings [13].
In order to put model-independent constraints on the Z-boson couplings to bottom and charm quarks, we consider in the following the ratio R b/c = A b /A c of the asymmetry in bb and cc production restricted to the mass bin [75,105] GeV. This ratio can be predicted to high accuracy in the SM [14], since many uncertainties cancel between numerator and denominator. For pp collisions at √ s = 13 TeV and employing the standard LHCb cuts (2.2), we find within the SM the result The given central value corresponds to the reference scale choice µ 0 = m QQ and the quoted uncertainty of around 5% includes scale variations as described in section 2.4 and PDF uncertainties. The dominant source of uncertainty arises from missing higher-order QCD corrections, meaning that the stated total uncertainty is in principle improvable by including NNLO corrections. We add that using instead the scale setting µ 0 = E T,Q leads to

JHEP03(2019)166
a central value of R SM b/c that agrees within errors with that in (5.1) and to a comparable total uncertainty. It is worth noting that the above prediction for R b/c corresponds to a specific bin around the Z-boson resonance, and that the prediction is sensitive to the exact location and width of the bin. It will therefore be important to carefully assess the impact of potential bin-to-bin migration effects as part of the experimental measurement of the ratio R b/c . Our code can be made available upon request, and can be used to assess the related systematic uncertainty.
The experimental measurements of the EW Z-pole observables obtained at LEP can be used to precisely extract the Z-boson couplings to all SM quarks but the top quark. In the case of the Zbb and Zcc couplings the combined results are [1]  In figure 9 we show the relative deviations of R b/c in four different planes of Zbb and Zcc couplings. Overlaid in green are the 68% CL regions that follow from the LEP measurements (5.2) and (5.3). The SM and best-fit points are indicated by black dots and black crosses in the figure. Numerically, we find that the best fits lead to relative deviations of −2.4%, 1.1%, −0.4% and −0.
L (lower left) and g b R -g c R (lower right) plane, respectively. These numbers indicate that for future LHCb measurements of the bb and cc asymmetries to reach the sensitivity of the existing LEP constraints on the Zbb and Zcc couplings, determinations of the ratio R b/c at the percent level are needed. Notice that to reach this goal not only the experimental precision but also the theoretical accuracy of the SM predicition (5.1) needs to be improved. Such a theoretical improvement would require the inclusion of NNLO QCD corrections in the prediction of R b/c , which is technically viable in view of [43].

Constraints on new light gauge bosons
Precision measurements of the gauge sector have shown agreement with expected SM properties at the permille level. Among the many observables, the bottom-quark forwardbackward asymmetry A b FB measured at LEP however presents a 3σ deviation with respect to the values expected in the SM [1]. While this deviation could be a result of statistical fluctuations, it is intriguing since it also could be associated with a large modification of the right-handed Zbb coupling cf. (5.2) and (5.4) , which can for instance arise if the Z boson is mixed with additional neutral gauge bosons.

JHEP03(2019)166
A recently proposed model [68] that aims at explaining the long-standing LEP anomaly of A b FB contains a new U(1) D boson (the corresponding mass eigenstate will be called Z in what follows), which couples with opposite charges to the right-handed components of the bottom and charm quarks. The low-energy spectrum of the model also includes two Higgs doublets, a singlet, and a charged and a neutral vector-like singlet, the latter being a good dark matter candidate. The phenomenology of the Z and the Z bosons is fully described by the masses M Z and M Z of the two gauge bosons, the new coupling constant g D , the sine s α of the mixing angle α of the neutral weak eigenstates and the mixings of the bottom (charm) quark with a heavy vector-like bottom (charm) partner, parameterised by the four variables s b,L , s b,R , s c,L and s c,R .
The following values of the U(1) D gauge coupling and the mixing parameters These couplings lead to a χ 2 /dof of 1.6, which represents a visible improvement compared to the χ 2 /dof value quoted after (5.4). Constraints on the Peskin-Takeuchi parameter T now put a bound on the size of the allowed mass splitting M Z − M Z [68]. For the s α value given in (5.4), one finds that the constraint T = 0.07 ± 0.12 [70] translates into the following 95% CL limit on the mass of the new gauge boson M Z ∈ [91.2, 174] GeV . (5.7) As pointed out in [68], the presence of a Z boson in this mass range is subject to several constraints. The first constraint comes from the CMS search [71] for narrow spin-1 resonances decaying to a qq pair in association with a high-transverse momentum jet from ISR. Other relevant constraints arise from the Z → + − searches [72,73]. For the benchmark parameters (5.5) the combination of the constraints [71][72][73] can however be shown to be fulfilled for most of the Z -boson masses in (5.7). In fact, the search [71] features an 2.9σ local excess for dijet invariant masses around 115 GeV, which has been interpreted in [68] as a Z boson in the U(1) D model with M Z 115 GeV and (5.6). In the following, we point out that Z boson with the properties (5.6) and (5.7), can also be probed by the LHCb measurement of Z → bb production in the forward direction [17]. To this purpose, we show in figure 10 four different dijet mass (m jj ) distributions predicted in the U(1) D model (blue curves). The chosen parameters are given in (5.5). The background-subtracted dijet mass distribution 6 as measured by the LHCb collaboration in [17] (black error bars), the SM prediction (red curves) and the one standard deviation total uncertainty band (grey bands) is also displayed. Δχ 2 Figure 11. ∆χ 2 distribution in the U(1) D model as a function of the Z -boson mass following from the LHCb measurement [17]. The benchmark parameters (5.6) have been used to obtain the shown predictions. The dashed black line corresponds to ∆χ 2 = 3.84, i.e. the 95% CL for a Gaussian distribution.
the CMS measurement [71] observed for dijet masses around 115 GeV. The corresponding χ 2 /dof is 1.9 and thus slightly better than the SM fit which leads to χ 2 /dof = 2.0. While the finding that both the CMS and LHCb measurement may indicate that a light new gauge boson is hiding in the data is probably accidental, we emphasise that the two-sided bound (5.8) on the mass of the Z boson is stronger than the limit that derives from a combination of the searches [71][72][73].
The above applications of our SM results presented in sections 3 and 4 show that the LHCb experiment can provide unique probes of new physics in heavy-quark production due to its efficient triggering, excellent vertexing and accurate event reconstruction. While already a handful of similar proposals of such "exotic" new-physics searches exist that specifically exploit the remarkable capabilities of LHCb (see e.g. [6,13,67,[74][75][76][77][78][79][80]), we believe that further research in this rich and developing field can turn out to be potentially very rewarding.

A Analytic results
In this appendix, we provide predictions for the hadronic process pp → QQX assuming the factorisation theorem of the form where the universal PDFs f i/p (x, µ F ) describe the probability of finding the parton i in the proton with longitudinal momentum fraction x and µ F is the factorisation scale. In analogy to (2.1), the differential partonic cross sections dσ ij may be written as an expansion in terms of the electromagnetic and strong coupling constant according to where dσ (n, m) ij denotes the differential partonic cross sections for the initial state ij with the α and α s factors stripped off.
The differential cross sections can be further decomposed to isolate the contributions which are asymmetric under interchange of the final state heavy quark and antiquark according to Here the notation indicates that in the process labelled by ij → QQX (ij →QQX) the angle θ corresponds to the scattering angle of the heavy quark (heavy antiquark) in the partonic centre-of-mass (CM) frame. The benefit of this decomposition is that the symmetric and asymmetric contributions to the cross sections can be numerically integrated separately, which substantially improves the efficiency of the numerical evaluation. The results in this work have been obtained at NLO accuracy (n + m = 3) for both symmetric and asymmetric contributions. Relevant results for the differential cross sections to this accuracy have previously been obtained in [20][21][22][23][24][25][26][27][28][29]81], and in many cases the analytic results have been given. The purpose of this appendix is to collect the analytic results for asymmetric bottom-and charm-quark pair production, suitable for direct numerical integration. We will provide analytic expressions for both the (renormalised) virtual and real emission contributions to the partonic cross section for various subprocesses, which contain explicit and implicit divergences, respectively. A number of these contributions contain only soft divergences, and so we have also included a soft function which describes the radiation of soft gluons integrated in phase space up to a cut E cut in the gluon energy. This function is suitable for applying the technique of phase-space slicing [41], which is found to be stable for these types of processes. In cases where both soft and collinear divergences are present, we have also regularised the numerical integration using dipole subtraction [42]. In the next subsection, we introduce our notation for describing the kinematics of two-and three-body partonic final states, and then list the relevant formulas ordered by their powers in the expansion (A.2).

A.1 Kinematics and notation
The partonic cross section for heavy-quark pair-production receives Born-level contributions from gluon fusion, quark annihilation as well as gluon-photon and photon-photon scattering. As an example of a partonic 2 → 2 subprocess, we consider quark annihilation where the four-momenta p 1,2 of the initial-state partons can be expressed as the fractions x 1,2 of the four-momenta P 1,2 of the colliding protons. The partonic cross section is a function of the kinematic invariantŝ and momentum conservation implies thatŝ +t Q +û Q = 0. In addition toŝ,t Q andû Q , we also use the velocity of the heavy quark and scattering angle to write our results, where θ denotes the angle between p 1 and p 3 in the partonic CM frame. Notice thatt which implies that c = (t Q −û Q )/ŝ, and as a result the variable c is strictly speaking not needed when writing our results. In some cases we will however use c, because the obtained expressions turn out to be more compact than those written in terms ofŝ,t Q andû Q . In addition to these kinematic variables, it also useful to introduce the following mass variables The complex squared-mass µ Z is introduced as the Z boson is treated as an unstable particle throughout our calculation. This is necessary when describing bottom-and charmquark pair production in the vicinity of the Z-boson resonance. The NLO corrections also involves the evaluation of 2 → 3 real emission processes of the form where again we have used the quark-annihilation subprocess as an example. The analytic formula for these processes are provided in terms of the following dimensionless variables All 2 → 3 processes can be characterised by five independent scalar quantities [22]. We emphasise that below we will write the 2 → 3 results such that the obtained expressions become as short as possible, and as a result our formulas will involve more than five independent y ij parameters.

A.2 O(α 2 ) contributions
The asymmetric O(α 2 ) effects arise from the interference between the partonic processes qq → Z/γ → QQ. The relevant diagram with s-channel Z-boson exchange is shown on the left-hand side in figure 12. In agreement with [82], we obtain for the corresponding asymmetric differential cross section the following compact expression are the axial-vector and vector coupling of the Z boson to a fermion f . These couplings depend on the third component T f 3 = ±1/2 of the weak isospin, the electric charge e f , and the sine s w and cosine c w of the weak mixing angle.

A.3 O(αα s ) contributions
In order to obtain the O(αα s ) contributions one has to consider interference contributions between t-channel W -boson (and would-be Goldstone boson) exchange and s-channel gluon JHEP03(2019)166 exchange. A diagram of this kind is given on the right in figure 12. In the case of cc production from a dd initial state, we arrive at where V cd denotes the relevant CKM matrix element. Our result (A.14) agrees with the expression given in [13]. In the case of asymmetric bb production, the W -boson mediated tchannel contributions are strongly suppressed either by the small CKM element V ub or by a bottom-quark PDF. Consequently, we include (A.14) in our numerical analysis only in the case of charm-quark pair production. As these corrections are numerically small, we have included neither the QCD nor the QED/weak correction to this process in our predictions.

A.4 O(α 3 s ) contributions
There is no asymmetric contribution to the production of heavy-quark pairs at O(α 2 s ). Starting at O(α 3 s ), however, quark annihilation qq → QQ(g) as well as flavour excitation qg → QQq receive charge-asymmetric contributions. The gluon-fusion gg → QQX subprocess must be convoluted with a symmetric initial state to provide a hadronic cross-section prediction, and so does not lead to an observable asymmetry.
Charge conjugation invariance can be invoked to show that, as far as the virtual corrections to qq → QQ are concerned, only the interference between the lowest-order and the QCD box graphs contributes to the asymmetry at O(α 3 s ). An example of a Feynman diagram that furnishes a contribution is shown on the left-hand side in the upper row of figure 13. The corresponding virtual corrections can be written as 15) with N c = 3 and d 2 abc = 40/3 colour factors. Here d abc = 2Tr T a , T b T c , while T a are the colour generators normalised such that Tr T a T b = δ ab /2. The one-loop function appearing in (A.15) is given by  Figure 13. Representative Feynman diagrams that contribute to the asymmetric production cross section of heavy-quark pairs at O(α 3 s ). Upper row: the two-particle cut (right) describes the interference of the one-loop box diagram with the tree-level graph, while the three-particle cut (left) corresponds to the interference of final-state with initial-state gluon corrections. Lower row: threeparticle cuts that represent a production of heavy quarks via flavour excitation.
where our definition of the Passarino-Veltman scalar integrals B 0 , C 0 and D 0 follows that used in FormCalc. We have verified that the formulas (A.15) and (A.16) recover the analytic results given in [83]. where is the relevant soft function. Here = (4 − d)/2 arises from dimensional regularisation in d dimensions, µ denotes the corresponding renormalisation scale, and we have furthermore introduced and (A.20) can be shown to agree with the expressions presented in [83]. Note that the IR 1/ pole in (A. 19) cancels against that in the virtual corrections (A.15) so that the sum of the virtual and soft contributions is IR finite and can be numerically integrated in four dimensions.
The asymmetric O(α 3 s ) contribution to the heavy-quark production cross section that is associated to the flavour excitation process can be obtained from the result (A.17) by crossing, i.e. interchanging the indices 2 ↔ 5 in the variables y ij as defined in (A.10). Examples of relevant Feynman diagrams are given in the lower row of figure 13. Noting a difference in the colour factor for averaging over the initial-state gluon with respect to [ The same result also holds in the case of the partonic reactionqg → QQq. In contrast to the asymmetric contribution from qq → QQg, the flavour excitation processes qg → QQq andqg → QQq are IR finite.  Figure 14. Left: possible two-and three-particle cuts that contribute via qq → QQ, qq → QQg and qq → QQγ to the asymmetric cross section in heavy-quark production at O(αα 2 s ). Right: example of a Feynman diagram that leads to a photon-initiated production asymmetry of heavy quarks at O(αα 2 s ).
Here the factor 12/5 arises from the ratio of QED and QCD colour factors (N 2 c − 1)/4 = 2 and d 2 abc /16 = 5/6. In the case of quark annihilation, we find the following relation between the corrections of O(αα 2 s ) and O(α 3 s ) to the differential asymmetric cross sections irrespectively of whether the contributions are virtual, real or soft. The additional overall factor of 3 reflects the three possible attachments of the photon in diagrams like the one shown on the left-hand side in figure 14. In the case of the qg-initiated transition only two different photon attachments are possible so that and similarly forqg → QQq. Our formulas (A.23) and (A.24) agree with the findings of the article [83]. Finally, in the case of the photon-initiated process qγ → QQq, which receives contributions from Feynman diagrams such as the one displayed on the right in figure 14, 25) and the same result applies to theqγ initial state. Notice that the factor of N 2 c − 1 = 8 arises from averaging over the photon in the initial state rather than the gluon.
As in the case of pure QCD, the O(αα 2 s ) corrections associated to Z-boson exchange receive contributions from both virtual and real corrections. For the interference contributions of box diagrams with tree-level s-channel exchange graphs, we obtain the following expression for the asymmetric heavy-quark pair production cross section. The loop function A(v, w) has already been given in (A. 16). It arises in the context of (A.26) from two-particle cut JHEP03(2019)166 q q Z Z q q g g g g qq q g g q q q Q Q Figure 15. Examples of graphs that affect the production asymmetry of heavy quarks at O(αα 2 s ) and involve the exchange of a Z boson. Consult the main text for additional explanations.
diagrams like the one depicted on the left in figure 15. The function B(v, w) is instead related to two-particle cuts in graphs of the type shown on the right-hand side of the same figure. In terms of standard Passarino-Veltman scalar integrals B 0 , C 0 and D 0 , we arrive at the result Notice that the these expressions are a function of the complex squared-mass µ Z , which is necessary to account for the width effects in region close to the Z-boson resonance (this is not required in the case of top-quark pair production since 2m t > m Z ). If the calculation is performed in the complex-mass scheme, the couplings v f become complex as discussed for instance in [60].
The real emission and soft contributions of O(αα 2 s ) can again be obtained by rescaling the corresponding QCD results. We define

JHEP03(2019)166
In terms of (A.17) and (A.28), we find for the real corrections (A.29) We emphasise that this contribution corresponds to the three-particle cuts displayed in figure 15, while real Z-boson emission of is not included in (A.29) as the Z boson is considered unstable in our calculation.
The corresponding soft function is given by where the relevant QCD results can be found in (A.18) and (A. 19).
The O(αα 2 s ) contribution to asymmetric heavy-quark production arising from the flavour excitation process can be obtained from (A.29) by crossing. Explicitly, we have (A.31) where the expression for QCD term has already been given in (A.21).

A.6 O(α 2 α s ) contributions
We finally consider the corrections of O(α 2 α s ), which correspond to QCD corrections to the contributions of O(α 2 ) provided in (A.12). Due to the colour structure of these corrections, they can be separated into those to either the massive final-state quarks or the massless initial-state quarks. Example diagrams are depicted in figure 16.
The relevant results for the corrections associated to final-state radiation (FSR) have been provided in [81], and may be written as dσ qq,A  Figure 16. Possible two-and three-particle cuts that contribute to asymmetric heavy-quark production at O(α 2 α s ). Graphs with photon exchange in the s-channel also give a contribution but are not explicitly shown. See text for additional details. heavy quark by where F soft O(α 2 αs) = α s C F 2π 2 + 1 + β 2 β ln 1 − β 1 + β 1 + 2 ln 2E cut µ (A.35) with E cut denoting the upper limit on the gluon energy. For the real emission corrections from the heavy quark lines, the differential cross section reads dσ qq,A where the kinematic function f Q is defined as The QCD corrections to the massless initial-state quark lines contain soft and/or collinear divergences. In this case we have chosen to provide an implementation of the O(α 2 α s ) corrections using the technique of phase-space slicing [41], and performed a crosscheck using dipole subtraction [42]. Rather than repeating the necessary details of both techniques, we instead refer the reader to section D of [41] for phase-space slicing, and appendix D of [42] for dipole subtraction. In the latter case, the relevant formula for the virtual correction is given in (D.9), while the general formula for the operator insertions are collected in appendix C. These relative O(α s ) corrections can be applied to our result for the Born-level cross section of O(α 2 ) provided in (A.12).
In addition to this, we provide the result for the real emission contributions to the differential cross section. They read dσ qq,A

JHEP03(2019)166
The cross section for the qg-initiated contributions can be obtained from crossing, and by additionally adjusting the colour averaging over the initial state.

A.7 Slicing parameter dependence
To conclude this appendix, we perform a numerical study of the dependence on the slicing parameter E cut used in the phase-space slicing technique. We do this by computing the asymmetric cross section for b-quark pair production within the LHCb acceptance (2.2) for pp collisions at √ s = 13 TeV. The calculation is performed with the input parameters given in section 2.4, with factorisation and renormalisation scales set to µ F = µ R = m bb . The invariant mass of the b-jet pair is furthermore restricted to m bb ∈ [75,105] GeV. In the left panel of figure 17 the contribution to the asymmetric cross section is shown for both the twoand three-body O(α 3 s ) contributions as well as their sum. This type of correction has been discussed in appendix A.2. In the lower panel, the y-axis is zoomed into the region around the sum of the O(α 3 s ) contributions, where the shown uncertainties are due to the accuracy of the numerical integration. A similar study is presented for the O(α 2 α s ) FSR corrections discussed in appendix A.6. The corresponding results are given on the right-hand side in figure 17. The numerical results of this work employ the choice E cut = 5 · 10 −5 GeV, and the NLO coefficients are obtained with a relative precision of around 1% to 2% in this case.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.