Hard diffraction in photoproduction with Pythia 8

Abstract

We present a new framework for modeling hard diffractive events in photoproduction, implemented in the general purpose event generator Pythia 8. The model is an extension of the model for hard diffraction with dynamical gap survival in $\mathrm {p}\mathrm {p}$ and $\mathrm {p}\overline{\mathrm {p}}$ collisions proposed in 2015, now also allowing for other beam types. It thus relies on several existing ideas: the Ingelman–Schlein approach, the framework for multiparton interactions and the recently developed framework for photoproduction in $\gamma \mathrm {p}$, $\gamma \gamma$, $\mathrm {e}\mathrm {p}$ and $\mathrm {e}^+\mathrm {e}^-$ collisions. The model proposes an explanation for the observed factorization breaking in photoproduced diffractive dijet events at HERA, showing an overall good agreement with data. The model is also applicable to ultraperipheral collisions with $\mathrm {p}\mathrm {p}$ and $\mathrm {p}{\mathrm {Pb}}$ beams, and predictions are made for such events at the LHC.

1 Introduction

Diffractive excitations represent large fractions of the total cross section in a wide range of collisions. A part of these has been seen to have a hard scale, as in e.g. the case of diffractive dijet production. These hard diffractive events allow for a perturbative calculation of the scattering subprocess, but still require some phenomenological modeling. This includes modeling of the Pomeron, expected to be responsible for the color-neutral momentum transfer between the beam and the diffractive system X. In the framework of collinear factorization, a diffractive parton distribution function (dPDF) may be defined. This can further be factorized into a Pomeron flux and a PDF, describing the flux of Pomerons from the beam and the parton density within the Pomeron, respectively.

Here we focus mainly on photoproduced diffractive dijets in $\mathrm {e}\mathrm {p}$ collisions. This scattering process can be separated into different subsystems, visualized in Fig. 1. The initial state consists of an electron and a proton, with the former radiating off a (virtual) photon. If the photon is highly virtual, we are in the range of deep inelastic scattering (DIS), while a photon with low enough virtuality can be considered (quasi-)real. This is the photoproduction regime. No clear distinction between the two regimes exists, however, and photons of intermediate virtuality require careful consideration to avoid double-counting. A special feature in the photoproduction regime is that there is a non-negligible probability for the photon to fluctuate into a hadronic state. These resolved photons open up for all possible hadron–hadron processes, including diffractive ones.

The next subsystem shown in Fig. 1 is the photon–proton scattering system. Here, diffraction could in principle occur on both sides if the photon is resolved. In direct photoproduction (and in DIS) the diffractive system can only be present on the photon side, as no Pomeron flux can be defined for point-like photons. In this article the emphasis will be on Pomeron emission from the proton.

The final subsystem is the hard scattering generated inside the diffractive system X. For direct photoproduction (and DIS) this includes the photon as an incoming parton, see Fig. 1a. In the resolved case, Fig. 1b, a parton is extracted from the hadronic photon, which then proceeds to initiate the hard scattering along with a parton extracted from the Pomeron. In both cases a beam remnant is left behind from the Pomeron, while resolved photoproduction also gives rise to a beam remnant from the hadronic photon. Multiple scatterings or multiparton interactions (MPIs) are expected between the remnants, but also in the larger photon–proton system. The particles produced by the latter type of MPIs may destroy the diffractive signature, the rapidity gap between the diffractive system and the elastically scattered proton (or meson, depending on the side of the diffractive system).

Fig. 1

Leading order Feynman diagrams for diffractive dijet production with photons in $\mathrm {e}\mathrm {p}$ collisions. Either the photon participates directly in the hard scattering matrix element (a) or a parton from the resolved photon participates (b)

The model for photoproduced diffractive dijets presented here is based on the general-purpose event generator Pythia 8 [1]. It combines the existing frameworks for photoproduction and hard diffraction, the latter originally introduced for purely hadronic collisions. The new model thus allows for event generation of photon-induced hard diffraction with different beam configurations. The model is highly dependent on the components of Pythia 8. The relevant ones – the model for MPIs, photoproduction and hard diffraction – are described in the following sections.

The first measurements of diffractive dijets was done by the UA8 experiment at the S$\mathrm {p}\overline{\mathrm {p}}$S collider at CERN [2]. Later on, similar events have been observed in $\mathrm {e}\mathrm {p}$ collisions at HERA [3], in $\mathrm {p}\overline{\mathrm {p}}$ collisions at the Tevatron [4], and nowadays also in $\mathrm {p}\mathrm {p}$ collisions at the LHC [5]. Similarly, diffractively produced $W^{\pm }$ and $Z^0$ bosons have been observed at the Tevatron [6]. All of these processes are expected to be calculable within a perturbative framework, such as the Ingelman–Schlein picture [7]. A model for such hard diffractive events was included in Pythia 8 [8], based on the Ingelman–Schlein approach and the rapidity gap survival idea of Bjorken [9]. The model proposed an explanation of the observed factorization breaking in hard diffractive $\mathrm {p}\overline{\mathrm {p}}$ collisions – the observation that with the Pomeron PDFs and fluxes derived from HERA DIS data, the factorization-based calculation was an order of magnitude above the measurement. The suppression factor required on top of the dPDF-based calculation, was dynamically generated by requiring no additional MPIs in the $\mathrm {p}\mathrm {p}$ (or $\mathrm {p}\overline{\mathrm {p}}$) system. The model predicted production rates in agreement with $\mathrm {p}\mathrm {p}$ and $\mathrm {p}\overline{\mathrm {p}}$ measurements, albeit some differential distributions did show room for improvement when comparing to Tevatron data. The latest preliminary analysis on diffractive dijets by CMS [10] finds a very good agreement between the model and data in all differential distributions.

First evidence of factorization breaking for diffractive dijets in $\mathrm {e}\mathrm {p}$ collisions was observed by an H1 measurement [11], where a suppression factor of 0.6 was required to describe the dijet data in the photoproduction region, whereas the analysis for the DIS region was, by construction, well described by the factorization-based model without a corresponding suppression factor. Advances in the formulation of the dPDFs improved the description of data in the DIS regime, but the discrepancies remained in the photoproduction limit. Several analyses have been performed by H1 and ZEUS for diffractive dijet production [12,13,14,15,16], all requiring a suppression factor between 0.5 and 0.9 in order for the factorization-based calculations to describe data.

The extension of the hard diffraction model in this article, to collisions with (intermediate) photons, makes it possible to explain the factorization-breaking in the photoproduction regime. The model is also applicable to the DIS regime, but here no further suppression is added since the highly virtual photons do not have any partonic structure that would give rise to the MPIs. Furthermore, the framework can also be applied to diffractive photoproduction in purely hadronic collisions, usually referred to as ultra-peripheral collisions (UPCs) [17]. The model predicts a substantial suppression for diffractive dijets in UPCs at the LHC.

The article is structured as follows: After the introduction in Sect. 1, we briefly describe in Sect. 2 the event generation procedure in Pythia 8. We then proceed in Sect. 3 to the photoproduction framework available in Pythia 8 and continue to a short description of the hard diffraction model in Sect. 4. We present results with our model compared to data from HERA on diffractive dijets in photoproduction in Sect. 5, and show some predictions for photoproduction in UPCs at the LHC in Sect. 6. We end with Sect. 7 where we summarize our work and provide an outlook for further studies.

2 Event generation with Pythia 8

Recently, Pythia 8 has undergone a drastic expansion. Where the earlier version, Pythia 6 [18], was designed to accommodate several types of collisions (lepton–lepton, hadron–hadron and lepton–hadron, excluding nuclei), the rewrite to C++ focused mainly on the hadronic physics at the Tevatron and the LHC. While the LHC will run for years to come, there are several future collider projects under consideration. A common feature between the projected colliders is that they will be using lepton beams either primarily (linear $\mathrm {e}^+\mathrm {e}^-$ colliders: CLIC and ILC [19, 20] or electron–ion collider (EIC) [21]), or as a first phase towards a hadronic collider (FCC [22, 23]). To enable studies related to these future colliders, Pythia 8 has been extended to handle many processes involving lepton beams. Another major facility has been the extension from $\mathrm {p}\mathrm {p}$ to $\mathrm {p}\mathrm {A}$ and $\mathrm {A}\mathrm {A}$ collisions with the inclusion of the Angantyr model for heavy ion collisions [24]. Combining the heavy-ion machinery with the recent developments related to lepton beams will also allow simulations of $\mathrm {e}\mathrm {A}$ collisions and ultra-peripheral $\mathrm {A}\mathrm {A}$ collisions. Work in this direction has been started within the Pythia collaboration.

The Pythia 6 description of lepton–lepton and lepton–hadron collisions included a sophisticated model for merging of the DIS regime (high-virtuality photons) and the photoproduction regime (low-virtuality photons) [25]. This, however, created upwards of 25 different event classes, each of which had to be set up differently. The model for the transition from photoproduction to DIS turned out not to agree so well with data, and the division of the different event classes was somewhat artificial. The aim for the Pythia 8 implementation of these processes has been to reduce the number of hard-coded event classes and increase robustness. The present framework, however, does not yet include a smooth merging of the high- and low-virtuality events and therefore the events with intermediate virtualities are not addressed. Work towards such a combined framework is currently ongoing. In addition, there is progress towards improving the parton showers for DIS events (see e.g. [26, 27]). In this paper we focus on the photoproduction regime, which is mature and well tested for hard-process events with virtuality $\lesssim 1~\mathrm {GeV}$ against LEP and HERA data [28,29,30].

The generation of non-diffractive (ND) $\mathrm {p}\mathrm {p}$ or $\mathrm {p}\overline{\mathrm {p}}$ events proceeds with the following steps. First, the incoming beams are set up with (possible) PDFs at a given (user-defined) energy. Then the hard scattering of interest is generated based on the matrix element (ME) of the process and the PDFs. The generated partonic system is then evolved with a parton shower (PS), in Pythia 8 using the interleaved evolution of both initial and final state showers (ISR, FSR) [31] and MPIs [32]. The splitting probabilities for the FSR and ISR are obtained from the standard collinear DGLAP evolution equations. The ISR probabilities also depend on the PDFs of the incoming beams, as the evolution is backwards from a high scale, set by the hard process, to a lower scale. Similarly, the MPI probabilities depend on the PDFs of the incoming beams, and these have to be adjusted whenever an MPI has removed a parton from the beam. Colour reconnection (CR) is allowed after the evolution to mimic the finite-color effects that are not taken into account in the infinite-color PS. After the partonic evolution, a minimal number of partons are added as beam remnants in order to conserve color, flavor and the total momentum of the event. Lastly, the generated partons are hadronized using the Lund string model [33] along with decays of unstable particles.

In $\mathrm {e}\mathrm {p}$ events, Pythia 8 operates with two regimes: the DIS regime, where the electron emits a highly virtual photon ($Q^2\gg 1$ $\hbox {GeV}^2$), and the photoproduction regime, where the photon is (quasi-)real ($Q^2\lesssim 1$ $\hbox {GeV}^2$). Currently no description is available for intermediate-virtuality photons. In DIS events, the hard scattering occurs between the incoming lepton and a parton from the hadron beam by an exchange of a virtual photon (or another EW boson). The photon can thus be considered devoid of any internal structure. In the photoproduction regime, the photon flux can be factorized from the hard scattering, such that the intermediate photon can be regarded as a particle initiating the hard scattering. In this regime, both point-like and hadron-like states of the photon occur. This significantly increases the complexity of the event generation, thus the photoproduction regime is thoroughly described in the next section.

3 The photoproduction framework

The (quasi-)real photon contains a point-like, direct part without substructure as well as a hadron-like part with internal structure. The latter part, the resolved photon, dominates the total cross section of the physical photon. The total cross section is expected to contain all types of hadronic collisions, including elastic (el), single- and double diffractive (SD, DD) and inelastic ND collisions. The ND collisions contain both hard and soft events, where the former can be calculated perturbatively, while the latter are modeled using the MPI framework in Pythia 8 [34]. Elastic and diffractive collisions require a phenomenological model for the hadronic photon.

The ND processes were first introduced in Pythia 8.215 [30], with a cross section given as a fraction of the total cross section, $\sigma _{\mathrm {ND}}=f\sigma _{\mathrm {tot}}$, $f<1$. The framework for photoproduction has since been expanded to include all soft QCD processes using the Schuler–Sjöstrand model [35] in Pythia 8.235, and with this the cross sections for each of the event classes is calculated separately. The full description of these event classes is postponed to a forthcoming paper [30], as we here concentrate on diffractive processes with a hard scale. Between the two versions, the $\gamma \mathrm {p}$ and $\gamma \gamma$ frameworks were extended to $\mathrm {e}\mathrm {p}$ and $\mathrm {e}^+\mathrm {e}^-$ by the introduction of a photon flux within a lepton, now giving a complete description of all photoproduction events in $\gamma \mathrm {p}$, $\gamma \gamma$, $\mathrm {e}\mathrm {p}$ and $\mathrm {e}^+\mathrm {e}^-$ collisions in the latest release, 8.240. Furthermore, an option to provide an external photon flux has been included, allowing the user to study photoproduction also in UPCs, where the virtuality of the intermediate photon is always small and thus the photoproduction framework directly applicable. An internal setup for these cases is under way.

The resolved photon is usually split into two: one describing a fluctuation of the photon into a low-mass meson and the other describing a fluctuation into a $\mathrm {q}\overline{\mathrm {q}}$ pair of higher virtuality. The former is usually treated according to a vector-meson dominance (VMD) model [36, 37], where the photon is a superposition of the lightest vector mesons (usually $\rho$, $\omega$ and $\phi$), whereas the latter, the anomalous part of the photon, is treated as “the remainder”, $\sigma _{\mathrm {anom}}=\sigma _{\mathrm {tot}} - \sigma _{\mathrm {direct}} - \sigma _{\mathrm {VMD}}$. A generalization of the VMD exists (the GVMD model) which takes into account also higher-mass mesons with the same quantum numbers as photons [38]. Note, however, that if the resonances are broad and closely spaced, they would look like a smooth continuum.

The event generation for the direct photons begins by sampling the hard scattering between the incoming photon and a parton (or another direct photon in case of $\gamma \gamma$), e.g. $\mathrm {q}\gamma \rightarrow \mathrm {q}\mathrm {g}$. The subsequent parton-shower generation always include FSR and in $\gamma \mathrm {p}$ case also ISR for the hadronic beam. The whole photon momentum goes into the hard process, $x_{\gamma }\sim 1$, as direct photons do not have any internal structure. Hence there is no energy left for MPIs and no photon remnant is left behind. The hadronization is then performed with the Lund string model as usual.

For resolved photons, a model for the partonic content of the hadronic photon, the photon PDF, needs to be taken into account. This PDF includes both the VMD and the anomalous contributions, the latter being calculable within perturbative QCD, the former requiring a non-perturbative input. As in the case of protons, the non-perturbative input is fixed in a global QCD analysis using experimental data. There are several PDF analyses available for photons [39,40,41,42] using mainly data from LEP, but some also exploiting HERA data to constrain the gluonic part of the PDF [43]. Ideally one would have a PDF for each of the VMD states, in practice one uses the same parametrization for all – or approximates these with pion PDFs.

After the setup of the photon PDFs, the hard collision kinematics has to be chosen. Here, a parton from the photon PDF initiates the hard process, carrying a fraction of the photon momentum, $x_{i}<1$, with parton i being extracted from the photon. Thus energy is still available in the fluctuation after the initial hard process, opening up for additional MPIs along with ISR and FSR in the subsequent evolution. As with other hadronic processes, a remnant is left behind, with its structure being derived from the flavor content of the original meson or $\mathrm {q}\overline{\mathrm {q}}$ state and the kicked-out partons.

As in pp collisions, the PS splitting probabilities with resolved photons are based on the DGLAP equations. The DGLAP equation governing the scale evolution of resolved photon PDFs can be written as [44]

$$\begin{aligned} \frac{\partial f_{i/\gamma }(x_{i}, Q^2)}{\partial \log (Q^2)}&= \frac{\alpha _{\mathrm {em}}(Q^2)}{2\pi }e_i^2P_{i\gamma }(x_{i}) \nonumber \\&\quad +\frac{\alpha _{\mathrm {s}}(Q^2)}{2\pi }\sum _j\int _{x_i}^1\frac{\mathrm {d}z}{z} P_{ij}(z)f_{j/\gamma }\left( \frac{x_{i}}{z}, Q^2\right) , \end{aligned}$$

(1)

where $f_{i(j)/\gamma }$ corresponds to the PDF of the photon, $x_{i}$ the fractional momenta of the photon carried by the parton i, $\alpha _{\mathrm {em}}$, $\alpha _{\mathrm {s}}$ the electromagnetic and strong couplings, $e_i$ the charge of parton i and $P_{ij}$, $P_{i\gamma }$ the DGLAP and $\gamma \rightarrow \mathrm {q}\overline{\mathrm {q}}$ splitting kernels, respectively. The term proportional to $P_{i\gamma }$ gives rise to the anomalous part of the photon PDF. In Pythia 8 the separation into VMD and anomalous contributions is not explicitly performed. By the backwards evolution of ISR, however, a resolved parton can be traced back to the original photon by a $\gamma \rightarrow \mathrm {q}\overline{\mathrm {q}}$ branching at some scale $Q^2$. Post facto, an event where this happens for $Q^2>Q_0^2$ can then be associated with an anomalous photon state, and where not with a VMD state. The dividing scale $Q_0$ is arbitrary to some extent, but would be of the order of the $\rho ^0$-meson mass. In the interleaved evolution of the parton showers and MPIs, additional MPIs and ISR splittings on the photon side become impossible below the scale where the photon became unresolved. This reduces the average number of MPIs for resolved photons compared to hadrons, and therefore has an impact also for the hard diffraction model as discussed in Sect. 3.1.

3.1 MPIs with photons

When the photon becomes resolved it is possible to have several partonic interactions in the same event. MPIs in Pythia 8 are generated according to the leading-order (LO) QCD cross sections, albeit being regularized by introducing a screening parameter $p_{\perp 0}$ [32],

$$\begin{aligned} \frac{\mathrm {d}\sigma }{\mathrm {d}p_{\perp }^2}\sim \frac{\alpha _{\mathrm {s}}^2(p_{\perp }^2)}{p_{\perp }^4} \rightarrow \frac{\alpha _{\mathrm {s}}^2(p_{\perp 0}^2+p_{\perp }^2)}{(p_{\perp 0}^2+p_{\perp }^2)^2}. \end{aligned}$$

(2)

Note here that $p_{\perp 0}$ can be related to the size d of the colliding objects, $p_{\perp 0}\sim 1/d$, thus a different value of the screening parameter could be motivated if the photon has a different size than the proton. Further, one could imagine working with different matter profiles for both the proton and the photon, and possibly also for each of the components of the photon. For now the shape is kept common for all systems, but possibly with different scale factors, i.e. average radii.

The screening parameter is allowed to vary with center-of-mass energy $\sqrt{s}$,

$$\begin{aligned} p_{\perp 0}(\sqrt{s})=p_{\perp 0}^{\mathrm {ref}} \left( \frac{\sqrt{s}}{\sqrt{s_{\mathrm {ref}}}}\right) ^p, \end{aligned}$$

(3)

with $p_{\perp 0}^{\mathrm {ref}}$, p tunable parameters and $\sqrt{s_{\mathrm {ref}}}$ a reference scale. Thus both the parameters from the matter profile and the parameters related to $p_{\perp 0}$ require input from data. These parameters can be fixed by a global tune, with the Monash tune [45] being the current default. The MPI parameters in this tune, however, are derived using only data from $\mathrm {p}\mathrm {p}$ and $\mathrm {p}\overline{\mathrm {p}}$ collisions. As the partonic structure and matter profile of resolved photons can be very different from that of protons, the values for the MPI parameters should be revisited for $\gamma \gamma$ and $\gamma \mathrm {p}$ collisions. The limitation is that there are only a few data sets sensitive to the MPIs available for these processes, and therefore it is not possible to perform a global retune for all the relevant parameters. Thus we have chosen to use the same form of the impact-parameter profile as for protons and study only the $p_{\perp 0}$ parameters (which allow for different scale factors).

For $\gamma \gamma$ collisions, LEP data is available for charged-hadron $p_{\perp }$ spectra in different $W_{\gamma \gamma }$ bins, allowing studies of the energy dependence of $p_{\perp 0}$ as shown in [28]. In the $\gamma \mathrm {p}$ case the HERA data for charged-hadron production is averaged over a rather narrow $W_{\gamma \mathrm {p}}$ bin. Hence a similar study of the energy dependence is not possible for $\gamma \mathrm {p}$, and it becomes necessary to assume the same energy dependence for $p_{\perp 0}$ in $\gamma \mathrm {p}$ as for $\mathrm {p}\mathrm {p}$ collisions. The value of the $p_{\perp 0}$-parameter, however, can be retuned with the available data. As discussed in [29] a good description of the H1 data from HERA can be obtained with a slightly larger $p_{\perp 0}^{\mathrm {ref}}$ in $\gamma \mathrm {p}$ than what is used in the $\mathrm {p}\mathrm {p}$ tune, $p_{\perp 0}^{\mathrm {ref}}(\gamma \mathrm {p}) = 3.00~\mathrm {GeV}$ versus $p_{\perp 0}^{\mathrm {ref}}(\mathrm {p}\mathrm {p})=2.28~\mathrm {GeV}$. Thus the photon-tune is consistent with a smaller size of the photon, i.e. that the photon does not quite reach a typical hadron size during its fluctuation.

The rule of thumb is that a larger screening parameter gives less MPI activity in an event, thus a smaller probability for MPIs with resolved photons is expected compared to proton–proton collisions. As the model for hard diffraction is highly dependent on the MPI framework, we expect that the increased screening parameter gives less gap-suppression in photoproduction than what was found in the proton–proton study. This is simply because there is a larger probability for the event to have no additional MPIs when the $p_{\perp 0}^{\mathrm {ref}}$-value is larger. Furthermore, since the ISR splittings may collapse the resolved photon into an unresolved state and, by construction, the direct-photon induced processes do not give rise to additional interactions, the role of MPIs is suppressed for photoproduction compared to purely hadronic collisions. Also, the invariant mass of the photon–proton system in the photoproduction data from HERA is typically an order of magnitude smaller than that in previously considered (anti-)proton–proton data, which further reduces the probability for MPIs. Anticipating results to be shown below, this is in accordance with what is seen in diffractive dijet production at HERA, where the suppression factor is much smaller than that at the Tevatron.

3.2 Photon flux in different beam configurations

In the photoproduction regime one can factorize the flux of photons from the hard-process cross section. For lepton beams a virtuality-dependent flux is used,

$$\begin{aligned} f_{\gamma /\mathrm {e}}(x,Q^2) = \frac{\alpha _{\mathrm {em}}}{2 \pi } \frac{1 + (1-x)^2}{x} \frac{1}{Q^2}, \end{aligned}$$

(4)

where x is the momentum fraction of the photon w.r.t. the lepton. Integration from the kinematically allowed minimum virtuality up to the maximum $Q^2_{\mathrm {max}}$ allowed by the photoproduction framework, yields the well-known Weizsäcker–Williams flux [46, 47]

$$\begin{aligned} f_{\gamma /\mathrm {e}}(x) = \frac{\alpha _{\mathrm {em}}}{2 \pi } \frac{1 + (1-x)^2}{x} \log \left[ \frac{Q^2_{\mathrm {max}} (1-x)}{m_{\mathrm {e}}^2 x^2} \right] , \end{aligned}$$

(5)

where $m_{\mathrm {e}}$ is the mass of the lepton.

In $\mathrm {p}\mathrm {p}$ collisions the electric form factor arising from the finite size of the proton, or equivalently that the proton should not break up by the photon emission recoil, needs to be taken into account. A good approximation of a $Q^2$-differential flux is given by

$$\begin{aligned} f_{\gamma /\mathrm {p}}(x,Q^2) = \frac{\alpha _{\mathrm {em}}}{2 \pi } \frac{1 + (1-x)^2}{x} \frac{1}{Q^2} \frac{1}{(1+Q^2/Q^2_0)^4}, \end{aligned}$$

(6)

where $Q^2_0 = 0.71~\mathrm {GeV}^2$. Integration over the virtuality provides the flux derived by Drees and Zeppenfeld [48],

$$\begin{aligned} f_{\gamma /\mathrm {p}}(x)&= \frac{\alpha _{\mathrm {em}}}{2 \pi } \frac{1 + (1-x)^2}{x} \nonumber \\&\quad \times \left[ \log (A) - \frac{11}{6} + \frac{3}{A} - \frac{3}{2A^2} + \frac{1}{3A^3} \right] , \end{aligned}$$

(7)

where $A = 1 + Q_0^2/Q_{\mathrm {min}}^2$ and $Q_{\mathrm {min}}^2$ is the minimum scale limited by the kinematics of a photon emission. Due to the form factor the photon flux drops rapidly with increasing virtuality and becomes negligible already at $Q^2\sim 2~\mathrm {GeV}^2$. This ensures that the photons from protons are well within the photoproduction regime and there is no need to introduce any cut on maximal photon virtuality.

In case of heavy ions it is more convenient to work in impact-parameter space. The size of a heavy nucleus is a better defined quantity than it is for protons, so the impact parameter b of the collision can be used to reject the events where additional hadronic interactions would overwhelm the electromagnetic interaction. Simply rejecting the events for which the minimal impact parameter, $b_{\mathrm {min}}$, is smaller than the sum of the radii of the colliding nuclei (or colliding hadron and nucleus for $\mathrm {p}\mathrm {A}$) provides a b-integrated flux,

$$\begin{aligned} f_{\gamma /\mathrm {A}}(x)&= \frac{\alpha _{\mathrm {em}}Z^2}{\pi x} \nonumber \\&\quad \times \left[ 2 \xi K_1(\xi ) K_0(\xi ) - \xi ^2(K^2_1(\xi ) - K^2_0(\xi ))\right] , \end{aligned}$$

(8)

where Z is the charge of the emitting nucleus, $K_i$ are the modified Bessel functions of the second kind and $\xi =b_{\mathrm {min}} \,x\, m_N$, where x is a per-nucleon energy fraction and $m_N$ a per-nucleon mass. The downside of working in the impact-parameter space is that the virtuality cannot be sampled according to the flux, as virtuality and impact parameter are conjugate variables. For heavy-ions, however, the maximal virtuality is very small (of the order of $60~\mathrm {MeV}$ [17]), and can be safely neglected for the considered applications. The different photon fluxes are shown in Fig. 2.

Fig. 2

The photon fluxes used for different beam types. Here $f_{\gamma /b}$ is the photon flux obtained from the beam b

When extending the photoproduction regime from pure photon-induced processes to collisions where the photon is emitted by a beam particle, some additions are needed. In direct photoproduction, the partonic processes can be generated by using the photon flux directly in the factorized cross-section formula, similar to what is done with the PDFs in a usual hadronic collision. In resolved photoproduction, a PDF for the partons from the photons emitted from the beam particle is needed. This can be found by convoluting the photon flux from the beam particle b, $f_{\gamma /b}(x)$, with the photon PDFs, $f_{i/\gamma }(x_{\gamma },Q^2)$, where $Q^2$ refers to the scale at which the resolved photon is probed. This scale can be linked to the scale of the hard(est) process, e.g. the $p_{\perp }$ of the leading jet in jet-production processes. The convolution yields

$$\begin{aligned} x_{i}f_{i/b}(x_{i},Q^2) = \int _{x_{i}}^1 \frac{\mathrm {d}x}{x} x f_{\gamma /b}(x) \frac{x_i}{x}f_{i/\gamma }\left( \frac{x_i}{x},Q^2\right) , \end{aligned}$$

(9)

with $x_{i}$ the energy fraction of beam particle momentum carried by parton i and x the energy fraction of the photon w.r.t. the beam. In practice the intermediate photon kinematics is sampled according to the appropriate flux during the event generation, thus taking care of the convolution on the fly.

4 Hard diffraction in Pythia 8

The Pythia model for hard diffractive events in $\mathrm {p}\mathrm {p}$ collisions was introduced as an explanation for the factorization breaking between diffractive DIS at HERA and the Tevatron [8]. The model can be applied to any process with sufficiently hard scales, including production of dijets, $Z^0,W^{\pm },H$ etc. It begins with the Ingelman–Schlein picture, where the diffractive cross section factorizes into a Pomeron-particle cross section and a Pomeron flux. Based on this ansatz a tentative probability for diffraction is defined as the ratio of diffractive PDF (dPDF) to inclusive PDF, as it is assumed that the proton PDF can be split into a diffractive and a non-diffractive part,

$$\begin{aligned} f_{i/\mathrm {p}}(x_i, Q^2)&= f_{i/\mathrm {p}}^{\mathrm {ND}}(x_i, Q^2) + f_{i/\mathrm {p}}^{\mathrm {D}}(x_i, Q^2)\nonumber ,\\ f_{i/\mathrm {p}}^{\mathrm {D}}(x_i, Q^2)&= \int _{x_i}^1 \frac{\mathrm {d}x_{\mathbb {P}}}{x_{\mathbb {P}}}f_{\mathbb {P}/\mathrm {p}}(x_{\mathbb {P}}) f_{i/\mathbb {P}}\left( \frac{x_i}{x_{\mathbb {P}}}, Q^2\right) \nonumber ,\\ P^{\mathrm {D}}_{A}&= \frac{f_{i/B}^{\mathrm {D}}(x_i, Q^2)}{f_{i/B}(x_i,Q^2)}\nonumber ,\\ P^{\mathrm {D}}_{B}&= \frac{f_{i/A}^{\mathrm {D}}(x_i, Q^2)}{f_{i/A}(x_i,Q^2)}, \end{aligned}$$

(10)

with $f_{i/\mathrm {p}}$ describing the PDF of the proton, $f_{i/\mathrm {p}}^{\mathrm {D}}$ being the diffractive part of the proton PDF defined as a convolution of the Pomeron flux in a proton ($f_{\mathbb {P}/\mathrm {p}}$) and the Pomeron PDFs ($f_{i/\mathbb {P}}$). The probabilities for side A, B to be the diffractive system are given as $P^{\mathrm {D}}_{A,B}$ and each relies on the variables of the opposite side.

This tentative probability is then used to classify an event as preliminary diffractive or non-diffractive. If non-diffractive, the events are handled as usual non-diffractive ones. If diffractive, the interleaved evolution of ISR, FSR and MPIs is applied, but only events surviving without additional MPIs are considered as fully diffractive events. The reasoning behind this is that additional MPIs in the $\mathrm {p}\mathrm {p}$ system would destroy the rapidity gap between the diffractive system and the elastically scattered proton. The gap survives if no further MPIs occur, and the event can be experimentally quantified as being diffractive, with e.g. the large rapidity gap method. This no-MPI requirement suppresses the probability for diffraction with respect to the tentative dPDF-based probability, and can thus be seen as a gap-survival factor. Unlike other methods of gap survival (e.g. [9, 49,50,51]) this method is performed on an event-by-event basis, thus inherently is a dynamical effect. Furthermore, it does not include any new parameters, but relies solely on the existing and well tested (for $\mathrm {p}\mathrm {p}/\mathrm {p}\overline{\mathrm {p}}$) MPI framework. Once the system is classified as diffractive, the full interleaved evolution is performed in the $\mathbb {P}\mathrm {p}$ subsystem. Here the model does not restrict the number of MPIs, as these will not destroy the rapidity gap between the scattered proton and the Pomeron remnant.

4.1 Hard diffraction with photons

In this article we extend the hard diffraction model to collisions involving one or two (intermediate) photons. The extension is straightforward. Changing the proton PDF in Eq. (10) to a photon PDF on one side, it is possible to describe hard diffraction in $\gamma \mathrm {p}$ interactions. Changing on both sides, the model is extended to $\gamma \gamma$ collisions. Thus Eq. (10) is valid in events with (intermediate) photons with the change $\mathrm {p}\rightarrow \gamma$. Connecting the event generation with an appropriate photon flux allows to study hard diffraction in both $\mathrm {e}\mathrm {p}$ and $\mathrm {e}^+\mathrm {e}^-$ collisions as well as in ultra-peripheral collisions of protons and nuclei. The differential cross section of the hard scattering ($X_{\mathrm {h}}$) in a diffractive system X, e.g. the dijet system within the diffractive system, for direct (dir) and resolved (res) photoproduction can then schematically be written as,

$$\begin{aligned} \mathrm {d}\sigma ^{AB\rightarrow X_\mathrm {h}B}_{\mathrm {dir}}&=f_{\gamma /A}(x)\otimes f_{\mathbb {P}/p}(x_{\mathbb {P}},t)\otimes f_{j/\mathbb {P}}(x_j,Q^2) \nonumber \\&\quad \otimes \, \mathrm {d}\sigma ^{\gamma j\rightarrow X_{\mathrm {h}}} , \nonumber \\ \mathrm {d}\sigma ^{AB\rightarrow X_\mathrm {h}B}_{\mathrm {res}}&=f_{\gamma /A}(x)\otimes f_{i/\gamma }(x_{\gamma }, Q^2)\otimes f_{\mathbb {P}/B}(x_{\mathbb {P}},t)\nonumber \\&\quad \otimes \, f_{j/\mathbb {P}}(x_j,Q^2)\otimes \mathrm {d}\sigma ^{i j\rightarrow X_{\mathrm {h}}}, \end{aligned}$$

(11)

with beam A emitting a photon, beam B emitting a Pomeron, and $AB\rightarrow X_{\mathrm {h}}B$ denoting that the diffractive system is present on side A. Changing $A\rightarrow B$ in Eq. (11) thus results in a diffractive system on side B. In the above, $f_{\gamma /A}$ denotes the photon flux from beam A, $f_{i/\gamma }$ the photon PDF, while $f_{\mathbb {P}/B}$ and $f_{j/\mathbb {P}}$ are the Pomeron flux and PDF, respectively. $\mathrm {d}\sigma ^{\gamma (i) j\rightarrow X_{\mathrm {h}}}$ are the partonic cross sections calculated from the hard scattering MEs. The full diffractive system X also contains partons from MPIs and beam remnants that also have to be taken into account, thus Eq. (11) only represent the hard subprocess part of the diffractive system. Presently, neither the double diffractive process $AB\rightarrow X_{\mathrm {h}}^AX_{\mathrm {h}}^B$ nor the central diffractive process $AB\rightarrow AX_{\mathrm {h}}B$ are modelled, and the Pomeron can only be extracted from protons and resolved photons. As the model is based on dPDFs and the dynamical gap survival derived from the MPI framework inside Pythia 8, the extension does not require any further modelling or parameters.

The dynamical gap survival is present only in the cases where the photon fluctuates into a hadronic state. Hence the tentative probability, Eq. (10), equates the final probability for diffraction in direct photoproduction and in the DIS regime, where no MPIs occur. In resolved photoproduction, the dynamical gap survival suppresses the tentative probability for diffraction, offering an explanation for the discrepancies between next-to-leading order (NLO) predictions for dijets in photoproduction compared to measured quantities at HERA, see e.g. [11, 13, 16]. The observed factorization breaking is not as striking as in $\mathrm {p}\mathrm {p}$ collisions, but the factorization-based calculation still overshoots the latest H1 analysis by roughly a factor of two [16].

Fig. 3

The two diffractive systems available for resolved photoproduction: either the proton is elastically scattered and the photon side contains the diffractive system (a), or the vector meson is elastically scattered and the proton side contains the diffractive system (b)

It should be noted that this extension allows for diffraction on both sides, i.e. the Pomeron can be extracted from the hadronic photon and/or the proton, see Fig. 3. Typically, the experiments only considered diffractive events where the diffractive system consists of a photon and a Pomeron, with a rapidity gap on the proton side (and a surviving proton, whether observed or not). The option to generate diffractive events on only one of the sides exist in Pythia 8, such as to avoid needless event generation.

4.2 Recent improvements in dPDFs

Since the publication of the hard diffraction model for $\mathrm {p}\mathrm {p}/\mathrm {p}\overline{\mathrm {p}}$, several improvements have been made for dPDFs. Work has been put into the inclusion of NLO corrections to the splitting kernels describing the evolution of the partons inside the Pomeron. Other work includes more recent fits to combined HERA data, or includes additional data samples into experiment-specific fits, so as to constrain some of the distributions in the dPDFs. A subset of these new dPDFs have been added to Pythia 8 recently and are briefly introduced below.

Fig. 4

The GKG18 LO Fit A, B and ZEUS SJ fluxes on a linear (a) and logarithmic (b) scale in $x_{\mathbb {P}}$. Note that t has been integrated over its kinematical range, $f(x_{\mathbb {P}})=\int \mathrm {d}t f(x_{\mathbb {P}},t)$

Fig. 5

The GKG18 LO Fit A, B and ZEUS SJ dPDFs on a linear (a, c) and logarithmic (b, d) scale. The upper figures shows the light quark content, the lower the gluonic content

Specifically two new sets of dPDFs have been introduced, along with the Pomeron fluxes used in these fits. The first set, the GKG18 dPDFs by Goharipour et al. [52], consists of two LO and two NLO dPDFs fitted to two different combined HERA data sets available, using the xFitter tool [53] recently extended to dPDFs. In addition, we consider an analysis released by the ZEUS collaboration offering three NLO dPDFs fitted to a larger sample of data. One of these, denoted ZEUS SJ, includes also diffractive DIS dijets from [54] in order to have better constraints for the gluon dPDF [14]. Using PDFs derived at NLO is not perfectly consistent with the LO matrix elements available in Pythia 8, but since the ZEUS SJ dPDF analysis is the only of the considered dPDF analyses including dijet data,¹ it is interesting to compare the results to other dPDFs.

Both the GKG18 and the ZEUS SJ fits uses the following parametrization for the Pomeron flux,

$$\begin{aligned} f_{\mathbb {P}}(x_{\mathbb {P}}, t) = A_{\mathbb {P}}\frac{\exp (B_{\mathbb {P}} t)}{x_{\mathbb {P}}^{2\alpha _{\mathbb {P}}-1}}, \end{aligned}$$

(12)

with $\alpha _{\mathbb {P}}=\alpha _{\mathbb {P}}(0) + \alpha _{\mathbb {P}}'t$ and A, B being parameters to be included in the fits. The dPDFs are typically parametrized as

$$\begin{aligned} zf_{i}(z, Q_0^2) = A_iz^{B_i}(1-z)^{C_i}, \end{aligned}$$

(13)

again with $A_i,B_i,C_i$ being parameters to be determined in the fits. The dPDFs are then evolved using standard DGLAP evolution [56,57,58,59,60,61,62] to higher $Q^2$. Different schemes for the inclusion of heavy quarks were invoked in the two fits; see the original papers for details. In both dPDFs the light quarks ($\mathrm {u}$, $\mathrm {d}$, $\mathrm {s}$) have been assumed equal at the starting scale, while heavy quarks ($\mathrm {c}$, $\mathrm {b}$) are generated dynamically above their mass thresholds. We show the new Pomeron PDFs and fluxes in Figs. 4 and 5, along with the H1 Fit B LO PDF [63] used as a default in Pythia 8. The GKG18 dPDFs are available with Pythia 8.240, while the ZEUS SJ set is expected in a forthcoming release.

5 Diffractive dijets in the photoproduction range

The production of dijets in a diffractive system is particularly interesting, as it provides valuable information on the validity of factorization theorems widely used in particle physics. These factorization theorems are not expected to hold in the case of diffractive dijets arising from resolved photoproduction, as this process essentially is a hadron–hadron collision, where the hard scattering factorization fails.

Both H1 and ZEUS have measured the production of diffractive dijets in both the photoproduction and DIS range. We here limit ourselves to showing results from two analyses, the H1 2007 and ZEUS 2008 analyses on diffractive dijets [12, 13]. Other analyses have been presented, including several ones examining only the DIS regime, but as the analysis codes or even the data itself have not always been preserved, we limit ourselves to reconstructing only a subset of these analyses. We aim to validate and provide the analyses used in this article within the Rivet framework [64].

Table 1 Kinematical cuts used in the experimental analyses by H1 [12] and ZEUS [13]. An asterisk ($^*$) indicates that the observable is evaluated in the photon–proton rest frame. $x_{\mathbb {P}}, M_Y, t$ are found in the rest frame of the hadronic system X, while the remaining are found in the laboratory frame

Both experiments have data on $\mathrm {e}\mathrm {p}$ collisions at $\sqrt{s}=318$ GeV using 27.5 GeV electrons and 920 GeV protons, with the proton moving in the $+z$ direction. Both use the large rapidity gap method for selecting diffractive systems. The experimental cuts in the two analyses are shown in Table 1. In the H1 analysis we concentrate on the differential cross sections as a function of four variables: invariant mass of the photon–proton system (W), transverse energy of the leading jet ($E_{\perp }^{*\,\mathrm {jet}\,1}$) and momentum fractions $z_{\mathbb {P}}^{\mathrm {obs}}$ and $x_{\gamma }^{\mathrm {obs}}$, both constructed from the measured jets as

$$\begin{aligned} x_{\gamma }^{\mathrm {obs}}&=\frac{\sum _{i=1}^2 \left( E^{\mathrm {jet},i} - p^{\mathrm {jet},i}_{z}\right) }{2 y E_e}, \nonumber \\ z_{\mathbb {P}}^{\mathrm {obs}}&=\frac{\sum _{i=1}^2 \left( E^{\mathrm {jet},i} + p^{\mathrm {jet},i}_{z}\right) }{2 x_{\mathbb {P}} E_p}, \end{aligned}$$

(14)

where $E_e$ ($E_p$) is the energy of the beam electron (proton) and the summation includes the two leading jets, i.e. the two with highest $E_{\perp }$. The inelasticity y and Pomeron momentum fraction w.r.t. the proton $x_{\mathbb {P}}$ are determined from the hadronic final state. In the ZEUS analysis the momentum fractions $z_{\mathbb {P}}^{\mathrm {obs}}$ and $x_{\gamma }^{\mathrm {obs}}$ are defined in terms of transverse energy and pseudorapidity of the jets,

$$\begin{aligned} x_{\gamma }^{\mathrm {obs}}&=\frac{\sum _{i=1}^2 E_{\perp }^{\mathrm {jet},\,i}\exp \left( -\eta ^{\mathrm {jet},i}\right) }{2 y E_e},\nonumber \\ z_{\mathbb {P}}^{\mathrm {obs}}&=\frac{\sum _{i=1}^2 E_{\perp }^{\mathrm {jet},\,i}\exp \left( \eta ^{\mathrm {jet},i}\right) }{2 x_{\mathbb {P}} E_p}, \end{aligned}$$

(15)

equivalent to the definitions in Eq. (14), if the jets are massless. In a LO parton-level calculation these definitions would exactly correspond to the momentum fraction of partons inside a photon ($x_{\gamma }$) and Pomeron ($z_{\mathbb {P}}$). Due to the underlying event, parton-shower emissions and hadronization effects, however, the connection between the measured $z_{\mathbb {P}}^{\mathrm {obs}}$ and $x_{\gamma }^{\mathrm {obs}}$ and the actual $x_{\gamma }$ and $z_{\mathbb {P}}$ is slightly smeared, but still Eqs. (14) and (15) serve as decent hadron-level estimates for the quantities. In place of W the ZEUS analysis provides the differential cross section in terms of invariant mass of the photon-Pomeron system, $M_X$.

There are several theoretical uncertainties affecting the distributions of the diffractive events. Here we focus on the most important ones:

• Renormalization- and factorization-scale variations, estimating the uncertainties of the LO descriptions in Pythia 8.

• dPDF variations affecting especially the $z^{\mathrm {obs}}_{\mathbb {P}}$ distribution and indirectly the number of events through the cuts on the squared momentum transfer, t, the momentum fraction of the beam carried by the Pomeron, $x_{\mathbb {P}}$ and the mass of the scattered (and possibly excited) proton, $M_Y$.

• $p_{\perp 0}^{\mathrm {ref}}$-variations, affecting the gap survival factor.

Fig. 6

The model with (solid lines) and without (dashed lines) gap survival compared to ZEUS data on $M_X$ (a), $z^{\mathrm {obs}}_{\mathbb {P}}$ (b), $x^{\mathrm {obs}}_{\gamma }$ (c) and $E_{\perp }^{\mathrm {jet\,1}}$ (d)

Fig. 7

The model with (solid lines) and without (dashed lines) gap survival compared to H1 data on W (a), $z^{\mathrm {obs}}_{\mathbb {P}}$ (b), $x^{\mathrm {obs}}_{\gamma }$ (c) and $E_{\perp }^{*\,\mathrm {jet\,1}}$ (d)

Other relevant parameters and distributions have also been varied, showing little or no effect on the end distributions. Remarkably, one of these was the choice of photon PDF. Pythia 8 uses the CJKL parametrization [42] as a default both in the hard process and in the shower and remnant description. As the MPI and ISR generation in the current photoproduction framework require some further approximations for the PDFs that are not universal and thus cannot be determined for an arbitrary PDF set, only the hard-process generation is affected by a change of photon PDF. Thus the effect of a different photon PDF on the various observables is not fully addressed with the present framework. The hard-process generation should, however, provide the leading photon PDF dependence. We find only a minimal change to the final distributions when changing to either the SaS [41], GRV [39] or GS-G [65] provided with LHAPDF5 [66]. There are two reasons for the weak dependence on photon PDFs. Firstly, the cuts applied by the experimental analyses presented here forces $x_{\gamma }$ to be rather large, where the photon PDFs are relatively well constrained by the LEP data. Secondly, the no-MPI requirement rejects mainly events from the low-$x_{\gamma }$ region, where the differences between the mentioned photon PDFs are more pronounced.

Two other analyses from HERA [11, 16] have also been used to check the current framework, giving results similar to the analyses presented here. For our baseline setup we show comparisons to both the H1 and ZEUS analyses, while for the more detailed variations we focus on comparisons to ZEUS.

5.1 Baseline results

In Figs. 6 and 7 we show the results obtained with Pythia 8 along with the experimental measurements. We show two simulated samples, one based on dPDFs solely without the dynamic gap survival (the “PDF” sample, dashed lines), and one including the dynamic gap survival (the “MPI” sample, solid lines). The results show that the “PDF” sample is too large compared to data in all distributions except for $x_{\gamma }$, thus showing evidence of factorization breaking. The “MPI” sample, however, seems to give a reasonably good description of data as the ratio of MC/data is smaller for the “MPI” sample than the “PDF” sample, thus hinting that it is the additional probability for multiparton interactions between the photon remnant and the proton that causes the factorization breaking.

A $\chi ^2$-test have been performed in order to quantify which of the models do better. Here, we have performed three different tests; using only either of the H1 or ZEUS datasets, or using both, Table 2. It is evident that the “MPI” model including the gap survival effect does a better job than the “PDF” model without it, within our baseline setup. The calculation of the $\chi ^2$ values include all differential cross sections provided by the experimental analyses, excluding the additional $x^{\mathrm {obs}}_{\gamma }$-binned distributions in ZEUS analysis to avoid counting the same data twice. Error correlations are not provided and so not considered.

Table 2 $\chi ^2$ tests using three different datasets

Fig. 8

The model with (solid lines) and without (dashed lines) gap survival compared to ZEUS data on $M_X$ in the direct-enhanced (a) and resolved-enhanced (b) regions

In general, most distributions are well described by the model including dynamical gap survival. The invariant mass distributions for the photon-Pomeron system ($M_X$) and for the photon–proton system (W) in Figs. 6 and 7a are both sensitive to the form of the photon flux from leptons. Both data sets are well compatible with the “MPI” samples, indicating that the standard Weizsäcker–Williams formula provide a good description of the flux.

It is, however, evident that in some observables the shape of the data is poorly described. Examples are $z^{\mathrm {obs}}_{\mathbb {P}}$ and $x^{\mathrm {obs}}_{\gamma }$, Figs. 6, 7b, c. The former is sensitive to the dPDFs used in the event generation. The baseline samples use the LO H1 Fit B flux and dPDF, fitted to data that is mainly sensitive to quarks. As the Pomeron is assumed to be primarily of gluonic content, it is expected that the vast majority of the dijets arise from gluon-induced processes. Thus a poorly-constrained gluon dPDF is expected to give discrepancies with distributions sensitive to this parameter, such as $z_{\mathbb {P}}$. In both the H1 and ZEUS analyses $z^{\mathrm {obs}}_{\mathbb {P}}$ is overestimated in the low end, while being underestimated in the high-$z^{\mathrm {obs}}_{\mathbb {P}}$ end. If the measured jets are dominantly gluon-induced, then it is expected that changing from the H1 LO Fit B dPDF to the ZEUS SJ fit should improve on the $z^{\mathrm {obs}}_{\mathbb {P}}$-distribution, as the low-$z^{\mathrm {obs}}_{\mathbb {P}}$ gluons are suppressed in this dPDF.

The latter observable, $x^{\mathrm {obs}}_{\gamma }$, is similarly underestimated in the low end and overestimated in the high end. The tight cut on $x_{\mathbb {P}}$ together with the requirement of high-$E_{\perp }$ jets reduces the contribution from lower values of $x^{\mathrm {obs}}_{\gamma }$. This suppresses the resolved contribution and therefore increases the relative contribution from direct processes, which typically are close to $x^{\mathrm {obs}}_{\gamma }=1$. The additional no-MPI requirement further suppresses the already low resolved contribution, and we end up with not being able to describe the shape of $x^{\mathrm {obs}}_{\gamma }$. As already discussed, the discrepancy cannot be explained with the uncertainties in the photon PDFs, as the sensitivity to different PDF analyses was found to be very low. The issue seems to be a problem with the relative normalizations of the direct and resolved contributions. This is evident from Fig. 8, where the ZEUS analysis conveniently splits the data into two regions, a direct- and a resolved-enhanced region with the division at $x^{\mathrm {obs}}_{\gamma }=0.75$. Here, the model underestimates the resolved-enriched part of the cross section and overestimates the direct-enriched part, confirming what we already observed in Figs. 6, and 7c.

Future measurements could shed more light on this issue, especially experimental setups in which the events passing the kinematical cuts would not be dominated by the direct contribution. In the experimental analyses considered here, a similar observation was made when comparing to a NLO calculation: the shape of $x^{\mathrm {obs}}_{\gamma }$ was well described by the NLO calculation (corresponding to our PDF selection) in the direct-enhanced region, but applying a constant suppression factor for the resolved contribution undershot the data at $x^{\mathrm {obs}}_{\gamma } < 0.75$, similar to what we observe. It is worth pointing out that both poorly-described distributions, $x^{\mathrm {obs}}_{\gamma }$ and $z^{\mathrm {obs}}_{\mathbb {P}}$, are constructed from the jet kinematics. Therefore further studies on jet reconstruction and their $\eta$ distributions could also offer some insights for the observed discrepancies.

Fig. 9

The model along with the uncertainty bands arising from varying the renormalization- and factorization scales compared to ZEUS data on $M_X$ (a), $z^{\mathrm {obs}}_{\mathbb {P}}$ (b), $x^{\mathrm {obs}}_{\gamma }$ (c) and $E_{\perp }^{\mathrm {jet\,1}}$ (d)

The jet variable $E_{\perp }$ can be used to check if the amount of activity within the diffractive system is properly described. As this system contains a Pomeron, it might very well be that the MPI parameters here could be different from the MPI parameters in the $\gamma \mathrm {p}$-system. It seems that using the same parameters for the $\gamma \mathbb {P}$ system as for $\gamma \mathrm {p}$ slightly overestimates the high-$E_{\perp }$ tail. This indicates that there might be too much MPI activity in the events, thus requiring a slightly larger $p_{\perp 0}^{\mathrm {ref}}$ value in the diffractive system than in the $\gamma \mathrm {p}$ system. The argument for a different $p_{\perp 0}^{\mathrm {ref}}$-value for $\gamma \mathrm {p}$ as compared to $\mathrm {p}\mathrm {p}$ can also be applied here: if the Pomeron has a smaller size than the proton, then the $p_{\perp 0}^{\mathrm {ref}}$-value can be increased. Having too much MPI activity in the $\gamma \mathbb {P}$-system may also push the $x^{\mathrm {obs}}_{\gamma }$ distribution towards higher values, as the $E_{\perp }$ of the jets may increase due to the underlying event. A full discussion of the MPI parameters in the diffractive system in $\mathrm {p}\mathrm {p}$ collisions has been provided in [8], but have not been pursued further here.

5.2 Scale variations

To probe the uncertainties in the choice of renormalization and factorization scales, $\mu _R$ and $\mu _F$, we employ the usual method of varying the scales up and down with a factor of two. Each is probed individually, such that one scale is kept fixed while the other is varied. Only the scales at matrix-element level are varied, thus the shower and MPI scales have been excluded from these variations. Each variation gives rise to an uncertainty band, and in Fig. 9 we show the envelope using the maximal value obtained from either of the two uncertainty bands. The envelope is dominated by the renormalization scale, giving the largest uncertainty in most of the figures shown – not unusual in a LO calculation. Note, however, that the scale uncertainty in the high-$x^{\mathrm {obs}}_{\gamma }$ bin actually reaches the upper error of the data point, essentially hinting that the model is able to describe the direct-enhanced region within theoretical uncertainties. The resolved region, however, cannot be fully accounted for within these theoretical uncertainties.

5.3 Variations of the dPDFs

As explained above, the considered observables are sensitive to the dPDFs, especially the fractional momentum carried by the parton from the Pomeron, $z^{\mathrm {obs}}_{\mathbb {P}}$. We here investigate if the increased amount of diffractive DIS data in the GKG LO dPDFs will provide a better description of the data than the less constrained H1 LO Fit B dPDF. We also show results obtained when using the NLO dPDF and flux from ZEUS SJ, as this dPDF includes data on diffractive dijets that is directly sensitive to the gluon distributions. Note, however, that a combination of NLO PDFs and LO matrix elements is still only accurate to LO and mixing different orders may result in different results compared to a situation where the matrix elements and PDF determination are consistently at the same perturbative order.

Fig. 10

The model without gap suppression using three different dPDFs: H1 LO Fit B (blue lines), GKG LO Fit A (green lines) and ZEUS NLO SJ (red lines) compared to ZEUS data on $M_X$ (a), $z^{\mathrm {obs}}_{\mathbb {P}}$ (b), $x^{\mathrm {obs}}_{\gamma }$ (c) and $E_{\perp }^{\mathrm {jet\,1}}$ (d)

Fig. 11

The model without gap suppression using the three dPDFs: H1 Fit B LO (blue lines), GKG LO Fit A (green lines) and ZEUS SJ (red lines) compared to ZEUS data on $M_X$ in the direct-enhanced (a) and resolved-enhanced (b) regions

In Fig. 10 we show results using two of the new dPDFs, ZEUS NLO SJ [14] and GKG LO Fit A [52] without the gap suppression factor. At first glance, the new dPDFs improve the overall description of data without a further need for suppression. Overall the new dPDFs seem to suppress the distributions as compared to H1 Fit B LO dPDF, with the ZEUS SJ dPDF performing slightly better than GKG LO Fit A as seen e.g. in the $z^{\mathrm {obs}}_{\mathbb {P}}$ distribution. Here, the ZEUS SJ dPDF flattens out at high $z^{\mathrm {obs}}_{\mathbb {P}}$ as compared to the GKG and H1 dPDFs, having a slightly larger $x_g$-distribution in this regime.

The distributions that the baseline study did not fully describe, also the new dPDFs fail to describe. Especially the $x^{\mathrm {obs}}_{\gamma }$ distribution is still underestimated at $x^{\mathrm {obs}}_{\gamma }<0.75$, which underlines the discrepancies with the relative normalization between the direct and resolved contributions. The $E_{\perp }^{\mathrm {jet\,1}}$ distribution is now well described with the GKG set. With the ZEUS SJ set the normalization is improved compared to the H1 Fit B but the shape of the distribution is similarly off.

A separation of $M_X$ into the two regimes, Fig. 11, shows that the direct-enhanced region is well described with the ZEUS SJ dPDFs. The GKG set improves the normalization but the shape of the distribution is still not compatible. The resolved region, however, is too suppressed with both of these, so the relative normalizations of the two contributions remain as an unresolved issue. Adding the gap suppression factor on top of this, Fig. 12, further suppresses the already suppressed resolved-enhanced region, worsening the agreement with the data in this regime. Little effect is seen in the direct-enhanced region, as expected.

Fig. 12

The model with gap suppression using the three dPDFs: H1 Fit B LO (blue lines), GKG LO Fit A (green lines) and ZEUS SJ (red lines) compared to ZEUS data on $M_X$ in the direct-enhanced (a) and resolved-enhanced (b) regions

These results thus puts forth the question whether the gap suppression is necessary if the dPDFs are refined and improved with additional diffractive data. The improvements seen especially with the ZEUS SJ dPDF in both the $x^{\mathrm {obs}}_{\gamma }$ and $z^{\mathrm {obs}}_{\mathbb {P}}$ distributions might hint towards this. As discussed earlier, this might partly follow from the tight cuts applied in the ZEUS analysis which does not leave much room for MPIs in the $\gamma \mathrm {p}$ system. Also, one should keep in mind that using NLO dPDFs with LO matrix elements might lead to different results compared to a full NLO calculation.

5.4 Variations of the screening parameter

The gap suppression method used here is highly sensitive to the model parameters of the MPI framework. Here we especially look at the screening parameter, $p_{\perp 0}^{\mathrm {ref}}$, as the value of this parameter differs between tunes to $\mathrm {e}\mathrm {p}$ and to $\mathrm {p}\mathrm {p}$ collisions. Changing the value of $p_{\perp 0}^{\mathrm {ref}}$ have only a small effect on the “PDF” samples. The “MPI” samples, however, are affected by the value of the screening parameter. A smaller value of $p_{\perp 0}^{\mathrm {ref}}$ results in more MPIs, thus we expect that the gap suppression will be larger if we decrease $p_{\perp 0}^{\mathrm {ref}}$ to its $\mathrm {p}\mathrm {p}$ value, as a smaller fraction of the events will survive the MPI-selection.

Fig. 13

The model with (solid lines) and without (dashed lines) gap suppression using two values of $p_{\perp 0}^{\mathrm {ref}}$: The $\mathrm {p}\mathrm {p}$-tune, $p_{\perp 0}^{\mathrm {ref}}=2.28~\mathrm {GeV}$ (red lines) and the $\mathrm {e}\mathrm {p}$-tune, $p_{\perp 0}^{\mathrm {ref}}=3.0~\mathrm {GeV}$ (blue lines). Again we show the samples in the observables $M_X$ (a), $z^{\mathrm {obs}}_{\mathbb {P}}$ (b), $x^{\mathrm {obs}}_{\gamma }$ (c) and $E_{\perp }^{\mathrm {jet\,1}}$ (d)

This effect is exactly what is seen in Fig. 13. The “PDF” samples are not affected, but the $\mathrm {p}\mathrm {p}$-tuned $p_{\perp 0}^{\mathrm {ref}}$ value in red causes a stronger suppression, best seen in the ratio plots where the solid red curves, the “MPI” sample with $p_{\perp 0}^{\mathrm {ref}}=2.28~\mathrm {GeV}$, is lower than the solid blue curves with $p_{\perp 0}^{\mathrm {ref}}=3.00~\mathrm {GeV}$. The value of $p_{\perp 0}^{\mathrm {ref}}$ has some effect on the shape of the distributions, mainly because a higher $M_X$ allows for more MPI activity, and thus a smaller fraction of events survive the no-MPI requirement. This means that the gap suppression increases with increasing energy available in the system, i.e. with increasing $M_X$, seen in Fig. 13a, where ratio-plot shows a suppression factor of approximately 0.9 in the low $M_X$ bin and 0.6 in the high $M_X$ bin.

5.5 Gap suppression factors

Several models have been proposed to explain the factorization breaking in diffractive hadronic collisions. Many of these employ an overall suppression factor, often relying primarily on the impact-parameter of the collision, see e.g. [49,50,51]. Some also include a suppression w.r.t. a kinematical variable, such as the $p_{\perp }$ of the diffractive dijets. But to our knowledge, the model of dynamical gap survival is the first of its kind to evaluate the gap survival on an event-by-event basis. This means it takes into account the kinematics of the entire event, and is thus also able to provide a gap suppression factor differential in any observable. In the model presented here, the ratio of “PDF” to “MPI” samples equates the gap survival factor, as the two samples only differ by the no-MPI requirement that determines the models definition of a fully diffractive event.

Fig. 14

The predicted gap suppression factors as a function of $M_X$ (a) and $E_{\perp }^{\mathrm {jet\,1}}$ (b) compared to the ZEUS analysis

Fig. 15

The predicted gap suppression factors as a function of W (a) and $E_{\perp }^{\mathrm {*\,jet\,1}}$ (b) compared to the H1 analysis

The theoretical uncertainties not directly related to MPI probability (e.g. scale variations) are expected to cancel in such a ratio. Even though many experimental analysis present similar ratios by using a NLO calculation as a baseline, such a ratio is not a measurable quantity, as it always require a theory-based estimation for the unsuppressed result. These ratios, however, are useful for demonstrating the effects arising from different models such as our dynamical rapidity gap survival. In order to estimate the factorisation-breaking effect in data w.r.t. our model, we show also the ratio between the data and the “PDF” sample.

In Fig. 14 we show the gap suppression differential in the observables $M_X$ and $E_{\perp }^{\mathrm {jet\,1}}$ from the ZEUS analysis and in Fig. 15 we show the gap suppression differential in the observables W and $E_{\perp }^{*\,\mathrm {jet\,1}}$ from the H1 analysis. These distributions demonstrate some of the main features of our dynamical rapidity gap survival model. We show the ratio of data to “PDF” sample (black dots) and the ratio of “MPI” to “PDF” sample (solid blue curve). This latter ratio is exactly the gap suppression factor predicted by the model. The shapes of the gap suppression factors agree reasonably well with the suppression factors derived from the data (the black dots), albeit the shape of Fig. 14b is off in the high-$E_{\perp }$ end, as already mentioned in the baseline results.

The model predicts a slowly decreasing suppression in $E_{\perp }^{(*)\,\mathrm {jet\,1}}$, while the suppression increases towards larger $M_X$ and W. This increase follows as the larger diffractive masses are correlated with larger invariant masses of the $\gamma \mathrm {p}$-system, where there is more room for MPIs at fixed jet $E_{\perp }$. This results in a larger fraction of the events having additional MPIs, thus a smaller fraction of the events survive as diffractive. Similarly, high-$E_{\perp }$ jets takes away more momentum than low-$E_{\perp }$ jets, again leaving less room for MPIs to take place. Thus we do not predict a flat overall suppression, as has often been applied in the experimental analyses.

Suppression factors in the range 0.7–0.9 are predicted in the shown observables. Given the uncertainty on the “PDF” sample, this is in agreement with the suppression factors of approximately 0.5–0.9, as observed by H1 [11, 12, 15, 16] and ZEUS [13]. A somewhat contradictory result was observed in Ref. [14], in which the ZEUS dijet data from Ref. [13] was found consistent with the purely factorization-based NLO calculation when using the ZEUS SJ dPDFs.

The experimental cuts applied in the ZEUS analysis, as compared to the analysis from H1, forces $x^{\mathrm {obs}}_{\gamma }$ to very large values, where the suppression from the MPIs does not have a large effect. Thus the ZEUS measurement requires less suppression than what is needed in the H1 measurement. The shown distributions, however, are still marred by the large theoretical uncertainties. One way to reduce these theoretical uncertainties would be to consider the ratio of photoproduced dijets to ones from DIS, as done e.g. in the recent H1 analysis [16]. The kinematic domain is slightly different due to different virtualities, but this would still greatly reduce dependency on dPDFs and scale variations, leaving only the mild photon PDF dependence in addition to the factorization breaking effects, that would be pronounced in this ratio. Unfortunately the current Pythia 8 description of DIS events at intermediate virtualities is not adequate to describe the inclusive DIS dijet data, so such a comparison is a project for the future.

6 Photoproduction in ultra-peripheral collisions

Because of the more than an order of magnitude larger $\sqrt{s}$ at the LHC, the accessible invariant masses of the $\gamma \mathrm {p}$ system are much larger than what could be studied at HERA. This allows us to study the factorization-breaking effects in hard-diffractive photoproduction in a previously unexplored kinematical region. Such measurements would fill the gap between the rather mild suppression observed at HERA and the striking effect observed in $\mathrm {p}\overline{\mathrm {p}}$ and $\mathrm {p}\mathrm {p}$ collisions at Tevatron and the LHC. This would provide important constraints for different models and thus valuable information about the underlying physics. Besides the results we present here, predictions for these processes have been computed in a framework based on a factorized NLO perturbative QCD calculation with two methods of gap survival probabilities, one with an overall suppression and one where the suppression is only present for resolved photons [67]. The authors here expect that the two scenarios can be distinguished at LHC, especially in the $x_{\gamma }^{\mathrm {obs}}$-distribution. The model presented in this work should thus be comparable to the latter suppression scheme from [67]. Another work considering similar processes is presented in Ref. [68].

In principle these measurements could be done in all kinds of hadronic and nuclear collisions, since all fast-moving charged particles generate a flux of photons. There are, however, some differences worth covering. In $\mathrm {p}\mathrm {p}$ collisions, the photons can be provided by either of the beam particles with an equal probability. The flux of photons is a bit softer for protons than with leptons, but still clearly harder than with nuclei. Experimentally it might be difficult to distinguish the photon-induced diffraction and “regular” double diffraction in $\mathrm {p}\mathrm {p}$, since both processes would leave a similar signature with rapidity gaps on both sides. In $\mathrm {p}{\mathrm {Pb}}$ collisions the heavy nucleus is the dominant source of photons, as the flux is amplified by the squared charge of the emitting nucleus, $Z^2$. Thus the photon-induced diffraction should overwhelm the QCD-originating colorless exchanges (Pomerons and Reggeons). Similarly, in ${\mathrm {Pb}}{\mathrm {Pb}}$ collisions the photon fluxes are large and thus would overwhelm the Regge exchanges. The latter type is currently not possible to model with Pythia 8, however, as in addition to regular MPIs, one should also take into account the further interactions between the resolved photon and the other nucleons, that could destroy the rapidity gap. Since these are currently not implemented in the photoproduction framework, we leave the ${\mathrm {Pb}}{\mathrm {Pb}}$ case for a future study.

Table 3 Kinematics for the UPC analyses

6.1 pPb collisions

The setup for the photoproduction in $\mathrm {p}{\mathrm {Pb}}$ collisions is the same as our default setup for $\mathrm {e}\mathrm {p}$, albeit the photon flux is now provided by Eq. (8). We here neglect the contribution where the proton would provide the photon flux, such that all photons arise from the nucleus. The jets are reconstructed with an anti-$k_{\mathrm {T}}$ algorithm using $R=1.0$ as implemented in FastJet package [69]. The applied cuts are presented in Table 3 and are very similar to the ones used by HERA analyses. The experimentally reachable lower cut on $E_{\perp }$ is not set in stone, however. This depends on how well the jets can be reconstructed in this process. On one hand, the underlying event activity is greatly reduced in UPCs as compared to $\mathrm {p}\mathrm {p}$ collisions, thus possibly allowing for a decrease of the reachable jet $E_{\perp }$. On the other hand, the increased W might require an increase of the minimum $E_{\perp }$ w.r.t. the HERA analyses. Feasibility of such a measurement has been recently demonstrated in a preliminary ATLAS study [70] which measured inclusive dijets in ultra-peripheral ${\mathrm {Pb}}{\mathrm {Pb}}$ collisions at the LHC.

The resulting differential cross sections for diffractive dijets from UPCs in $\mathrm {p}{\mathrm {Pb}}$ collisions are presented in Fig. 16. Similar to Sect. 4 we show the results differential in W, $M_X$, $x_{\gamma }^{\mathrm {obs}}$ and $z_{\mathbb {P}}^{\mathrm {obs}}$. The “PDF” samples (dashed lines) are without the gap suppression and “MPI” samples (solid lines) are with the gap suppression. The lower panels show the ratio of the two, corresponding to the rapidity gap suppression factor predicted by the model. As discussed earlier, the energy dependence of the $p_{\perp 0}$ screening parameter in $\gamma \mathrm {p}$ collisions was constrained by HERA data in a narrow W bin around $200~\mathrm {GeV}$. As the UPC events at the LHC will extend to much higher values of W, the poorly-constrained energy dependence of $p_{\perp 0}$ will generate some theoretical uncertainty for the predictions. To get a handle on this uncertainty we show samples with both the $\mathrm {p}\mathrm {p}$-tuned (red lines) and $\mathrm {e}\mathrm {p}$-tuned (blue lines) values for $p_{\perp 0}$.

The predicted gap suppression factor is rather flat as a function of $z_{\mathbb {P}}^{\mathrm {obs}}$ at around $\sim 0.7$. The suppression factor is, however, strongly dependent on W and $M_X$, also observed in the HERA comparisons. It is more pronounced at the LHC thanks to the extended range in W, with an average suppression being roughly two times larger than at HERA. A similar strong dependence is also seen in $x_{\gamma }^{\mathrm {obs}}$. As concluded earlier, the increasing suppression with W follows from the fact that the probability for MPIs is increased with a higher W, due to the increased cross sections for the QCD processes. Thus a larger number of tentatively diffractive events are rejected due to the additional MPIs. Similarly, decreasing $x_{\gamma }^{\mathrm {obs}}$ will leave more room for the MPIs to take place, since the momentum extracted from the photon to the primary jet production is decreased.

A reduction of the $p_{\perp 0}^{\mathrm {ref}}$-value from 3.00 to $2.28~\mathrm {GeV}$ increases the MPI probability, thus having a twofold effect. Firstly, it increases the jet cross section in the “PDF”-sample, as the additional MPIs allowed with the lower reference value increase the energy inside the jet cone. Secondly, the enhanced MPI probability rejects a larger number of tentatively diffractive events, thus giving a larger gap suppression effect. Collectively, these effects lead to 20–30% larger gap-suppression factors as compared to the $\gamma \mathrm {p}$ value for $p_{\perp 0}^{\mathrm {ref}}$.

Fig. 16

Cross section for diffractive dijets in ultra-peripheral $\mathrm {p}{\mathrm {Pb}}$ collisions for observables W (a), $M_X$ (b), $x^{\mathrm {obs}}_{\gamma }$ (c) and $z^{\mathrm {obs}}_{\mathbb {P}}$ (d). Vertical bars denote the statistical uncertainty in the MC generation

Fig. 17

Cross section for diffractive dijets in ultra-peripheral $\mathrm {p}\mathrm {p}$ collisions for observables W (a), $M_X$ (b), $x^{\mathrm {obs}}_{\gamma }$ (c) and $z^{\mathrm {obs}}_{\mathbb {P}}$ (d). Vertical bars denote the statistical uncertainty in the MC generation

6.2 pp collisions

The kinematical cuts applied in $\mathrm {p}\mathrm {p}$ equals those from $\mathrm {p}{\mathrm {Pb}}$. Due to the increased $\sqrt{s}$ and the harder photon spectrum from protons compared to heavy ions, the W range probed is extended to even larger values. When keeping jet kinematics fixed this leaves more room for MPIs in the $\gamma \mathrm {p}$-system, while also increasing the relative contribution from resolved photons. Thus the predicted gap-suppression factors are further increased here, as compared to $\mathrm {p}{\mathrm {Pb}}$ and $\mathrm {e}\mathrm {p}$ case, cf. Fig. 17. At extreme kinematics – high-$M_X$, low-$x_{\gamma }^{\mathrm {obs}}$ – the gap-suppression factors are almost as large as what have been found in hadronic diffractive $\mathrm {p}\mathrm {p}$ events. The $\mathrm {p}\mathrm {p}$ suppression factors should provide an estimate of the upper limit for photoproduction, as the latter includes the (unsuppressed) direct contribution. The suppression factors show a similar sensitivity to the value of $p_{\perp 0}^{\mathrm {ref}}$ as in $\mathrm {p}{\mathrm {Pb}}$ collisions, such that the lower value gives more suppression. Notice that the cross sections are calculated assuming that the photon is emitted from the beam with positive $p_z$.

A particularly interesting observable is the $x_{\gamma }^{\mathrm {obs}}$ distribution. Due to the extended W reach, the dijet production starts to be sensitive also to the low-x part of the photon PDFs. Here, the photon PDF analyses find that gluon distributions rise rapidly with decreasing x, the same tendency as seen in proton PDFs. This generates the observed rise of the cross section towards low values of $x_{\gamma }^{\mathrm {obs}}$ when the MPI rejection is not applied. However, the contribution from the low-$x_{\gamma }^{\mathrm {obs}}$ region is significantly reduced when the rejection is applied, as these events have a high probability for MPIs. Note, however, that there are large differences in the gluon distributions between different photon PDF analyses in this region. Thus here a variation of the photon PDF in the hard scattering could have some effect on the predicted gap-suppression factor, even though only very mild impact was seen in the HERA comparisons. But as most of these events with a soft gluon in the hard scattering will be rejected due to the presence of additional MPIs, the predicted cross-sections shown in Fig. 17 is expected to be rather stable against such variations. Further uncertainty again arises from the dPDFs. But as the purpose of the shown UPC results is to demonstrate the gap survival effects, we do not discuss the sensitivity to dPDF variation here explicitly.

7 Conclusions

In this paper we present a model for explaining the factorization-breaking effects seen in photoproduction events at HERA. The model, implemented in the general purpose event generator Pythia 8, is an extension of the hard diffraction model to photoproduction. It is a novel combination of several existing ideas, and it is the first model of its kind with a dynamical gap suppression based on the kinematics of the entire event.

The starting point is the Ingelman–Schlein approach, where the cross section is factorized into a Pomeron flux and a PDF, convoluted with the hard scattering cross section. The Pomeron flux and PDF are extracted from HERA data, but if used out-of-the-box these give an order-of-magnitude larger cross sections in pure hadron–hadron collisions, while the differences in photoproduction are around a factor of two at most. Thus, factorization was observed to be broken in diffractive events with a hard scale. The dynamical model extended here, explain this factorization breaking with additional MPI activity filling the rapidity gap used for experimental detection of the diffractive events. Thus the MPI framework of Pythia 8 is used as an additional suppression factor on an event-by-event basis, giving cross sections very similar to what is seen in data, both in $\mathrm {p}\mathrm {p}$ and $\mathrm {p}\overline{\mathrm {p}}$ events. As low virtuality photons are allowed to fluctuate into a hadronic state, MPIs are also possible in these systems. Thus the same mechanism is responsible for the factorization breaking in photoproduction events in $\mathrm {e}\mathrm {p}$ collisions, and also here the model predicts cross sections similar to what is seen in $\mathrm {e}\mathrm {p}$ data.

We present results obtained with the model compared to experimental data from H1 and ZEUS for diffractive dijet photoproduction. The agreement with the data is improved when the MPI rejection is applied, supporting the idea behind the factorization-breaking mechanism. However, the kinematical cuts applied by the experiments reduce the contribution from resolved photons, so the observed suppression is rather mild with the HERA kinematics, especially for the ZEUS data. The improvements in the dPDFs raises the question if such a suppression is actually needed, as the new dPDFs seem to describe data fairly well without, especially in the direct-enhanced region of phase space. Furthermore, there are several theoretical uncertainties that hamper the interpretation of the data, and the description is far from perfect for all considered distributions. Many of these theoretical uncertainties could be reduced by considering ratios of diffractive dijets in DIS and photoproduction regimes, but have not been pursued here as the description for DIS in Pythia 8 is not yet complete.

As an additional example for the range of the model, we present predictions for diffractive dijets in ultra-peripheral $\mathrm {p}\mathrm {p}$ and $\mathrm {p}{\mathrm {Pb}}$ collisions at the LHC. In these processes a quasi-real photon emitted from a proton or nucleus interacts with a proton from the other beam. Due to the larger invariant masses of the $\gamma \mathrm {p}$ system in these processes, the contribution from resolved photons is significantly increased. Thus UPCs is an excellent place to study the gap suppression in photoproduction. The results demonstrate that a measurement of photoproduced diffractive dijet cross sections in $\mathrm {p}\mathrm {p}$ collisions would provide very strong constraints on our dynamical rapidity gap survival model, as the effects are much more pronounced than with HERA kinematics. The distinct features of the model are well accessible within the kinematical limits for UPCs at LHC. If such a measurement is not feasible due to the pure QCD background, a measurement in $\mathrm {p}{\mathrm {Pb}}$ collisions would be sufficient to confirm the factorization breaking in diffractive photoproduction and provide constraints on the underlying mechanism.

Future work consists of opening up for different photon PDFs in the photoproduction framework, improving the DIS description in Pythia 8 and merging the two regimes in a consistent manner. The first allows for probing additional theoretical uncertainties of the photoproduction framework, the second allows for probing the double ratios of photoproduction to DIS cross sections for diffractive dijets. The merging of the two regimes would allow for full event generation of all photon virtualities needed for future collider studies. Similarly a combination of the current model and the Angantyr model for heavy ions is planned, such that $\mathrm {e}\mathrm {A}$ and ultraperipheral UPCs in $\mathrm {A}\mathrm {A}$ collisions could be probed as well.