Double parton scattering in four-jet production in proton-proton collisions at the LHC

We study the contribution from double parton scattering (DPS) to four-jet production at the LHC, both at the leading order accuracy and after incorporating the effects of QCD radiation. Apart from DPS, we also include and discuss the contribution from single parton scattering (SPS). We find that the QCD radiation impacts theoretical predictions significantly, with DPS contributions more affected than the SPS. We also examine a number of observables in regard to their effectiveness for discrimination between DPS and SPS events and propose sets of kinematical cuts to improve the prospects of measuring DPS in four-jet production.


Introduction
A composite nature of hadrons leads to a complicated structure of the underlying event in hadronic collisions. In particular it gives rise to a possibility of several interactions per one collision, a phenomenon referred to as multiple parton interactions (MPI). Rapidly increasing fluxes of partons in hadrons at small momentum fractions x make their occurence more frequent at higher collision energies or, alternatively, at lower invariant masses of the measured system. A particular subset of MPI with two hard interactions per single hadron-hadron collision is called double parton scattering (DPS). It is a simplest possible MPI system. Provided the final state carries enough transverse momenta, it also is a relatively clean system experimentally. DPS processes have been first observed by the AFS [1] and the UA2 [2] collaborations. After these pioneering works a series of measurements were performed at the Tevatron collider [3,4,5,6,7,8,9,10] and at the LHC [11,12,13,14,15,16,17,18,19,20,21,22,23].
Among different DPS production channels the four-jet DPS production takes a special place due to high abundance of the multi-jet events in hadronhadron collisions and correspondingly, a large DPS cross section. The four-jet DPS production was measured in multiple experiments at various colliders: ISR [1], SPS [2], Tevatron [3] and the LHC [16,18,23]. The theoretical efforts to describe (four-) jet production from multiple scatterings have a long history [24,25,26,27,28,29]. In more recent years studies focused on, among other, including modelling of radiation effects [30] and LHC's potential to gain new information on the properties of two-parton distribution functions in the transverse plane [31]. The importance of four-jet production in the context of DPS studies, especially in the back-to-back regime, was further highligted in [32,33,34]. In this series of papers, concurrently to other efforts at that time, a theoretical DPS framework that accounts for perturbative splittings of a single parton into two was devel-oped. It was later followed by a proposal to implement the new approach in the PYTHIA [35,36] event generator through modifying its MPI model with the help of a dedicated tune [37]. Furthermore, the DPS observables and kinematics of four-jet production [38,39] as well as calculations in the high-energy factorization approach [40,41] were explored. In parallel, there has been an enormous progress in theoretical understanding and decription of DPS on a more fundamental level [33,34,42,43,44,45,46,47,48,49], see also recent review in [50]. Concurrently to the aforementioned advances in DPS description, development of Monte Carlo DPS and MPI models progressed significantly [51,52,53,54,55,56,57,58,59].
The DPS four-jet production unavoidably occurs simultaneously with the four-jet production in single parton scattering (SPS). Theoretical predictions for this process are known at the next-to-leading order (NLO) accuracy [60,61]. The di-jet production, which is a major building block of double scattering resulting in four jets, has been studied up to the next-to-next-toleading order (NNLO) in perturbation theory [62]. Although a significant progress has been made recently towards DPS calculations at NLO [49,63] no DPS process has been described at this accuracy yet.
In this work, we revisit the theoretical predictions for four-jet DPS production in pp collisions at the LHC. We first review the LO calculations at the partonic level, both for DPS and SPS, and discuss selected differential quantities as well as their uncertainties. We also investigate the impact of longitudinal parton correlations present in double parton distribution functions, as defined below. The LO analysis is then extended by studying the effect of initial and final state radiation on both the SPS and DPS total cross sections and distributions. The radiation is simulated with the parton shower algorithm as implemented in the PYTHIA event generator.
In this sense, the work presented here goes beyond parton level LO calculations in the collinear factorization [38,39], or adding to it one (undetected) real emission, as done in [30]. Our studies can be also seen as parallel to the effort of including higher-order effects in the DPS cross sections using the framework of high-energy factorization of [40,41]. The developments presented here are also different from the option of second hard scattering built in PYTHIA. While PYTHIA constructs two-parton parton distribution functions from single parton distributions in a very specific way which cannot be changed by a user, and requires additional adjustments to obtain correct normalization of DPS cross sections, our approach allows to specify how double parton distributions are modelled. Besides, in earlier versions of PYTHIA (before 8.240) the DPS events are dependent on the ordering of the hard interactions constituting the DPS event. Furthermore, our approach offers another way to study effects of showering on DPS predictions in nuclear collisions [64], different from the Angantyr model of PYTHIA [65] 1 .
The paper is structured as follows: in Chapter 2 we briefly review the DPS framework used in our studies. This is followed in Chapter 3 by the description of the implementation and details of various technical aspects of the simulation. This chapter also contains our numerical results and their discussion. We conclude and present an outlook for future work in Chapter 5.

Phenomenology of double parton scattering
We begin with a description of the framework in which four-jet DPS production is studied. Its origins go back to the work of Paver and Treleani [25] and Mekhfi [67]. The results of [25], [67] were later generalized and extended to the case of n-hard interactions [42,43].
The total DPS cross section is given by the following expression where A and B denote final states in i k → A and j l → B processes and objects Γ ij/h are called generalized two parton distribution functions which can be in a first approximation thought of as a probability to find two partons i and j with longitudinal momentum fractions x 1 and x 2 separated by distance |b| in a transverse plane of a hadron h.
In the following, we assume factorization of Γ ij/h into a product of longitudinal and transverse dependent pieces Using Eq. 2 one can write a total DPS cross section for pp collisions as where we have dropped the subscript h in order to simplify notation. The prefactor 1/ (1 + δ AB ) in Eq. 3 was introduced to reflect the fact that in case of production of two indistinguishable finals states A and B one has to divide a total DPS cross section by 2. In the following we refer to the distributions D ij in Eq. 3 as to double parton distribution functions (dPDFs). The quantity σ eff in Eq. 3 is given by and can be interpreted as an effective interaction area. It should be noted, however, that the factorization of longitudinal and transverse pieces is an assumption driven by a practical need of modelling of two-parton parton distribution functions for DPS production in pA collisions. phenomenology purposes. The generalized distributions evolve differently in the position and in momentum space [42,43], and consequently at b = 0 and b = 0 [44]. This is inconsistent with the factorization assumption, which should hold for all values of b. As discussed in e.g. [49], the cornerstone theoretical expression for the DPS cross section is given by Eq. (1), with Γ ij/h evolving according to a homogenous double DGLAP equation [42,43]. Assuming no partonic correlations in x-space one can substitute in Eq. 3 which gives us the "pocket formula of DPS" Such factorization violates momentum and number (flavour) dPDF sum rules proposed by Gaunt and Stirling [68]. One can avoid unphysical contributions by multiplying a factorized product of PDFs by an appropriate cutoff function, for example where θ(1 − x 1 − x 2 ) excludes unphysical region where x 1 + x 2 > 1. However, one has to keep in mind that Eq. 7 still violates momentum and number sum rules for dPDFs. Notwithstanding all the above concerns, Eqs. 3 and 7 are still commonly used for a phenomenological modelling of DPS. In this paper we refer to the dPDFs approximated according to Eq. 7 as to "naive" dPDFs. We also should notice that "naive" dPDFs, as in Eq. 7, do not allow to reduce Eq. 3 to the "pocket formula" Eq. 6. Instead, by substituting Eq. 7 into Eq. 3, we get In our analysis we use two different models of dPDFs. Namely, we use "naive" dPDFs, as in Eq. 7, constructed out of MSTW2008 [69] or CT14 [70] leading order (LO) PDFs. In order to estimate the impact of partonic correlations in xspace we compare between results obtained according to Eq. 3 and supplemented with either "naive" dPDFs or with GS09 dPDFs [68]. The latter are LO dPDFs, which evolve according to the double DGLAP equation and obey the momentum and number sum rules. Their initial parametrization is predominantly based on MSTW2008LO as input single PDFs. At the root of the difference in the evolution of generalized distributions in the position and in momentum space lies perturbative splitting of one parton into two partons [42,43], a mechanism which contribution to double parton distributions needs to be accounted for. The corresponding 1 → 2 term in the Γ ij/h has a 1/b 2 behaviour at small b, and renders the cross section in Eq. 1 UV-divergent. The UV divergence is an artefact of using the DPS description outside of its region of validity, where the SPS picture is more suitable. A consistent scheme which treats the problem of UV divergencies, as well as a closely related problem of double counting between DPS and SPS, was proposed in [49]. Earlier approaches and discussion of the problems can be found in [33,34,45,46,47,48]. According to the scheme of [49], both DPS cross section and SPS cross section (at the same order in perturbation theory as DPS) contribute to the total cross section, while the double counting and the UV divergencies are removed by suitably designed subtraction terms. However, in the case of four-jet production it is not possible to apply this prescription since the corresponding SPS calculation at NNLO is technically out of reach. Instead, we rely on the observation made in [49] that the SPS contribution and the associated subtraction terms lose their relevance if the considered system probes sufficiently low values of x in two-parton distribution functions (see also [34]). The four-jet production at relatively low p T can be then seen as a promising setup to study DPS, in addition to e.g. same-sign W-pair production [49,71,72,73,74].
Observables used in DPS phenomenology can be often characterized as "transverse" and "longitudinal", depending on the direction of momenta of the final state particles which are predominantly probed. Both types of observables exploit the fact that the individual partonic collision in a DPS event must obey momentum conservation, i.e. make use of correlations between final state particles either in DPS or SPS topologies. In addition, the rapidity-based observables not only exploit differences between SPS and DPS topologies [75,76], but also are directly sensitive to the correlations among the incoming momenta fractions for the parton coming from the same hadron. It has been observed in e.g. [75] that the DPS observables of the transverse type are particularly sensitive to real radiation effects. In particular, in many cases the DPS-sensitive variables based upon p T imbalance and angular correlations between produced jets are trivial at the partonic level. Consequently, adding the real radiation results is needed in order to obtain a realistic description of these observables.

Predictions for four-jet production in pp collisions
The numerical calculations of DPS production at the parton-level are performed using an in-house built Monte Carlo programme, which calculates LO matrix elements for 2 → 2 ⊗ 2 → 2 scattering. In order to add initial and final state radiation (ISR and FSR) effects to our parton-level simulations we use the PYTHIA event generator [36,35] with modifications necessary to read and "shower" output of our DPS code [77]. In the first step we output the DPS events using modified ("double") Les Houches Event (LHE) file records [78], see Appendix A for more technical details. The double LHE files are then supplied to PYTHIA for showering using SecondHard:generate = on and PartonLevel:MPI = off settings. While generating partonic events, we use the original Lund Monte Carlo algorithm proposed by Bengtsson [79] which allows to generate colour charges of the initial and final state partons within the leading colour approximation 2 . The generation of the colour charges is required to take into account colour coherence effects [80], [81]. For each individual 2 → 2 scattering, all 8 subprocesses involving combinations of (anti-)quarks and/or gluons are taken into account. The scales of the two partonic collisions are treated as independent and chosen equal to the value of the transverse momentum partons in the two di-jets. In each case, the central value of the factorization scale is set to be equal to renormalization scale, where p i is the momentum of one of the partons in the collision i and i ∈ {1, 2} with 1, 2 indicating the two scatterings. The SPS predictions are obtained using the MadGraph package [82]. In order to account for the NLO effects in the normalization of the SPS total cross section,we apply an effective K-factor [60,61] to the LO predictions, similarly to approach of Ref. [39,41]. The scale of the SPS process is chosen as [60,61]. In the analysis at the partonic level we produce four partons and apply to them a jet separation criterium of R ij = (y i − y j ) 2 + (φ i − φ j ) 2 > 0.4, unless otherwise stated. In the parton shower analysis final state particles are clustered by means of the Fast-Jet [83] package with the anti-k t jet clustering algorithm [84] and jet radius taken to be R = 0.4. Whenever comparing the SPS predictions to the showered DPS results, the SPS results are also showered using the same version of PYTHIA and the same method of clustering of partons into jets is applied. Unless otherwise stated, we consider a rapidity coverage of −4.7 < y j < 4.7 for all jets j. Since the generation of DPS events with jets carrying high transverse momentum is suppressed due to low fluxes of incoming parton pairs and, correspondingly, high-p T jets originate predominantly from SPS, we apply an upper cut on p T,j for all jets.
A number of checks at various stages of the calculation have been performed. The MadGraph set-up for the computation of SPS cross sections have been crosschecked against the outcome of the ALPGEN code [85] for the process pp → jjgg. The cross section for results for the individual 2 → 2 cross sections, i.e. the building blocks of the DPS cross sections, have been also checked against MadGraph. The procedure to assign colour to initial and final state partons in LHE files was checked by comparing PYTHIA-showered results against results based on MadGraph LHE events, also showered with PYTHIA. Our results for the ratios of the DPS to SPS+DPS total cross sections agree with results of LO partonic results of Ref. [39] within a few percent, depending on the choice of the value of R ij . We have also checked that with our set-up for SPS calculations we can reasonably well reproduce the LO and LO matrix element matched to parton shower (ME+PS) four-jet cross sections from [60].

Parton-level analysis
We begin our studies by investigating the DPS and SPS four-jet production at the partonic level. In order to judge the reliability of the predictions, in Fig. 1 we check how well they fare against CMS measurement of four jet production at √ S = 7 TeV [16], which uses relatively low cuts on jets transverse momentum: p T > 50 GeV for the two most leading jets and p T > 20 GeV for the third and fourth jet. In accordance with experimental analysis, we use here R jj = 0.5.
The uncertainty bands for DPS and SPS events are estimated independently and combined together afterwards. For the SPS events we start by generating central value predictions at LO using MSTW2008 LO PDFs. Since NLO corrections typically lead to the change in the normalization of the cross section predictions, similarly to the approach of [39], we multiply the LO predictions by a K-factor of 0.5. Such value of the K-factor follows from the fixed-order four-jet NLO in [61] and has been obtained for the MSTW PDFs. After producing central value predictions we obtain the conservative estimate of uncertainties due to the scale variation by performing simultaneous changes of factorization and renormalization scales up and   The theoretical predictions are obtained as a sum of the SPS result (including the K-factor of 0.5) and the DPS result, with their corresponding error bands. While we are not in position to calculate the DPS prediction at NLO, nor have additional information on the NLO effects 3 , we take the spread provided by the scale uncertainty and variation in the parameter of σ eff = 15±5 mb for two different PDF sets, MSTW2008LO and CT14 LO, used in the "naive" dPDFs construction, Eq. 7. These effects are then combined using the envelope method [86,87]. Finally, the central values of the DPS and SPS predictions as well as the uncertainty estimates are combined. We observe an overall agreement within uncertainty with the data for both each jet's p T and y, see Fig. 1, with tendency for theoretical results to overestimate the data.
In Table 1, we compare central values of the total cross cross sections for the LHC collision energy of 7 and 13 TeV and two sets of cuts on jets' p T : a stricter 35 GeV < p T,j < 100 GeV and a looser one 20 GeV < p T,j < 100 GeV. The upper cut on jet's p T is inspired by the well-known fact that DPS is enhanced at lower partonic energies, and by the earlier literature [39]. In the DPS computations we use "naive" dPDFs constructed out of MSTW2008 LO PDFs. Correspondingly, the SPS results are also obtained with MSTW2008 PDFs. Apart from the least DPS-favourable case of √ S = 7 TeV and 35 GeV < p T,j < 100 GeV, in all other cases the upper cut on maximal p T,j leads to higher total DPS cross sections than SPS cross sections. As expected, the DPS cross sections increase greatly when the minimal p T,j cut is lowered, and the ratio of DPS to SPS cross sections improves. Due to growing parton fluxes at low fractions of momenta x, an increase in the ratio is also observed at higher LHC energies.    There we show predictions for the leading jet p T and the rapidity difference ∆Y ≡ max|y j − y k | distributions [75] for the same values of LHC energy and the same sets of cuts on jets' p T,j as in Table 1. Note that, unlike distributions in Fig. 1, we do not combine DPS and SPS predictions together. Moreover, we plot DPS distributions generated with MSTW2008 and CT14 PDFs separately, such that the effect due to the choice of the PDF set remains clearly visible. Shaded areas in all plots show the size of scale variation error, obtained by varying the central renormalization and factorization scales by a factor of 2 and 1/2 simultaneously. Hatched areas correspond to the additional variation of the σ eff parameter in the 15 ± 5 mb range, corresponding to the measurement in four-jet final state by the ATLAS collaboration [18]. Histograms in Figs. 2-5 confirm naive expectations that lowering the cut on minimal p T of jets leads to an increase in the DPS scattering vs. SPS scattering at low p T of the leading jet, as well as at high ∆Y 4 . The shapes of the distributions, shown in lower panels of Figs. 2-5, also become more steep as the minimal p T -cut decreases. Similar qualitative changes are observed as the energy of the hadronic collision increases.

Cuts and collision energy
In the remaining part of this paper we will present the DPS predictions obtained only with dPDFs built on the basis of the MSTW2008LO PDFs. Firstly, the only publicly existing dPDFs package, GS09 [68], is built on the basis of the MSTW2008 LO parametrization. Secondly, Figs. 2-5 clearly demonstrate that although the DPS predictions come with a large theoretical uncertainty, qualitative behaviour of the results is similar, independently of the LO PDFs used. Furthermore, the SPS predictions we refer to [60,61] are obtained using the MSTW2008 PDFs. Finally, most of our analysis presented in the following is concerned with shapes of the distributions, where the PDF effects mostly cancel out.          One of the observables which is considered to have a good discriminating power between DPS and SPS is the azimuthal difference between the two most remote rapidity jets, ∆φ jj [39]. We show this distribution for √ S = 13 TeV and 35 GeV < p T,j < 100 GeV. Indeed, as Fig. 6 demonstrates, while away from ∆φ jj = π the DPS events are expected to be almost equally distributed, the SPS events prefer back-to-back configurations for the two most distant jet in rapidity, leading to a strong depletion of the cross section at small ∆φ jj . The small dents in ∆φ jj DPS distribution are due to jet clustering, since the jet separation obviously cuts out a certain number of events where the partons have small angular separation.  Our results given in Figs. 2-5 agree with the results of [39] at qualitative level. We shall note, however, that the qualitative comparison between ∆φ jj distributions in [39] and our results in Fig. 6 shows some differences at small and large values of ∆φ jj . More precisely, the ∆φ jj distribution given in [39] does not have a peak at ∆φ jj = π as well as curvatures at small and large values of ∆φ jj caused by the In a next step we investigate the impact of longitudinal correlations within the pairs of partons taking part in the double scattering, as implemented in the GS09 package, on the predictions. To this end, in Figs. 7 and 8 we compare the leading jet p T distributions at √ S = 13 TeV obtained using the GS09 dPDFs and using "naive" dPDFs for the same two sets of cuts on p T,j as above. For the whole range of the p T considered, as well as for the practically accessible values of ∆Y, the  Leading jet p T distribution in four-jets DPS production at √ S = 13 TeV, calculated using "naive" dPDFs based on MSTW2008LO and GS09 dPDFS. Upper two panels show results obtained with the 35 GeV < p T,j < 100 GeV cut, lower two panels are for the 20 GeV < p T,j < 100 GeV cut. difference between the two distributions is not more than 10%. At very high values of ∆Y, where large x values are probed, this difference grows bigger. However, given that the effect is overall small, we will not consider it in further studies and from now on only use the "naive" dPDF modelling.

Impact of QCD radiation
We now move on to discussion of the impact of the initial and final state radiation on the DPS and SPS four-jet predictions. We use LHE files with events generated at are obtained using the "naive" dPDFs based on the MSTW2008LO set. Since the effect of adding radiation is expected to be universal, i.e. independent on the form of the dPDFs, this is sufficient for the discussion to follow. The LHE files with the DPS events are generated by our in-house C++ code, while the SPS events are generated with the MadGraph event generator [82]. In the following, the SPS events are not longer multiplied by a K-factor.
In order to add radiation effects to our simulations we supply SPS and DPS LHE files to the PYTHIA event generator. To shower the SPS events we follow a standard procedure where we first create LHE files using the MadGraph package and supply them to PYTHIA to add ISR and FSR to our simulations. The DPS case, however, is more involved: since the standard LHE event record contains information only about one hard interaction, one has to adapt the LHE files to contain also information about the second hard interaction in a DPS process, see Appendix A. The resulted "double" LHE files are passed to the PYTHIA code which was modified for this purpose [77]. In this way, the information on two individual hard scatterings propagates through the shower.
The DPS events are then showered in the interleaved way in accordance with the MPI model of PYTHIA8 [94]. The approach of [94] has several advantages comparing to the simplified case where two hard processes constituting a DPS event are showered independently from each other. The DPS showering performed in accordance with [94] implies that ISR and FSR cascades satisfy constraints imposed by energy and momentum conservation at each evolution step. In particular, in this way it is ensured that the kinematic constraint on Bjorken x-es given by Eq. 7 is preserved by the ISR cascades at each step of the ISR evolution. Another advantage of the approach of [94] is the fact that the ISR cascade which first reaches the proton reduces amount of energy and momentum available for the second ISR cascade. Finally, we shall note that the interleaved ISR evolution accounts for the changes in the parton content of the beam remnants which introduces non-trivial correlations in flavour of interacting partons [53]. Therefore, our simulations provide a better treatment of the beam remnants comparing to the case where one "showers" two hard processes completely independently. 6 Next, we discuss our results for the 4-jet production. Apart from showing the absolute value of the differential cross sections, we also study their shapes. In all plots in this section, we show an estimate of the statistical error on the (nominal or normalized) cross section in each bin.
It is well-known that if the same p T cut is applied on the two observed jets the higher-order calculations for di-jet production become unstable due to restricted phase space available for soft gluon radiation [95]. Although the four-jet production is not affected by this problem, the contributions to the DPS production which originate from double di-jet production might be. For that reason, apart from applying symmetric cuts of p T,j > 35 GeV on all four jets, we additionally introduce asymmetric cuts by requesting p T,1 > 55 GeV for the leading jet and p T,j > 35 GeV for the remaining jets. 7 In Figs. 9 and 10 we present distributions in the leading jet p T and ∆Y studied in the previous section, now including the radiation effects. We find that adding radiation decreases the total cross section for both the SPS production and DPS production, see also the first row of Table 2. The exact value of the reduction factor depends on the set-up of the calculations, in particular the chosen set of kinematical cuts, but in general the DPS predictions are more impacted by the radiation. As can be seen from the leading jet p T distribution discussed in the previous section, most of the partonic DPS events happen to have jets with very low p T , just passing the cut. This makes the DPS distribution much more vulnerable to the effects of radiation, meaning the adjustment of jet's p T due to radiation can cause a substantial proportion of the events not to pass the cuts. The same is not true for SPS, where the peak of the leading jet's p T distribution is at much higher values of p T . Having studied origin of jets passing the selection cuts, it 6 We shall note that the PYTHIA approach to the ISR evolution of the DPS (MPI) events does not take into account cases when two ISR cascades merge into a single ISR cascade at certain step of backward evolution (so called "1v2" events). Whereas the first results on implementation of such effects in to the parton-shower framework for the same-sign W W production were recently reported [74], the implementation of the aforementioned effects for the four-jet DPS production is highly non-trivial [49] and is beyond the scope of this work. 7 However, one should note that in order to make sure no instability of that sort affects the DPS calculations entirely, fully asymmetric cuts, i.e. with different minimal pT for each jet, would be required. Given that the cuts need to be sufficiently separated, a measurement of DPS based on events selected in that way might be very difficult due to low statistics.  also appears the SPS events have a higher partonic center of mass energy √ŝ , and correspondingly higher values of Bjorken's x carried by the partons participating in a collision, in comparison with individual DPS collisions. As a consequence, radiation in the SPS events more often gives rise to an extra jet passing the cut. In this way some SPS events which at partonic level would have not passed the cut are accepted now. The same effect is much suppressed for the DPS due to lower x values of incoming partons. The normalized distributions, on the contrary, seem to preserve their main features first observed at the partonic level, cf. Fig. 9. For example, the DPS leading jet's p T distribution, though it gets flattened out, remains more peaked at small p T than the SPS one, or the ∆Y DPS distribution stays flatter at very high ∆Y. It is also worthy noticing, that showering alters the leading jet p T distributions more than the ∆Y distribution, in accordance with the observation of higher impact of radiation on the observables of the "transverse" type. Introducing the asymmetric set of cuts does not lead to significant changes in the behaviour of the distributions.
The whole analysis shows clearly that accounting for radiation can dramatically change the size of the four-jet cross sections calculated under assumptions of certain sets of cuts, and conclusions derived from analysis of the partonic level do not hold after a more realistic simulation of the production process is employed. It is interesting to note that a similar dampening effect of radiation on DPS predictions have been observed in the k T factorization framework [40,41]. However, a comparison between Fig. 4 and Figs. 9, 10 indicates that the radiation effects can also influence the differences between shapes of DPS and SPS distributions. Hence, even if only the information on normalized distributions is used for distinguishing between SPS and DPS, the impact of radiation should be taken into account.
Next we study the impact of radiation on the differential cross section in ∆φ jj , discussed at the partonic level in the previous section, see Fig. 11. While in that case DPS distribution was higher than the SPS for smaller angular differences ∆φ jj , adding radiation turns the SPS production mechanism into a dominant one everywhere. As expected, the typical peak in the DPS distributions calculated at LO due to back-to-back configurations, observed here at maximal value of ∆φ jj , vanishes once adtional radiation is introduced. Furthermore, the radiation also changes the shape of the DPS distribution from flat to rising at higher ∆φ jj , making it almost indistinguishable from the shape of the SPS distribution. This indicates that the ∆φ jj observable does not deliver an efficient discrimination between DPS and SPS mechanisms, in this setup of cuts. However, as we show later, the discrimination power of various observables can be improved on if additional cuts are introduced.
In the following we are going to use different combinations of commonly used DPS-sensitive observables constructed to increase the fraction of DPS events. The general idea behind construction of such observables relies on different kinematics of jets produced via DPS and SPS mechanisms. Since the DPS events are essentially constructed out of two di-jet events, we expect them to give rise to two di-jet topologies well separated in rapidity and azimuthal angle φ and also well balanced in p T .
Using the aforementioned consideration as a guiding line we can study discriminating power of different DPS sensitive observables. In Figs. 12, 13, 14 and 15 we show the absolute and normalized DPS and SPS differential cross sections in other commonly studied observables in this context, such as the momentum imbalances as in four-jet DPS ATLAS [18] measurements, and the closely related observables defined in [37] and [30] respectively. We also study the azimuthal differences as in [18], and the invariant mass m ij of the two jets, i and j, with lowest ∆ p T ij , see Figs. 16, 17, and 18, respectively. We also shall note that in Eqs. 9 -11 and Eqs. 13 -14 we assume that | p T 1 | > | p T 2 | > | p T 3 | > | p T 4 | whereas in Eq. 12 we consider all non-equivalent combinations of jets i, j, k, l. For all observables, there is no value at which the DPS production would prevail over the SPS. We also find that some shapes of the DPS distributions at the parton level do not survive after showering. The case in point is the S ⊥ DPS distribution, showing distinct peaks at the minimal and maximal value of S ⊥ at the parton level [30] which get washed out by the radiation effects, cf. Fig. 15. On the positive side, for some observables the shape of the DPS and SPS distributions can differ substantially, as shown in the lower parts of Figs. 9-18. In particular, compared to SPS, we observe narrower distributions for DPS at small values of ∆ p T 12 and ∆ p T 34 . These regions correspond to back-to-back configurations of the two di-jets. The observed enhancements are then expected from theoretical considerations [32,33], if we identify the two leading jet pair and the subleading jet pair as coming from two separate hard collisions, which can happen rather often. We also note that in the back-to-back region the contributions to the cross section from the perturbative splitting mechanism are less relevant [34], meaning our predictions should not be much affected if these contributions are taken into account in the simulation. Similarly to the leading jet p T and ∆Y distributions, also for the other observables shown here we do not see a significant differences between results obtained with our symmetric and asymmetric cuts.
The information on shapes of the distributions can be harvested to improve discrimination power between DPS and SPS by imposing additional cuts. The choice of the cut parameters is driven by the shape of the distributions. In Table 2, we list values of the DPS and SPS total cross sections for various sets of cuts before and after radiation is included. Additionally, we provide percentage of DPS contributions to the total cross sections. We observe that, for our choice of the cut parameters, cuts on ∆ p T 12 and ∆ p T 34 increase the DPS signal in the most efficient way, even yielding some regions of phase-space where the DPS signal dominates, see Table 2 and Figs. 19, 20. A further significant improvement can be achieved by          (1/σ) dσ/dp T (1/GeV)  Figure 19: Same as in Fig. 9 but after imposing additional cuts ∆ p T 12 < 0.2 and ∆ p T 34 < 0.2.
combining these cuts with cuts on the invariant mass m ij of a di-jet pair with the smallest value of transverse momentum imbalance ∆ p T ij ,see Figs. 21, 22, 23. At this level, the difference in the number of events, i.e. lower statistics due to more stringent asymmetric set of cuts compared to symmetric set becomes visible. In this case, it might be advisable to decrease the cuts on the value of p T for all jets in the asymmetric set.  Table 2: Impact of parton shower (PS) effects on the DPS and SPS cross sections, √ S = 13 TeV, |y| < 4.7. The parton-level SPS cross sections are multiplied by K-factor equal to 0.5. All cross sections are given in nb. The DPS fraction is given at the 1% precision level.
Before finishing this section we would like to note that the different DPSsensitive observables considered in this paper were already used in earlier publi- (1/σ) dσ/dp T (1/GeV)  Figure 21: Same as in Fig. 9 but after imposing additional cuts ∆ p T 12 < 0.2, ∆ p T 34 < 0.2 and m ij < 100 GeV.   cations [18,30,37,39,72]. However, to the best of our knowledge, a combination of cuts on the aforementioned observables listed in Table 2, and the corresponding distributions in Figs. 19 -23, are new. We hope that the new combinations of cuts proposed in this work would allow to achieve better separation between DPS signal and SPS background in the four-jet production processes and, as a consequence, to improve the estimate of value of σ eff .

Summary and discussion
In this work, we have studied double parton scattering in four-jet events at the LHC, both at the partonic level and incorporating the effects of QCD radiation.
To this end, we have developed a parton-level Monte Carlo simulation outputting modified LHE files event records, which then can be showered by the PYTHIA event generator to which additional modifications are applied. Apart from studying the impact of various cuts and collider energies on the DPS and SPS contributions to the cross sections and estimating their uncertainties at the partonic level, we have also investigated the effect of longitudinal correlations as implemented in the GS09 package. After adding QCD radiation, we find that it can affect DPS and SPS predictions significantly, with DPS contributions being modified more. We have also examined a number of observables in regard to their discriminating power between DPS and SPS. Applying just a basic set of cuts, we observe that many of these observables substantially differ in shape in a specific range of values. This information can be then used to propose more elaborated sets of cuts which increase the efficiency of the discrimination. In particular, we find that a combination of cuts on jets' p T and transverse momentum imbalance ∆ p T ij as well as invariant mass m ij of a di-jet pair with the smallest value of transverse momentum imbalance ∆ p T ij provides a very promising way for selecting DPS contributions in four-jet events.
In this foray towards including higher order effects in DPS predictions, we have implemented a simple model of dPDFs which only partially accounts for correlations between partons in the same proton. In particular, this model neglects contributions from a perturbative splitting of one parton into two. We have, however, checked that at the parton level the predictions obtained with a publicly available dPDF package GS09, which under assumption of transverse-longitudinal factorization of double distributions accounts for longitudinal correlations resulting from 1 → 2 splittings, differ only minimally from predictions obtained with our naive dPDF model. As discussed above, a consistent theoretical framework according to [49] would require generalized parton distributions dependent on the impact parameter. Their modelling involves an "intrinsic" and a "splitting" part, which both evolve according to the homogenous double DGLAP equation. Provided such set of double parton distributions, our parton shower calculations can easily accommodate them. We note that a complementary approach, relying on altering the parton shower to include 1 → 2 splittings was reported recently in [74,96].
Naturally, in this exploratory study, we could only explore a few particular setups for the calculations. It is conceivable that e.g. other kinematical cuts could lead to a better discrimination between DPS and SPS. Additionally, lowering a minimal cut on jets' p T would select proportionally more DPS events. On the other hand however, it would inevitably lead to bigger contamination from background. Studies of various set-ups are beyond the scope of this work, but we hope that our work and simulation tool, which is available on request, can be used for optimizing the measurement of DPS in four-jet events at the LHC in the future. Figure 24: An extension of the Les-Houches version="1.0" standard to the DPS events. In order to ease the reading we keep only one digit after comma for the components of the four-momenta. version = "1" files has been adopted for the DPS events, that's not yet the case for the LHE version = "3" files. Therefore, the usage of the tag <LesHouchesEvents version="3.0"> will invoke calling routines from the file LHEF3.cc, leading then to a wrong assignment of the mother-daughter labels and MPI, ISR and FSR scales. [23] CMS Collaboration. Study of double-parton scattering in the inclusive production of four jets with low transverse momentum in proton-proton collisions at √ s = 13 TeV. 2021.