Precision study on $ZZ\gamma$ production including $Z$-boson leptonic decays at the ILC

We report on the precision predictions for the $e^+e^-\to ZZ\gamma$ process including $Z$-boson leptonic decays at the ILC in the standard model (SM). The calculation includes the full next-to-leading (NLO) electroweak (EW) corrections and high order initial state radiation (h.o.ISR) contributions. We find that the NLO EW corrections heavily suppress the LO cross section, and the h.o.ISR effects are notable near the threshold while become small in high energy region. We present the LO and the NLO EW+h.o.ISR corrected distributions of the transverse momenta of final $Z$-boson and photon as well as the $Z$-pair invariant mass, and investigate the corresponding NLO EW and h.o.ISR relative corrections. We also study the leptonic decays of the final $Z$-boson pair by adopting the {\sc MadSpin} method where the spin correlation effect is involved. We conclude that both the h.o.ISR effects and the NLO EW corrections are important in exploring the $ZZ\gamma$ production at the ILC.


INTRODUCTION
The standard model (SM) is the most successful particle physics model until now, since its theoretical predictions are consistent with high energy experimental results excellently. In the SM, electroweak symmety breaking (EWSB) is achieved by introducing the Higgs mechanism, which gives masses to the elementary particles and implies the existence of a SM Higgs boson. A giant step of particle physics was made with a new boson with mass of about 126 GeV observed by both ATLAS and CMS collaborations at the Large Hadron Collider (LHC) in July 2012 [1,2]. The properties of this new boson are very well compatible with the SM Higgs boson but leave the room for new physics. One of the next important tasks is to investigate and measure the nature of this new particle, and finally to determine whether it is really the SM Higgs boson. The International Linear Collider (ILC) is an ideal machine to complete this task.
The H → Zγ channel is a rare decay mode but remarkable to discover the nature of the Higgs boson. This channel does not occur at tree level but is induced by loop diagrams mediated by a heavy quark loop [3] or a W -boson loop [4], just like the H → γγ decay. The H → Zγ decay mode may provide hint for new physics beyond the SM, because if new particles circulate in the loop or H is a non-SM scalar boson, the H → Zγ decay rate will be different from the SM prediction. At the ILC, the Higgs-strahlung, e + e − → ZH, is the predominantly process for Higgs production. Therefore, the most serous and irreducible background arises from the e + e − → ZZγ process in the search of Higgs boson via H → Zγ decay channel at the ILC.
The multiple gauge boson production plays an important role in probing the gauge self-couplings and would be helpful for verification of the non-Abelian gauge structure of the SM. The direct measurements of the quartic gauge boson couplings (QGCs), which can provide a connection of the mechanism of EWSB, require the theoretical study of triple gauge boson production. In the last few years, the calculations up to the QCD next-to-leading order (NLO) on triple gauge boson production at hadron colliders in the SM have been completed [5][6][7][8][9][10][11][12][13]. The NLO EW corrections to W W Z and W ZZ productions at the LHC have been provided in Refs. [14,15]. The NLO EW corrections to W W Z, ZZZ, W W γ and Zγγ productions at the ILC including the h.o.ISR contributions are calculated in Refs. [16][17][18][19][20].
The ZZγ production at the ILC is not only an important background process to study the nature of Higgs boson, but also a signal process to explore the ZZZγ and ZZγγ anomalous QGCs. These anomalous QGCs vanish at tree level in the SM and might provide a clean signal of new physics if any deviations from the SM predictions are observed. The effects of anomalous QGCs in ZZγ production at the LEP, ILC and CLIC were phenomenologically studied in Refs. [21][22][23][24]. The theoretical leading order (LO) predictions in the SM were provided in Ref. [23], while the EW corrections to e + e − → ZZγ, which would be necessary to match the ILC experimental accuracy, are still missing.
In this paper, we investigate the full NLO EW corrections and the h.o.ISR contributions to the ZZγ production at the ILC in the SM. The rest of the paper is organized as follows: In the following section we present the NLO EW and h.o.ISR analytical calculations for the e + e − → ZZγ process. The numerical results and discussions are given in Section III. Finally, we give a short summary.

CALCULATION SETUP
In our calculation, we use the 't Hooft-Feynman gauge and apply FeynArts-3.7 package [25] to generate automatically the Feynman diagrams. The corresponding amplitudes are subsequently reduced by using FormCalc-7.4 program [26]. Because of the smallness of the electron mass, we neglect the contributions from the Feynman diagrams which involve the Higgs/Goldstone-electron-positron Yukawa couplings. We adopt a mixed scheme for the electromagnetic coupling, just as our previous work [19], i.e., the couplings related to the external photons are fixed with α = α(0) in the α(0)-scheme and the others with in the G µ -scheme. Then the LO and NLO EW corrected cross sections for e + e − → ZZγ are of O(α 2 Gµ α(0)) and O(α 2 Gµ α(0) 2 ), respectively.

Virtual correction
The NLO EW correction to the e + e − → ZZγ process consists of virtual loop correction and real photon emission correction. The virtual EW correction involves 1382 diagrams which can be classified into self-energy (64), triangle (722), box (507), pentagon (56) and counterterm (33) graph groups. The amplitude for these one-loop Feynman diagrams contains both ultraviolet (UV) and infrared (IR) singularities. We adopt the dimensional regularization scheme to regularize the UV divergences in loop integrals. The relevant renormalization constants are detailed in Refs. [27,28] by using the on-mass-shell conditions. After performing the renormalization procedure, the UV divergences are canceled exactly. Since the logarithmic divergences contributed by light quarks to the photon wave-function renormalization constant have been absorbed in α Gµ , we should subtract these contributions from the virtual correction in order to avoid double counting, i.e., the electric charge renormalization constant is given in the G µ -scheme as where δZ α(0) e is the electric charge renormalization constant in the α(0)-scheme and ∆r comes from the weak corrections to muon decay [29]. The IR singularities are regularized by using infinitesimal fictitious photon mass. After adding the contribution of real photon emission process, the IR singularities are canceled. Finally, a UV-and IR-finite EW correction can be obtained.
We adopt the LoopTools-2.8 package [26] for the numerical calculations of the scalar and tensor integrals, in which the n-point (n ≤ 4) tensor integrals are reduced to scalar integrals recursively by using Passarino-Veltman algorithm and the 5-point integrals are decomposed into 4-point integrals by using the method in Ref. [30]. In our previous work [19,31], we addressed the numerical instability originating from the small Gram determinant (detG) scalar 4-point integrals [18]. In order to solve this problem, we developed the LoopTools-2.8 package and checked the results with those by using OneLoop package [32] to verify the correctness of our codes.
The one-loop Feynman diagrams for e + e − → ZZγ with possible Higgs resonance are demonstrated in Fig.1. The interference between the amplitudes for these one-loop Feynman diagrams and the LO diagrams leads to a propagator factor of H +iMH ΓH ) . The contribution of the imaginary part is so tiny that can be ignored in the total NLO EW correction.

Real photon emission correction
Technically, to extract the IR singularities from the real photon emission correction and combine it with the virtual contribution, we employ the dipole subtraction (DS) method [33]. In the DS approach, the IR-finite real correction is obtained by subtracting an auxiliary function from the squared amplitude for the real photon emission process before integrating over phase space. Because the subtraction function have the same singularity structure as the squared amplitude, this phase space integration can be performed numerically. In order to obtain an unchanged result, the subtracted term is added again after analytical integration over the bremsstrahlung photon phase space. The dipole subtraction formalism is a process independent approach, which was first presented by Catani and Seymour for QCD with massless unpolarized partons [34][35][36] and subsequently was generalized to photon radiation off charged fermions with arbitrary mass by Dittmaier [37]. In our calculations, we follow the approach of Ref. [37] and check the independence on the parameter α ∈ (0, 1], which is introduced to control the size of dipole phase space [38,39].

High order initial state radiation
The emission of photon collinear to the incoming electron or positron, i.e., initial state radiation (ISR), induces the collinear IR quasi-singularities due to the smallness of the electron mass. The virtual correction can partially cancel the ISR collinear IR quasisingularities. The left ones would lead to large radiative corrections in the form α n log n (m 2 e /Q 2 ) at the leading logarithmic (LL) level. In order to achieve an accuracy at the few 0.1% level, the high order contributions from this part beyond NLO have to be taken into account. According to the mass factorization theorem, the LL initial state QED correction can be expressed as a convolution of the LO cross section with structure functions by using the structure function method [40,41], where x 1 and x 2 denote the momentum fractions carried by the incoming electron and positron just before the hard scattering, Q 2 is the typical scale at which the hard scattering occurs chosen as the colliding energy √ s in our calculation, and the LL structure functions Γ LL ee (x, Q 2 ) are detailed in Ref. [41] up to O(α 3 ). In summing the contribution from Eq.(2) with the NLO EW corrected result, we must subtract the LO and one-loop contributions to avoid double counting. The explicit expression for the subtracted terms are presented in Ref. [40]. In the following, the subtracted ISR effect is called the high order ISR (h.o.ISR) contribution beyond O(α). We define the summation of the NLO EW corrected cross section and the h.o.ISR contribution as the total EW corrected result.

Input parameters and event selection criterion
The relevant input parameters used in our calculation are taken as [42]: where the current mass values of light quarks (all quarks except t-quark) can reproduce the hadronic contribution to the photonic vacuum polarization [43]. The Higgs boson mass is taken as M H = 125.09 GeV [44], and its decay width is obtained by using the HDECAY program [45]. The Cabibbo-Kobayashi-Maskawa (CKM) matrix is set to be unit matrix.
There is only one photon in the final state of the process e + e − → ZZγ at the LO, while at most two photons up to EW NLO. We apply the Cambridge/Aachen (C/A) jet algorithm [46] to the photon candidates. For a two-photon event originating from the real emission correction, we merge them into one new photon with four-momentum p ij,µ = p i,µ + p j,µ when the separation R = ∆y 2 + ∆φ 2 of the two final photons satisfies the condition of R < 0.4, where ∆y and ∆φ are the differences of rapidity and azimuthal angle between the two photons, then we call this event as a "one-photon" event, otherwise it is considered as a "two-photon" event. In our calculations, we collect all the "one-photon" and "twophoton" events, and at least one photon is required to satisfy the constraints Thereby we exclude the inevitable infrared (IR) singularity in the LO calculation. In the "two-photon" event, when both the two photons pass the cuts in Eq.(4) we call the photon with the largest transverse energy as the leading photon and the another one as the sub-leading photon.

Total cross section
In Fig.2(a), we present the LO and total EW corrected integrated cross sections as the functions of the colliding energy √ s for the e + e − → ZZγ process in the SM, and in Fig.2 Fig.2(b), we can see that the NLO EW relative corrections near the threshold are very large. That is because for energies near the threshold the virtual EW corrections are enhanced by Coulombic photon exchange between the electron and positron in the initial state. When the photon momentum approaches to zero a so-called Coulomb singularity occurs. At high energy range the NLO EW relative correction is also notable and increases slowly with the increment of √ s. That is due to the Sudakov logarithms from the virtual exchange of soft or collinear massive weak gauge bosons which dominate the weak corrections at high energies. As indicated in Fig.2(b), the h.o.ISR effect beyond O(α) is significant near the threshold whose relative correction can reach 11.81% at √ s = 200 GeV, but goes down and approaches the value lower than 1% when √ s > 250 GeV. In order to show the results more explicitly, we present some representative numerical results of the LO and total EW corrected cross sections, and the corresponding NLO EW and h.o.ISR relative corrections in Tab.I.

Kinematic distributions
We discuss the kinematic distributions of final produced particles for the ZZγ production process in this subsection. The LO and total EW corrected transverse momentum distributions of the final Z-boson (i.e., dσLO dp Z T and dσEW dp Z T ) at the √ s = 500 GeV ILC are provided in Fig.3(a). There we pick the p Z T of each   of the two Z-bosons as an entry in this histograms, then the final result of the differential cross section is obtained by multiplying factor 1/2. From this figure, we can see that the LO differential cross sections are enhanced by the total EW correction when p Z T ≤ 45 GeV, while suppressed in the rest of plotted region. The transverse momentum distributions of the leading photon are plotted in Fig.4(a). It shows that both the LO and total EW corrected p γ distributions for the leading photon decrease continuously with the increment of p γ T , and the LO p γ T distributions are always suppressed by the total EW correction. We also present the corresponding NLO EW and h.o.ISR relative correction distributions of p Z T and p γ T in Fig.3(b) and Fig.4(b), respectively. From these figures we can see that the h.o.ISR relative correction is usually very small, but when p Z T (p γ T ) approaches to 240 (225) GeV, the h.o.ISR relative correction can reach 8.3% (8.8%). Due to the Sudakov effect, the absolute size of the NLO EW relative corrections continuously grow up with the increments of p Z T and p γ T at high p T regions. In Fig.5(a), we depict the LO and total EW corrected distributions of the Z-pair invariant mass M ZZ . The corresponding NLO EW and h.o.ISR relative correction of M ZZ distribution are also presented in Fig.5(b). The differential cross sections of M ZZ are drawn in the range of M ZZ ∈ [2M Z , 485 GeV], where the upper limit on M ZZ is determined by the colliding energy and the transverse momentum lower cut on the final photon. As displayed in Fig.5(a), the LO distribution is enhanced by the total EW correction in small M ZZ region, while suppressed when M ZZ > 370 GeV. The obvious NLO EW correction in large M ZZ region, which can amount up to −61.5% for M ZZ ≃ 485 GeV, can be attributed to the Sudakov logarithms from the virtual corrections. Now we consider the leptonic decays of the final Z-boson pair and study the e + e − → ZZγ → l + 1 l − 1 l + 2 l − 2 γ(+γ) (l 1 , l 2 = e, µ) process by adopting the following two methods within narrow width approximation (NWA) to generate the subsequent decays.
(1) The naive NWA method. In this method, Z-boson is treated as an on-shell particle and its spin information is dropped.
(2) The MadSpin method. In this method, the improved NWA is adopted and implemented in the MadSpin package, which is part of Mad-Graph5 aMC@NLO [49] and can be used to generate the heavy resonance decay taking into account the spin correlation and finite width effects to a very good accuracy, just as demonstrated in our previous papers [15,50].
In order to display the spin correlation effect and the accuracy of Madspin approximation, we calculate the full amplitude for the signal process e + e − → ZZγ → l + 1 l − 1 l + 2 l − 2 γ at LO by taking full off-shell effects into account. In Fig.6(a), we present the LO distributions of the transverse momentum of the final negative charged lepton after Z-boson decays by applying both the naive NWA and MadSpin methods separately, which are labeled as dσNWA dp l − T and dσ MadSpin dp l − T , respectively. The distributions are depicted by entering for each p l − T and then normalising by a factor 1/2. We also plot the same kinematic distribution for the signal process by taking full off-shell effects ( dσ full dp l − T ) into account in Fig.6(a). The relative deviations are given in Fig.6(b). From these figures, we can see that the results of the naive NWA method sharply deviate from the full off-shell results (at most 25%), while the MadSpin method is less deviated (less than 5%). We see clearly that the spin correlation effect is manifested obviously in the LO p l − T distribution, and the MadSpin program is an efficient tool in handling the spin-entangled decays of resonant Z-boson in an accurate way.
We present the LO and total EW corrected distributions of the final negative charged lepton transverse momentum in Fig.7(a), and the corresponding NLO EW and h.o.ISR relative corrections in Fig.7(b). We can see that both the LO and total EW corrected p l − T distributions reach their maxima at the position of p l − T ∼ 30 GeV and the LO distributions are suppressed significantly by EW correction in the whole plotted region. In above discussion we didn't mentioned about the QED final state radiation (FSR) from the outgoing charged leptons, which could enhanced the results in the "bare" lepton scheme in measuring final leptons due to large logarithms form α n ln n (m 2 l /s) terms induced by the small lepton mass. But in the "dressed" lepton scheme we obtain the invariant mass m ll and transverse momentum p Z T by recombining the final state leptons with radiated photons within a cone, e.g.,∆R < 0.1 (called "dressed" leptons) [51]. Normally, final electrons are detected in calorimeters, photon recombination is automatically involved in the reconstruction from electromagnetic showers. While, muons can be observed as "bare" leptons from their tracks in the muon chambers, but in order to reduce FSR corrections we can sometimes reconstruct the observed muons as dressed leptons via photon recombination [52]. In this "dressed" lepton scheme, the mass singular FSR effects vanish and the resulting cross section does not depend on the mass of charged lepton, i.e. the reconstructed lepton looks universal (at least electrons and muons). In this work we study the e + e − → ZZγ → l + 1 l − 1 l + 2 l − 2 γ(+γ) process by using the "dressed" lepton scheme.

SUMMARY
The e + e − → ZZγ process is very important for understanding the nature of the Higgs boson and exploring the ZZZγ and ZZγγ anomalous QGCs. In this paper, we report on the full NLO EW corrections and the h.o.ISR contributions to the ZZγ production at the ILC in the SM. We analyze the EW quantum effects on the total cross section and find that both the NLO EW and total EW corrections suppress the LO cross section when √ s goes up from 200 GeV to 1 TeV. Due to the Coulomb singularity effect, the NLO EW relative correction is very large near the threshold (e.g., δ N LO = 42.72% with √ s = 200 GeV). The h.o.ISR effect beyond O(α) is also distinct near the threshold (e.g., the relative correction is 11.81% when √ s = 200 GeV), while becomes small at the high colliding energy region. We provide the LO and total EW corrected p Z T , p γ T and M ZZ distributions and investigate the corresponding NLO EW relative correction and h.o.ISR relative correction. We find that the h.o.ISR effect becomes notable near the maximum value of p Z T , p γ and M ZZ . We also investigate the leptonic decays of the final Z-boson pair by adopting the MadSpin program to include the spin correlation effect and find that the LO p l − T distributions are suppressed in the whole plotted region. Ascribed to the Sudakov effect, the NLO EW correction becomes very large with the increment of these kinematic variables. We conclude that both the h.o.ISR and the NLO EW corrections are worth being taken into account in precision measurement of the ZZγ production at the ILC.