Merging high energy with soft and collinear logarithms using HEJ and PYTHIA

We present a method to combine the all-order treatment of the High Energy Jets exclusive partonic Monte Carlo (HEJ) with the parton shower of Pythia8, while retaining the logarithmic accuracy of both. This procedure enables the generation of fully realistic and hadronised events with HEJ. Furthermore, the combination of the two allorder treatments leads to improvements in the quality of the description of observables, in particular for regions with disparate transverse scales.


Introduction
The analyses of collider data collected at both the Tevatron [1] ( √ s = 1.96 TeV) and the LHC [2][3][4][5][6][7] ( √ s = 7, 8 TeV) indicate that perturbative terms beyond fixed order are required for the description of observables in processes involving at least two jets, in the region of large partonic centre-of-mass energy √ŝ compared to the typical transverse momentum scale p ⊥ , √ŝ /p ⊥ > 5. This corresponds to the jets spanning more than 3 units of rapidity. It is of course well-known and indeed not surprising that the convergence of the perturbative series requires input beyond fixed-order perturbation theory in certain regions of phase space. The dominant and large corrections in this particular region of phase space is the focus of one of the applications of the theory of Balitsky-Fadin-Kuraev-Lipatov (BFKL) [8][9][10][11]. The study of such perturbative effects has received renewed interest, not only because the increasing energy of colliders allows for detailed study of observables in these regions, but also because some measurements specifically concentrate on experimental signatures with large rapidity spans [3][4][5][6][7][12][13][14].
The BFKL formalism can provide a systematic description of large perturbative corrections in two separate kinematic limits. Firstly, in the small-x limit ofŝ/s 1 (where

JHEP09(2018)074
would not contain the correct amount of collinear radiation. Although the algorithm properly prevented the double-counting of such soft emissions, there was no mechanism in place to account for the probability that the parton shower might preferentially have inserted collinear emissions at an earlier stage.
In this paper we therefore present a new method for combining the effects of soft and collinear logarithms with those of the all-order summation of HEJ based on the advances made in the merging of parton showers with fixed-order matrix elements. A crucial feature of our approach is that the exclusive n-parton events generated according to the HEJ all-order matrix elements will be reweighted using properly subtracted collinear Sudakov factors, and moreover the parton shower will be able to insert collinear emissions where it is appropriate to do so. This has been implemented for the interleaved parton shower of PYTHIA8 [37], allowing for the inclusion of multiple partonic interactions as well as the subsequent hadronisation of the event.
The outline of the paper is as follows. The all-order calculation of HEJ is described in section 2. This is followed in section 3 by a brief description of the relevant parts of PYTHIA8 and the tree-level matrix element merging procedures from which our algorithm is inspired. Section 4 will present the method for combining the calculations of PYTHIA8 with HEJ. In section 5 we examine the performance our merged description, firstly by demonstrating the capacity of the new approach to describe jet profiles. Secondly we compare to data for a set of observables that measure additional radiation in inclusive dijet events. We note that in this paper we restrict our focus to pure dijet studies, despite the relevance to Higgs phenomenology. The reasons for this are two-fold. Firstly, the observables of interest have not yet been measured for Higgs plus dijet processes, and secondly it is preferable to test newly developed tools in cleaner environments where there is no expectation of new physics. Nevertheless the method we present should be easily applicable to other processes. Finally we present the conclusions and outlook in section 6.

The High Energy Jets formalism
The framework of High Energy Jets (HEJ) [28][29][30] provides an approximation to the perturbative hard scattering matrix elements for jet production to any order in the strong coupling. The results are exact in the limit of large invariant mass between all particles. The formalism is inspired by the high energy factorisation of matrix elements (as pioneered by BFKL [8][9][10][11]), and obtains a power series inŝ for the square of the scattering matrix elements. Within HEJ, approximations are only applied to the matrix elements. This is different to the framework of BFKL, where numerous kinematic approximations are applied in order to cast the cross section in the form of a two-dimensional integral equation. The highest power inŝ/p 2 t from the square of the matrix element gives the leading-logarithmic contribution (inŝ/p 2 t ) to the cross section. Logarithmic corrections additionally arise from virtual corrections. Recently it was shown that some next-to-leading contributions may be reached within HEJ [38] by including so-called unordered emissions, which have the square JHEP09(2018)074 of the matrix elements suppressed by one power inŝ compared to the leading flavourconfiguration for the same rapidity-ordered momenta. However, in the present study we consider only the leading-logarithmic contributions to the cross section, where only certain Fadin-Kuraev-Lipatov [9] (FKL) partonic configurations contribute.
We will now discuss in more detail the features of HEJ that are relevant to the construction of an algorithm for merging with a parton shower. The all-order perturbative treatment of pp → jj scattering in HEJ starts with an approximation to the tree-level amplitude for the scattering process f 1 f 2 → f 1 g · · · gf 2 , where the final-state particles are listed according to their ordering in rapidity, and f 1 , f 2 can be quarks, antiquarks or gluons. These are the FKL configurations that give rise to the leading contribution to the inclusive n-jet cross section in the Multi-Regge-Kinematic (MRK) limit (see ref. [38] for a recent discussion of the power-suppression of other partonic contributions to the same multi-jet process). The MRK limit can be specified as the limit of large rapidity separations between all particles, for fixed transverse momentum scales: ∀i : y 1 · · · y i−1 y i · · · y n ; p ⊥i ∼ p ⊥ (2.1) It should be noted that the existence of large transverse momentum hierarchies is not compatible with the MRK limit, which will be of importance later. The 2 → n scattering amplitude is approximated at lowest order by the following expression [30]: where S f 1 f 2 →f 1 f 2 2 denotes the square of a pure current-current scattering, K f 1 , K f 2 are flavour-dependent colour-factors (which can depend also on the momentum of the particles of each flavour f 1 , f 2 , see ref. [30] for more details); q i are the momenta of the colouroctets exchanged in the t-channel, and t i = q 2 i . The leading-logarithmic contribution to jet production beyond the first two jets is given by gluon emission from the underlying 2 → 2 process f 1 f 2 → f 1 f 2 , and the effective vertex for gluon emissions takes the form [28]: This form of the effective vertex is fully gauge invariant; the Ward Identity, p j · V = 0 (j = 2, . . . , n − 1) can easily be checked, and is valid for any values for the outgoing momenta p j (that is, not just in the MRK limit).

JHEP09(2018)074
The virtual corrections to the amplitude for each multiplicity are approximated in D = 4 + 2ε dimensions with the Lipatov ansatz [11] for the t-channel gluon propagators (see ref. [28] for more details). This is obtained by the simple replacement in eq. (2.2), where y i are the rapidities of the outgoing partons and is the Regge trajectory which is regulated in D = 4 + 2ε dimensions, in which q 2 i is the Euclidean square of the transverse components of q i . The cancellation of the poles in ε between the real and virtual corrections is organised using a subtraction term, such that the regulated matrix elements used in the all-order summation of HEJ are given by [30]: where and λ is a regularisation parameter describing the extent of the subtraction terms in the real emissions phase space.
Here α s is evaluated using a renormalisation scale µ R , which typically is chosen to reflect the momenta of the final-state partons. Possible choices include half the scalar sum of transverse momenta (H T /2) and the maximum jet transverse momentum (p T max ). Since the matrix elements have been regulated, this allows for a finite numerical approximation to the all-order scattering amplitude to be constructed, and for this to be integrated over all of phase space using a Monte Carlo approach (allowing for the application of arbitrary phase space cuts). Just as for perturbative fixed-order calculations, the parton momenta in eq. (2.6) are interpreted as arising from identifiable partons. An NLO calculation of the production of dijets would deliver the exclusive dijet cross-section to order α 3 s and the inclusive trijet cross-section at the same order in α s . The perturbative result in eq. (2.6) contains real and virtual corrections to any order, and the momenta and multiplicities should all be considered exclusive (to the logarithmic accuracy of HEJ).

JHEP09(2018)074
Indeed, the all-order dijet cross section is simply calculated by explicitly summing the exclusive n-parton cross sections (calculated by numerically integrating the matrix elements squared from eq. (2.6) over all of phase space) over all numbers of gluon emissions from the initial scattering f 1 f 2 → f 1 f 2 . In addition, matching to tree-level matrix elements is performed by reweighting each exclusive m-jet event with the factor: This is just the ratio of the square of the full tree-level matrix element (evaluated using MADGRAPH5 aMC@NLO [39]) to the approximation of this in eq. (2.2), both evaluated on a set of shuffled momenta p new J l ({p i }) derived from the hard jets only. This procedure is summarised in the following formula: where n is the partonic multiplicity of the final state, and the operator O e mj returns one if there are exactly m jets, and zero otherwise. We also define the inclusive dijet operator O 2j = ∞ n=2 O e mj , and require that the extremal partons from HEJ are members of the extremal jets, in order to ensure that the partonic configuration matches the situation for which the HEJ framework was developed. The last line in eq. (2.9) corresponds to the inclusion of parton density functions (PDFs) and the momentum-conserving deltafunctional. Finally we note that while the sum in the first line of eq. (2.9) is over all numbers of final-state partons, 2 ≤ n < ∞, in practice the sum needs to include only a finite number of terms: for finite rapidities and collider energies, the contribution beyond a certain number of gluons is perturbatively suppressed. The upper bound N is chosen such that the results are insensitive to further emissions. This check is performed by simply keeping track of the contribution from each term in the series, and N = 22 is sufficient for this study. Nevertheless, for completeness this choice will enter into the merging algorithm described in section 4.
The matching of HEJ to tree-level accuracy is currently performed up to four jets. The limit on the multiplicity is determined by the time taken to evaluate the full expressions. In addition, the partonic configurations not conforming to the ordering described above are included in HEJ by simply adding the contributions order by order (again using MADGRAPH5 aMC@NLO [39]), but no all-order summation is performed for these non-FKL configurations. In the current study, we will focus on the FKL configurations, since this is where special attention is needed in order to avoid double-counting. This is unlike Examples of a colour flow (left) which contributes in the limit of wide angle, hard radiation, and (right) a configuration which is suppressed in the same limit. In these diagrams, the final-state gluons (on the right of each picture) are ordered according to their rapidity. the challenge addressed by typical fixed-order merging algorithms, because the description in HEJ goes beyond approximating leading-order matrix elements. As previously discussed, application of the Lipatov ansatz through eq. (2.4) is used to sum to all-orders the leadinglogarithmic virtual corrections to the t-channel poles. Although the approaches of HEJ and PYTHIA8 are complementary and calculate different all-order contributions to the perturbative series, they cover overlapping regions of phase space, and the combination of HEJ and PYTHIA8 therefore requires a new merging algorithm. We conclude this overview on HEJ by reiterating that a parton shower framework such as PYTHIA8 [37] is necessary in order to evolve the partonic state of HEJ to the state of hadronisation, primarily by populating the partonic state with further soft and collinear radiation. In order to obtain the logarithmic accuracy of the shower, it should also populate (with the appropriate probability) any region between disparate transverse scales, which might be generated by HEJ. Since the shower, as well as the subsequent string hadronisation, relies on well-defined colour connections between partons, we now briefly discuss the colour connections arising in HEJ.

The colour connections of High Energy Jets
The colour-ordered Parke-Taylor amplitudes [40] for tree-level gg → g · · · g scattering allow for a very neat analysis [41,42] of the dominant colour configurations in the limit of widelyseparated, hard gluons. The conclusion, as presented in references [41,42], is that the leading contribution in the MRK limit is provided by the colour configurations which can be untwisted into two non-crossing ladders that connect the rapidity-ordered gluons. Figure 1 (left) contains an example of a configuration contributing in the MRK limit, and one (right) which is suppressed. The numbering of the final-state partons is assigned according to their rapidity; as drawn their vertical ordering also coincides with their ordering in rapidity.
The colour connections in figure 1 (left) can be summarised as a134b2a. It is possible to arrange the final-state partons such that no colour lines cross without modifying the vertical order of the final-state particles, namely by moving particle 2 to the left side of the same plot. Since the vertical ordering is unchanged, the rungs of the resulting un-crossed JHEP09(2018)074 ladders are also ordered in rapidity. Such manoeuvres are always possible when the order of the particles in the colour connection string between the two initial-state gluons a . . . b and b . . . a also reflects their order in rapidity, as in the case of {134} and {2} in a134b2a.
The colour connections in figure 1 (right) can be summarised as a1324ba; in this case the string {1324} between a and b is not ordered in rapidity. The only manoeuvre which will untangle the colour connections requires flipping the vertical arrangement of particles 2 and 3 such that their vertical ordering is no longer equivalent to their ordering in rapidity. This configuration is therefore suppressed in the MRK limit, because the two un-crossed ladders are not rapidity-ordered.
Furthermore, the study of references [41,42] shows that all the leading configurations each have the same limit in the MRK limit, resulting in a colour factor C A for every final-state gluon. The limit agrees with that predicted by the amplitudes of Fadin-Kuraev-Lipatov (FKL) [9]. When we pass an event from HEJ to PYTHIA8, we choose a colour configuration at random from the set of colour connections which are leading in the MRK limit, and pass the event using an interface conforming to the Les Houches accord [43]. This method is identical to that applied in ref. [35].

PYTHIA8 and CKKW-L
There are several reasons for using PYTHIA8 to handle the collinear resummation rather than ARIADNE as was done in ref. [35]. First of all, the handling of initial-state radiation in ARIADNE is somewhat peculiar [44] and does not quite fit into a conventional resummation scenario. Furthermore, PYTHIA8 has a much more advanced infrastructure for handling matching and merging. Finally, PYTHIA8 has a very advanced model for multiple partonic scattering (based on ref. [45]) which is needed to have a realistic description of the underlying event.

The interleaved shower in PYTHIA8
PYTHIA8 implements a transverse-momentum-ordered shower [36], which includes not only initial-and final-state emissions, but also interleaves these with multiple partonic scatterings. The general philosophy is that emissions (or sub-scatterings) with high transverse momentum should always be performed before those with lower transverse momentum.
As in all parton shower algorithms the ordering is used to ensure that the probability for any emission remains finite, and that the whole shower process is unitary. Even though the splitting functions for an emission diverge for small transverse momenta according to P (k 2 ⊥ , z) ∝ 1/k 2 ⊥ , at each step of the shower the basic splitting probabilities are amended by the probability that no splittings with larger transverse momenta had happened before. The probability that the hardest emission occurs at the scale k 2 ⊥ with an energy splitting z is given by:

JHEP09(2018)074
Here the Sudakov factor ∆(k 2 ⊥max , k 2 ⊥ ) corresponds to the no-emission probability, ensuring that there were no other emissions between the maximum scale k 2 ⊥max and k 2 ⊥ . A lower cutoff, k 2 ⊥cut , is still needed but can be taken very small and still result in probabilities below unity. Formally, eq. (3.1) resums the leading double-logarithms of k 2 ⊥max /k 2 ⊥ in the soft-collinear limit in the leading colour (large N c ) approximation. It should be noted however that many formally subleading contributions, such as momentum conservation, which in practice give rise to large effects are also included.
The no-emission probabilities are fairly easily implemented using the Sudakov veto algorithm [24], and has simple factorisation properties if several different types of emissions are possible, due to the exponential form. The ordering variable, k ⊥ , used in the evolution is not necessarily the actual transverse momentum of an emission in any Lorentz frame, and it is defined slightly differently depending on the class of emission in the interleaved shower. For final-state radiation -FSR -(or time-like splittings) it is defined as k 2 ⊥FSR = z(1 − z)Q 2 , where Q 2 is the invariant mass of the two final-state partons. For initial-state radiation -ISR -(or space-like splittings) we instead have k 2 ⊥ISR = (1 − z)Q 2 , where now Q 2 is the virtuality of the incoming parton entering the hard sub-system after the emission. Finally for multi-parton interactions (MPI ), k 2 ⊥MPI is simply defined as the transverse momentum in the lab system for the 2 → 2 scattering.

Mergingà la CKKW(-L)
The partonic states generated by a parton shower are exclusive; in other words, the probability to produce an n-parton state in the parton shower is approximately given by the exclusive cross section for exactly n partons. This is in contrast to n-parton states generated by a matrix element generator, where the state is exactly given by the inclusive cross section for having at least n partons. The main principle of algorithms that merge matrix elements with parton showers is therefore to take several inclusive samples with different numbers of partons from a matrix element generator and reweight them with no-emission probabilities to make them exclusive. This allows the samples to be safely added and subsequently showered without any double-counting.
The general idea in this paper is to use HEJ as a matrix element generator and add emissions from PYTHIA8 in a consistent way. In doing so we will use ideas from the CKKW-L merging algorithm [46][47][48], but with some important modifications which will be described in section 4. Here we shall review the pertinent features of the CKKW-L method.
Similarly to merging algorithms such as CKKW [49] and MLM [50], the CKKW-L method takes matrix-element-generated states and tries to reconstruct a sequence of emission scales from which the no-emission probabilities are calculated. While some merging procedures use jet clustering algorithms to do this, CKKW-L looks at the partonic states and tries to answer the question "How would my parton shower have generated this state?", and then reconstructs the full kinematics of the corresponding sequence of emissions in the parton shower. Often there is more than one sequence of emissions possible, in which case one sequence is chosen at random with relative weights given by the product of the values of the corresponding splitting functions. The sequence chosen is referred to as the parton shower history and will comprise of a complete set of intermediate states, {S 0 , . . . , S n }

JHEP09(2018)074
(where S 0 is the lowest multiplicity state and S i has i additional partons) and a series of n parton shower emissions. Each emission i is characterised by an ordering scale k 2 ⊥i , a splitting fraction z i , and an azimuthal angle, φ i . This procedure differs from the standard CKKW algorithm, where the intermediate states are not needed and instead only the emission scales are calculated by the k ⊥ jet-clustering algorithm. Formally this difference only affects sub-leading logarithms.
The no-emission probabilities are then calculated by generating trial emissions from each intermediate state in turn, starting at S 0 . The emission generated from S i will have a maximum scale given by k 2 ⊥i . The probability that this emission has a scale above k 2 ⊥i+1 is exactly the no-emission probability ∆ i (k 2 ⊥i , k 2 ⊥i+1 ). Giving the matrix-element-generated state a weight zero if a trial emission from a given state S i has a scale above k 2 ⊥i+1 is therefore equivalent to reweighting the cross section by the no-emission probability: Here k 2 ⊥0 is the maximum possible scale and corresponds to the scale of the Born level process; k 2 ⊥n+1 ≡ k 2 ⊥M is the merging scale which is given by the cut used in the matrix element generator, and is used to isolate the region of soft and collinear divergences where the parton shower is assumed to give a better description.
We can now freely add more partons below the merging scale with our parton shower. For the case that n = N is the maximum multiplicity of the matrix element samples to be merged, trial emissions from S N are not performed and the last factor ∆ N (k 2 ⊥N , k 2 ⊥N +1 ) is omitted. (This is because there is no possibility of double-counting with states of higher multiplicity, n > N .) Consequently the shower is instead started from k 2 ⊥N . In addition to the reweighting of the cross section by the no-emission probability, there is also a reweighting of the value for α s used in the matrix element, typically evaluated at some fixed renormalisation scale µ R characteristic of the Born level process. For a parton shower resummation, however, it can be shown that true collinear logarithms are better reproduced if α s is evaluated at the scale of the individual shower splittings. The states are therefore reweighted by the factor: Additionally there is a reweighting with PDFs, related to the fact that the no-emission probabilities contains PDF ratios, as explained in more detail in [47]. The end result is that the merged event sample will be constructed by exclusive partonic states where the N hardest emissions above the merging scale are given by the full tree-level matrix element, and all softer emissions are given by the shower. So far we have only considered initial-and final-state showers, but to get a realistic description we also need to consider the underlying event. This cannot be described by a tree-level matrix element, but it may be accounted for in an interleaved parton shower using MPI. This means that we also want to incorporate the MPI "emissions" from the shower

JHEP09(2018)074
in the merged event sample. The underlying event may actually contain hard jets, and it is impossible to separate these from the jets given by the matrix element generator. Therefore we cannot blindly include MPI emissions in the CKKW-L no-emission probabilities above, but the procedure is modified [48] as follows.
As before we reconstruct a parton shower history for every matrix element state. We make trial emissions including MPI from each intermediate state S i for i < n, giving the event a weight zero if the emission scale is above k 2 ⊥i+1 . The last state S n is treated separately. When an emission is generated above k 2 ⊥M , if it corresponds to either ISR or FSR we still give the event a weight zero; if however an MPI is generated we will accept the generated state and continue the shower below the emission scale rather than below the merging scale. The end result is thus changed such that the merged event sample will now consist of exclusive partonic states where the N hardest emissions above the merging scale that are not from an MPI are given by the full tree-level matrix element, and all softer emissions are given by the shower.
We note that merging procedures such as CKKW-L are not necessarily unitary, in that the inclusive lowest multiplicity Born level cross section is not preserved, as is the case in parton showers. This is because of the mismatch between the ratio of full matrix element describing the addition of a parton and the splitting function used in the no-emission probabilities. This is in contrast to matching procedures (see e.g. [51]) where it is the matrix element ratios that are exponentiated in the no-emission probabilities.

Modified merging for HEJ
In merging matrix elements with parton showers there are two primary challenges encountered, which we recapitulate here so as to compare the corresponding challenges in merging HEJ with a parton shower. The first challenge is to ensure there is no double-counting between the fixed order matrix elements and the parton shower. In fixed order merging algorithms, this is achieved through the merging scale, which provides a clear partition of phase space. Above the merging scale, the multiplicity of hard jets should not be increased by the parton shower, and the distribution of hard jets should be determined by the fixed order matrix elements. Below the merging scale, soft and collinear radiation from the parton shower is added, smearing the energy of the original hard partons, but leaving the original jets' energies largely unchanged.
We want the merging of HEJ and PYTHIA to obtain the logarithmic accuracy of both. Therefore, the parton shower should not change the jet multiplicity relative to HEJ in the MRK limit (namely, at large rapidities with no transverse momentum hierarchies). The parton shower should however be able to add collinear radiation inside the jet cone. One could envisage using a phase space slicing mechanism such that regions populated by HEJ and parton shower are not allowed to overlap. However, in combining HEJ with a parton shower we are aiming to correctly model the amount of radiation (for example, the multiplicity of jets) in regions of phase space sensitive to both high energy and soft-collinear logarithms, which is hard to achieve with a strict partition. Instead we will allow both formalisms to populate their respective overlapping phase spaces and define a subtraction JHEP09(2018)074 term for the splitting functions and corresponding no-emission probabilities in the shower. Double-counting is then avoided by reducing the probability of producing a certain emission in the shower by the probability that HEJ had already performed that emission.
The second challenge in fixed order merging is to avoid double-counting between the inclusive event samples that are combined, which is resolved by making those event samples exclusive through reweighting with Sudakov factors. The picture for merging HEJ with a parton shower is slightly different, because a given n-parton event generated by HEJ is already exclusive. (This is ensured by the inclusion of virtual corrections at the level of the matrix elements to all orders in eq. (2.6).) It would therefore be inappropriate to naïvely reweight events with the Sudakov factors in eq. (3.2) whose kernels correspond to the full Altarelli-Parisi splitting functions. Instead we have devised a procedure where collinear emissions from the shower are added to states produced by HEJ in a way such that the corresponding collinear Sudakov form factors only change the relative weight of different HEJ multiplicities, and retain the inclusive cross section. It does so by inserting emissions also in early stages of the reconstructed parton shower history to avoid under -counting of collinear emissions due to phase space limitations set by the full generated HEJ state, which was the main drawback of the previous approach in [35].
In section 4.1 we will outline the merging procedure without specifying the particulars of how the division of phase space between HEJ and parton shower is achieved. We simply assume that there exists a consistent way to classify a given emission as being belonging to the either the HEJ or parton shower regimes. We use the jet cone radius as an example of a cut to highlight some features of our algorithm. Interpreting this statement in the language of the parton shower implies that it is possible to define both HEJ and collinear (or subtracted) Sudakov factors. We will then develop these ideas, in particular the definition of and procedure to calculate these subtracted Sudakov factors in section 4.2. Finally, the algorithm will be disclosed in full in section 4.3.

CKKW-L and HEJ
A prescription for dressing HEJ events with collinear radiation may be obtained in a manner analogous to how MPI were added to samples of tree-level events in CKKW-L. To understand how this works, we first note that the MPI algorithm could have been reformulated such that one first does a normal reweighting with the no-emission probabilities excluding MPI in the trial emissions, and then go through the reconstructed states a second time making trial emissions using only MPI. In this second round, starting from S 0 , as soon as an MPI trial emission from S i with a scale k ⊥ > k ⊥i+1 is found one simply replaces the original S n state with the S i state plus the additional generated MPI emission. The shower subsequently evolves from the MPI emission scale k ⊥ .
In the analogous procedure for HEJ, since the states are already exclusive we completely skip the first round of reweighting, and proceed directly to the adding of collinear splittings and MPI. As before this is done by first constructing the parton shower history, however the reconstructed states should only correspond to configurations which HEJ could have generated. We will define such 'HEJ states' more precisely in the next section. This is JHEP09(2018)074 followed by the generation of trial emissions from each reconstructed state which HEJ could not have done, namely the collinear emissions and MPI.
If a trial emission from state S i is generated that has a scale k ⊥ > k ⊥i+1 , the original S n state is replaced by the reconstructed state S i with the additional trial emission. If the original event is replaced, the shower is allowed to evolve freely from the scale k ⊥ . In such a prescription, the N hardest emissions that are neither collinear nor correspond to an MPI are generated by HEJ; everything else is generated by the shower. We also skip the reweighting of α s in eq. (3.3), since as discussed in section 2 the scale used in HEJ has been chosen to be characteristic of the event topology. We will however still use α s (k 2 ⊥ ) in the addition of collinear emissions from the shower.
To illustrate how the algorithm works quantitatively we will in the following assume that there is a clean separation in phase space between the HEJ states and the region where we want to dress the jets with collinear emissions from PYTHIA. For instance, we could imagine a simple phase space cut, where HEJ states are required to have a ∆R between any two partons larger than some value, and the PYTHIA splitting functions are set to zero if they result in such states.
We start by reformulating the n-parton state produced by HEJ within the showerformalism, written as a basic two-parton inclusive cross section multiplied by a series of 'HEJ splittings' with decreasing values of the scale reconstructed by the merging algorithm. The fact that the HEJ states are exclusive means we can write the cross section for an n-parton state produced by HEJ (that is, prior to merging) as: Here P H i (k 2 ⊥ ) is the splitting function for emitting a parton at the scale k 2 ⊥ from the state i according to HEJ, integrated over the energy fraction z; ∆ H i (k 2 ⊥i , k 2 ⊥ ) is the probability that there were no 'HEJ-like' emissions from the state i between the scales k 2 ⊥i and k 2 ⊥ . Finally, dσ 2 is the inclusive differential cross section for the initial two-parton state.
The shower-merging will have to construct all shower-histories, which could have produced a given n-parton state from HEJ. This ends up being the time-consuming step for the high-multiplicity states produced by HEJ. These states can be of much higher multiplicity than the current limit experienced with fixed-order matchings, where the shower-histories are also reconstructed. In order to reduce the complexity of the shower history reconstruction, we trim the parton-content of the high-multiplicity states from HEJ before they are passed to the shower. This is done removing any parton with a transverse momentum smaller than a scale k ⊥M from the event record (and reshuffling the remaining momenta to absorb the transverse momentum thus removed). The effect of introducing the trimming though is that the event record contains no partons with transverse momenta less than k ⊥M . After the trimming, this phase space is therefore left completely for the shower to populate, and the trimming scale is thus the final scale for the last Sudakov factor in eq. (4.1).

JHEP09(2018)074
The inclusion of trimming can speed up the merging significantly; however, it should be emphasised that formally we should consider the limit where k ⊥M → 0. Nevertheless, as the weight of the event is kept unchanged, as long as k ⊥M is smaller than the scale of the jet threshold, any observable based on jet momenta is only weakly dependent on this trimming if at all. We will later (in figure 8) investigate directly the numerical impact of the transverse scale used in the trimming of the event record in passing events from HEJ to PYTHIA, which indeed is found to be insignificant even on the observables which are very sensitive to the jet multiplicities of the events.
Before continuing, some comments may be needed to clarify eq. (4.1): • At this stage we do not need to know anything about P H . The fact that we have a sequence of emissions means that we can describe it as a product of splitting functions accompanied by corresponding no-emission probabilities which are of the form given in eq. (3.1), even if the states were not produced that way by HEJ.
• In rare cases it is not possible to reconstruct an ordered history of shower emission. Such cases are handled by joining two (or more) subsequent steps into one, as described in [48]. Such unordered paths are by definition far the parton shower resummation regions and does not affect the logarithmic accuracy of the procedure.
• In other rare cases, it is not possible to find intermediate states corresponding to HEJ states. Again we treat these by joining several steps into one, so that the trial emissions always come from HEJ-like states.
• The total inclusive dijet cross section is given by σ 2 and is not the basic tree-level 2 → 2 cross section. This is because HEJ includes non-unitary corrections beyond leading order, and these we would like to preserve.
We shall now consider the different possibilities which may arise in the merging procedure described above. In the first case, the original state generated by HEJ is not replaced by one generated by the shower. This will occur only if at each reconstructed state in the history, no trial (non-HEJ-like) emission is generated above the scale k ⊥i+1 of the next reconstructed state. This corresponds to multiplying eq. (4.1) by the following product of no-emission factors: where ∆ C i is a suitably modified (collinear) Sudakov factor, for example, corresponding to the exponentiation of a PYTHIA splitting function with the ∆R cut assumed above. Defining ∆ M i = ∆ H i ∆ C i , such events will contribute to the cross section according to: Furthermore it is clear that we can freely dress these states with full PYTHIA splittings below the merging scale.

JHEP09(2018)074
If instead a trial emission is generated from a reconstructed state m ≤ n at the scale k 2 ⊥C > k 2 ⊥M (and above the scale of the next reconstructed state) the original n-parton state from HEJ will be replaced by the reconstructed m-parton state plus the accepted trial emission. Calculating the contribution to the cross section from such states requires summing and integrating over all possible reconstructed HEJ emissions below k 2 ⊥C , where P C m is the ∆R-truncated PYTHIA splitting function. For n = m there is no integral and we just get the probability that there are no extra emissions. For n = m + 1, we get the probability that there is exactly one extra emission, for n = m + 2 exactly two, etc. The sum of all these must necessarily sum up to unity, and the final result is just the two first lines of eq. (4.4).
In practice there is an upper cut, N , on the parton multiplicity in HEJ. Assuming that the corresponding cross section is inclusive over the last emission, the last integral in the third line of eq. (4.4), becomes where we have used the property of no-emission probabilities that its derivative is simply itself times the splitting function, d . When adding this to the N − 1 contribution, this will explicitly cancel the last no-emission factor there, and we can do the last integral in the N − 1 contribution in the same way, and so on, until we cancel also the last no-emission factor in the m = n contribution.
Adding full PYTHIA shower splittings below k 2 ⊥C , we can now write the exclusive probability that we have exactly n partons above the merging scale as where P P is now the full PYTHIA splitting function (possibly also including MPI) and ∆ P the corresponding no-emission probability. Comparing this with what PYTHIA would give JHEP09(2018)074 on its own, we see that for n partons above the merging scale the m hardest ones will always be produced by HEJ, and if there are partons from PYTHIA above the merging scale the hardest one of these will always be a collinear splitting. We also see that the procedure is unitary, in that the inclusive jet cross section is still given by σ 2 as calculated by HEJ. All we have done is to add (unitary) parton shower emissions and, in some cases where these are harder than the HEJ ones, reclustered the original HEJ state into a lower multiplicity state, and then added the parton shower. We note that the action of multiplying by eq. (4.2) was not present in the algorithm presented in [35] for matching HEJ with ARIADNE. That is to say, the probability that the parton shower might have produced a collinear emission at an earlier stage in the reconstructed history was not taken into account. It was the lack of this step which allowed the inclusion of soft gluons from HEJ that interfered with the ordering of the parton shower and prevented a proper parton shower evolution in the full phase space. Furthermore these collinear emissions which according to the parton shower should have occurred were not inserted; instead such emissions could only be included below the matching scale. We emphasise that in this regard the approach we take here is fundamentally different from how HEJ was matched with ARIADNE.

The subtracted shower
In the previous subsection it was assumed that we could make a simple phase space cut between collinear splittings to be described by PYTHIA and large angle splittings from HEJ. However, the MRK-limit in HEJ does not take into account large logarithms that arise in case we have large transverse momentum hierarchies between (possibly widely separated) jets. Such logarithms are included in the parton shower, and we would like to include them in our merging.
To accomplish this we go beyond a simple phase space cut and use subtracted splitting functions instead. The idea is that where resummation is important the splitting functions are large, and we could naively say that where the PYTHIA splitting function is less than the HEJ one we set it to zero, and vice versa. This would still correspond to a simple phase space cut, albeit more complicated than the ∆R cut assumed above. However, this would be fairly wasteful as we would throw away many of the jets produced by HEJ. Instead we have introduced a procedure where the HEJ splitting function is subtracted from the PYTHIA one.
To do this it is now necessary to obtain an explicit definition of the splitting functions and no emission probabilities for HEJ. Although no such expressions appear explicitly within the HEJ formalism, we note that the Altarelli-Parisi splitting functions may be derived as the soft and collinear limit of a ratio of matrix elements [52]:

JHEP09(2018)074
This is just the normal universal behaviour of matrix elements in the soft and collinear limit. The Altarelli-Parisi splitting functions precisely capture the soft and collinear singularities which must be exponentiated to calculate the leading DGLAP logarithms in the parton shower no emission probabilities. If instead we replace the full matrix elements by the HEJ ones, this will no longer contain any collinear singularities, but only the soft singularities. Such a function is precisely what is needed to define a subtraction term for the parton shower. Of course, we could take the MRK limit of this and retain the same logarithmic accuracy, but by using the full matrix elements we retain more of the HEJ accuracy. Therefore, as in the approach of [35], we define the HEJ splitting function as a ratio of HEJ matrix elements given by eq. (2.6) corresponding to an event before and after the insertion of an emission as generated by the parton shower. Of course as noted in section 2, these matrix elements are only valid for FKL configurations, but there is no restriction upon the kind of configuration which may be generated by the parton shower. We must therefore assert that the following criteria define a 'HEJ state': 1. The most forwards outgoing parton should have the same flavour as the parton incoming along the positive z axis.
2. The most backwards outgoing parton should have the same flavour as the parton incoming along the negative z axis.
3. All other outgoing partons must be gluons. 4. It must be possible to untangle the colour connections into two 'ladders' of rapidityordered partons.
5. The outgoing partons must cluster into at least two jets.
6. Each extremal (most forwards or backwards) parton must be a member of the corresponding extremal jet.
7. Each parton must have a transverse momentum above the merging scale k ⊥M . Criteria 1-3 simply define an FKL configuration; criterion 4 is required since a given set of colour connections is chosen (as described in section 2.2) and this has an impact upon which dipoles arise in the PYTHIA parton shower; criteria 5-6 are kinematic constraints on HEJ events. Finally, although not strictly necessary for the purpose of efficiency events generated by HEJ containing soft emissions below k 2 ⊥M are reclustered (in a manner that does not alter the rapidities of the resulting jets) and this is therefore reflected by the requirement given in condition 7. The reclustering reduces the complexity of constructing all possible shower histories for the state passed to PYTHIA. This reduces the CPU time needed to obtained merged predictions, and explicit tests indicate that the impact of the reclustering is unnoticeable on observables based on hard jets, such as those studied in this paper. Now, for any emission resulting in a configuration not corresponding to a HEJ state we set P H = 0 (because there is nothing to subtract for non-FKL states), and otherwise we can define: (4.9)

JHEP09(2018)074
The factor of 1 2 accounts for the fact that the matrix elements are summed over all possible colour connections, but for each parton shower emission we wish to calculate the splitting function for one of two possible choices (each of which contribute equally in the MRK limit). This expression however only accounts for time-like emissions. For a space-like branching i → jk, where parton j is evolved backwards to parton i with a higher momentum fraction x i = (1/z)x j , we instead define: . (4.10) where the PDFs f i,j should be evaluated at an appropriately chosen factorisation scale µ F . Our effective Sudakov factor ∆ H from the previous section would then simply be the exponentiation of this splitting kernel, however, we will not need to compute this explicitly in the numerical implementation of the algorithm. For completeness we shall also write down the PYTHIA8 splitting functions, evaluated as a function of the evolution variables k 2 ⊥ and z: where P (z) is the appropriate unregulated Altarelli-Parisi splitting function. There is an additional factor of 1/(2π) to average over azimuthal angle, since the matrix elements in eq. (4.9) will be evaluated for a given choice of azimuthal angle for the generated emission.
In the case of a space-like branching we modify eq. (4.11) to include the ratio of PDFs for the branching, as in eq. (4.10).
To illustrate the differences between the HEJ and PYTHIA splitting functions, in figure 2 we show their typical behaviour as a function of (a) the angular distance between the emitted gluon and the nearest parton, ∆R, and (b) the transverse momentum of the emitted parton in the lab frame, p ⊥ . What is shown is the average value of the splitting functions in the first emission in HEJ-generated qQ → qQ events, excluding factors of α S and ratios of PDFs. In figure 2a we average over emissions with p ⊥ > 10 GeV. The discontinuity in the HEJ splitting function at the jet radius is an artefact of the regularisation procedure used for emissions inside the jet cone of extremal jets [30]. This is in any case the phase space region we want to populate with collinear shower splittings.
For the p ⊥ plot we average of all emissions with ∆R > 0.6. We clearly see that for small p ⊥ the PYTHIA splitting function exceeds the HEJ one (and also the analytic MRKlimit splitting shown for comparison). This is the region of large transverse momentum hierarchies, where the MRK approximation fails to properly resum the corresponding logarithms. Such large logarithms are present and are resummed by PYTHIA, and we would therefore like to add such splittings even if they are far away from the collinear region.
From this we see that it makes sense to use a simple phase space cut based on the relative sizes of the splitting functions. However, instead we can go one step further and in regions where P P > P H we subtract the HEJ splitting function from the PYTHIA one. There we define the subtracted PYTHIA splitting function as P S (k 2 ⊥ , z) = max P P (k 2 ⊥ , z) − P H (k 2 ⊥ , z), 0 . Here the arguments of P H are intended to be schematic. This is intended to denote that having generated an emission with corresponding evolution variables k 2 ⊥ and z, and having inserted this into the event with an appropriate recoil strategy, the matrix element containing n + 1 partons should be evaluated with the resulting set of n + 1 final-state (recoiled) momenta. With this notation we can now define the subtracted Sudakov factor: It should be clear that we then want ∆ M = ∆ S ∆ H , and in order that such a Sudakov factor might be employed during a trial shower, it is sufficient to generate emissions according to the full PYTHIA splitting function P P (k 2 ⊥ , z), but veto emissions with probability (4.14) in accordance with the Sudakov veto algorithm.
Armed with this we can go through the steps in section 4.1 again and arrive at exactly the same formulae except with P C and ∆ C replaced by P S and ∆ S . The net result is that in phase space regions where P H > P P , where we believe HEJ is doing a good job, we never add any PYTHIA splittings, while in the complementary region emissions are added in proportion to the subtracted splitting function so that in total they will correspond to populating that region only with PYTHIA splittings.

The merging algorithm
For completeness we now explicitly disclose the full algorithm for merging HEJ with PYTHIA as follows: 1. Generate samples of n-parton HEJ states with n ≤ N . Recluster any partons that have momenta above k ⊥M in such a way that the rapidities of the resulting jets is unchanged.
2. For each n-parton state from HEJ (2 < n ≤ N ), reconstruct all possible PYTHIA shower histories where each clustering has the reconstructed scale k 2 ⊥i , and set k 2 ⊥n+1 = k 2 ⊥M . If n = 2 calculate the scale k 2 ⊥2 and continue to step 3, and otherwise continue as follows.
(a) Throw away all histories that do not correspond to a sequence of HEJ states.
(b) If there is at least one history that is correctly ordered in k 2 ⊥i , throw away every other history.
(c) Give each history that is left a weight proportional to the HEJ matrix element squared for the lowest multiplicity (HEJ) state, times the product of PYTHIA splitting functions for the sequence of emissions that gives the original n-parton state. Pick a history at random according to its relative weight.
(d) Starting from the most clustered state in the history, make a trial emission from each intermediate state in the selected history starting from k 2 ⊥i . i. If the emission scale is below the reconstructed scale of the next state in the history, k 2 ⊥i+1 , continue to the next state in the history. If this is the original event we started with, continue to step 3.
ii. If the emission scale is above the reconstructed scale of the next state in the history, but has produced a HEJ state, veto the emission with probability P H /P P . If the emission is vetoed, generate a new trial emission starting from the current emission scale, and return to 2(d)i. If the emission is not vetoed replace the original event with this state and continue to step 3. iii. If the emission scale is above the reconstructed scale and has not produced a HEJ state, we substitute the original event with this state and continue to step 3.

(a)
If in the previous step we replaced the original event with one that could not have been produced by HEJ, continue the shower from the emission scale of the new state without restriction.
(b) If this is the original event and we have n < N start the shower from the reconstructed scale k 2 ⊥n and check the first emission. If it gives a new HEJ state, discard the emission with probability P H /P P and continue generating the first emission starting from the scale k 2 ⊥n . Once a first emission is accepted, the shower continues from the emission scale, radiating freely. 4. Once the parton shower has evolved below the cut-off scale, hadronise the event.

JHEP09(2018)074
This method represents one of the possibilities for merging HEJ with a parton shower. In particular, it retains the dijet cross section and logarithmic accuracy of HEJ: indeed, each event configuration and weight is first generated by HEJ, and all of phase space is thus covered. In the MRK limits of similar transverse scales for all emissions, the Sudakov factors introduced in equation (4.3) all evaluate to unity, since the scale used in the evolution of PYTHIA in the MRK limit tends to the transverse scale of the lab frame. Since the MRK phase space is populated, and the matrix elements are unchanged in this limit, the merging maintains the logarithmic accuracy of HEJ.
The logarithmic accuracy of PYTHIA is ensured since the full allowed phase space in PYTHIA is covered, and the appropriate Sudakov factors between emissions are applied in the shower evolution, with the possibility of generating further emissions from the shower evolution. Such emissions are then vetoed with a probability that said emissions were also generated from HEJ, such that double-counting is avoided.
For completeness we here mention three potential issues with the algorithm. One limitation of the proposed method is that only the hardest emissions (as ordered by PYTHIA) will be merged to HEJ, which is not itself ordered in hardness: it is possible for the parton shower to modify a state classified as non-FKL (according to the momenta above the merging scale) to a FKL state (accounted for in HEJ) through an emission, and such emissions will not have their splitting kernels subtracted. However, the non-FKL configurations account for a logarithmically suppressed part of the cross section, quickly diminishing with increasing rapidity [38]. Furthermore, future accounting for next-to-leading logarithmic contributions in HEJ will decrease further the significance of the parton shower changing non-HEJ to HEJ states.
In addition, we reiterate that the method we are presenting is currently only applicable to FKL configurations. The impact of non-FKL corrections on the observables presented in this study is relatively small and within the indicated scale variations of the FKL results. To include non-FKL configurations, it would be necessary to extend the definition of what constitutes a HEJ state, and ensure that the appropriate tree-level matrix elements are used when calculating the veto probability for non-FKL states. This is so that no problems arise from double-counting. Primarily such changes would affect what states may be included in the parton shower history, and which states may be inserted by the parton shower.
Finally we note that the factor one half in eq. (4.9) is based on the fact that the colour flows which have a leading logarithmic contribution will contribute equally to the coloursummed matrix element squared in the MRK limit. This means that there will always be just two possible colour flows for inserting a gluon in the exchange, and they will have the same leading kinematic term in the MRK limit. While it is relevant to take into account the different kinematic contribution from each possible colour connection when matching full matrix elements to the parton shower, the fact that the collinear divergences are absent from the formalism of HEJ means that the kinematic contributions from different colour connections differ far less than in the full theory. While it would be possible to account for the colour flow dependence in the contribution from the sub-leading (and non-divergent) terms introduced in HEJ compared to BFKL, we choose in this study to assign an equal weight to the each of the possible colour flows, just as will be the case in the MRK limit. This allows for the simple attribution of 1 2 in eq. (4.9). HEJ partons HEJ+Pythia particles Figure 3. A Lego-plot of the momenta of partons arising from a single event from HEJ (blue) and average of the momenta arising from particles from that event showered 10,000 times with PYTHIA8 using the merging of HEJ+PYTHIA. For this event configuration from HEJ, which contains partons of similar transverse momenta, the effect of the showering is mostly to distribute physical particles around the partons of HEJ.

Results
In this section, we will present the results of the formalism developed, and contrast it with experimental data. We start however by qualitatively examining the effect of the parton shower and hadronisation on a specific partonic event from HEJ. This is presented in figure 3, which shows a LEGO-plot of the average transverse momentum deposit (greater than 100 MeV) in bins of 0.2 × 0.2 units of rapidity and azimuthal angle. A single HEJ event with 5 partons (and with fixed colour connections), of which 4 are sufficiently hard to form individual jets of p T > 30 GeV, is shown in blue. The average result of passing this event 10,000 times to PYTHIA8 is shown in red.
The effect of the shower on average is to spread out the momentum of each HEJ parton over an area with radius R ∼ 1 around that parton. Indeed, for events similar to the chosen one, the effect of the shower seems to be limited to filling the jet cones, and in section 5.1 we study in more detail the accuracy with which the jet cones are filled. In section 5.2 we will study multi-jet observables, some of which probe large hierarchies in transverse momentum, and in such regions PYTHIA8 can additionally supplement the jets produced by HEJ.
We will mainly look at LHC analyses especially designed to probe effects of high energy logarithms. It should be noted however, that these have so far employed a relatively soft definition of hadronic jets, typically requiring a transverse momentum of less than 40 GeV. This results in broad shower profiles, where the description of the spill-over outside the jet JHEP09(2018)074 cones is necessary for accurate results. Furthermore, it was noted in ref. [53] that these analyses often use cuts that also enhance soft and collinear logarithms. In section 5.2 we therefore propose to use a slightly harder threshold for jets to reduce the dependence on shower and hadronisation effects, and crucially, investigate the full rapidity range of the hard event rather than just the region in-between the two hardest jets. This allows for a much cleaner probe of the high energy logarithms.
We note that there exist many parameters in PYTHIA8 that control non-perturbative effects, and which are fixed by tuning to measurements of certain soft observables. We investigate an example of such an observable in the next section. As we will shortly see, the combination of HEJ and PYTHIA8 obtains a very similar description to PYTHIA8 alone. Therefore, in the results that follow we do not retune PYTHIA8 for use with HEJ, even if this might further improve the agreement of HEJ+PYTHIA with data; instead we use the default Monash 2013 tune [54] for both PYTHIA8 and HEJ+PYTHIA.

The description of the profile of jets
The jet profiles were measured at the LHC in early 7 TeV runs, for example by ATLAS in ref. [55], accepting events with just one primary vertex (no pile-up) and at least one jet with transverse momentum p ⊥ > 30 GeV and rapidity |y| < 2.8. For such events, the differential jet profile ρ(r) as a function of the distance r = ∆y 2 + ∆φ 2 to the jet axis is defined as the average fraction of the jet transverse momentum in an annulus between r − ∆r/2 and r + ∆r/2 around the jet axis in the (y, φ)-plane. As such, ρ(r) = 1 ∆r 1 N jets jets p ⊥ (r − ∆r/2, r + ∆r/2) p ⊥ (0, R) , where p ⊥ (r 1 , r 2 ) is the summed p ⊥ in the annulus of the two circles of radii r 1 and r 2 , and N jets is the number of jets. The measurement of ref. [55] used ∆r = 0.1 and the anti-k T jet algorithm [56] with R = 0.6. This analysis is implemented in recent versions of Rivet [57], which we use to analyse generated events (both here and in section 5.2). The HEJ formalism captures the logarithms associated with wide-angle emissions, but not those associated with the collinear emissions. HEJ is thus not expected to fill the jet cones with radiation, and it is expected that the results of HEJ +PYTHIA8 for the jet shapes is similar to those of pure PYTHIA8 (since the merging procedure should produce results similar to PYTHIA8 in regions where HEJ does not radiate). In figure 4 we compare the predictions of PYTHIA8 and HEJ+PYTHIA for ρ(r) in slices of jet transverse momentum to data [55]. While HEJ alone would primarily have filled just the first bin in each distribution, HEJ+PYTHIA gives the same very good description of the jet shapes as the parton shower of PYTHIA8. The merging procedure has therefore performed a perfect job of populating the jet areas (through collinear emissions), which are mostly empty in the pure partonic description of HEJ-and has of course furthermore fully hadronised the partonic states. The ability to describe this observable represents an improvement relative to the matching of HEJ + ARIADNE.   Figure 4. The data and predictions for the differential jet profile as defined in eq. (5.1). The parton shower of PYTHIA8 gives a very good description of data, which is inherited by HEJ+PYTHIA.

Impact on multi-jet observables
In this section we will investigate the impact of the merging on observables which depend only on the identified hard jets of the event. We shall make comparisons between pure HEJ, PYTHIA8 and HEJ+PYTHIA. Events in HEJ (both with and without showering) were generated with the PDF set CT14nlo [58,59], and the renormalisation and factorisation scales were taken to be the maximum jet transverse momentum µ R = µ F = p T max . In the case of pure HEJ the scale uncertainties were estimated by varying µ R and µ F independently between twice and half the central scale choice, and these uncertainties will be denoted as a band around the central predictions. The vertical lines indicate the statistical uncertainty on the results. As before, PYTHIA8 predictions were generated using the Monash 2013 tune. The expectation is that there should be little impact on the results of HEJ in phase space regions where the jets have similar transverse momenta but are widely separated in rapidity. This is the region where HEJ should already control the dominant logarithms to all orders. The parton shower should therefore not introduce sizeable corrections. On the other hand, as mentioned in section 2.1 regions with large disparate transverse scales are not encompassed by the kinematic assumptions of the HEJ formalism, and should therefore receive additional hard emissions from the parton shower. We will first consider two ATLAS analyses [2,3] that measure the amount of additional radiation in inclusive dijet events. Dijet systems are of course simple at the Born level, characterised by two jets of equal transverse momenta that are back-to-back in the azimuthal plane; however, this simple topology is in general spoiled by radiative corrections. The analyses in question both require the existence of a dijet pair above some transverse momentum cut, defining the tagging jets; in what follows the tagging jets are identified as the two hardest (leading) jets in the event. The number of jets in the rapidity interval between the tagging jets, each having a transverse momentum above a given veto scale Q 0 , is then measured. This allows the definition of two observables: first, the gap fraction, and secondly the average number of jets in the rapidity interval N jet . Events having no jets above the veto scale in the rapidity interval between the tagging jets are classified as gap events. The gap fraction as defined by ATLAS [2,3] is then simply the ratio of the contribution to the cross section from these gap events to the inclusive dijet cross section.
We start with the ATLAS analysis presented in ref. [2], in which jets were defined using the anti-k T jet algorithm with R = 0.6 and having rapidity |y j | < 4.4. In figure 5 we show a plot of gap fraction as a function of the rapidity interval between the tagging jets |∆y|, where the veto scale was taken to be Q 0 = 20 GeV. This is shown in bins of the average transverse momentum of the tagging dijets p T , from 70 GeV-90 GeV to 240 GeV-270 GeV. By construction, the gap fraction will be 1 at |∆y| = 0, since the phase space where a third jet would be counted is vanishing (since only jets in-between the two hardest jets are counted, and the rapidity difference between the two hardest jets is zero). The variation between the predictions is small, and discernible only for the bins with the largest p T . Here there is a large hierarchy between the scale of the tagging jets and the scale of any additional jets (which are characterised by the veto scale). It is therefore not surprising that HEJ predicts too few additional jets in this region instead requiring the DGLAP resummation of the parton shower. Moreover the combination of HEJ+PYTHIA results in a description that at least as good as, or better than, PYTHIA8 or HEJ individually. The average number of hard jets is a potentially better discriminant between the predictions than the gap fraction, simply because the average number of jets has a larger range of variation. We now consider this observable as measured by ATLAS in ref. [3], where the hardest and second hardest jets (also defining the tagging jets) were required to have transverse momenta above 60 GeV and 50 GeV respectively. 1 Jets were again defined using the anti-k T algorithm with R = 0.6. In figure 7a the average number of jets in the interval between the tagging jets is shown as a function of the rapidity interval between the tagging jets (with Q 0 = 20 GeV). While the differences in the predictions are again small, we observe that although the data from ATLAS lie within the scale uncertainty band for pure HEJ, the central line for HEJ nevertheless underestimates the number of additional jets. Meanwhile, the predictions of HEJ+PYTHIA are better in line with data, and are of a similar quality to that of PYTHIA8.  [3] for the average number of jets in the rapidity interval between the two tagging jets, as a function of (a) the rapidity interval between the two tagging jets, and (b) the average transverse momentum of the two tagging jets.
In figure 7b the average number of jets in the interval between the tagging jets is instead shown as a function of the average transverse momentum of the tagging jets (with Q 0 = 30 GeV), and where the dijets were required to be separated by at least one unit of rapidity. As the average transverse momentum of the two hardest jets increases to 1 TeV, the number of 30 GeV jets is unsurprisingly no longer well-described without a shower resummation. Indeed, for increasing p T , the predictions of pure HEJ rises from 0.12 additional jets to 0.3, whereas data rises from 0.15 to 0.5. Both PYTHIA8 and HEJ+PYTHIA give a good description of this distribution. For such large ratios of transverse scales, the effect of the shower resummation is large, and therefore the results for HEJ+PYTHIA are outside the scale uncertainty band for pure HEJ.
It should be apparent at this stage that in distributions that probe large differences in transverse momentum such as figure 7b a parton shower is necessary for an accurate description, and therefore the addition of PYTHIA8 to HEJ gives rise to a notable improvement relative to HEJ. Likewise, in distributions that probe large rapidity spans, one might have expected that HEJ (and hence HEJ+PYTHIA) would provide a superior description relative to PYTHIA8. Indeed, it is perhaps surprising that the predictions of HEJ and PYTHIA8 are so similar for the rapidity distributions studied so far. (In fact, we note that in some cases the description of PYTHIA8 is closer to data than the predictions of PYTHIA +POWHEG [61][62][63] which were used in the original analyses [2,3]. This could be an effect of the later tunings of the non-perturbative parameters of PYTHIA, and reiterates the possible benefits of performing similar analysis with much harder jet scales, such that the sensitivity to the tunings of the MPI and non-perturbative effects are reduced). Firstly, the restrictive definition of the chosen observables prevents much variation in their values. Also, as discussed in ref. [53], the softness of the veto scale relative to the tagging dijets' transverse momentum results in event samples that are influenced by both high energy and soft-collinear logarithms, spoiling the applicability of the HEJ formalism.
Simpler event samples were suggested in ref. [53] to disentangle the two sources of logarithmic corrections, together with more inclusive observables that better expose the differences in the description of a fixed-order calculation, a parton shower and HEJ. The analysis considered inclusive dijet events, requiring at least one jet with transverse momentum above 45 GeV, and with all other jets required to have transverse momentum above 35 GeV. Jets are defined using the anti-k T algorithm with R = 0.5 and with rapidities |y j | < 4.7. Comparisons between PYTHIA8, HEJ, and HEJ+PYTHIA were made for this analysis and the results for the average total number of jets are shown in figure 8. We emphasise that the additional jets are no longer required to be in the rapidity interval between the two hardest jets, and there is no longer a significant disparity between their transverse momenta and that of the two hardest jets. This results in a greater number of jets passing the selection cuts, and consequently the potential for variation between different predictions is slightly higher.
Also shown in figure 8 as a shaded red band around the central predictions for HEJ+PYTHIA are variations of the merging scale k ⊥M (with a central scale of 15 GeV) between 7.5 GeVand 30 GeV. k ⊥M should be set to a value below the minimum jet transverse momentum used in the analysis, which in this case is 35 GeV. We see that even for these very exclusive multiplicity-dependent observables, allowing the merging scale to get very close to the analysis scale leads to only modest variations, and do not exceed the JHEP09(2018)074 size of the HEJ renormalisation and factorisation scale uncertainties. As this plot is most sensitive to differences between HEJ, PYTHIA8 and HEJ+PYTHIA, we expect the merging scale dependence in other plots to be comparable or smaller than that observed here.
In figure 8a the average number of jets is shown as a function of the rapidity interval between the most forward and backward jets ∆y f b ; we expect this to be particularly sensitive to the logarithms inŝ/|t| ∼ e ∆y summed by HEJ. The predictions of PYTHIA8 are significantly lower than those of HEJ and HEJ+PYTHIA, and moreover are outside the scale variation band for pure HEJ beyond ∆y f b > 4.5. This implies that in this regime, the merging of HEJ with PYTHIA8 increases the number of wide-angle jets relative to PYTHIA8 alone, as we should expect. Such differences should be even more pronounced with a larger centre-of-mass energy than the choice of √ s = 7 TeVwhich was used for these comparisons.
It is interesting to note that not only is the spread of predictions significantly larger in figure 8a than in figure 7a, but also that the prediction of HEJ+PYTHIA is lower than that of HEJ alone. There are several possible explanations for this. Firstly the addition of a parton shower extends the shower profile beyond the jet radius, such that potentially fewer of the jets from the partonic calculation pass the relevant criteria. Secondly, at ∆y = 10 two partonic jets of 45 GeV transverse momentum would take up all the energy available at a 7 TeV collider, and all predictions for the average number of jets would therefore have to return to 2 at this point. Since the parton shower uses some of the available collider energy in (for example) the description of the underlying event, the turnover of the average number of jets will have to happen earlier than in the pure partonic prediction. Alternatively it could be that too many non-FKL configurations of lower multiplicity are being inserted by the merging algorithm, an issue that could be resolved by the extension of this method to merge non-FKL events.
Finally, in figure 8b the average number of jets is shown as a function of the scalar sum of transverse momentum H T ; we expect this observable to be sensitive to the double logarithms in transverse momentum summed by the parton shower. PYTHIA8 now adds further hard radiation to that of HEJ, which is both as expected and is consistent with the previous results.
The choice of more inclusive observables and simpler selection cuts leads to a clearer separation of the effects of the logarithms included in the parton shower and those of the all-order summation of high energy logarithms in BFKL or HEJ. A simple experimental investigation with a similar set of cuts and distributions would be extremely interesting in exposing the shortcomings of either predictions, and the benefits of the combined formalism presented in this paper. Such an experimental analysis would further aid the development of predictions valid for the separation of the VBF and GF contribution to Higgs-boson production in association with dijets.

Summary and outlook
We have introduced a new CKKW-L-inspired merging algorithm for combining the allorder summation of high energy logarithms in HEJ with a parton shower. For the first time HEJ events have been fully evolved down to particle level using the modern parton shower,

JHEP09(2018)074
hadronisation and modelling of MPI in PYTHIA8. The merging algorithm systematically combines the dominant perturbative corrections due to hierarchies of transverse momenta (i.e. of soft and collinear origin) from the parton shower with those due to large invariant masses between jets of similar transverse momenta, as implemented in HEJ or BFKL.
The performance of the merging algorithm was assessed by considering observables which measure the additional radiation in the rapidity interval between two tagging jets. Many of the observables measured so far have (intentionally or not) a hierarchy of transverse scales induced, and so require a systematic resummation of the logarithms from the parton shower in order to arrive at a satisfactory description. For such observables we find that the description of HEJ+PYTHIA is consistent with standalone PYTHIA and data. The improvement upon HEJ in such distributions is notable. In addition, an investigation of related observables but with more inclusive cuts demonstrated that the jet multiplicity for large rapidity intervals is increased relative to PYTHIA through merging. A measurement of such clean observables can serve as test of high energy evolution. These results demonstrate that we have combined effects originating from both parton shower and from HEJ, providing a proof of concept for this method.
Notwithstanding what has been so far achieved, what has been presented constitutes a first attempt at merging HEJ with a parton shower. We envisage numerous refinements that can be made. There is a need to implement a prescription for incorporating full fixed-order matching into the merging procedure and the inclusions of sub-leading partonic channels (non-FKL) to the HEJ resummation [38]. In particular this will have an impact upon which states may be inserted by the parton shower. A systematic inclusion of such events in the prescription would not require dramatic changes to the algorithm. Firstly, the definition of what constitutes a HEJ state would need to be extended; secondly, the appropriate tree-level matrix elements should be used when calculating the veto probability of trial emissions. Nevertheless, the impact upon the observables discussed in this paper should be relatively modest; this was assessed by studying the relative size of the contributions of fixed-order non-FKL events in pure HEJ.
As discussed in section 4.3, a limitation to the method is that only the hardest emission of the shower received subtractions in their associated splitting kernel. This limitation could be addressed by re-inserting HEJ emissions in those events that were modified by PYTHIA above the merging scale, at the appropriate evolution scales reconstructed during the parton shower history. However, such a procedure has several ambiguities, such as where in the (modified) colour flow the emission should be inserted, and precisely how the recoils should be performed. Preliminary studies indicate the effect of reinserting HEJ emissions has a small effect, even upon the most sensitive observables shown in figure 8. However, we postpone a systematic study of these effects to a future publication.
Finally, also noted in section 4.3, a more advanced treatment for the weighting of colour flows in HEJ events that is informed by the parton shower may be necessary. The impact of this last effect is not obvious, and its resolution will require further study.
In this paper we considered the effects of our merging algorithm in pure dijet analyses. Partially this was due to the availability of data; in addition it is worthwhile to first consider the effects of a new method in a simpler environment where there is no expectation of new JHEP09(2018)074 physics. Nevertheless it is also important to apply this method to processes other than those which are purely QCD. Since one of the primary motivations was to assess the impact of jet vetoes that are relevant for Higgs plus dijets studies, this process is the next natural arena for study. We emphasise that this should not require any significant modifications to the method; the task is primarily an exercise in software development, rather than a theoretical challenge.
Finally, although we chose to implement this method for PYTHIA, in principle it should be possible to implement for other parton showers. It would be interesting to compare the effect of merging HEJ with different choices of parton shower. It would also be informative to perform the jet analysis with a much harder jet threshold, such that the sensitivity to the tunings of the non-perturbative effects are reduced. This would result in a much cleaner comparison of the perturbative merits.
Although we have been able to draw many positive conclusions by comparing with experimental data, the cuts that were chosen are not conducive to examining the effect of high energy logarithms. Both these points entail that it is difficult to discriminate between theoretical predictions that model different physics and should be expected to differ. We hope that as more data is collected, future analyses will examine a similar set of observables but with more inclusive cuts, as discussed.
This work has reinforced the notion that the interplay between different types of logarithms is not necessarily straightforward, and that there are circumstances under which the combination of two all-order summations is necessary. We hope that the merging algorithm we have developed may be used in future as a tool to inform analyses what selection of cuts and observables are sensitive to parton shower effects, high energy effects, or both.