Diphoton production in vector-boson scattering at the LHC at next-to-leading order QCD

In this paper, we present results at next-to-leading order (NLO) QCD for photon pair production in association with two jets via vector boson scattering within the Standard Model (SM), and also in an effective field theory framework with anomalous gauge coupling effects via bosonic dimension-6 and 8 operators. We observe that, compared to other processes in the class of two electroweak (EW) vector boson production in association with two jets, more exclusive cuts are needed in order to suppress the SM QCD-induced background channel. As expected, the NLO QCD corrections reduce the scale uncertainties considerably. Using a well-motivated dynamical scale choice, we find moderate $K$-factors for the EW-induced process while the QCD-induced channel receives much larger corrections. Furthermore, we observe that applying a cut of $\Delta \phi_{j_2 \gamma_1}^{\text{cut}}<2.5$ for the second hardest jet and the hardest photon helps to increase the signal significance and reduces the impact of higher-order QCD corrections.


I. INTRODUCTION
Vector boson pair production in association with two jets, denoted as V V jj in the following, has been studied in detail in recent years, both from the theoretical and the experimental sides of the Large Hadron Collider (LHC) physics community, since it provides information on weak boson scattering and because it is sensitive to beyond Standard Model (SM) physics via anomalous gauge boson couplings.
From the theoretical point of view, γγjj events can be produced at leading order (LO) via electroweak (EW)induced channel of order O(α 4 ), QCD-induced channel of order O(α 2 s α 2 ), and the interference between them of order O(α s α 3 ). The EW-induced channel is considered here to be the signal since it is sensitive to the weak boson scattering and to EW quartic gauge couplings, Fig. 1. The QCD-induced channel is an irreducible background. Despite the apparent O(α 2 /α 2 s ) suppression factor of the EW mechanism, once appropriate kinematical cuts are applied, both mechanisms yield integrated cross sections of the same order, and they provide distinct differential distributions in selected observables. Therefore, the EW mechanism turns out to be an excellent channel to test the SM and search for hints of beyond SM physics.
Using the vector boson scattering (VBS) cuts defined in Section III, the relative contributions at LO of the EW, QCD and interference contributions are 46%, 53%, and 1%, respectively (see Table I). At next-to-leading order (NLO) QCD, the distinction between EW-and QCD-induced channels becomes more blurry as new interference terms of the order O(α 2 s α 3 ) occur, mixing the two mechanisms together. These corrections have been found to be completely negligible in Ref.
[1] for the W + W + jj process with a looser VBS-cut setup. It can therefore be safely assumed that all interference effects between the EW-and QCD-induced mechanisms can be neglected at NLO QCD, given the fact the NLO QCD scale uncertainties are at the level of 20% on the QCD-induced cross section with tight VBS cuts, see Table IV.
The EW-induced mechanism can be further classified into s-channel contributions, which can be considered as a triple EW boson production with a subsequent hadronic decay of a vector boson, i.e. γγV → γγjj, and the t/u-channel VBS. In the VBS approximation defined in Ref. [2], which we will use to calculate the VBS signal at NLO QCD, the s-channel and its interference with the t/u-channels will be neglected because these effects are very small compared to the scale uncertainties of the QCD background when VBS cuts are applied. To further reduce the impact of the s-channel contributions, we will remove a window of 15 GeV around the V → jets (M V = (M W + M Z )/2) resonances. Moreover, since the interference between the t and u channels, occurring for sub-processes with identical quark lines, is negligible, the VBS approximation includes the contributions from the t and u channels independently.
Measuring the VBS signal is now an active field of research using Run-2 data at the LHC. Recently, 13 TeV results from ATLAS and CMS measuring the EWinduced channels involving two massive vector bosons with leptonic decays have been published for the samesign W W jj [3,4], W Zjj [5, 6], and ZZjj [7], showing agreement with the SM predictions. In this context, we would like to mention that NLO QCD corrections to the massive EW channels have been known for a decade in the VBS approximation [8][9][10][11][12], and have also been re-cently more precisely calculated beyond the VBS approximation with full off-shell and interference effects taken into account in Ref.
The EW induced Zγjj [30][31][32] and W ± γjj [33] channels have also been measured by the ATLAS and CMS experiments. Again, no deviation with the SM predictions was found. The NLO QCD corrections to the EW processes were calculated in Refs. [34,35] in the VBS approximation. For the QCD channels, the NLO QCD corrections were computed in Refs. [36,37]. The s-channel NLO QCD predictions are also available in the VBFNLO package. They were first computed in Refs. [20,21] for the leptonic decay modes and in Ref. [15] the hadronic decays were included.
EW di-photon production in association with two jets, γγjj, is the only process for which the NLO QCD predictions have not been studied and measurements are not available. It is an important process providing additional information on the EW boson scatterings and to beyond standard model physics via anomalous gauge boson couplings. Results at NLO QCD for the QCD-induced process have already been calculated in Refs. [38][39][40] In this paper, we present results at NLO QCD for the VBS t/u-channel for the process in the VBS approximation. Representative diagrams are shown in Figure 1. In the upper left diagram, the sensitivity of the process to vector boson scattering and to weak quartic gauge boson couplings is manifest. Using an effective field theory (EFT) approach, we have also included at NLO QCD the effect of dimension-6 and 8 operators involving the EW gauge bosons and the Higgs field. While dimension-6 operators also contribute here they are better investigated in vector-boson-pair production with much higher statistics. Besides, we have also computed the SM NLO QCDinduced channel, which is used here as the background process. Both EW-and QCD-induced channels have been implemented in the VBFNLO package, which will be available in the new release of the program or upon request. In addition, for a set of selected processes such as the complete list of EW V V jj processes, it can be linked at NLO to the Herwig [41] event generator to study shower and hadronization effects.
The paper is organized as follows: Section II describes the method used to compute the cross sections. Section III presents phenomenological results at the integrated cross section level and for differential distributions. Furthermore, as an illustrative example, the differential distribution of the diphoton invariant mass including anomalous couplings effects from an dimension-8 operator is shown. Finally, conclusions are presented in Section IV.

II. CALCULATIONAL SETUP
In this section, we shortly describe the method used to calculate γγjj production, both for the EW VBS channel and the QCD mechanism. We closely follow the strategy of other similar EW and QCD processes implemented in VBFNLO. The codes are based on a simplification of the VBS and QCD + − γjj processes (called EW and QCD Z γjj for simplicity from now on).
We work in the five-flavor scheme and top-quark loops are taken into account in the QCD mechanism. We note that sub-processes with external bottom quarks are included as long as there is no external top quark involved. This means that for the EW-induced channel, the bottom quark is included in the neutral currents but not in the charged currents. For the propagator of the massive vector bosons, a fixed width is used for both resonant and non-resonant propagators. The weak-mixing angle and other coupling constants are kept at real values, i.e. complex masses are not used to calculate the weak-mixing angle. Photon fragmentation functions are not included and we use instead the Frixione smooth-cone isolation criteria [42] for the external photons.
For the QCD γγjj induced channel, the gauge invariant set of diagrams with the photons attached to closed-quark loops are discarded due to its per mille level contribution effect [27]. Closed-quarks loops with three or two gluons are included, however.
For the EW process, we work in the VBS approximation and consider only the t/u-channels Feynman diagrams. Interference effects between the t-and u-channel diagrams as well as with the s-channel contributions are neglected. A detailed study on the validity of this approximation for the same-sign W W jj production channel can be found in Ref. [13] -the accuracy of the VBS approximation once tight VBS cuts are applied should be enough for present and near future experiments, being around the few-percent level. Consistently, virtual corrections with a gluon exchange between the two quark lines are not considered since, due to the color structure of the amplitudes, they are only non-vanishing for the interference of the t-and u-channel diagrams, which are phase-space suppressed and neglected in the VBS approximation.
In the following, we give a brief description of the technical details implemented to generate the code for the EW process -the QCD induced channel has been obtained using a similar procedure. As mentioned above, our code is based on a simplification of the VBS + − γjj amplitudes already programmed in VBFNLO. We therefore first discuss the implementation of this process. We use the effective current approach and the spin-helicity formalism of Ref. [43,44], which allows to factorize the EW-dependent leptonic tensors from the QCD amplitudes. The lepton pair can either originate from the decay of an intermediate boson, V 1 = Z/γ * → + − or V = Z/γ * → + − γ, which is radiated off the quark lines, or it can stem from the scattering of the t-channel vector bosons, e.g. V V /W + W − → + − . We denote the latter as leptonic tensor T and we define the tensors T γ and T γ in a similar way. Using this notation, the + − γjj amplitude can be written as the sum of the generic qq → V 1 γqq, qq →V qq, qq → V 1 T γ qq, qq → T γqq, and qq → T γ qq contributions. All spin correlation and off-shell effects are taken into account in the definition of the leptonic tensors and decay currents.
With this approach, it is trivial to obtain the γγjj code. We simply select the same generic amplitudes then replace V 1 = γ 2 andV = 0 in the effective currents. In the amplitudes containing the leptonic tensors, we set T = T γ2 , T γ = 0. We use the sub-index 2 to distinguish the final state photons. Finally, we define new leptonic tensors V V /W + W − → γγ ≡ T γγ , and use it in the qq → T γγ qq amplitude.
We have performed several tests to validate our code. The LO and real radiation matrix elements have been cross-checked with Madgraph [45] at the amplitude level and with Sherpa [46,47] for integrated cross sections, finding agreement at the machine precision and per mille level, respectively. For the EW-induced process, we have subtracted the s-channel contributions from the real matrix elements in Sherpa, otherwise, percent-level agreement is found for typical VBS cuts. The impact of the neglected s-channel is only noticeable in the real correc-tions, thus, its effect at the total NLO QCD cross section should be below the percent level. Additionally, for the virtual contributions the factorization of the poles, gauge and parametrization invariance [44] have been proved at the machine precision level and the convergence of the Catani-Seymour subtraction algorithm at the per mille level or better. Finally, the numerical stability of the code is controlled via Ward identities. For the EW process the amplitude is set to zero, if the identities are not satisfied at the per mille level using double precision. For the QCD process, we have a rescue system in quadruple precision, and the amplitude is only set to zero for the points that do not satisfy the Ward identities at the per mille level in quadruple precision. The fraction of points for which this takes place is at the per mille level in the EW process. Since the contribution of the virtual corrections, after the cancellation of the infrared divergences, is about a few percent, the error induced by this procedure is irrelevant. For the QCD channel, the fraction of rejected points is well below the per mille level, thus, completely negligible.
Additionally, for the QCD-induced mechanism, using the setup described in Ref. [39], we obtain [σ LO = 2.045(1), σ NLO = 2.714(3)] pb, to be compared with [2.046(2), 2.691(7)] pb of Ref. [39], which includes the fermion loops with one or two photons coupling to it. For the sake of comparison, we have estimated this fermionloop contribution separately and obtained σ γ,γγ fer. loop ≈ −7.8 fb, giving an agreement at the 2 standard-deviation level for the NLO cross section. Furthermore, the topquark loop contribution in the gluon self-energies and three-gluon vertices is included in our calculation while being omitted in Ref. [39]. This small effect may contribute to the above small discrepancy.
Concerning the anomalous-gauge-coupling implementation, new physics effects only occur in the photonic tensors T γ , T γ2 , and T γγ . For dimension-6 and 8 operators, which have also been implemented for the other EW V V jj processes in the VBFNLO program [17,48,49], we have crosschecked our implementation at the LOamplitude level against Madgraph with the FeynRules [50,51] model file EWdim6 [49,52,53] (for dimension-6) and with the FeynRules model files for quartic-gauge couplings [54,55] (for dimension-8). Agreement at the machine-precision level has been found at random phasespace points.

III. PHENOMENOLOGICAL RESULTS
We use the following SM input parameters [56] from which the electromagnetic coupling is calculated /π and the widths as Γ Z = 2.507426 GeV, Γ W = 2.096211 GeV. The mass of all the other light fermions are neglected as the results are insensitive to them. The top-quark mass dependence occurs via the fermion-loop corrections in the QCD-induced channels. This contribution is known to be very small, hence the results depend very weakly on the top-quark mass. The Cabibbo-Kobayashi-Maskawa matrix is set to unity in our calculations.
Concerning kinematic cuts, we require where jets are reconstructed from massless partons satisfying |y parton | < 5 using the anti-k t algorithm [57] with the radius parameter R = 0.4, y denoting the rapidity. In this paper we define the R-separation between two particles a and b as ∆R ab = (∆y ab ) 2 + (∆φ ab ) 2 where ∆y ab = |y a − y b | (here and the following) and ∆φ ab = |φ a − φ b | ≤ π. In order to isolate prompt photon events and minimize the parton-to-photon fragmentation contribution, we use Frixione's smooth-cone isolation criteria [42]. Events are accepted if (4) where the index j runs over all partons, δ 0 = 0.4 is the cone-radius parameter, and is the efficiency. The notation R γj = ∆R γj has been used for shortness. We choose as default = 0.05, following the recommendation of tight-isolation cuts in Ref. [58]. Note that Eq. (4) must be applied independently for the two photons. We see clearly from Eq. (4) that soft gluons are accepted while a hard quark exactly collinear to a photon is rejected, thereby ensuring IR safety while removing the collinear contribution. The cut ∆R jγ > 0.8 in Eq. (3) helps to isolate the photons further from the jets. In our calculation, since the quark-photon collinear events have been discarded, no fragmentation contribution occurs. In experiments, the above smooth-cone isolation criteria cannot be exactly implemented due to the finite resolution of the detector. However, using a tight-isolation cut both in theoretical calculations and measurements (with either a standard-cone cut or a discretized version of the smoothcone criteria) is expected to produce a very good agreement at the few-percent level, according to the study in Ref. [58] for pp → γγ. This is because the tight cut suppresses the fragmentation contribution, where a photon usually lies inside a hadronic jet.
To enhance the signal, we employ further the following VBS cuts m j1j2 > 800 GeV, |y j1 − y j2 | > 3, y j1 y j2 < 0, where j 1 and j 2 are the two tagging jets ordered by p T with j 1 being the hardest jet. Furthermore, the photons are required to fall inside the rapidity gap of the two tagging jets, i.e. y min < y γi < y max with y min = min(y j1 , y j2 ), y max = max(y j1 , y j2 ). Additionally, to remove the pp → V γγ (V → jj) contribution with V = W, Z, we accept only events satisfying where m jets can be the mass of any massive jet or the invariant mass of any combination of two or more jets. Note that we have m 3jets > 800 GeV due to the invariant mass cut in Eq. (5) in the default setup. However, when the value of the m j1j2 cut becomes smaller than (M W + M Z )/2 + 15 GeV as can happen in the scan shown in Figs. 4 and 10, then the three-jet contribution is affected by the cut in Eq. (6). The ∆y j1j2 cut choice together with this last cut should guarantee the validity of the VBS approximation at the percent level, as demonstrated in EW-Hjjj production [59], independently of the di-jet invariant mass cut used. To calculate hadronic cross sections, parton distribution functions (PDF) and the strong coupling constant α s (µ R ) are calculated using the LHAPDF6 program [60] with the PDF4LHC15 nlo 100 set [61][62][63][64]. The same PDF set is used both for the LO and NLO results. Our default choice for the renormalization and factorization scales is for both EW and QCD-induced contributions (see the scale-dependence discussion in the next section). In the following, SM results for the LHC at √ s = 13 TeV will be first presented. After that we will show briefly results with anomalous gauge couplings from dimension-8 operators as an illustration of the capabilities of our computer program.  With the above setup, we first show in Table I the full LO cross section with full EW amplitudes (t, u, s channels included), QCD amplitudes and their interference. We see that the EW-QCD interference is 1.1%, completely negligible compared to the scale uncertainties of the QCD cross section discussed below. In addition, the VBS approximation is also provided, showing an excellent agreement with the full EW result.
From now on we will use the notation EW to denote the VBS approximation. Within the default VBS cuts, we expect that the difference between the EW and the full EW results is completely invisible also at NLO QCD, given the fact that the QCD corrections are very small. 40 We now discuss the scale dependence of the integrated cross sections. Setting µ F = µ R = ξµ 0 , the scale dependence is shown in Fig. 2 for the EW-induced process. The scale dependence of the QCD-induced process has already been provided in Refs. [38][39][40], showing that H T /2 is a good scale choice. For the sake of comparisons, results with a fixed scale choice centered around µ fix 0 = M Z and with other dynamical scale choices √ p T,j1 p T,j2 and Q i are also shown for the EW case. Q i is the momentum transfer from quark line 'i' and is set independently for both quark lines. At LO, there is no µ R , hence the dependence comes from µ F via the PDFs. The scale choice H T /2 gives the smallest correction around the central scale µ 0 when moving from LO to NLO and, in this view, the scale √ p T,j1 p T,j2 comes second, being more consistent than the Q i scale. Moreover, the dependence on ξ at NLO shows that H T /2 and √ p T,j1 p T,j2 provide the most stable behavior. From this evidence, we conclude that H T /2 is a good scale choice for calculating the EW cross sections in the present setup. The scale dependence for EW pp → Zγjj presented in Ref. [35], for a similar setup, also shows that H T /2 gives the most stable behavior at NLO. We next study the dependence of the cross sections on the photon-isolation parameters. Results are shown in Fig. 3 and Table II. In Fig. 3, we show the dependence of the EW and QCD cross sections on the parameter defined in Eq. (4) for two cases δ 0 = 0.4 and 0.7. Note that, for all cases the LO cross section is independent of δ 0 and of because of the ∆R jγ > 0.8 cut. At NLO, one additional partonic radiation occurs. This radiation is included or excluded depending on the values of δ 0 and . Only events with at least a parton in the vicinity of a photon satisfying ∆R γ,parton < δ 0 can be rejected. This   explains why the NLO cross section decreases as δ 0 increases. It also explains why the cross section is more sensitive to when the cone-radius δ 0 is larger. Numerically, we find that the EW cross section increases about 2% (6%) when varying ∈ (0.01, 1) for δ 0 = 0.4 (0.7). For the QCD channel, we have 35% (97%), correspondingly. We observe that the dependence on is significantly milder for the smaller cone radius. However, the K-factors defined as σ N LO /σ LO are larger, in particular for the QCD channel, when the cone radius is decreased.
In experimental analyses, it is important to find out an optimal set of cuts to enhance the EW-induced channel. For this purpose, we show in Fig. 4 and Table III the      pendence of the significance defined as S = EW/ √ QCD, where EW and QCD represent the number of events of the two production mechanisms, calculated at NLO with luminosity L = 1 fb −1 , using the cuts m j1j2 > m cut j1j2 and ∆y j1j2 > ∆y cut j1j2 . For the other cuts, default values are used. We see that the maximal significance region of S > 4.25 is m cut j1j2 ∈ (1.1, 1.6) TeV and ∆y cut j1j2 < 3.5. If we require S > 3.75 then the region becomes m cut j1j2 ∈ (0.64, 2.2) TeV and ∆y cut j1j2 < 4.8. We note that for arbitrary luminosities L, the significance can be calculated as S(L) = S(L = 1 fb −1 ) · √ L fb. It is also interesting to compare the γγjj process with similar ones of + − γjj and ν γjj with = e, µ. This comparison is shown in Table IV. The default kinematic cuts used for the γγjj process are also applied to the other ones. Additionally, to select the charged leptons we use p T, > 30 GeV, |y | < 2.5, ∆R j, > 0.4, ∆R γ > 0.8.
We further require, for the + − γjj process, to suppress the γ → + − and Z → + − γ contributions, respectively. The value of 120 GeV was recommended in Ref. [37]. For the ν γjj processes, following Refs. [36,65], we use with E T = m 2 γ + p 2 T, γ + p T,ν to suppress the W → ν γ contribution, which is, as the Z → + − γ one, a background to the anomalous-quartic-gauge-coupling measurements. Concerning the scale choice, similar to Eq. (7) for the γγjj case, we use where E T,V = m 2 V + p 2 T,V with m V being the reconstructed mass. With this setup, we have used VBFNLO version 3.0.0 beta 4, where the calculations of + − γjj EW [35] and QCD [37] processes, as well as ν γjj EW [34] and QCD [36] processes are included, to produce the NLO cross sections. We see that the significance for the γγjj process is largest.
We now turn to differential cross sections. To understand the energy scale of the final state particles, we show in Fig. 5 the transverse momentum distributions of the tagging jets and the photons, individually, for the EW-induced channel at LO and NLO and for the QCD-induced process at NLO. Given the above cuts, the EW cross section is largest at p T,j1 ≈ 110 GeV, p T,j2 ≈ 40 GeV, p T,γ1 ≈ 60 GeV, p T,γ2 ≈ 35 GeV at NLO. For the QCD background, the maximal position is at p T,j1 ≈ 70 GeV, p T,j2 ≈ 40 GeV, p T,γ1 ≈ 60 GeV, p T,γ2 ≈ 35 GeV. We observe that, for the photon distributions, the EW and QCD processes have the same shapes. However, for the jets, the QCD distribution falls faster than the EW one. On the small panels, the K-factors are shown for both EW and QCD processes. On all panels, the scale-uncertainty bands calculated from the maximum and minimum of [dσ(H T /4), dσ(H T /2), dσ(H T )]   are plotted, where both scales are set equal. The Kfactor bands are calculated from these maximum and minimum with a common normalization to the central LO cross section dσ LO (H T /2). As expected, we see that, for the EW process, the NLO bands are much shrunk compared to the LO ones. We also see that the scale uncertainties on the QCD process are significantly larger than on the EW one. The K-factors of the QCD channel are also much larger. In the low p T region, the K-factors reach very large values, e.g. at p T,j1 ≈ 40 GeV we get K QCD ≈ 3.8, K EW ≈ 1.3. In the maximal cross section region, where p T,j1 ≈ 100 GeV, those values change to 1.9 and 1.0, respectively. This result is possibly caused by the large separation of scales, e.g. p T,j and m jj , and it indicates that a fixed-order NLO calculation is not sufficient for reliable predictions when applying VBF cuts. Since the calculation of the next-to-next-to-leadingorder (NNLO) prediction of this process is beyond the current reach of higher order calculations, we propose the merging of NLO predictions with various jet multiplicities within a parton-shower framework to improve the accuracy of the prediction of the QCD channel, which  we leave for a future work. For the EW signal, an NLO prediction is sufficient.
In order to see how the events look like, we show in Fig. 6 the rapidity and azimuthal-angle separation between the two tagging jets and between the photons. The ∆y j1j2 plot is surprising. Normally, we expect that the most likely rapidity separation for the EW process is larger than for the QCD one, as shown in Ref. [26] for the case of W + W + jj production and in Ref. [35] for Zγjj. However, the plot in Fig. 6 (top left) shows the contrary: the QCD cross section peaks at ∆y j1j2 ≈ 5.5 while the EW one at about 5.3. This means that a large ∆y j1j2 cut is not efficient to enhance the signal-over-background ratio in this case. We have chosen the cut ∆y j1j2 > 3 as a default setting in this paper. However, the ∆y j1j2 distributions show that using a looser cut can be good as well. The difference between the QCD and EW channels is also very pronounced in the ∆φ j1j2 distributions. Although both processes peak at ∆φ j1j2 ≈ π as expected, because the hardest jet is recoiling against other particles, the QCD cross section is more uniformly distributed than the EW one. Similar distributions for the two photons are also presented. The ∆y γγ plot shows that a small rapidity separation is the preferred configuration   for both EW and QCD processes. The EW distribution has one local maximum at ∆y ≈ 0.45 then drops gently as the separation increases. The QCD distribution however has two local maximums at ∆y ≈ 0.1 and at 0.45, with the first peak a bit higher. After the second peak, the QCD cross section drops more rapidly than the EW one. The ∆φ γγ distributions have a jump at ∆φ ≈ 0.4 because of the ∆R γγ > 0.4 cut. The results show that a large ∆φ separation is slightly preferred for both channels. Compared to the two-jet case, it is more uniformly distributed.
The distribution of the invariant mass of the two tagging jets and of the two photons are shown in Fig. 7. As expected, we see that the QCD-induced cross section drops more rapidly with increasing di-jet invariant mass than the EW-induced one. This justifies the high value of the m j1j2 cut used in this paper. The EW K-factor is close to unity for a large range of the invariant mass, up to 3 TeV, suggesting that this distribution is perturbatively well behaved even at very high energies for the signal. For the QCD background, the K-factor is also rather constant, but much larger, slowly increasing from 1.7 to 2.0. This large K-factor together with a large uncertainty band again signal the importance of predictions beyond the fixed-order NLO accuracy for the QCD process. The m γγ plot on the right shows that the EW NLO distribution peaks at 120 GeV while the QCD at 90 GeV. We note that the Higgs contribution is not included in the EW channel as it is beyond the fixed-order corrections considered here. Concerning the K-factor, it is very close to unity for the EW process. For the QCD channel, it is large at small invariant masses, decreasing from 2.7 at m γγ ≈ 30 GeV, reaching 1.0 at around 600 GeV then being rather constant after that.
To understand the jet-photon separations, we show in Fig. 8 the distributions of the z γi with i = 1, 2 (top row) defined as z X = y X − (y j1 + y j2 )/2 |y j1 − y j2 | , X ∈ (γ 1 , γ 2 , j 3 ), (12) and of the R-separation between the hardest photon and the tagging jets (bottom row). The z γi distributions show the distance of the photon with respect to the tagging jets with values of 1/2 and −1/2 when the photon equals the rapidity of jet 1 and 2, respectively. Due to the cuts imposed, we observe as expected that in both EW and QCD induced channels the photons are nearly homogeneously distributed in the center between the tagging jets. This has to be compared with the differential distribution of z j3 , computed here only at LO and not shown, where the third jet aligns in the EW-induced channel with either of the two leading jets while in the QCD-induced process, the distributions have a pronounced peak at z j3 = 0. Note that, the two z γi plots in Fig. 8 are not identical because the photons are ordered by p T . As expected, the distributions are flatter for the softer photon. The peaks in the ∆R j1γ1 distributions show that the preferred configuration is when the hardest jet and the hardest photon are back-to-back (i.e. ∆φ j1γ1 = π ) for both channels. For small separations, the K-factors are large with K ≈ 2.2 (4.8) for EW (QCD) induced mechanisms at ∆R ≈ 1. The K-factors then decrease steadily before reaching a constant value of K ≈ 1(2) for ∆R ≥ 1.4. The large values of the K-factors at small separations can be understood as follows. For this configuration to exist, the j 1 −γ 1 system is mostly recoiling against the j 2 − j 3 − γ 2 system. Thus, the large K-factors are due to this three-jet contribution and almost-vanishing two-jet cross sections. Note that, the three-jet contribution is only calculated at LO. This also explains why the scale-uncertainty bands are large at small separations. The ∆R j2γ1 distributions in the bottom-right plot show that the events are more uniformly distributed compared to the ∆R j1γ1 distributions and this happens already at LO. As a consequence, the K-factors are more moderate. For the EW-induced channel, the K-factor increases steadily from 0.9 to 1.2 for ∆R ∈ (1, 6). For the QCD case, the K-factor increases more rapidly from 0.5 to 2.0 for ∆R ∈ (1, 3), then suddenly changes its behavior to be rather constant varying from 2.0 to 2.5 for ∆R ∈ (3, 6). This sudden change at ∆R ≈ 3 is also visible in the ∆R j1γ1 case, but to a much lesser extent.
These findings are confirmed in Fig. 9 where the ∆φ jiγ1 -separation (top row) and ∆y jiγ1 -separation (bottom row) of the hardest photon and tagging jets are shown. One can observe that the region of ∆R > 3 is mostly explain by ∆y while the ∆R < 3 region is an interplay of the two contributing observables ∆y and ∆φ. It is clearly visible in the upper plots the three-jet contribution being maximal at around ∆φ j1γ1 = 0 and ∆φ j2γ1 = π. For the QCD (EW)-induced mechanism, the K-factor is 16(4) at ∆φ j1γ1 = 0.5, reaching values larger than 50(24) for ∆φ j1γ1 < 0.25. In the upper right plot, we observe in the ∆φ j2γ1 distribution that the Kfactors are more moderate with values up to 4.3 (1.6) at π, however, the relative relevance is higher since the maximum K-factor is reached for the QCD-induced sample in the dominant region of the differential distribution. This region is completely dominated by three jet events, computed only at LO, explaining the larger scale uncertainties. These large K-factors highlight the relevance of further radiation which can be studied at parton level at NNLO or including parton-shower effects and merging different jet multiplicities at NLO. We leave for future work the study of parton-shower effects in the framework of VBFNLO and Herwig. Additionally, this observable discriminates the QCD-and EW-induced processes in the region of ∆φ j2γ1 > 2. To study the discriminant capacities of the observable, in Fig.10, we show the NLO significance EW/ √ QCD calculated with luminosity L = 1 fb −1 on the cuts m j1j2 > m cut j1j2 and ∆φ j2γ1 < ∆φ cut j2γ1 . On top of this, the other default cuts are used. We observe a maximum of 4.43 at around ∆φ j2γ1 ≈ 2.51 and m j1j2 ≈ 1320 GeV, which has to be compared with the maximal value of 4.31 without the ∆φ j2γ1 cut in Fig. 4. cut is applied. While the significance do not increase considerably, we consider that this cut should be applied in the search for new physics since the region of the phase space removed is mainly dominated in the QCD-induced mechanism by three or more jet events as previously discussed, and, thus, more sensitive to higher order corrections. We note that after applying a cut of ∆φ j2γ1 < 2 the integrated QCD K-factor decreases from 1.78 to 1.43 and at the differential level from 4.21 at π to 1.93 at ∆φ j2γ1 = 2. This effect is much less pronounced in the EW-induced mechanism with the changes of K-factor from 0.99 to 0.95 and from 1.66 to 0.99 at the integrated and differential cross-section level, respectively.
Finally, we show in Fig. 11 the m γγ distribution in the presence of anomalous gauge couplings. As explained in Refs. [17,48], dimension-6 and 8 operators are included for the set of EW-induced pp → V V jj processes in the VBFNLO program using the effective Lagrangian approach. This Lagrangian reads where the operators have been defined in Refs. [48,54,66,67]. For the illustration in Fig. 11 we turn on only the following operator whereB µν = ig (∂ µ B ν − ∂ ν B µ )/2 with g and B µ being the coupling and the gauge field associated with the U (1) Y group, as defined in Ref. [48]. We note that other operators defined in Refs. [17,48] are included as usual in the EW-induced γγjj process. However, not all dimension-6 and 8 operators are included. For example, With the values of f T8 /Λ 4 = 600 and 1200 TeV −4 satisfying the current experimental bound [30] (without introducing form factors), we see in Fig. 11 that the di-photon invariant mass distribution is very sensitive to this operator at high energies.

IV. CONCLUSIONS
Results at NLO QCD for photon pair production in association with two jets via vector boson scattering has been presented for the first time in this paper. We also showed results for the QCD-induced mechanism finding good agreement with the previous calculations. In order to guarantee the validity of the VBS approximation we remove the s-channel contributions using the cut of Eq. 6 on top of a tight VBS-cut setup. With these cuts, the VBS approximation is almost identical to the full EW result and the EW-QCD interference is negligible. We have investigated the dependence of the cross sections on the photon-isolation parameters. Our results show that choosing the cone radius δ 0 = 0.4 and = 0.05 as suggested in Ref. [58] in the context of inclusive diphoton production is also good for the γγjj channel.
To increase the signal versus background ratio S = EW/ √ QCD, we studied the dependence of the cross sections for the two production mechanisms on the two typical VBS cuts m j1j2 and ∆y j1j2 . We find that tight VBS cuts, in particular a large m j1j2 cut, are needed in order to optimize the significance, which turns out to be higher than for other VBS processes. We used as a default setup m j1j2 > 800 GeV, |y j1 − y j2 | > 3 and the two photons are required to be between the rapidity gap of the two tagging jets to obtain a significance of about 4. A higher significance of 4.3 is found for m j1j2 > 1300 GeV and ∆y j1j2 > 3 with an integrated cross section of σ EW NLO = 13.35 fb. A big plateau with S > 4.25 is identified for m cut j1j2 ∈ (1, 1.6) TeV and ∆y cut j1j2 < 3.5. Furthermore, we observe that applying a cut of ∆φ cut j2γ1 < 2.5 not only increases the significance a little bit to 4.43, but, as well, reduces drastically the impact of higher-order QCD corrections.
For the EW-induced process, using the default cuts, different scale choices have been studied and we find that the scale dependence is significantly smaller for all choices when the NLO corrections are included -from about 15% at LO to a few percent or less at NLO QCD when varying both scales simultaneously µ F = µ R by a factor of 2 around the central scale. Additionally, at the central scale, we observe that the integrated NLO QCD predictions for different scale choices are as well consistent among each other at the percent level, while at LO larger differences up to 10% level are visible. This highlights the relevance of the NLO QCD predictions. We find that the scales H T /2 and √ p T,j1 p T,j2 provide the most stable results with K-factors close to 1 at the integrated cross section level. With our default scale, H T /2, we studied the effect of the NLO QCD corrections at the differential cross section level. With our VBS default cuts, for the VBS channel, corrections are small, at the a-few-percent level, for EW observables in the whole spectrum while larger corrections up to 50% are visible for jet observables. These large corrections occur in the region of phase space where the LO cross section is suppressed due to kinematic reasons or when the p T of the tagging jets are small. For the QCD-induced mechanism, corrections can be much larger reaching K-factor of 10 or higher (see the ∆φ j1γ1 distribution).
In addition, we have included in our code dimension-6 and 8 operators involving EW gauge bosons and the Higgs field using an effective Lagrangian framework. We have shown, as an illustrative example, in Fig. 11 the capability of our program to study anomalous-quarticgauge couplings at NLO QCD. The code will be available in the next release of the VBFNLO program or upon request.
With the results obtained, we see that one of the main challenges to achieve a precise EW γγjj measurement is to reduce the theoretical uncertainties of the QCD process. Since the NNLO corrections for this 2 → 4 process is unlikely to be available in the near future, further studies including parton-shower effects or using the all-order resummation discussed in Ref. [68] for the case of Hjj can be valuable.