Precision study of ZZγ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ZZ\gamma $$\end{document} production including Z-boson leptonic decays at the ILC

We report on the precision predictions for the e+e-→ZZγ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^+e^- \rightarrow Z Z\gamma ~$$\end{document}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+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$+$$\end{document}h.o. ISR corrected distributions of the transverse momenta of final Z-boson and photon as well as the Z-pair invariant mass, and we 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 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γ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ZZ\gamma $$\end{document} 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 symmetry 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 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 a e-mail: dpffqs@mail.ustc.edu.cn boson. The International Linear Collider (ILC) is an ideal machine to complete this task.
The H → Z γ channel is a rare decay mode but suitable 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 − → Z H, is the predominantly process for Higgs production. Therefore, the most serous and irreducible background arises from the e + e − → Z Zγ 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 Z Z productions at the LHC have been provided in Refs. [14,15]. The NLO EW corrections to W W Z, Z Z Z, W W γ , and Z γ γ productions at the ILC including the h.o. ISR contributions are calculated in Refs. [16][17][18][19][20].
The Z Zγ 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 Z Z Zγ and Z Zγ γ 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 Z Zγ 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 − → Z Zγ , 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 Z Zγ 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 − → Z Zγ process. The numerical results and discussions are given in Sect. 3. 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 Then the LO and NLO EW corrected cross sections for e + e − → Z Zγ are of O(α 2 G μ α(0)) and O(α 2 G μ α(0) 2 ), respectively.

Virtual correction
The NLO EW correction to the e + e − → Z Zγ 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 onmass-shell conditions. After performing the renormalization procedure, the UV divergences are canceled exactly. Since the logarithmic divergences contributed by the 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 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 (det G) scalar 4point 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 − → Z Zγ 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 In order to regulate it, we replace The contribution of the imaginary part is so tiny that it 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 has 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 a 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 quasi-singularities. The ones left 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 expressions 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] G μ = 1.1663787 × 10 −5 GeV 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 − → Z Zγ 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 i j,μ = 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 a "onephoton" event, otherwise it is considered a "two-photon" event. In our calculations, we collect all the "one-photon" and "two-photon" events, and at least one photon is required to satisfy the constraints p γ T ≥ 15 GeV, |y γ | ≤ 2.5.
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 the leading photon and the another one the sub-leading photon.

Total cross section
In Fig. 2a, we present the LO and total EW corrected integrated cross sections as functions of the colliding energy √ s for the e + e − → Z Zγ process in the SM, and in Fig. 2b we show the corresponding h.o. ISR and NLO EW relative corrections, defined as δ ≡ σ −σ LO σ LO . From these figures, we find that all the LO and total EW corrected integrated cross sections are sensitive to the colliding energy, and they reach their maxima at the position of √ s ∼ 350 GeV. The LO cross sections are suppressed by both the NLO EW and the total EW correction in the whole plotted √ s region. From Fig.  2b, 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 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. 2b, the h.o. ISR effect beyond O(α) is significant near the threshold whose relative correction can reach 11.81 % at √ s = 200 GeV, but it falls 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 Table 1.

Kinematic distributions
We discuss the kinematic distributions of final produced particles for the Z Zγ production process in this subsection.   Fig. 3a. 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. 4a. It shows that both the LO and the total EW corrected p (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 [47,48], which is part of MadGraph5_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 − → Z Zγ → l + 1 l − 1 l + 2 l − 2 γ at LO by taking full off-shell effects into account. In Fig. 6a Fig. 6b. 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 In the above discussion we did not mention 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 the 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, and 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, 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 − → Z Zγ → l + 1 l − 1 l + 2 l − 2 γ (+γ ) process by using the "dressed" lepton scheme.

Summary
The e + e − → Z Zγ process is very important for understanding the nature of the Higgs boson and exploring the Z Z Zγ and Z Zγ γ anomalous QGCs. In this paper, we report on the full NLO EW corrections and the h.o. ISR contributions to the Z Zγ 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 the 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., at the high colliding energy region. We provide the LO and total EW corrected p Z T , p γ T , and M Z Z 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 γ T , and M Z Z . We also investigate the leptonic decays of the final Zboson 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 of being taken into account in precision measurement of the Z Zγ production at the ILC.