Generation of central exclusive final states

We present a scheme for the generation of central exclusive final states in the Pythia 8 program. The implementation allows for the investigation of higher order corrections to such exclusive processes as approximated by the initial-state parton shower in Pythia 8. To achieve this, the spin and colour decomposition of the initial-state shower has been worked out, in order to determine the probability that a partonic state generated from an inclusive sub-process followed by a series of initial-state parton splittings can be considered as an approximation of an exclusive colour- and spin-singlet process. We use our implementation to investigate effects of parton showers on some examples of central exclusive processes, and find sizeable effects on di-jet production, while the effects on e.g. central exclusive Higgs production are minor.


Introduction
Compared to the fairly clean environment of e + e − annihilation, proton collision events are in general very messy, especially at the LHC at high luminosity. Even at lower luminosity where pile-up events are absent, the existence of multiple soft interactions and initial-state parton showers means that any hard sub-process of interest will be obscured by soft and semi-hard hadrons x 1 x 1 ′ x 2 ′ Fig. 1 The basic diagram for a general central exclusive process pp → p + X + p.
smearing the measurements. However, for some rare events, a central colour-singlet hard sub-process may appear in complete isolation, with rapidity gaps on both sides stretching all the way out to the (quasi-) elastically scattered protons, giving a nice and clean environment to study its properties.
Such Central Exclusive Processes (CEPs) have been extensively studied in the so-called Durham formalism, first described by Khoze, Martin and Ryskin in [1] and reviewed in detail in [2]. The simplest such process is Higgs production, where two gluons in a colour-and spin-singlet state fuse together via a top-quark loop into a Higgs particle as outlined in figure 1. With an additional virtual exchange of a (semi-hard) gluon, the net colour exchange between the colliding protons can be zero and we may end up with a very simple and clean final state consisting only of two (quasi-) elastically scattered protons along the beam pipe and the Higgs decay products in the central rapidity region. 2 The formalism can be generalised to any coloursinglet hard sub-process, and the main ingredients to construct the amplitude is the matrix element for this sub-process and the so-called off-diagonal unintegrated parton densities. The latter can be interpreted as the amplitude related to the probability of finding gluons in a proton with equal but opposite transverse momentum, q ⊥ , and carrying energy fractions x and x ′ each, one of which is being probed by a hard scale µ 2 . These densities also include a Sudakov form factor describing the probability that there is no additional initial-state radiation from the incoming gluon between the scales q ⊥ and µ, which could destroy the rapidity gaps. Additional emissions below q ⊥ are then suppressed since they cannot resolve the individual colours of the two gluons. A third ingredient is the so-called soft survival probability which gives the probability that there are no additional soft or semi-hard interactions between the colliding protons which could destroy the rapidity gap.
Implementing the Durham formalism for CEP in an event generator is fairly straightforward since the final states are quite simple and clean. The cross section for any sub-process can be decomposed in a central exclusive luminosity function which is folded with a coloursinglet matrix element in a specific spin state. Several implementations have been made [3,4,5], and in this paper will present yet another.
Our implementation is provided as an add-on 1 to PYTHIA8 [6] and is inspired by the observation that the Sudakov form factors in the off-diagonal unintegrated parton densities used within the Durham formalism can be interpreted in terms of no-emission probabilities in the parton shower language of PYTHIA8.
In this way, we can reformulate the cross section for producing a CEP event in terms of a probability that a standard inclusive sub-process generated by PYTHIA8 at some scale during the parton shower evolution is converted to a colour-singlet, and thereafter be considered a CEP event disallowing further initial-state shower splitting. The main advantage of this approach is that we actually are allowed to include initial-state shower splittings, and thus can approximately model higher order corrections to the original sub-process. In addition, we have the option of using the multiple interactions machinery of PYTHIA8 to directly model soft survival probability as suggested in [7].
Consider, e.g., the central exclusive production of di-jets. We would start by generating the basic 2 → 2 hard partonic scattering from the inclusive matrix element. We would then generate an initial-state parton emission from each of the incoming partons. This im- 1 The code uses the PYTHIA8 UserHooks machinery and is available on request from the authors.
plicitly includes the probability that no emission has been made at a higher scale than the two generated splittings. If the colour and spin state of the original 2 → 2 is consistent with a CEP, we basically take the ratio of the corresponding exclusive cross section and the one calculated with the inclusive matrix element, using the generated scales as factorisation scale. This gives us the probability to discard the generated splittings and continue the event as a CEP, or to keep the hardest emission and continue as a normal inclusive event.
If we continue the event as inclusive, we keep the hardest splitting and again generate one initial-state splitting from each side. Now that we have a threeparton final state which must be checked if it can be a candidate for CEP, but otherwise the procedure is repeated. It should be noted that PYTHIA8 in general does not assign spin states to particles, so the procedure here also involves a spin decomposition of the parton splitting probabilities and the matrix elements to correctly get the probability for this to be a CEP. This will become cumbersome when we go up in parton multiplicity, but is still fairly straightforward.
The fact that we can stop the parton shower at any stage and check if we can convert the generated state exclusive, does not only mean that we can approximate higher order contributions from initial-state radiation. If we continue the parton shower evolution to low scales we are also able to investigate the transition region between the Durham formalism and the resolved Pomeron formalism [8] which may produce similar final states.
The outline of this paper is as follows. First, we recapitulate the main features of the Durham formalism in section 2. Then we describe the different parts of our implementation in the PYTHIA8 program, starting in section 3 with the reinterpretation of the exclusive luminosity function in terms of the parton shower no-emission probabilities, and followed by a description (section 4) of the spin and colour decomposition of a given partonic state generated by a parton shower from an inclusive hard matrix element. In section 5 we then present some proof-of-concept results for some sample processes before we conclude with a summary and outlook in section 6.

The Durham formalism
Within the Durham model, the amplitude, A, of the central exclusive process in a pp collision, can be written as , (2) where the integration runs over the two-dimensional transverse momentum of the screening gluon q (Fig. 1). The transverse momenta of the outgoing protons are denoted as p 1 and p 2 . The scales Q 2 1 and Q 2 2 are within Durham model [9] defined as i.e. as the smaller of the virtualities of the screening gluon and the fusing gluon momenta related to the particular proton. The proton form factors, F N , are absorbed into offdiagonal unintegrated PDFs f g . These PDFs are assumed to factorise as: where t 1,2 ≈ −p 2 1,2 . The proton form factors in the simplest approach are F N (t) = e bt/2 , with b ≈ 4 GeV −2 . These off-diagonal unintegrated densities depend on the momentum fraction of the fusing (screening) gluon x (x ′ ) and on two scales: the scale of the hard sub-process µ 2 ; and the scale corresponding to the screening gluon transverse momentum Q 2 .
The kinematic regime relevant for CEP is Q/ √ s ∼ x ′ ≪ x ∼ M X / √ s, which allows to integrate out the x ′ -dependency of f g and express them using generalised gluon PDFs H g and the Sudakov factor T M : The Sudakov factor T M describes the probability of no emission from the fusing gluons between scales Q 2 and µ 2 . It resumes singularities from virtual diagrams with soft or collinear emissions up to (modified) nextto-leading logarithmic accuracy and ensures that the integral (2) is finite, as the Sudakov factors exponentially suppress low q contribution.
In this expression both the splitting functions P gg , P qg and the running of α s are in the leading order form. The upper bound of the z integration 2 depends on the mass 2 Usually simply ǫ(k/M X ) = k/M X . of the exclusive system which makes the Sudakov factor T M also M X -dependent as indicated by the subscript, M .
The generalised PDF H g [10] can be approximately calculated from the ordinary parton distribution function of gluon g(x, Q 2 ) using relation: In a much used approximation the generalised PDF, H g , is simply proportional to the conventional one where the constant R g is about 1.3 for LHC energies [10]. The sub-process amplitude M depends on transverse momenta of the fusing gluons q 1 and q 2 and on the gg → X vertex V ab ij which is averaged over identical colour indexes a = b.
The sub-process amplitude is typically known in the helicity basis. In this basis the vertex term q i 1 q j 2 V aa ij takes form: where we have introduced the kinematic spin factors S Jz (q 1 , q 2 ) for future reference with q i 1 = q i − p i 1 and q i 2 = −q i − p i 2 (i = x, y), and the amplitudes, e.g. A ++ , depend on M X , the momenta of outgoing particles, the helicities of outgoing particles as well as the helicities of incoming gluons (here both +1). The amplitudes also depend on the colours of the particles in the sub-process. Here, the amplitudes A are the result of averaging over colour indexes of the incoming partons in such a way that the exclusive system is a colour singlet.
It is useful to note that in the collinear limit, where q 1 = −q 2 = q the kinematic factors are simply: where the azimuthal angle of q was labelled as φ. The state J z = 0 − is trivially zero as a consequence of cross product of two collinear vectors. States |J z | = 2 become zero after integration over φ because the remaining part of (2) is φ independent in the collinear limit.
Beyond such collinear limit also non spin-singlet states contribute to the cross section, but are suppressed as q 2 2 / p 2 1,2 2 ∼ 0.01 with respect to the spin singlet term [5] (if all helicity amplitudes, A, are of the same size).
Note that the formula (2) does not include a soft survival probability which cannot be calculated in a perturbative way. It reduces the CEP cross section at LHC by typically two orders of magnitude and can, e.g., be determined using the eikonal model [10,11,12].
3 Reinterpretation of the exclusive cross section

Prerequisites
In the following discussion it will be beneficial to reformulate the CEP amplitude of the Durham model into more straightforward form and rather work directly with the exclusive cross section σ exc The variables y and M denote the rapidity and mass of the exclusive system and are related to the momenta fractions x 1 and x 2 by formulas The transverse momenta p 1 and p 2 of the scattered protons are assumed to be distributed according the simple one-channel model with the slope of the exponential equal to b ∼ 4 GeV −2 . The integration inside the absolute value is performed over transverse momentum q of the "screening" gluon. The other variables are recognised from the previous section. For completeness the dependency of the off-diagonal generalised PDFs on mass M (via a Sudakov form factor) as well as the arguments of the V aa ij are written explicitly. The kinematics of all outgoing particles of the exclusive system X is denoted by w and dΣ w is the corresponding phase space element. The whole expression is integrated over phase space of these outgoing momenta w which satisfy imposed kinematic cuts.
Inspired by a similar form of the derivative of the exponential function, we factorise the Sudakov factor in expression (6) in front of the bracket which leads to where the newly introduced modified PDF φ M is defined as: In the text below, for simplicity, only the dominant spin singlet part will be considered and the following abbreviation is introduced: where the indexes 1, 2 indicate the dependency of the differential on transverse momenta p 1 and p 2 and the dot represents the scalar product of two-dimensional vectors. In the collinear limit 3 this differential simplifies to −d 2 q/q 4 . Using these notations and assumptions the formula (14) takes form: where q ′ is the transverse momentum of the screening gluon in the complex conjugate amplitude.
In the next step, the exclusive cross section (19) will be expressed as a product of the exclusive luminosity and the exclusive cross section. We define the spinsinglet colour-singlet cross section of hard sub-process 5 as: where the factor 1/4 follows from the probability to have some particular helicity configuration of the incoming gluons and the coefficient 1/64 has an analogous meaning for the colours. The term 1/2M 2 represents the "flux factor". The symmetrisation factor, N s , important if identical particles occur in the final state, is assumed to be incorporated in the amplitudes A, i.e. amplitudes are scaled by 1/ √ N s . Letters λ and j denote all possible helicity and colour configurations of the exclusive system.
The exclusive luminosity which corresponds to the exclusive cross section of the hard sub-process (20) is: In the collinear limit, where p 1 and p 2 are neglected with respect to q and q ′ , the luminosity can be integrated over p 1 and p 2 and over the azimuthal angles of q and q ′ which leads to: The integrand in (23) is fairly flat if considered as a function of 1/q 2 and 1/q ′2 which makes these variables suitable for Monte Carlo integration 4 , especially be- The integrand in (24) is shown in Fig. 2 for ln q 2 as an integration variable. For calculating of the exclusive luminosity the upper limit of integration is set to the factorisation scale µ F , nevertheless the integrand is typically negligible in the high q region. Whereas it dominates for q of 2 − 3 GeV in case of LHC Higgs production and for even smaller values (1 − 2 GeV) in case of di-jet production at Tevatron. Values below the starting scale q 0 = 1 GeV of MMHT2014 LO PDF [14] can be extracted using a backward DGLAP evolution [15]. In reality, it was argued in [16] that for q 0.85 GeV the gluon propagator would be modified by non-perturbative dynamics which effectively suppress such low q contributions. Rather than a sharp cut-off, the damped gluon PDF [5] is used to calculate φ M below q 0 in order to suppress the region of low transverse momentum: The coefficients a 1,2 are chosen in such a way that the function is smooth in q 0 up to the second derivative. Finally, the CEP cross section expressed as a convolution of the exclusive luminosity (24) and the exclusive cross section of the hard sub-process (20) takes from in analogy with the formula formula for the inclusive cross section where the inclusive luminosity L inc is 3.2 Screening gluons in the PYTHIA8 interleaved parton shower The way parton showers are included in PYTHIA8 is through an interleaved process where we normally have three competing processes, which we will denote ISR, FSR and MPI. ISR is an initial-state splitting where one of the incoming partons to the hard sub-process is evolved to lower scale and higher energy fraction, by emitting a parton into the final state; FSR is the final-state splitting of a parton in the final state; while MPI is the appearance of an additional parton-parton interaction. All of these occur at decreasing scale, where the highest scale is given by the kinematics of the hard sub-process. In each step in the shower we then pick a process which has a lower scale than the previous one, and the factorisation property of the no-emission probability means that we can generate one of each of the possible processes independently and simply pick the one which yielded the highest scale in each step.
For simplicity we will here only consider the ISR, concentrating on the initial-state g → gg splittings, and show how we can reinterpret the formula for CEP as an extra process in the interleaved shower, which transforms an inclusive event into an exclusive one.
First, let's consider the inclusive processes with no emission from the space-like shower between scales µ 2 and µ 2 F . The scale µ 2 F is considered as a starting scale of the backward space-like parton shower. The cross section for such processes take form: M is the no-emission probability, quantifying the probability that no extra emission from parton are present between two given scales, if the higher scale is taken as a reference. This no-emission term is used in the backward evolution of the space-like showers in Pythia and, for the case of an incoming gluon, it is defined as: where the sum runs over gluons and all possible flavours of quarks and anti-quarks. It can be shown that these no-emission probabilities are linked to the standard Sudakov form factors by the relation [17,18]: The cross section σ inc , differential in the variable ln µ 2 , corresponding to the scale of the first parton shower emission is: Employing the relation (30) the derivative of the no-emission probabilities can be expressed using the Sudakov factors only: where the newly defined distribution functiong is: 5 The functiong, resembling φ M , depends also on mass M .
It allows to re-express the differential cross section (31) using the standard Sudakov form factors: The inclusive cross section is then simply: where the lower integration limit is assumed to be so small that the Sudakov factor between this limit and µ 2 F is close to zero. We can now apply a similar procedure to the exclusive cross section where the variable µ 2 is now interpreted as the minimal transverse momentum of the screening gluons exchanged during the interaction. Let's define: Then the derivative of the exclusive cross section according to this variable is simply: The ratio of the integrands in eqs. (37) and (35) defines the exclusive probability This means that we now have for each step in the interleaved shower a probability for a given partonic state to become exclusive by the exchange of a screening gluon. The physical picture is the same as in the Durham model, in that a screening gluon with low transverse momentum cannot affect the colour structure of an emission at a higher scale. It also has the nice property that we become somewhat insensitive to the low transverse momentum behaviour of the parton densities in the integral of the exclusive luminosity.
There is, however, a problem with this approach, in that the p exc is very peaked at small µ, making the generation of the exclusive events very inefficient. Although a weighting procedure could be applied, it would be difficult to incorporate into the current framework of PYTHIA8. For the purpose of this paper, we have therefore chosen to implement a simpler procedure.

A simpler approach
Instead of adding the exchange of a screening gluon as an extra process in the interleaved shower, we simply calculate before each shower step, if given state should be made exclusive, using the probability where µ 2 is the scale of the latest emission. We note that the integration in L exc now goes down to very low transverse momenta, but it turns out that the results are the same as in the more complicated approach above.
The modified cross section σ ′i is defined as: where the P i splitting functions in principle could be either initial-or final-state splittings. For modified singlet sub-process cross section the definition is the same, only the formula is corrected for the fact that the incoming gluons are in colour and spin singlet state which, for example, makes H + jet exclusive cross section equal to zero. Note that for calculation of the modified singlet cross section, not only classical matrix element squared but also amplitudes for all possible helicity configuration must be known. The procedure for calculating σ ′ s will be provided in the next section. The full procedure would be to start the generation of an inclusive process in PYTHIA8 and calculate the probability in (39) of that process being exclusive. Then, after each ISR or FSR step in the interleaved shower, we would again check if the current state can be made exclusive by (39). If the event is to be considered to be exclusive we would rearrange the colour flow accordingly and insert the quasi-elastically scattered protons, but let the shower continue without the ISR process.
Note that the soft survival probability, which we have so far left out of the exclusive luminosity function, corresponds exactly to the probability of having no additional multi-parton interactions, so any MPI in the interleaved shower (before or after the event has 8 been made exclusive) will mean that the event will stay inclusive.
To make the generation of exclusive events more efficient we have here decided to simplify the procedure even more. A final-state emission does not modify the parton densities used in the luminosity functions, and although they may affect the exclusive cross section, we have here decided to leave them out and generate them separately. Also the calculation of the soft survival probability by vetoing any exclusive event in case the shower gives a MPI is extremely inefficient, and we instead calculate that separately and simply multiply the inclusive cross section with that factor 6 . It should be noted that using the MPI model in PYTHIA8 presented in [7] has not been carefully investigated. It has the interesting feature that the probability for additional scattering depends on the hardness of the primary scattering, since harder processes have larger overlap (smaller impact parameter). It also has a natural dependence on the collision energy. The downside is that it is sensitive to the soft behaviour of MPI model and may vary strongly between different tunings of the parameters. In this paper we will simply use the default tune in PYTHIA8 and postpone a proper investigation of the procedure to a future publication.
In the end, the procedure will look as follows: 1. Generate the hard sub-process of interest in PYTHIA8, use the standard inclusive cross section. 2. Use only the initial-state shower in PYTHIA8. 3. Before generating the next initial-state emission, make the event exclusive with the probability in (39). 4. If the event is made exclusive, switch off the initialstate cascade, rearrange colours, remove the proton remnants, insert the scattered protons and continue with final-state radiation from the exclusive state. 5. If the event stays inclusive, generate the next initialstate emission (continue with step 3).
As an alternative, to enable detailed studies of the exclusive processes we will below also use a procedure where we study a specific number of initial-state splittings or only initial-state splittings above a certain scale, µ exc , in which case we run the initial state shower without modification to get the desired states and only afterwards decide if the states should become exclusive.
This approach is efficient if the ratio σ ′s σ ′i does not heavily depend on the number of emissions. This is normally the case for higher emission multiplicities in contrast the first few emissions where the interference effects play a role. In the program both approaches are implemented and we use each of them in such a frequency so that the event weights variation is minimal.

Modified luminosity for all possible helicity configurations
Before proceeding to calculate σ ′s in (39) we have to generalise the expression to include all possible helicity combinations contributing to the exclusive production. First, we define four kinds of the "luminosity amplitudes", L (Jz) , for J z = {0 + , 0 − , +2 + , −2 + }, using the notation in (12): For the future use, it is more convenient to work with L directly related to some particular helicity configuration of the partons entering to the hard sub-process. Let's define These relations allow us to rewrite the CEP cross section (14) as: where j and λ denotes the colour state and helicity state of all final state particles of the hard sub-process. The index i denotes the colour of the fusing gluons. Finally, it is possible to formally factorise the cross section formula into process independent luminosity part and the cross section part: where there is an implicit summation in the last expression over all four helicity indices, and the generalised colour singlet cross section σ s incorporates the normalisation factor 1 512M 2 . The origin of this factor is explained below equation (20). Notice, that the cross section formula (45) is a direct generalisation of relation (25), where only the spin singlet component was considered.
A more general expression for the probability p exc of event being exclusive is then where the modified inclusive cross section σ ′i which incorporates the splitting functions is defined by formula (40), whereas the corresponding cross section σ ′s will be derived in the next section. Both luminosities are evaluated at scale of the latest initial-state emission.

Approximation of matrix elements using shower splittings
To calculate the central exclusive cross section, the amplitude for every helicity and colour combination of the studied sub-process must be known, rather than the spin-and colour-averaged sub-process cross section. For 2 → 2 processes the amplitude A λ l1 λr1→λ3λ4 a l1 ar1→x3x4 depends on the helicities of the incoming (λ l1 , λ r1 ) and outgoing particles (λ 3 , λ 4 ) as well as on the colours of the corresponding particles a l1 , a r1 , x 3 , x 4 . The additional dependency on particle momenta is not written out explicitly.
It is useful to introduce the "generalised cross section" σ of the hard sub-process: which is summed over helicities and colours of the final state particles but not over initial-state one. Moreover, the colour and helicity indexes of the incoming particles are in general considered to be different for the amplitude and its complex conjugated. The generalisation of this cross section to 2 → n processes is straightforward. Knowing the generalised cross section, the inclusive cross section takes the form: where the nominal and complex conjugate indexes were put to be equal by means of delta functions.
The colour singlet spin singlet cross section σ S(0 + ) is obtained by an analogous formula (imposing "left" and "right" colours and helicities to be identical): Or using another notation: where the summation over helicities is made explicit and the generalised colour singlet cross section σ s , firstly used in (45), is: The aim of the following sections is to derive an approximative form of the generalised cross section for the case where we have initial-state parton shower splittings from the left and right incoming partons. These emissions are assumed to be strongly ordered in their p ⊥ .
To be specific, let's consider the gluon emission from the "left" incoming parton. The new generalised cross section will then take the form: where the splitting g → gg was considered. Einstein summation convention is employed, in particular it is summed over helicity λ l e1 and colour e 1 of the emitted gluon.
For splittings that includes quarks, the SU(3) structure constants, f abc , must be replaced by the Gell-Mann matrices, T a ij , and the splitting amplitude, P gg , by the corresponding one.
It is obvious that the spin-momentum and colour parts describing emissions factorise and the resulting generalised cross section after n − 1 emission from the left parton and m − 1 emissions from the right parton is given by: where the colour emission tensor T em depends on type and order of the splittings and its indexes correspond to nominal (1 − 3) or adjoint (1 − 8) SU (3) representations depending on the type of the first and last parton in the shower. The space-like emission tensor P em depends on type, order and kinematics of the splittings. Both of these tensors will be discussed in the following.

The colour emission tensor
In case the first parton in the parton shower as well as the parton entering to the hard sub-process are gluons, the general form of the colour emission tensor can be expressed as: where Tr (T x1 T x2 ) = 1 2 δ x1x2 . In our notation the trace prescription is used for gluon colour indexes and delta function for quark colour indexes that means especially that Tr (T x1 T x1 ) = 4 and δ x1x1 = 3. Consequently, any colour emission tensor between gluon and gluon is fully determined by 9 real numbers irrespectively on the parton shower composition between these two gluons. Especially, the colour emission tensor representing no emissions in the parton shower has only one nonzero coefficient A 12 = 4. A colour tensor with one extra gluon emission is related to the original one by linear transformation (provided in appendix Appendix A) applied to the coefficients A ij , B ijk .
In an analogous way, the colour emission tensor can be constructed for the transition between (anti-)quark and (anti-)quark: The case of no emissions will here corresponds to K 3 = 1, K 1,2 = 0 .
The last alternative is the transition between quarks and gluons which is given by where in the last expression the coefficients, C 1,2 are zero if the first parton is a quark and C 1c,2c are zero if the first parton is an anti-quark.
In appendix Appendix A we provide the complete list of linear transformations, necessary for calculating the colour emission tensor with an extra emission in the beginning of the shower. Considering the initialstate parton shower described by the tensor 7 T em p f p h , if the first parton p f is gluon it can be evolved backward to a quark, anti-quark, or a gluon, making this parton the starting one. A quark can be evolved backward to a quark or a gluon and, finally, an anti-quark can be evolved to an anti-quark or a gluon. This gives in total 7 possibilities. The parton p h attached to the hard subprocess can be quark, anti-quark, or gluon making the overall number of possible linear transformation equal to 7 × 3 = 21. In reality, some of these transformations are independent of whether the parton is quark or antiquark which reduces the number of non-identical linear transformations to 13. In the procedure we have developed, the two colour emission tensors are constructed, one for "left" side and one for "right". Before starting the backward parton shower evolution these tensors describe the shower with zero emissions and the particular type (T em gg or T em qq ) is chosen according to the "left" and "right" type of the parton entering to the hard sub-process. After every step in the backward evolution, one of these tensors is modified using the appropriate transformation.

The spin emission tensor
The leading-order amplitudes corresponding to the possible splitting will depend on the helicities of incoming, outgoing and emitted parton, as well as on the momentum fraction z and the azimuthal angle φ of the emission. These can all be found in the literature [19]. Using these splitting amplitudes, the spin emission tensor can be defined, and in the particular case of only one emission, this tensor has a form: which depends on the helicity and "conjugate" helicity of the incoming and outgoing parton. If two emissions are considered this tensor is determined by summing over the intermediate helicity and "conjugate" helicity 7 The p f denotes flavour of parton which initiates the parton shower and p h is the flavour of parton entering to the hard sub-process. Both partons are assumed to be on the "left" or on the "right" side.

Introducing the helicity indexλ which takes integer values between 1 and 4 and incorporates information about helicity and "conjugate" helicity index the form of the last equation
resembles simple matrix multiplication. To add the emission one therefore only need to multiply the current spin tensor (matrix) with matrix corresponding to the particular emission. The different forms of these matrices for all possible splittings provided in appendix Appendix B.

Amplitude definition
Within the process library the amplitude of each process is defined in the colour trace basis, in a form similar to what used for example in MadGraph [20]. For the gg → gg process, which has the most complicated colour topology, the colour trace basis has 6 terms. 8 It means that 6 amplitudes depending on Mandelstam variables s, t, u must be provided for every helicity combination (together 6 × 16 = 96 amplitudes). Fortunately, most of these amplitudes are equal to zero or are identical to each other, e.g. due to parity invariance. The general form of the amplitude decomposition in the colour trace basis is where the B i are the colour basis vectors. The amplitudes A λ i depend on Mandelstam variables and helicities of the particles but not on their colours.
To calculate the generalised cross section,σ, defined above, the product of every two basis vectors M ij must be known: This product does not depend on the final state colours but only on the colours and "conjugate" colours of the incoming particles. In analogy with the procedure used for colour emission tensor, each of these products can be expressed as a linear combination of few colour tensors B α creating basis. 8 The colour basis of, for example, qq → gg process has two terms only.
If the partons entering to the sub-process are gluons then the tensors B α can be one of the following: T ar1 ) and for sub-processes with the incoming quarks or antiquarks the possible colour tensors B α are: . (63) To be able to calculate the exclusive cross section, the amplitudes A λ l1 λr1→λ3λ4 i and the linear decompositions of the M ij into colour tensor basis listed above must be provided for each sub-process (see appendix Appendix C for an example). If the coefficients in the decomposition of M ij are denoted as M ij α then the generalised sub-process cross section,σ, takes the form: and the generalised colour-singlet cross section, σ s , defined in (51), can be calculated (if initial-state radiation is present) using the following formula: The term in the last line does not depend on the colours but only on the helicities and "conjugate" helicities entering into the hard sub-process. To calculate this term, first the contractions between colour emission tensors T em and all members of the B α basis must be calculated. Using these numbers the coefficients M ij (the colour indexes were contracted) are evaluated and consequently the whole term in the last line of (65). It represents 16 values corresponding to all possible helicity combinations. The full generalised colour singlet cross section, σ s , is finally obtained by multiplying by the "left" and "right" spin emission matrices P em l and P em r .

Sample results
In this section we present a few sample results from our implementation of the Durham formalism. We will focus the discussion on the unique feature of our implementation, i.e. possibility of generation the exclusive states with higher particle multiplicities. Currently, the process library includes all hard QCD 2 → 2 processes; Higgs boson production via gg → H; the single Z 0 production via qq → Z 0 ; and two photon production via gg → γγ and qq → γγ. The program is modular, however, and new processes can be easily added. All presented calculations of the CEP cross sections are made for pp collisions at √ s = 13 TeV and are based on MMHT2014 LO PDF [14]. They incorporate hadronisation as well as the final-state radiation. The initialstate shower is evolved down to µ exc = 1.5 GeV if not said otherwise. All predictions incorporate the soft survival probability estimated using veto on MPI. This probability is around 0.06 with little kinematic dependency.

Di-jet production
We start by studying the properties of our Monte Carlo model for di-jet production at the LHC. In contrast to other implementations of the Durham formalism, our program allows for the generation the exclusive di-jet event from any 2 → 2 QCD hard sub-process, as long as the partons which initiate the space-like parton shower are gluons that can be in a colour singlet state. 9 The variable which describes the size of the phase space available for the space-like parton shower is µ exc which is the lowest allowed transverse momentum of the emission. The maximal allowed p ⊥ of an emission is set to be equal to the hard scale of the sub-process, given by the transverse momentum of the leading jet. For ordinary inclusive events the cut-off scale for the initialstate radiation (ISR) in PYTHIA8 is around 2 GeV. In our discussion, we study the events with transverse momenta of the emissions starting at 1.5 GeV.
The inclusion of possible initial-state splittings in our approach will naturally increase the exclusive cross section for di-jet production, and cause a smearing towards low values in the distribution of M 12 /M X . The M 12 /M X observable can be seen as an experimental measure of the "exclusivity" of the particular event. M 12 is here the invariant mass of the two leading jets, and the total mass, M X , of the exclusive system X can, in principle, be calculated from the outgoing protons relative momentum loss, ξ, as √ ξ 1 ξ 2 s. Without any parton showers this ratio equals to 1 on the parton-level. The final-state radiation and hadronisation can smear this distribution, especially if the jet radius of the jet algorithm is small, since a final state parton may radiate outside the jet cone, giving to the smaller value of the invariant di-jet mass M 12 .
We have here used the "anti-k ⊥ " jet algorithm [21] with R = 0.7, a minimum transverse momentum of the jets of 40 GeV, and the absolute value of the pseudorapidity of jets smaller than 2.5. As seen in figure 3 there is indeed a smearing from the final-state radiation (blue curve), but the smearing increases significantly if initial-state radiation is included (black curve). The distribution with initial-state radiation resembles what one would expect form the double Pomeron scattering and the final states generated by these two mechanisms overlap. However, the physical nature of both processes are different since, in DPE, where the Pomeron in the simplest approximation is a gg object, the colour neutralisation of the hard system comes from the Pomeron remnant gluons, while for CEP it is due to an additional gluon exchange. Despite different pictures, the final states could still be indistinguishable, as low-p ⊥ initial state emission on either side of the hard scattering in CEP could look exactly like Pomeron remnants.
The exclusive cross section is less sensitive to the space-like emissions if only the events with, for example M 12 /M X > 0.8 are accepted as is demonstrated in figure 4.
Here, the left plot shows that the di-jet cross section consists of events either with p last ⊥ ∼ 40 GeV where there was typically no space-like emission and p last ⊥ was identified with the hard scale of the process, or events with p last ⊥ ∼ 3 GeV. In this case, there are usually many emissions and p last ⊥ denotes the transverse momentum of the latest one with the smallest p ⊥ .
It can be seen the di-jet cross section differential in the p ⊥ of the last ISR emission peaks for p last ⊥ ∼ 3 GeV and decreases for lower transverse momenta.
More comprehensive picture of the situation provides the two-dimensional plot (Fig. 5), where the correlation of the mass ratio and the p last ⊥ for a particular event is shown. The depletion of the emissions for p last ⊥ slightly below 40 GeV is partially due to the small gg → ggg colour singlet cross section and partially just a statistical effects given by the low probability of having no further emission below such high p ⊥ . The treelevel gg → ggg spin-singlet colour-singlet cross section in the analytic form is provided in [22]. This cross section is zero in the "parton-shower" limit where one of the outgoing gluons has small p T compared to the remaining two andt =û = −ŝ/2, where the Mandelstam variables are derived from the two hardest gluons whereas the softest one is supposed to be part of  is identified with the scale of the sub-process. The total exclusive cross section as a function of the cut-off scale µexc is shown on the right. This scale means that that only events with p last ⊥ > µexc are accepted. As is seen from the right plot the incorporation of the ISR with µexc = 1.5 GeV increases the CEP cross section from 19 pb to 63 pb. The phase space is defined by p jet1,2 ⊥ > 40 GeV, |η jet1,2 | < 2.5, M 12 /M X > 0.8 and ξ 1,2 < 0.03. the shower. Such behaviour agrees with our calculations based on procedure introduced in section 4.
The cut-off parameter µ exc can be understood as a variable which describes the transition between perturbative and non-perturbative region. Not only due to the possible overlap with the double Pomeron exchange process but also because the scale µ exc denotes the minimal allowed p ⊥ of the ISR emission and the p ⊥ of the softest emission is simultaneously the highest allowed transverse momentum of the screening gluon. Choosing small µ exc leads to low p last ⊥ and consequently the main contribution to the exclusive luminosity given by integral (24) stems from small transverse momenta, i.e. smaller than 1 GeV, where the perturbative QCD is not justified [16].
Finally in table 1 we show the contribution to the exclusive cross section from the different possible hard sub-process. Table 1 The table demonstrates how the particular hard sub-processes in PYTHIA8 contribute to the total exclusive cross section of the di-jet production at LHC ( √ s = 13 TeV). The hard processes are defined using the Pythia convention and are accompanied by the Pythia process Id [23]. The letter q denotes any light quark flavour, therefore, e.g. gg → qq represents the sum of gg → uū, gg → dd and gg → ss cross sections. The jets in the di-jet system are required to have p jet1,2 ⊥ > 40 GeV and |η jet1,2 | < 2.5. In addition the leading protons momentum loss ξ must be ξ 1,2 < 0.03. The σ n Em =0 exc are the exclusive cross sections with no initial-state radiation. The σexc and σ M12/M X >0.8 exc are the exclusive cross section with allowed initial-state radiation down to 1.5 GeV; the last one has an additional constrain M 12 /M X > 0.8.

Id
Process  The exclusive di-jet cross section for pp collisions at √ s = 13 TeV double differential in M 12 /M X and ln (p last ⊥ ) 2 . For sake of clarity the vertical axis is denominated directly in p last ⊥ . The phase space is defined by p jet1,2 ⊥ > 40 GeV, |η jet1,2 | < 2.5 and ξ 1,2 < 0.03.
One can see that even with space-like parton showers enabled, the gg → gg sub-process dominates. The second largest cross section is given by the qg → qg process which is forbidden without ISR. Consequently, the fraction of di-jet events where at least one of them is quark-induced with respect to the total exclusive dijet cross section is much higher than ∼ 10 −4 predicted in [2]. This fact makes it problematic to use the CEP as a pure source of gluonic jets.
Within the collinear approximation the gg → qq cross section is predicted to be suppressed as m 2 q /s with respect to the gg → gg cross section. This is well-known consequence of the spin singlet selection rule. It is interesting that without using such collinear approximation the exclusive production of light flavour qq jets is not so heavily suppressed since the |J z | = 2 contribution, absent in collinear case, has a similar size and is quarkmass independent [2]. This effect is even stronger if the ISR is included.
Nevertheless, for higher M 12 /M X the fraction of the heavy flavours jets with respect to the whole CEP's dijet sample is still predicted to be lower compared to the DPE which makes such quantity a vital experimental variable for studying the transition region between CEP and DPE as was first done at the Tevatron [24].
In figure 6 we have tried to compare the results from our CEP program for the distribution in M 12 /M X with data published by the CDF collaboration [24]. The comparison is a bit uncertain as the data has not been corrected to the hadron-level, and the acceptance in different regions of phase space is difficult to disentangle. Nevertheless we have checked that our implementation of the DPE gives results similar to what was published in [24] for normalised distributions.
Our DPE implementation uses diffractive parton densities, as measured e.g. by HERA. Specifically, we will here use the HERA H1 2006 Fit B DPDFs [25] with, for simplicity, the same soft survival probability as for the CEP process, although we are aware that the soft survival probability may very well be different for the DPE process as compared to the CEP one. Technically the DPE simulation is done in Pythia (with the same setting as for the CEP) by colliding two hadrons, the  Fig. 6 The di-jet event counts binned in M 12 /M X variable as measured by the CDF collaboration [24] in pp collisions at energy √ s = 1.96 TeV. These measured event counts (not corrected for the detector related effects) are compared with the DPE contribution (black dashed line) and CEP contribution (red line) and their sum (black solid line). The normalization of these curves is fixed in such a way that the total predicted event count is the same as in data. In figures c and d the CEP contributions were scaled down by a factor of 0.25 compared to their nominal value.
Pomerons, with energies 1 2 ξ 1 √ s and 1 2 ξ 2 √ s and with parton densities described by the HERA DPDFs. 10 In figure 6a and b we show the results of simply adding our CEP generated events for two different selection cuts (two jets above 10 and 25 GeV respectively and no third jet above 5 GeV). Further selection criteria, identical for both phase spaces are given in [24]. We see that the addition of CEP severely overshoots the data in the exclusive region of high M 12 /M X . There are, however, many uncertainties, especially when it comes to the soft survival probability, both for the CEP and DPE contribution. As a demonstration we show in figure 6c and d the effect of introducing a relative normalization factor of 0.25 between the CEP and DPE contribution, which gives a quite reasonable description of the data. We note that our CEP, as expected, contributes quite noticeably also away from the purely exclusive region. A more detailed study of the differences between our new CEP procedure and the DPE one, especially in the regions of the Pomeron remnants in DPE, may result in observables that could further improve the experimental separation between the two processes. 16

Higgs production
The possibility to measure the Higgs boson in the central exclusive production was studied extensively in the last decade [27,28]. The discussion was mainly focused on the dominant decay channel gg → H → bb with a standard model branching ratio of 59%.
The main advantage of this production mechanism is a huge suppression of the irreducible standard model background from gg → bb due to the J z = 0 selection rule in CEP. Furthermore, the scalar nature of the Higgs boson means that the ratio of exclusive to inclusive cross sections is relatively enhanced as compared to the background, 11 as the spin and colours have to match also in the inclusive sub-process. Both these effects improve the signal/background ratio for the Higgs boson production compared to the inclusive production.
The main background to the exclusive Higgs boson production comes from the gg → gg sub-process, which can be substantially suppressed using b-jet tagging techniques. The other experimental challenge is the detection of the scattered protons in the forward detectors in a high pile-up environment where protons from several interactions can simultaneously hit the forward detector within one bunch crossing 12 .
Note, that the cross section of the Higgs boson production in CEP is only around 2 fb, including the calculated soft survival probability of 0.06, which is about four order lower than the inclusive Higgs cross section ∼ 20 pb. The signal event's count is further reduced due to selection criteria and inefficiency of the b-jet tagging. In particular, the QCD background must be suppressed by selecting only high p ⊥ b-jets (comparable to M H /2) because the Higgs boson decay is isotropic whereas the QCD jet production is suppressed at high p ⊥ at least as 1./p 4 ⊥ . We included the Higgs boson production in the process library of our program to study the production rates compared to the background processes. The simulation incorporates the parton showers as well as hadronisation of the resulting partons into "stable" particles, where the particles with lifetime higher than 0.01 mm/c are considered to be stable. The inclusion of initialstate showers have negligible effect on the exclusive Higgs cross section but can substantially increase the s ≈ 0.02. 12 The background protons typically originate from single diffractive excitation. Two such soft single diffractive events together with one inclusive can fake the CEP topology. Fortunately, this kind of experimental background is suppressed for higher masses of the exclusive system (higher ξ).  Fig. 7 The differential distribution of the invariant mass of two leading b-jets for pp CEP at LHC ( √ s = 13 TeV). To enhance the signal fraction the additional cuts p jet1,2 T > 50 GeV and M 12 /M X > 0.9 were applied. The differential cross section stemming from the Higgs decay is given by the red solid curve and the QCD background from gg → bb process is given by the black curve. In addition, the dashed curves indicate the corresponding cross section if initial-state showers are not considered.
gg → bb background and spoil the signal significance (see table 1).
We have simulated Higgs production at the LHC at √ s = 13 TeV. The hard scale is set to be equal to the Higgs mass for the signal, and to p ⊥ of the leading jet for the background, the cut-off for the space-like showers is 1.5 GeV in both cases. To pass the selection cuts, the events are required to contain at least two b-jets with p ⊥ higher than 50 GeV and, additionally, the ratio M 12 /M X must be higher than 0.9. As before, the jets are identified using anti-k ⊥ jet algorithm with R = 0.7 and a jet is tagged as a b-jet if it contains at least one bottom hadron. For now, the kinematics of the scattered protons is not constrained. The result of these calculations is presented in figure 7, where the dotted lines indicate the fraction of the cross section without spacelike emissions. It is obvious that events with space-like emissions play a role only for the gg → bb process and their rate can be probably further reduced using more sophisticated selection techniques. Note that the signal peak is a little bit shifted towards lower values compared to the Higgs mass m H = 125 GeV, due to the fact that sometimes not all produced particles in the hadronisation of the b-quarks are incorporated into the b-jets.
In reality, the forward proton spectrometers installed to ATLAS and CMS have a limited acceptance in ξ, the lowest measurable value of ξ is projected to be around  Fig. 8 The cross section for the CEP production of two b-jets at LHC with p jet1,2 ⊥ > 50 GeV and invariant mass 116 < M 12 < 127 GeV differential in M 12 /M X . The red and black solid lines denote the cross section rising from the Higgs boson decay and the QCD background, respectively. The dashed lines denoted the corresponding cross sections with no initial-state radiation. In addition, the double Pomeron cross sections for Higgs production and the QCD gg → bb process are plotted by the dotted lines.
0.015 which restricts the minimal value of the exclusive system mass to M X = ξ √ s = 195 GeV. This acceptance limit makes the observation of single Higgs boson production without no other activity impossible. On the other hand, there is still a hope of the signal of the Higgs boson accompanied by jets originating from the space-like emissions 13 . To see the size of such cross section we plot the b-jets cross section (both of them must still have p ⊥ > 50 GeV) in the mass window between 116 and 127 GeV where the signal peak is expected. This cross section is shown in figure 8 as a function of the M 12 /M X ratio both for signal and background Monte Carlo sample. The ratio of these cross sections roughly matches the signal/background estimate. It is quite good for M 12 /M X > 0.9 which is the kinematic phase space shown in figure 7 whereas deteriorates for lower values of M 12 /M X .
To reach the acceptance of LHC forward detectors, the mass ratio must be lower than 0.6. Assuming 0.5 < M 12 /M X < 0.6 the signal cross section of the b-jets production is around 0.05 fb 14 and is around 200 times smaller than the QCD background. This small signal cross section and huge background contamination leads to a luminosity of ∼ 60, 000 fb −1 to reach 4-sigma precision. Although the possibility of measure the Higgs pro- 13 Due to the colour singlet nature of Higgs production, at least two emissions are needed. 14 Compare to 0.4 fb for M 12 /M X > 0. 9. duction in this experimental setup is rather academic, our framework allows to determine such cross sections as well as more realistically evaluate the contamination from the QCD background processes.
In figure 8 we also show the corresponding calculation from the DPE process, which becomes significant at low values of M 12 /M X both for the signal and background, but clearly does not give any increase in the significance.

Z 0 production
Considering the acceptance of the forward proton spectrometers of the ATLAS and CMS detectors, which is around 0.015 < ξ < 0.1, the mass of Z 0 resonance is much smaller than the acceptance limit M X ≈ 200 GeV. This makes the study of direct production 15 of the Z 0 within the CEP mechanism even more impossible than the Higgs. Moreover, the Z 0 is produced by qq → Z 0 sub-process, which cannot be handled directly in the standard implementations of the Durham model. However, our model allows for initial-state radiation from the partons entering to the hard sub-process, which can change the identity of incoming quarks to gluons, which can be then treated using the standard Durham exclusive luminosity. To do so, at least one g → qq emission from each side is needed.
Due to the colour singlet nature of the Z 0 and the fusing quarks, there will probably be a non-negligible cross section for no space like emissions and qq CEP luminosity with a screening quark as discussed briefly above. Here, we will make no attempt to evaluate such cross section although our model can, in principle, be extended to cover this production mechanism as well.
The other mechanism for central (semi-)exclusive production of a Z 0 is through DPE. To estimate such a cross section we use the procedure described in section 5.1, where we again we assume that the DPE soft survival probability is the same as in CEP. Contrary to the CEP where the M X mass is higher than M Z mostly due to space like emissions, in DPE both the space like emissions and the Pomeron remnants contribute to the mass.
We will look at semi-exclusive Z 0 → µ − µ + production at the LHC at √ s = 13 TeV, requiring a minimum transverse momentum of 30 GeV for the muons in the pseudorapidity region of |η| < 2.5. Both quasielastically scattered protons are required to have 0.015 < ξ 1,2 < 0.1 in accordance to the acceptance of forward proton spectrometers. The electroweak process qq → µμ includes both Z 0 exchange and γ exchange as well as the interference terms. We find that the resulting DPE cross section is about ten times higher than for CEP. Quantitatively the CEP cross section is around 3.5 fb compared to 40 fb for DPE. The Z 0 can be produced also via the exclusive photoproduction. The predicted cross section for this process, including Z 0 → µ − µ + branching ratio, is, however, about 0.3 fb [29] and would be even smaller if the selection criteria for muons p T and pseudorapidity had been applied.
The shapes of the M µµ and M µµ /M X distributions are compared in figure 9. The shapes of M µµ /M X are rather similar for both processes, with the DPE curve somewhat shifted towards lower values as compared to the CEP one which prefers more "exclusive" configurations.

Conclusions
In this paper we introduced a new Monte Carlo implementation of the Durham formalism to calculate the central exclusive processes in pp and pp collisions. Our model is based on PYTHIA8 generator, and naturally incorporates partonic showers and hadronisation, as well as multi-parton interactions.
The main advantage of our implementation is the possibility to study the effects of initial-state parton radiation on CEPs. This is done by allowing that any inclusively produced sub-process is converted to an exclusive at any stage in the shower. To do this we have implemented a colour and spin decomposition of the initial-state shower in PYTHIA8 which, together with a similarly decomposed (user supplied) matrix element, can be used to determine the probability that a given partonic state can be exclusive.
We have shown that this way of approximating higher jet multiplicities gives rise to to new, non-trivial, physical consequences. In particular, for exclusive di-jet production, it leads to event topologies with medium values of M 12 /M X which naturally fill the gap between double Pomeron exchange and pure central exclusive production. Moreover, the incorporation of the parton showers enables the generation of quark-initiated processes such as Z 0 production.
All predicted cross sections depend on the parameter µ exc , the scale related the transition between the perturbative and non-perturbative region in the parton shower. The actual value of µ exc will have to be determined from experiment. For the time being we set its value equal to 1.5 GeV.
The cross sections also depend on the soft survival probability used. Here we have used the MPI model in PYTHIA8 to simply estimate the probability of having no additional scatterings, equating this to the soft survival probability. Although this procedure was suggested long ago, it has not been properly investigated, and we intend to return with a detailed study of this model in a future publication.
Currently, the program process library includes QCD 2 → 2 processes, H production, Z 0 production and γγ production, but it can be easily extended. In particular, it would be interesting to add production of vector mesons (ρ, φ, . . . ) and/or quarkonia χ c,b . These processes have large cross sections which make them experimentally accessible even at low luminosities.
Our framework to treat colour and spin states within the partonic shower is rather general and can, in principle, be extended to simulate the central exclusive processes initiated by qq fusion, in addition to standard gg-initiated processes. Here a screening quark rather than screening gluon is exchanged to cancel the colour flow. Such processes would be especially interesting for e.g. central exclusive Z 0 production.
It should also be possible to extend our treatment of the colour and spin structure of the parton showers to treat final-state splittings. This would give an additional way of studying approximate higher order effects in the hard sub-process matrix elements.
These and other possible improvements will be discussed in a future publication.  Fig. 9 The differential cross section of pp → p + µ − µ + X ′ + p process as a function of the invariant mass of the µ − µ + pair (left plot) and of the ratio Mµµ/M X (right plot). The black lines show the CEP contribution whereas the DPE result is given by the red lines. The DPE cross sections are normalised by a factor of 0.1. In the CEP calculations the µexc cut-off parameter is set to 1.5 GeV.