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 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.


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, fixedorder calculations in perturbative QCD based on collinear factorization miss the effect of large energy logarithms, entera e-mail: papa@cs.infn.it ing 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][2][3][4] 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) [5,6] and for any fixed (not growing with energy) momentum transfer t and any possible two-gluon color state in the t-channel [7][8][9][10][11].
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 good agreement with the HERA data [12,13], 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 inves-tigate more exclusive observables, which can, in principle, only be 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 [14] 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 [15][16][17][18][19]. 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 [20][21][22], 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 [23][24][25][26][27][28][29][30][31][32][33][34][35] 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 [36,37] 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 [38,39]. Recently [40], 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. [41,42] and recently investigated with NLA BFKL accuracy in Ref. [43]. Its natural extension, the four-jet production process, has been proposed in Ref. [44] and studied in Ref. [45].
In a recent paper [46] we suggested a novel possibility, i.e. the inclusive dihadron production 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,

Fig. 1
Inclusive dihadron production process in multi-Regge kinematics 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 [47] within NLA: it was shown there that ultraviolet divergences are taken care of by the renormalization of the QCD coupling, soft and virtual infrared divergences cancel each other, whereas the surviving infrared collinear ones are compensated by the collinear counterterms related with the renormalization of PDFs for the initial proton and parton fragmentation functions (FFs) describing the detected hadron in the final state within collinear factorization. 1 Hence, infrared-safe NLA predictions for observables related with this process are amenable, thus making this process an additional clear channel to test the BFKL dynamics at the LHC. The reaction (1) can be considered complementary to Mueller-Navelet jet production, since hadrons 1 The identified hadron production vertex in the NLA was found within the shockwave approach (or Color Glass Condensate effective theory) in [48]. 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 also [49]. Unfortunately, the comparison between the results of [48] and those of [47] is not simple and straightforward, since the distribution of radiative corrections between the kernel and the impact factor is different in the shock-wave and the BFKL frameworks. Non-trivial 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.
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 analyzed by CMS [50,51] and ATLAS [52] 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. [46] 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 NLA corrections of the hadron vertices. It is well known that the inclusion of NLA terms has a large impact on the theory predictions for the Mueller-Navelet jet cross sections and the jet azimuthal angle distributions. Similar features are expected also for our case of inclusive dihadron production. As for Mueller-Navelet jets, the inclusion of full NLA effects in the process (1) is very important in order to have a control on the accuracy of predictions, in particular on effects related with the choice of the renormalization scale μ R and the factorization scale μ F .
The main aim of this paper is to extend and complete the analysis done in Ref. [46] 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 [53]. In BLM the renormalization scale ambiguity is eliminated by absorbing the non-conformal terms, proportional to the QCD β 0 -function, into the running coupling. Such an approach was successfully used, first in [28], for a satisfactory description of the LHC data on the azimuthal correlations of Mueller-Navelet jets [40], obtained by the CMS Collaboration.
As for the factorization scale, we chose either to fix it equal to the renormalization scale, μ F = μ R = μ BLM R , or we used 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 Sect. 1 we present the theoretical framework and sketch the derivation of our predictions; in Sect. 2 we show and discuss the results of our numerical analysis; finally, in Sect. 3, we draw our conclusions and give some outlook.

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 , the space part of the four-vector p 1 here being taken positive.
In QCD collinear factorization the cross section of the process (1) reads dσ dα 1 dα 2 d 2 k 1 d 2 k 2 = a,b=q,q,g where the a, b indices specify the parton types (quarks q = u, d, s, c, b; antiquarksq =ū,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 centerof-mass energy of the parton-parton collision subprocess.
In the BFKL approach the cross section can be presented (see Ref. [26] for the details of the derivation) as the Fourier sum of the azimuthal coefficients C n , thus dσ dy 1 dy 2 d| k 1 | d| k 2 |dφ 1 dφ 2 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 Hereᾱ is the first coefficient of the QCD β-function, where n f is the number of active flavors. We have being the leading-order (LO) BFKL characteristic function, c 1,2 (n, ν) are the LO impact factors in the ν-representation, which 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) and 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. [47]. It is well known [32] 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 well known [54][55][56][57] that in the BLM approach, applied to semi-hard 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 have with T = T β + T conf , T conf = 3 8 2.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 β 0dependent part in the expression for the observable of interest vanish. In [32] 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.

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 [40] of Mueller-Navelet jet production is much larger, k (jet) min = 35 GeV. In our calculations we use the PDF set MSTW 2008 NLO [58] with two different NLO parameterizations for hadron FFs: AKK [59] and HKNS [60] (see Sect. 3.2 for a related discussion). In the results presented below we sum over the production of charged light hadrons: 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, μ N = | k 1 || k 2 |, so that m R = μ BLM R /μ N , and look for the values of m R such that Eq. (15) is satisfied.
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. [46]; in the case (μ F ) 1,2 = | k 1,2 | they turn out to be generally lower than in the previous case (see Fig. 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 out 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 [32]): 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 [61] and its expression is given, e.g., in Eq. (23) of [26], whereasc 1,2 are the NLA parts of the hadron vertices [47]. 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.

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: MSTW 2008 [58], MMHT 2014 [62] (which is the successor of the MSTW 2008 one), and CT 2014 [63], and convolving them with the three following NLO FF routines: AKK [59], DSS [64,65], and HNKS [60]. Our tests have shown no significant discrepancy when different PDF sets are used in our kinematic range. In view of this result, in the final calculations we selected the MSTW 2008 PDF set (which was successfully used in various analyses of inclusive semihard processes at LHC, including our previous studies of Mueller-Navelet jets), together with the FF interfaces mentioned above. We do not show the results with the DSS routine, since they would be hardly distinguishable from those with the HKNS parametrization.
Specific CERN program libraries [66] were used to evaluate the azimuthal coefficients given in Eq. (17), which requires a complicated eight-dimensional numerical integration (the expressions forc (1) 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 [67] and Psi [68] 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 four-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 [66]. 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 Eqs. (8) and (10)), the one-dimensional integration over the longitudinal momentum fraction ζ in the NLO impact factor correction (see Eq. (10)), and the upper cutoff in the numerical integrations over | k 1,2 | and ν, are negligible with respect to the first one. For this reason the error bars of all predictions presented in this work are just those given by the Dadmul routine.

Discussion
In Fig. 3 we present our results for C 0 in the MS scheme (as implemented in Eq. (5)) for we already specified above the scale settings √ s = 7, 13 TeV, and in the two cases of Y ≤ 4.8 and Y ≤ 9.4. We clearly see that NLA corrections become negative with respect to the LLA prediction when Y grows. Besides, it is interesting to note that the full NLA approach predicts larger values for the cross sections in comparison to the case where only NLA corrections to the BFKL kernel are taken into account. It means that the inclusion into the analysis of the NLA corrections to the hadron vertices makes the predictions for the cross sections somewhat big-ger and partially compensates the large negative effect from the NLA corrections to the BFKL kernel.
The other results we presented below are obtained using BLM in the MOM scheme, as 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. [46] that changing the renormalization scheme produces a non-exponentiated extra factor in Eq. (17) proportional to T conf 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 R , while Y lies on a larger range, i.e. Y ≤ 9.4.
For comparison, we show in Figs. 10 and 11 the results for the same observables with the choice of (μ F ) 1,2 = | k 1,2 |. We clearly see that, in the case of larger rapidity intervals Y and with the natural choice for the factorization scale, the situation is different in comparison to the μ F = μ BLM R choice: the NLA corrections to the cross section C 0 are negative, while the pattern of C 1 /C 0 shows a somewhat unexpected "turn-up" at large Y , and these effects are more pronounced for the lower LHC energy, √ s = 7 TeV. Such a sensitivity to the factorization scale setting may be an indication of the fact that with the increase of Y values we are moving towards the threshold region, where the energy of detected dihadron system becomes comparable with √ s. In this situation the FFs and PDFs are probed in regions that are close to the endpoints of their definitions, where they exhibit large dependence on the factorization scale. From the physical side, in this kinematics the undetected hard-gluon radiation is getting  5)) at natural scales for μ R and μ F , √ s = 7, 13 TeV, and in the two cases of Y ≤ 4.8 and Y ≤ 9.4. Here and in the following figure captions "LLA" means pure leading logarithmic approximation, "NLA kernel" means inclusion of the NLA corrections from the kernel only, "NLA" stands for full inclusion of NLA corrections, i.e. both from the kernel and the hadron vertices restricted and only radiation of soft gluons is allowed. Softgluon radiation cannot change the kinematics of the hard subprocess, therefore one expects restoration of the correlation of the detected dihadrons in the relative azimuthal angle when we approach the threshold region. It is well known that in this situation large threshold double logarithms appear in the perturbative series, and such contributions have to be resummed to all orders. Resummation in the kinematics where both threshold and BFKL logarithms are important is an interest-ing task, but it goes well beyond the scope of the present study. Here we just note that pure BFKL predictions in the region of largest Y become rather sensitive to the choice of the factorization scale.
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- 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 . This originates from the increasing amount of hard undetected parton radiation in the final state allowed by the growth of the partonic subprocess energy.

Conclusions and outlook
In this paper we studied the inclusive dihadron production process at the LHC within the BFKL approach, giving the first complete phenomenological predictions for cross sections and azimuthal correlation momenta in the full NLA approximation. We implemented the exact version of the BLM optimization procedure, which requires the choice of renormalization scale μ R = μ BLM R such that it makes completely vanish the NLA terms proportional to the QCD βfunction. 2 This procedure leads to rather large values of the scale μ BLM R and it allows one to minimize the size of the NLA corrections in our observables. We considered two center-of-mass energies, √ s = 7, 13 TeV, and two different ranges for the rapidity interval between the two hadrons in the final state, Y ≤ 4.8 and Y ≤ 9.4, which are typical for the last CMS analyses. The first rapidity range we investigated, Y ≤ 4.8, may look to be not large enough for the dominance of BFKL dynamics. But we see, however, that in this range there are large NLA BFKL corrections, thus indicating that the BFKL resummation is playing here a non-trivial role. To clarify the issue it would be very interesting to confront our predictions with the results of fixedorder NLO DGLAP calculations. But this would require new numerical analysis in our semi-hard kinematic range, because the existing NLO DGLAP results cover the hard kinematic range for the energies of fixed target experiments; see for instance [69,70].
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, which 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 the 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 may be a physical reason for the observed strong dependence on the factorization scale choice in our pure BFKL approach, and it definitely deserves 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 [71]). The double-parton scattering contribution to the Mueller-Navelet jet production was considered in Refs. [31] and [72], 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 as regards the dynamics of strong interactions in the Regge limit.