1 Introduction

Photon-induced processes provide a rich testing ground for a wide range of physics effects. This is, on one hand, because photons will couple to any electromagnetically charged particle, resulting in a wide spectrum of accessible final states. On the other hand, photons have the quantum numbers of (neutral vector) mesons and their wave functions therefore have a sizeable hadronic component which lets them interact strongly, with correspondingly large cross sections. The production of low–multiplicity final states in \(\gamma \gamma \) collisions has been observed in many experiments [1,2,3,4], and yields interesting insights into the physics of hadrons and hadron resonances described by effective theories of the strong interactions. At increasing centre-of-mass energies of the colliding photons, new channels open up, and the production of jets has been studied, for example, at [5,6,7,8,9,10]. Similarly, the photo–production of various final states, including jets, has been analysed by the experiments [11,12,13,14,15,16,17].

The results obtained in these experiments have allowed to parameterise the parton content of the hadronic component of quasi-real photons in the form of parton distribution functions (PDFs), for example in the GRV [18], CJK [19, 20], SAL [21], and SaS [22, 23] sets.Footnote 1 In fact, based on these PDFs, satisfying agreement between data and calculations has been achieved, and a complete model for photon structure functions and high–energy photon interactions [24] has been encoded in [25]. There, the hard production of QCD final states at large scales, i.e. jets, is simulated in the usual way by dressing the hard parton–level matrix element with subsequent parton showers, the fragmentation of the resulting partons into hadrons during hadronization, possibly including an underlying event. Recently, the model was extended to also include the perpendicular component of the photon momentum [26].

Fig. 1
figure 1

Schematic sketch of the phase space mappings between the different steps in the initial states, i.e. the Equivalent Photon Approximation (EPA) and the Initial State Radiation (ISR), and the Matrix Element (ME). Each coordinates pair of Mandelstam-\(s^{\prime (\prime )}\) and rapidity \(y^{\prime (\prime )}\) is randomly sampled and the momenta are calculated as functions of these

Full size image

Pushing for higher accuracy, there have been a few predictions for inclusive jet-production at fixed-order at [27,28,29,30,31] and the [32], while more attention has recently been paid to exclusive meson production processes and photo-production at heavy-ion collisions [33,34,35,36,37,38].

Anticipating the increased precision requirements for successfully operating a possible future lepton collider such as or the planned electron–ion collider, , motivates to revisit the physics of photon–induced processes and to arrive at fully-differential predictions at Next-To-Leading Order (NLO) in QCD perturbation theory. We report here on the systematic inclusion of PDFs for quasi–real photons into the event generation framework [39, 40], the modelling of multiple-parton interactions, and the extenstion of the calculation to Next-to-Leading Order matched to the parton shower. The paper is organised as follows. In Sect. 2 we briefly discuss how combines the equivalent photon approximation (EPA) and photon PDFs, with some emphasis on the efficient integration over the resulting phase space, the matching to the parton shower at NLO and the multiple-parton interactions (MPI) model. In Sect. 3 we will show some first results of full Monte Carlo simulations, comparing results obtained with to data from both the and experiments at accuracy. We present comparisons to Leading Order and study the effect of the MPI. We summarise our findings in Sect. 4.

2 Equivalent photons and their PDFs

For the simulation of photo-production events in , we use its existing EPA interface, improved the phase space handling for the initial states for a more efficient integration, and added relevant photon PDFs to ’s internal PDF interface. The resulting code will be publically available as part of the upcoming release of 3.0; in the meantime it can be obtained from the authors upon request.

2.1 Phase space handling

In the following we detail the structures for efficient phase space sampling, using the most involved example of doubly resolved photon–photon collisions at lepton colliders, schematically depicted in Fig. 1; the cases of direct photons, the corresponding ISR terms do not need to be taken into account and the corresponding phase space integration is simplified.

The two incoming leptons have momenta \(p_1\) and \(p_2\), and a (beam) c.m.-system characterised by the c.m.-energy squared \(s_{12}\) and its rapidity \(y_{12}\) in the lab system. The momenta of the photons emitted by the leptons, \(p'_1\) and \(p'_2\), create a (photon) c.m.-system characterised by its c.m.-energy squared \(s'_{12}\) and rapidity \(y'_{12}\) with respect to the beam system. The partonic structure of the photons, as described by the PDFs, results in two partons with momenta \(p''_1\) and \(p''_2\) to finally enter the hard process which will result in final state particles with momenta \(q''_i\). The hard scattering is characterised by a c.m.-energy of \(s''_{12}\) and a rapidity \(y''_{12}\) w.r.t. the photon system. This structure requires two nested integrations for the two successive “initial states” (photons and partons): first an integration over \(s'_{12}\) and \(y'_{12}\), with factors given by the EPA spectra, and then an integration over \(s''_{12}\) and \(y''_{12}\), with factors given by the PDFs, before adding the integration over the final state phase space over the outgoing momenta \(q''_i\). Efficient integration over this complex phase space in is facilitated through the multi-channel method [41] with automatically generated integration channels that map out intrinsic structures such as s-channel resonances etc.

Table 1 Photon PDF libraries included in and their properties
Full size table

After the successful generation of a phase space point, the corresponding weight is calculated, given by the factors stemming from the EPA photon spectra and the PDF weights.

2.2 Equivalent photon approximation

The equivalent photon approximation encoded in the Weizsäcker-Williams formula [42,43,44] is based on the observation that quasi-virtual photons can be approximated through real photons for small virtualities \(Q^2 < \Lambda _\text {cut}^2\). As photo-production events are dominated by the interaction of low-virtuality photons, the differential cross-section can be substituted by \(\text {d}\sigma _{e X} = \sigma _{\gamma X}(Q^2=0) \text {d}n\). uses an improved version of the spectrum, following [45], which introduces the term proportional to \(m_e^2\) to the spectrum, and the photons are assumed to be collinear to the electron beam. The dependence of the photon spectrum on the photon virtuality is integrated out. This results in the following spectrum for electrons:

$$\begin{aligned} \textrm{d}n{} & {} = \frac{\alpha _\textrm{em}}{2 \pi } \frac{\textrm{d}x}{x} \left[ \left( 1 + {(1 - x)}^2 \right) \log \left( \frac{Q^2_\textrm{max}}{Q^2_\textrm{min}} \right) \right. \nonumber \\{} & {} \quad \left. - 2 m_e^2 x^2 \left( \frac{1}{Q^2_\textrm{min}} - \frac{1}{Q^2_\textrm{max}} \right) \right] \end{aligned}$$
(1)

Here, x denotes the ratio of photon to electron energy, \(\frac{E_\gamma }{E_e}\), and \(\alpha _\textrm{em}\) is the electromagnetic coupling constant. \(Q^2_\textrm{max}\) and \(Q^2_\textrm{min}\) denote the maximal and minimal photon virtuality, where the latter can be calculated from kinematic restrictions and is given by

$$\begin{aligned} Q^2_\textrm{min} = \frac{m_e^2 x^2}{1 - x} \,. \end{aligned}$$
(2)

The maximal virtuality is given by the experimental setup and the maximal deflection angle of the electron, \(\theta _\textrm{max}\), below which the hard process is still considered to be photon-induced. It is given by

$$\begin{aligned} Q^2_\textrm{max} = \textrm{min} \left( Q^2_\textrm{min} + E_e^2 (1 - x) \theta _\textrm{max}^2, Q^2_\textrm{max,fixed} \right) . \end{aligned}$$
(3)

Default choices are \(\theta _\textrm{max} = 0.3\) and \(Q^2_\textrm{max,fixed} = 3\ \textrm{GeV}^2\), but they can be overwritten by the user, cf. the manual [46].

2.3 Photon PDFs

To facilitate a comparison over different parameterisations, four PDF libraries have been included in , see Table 1 for a summary.

Currently, all PDFs are evaluated at virtuality \(Q^2 = 0\); the extension to virtual photons, taking also into account longitudinal polarisations, will be introduced in a later release. The extraction of a parton from the photon is complemented by the corresponding treatment of the remaining partons. For this a similar procedure is applied as in the remnant construction for a hadron, however, with a few simplifications. As there are no valence quarks present in the photon, the remnant is constructed as the anti-particle of the hard-scattering quark. In the case of the gluons, a quark-antiquark pair is constructed from one of the light quark generations. In the following processing, the photon remnants is treated analogously to a hadron remnant, i.e. its longitudinal momentum is distributed among the partons according to the kinematics, and their transverse momenta, given by some non-perturbative intrinsic \(k_T\) distribution, compensate each other.

2.4 NLO and matching to the parton shower

The procedure for a \(\textrm{NLO}\) calculation in photo-production has been presented in [47]. We will follow the same line of arguments for the matching to the parton shower, summarising the procedure in the following paragraphs. In contrast to the more familar case of hadronic PDFs, the evolution in the photon PDFs has an additional ”source” term, related to the photon splitting into a quark–anti-quark pair, which is proportional to \(\alpha _\textrm{em}\). This means that the evolution of the parton distributions \(f_i^{(\gamma )}\) for parton i in the photon follow [47]

$$\begin{aligned} \frac{\partial f_i^{(\gamma )}}{\partial \log \mu ^2} = \frac{\alpha _\textrm{em}}{2 \pi } P_{i\gamma } + \frac{\alpha _S}{2 \pi } \sum _j P_{ij} \otimes f_j^{(\gamma )}. \end{aligned}$$
(4)

The first term on the r.h.s. of the equation above gives rise to the ”anomalous component” mentioned before.

Fig. 2
figure 2

Differential dijet inclusive cross section with respect to \(\cos \theta ^*\), defined in Eq. (9), comparing results of our simulation at hadron-level with the LO simulation (top row) and with simulations at parton-level (bottom row) and with data from at an \(e^-e^+\) c.m.-energy of 198 GeV [5]. In the left and right panels the requirement \(x^\pm _\gamma <0.75\) and \(x_\gamma ^\pm >0.75\) are applied and enhance resolved and direct contributions, respectively

Full size image

At Next-to-Leading-Order, there will be additional collinear divergences in the real-corrections, stemming from the photon splitting. These divergences appear only in the direct-photon component, but they can be subtracted using the splitting kernel from the corresponding term in the photon PDF evolution.

To ensure an exact cancellation between these terms, the PDF has to use the same factorisation scheme as the subtraction scheme, which in the case of is the \(\overline{\textrm{MS}}\) scheme. This reduces the number of possible PDFs sets directly available for the NLO calculation in down to SAS1M and SAS2M from the SaS library. We note that additional PDFs can be made available easily by adding the respective factorization scheme dependent correction terms [48]. Apart from this subtlety in the choice of PDFs both NLO calculations and full simulations proceed in full analogy to the more familiar case of, e.g., proton–proton collisions.

3 Validation of the framework

We will now turn to the validation of our implementation, comparing results at precision with photo-production data from the and experiments. For each collider set-up and energy we generated samples of \(5 \cdot 10^6\) events per component at accuracy. For the calculation of matrix elements we used [49] and [50] for the tree-level matrix elements and subtraction terms [51] and [52] for the one-loop matrix elements. We added the parton shower [53] for the jet evolution, matched with the method [54, 55] as implemented in [56]. Underlying event effects have been included through an implementation of the Sjostrand-van Zijl model [57, 58] within ,Footnote 2 and the partons were hadronized with the cluster fragmentation of [59]. We consistently used the current default value for \(\alpha _S = 0.118\) with three-loop running. As described in 2.4, the event generation must currently use PDFs based on the \(\overline{\textrm{MS}}\) scheme, so the resolved-photon predictions were generated with, and averaged over, the SAS1M and SAS2M PDF sets. The factorisation scale and the renormalisation scale were both set to \(\mu _F = \mu _R = H_T/2\) and we kept the 7-point variation for the scale uncertainty estimate.

Fig. 3
figure 3

Distributions \(x_\gamma ^\pm \), collectively denoted as \(x_\gamma \), in different bins of average transverse jet energy: \({\bar{E}}_T\in [5\,\textrm{GeV}, 7\,\textrm{GeV}]\) (left), \({\bar{E}}_T\in [7\,\textrm{GeV}, 11\,\textrm{GeV}]\) (middle), \({\bar{E}}_T\in [11\,\textrm{GeV}, 25\,\textrm{GeV}]\) (right). Results of the simulation with accuracy are compared with results at LO and with data from at an \(e^-e^+\) c.m.-energy of 198 GeV [5]

Full size image

We used the [60] framework with the existing analyses implementing [5,6,7, 11, 12]. For each experiment, different components of the cross-section have to be summed over, for example,

$$\begin{aligned} \sigma _\text {tot} = \sigma _{\gamma \gamma } + \sigma _{j \gamma } + \sigma _{\gamma j} + \sigma _{jj} \end{aligned}$$
(5)

for and

$$\begin{aligned} \sigma _\text {tot} = \sigma _{\gamma j} + \sigma _{jj} \end{aligned}$$
(6)

for , where in both cases j denotes a photon or proton resolved through a PDF.

3.1 Comparison with data

The analysis of photon-induced di-jet production at 198 GeV c.m.-energy [5] offers the most differential observables, and we use it as the primary reference for our validation of photo-production at . To comply with experimental cuts, the \(k_T\) algorithm with a jet radius of \(R = 1.0\) is used and cuts of \(E_T > 4.5\ (2.5)\) GeV for the leading (subleading) jet are imposed.

$$\begin{aligned} x_\gamma ^\pm = \displaystyle { \frac{\sum \limits _{j=1,2} E^{(j)}\pm p_z^{(j)}}{\sum \limits _{i\in \textrm{hfs}} E^{(i)}\pm p_z^{(i)}}}, \end{aligned}$$
(7)

are used in this analysis to disentangle direct, singly, and doubly resolved production. In their definition, the sum in the numerator is over the two jets, and the sum in the denominator is over all hadronic final state particles, thereby distilling the energy-fraction the jets have w.r.t. the overall photon energies.

This is exemplified by, e.g., the distribution of events in

$$\begin{aligned} \cos \Theta ^* = \tanh \frac{\eta _{1}-\eta _2}{2}, \end{aligned}$$
(8)

an approximation of the angle between the two jets, and exhibited in Fig. 2. Apart from satisfying agreement with data, a few things are worth noticing here, which we will continue to observe also in the following: For \(x_\gamma ^\pm > 0.75\) the direct component dominates by about 1.5 orders of magnitude, with only small scale uncertainties, indicated by the pink band. Conversely, for \(x_\gamma ^\pm < 0.75\) the doubly-resolved component dominates with a significantly larger scale uncertainty, which, in this case, also includes factorization scale uncertainties. Intuitively, the singly-resolved component in both case ranges between the two other components. In addition we observe that hadronization effects reduce the cross section in the unresolved domanin, while the combination of hadronization and multiple parton scattering increases it in the doubly-resolved regime. The visible effect in the latter suggests that a careful retuning of the MPIs may further improve agreement with data.

We report that distributions in \(x_\gamma ^\pm \) for three different average di-jet transverse energies \({\bar{E}}_T = \sum _{j=1,2} E_T^{(j)} / 2\) experience a significant improvement in shape when going from Leading to Next-to-Leading Order, cf. Fig. 3. However, in the transition region between doubly resolved to unresolved events, we notice a clear difference in shape: While for \(x_\gamma ^\pm <0.6-0.7\) the prediction is relatively flat below the data, the underprediction at around \(x_\gamma ^\pm \approx 0.8\) persists at NLO. Apart from possibly insufficient photon PDFs – a point we will elucidate below – there are a number of possible explanations: First of all, as before, a retuning of MPIs may come to the rescue and fill up the gap.

Fig. 4
figure 4

Distributions of \(|\Delta \eta |\) (left), \(|\eta _\textrm{cntr}|\) (middle), and \(|\eta _\textrm{fwd}|\) (right), comparing and LO. Results of the simulation are compared with results from at an \(e^-e^+\) c.m.-energy of 198 GeV [5]

Full size image
Fig. 5
figure 5

Differential dijet inclusive cross section with respect to \(\cos \theta ^*\) for \(x^\textrm{obs}_\gamma <0.75\) (left) and \(x_\gamma ^\textrm{obs}>0.75\) (right), comparing results of our simulation with Run-1 data [11]

Full size image

Secondly, this drop around \(x_\gamma ^\pm \approx 0.7\) could be attributed to the missing QED splitting kernel in the evolution of the parton shower. Including this term would impact the backwards evolution of the photonic initial state radiation leading to a photon being reconstructed as the initial state also in the case of a resolved process. This again would lead to fewer radiation being generated, therefore shifting the distribution of the resolved process towards larger \(x_\gamma ^\pm \) values. The inclusion of this term in the evolution of the initial state showering is left for future work. Finally, we should stress that our singly resolved events are described by the \(2\rightarrow 2\) scattering of on-shell photons with partons from the resolved photons, an approximation which is probably not entirely correct as virtual photons would lead to a DIS-like scattering of the resolved photon, thereby inducing a somewhat different kinematics and scale choices.

Fig. 6
figure 6

Differential single-jet inclusive photo-production cross section with respect to the pseudo-rapidity of the leading jet, \(\frac{\textrm{d} \sigma }{\textrm{d} \eta }\), comparing results of our simulation at parton (left) and hadron-level (right) with Run-2 data [12]

Full size image
Fig. 7
figure 7

Single-jet inclusive transverse energy spectra for the leading jet in different pseudo-rapidity bins of the \(k_T\)-jets: \(0<\eta <1\) (top left), \(1<\eta <1.5\) (top right), \(1.5<\eta <2\) (bottom left), and \(2<\eta <2.5\) (bottom right), comparing results with Run 2 data [12]

Full size image

Figure 4 shows distributions of jet pseudo-rapidities and their differences. Again, the overall shape of the prediction is improved and the lowered NLO cross-section is countered by the inclusion of Multiple-Parton Interactions (MPIs).

3.2 Comparison with data

For the further validation of our implementation in electron–proton collisions we mainly rely on data [12] taken at Run 2. The kinematic cuts on the final states in the hard matrix element calculation were chosen to be a minimal transverse momentum of \(p_T > 11 (8)\) GeV for the (sub-)leading jets using the \(k_T\) clustering algorithm with radius \(R = 1.0\) to safely capture the phase space cuts of \(p_T^{(1)}>14\) GeV and \(p_T^{(2)}>11\) GeV used in the analysis of the Run-1 data [11], even after taking into account MPIs. For the Run-2 data [12], we chose \(p_T>13\) GeV to comply with the experimental cut of \(E_T>17\) GeV on the leading jet, with otherwise same settings.

Fig. 8
figure 8

Differential dijet inclusive cross section with respect to the transverse momentum of the leading jet, with \(0<\eta ^{(1)}<1\) (upper panel) and \(1<\eta ^{(1)}<2.4\) (lower panel) and in different bins for \(\eta ^{(2)}\), comparing results of our simulation at hadron-level incl. MPI effects with the LO simulation and with data [11] taken by at Run 1

Full size image

We use the same scale setting algorithm as before, and we also evaluate the theory uncertainties as combination of scale and PDF uncertainties like in the case of electron–positron colliders in the section above. Defining, in analogy to the case at lepton colliders above, Eqs. (7) and (8), respectively,

$$\begin{aligned} x_\gamma ^\textrm{obs}{} & {} = \frac{E_T^{(1)} e^{-\eta ^{(1)}}+ E_T^{(2)} e^{-\eta ^{(2)}}}{2yE_e} \;\;\;\text{ and }\;\;\; \nonumber \\{} & {} \cos \theta ^* = \tanh \frac{\eta ^{(1)}-\eta ^{(2)}}{2} . \end{aligned}$$
(9)

Here \(^{(1,2)}\) labels the leading and sub-leading jet, \(E_e\) is the lepton energy, and y is the energy fraction of the photon w.r.t. the lepton. As before we observe that the \(x_\gamma \) is excellently suited to disentangle unresolved and resolved photon interactions, cf. Fig. 5. In the resolved regime we again observe satisfactory agreement with data, indicating that the summation over the different components is correct. Interestingly, as suggested by the right panel of the figure, it appears as if the direct component is not sufficient to fully recover the experimental cross section. Possible explanations, as before, are related to the missing ”anomalous” \(\gamma \rightarrow q\bar{q}\) splitting in the backwards evolution, or, possibly more relevant here, a failure of the strictly on-shell approximation of the incident photons inherent to the treatment through EPA.

With this caveat in mind we will now turn into a more differential analysis of QCD final states in photo-production at . In Fig. 6 we compare the parton- against the hadron-level results and as previously observed, the shape improves significantly through the combined effect of hadronization and MPIs. This is most visible in the phase space of \(\eta < 0\), rendering the distribution of the simulation data flat compared to the experimental data. While there remains a discrepancy between simulation and data of around 10-20%, it might be explained by three observations. First, this analysis cuts requires only one jet in the final state which allows for contributions from the DIS region to leak into this measurement, for example if the scattered beam electron has not been correctly detected or identified. Secondly, the precise value of the strong coupling in the fit of the photon PDFs is not explicitly mentioned in the corresponding publications [22, 23]. An updated photon PDF fit would be performed with the current world average of \(\alpha _S\) and might further reduce the discrepancy. Lastly, the modelling of multiple-parton interactions for photon–proton interactions needs a fitting to the data for both proton–proton and photon–proton data. As neither have been tuned to data yet it can be suspected that the MPI will receive larger contributions, thereby improving overall agreement of simulation and data.

The spectrum of the jet transverse energies \(E_T\), displayed in Fig. 7 exhibits a slight shape in the distribution for forward jets and small \(E_T\). The same effect can be seen in Fig. 8, where the simulation describes data well for central leading jets, but worsens as both go more forward. As this part of the phase receives contributions from the MPIs, we again suspect that the naive parameter choice underestimate the amount of additionally generated radiation.

Fig. 9
figure 9

Distributions for \(\cos \Theta ^*\) at (left) and for \(\cos \theta ^*\) at (right), comparing ’s LO simulations with data from  [5] and  [11]. Scale uncertainties at LO are indicated by the pink band, while PDF uncertainties are shown with the blue hatched area

Full size image
Table 2 Inclusive cross-sections for one million events in a Run 1 dijet photo-production setup with two different proton PDFs and the same PDF for the photon
Full size table

3.3 Photon PDF for precision phenomenology

The predictions of photo-production cross sections and distributions in low-\(x_\gamma \) space exhibit large variations depending on the used photon PDF. In fact, these deviations can be as large or even larger in value than the estimate of higher-order corrections through the scale variations, especially for the simulation of photo-production at lepton colliders, see Fig. 9. There we show the very inclusive \(\cos \Theta ^*\) and \(\cos \theta ^*\) distributions, obtained at leading order, at and , respectively., indicating LO scale and PDF uncertainties separately. For the latter we present the full range of results from all available leading order PDFs, and we observe that the PDF uncertainties alone are of equal size as the LO scale uncertainties. This underlines the need for a comprehensive retuning of the photon PDFs with available data, and including higher-order calculations and simulations.

It is also worth noting that a consistent simulation of photo-production at hadron–lepton colliders necessitates a combined fit of the photon and proton PDFs. Depending on the proton PDF and its value for \(\alpha _S\) the inclusive cross-section in a Run 2 simulation gives deviations of about 20%, as can be seen in Tab. 2. This underlines the necessity for a systematic refit of photon PDFs to use in conjunction with modern proton PDFs. While no new data has been taken since the retiring of the collider, a consistent fit to all the data with the updated values for \(\alpha _S\) and including error estimates should increase the confidence in precision phenomenology for photo-production for the planned and other future colliders.

4 Summary

In this paper we reported on developments of the event generator towards the first hadron-level simulation of photo-production at NLO accuracy, to be published in our upcoming 3.0 public release. We described how the initial state phase space is treated to allow for a flexible yet efficient integration of different spectra and distributions, based on the equivalent photon approximation and the inclusion of a broad range of photon PDFs. We validated our model by comparison to data from the and experiments and found satisfactory agreement between the data and the simulation, noting that the non-perturbative aspects of the simulation, in particular intrinsic \(k_T\) and multiple-parton interactions, still require a comprehensive tuning to data. Apart from possible improvement in this sector of the simulation, our analysis of uncertainties underlined the necessity for a refitting of the photon PDFs, especially in view of the anticipated precision of the planned electron-ion collider. We believe that this step – a first refitting of photon PDFs after 20 years – based on higher-order calculations, modern simulation tools at least at NLO accuracy, and recent proton PDFs will siginifcantly improve the quality of our theoretical preprations for this new collider experiment.