Dihadron production at the LHC: full next-to-leading BFKL calculation

The study of the inclusive production of a pair of charged light hadrons (a"dihadron"system) featuring high transverse momenta and well separated in rapidity represents a clear channel for the test of the BFKL dynamics at the Large Hadron Collider (LHC). This process has much in common with the well known Mueller-Navelet jet production; however, hadrons can be detected at much smaller values of the transverse momentum than jets, thus allowing to explore an additional kinematic range, supplementary to the one studied with Mueller-Navelet jets. Furthermore, it makes it possible to constrain not only the parton densities (PDFs) for the initial proton, but also the parton fragmentation functions (FFs) describing the detected hadron in the final state. Here, we present the first full NLA BFKL analysis for cross sections and azimuthal angle correlations for dihadrons produced in the LHC kinematic ranges. We make use of the Brodsky-Lapage-Mackenzie (BLM) optimization method to set the values of the renormalization scale and study the effect of choosing different values for the factorization scale. We also gauge the uncertainty coming from the use of different PDF and FF parametrizations.


I. INTRODUCTION
Semi-hard processes in the large center-of-mass energy limit represent a unique arena to test strong interactions in kinematic regimes so far unexplored, the high luminosity and the record energies of hadronic processes at the Large Hadron Collider (LHC) providing with a wealth of useful data. In the kinematical regime s ≫ |t|, known as Regge limit, fixed-order calculations in perturbative QCD based on collinear factorization miss the effect of large energy logarithms, entering the perturbative series with a power increasing with the order and thus compensating the smallness of the coupling α s . The Balitsky-Fadin-Kuraev-Lipatov (BFKL) approach [1] serves as the most powerful tool to perform the all-order resummation of these large energy logarithms both in the leading approximation (LLA), which means all terms proportional to (α s ln(s)) n , and the next-to-leading approximation (NLA), which means all terms proportional to α s (α s ln(s)) n . In the BFKL formalism, it is possible to express the cross section of an LHC process falling in the domain of perturbative QCD as the convolution between two impact factors, which describe the transition from each colliding proton to the respective final state object, and a process-independent Green's function. The BFKL Green's function obeys an integral equation, whose kernel is known at the next-to-leading order (NLO) both for forward scattering (i.e. for t = 0 and color singlet in the t-channel) [2,3] and for any fixed (not growing with energy) momentum transfer t and any possible two-gluon color state in the t-channel [4][5][6].
The too low √ s, together with small rapidity intervals among the tagged objects in the final state, had been so far the weakness point of the search for BFKL effects. Furthermore, too inclusive observables were considered. A striking example is the growth of the hadron structure functions at small Bjorken-x values in Deep Inelastic Scattering (DIS). Although NLA BFKL predictions for the structure function F 2,L have shown a good agreement with the HERA data [7,8], other approaches can fit these data. The LHC record energy, together with the good resolution in azimuthal angles of the particle detectors, can address these issues: on one side larger rapidity intervals in the final state are reachable, allowing us to study a kinematic regime where it is possible to disentangle the BFKL dynamics from other resummations; on the other side there is enough statistics to define and investigate more exclusive observables, which can, in principle, be only described by the BFKL framework.
With this aim, the production of two jets featuring transverse momenta much larger than Λ 2 QCD and well separated in rapidity, known as Mueller-Navelet jets, was proposed [9] as a tool to investigate semi-hard parton scatterings at a hadron collider. This reaction represents a unique venue where two main resummations, collinear and BFKL ones, play their role at the same time in the context of perturbative QCD. On one hand, the rapidity ranges in the final state are large enough to let the NLA BFKL resummation of the energy logarithms come into play. The process-dependent part of the information needed to build up the cross section is encoded in the impact factors (the so-called "jet vertices"), which are known up to NLO [10][11][12][13][14]. On the other hand, the jet vertex can be expressed, within collinear factorization at the leading twist, as the convolution of the parton distribution function (PDF) of the colliding proton, obeying the standard DGLAP evolution [15], with the hard process describing the transition from the parton emitted by the proton to the forward jet in the final state.
A large number of numerical analyses [16][17][18][19][20][21][22][23][24][25][26][27][28] has appeared so far, which have been devoted to NLA BFKL predictions for the Mueller-Navelet jet production process. All these studies are involved in calculating cross sections and azimuthal angle correlations [29,30] between the two measured jets, i.e. average values of cos (nφ), where n is an integer and φ is the angle in the azimuthal plane between the direction of one jet and the direction opposite to the other jet, and ratios of two such cosines [31,32]. Recently [33], the CMS Collaboration presented the first measurements of the azimuthal correlation of the Mueller-Navelet jets at √ s = 7 TeV at LHC. Further experimental studies of the Mueller-Navelet jets at higher LHC energies and larger rapidity intervals, including also the effects of using asymmetrical cuts for the jet transverse momenta, are expected.
In order to uncover the dynamical mechanisms behind partonic interactions in the Regge limit, new observables, sensitive to the BFKL dynamics and less inclusive than the Mueller-Navelet ones, need to be proposed and considered in the next LHC analyses. An interesting option, the detection of three jets, well separated in rapidity from each other, has been proposed in Refs. [34,35] and recently investigated with NLA BFKL accuracy in Ref. [36].
Its natural extension, the four-jet production process, has been proposed in Ref. [37] and studied in Ref. [38].
In a recent paper [39] we suggested a novel possibility, i.e. the inclusive dihadron pro- when the two charged light hadrons: π ± , K ± , p,p with high transverse momenta and separated by a large interval of rapidity, together with an undetected hadronic system X, are produced in the final state (see Fig. 1 for a schematic view).
This process is similar to the Mueller-Navelet jet production and shares with it the underlying theoretical framework, the only obvious difference lying in the vertices describing the dynamics in the proton fragmentation region: instead of the proton-to-jet vertex, the vertex for the proton to identified hadron transition is needed. Such a vertex was considered in [40]  The identified hadron production vertex in the NLA was found within the shockwave approach (or Color Glass Condensate effective theory) in [41]. It was used there to study the single inclusive particle production at forward rapidities in proton-nucleus collisions; for recent developments of this line of research, see namics at the LHC. The reaction (1) can be considered complementary to Mueller-Navelet jet production, since hadrons can be detected at the LHC at much smaller values of the transverse momentum than jets, thus giving access to a kinematic range outside the reach of the Mueller-Navelet channel.
Note that the inclusive dihadron production was analysed by CMS [43,44] and AT-LAS [45] Collaborations at different LHC energies. The focus was put on two-particle azimuthal angle and rapidity correlations for charged hadrons at low and medium transverse momenta. Here we suggest to analyze this reaction in the region of larger transverse momenta, where data could be confronted with perturbative QCD predictions.
In Ref. [39] we gave the first predictions for cross sections and azimuthal angle correlations of the process (1) in an approximated way, since we neglected, for the sake of simplicity, the The main aim of this paper is to extend and complete the analysis done in Ref. [39] by giving full NLA predictions at √ s = 7, 13 TeV and considering two distinct ranges for the rapidity interval Y between the two hadrons: Y ≤ 4.8 and Y ≤ 9.4. It is well known that, even after the account of the NLA effects, predictions within BFKL resummation still suffer from large ambiguities in the choice of scales. As an idea for the renormalization scale choice setting we adopt the Brodsky-Lepage-Mackenzie (BLM) scheme [46]. In BLM the renormalization scale ambiguity is eliminated by absorbing the non-conformal, proportional to the QCD β 0 -function, terms into the running coupling. Such approach was successfully used, first in [21], for a satisfactory description of the LHC data on the azimuthal correlations also [42]. Unfortunately, the comparison between the results of [41] and those of [40] is not simple and straightforward, since the distribution of radiative corrections between the kernel and the impact factor is different in the shockwave and the BFKL frameworks. Nontrivial kernel and impact factor transformations are required for such a comparison. It certainly deserves a separate study, and the consideration of the process (1) within both the shockwave and the BFKL resummation schemes seems the best possibility to this purpose.
As for the factorization scale, we chose either to fix it equal to the renormalization scale, , or we use a scheme with two separate values of the factorization scale and fix them at the transverse momentum of one or the other of the two detected hadrons, (µ F ) 1,2 = | k 1,2 |, depending on which of the two vertices is considered.
The summary of the paper is as follows: in Section 1 we present the theoretical framework and sketch the derivation of our predictions; in Section 2 we show and discuss the results of our numerical analysis; finally, in Section 3, we draw our conclusions and give some outlook.

II. THEORETICAL FRAMEWORK
The process under investigation (see (1) and Fig. 1) is the inclusive production of a pair of identified hadrons featuring large transverse momenta, k 2 1 ∼ k 2 2 ≫ Λ 2 QCD and separated by a large rapidity interval in high-energy proton-proton collisions. The protons' momenta p 1 and p 2 are taken as Sudakov vectors satisfying p 2 1 = p 2 2 = 0 and 2(p 1 p 2 ) = s, so that the momentum of each hadron can be decomposed as In the center of mass system, the hadrons' longitudinal momentum fractions α 1,2 are connected to the respective rapidities through the relations y 1 = 1 2 ln , and y 2 = 1 2 ln k 2 2 α 2 2 s , so that dy 1 = dα 1 α 1 , dy 2 = − dα 2 α 2 , and Y = y 1 − y 2 = ln α 1 α 2 s | k 1 || k 2 | , here the space part of the four-vector p 1 being taken positive.
In QCD collinear factorization the cross section of the process (1) reads where the a, b indices specify the parton types (quarks q = u, d, s, c, b; antiquarksq = u,d,s,c,b; or gluon g), f a (x, µ F ) denotes the initial proton PDFs; x 1,2 are the longitudinal fractions of the partons involved in the hard subprocess, while µ F is the factorization scale; dσ a,b (ŝ) is the partonic cross section andŝ ≡ x 1 x 2 s is the squared center-of-mass energy of the parton-parton collision subprocess.
In the BFKL approach the cross section can be presented (see Ref. [19] for the details of the derivation) as the Fourier sum of the azimuthal coefficients C n , having so: where φ = φ 1 − φ 2 − π, with φ 1,2 are the two hadrons' azimuthal angles, while y 1,2 and k 1,2 are their rapidities and transverse momenta, respectively. The φ-averaged cross section C 0 and the other coefficients C n =0 are given by dν .
Hereᾱ s (µ R ) ≡ α s (µ R )N c /π, with N c the number of colors is the first coefficient of the QCD β-function, where n f is the number of active flavors.
is the leading-order (LO) BFKL characteristic function, c 1,2 (n, ν) are the LO impact factors in the ν-representation, that are given as an integral in the parton fraction x, containing the PDFs of the gluon and of the different quark/antiquark flavors in the proton, and the FFs of the detected hadron, and while c (1) are the NLO impact factor corrections in the ν-representation. The expressions for them can be derived from the last two lines of Eq. (4.58) in Ref. [40]. It is known [25] that contributions to the NLO impact factors that are proportional to the QCD β 0 -function are universally expressed in terms of the LO impact factors of the considered process, through the function f (ν), defined as follows: It is known [47], that in the BLM approach applied to semihard processes, we need to perform a finite renormalization from the MS to the physical MOM scheme, whose definition is related to the 3-gluon vertex being a key ingredient of the BFKL resummation. So, we with T = T β + T conf , T conf = 3 8 3439 and ξ is the gauge parameter of the MOM scheme, fixed at zero in the following. The optimal scale µ BLM R is the value of µ R that makes the β 0 -dependent part in the expression for the observable of interest vanish. In [25] some of us showed that terms proportional to the QCD β 0 -function are present not only in the NLA BFKL kernel, but also in the expressions for the NLA impact factor. This leads to a non-universality of the BLM scale and to its dependence on the energy of the process.
Finally, the condition for the BLM scale setting was found to be The first term in the r.h.s. of (15) originates from the NLA corrections to the hadron vertices and the second one (proportional to α MOM s ) from the NLA part of the kernel.

III. RESULTS AND DISCUSSION
A. Integration over the final state phase space In order to match the actual LHC kinematical cuts, we integrate the coefficients over the phase space for two final state hadrons, For the integrations over rapidities we consider two distinct ranges: 1. y 1,min = −y 2,max = −2.4, y 1,max = −y 2,min = 2.4, and Y ≤ 4.8, typical for the identified hadron detection at LHC; 2. y 1,min = −y 2,max = −4.7, y 1,max = −y 2,min = 4.7, and Y ≤ 9.4, similar to those used in the CMS Mueller-Navelet jets analysis.
As minimum transverse momenta we choose k 1,min = k 2,min = 5 GeV, which are also realistic values for the LHC. We observe that the minimum transverse momentum in the CMS analysis [33] of Mueller-Navelet jet production is much larger, k (jet) min = 35 GeV. In our calculations we use the PDF set MSTW 2008 NLO [48] with two different NLO parameterizations for hadron FFs: AKK [49] and HKNS [50] (see Section III B for a related discussion). In the results presented below we sum over the production of charged light hadrons: π ± , K ± , p,p.
In order to find the values of the BLM scales, we introduce the ratios of the BLM to the "natural" scale suggested by the kinematic of the process, The values for m R are not affected by the inclusion of the NLO corrections to the impact factor, therefore, for the region Y ≤ 4.8 and for the case µ F = µ BLM R , are exactly the same shown in Fig. 1 of Ref. [39]; in the case (µ F ) 1,2 = | k 1,2 | they turn to be generally lower than in the previous case (see Figure 2 for the summary of all determinations for m R in the region Y ≤ 4.8). In the rapidity region 4.8 < Y ≤ 9.4 we got values for m R similar to those shown in Fig. 2, except for n = 3, where m R turned to be four to five times larger than in the region Y ≤ 4.8.
Then we plug these scales into our expression for the integrated coefficients in the BLM scheme (for the derivation see [25]): The coefficient C 0 gives the total cross sections and the ratios C n /C 0 = cos(nφ) determine the values of the mean cosines, or azimuthal correlations, of the produced hadrons. In Eq. (17),χ(n, ν) is the eigenvalue of NLA BFKL kernel [51] and its expression is given, e.g.
in Eq. (23) of [19], whereasc 1,2 are the NLA parts of the hadron vertices [40]. As anticipated, we give predictions for C n by fixing the factorization scale µ F in two different ways: All calculations are done in the MOM scheme. For comparison, we present results for the φ-averaged cross section C 0 in the MS scheme (as implemented in Eq. (5)) for √ s = 7, 13 TeV and for Y ≤ 4. 8, 9.4. In this case, we choose natural values for µ R , i.e. µ R = µ N = | k 1 || k 2 |, and the option 2., i.e. (µ F ) 1,2 = | k 1,2 | for the factorization scale.

B. Used tools and uncertainty estimation
We performed all numerical calculations in Fortran, choosing a two-loop running coupling setup with α s (M Z ) = 0.11707 and five quark flavors. It is known that potential sources of uncertainty could be due to the particular PDF and FF parametrizations used.
For this reason, we did preliminary tests by using three different NLO PDF sets, expressly: 1,2 contain an additional longitudinal fraction integral in comparison to the formulas for the LLA vertices, given in Eqs. (8) and (9)). Furthermore, slightly modified versions of the Chyp [56] and Psi [57] routines were used to calculate the Gauss hypergeometric function 2 F 1 and the real part of the ψ function, respectively.
The most significant uncertainty comes from the numerical 4-dimensional integration over the two transverse momenta | k 1,2 |, the rapidity y 1 , and over ν. Its effect was directly estimated by Dadmul integration routine [55]. The other three sources of uncertainty, which are respectively: the one-dimensional integration over the parton fraction x needed to perform the convolution between PDFs and FFs in the LO/NLO impact factors (see Eq. (8) and (10) The other results we presented below are obtained using BLM in the MOM scheme, as it is given in Eq. (17). In Figs. 4 and 5 we present our results for C 0 and for several ratios C m /C n at √ s = 13 and 7 TeV, respectively; µ F is set equal to µ BLM R , while Y ≤ 4.8. It is worth to note that in this case the NLA corrections to C 0 are positive, so they increase the value of the φ-averaged cross section at all values of Y . This is the result of the combination of two distinct effects: on one side, we already saw in Ref. [39] that changing the renormalization scheme produces a non-exponentiated extra factor in Eq. (17) proportional to T conf , and that is positive. On the other side, we found that the C gg coefficient in Eq. (10) gives a large and positive contribution to the NLO impact factor. We see also that NLA corrections increase the azimuthal correlations: C 1 /C 0 , C 2 /C 0 , and C 3 /C 0 , while their effect is small with respect to LLA predictions in their ratios, C 2 /C 1 and C 3 /C 2 . The value of C 1 /C 0 for Y ≤ 2.75 in some cases exceeds 1. We consider this as an effect due to the fact that, at very small Y , which corresponds to the small values of partonic subenergiesŝ, we are crossing the applicability limit of the BFKL approach, which systematically neglects any contributions that are suppressed by the powers ofŝ.
For comparison, we show in Figs. 6 and 7 the results for the same observables with the choice of (µ F ) 1,2 = | k 1,2 |. The patterns we have found are very similar to the previous ones, but we see that the effect of having C 1 /C 0 larger than 1 at small Y is reduced. Furthermore, NLA corrections are negative for larger Y values. On the basis of this, we may conclude that, in the Y ≤ 4.8 kinematical regime, the choice of natural scales for µ F stabilizes the results.
In Figs. 8 and 9 we present our results for C 0 and for several ratios C m /C n at √ s = 13 and 7 TeV respectively; µ F is set equal to µ BLM To better assess the factorization scale dependence, we have considered also the case when µ F is varied around its "natural value" k 1 k 2 by a factor r taking values in the range 1/2 to four. In Fig. 12, as a selection of our results, we present the plots for C 0 and C 1 /C 0 at a squared center-of-mass energy of 7 and 13 TeV for the rapidity region Y ≤ 4.8 and the HKNS parametrization of the fragmentation functions.
At the end of this section it is worth to note that the general features of our predictions for dihadron production are rather similar to those obtained earlier for the Mueller-Navelet jet process. Although the BFKL resummation leads to the growth with energy of the partonic subprocess cross sections, the convolution of the latter with the proton PDFs makes the net effect of a decrease with Y of our predictions. This is due to the fact that, at larger values of Y , PDFs are probed effectively at larger values of x, where they fall very fast. For the dihadron azimuthal correlations we predict a decreasing behavior with Y . That originates from the increasing amount of hard undetected parton radiation in the final state allowed by the growth of the partonic subprocess energy.

IV. CONCLUSIONS AND OUTLOOK
In this paper we studied the inclusive dihadron production process at the LHC within the require new numerical analysis in our semihard kinematic range, because the existing NLO DGLAP results cover the hard kinematic range for the energies of fixed target experiments, see for instance [59,60].
As for the hadron's transverse momenta, we imposed the symmetrical lower cutoff: | k 1,2 | ≥ 5 GeV. Considering a region of lower hadron transverse momenta, say | k 1,2 | ≥ 2 GeV, would lead to even larger values of the cross sections. But it should be noted that in our calculation we use the BFKL method together with leading-twist collinear factorization, which means that we are systematically neglecting power-suppressed corrections. Therefore, going to smaller transverse momenta we would enter a region where higher-twist effects must be important.
The general features of our predictions for dihadron production are rather similar to those obtained earlier for the Mueller-Navelet jet process. In particular, we observe that the account of NLA BFKL terms leads to much less azimuthal angle decorrelation with increasing Y in comparison to LLA BFKL calculations. As for the difference between the Mueller-Navelet jet and dihadron production processes, we would mention the fact that, contrary to the jets' case, the full account of NLA terms leads in dihadron production to an increase of our predictions for the cross sections in comparison to the LLA BFKL calculation.
We considered the effect of using different parametrization sets for the PDFs and the FFs, that could potentially give rise to uncertainties which, in principle, are not negligible.
We did some preliminary tests devoted to gauge the effect of using different PDF routines, showing that it leads to no significant difference in the results. Then, we investigated the Y -behavior of our observables by using two different FF parametrizations. Our calculation with the AKK FFs gives bigger cross sections, while the difference between AKK and HKNS is small, since the FFs uncertainties are mostly wiped out in the azimuthal ratios.
We studied the effect of using two different choices for the factorization scale, µ F = µ BLM R and (µ F ) 1,2 = | k 1,2 |, whereas µ R = µ BLM R runs at BLM scales. We see some difference in predictions within these two approaches, especially for larger values of Y and at the smaller value of the energy √ s = 7 TeV. In this region, the kinematic restriction for the undetected hard gluon radiation may start to be important, requiring resummation of threshold double logs together with BFKL logarithms of energy. This issue maybe a physical reason for the observed strong dependence on the factorization scale choice in our pure BFKL approach, and it definitely deserves a further study.
The applicability border for our approach could be established either by comparing our predictions with future data or by confronting it with some other theoretical predictions which do include higher-twist effects. For the last point, one can consider an alternative, higher-twist production mechanism, related with multiparton interactions in QCD (for a review, see [61]). The double-parton scattering contribution to the Mueller-Navelet jet production was considered in the papers [24] and [62], using different approaches. It would be very interesting if similar estimates were done also for the case of dihadron production.
We plan to extend this study by investigating the effect of using asymmetrical cuts for the hadrons' transverse momenta as well as studying less inclusive processes where at least one light charged hadron is always tagged in the final state.
We encourage experimental collaborations to include the study of the dihadron production in the program of future analyses at the LHC, making use of a new suitable channel to improve our knowledge about the dynamics of strong interactions in the Regge limit.

V. ACKNOWLEDGMENTS
We thank G. Safronov and I. Khmelevskoi for stimulating and helpful discussions. This work was supported in part by the RFBR-15-02-05868.