ZA production in vector-boson scattering at next-to-leading order QCD

Cross sections and differential distributions for ZA production in association with two jets via vector boson fusion are presented at next-to-leading order in QCD. The leptonic decays of the Z boson with full off-shell effects and spin correlations are taken into account. The uncertainties due to different scale choices and pdf sets are studied. Furthermore, we analyze the effect of including anomalous quartic gauge couplings at NLO QCD.


I. INTRODUCTION
Electroweak photon production in association with two charged leptons and two jets is an important channel at the LHC since it provides information on weak boson scattering. It is also sensitive to beyond standard model (BSM) physics via anomalous gauge boson couplings.
At LO, there are two mechanism to produce + − γjj events at the LHC. The QCD-induced of order O(α 2 s α 3 ) and the EW-induced of order O(α 5 ) which is further classified into the s-channel contributions, given mainly by tri-boson production with a subsequent hadronic decay of one of the vector bosons, and the t/u-channel vector boson fusion (VBF) contributions.
The QCD-induced mechanism is considered to be an irreducible background of the EW mechanism due to the lack of weak boson scattering and quartic gauge boson couplings, V V → V V . Despite the α 2 s /α 2 enhancement, its contribution is comparable in typical VBF searches. The interference effects among the different mechanism/channels are expected to be small in LHC measurements. A dedicated study of these effects was carried out in Ref. [1] for same sign W W jj production, where interference effects are expected to be larger due to the absence of gluon initiated processes and the fixed chirality of the quark lines.
The NLO QCD corrections for the QCD-induced mechanism were given in Ref. [2]. The corrections are moderate, if adequate central scales are chosen, but phase space dependent and lead to a significant reduction of the scale uncertainty.
The NLO QCD corrections of the tri-boson production processes were first computed including the leptonic decays in Refs. [3][4][5] and afterwards including the hadronic decays in Ref. [6,7]. They turned out to be large, around 70%, and not covered by the scale uncertainties. This is due to logarithmically enhanced configurations [8,9] and new channels opening up at NLO.
Both the QCD-induced and tri-boson production processes are available in the VBFNLO package [7,10]. In this paper, we provide results at NLO QCD for the VBF t/u-channel for the processes pp → + − γjj + X, "Z γjj", (1) pp → ννγjj + X, "Z ν γjj", focusing mainly on the charged leptonic channel. Representative Feynman diagrams are shown in Fig. 1. The bulk of the cross section for both processes comes from regions of the phase space where the intermediate Z boson is approximately on-shell. For simplicity, in the following, we often refer to the processes by "Z γjj" or "Z ν γjj" production, although we will consider all off-shell effects, non-resonant diagrams and spin correlations. Furthermore, using an Effective Theory (EFT) approach, we will study, at NLO QCD, BSM effects due to anomalous quartic gauge couplings.
The processes considered here have been implemented in VBFNLO [7,10], a parton level Monte Carlo program which allows the definition of general acceptance cuts and distributions.
The rest of the paper is organized as follows: in section II the calculational setup as well as the checks performed to ensure the correctness of the calculation are given. In section III, results at the integrated cross section level, for differential distributions, as well as for anomalous couplings studies are presented. Finally, we present our conclusions in section IV.

II. CALCULATIONAL SETUP
The calculation method of the "Z γjj" and "Z ν γjj" production follows closely the one of pp → + 1 − 1 + 2 − 2 jj + X [11] (called from now on Z Z jj production for simplicity) and other VBF channels implemented in VBFNLO.
We work in the VBF approximation and, therefore, only the t/u-channel Feynman diagrams, neglecting the interference between them, are considered. The s-channel contributions at NLO are accessible in VBFNLO via "ZAV " production, with V ∈ (W, Z, γ) decaying hadronically. Their size as well as the interference effects are small once typical VBF cuts are applied.
We work in the five flavor scheme. For the potentially resonating massive vector bosons, we use a modified complex mass scheme as implemented in MadGraph [12].
Technically, to obtain the NLO QCD corrections to EW t/u-channel contributions of Z γjj and Z ν γjj production, we adapt the code with some modifications from the process pp → + 1 − 1 + 2 − 2 jj + X and pp → + − νν jj + X implemented in VBFNLO, respectively. Here, to be self-contained, we give a brief description of the method used for the Z γjj process. The Z ν γjj channel is calculated analogously.
To compute the amplitudes, we use the helicity formalism of Ref. [13]. This allows us to factorize amplitudes into a QCD part and electroweak factors. The latter contain not only the decay currents for V → + − andṼ → + − γ with V,Ṽ ∈ (Z, γ * ), but also "leptonic tensors" for V V /W + W − → + − /γ/ + − γ, containing the scattering of the t-channel vector bosons connecting the quark lines. The EW currents and tensors are calculated only once per phase space point using the routines of the HELAS package [14]. Afterwards, the full LO amplitudes for all sub-processes, q 1 q 2 → q 3 q 4 + − γ and crossing related ones, can be easily obtained using the pre-calculated structures. Similarly, we obtain the real emission amplitudes by adding an additional gluon emission to a quark line.
For the virtual corrections, we do not consider graphs with a gluon exchange between the two quark lines. Due to the color structure of the amplitude, these would only give non-vanishing contributions for the interferences of t-and u-channel diagrams, which are phase-space suppresses and neglected in the VBF approximation. Thus only corrections to a quark line (the upper or the lower one) with up to three bosons emitted have to be considered. We make use of the routines Boxline and Penline computed in Ref. [15], which respectively combine all the loop diagrams to a quark line with two or three bosons emitted in a fixed order permutation of the external bosons. The full amplitude is obtained by including all physically allowed permutations.
We use the Catani-Seymour formalism [16] to deal with the infrared divergences and we follow Ref. [17] to obtain individual factorization and renormalization scales for each of the two quark lines. This is possible because they each form a gauge invariant sub-set.
To deal with the real photon in the final state, we implement the Frixione smooth cone isolation cut [18], which preserves IR safety, without the need of introducing photon fragmentation functions. Additionally, we split the phase space integration into two separated regions to improve the convergence of the Monte Carlo integration. The regions are generated as double EW boson production with (approximately) on-shell Z → decay or as Z production with Z → γ three-body decay, respectively, and they are chosen according to whether m( γ) or m( ) is closer to M Z . The final result is obtained by adding the two integrals. The presence of intermediate off-shell photons, which are far from on-shell for typical lepton cuts, does not pose numerical problems and their contribution is integrated with the Z ones. For the Z ν γjj production channel, this phase space splitting is not necessary.
To ensure the correctness of our results, we have crosschecked our LO and real matrix elements with Madgraph [12] and we compared the integrated cross sections with Sherpa, finding agreement at the machine precision and at the per mille level, respectively. For the real emission contributions, we need to subtract the s-channel contributions explicitly which are included in Sherpa, otherwise, agreement at the few percent level is found for typical VBF cuts. Thus, the neglected s-channel contributions only give a noticeable contribution to the real emission, leading to an increase at the few per mille level of the total NLO QCD cross section. Hence, they can be safely neglected. In the following, we discard these contributions. However, they can be included easily within the VBFNLO package at NLO QCD, by adding the cross sections for triple electroweak boson production.
Furthermore, we have checked the convergence of the Catani-Seymour subtraction and, for the virtual contributions, the factorization of the poles, as well as gauge and parametrization invariance [15].
The numerical stability of the calculation of the virtual amplitudes is controlled with the use of Ward identities [15]. The amplitude is set to zero if they are not satisfied with a precision better than one per mille. The fraction of phase space points rejected is around the per mille level and thus, the error induced by this procedure is negligible since the total virtual contributions, after analytic cancellation of infrared divergences, are at the level of a few per cent.

III. PHENOMENOLOGICAL RESULTS
In the following, we present results for the LHC mainly for the Z γjj channel. Two lepton families are included in the results presented. Since we apply the same cuts for the first two families, we compute the process pp → e + e − γjj + X, and multiply the result by two. As EW input parameter, we use the Fermi constant as well as the Z and W mass and derive the remaining EW parameters  via tree-level relations, i.e., we use The widths of the bosons are taken as As default, we use the PDF4LHC15 nlo 100 [19] PDF set both at LO and NLO. The jets are defined with the antik t algorithm [20] with radius R = 0.4 and are required to have a transverse momentum p T,j > 30 GeV and rapidity |y j | < 4.5. The jets are ordered by transverse momenta and the tagging jets at NLO are defined as the two hardest jets. To simulate typical VBF searches and LHC detector capabilities, we use The last cut in Eq. (4) eliminates the singularity arising from a virtual photon, γ * → + − . Additionally, we impose the typical VBF cuts on the tagging jets, Furthermore, as a photon isolation cut, following the "tight isolation accord" presented in Ref. [21], events are accepted if they satisfy with δ 0 = 0.4 and efficiency, = 0.05. As default factorization and renormalization central scale, we chose the absolute value of the momentum transfer Q i from quark line "i" to the EW process, A. Total Cross Section In the following, we present results for the integrated cross section at the LHC at a center of mass energy of 13 TeV for different scale choices and PDF sets. In the left panel of Fig. 2, we vary independently the factorization and renormalization scale µ = ξµ 0 in the range ξ ∈ (0.1, 10) around the central scales µ 0 = M Z (orange line) and µ 0 = Q i (blue line). At ξ = 1, we find for µ 0 = Q i σ LO = 2.9378(7) −8% +9% fb and σ NLO = 2.837(1) −0.3% −1% fb with a K-factor, defined as the ratio of the NLO over the LO predictions, of K = 0.97. Correspondingly, for µ 0 = M Z we find σ LO = 3.083(2) −8% +10% fb and σ NLO = 2.848(4) +0.3% −2% fb with a K-factor of K = 0.92. The upper(sub)-scripts correspond to the scale uncertainty taken at ξ = 2(ξ = 0.5) being of order 10% at LO and few percent at NLO. The numbers in parenthesis quote MC statistical errors. Note that at leading order, we only have factorization scale uncertainties and that the differences of the LO and NLO predictions of about 5% and 0.5%, correspondingly, at the central scale for the two different scale choices are contained in the scale uncertainties. This is, however, not always the case at LO as we can see in the right panel of Fig. 2 where we have plotted additionally, setting µ F = µ R = ξµ 0 for simplicity, the curves for µ 0 ∈ (H T , H T * , E T ), which are often used in the corresponding QCD V V jj induced processes, and are defined as: with V i ∈ (Z, γ). m Vi denotes the invariant mass of the corresponding leptons (m Vi = 0 for on-shell photons) and y 12 = (y 1 + y 2 )/2 the average rapidity of the two hardest (or tagging) jets, ordered by decreasing transverse momenta. E T (jj) and E T (V V ) stand for the transverse energy of the two tagging jets and of the V V system, respectively. In the last two scale choices of Eq. (8), the first term interpolates between m jj and p T,jets for large and small ∆y jj = |y 1 − y 2 | values, respectively. First, we note that the LO predictions for the different central scale choices are not covered by the scale uncertainties. At, ξ = 1, we found at LO differences of about 40% for µ 0 = M Z and µ 0 = H T * , while the NLO differences are about 5% for µ 0 = H T and µ 0 = H T * . Also, note, that the K-factor varies from 0.92 for µ 0 = M Z to 1.24 for µ 0 = H T * , and greatly depends on the value of the LO predictions. Even larger differences can be seen in differential distributions, thus, pointing out the necessity of using NLO predictions for obtaining robust results.
The order 5% uncertainty found above for the more extreme scale choices is, in fact, a better, but still a low estimate for the error induced by missing NNLO corrections. As was found for VBF Higgs production, cross sections and various distributions within VBF cuts are lowered by another 5 to 10% by NNLO corrections as compared to NLO results [22]. This effect is due to a wider energy flow within NNLO quark jets as compared to the NLO approximation: the narrow R = 0.4 jets capture less of the jet energy at NNLO and thus have a harder time passing the m jj > 600 GeV requirement [23]. Since this effect should be universal for all VBF processes, we expect another order 5% reduction of the fiducial cross section at NNLO, which is not accounted for by the scale considerations.
In Fig. 3, we plot the scale variation, for equal factorization and renormalization scale, i.e. µ F = µ R = ξQ i , as well as the associated PDF uncertainty at NLO for the PDF4LHC15 nlo 100 and HERAPDF20 [24] (EIG) sets. In addition, we show the scale variation for the ABM11 [25] (5 flavors), CT14 [26], MMHT2014 [27], and NNPDF30 [28] (with α s (M Z ) = 0.118) sets used at the corresponding order, if available. For NNPDF30 and ABM11, we use the NLO sets for the LO curves and for the CT14 curve at LO, we use the LO CT10 [29] sets. At LO, we observe at the central point an 8% maximum difference between the MMHT2014 and the CT10 predictions, which is of the same order as the scale uncertainty reported previously, while the error of the PDFs is around the few percent level and, thus, do not cover the uncertainty observed. At NLO, at the central point, a 6% maximum difference between HERAPDF20 and AMB11 is found. It is neither covered by the scale uncertainty, varying only one of the particular scales, nor by the PDF uncertainties, which are 2% and 1% for these PDFs, respectively. When comparing our default PDF set (PDF4LHC15 nlo 100), with a 2% PDF uncertainty, with the HERAPDF20 predictions, they almost overlap with a difference of 4.7%. The size of the PDF uncertainties is of the same order in the whole range ξ ∈ (0.1, 10).
In Table I, one finds the values at NLO obtained with the different PDFs as well as the associated asymmetric error at ξ = 1. In Fig. 4, we plot the cross sections for different energies. In addition to our default PDF set, we plot the lines for the AMB11 and NNPDF30 sets, where larger differences of about 20% are seen at a center of mass energy of 100 TeV. The bands include the scale and PDF uncertainties, which are added in quadrature. The combined uncertainty at 100 TeV rounds the 8% for the AMB11 set, dominated by the PDF uncertainty of about 7%, and 3% for the NNPDF30 and our default PDF4LHC15 nlo 100 set.
In order to remove contributions not relevant for the study of anomalous coupling effects, in Fig. 5, we plot the integrated cross section depending on the minimum value required for the invariant mass, m Zγ , of the EW system both for the Z γjj and Z ν γjj channels, neglecting the corresponding cut specified in Eq. (4). For leptonic decays of the Z-boson, two lepton generations are counted while for Z → νν only one generation is considered, which gives roughly equal cross sections for easier comparison. When applying cuts, we treat the neutrinos  the same as charged leptons, such that differences of the two processes can be entirely associated with differences in the amplitude and coupling constants. For the Z γjj channel, we show lines at LO (red) and NLO (blue), including the scale uncertainty, while only the central value at NLO (green) is shown for the Z ν γjj channel for comparison. In the middle panel, the Kfactors are shown, while in the third panel, the ratios of the Z γjj and Z ν γjj channels are plotted. We observe that imposing a minimum value m Zγ > 90 GeV, the contributions from a photon radiated off the charged leptons (radiative Z decays), which is absent in the Z ν γjj channel, start to decrease significantly since configurations of the decay Z → − + γ can only be off-shell. For values of m Zγ > 120 GeV, these contributions are almost completely gone. This is confirmed in the third panel with a ratio almost constant beyond this value. Similar observations were found in the QCD induced process in Ref. [2].
Following the recommendation of the "tight isolation accord" of Ref. [21], we set the efficiency to = 0.05. In Fig. 6, we study the dependence of the cross section on the efficiency parameter plotting the cross section varying in the range ∈ (0.01, 1). We observe mild dependencies of around 3%(6%) in the whole range shown and of order 2%(3%) in ∈ (0.01, 0.05) for δ 0 = 0.4(0.7). For comparison, we also show results for the QCD induced process. In this case, we find larger differences of about 12%(23%), with 6%(13%) differences in the ∈ (0.01, 0.05) range with δ 0 = 0.4(0.7). The milder variation found in the EW channel is expected due to the characteristic signature of the VBF processes with two forward jets and the photon produced in the central region, where little hadronic activity is expected due to the formation of a rapidity gap in VBF processes.
Finally, in Table II we give the NLO cross sections of the EW and QCD production processes at different collider energies including scale uncertainties. The number in parenthesis represents the statistical Monte Carlo integration error, while the upper and the lower numbers represent the scale uncertainty error at ξ = 2 and ξ = 0.5, respectively. With the given VBF cuts, both mechanisms are of the same order, in spite of the (α/α S ) 2 suppression of the EW channel.  In the following, we show results for differential distributions at √ s = 13 TeV for the QCD and EW induced Z γjj processes using our default settings described in section II. In the top panels, we show the EW LO (red) and NLO (dark-blue) curves, including the scale uncertainties. In light-blue, we show the QCD induced process at NLO, including its scale uncertainties. In the bottom panels, we show the corresponding EW K-factor as well as the scale uncertainty band compared to the LO result at the central scale. PDF uncertainties of the EW process are shown as hatched bands. In Fig. 7, we show in the upper row the differential distribution of the invariant mass of the electroweak system (left) and the minimum R-separation between the photon and one of the leptons (right). In the lower row, we show jet observables for the tagging jets, the dijet invariant mass (left) and the rapidity separation (right). Given the appropriate scale choice, Q i , we observed modest K-factors, close to one, in the whole spectrum, with larger variation in the rapidity separation plot, ranging from 0.90-1.10, and a drastic reduction of the scale uncertainties. In the topright plot, we show in addition the curve (green) without the m ZA > 120 GeV cut. As expected, the cut only reduces events with photons emitted close to the charged leptons. In the bottom-left plot, one can observe clearly the distinct behaviour of the invariant mass distribution of the tagging jets for the EW vs QCD channels, with a steeper fall-off of the cross section for the QCD induced process. In Fig. 8, we show the normalized centralized rapidity distribution of the reconstructed Z boson system (left) and the photon (right) with respect to the tagging jests, Whereas for the EW process, the electroweak particles are nearly exclusively produced in the central region between the two tagging jets located at z * = ±1/2, they are produced in a broader rapidity range for the QCD process. Note that the distributions are not symmetric because the jets are p T -ordered. In particular for the QCD process, larger contributions can be found in the vicinity of the hardest jet. The particular shape of the QCD distributions can be explained by kinematic configurations, where the hardest jet recoils against the EW system, and a second jet, possibly stemming from gluon radiation, is produced at large separation to fulfill the VBF cuts.

C. Anomalous Couplings
In our implementation of the Z γjj production cross section, we allow for modified gauge couplings in the framework of an effective Lagrangian where the operators of dimension 6 and 8 have first been defined in Refs. [30][31][32]. Due to minor differences in the definition of the field strength tensors in VBFNLO, our conventions for the dimension 8 operators differ from the ones given in Ref. [32]. The exact definition, as well as the corresponding conversion rules can be found in Ref. [7]. While all operators given in Refs. [30][31][32] affect the Z γjj production cross section, the operators with are of particular interest for Z γjj production since they only involve neutral gauge bosons. Hence, they can first be constrained in vector boson scattering pp → V V jj or triboson production pp → V V V of neutral gauge bosons (V ∈ (Z, A)). Current experimental constraints on these operators can be found in Refs. [33,34].  the underlying vector boson scattering process, leading to unitarity violation for large invariant masses. Unitarity of the scattering amplitude can be restored by multiplying the amplitude with a form factor of the form A different approach to unitarize the amplitude, via Kmatrix unitarization, has been explored in Ref. [35], leading to a modification of the normalized eigen-amplitudes a IJ according to In the large s limit, K-matrix unitarization leads to a behavior similar to applying a modified, complex form factor [36], and in the following we compare this complex form factor with the conventional form factor defined in Eq. (14). The form factor scales Λ F F and Λ c F F are set according to the unitarity constraint, such that the helicity combination with the largest contribution to the zeroth partial wave fulfills the unitarity condition for all vector boson scatterings V V → V A and W W → V A (V ∈ (Z, A)).
In Fig. 9, we show the dependence of the Z γjj cross section on the invariant mass m Zγ of the electroweak system for different values of f T 8 . It can be seen that well below the form factor scale, the results using the complex form factor (dashed lines) closely follows the results where no unitarization is applied. Only close to and above the form factor scale, the modified form factor leads to a significant suppression, removing the unitarity violating tail of the distribution. The different values of f T 8 influence the shape of the distribution well below the form factor scale, but for very high invariant masses they lead to identical results, given by the unitarity condition. In contrast, the conventional form factor (dotted lines) also leads to a significant reduction of the cross section well below the form factor scale, reducing the effect of the anomalous coupling in a larger region of the phase space. While Fig. 9 illustrates the differences of the two form factors over a broad range of m Zγ , it is clear that the phase-space region above the form factor scale is determined by the unitarization procedure and shouldn't be used to constrain anomalous couplings. In Fig. 10 we fo-cus on smaller invariant masses and, in addition, we show the cross section differential in the final state transverse momenta. Similar to the m Zγ distribution, we observe that the conventional form factor (dotted lines) leads to a large suppression, whereas its complex version (dashed lines) leads to results much closer to the results without unitarization. However, also for the complex form factor, the deviations from the result without unitarization start already at small transverse momentum. In particular for the transverse momentum of the softer lepton, we obtain a large suppression already for very small values of p T,l2 .
Since anomalous gauge couplings lead to an increased cross section in the tails of the distributions, experimental limits are often obtained from a comparison of the observed event count with the cross section where a lower cut on the invariant mass of the electroweak system is applied. We therefore show the dependence of the cross section on this cut in Fig. 11, where it can be seen that a large fraction of the cross section results from contributions with invariant masses above the form factor scales. We want to point out that theoretical predictions in this phase space region highly depend on the unitarization procedure. A preferred procedure to compare experimental results with theoretical predictions should therefore be a comparison based on differential distributions, restricted to invariant masses of the electroweak system well below the form factor scale.

IV. CONCLUSIONS
In this article, we have reported results at NLO QCD for VBF Zγ production, including the leptonic decay of the Z boson with all off-shell effects and spin correlations taken into account. While, at LO, the results greatly depend on the scale choice, with up to 40% differences at the central value, the NLO results reduce considerably the scale uncertainties, to the few percent level. PDF uncertainties of individual sets have been studied yielding errors of a few percent, which are propagated quite homogeneously over the available phase space. However, the central values of the predictions for the different sets differ by up to 6%, which is not covered by the combined uncertainty associated to the pdf sets and scale variation. Furthermore, we have presented results for anomalous couplings and we pointed out the necessity to constrain anomalous couplings with data well below the form factor scale only.

V. ACKNOWLEDGMENTS
FC acknowledges discussions with J. Bernabeu and financial support by the Spanish Government and ERDF funds from the European Commission (Grants No. RYC-2014-16061, FPA2014-53631-C2-1-P , FPA2014-57816-P, and SEV-2014-0398). The work of DZ was supported by "BMBF Verbundforschung Teilchenphysik" under grant number 05H15VKCCA. . The line styles are as given in Fig. 10.