1 Introduction

Probing the mechanism of electroweak symmetry breaking (EWSB) is one of the most important tasks in particle physics. In the standard model (SM), symmetry breaking is achieved by introducing the Higgs mechanism, which gives masses to the elementary particles and implies the existence of an SM Higgs boson. Therefore, to uncover the origin of EWSB and to determine whether the SM Higgs boson really exists is one of the highlights of the Large Hadron Collider (LHC) physics program [1]. In July 2013, both the ATLAS and CMS collaborations at the LHC reported that they had observed a new neutral boson with mass of around \(126\) GeV [2, 3], and this particle is tentatively identified as a Higgs boson. The more precise measurements on its properties are still going on at the LHC, but in light of the current data, its properties are very well compatible with the SM Higgs boson. However, it has been understood for a long time that there are intrinsic limitations from the ability of hadron colliders in precision measurement. The International Linear Collider (ILC) is an ideal machine to address this problem [4]. One of the major aspects of the physics program of the ILC is to make detailed precision measurements as regards the nature of the Higgs boson discovered at the LHC [4, 5].

For any observed new particle, the determination of its fundamental properties will be a primary goal. The measurement of the branching fraction of the Higgs boson decaying into two photons, \(Br(H\rightarrow \gamma \gamma )\), turns out to be an absolutely necessary ingredient in extracting the total width [6]. Besides, this measurement may possibly provide hints for new physics if the deviation from the SM prediction is larger than the measurement accuracy. At the ILC the Higgs boson is predominantly produced by the Higgs-strahlung process \(e^+e^-\rightarrow ZH\). The most serious and irreducible background for Higgs search via the \(H \rightarrow \gamma \gamma \) decay channel arises from the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process, which it is hard to get rid of and needs to be explored in depth [7].

The precision measurement of the quartic gauge boson coupling (QGC) can provide a connection to the mechanism of electroweak symmetry breaking. The anomalous QGC, such as \(ZZ\gamma \gamma \), vanishes in the SM at the tree level and might provide a clean signal of new physics, since any deviation from the SM prediction might be connected to the residual effect of electroweak symmetry breaking. The effect of \(ZZ\gamma \gamma \) coupling has been theoretically investigated at the LEP and the ILC [813]. The measurement of the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process at LEP2 by L3 Collaboration [14] shows that the anomalous \(ZZ\gamma \gamma \) coupling leads to a negligible effect at LEP energy, while it might be detectable at the ILC with higher colliding energy. Since this effect can be small and subtle, theoretical predictions for the cross section with high precision are mandatory.

At the ILC, the accuracy of the cross section measurement for the triple gauge bosons production process could reach per mille level. It is necessary to presume upon an accurate theoretical calculations to match the experimental accuracy. Thus good theoretical predictions beyond leading order (LO) are indispensable. In the last few years, a lot of work contributed to the phenomenological studies in the SM up to the QCD next-to-leading order (NLO) on triple gauge boson production processes at hadron colliders [1523]. Most recently, the calculation of the NLO electroweak (EW) correction to the \(W^+W^-Z \) production at the LHC was present in Ref. [24]. The NLO EW calculations to the \(W^+W^-Z \) and \(ZZZ\) productions at the ILC were provided in Refs. [2527], while a prediction for the \(Z\) production associated with two photons at the ILC in the NLO EW precision is still missing.

In this paper, we investigate the complete NLO EW corrections to the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process at the ILC in the SM. The rest of the paper is organized as follows: In the following section we present the LO and NLO EW analytical calculations for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process. The numerical results and discussions are given in Sect. 3. Finally, we will give a short summary.

2 Analytical calculations

The LO and NLO EW calculations for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process in the SM are presented in this section by using the ’t Hooft–Feynman gauge. We apply FeynArts-3.7 package [28] to automatically generate the Feynman diagrams and the FormCalc-7.2 program [29] to algebraically simplify the corresponding amplitudes. In our calculations we neglect the contributions of the Feynman diagrams which involve the \(H\)\(e\)\(\bar{e}, G^0\)\(e\)\(\bar{e}, G^+\)\(e\)\(\bar{\nu }_{e}\) or \(G^-\)\(\nu _{e}\)\(\bar{e}\) vertices, because the Yukawa coupling strength of Higgs/Goldstone to fermion pair is proportional to the fermion mass.

We denote the process

$$\begin{aligned} e^+(p_1)+e^-(p_2)\rightarrow Z(p_3)+\gamma (p_4)+\gamma (p_5), \end{aligned}$$
(2.1)

where \(p_i(i=1,5)\) represent the four-momenta of the initial and final particles. There are six generic tree-level Feynman diagrams for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process, and some of them are depicted in Fig. 1. The LO cross section for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process can be obtained as

$$\begin{aligned} \sigma _\mathrm{LO}=\frac{1}{2!}\frac{(2\pi )^4}{2s}\int \mathrm{d}\Phi _{3} \overline{\sum _\mathrm{spin}} |\mathcal{M}_\mathrm{LO}|^2, \end{aligned}$$
(2.2)

where \(\mathcal{M}_\mathrm{LO}\) is the LO amplitude, the factor \(\frac{1}{2!}\) comes from the two identical final photons and the bar over summation recalls averaging over the initial spins. The phase space element of the final three particles is defined as

$$\begin{aligned} {\mathrm{d}\Phi _3}=\delta ^{(4)} \left( p_1+p_2-\sum _{i=3}^5 p_i \right) \prod _{i=3}^5 \frac{d^3 \vec {p}_i}{(2 \pi )^3 2 E_i}. \end{aligned}$$
(2.3)
Fig. 1
figure 1

The tree-level Feynman diagrams for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process. The graphs with the exchange of the final two photons are not drawn

The virtual EW correction to the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process at \(\mathcal{O}(\alpha ^4)\) involves \(1003\) diagrams, including \(36\) self-energy diagrams, \(472\) triangles, \(418\) boxes, \(47\) pentagons and \(30\) counterterm graphs. The most complicated topology involved in the EW one-loop contribution contain \(5\)-point integrals up to rank \(4\), which are deduced by using the reduction method in Ref. [30]. The numerical calculations of \(n\)-point \((n \le 4)\) tensor integrals are implemented by using the Passarino–Veltman reduction algorithm [31]. We adopt mainly the LoopTools-2.8 package [29] for the numerical calculations of the scalar and tensor integrals. In order to avoid instability in the numerical calculations of the 5-point tensor integrals of rank \(4\), we developed the program coded in Fortran77 with quadruple precision for the numerical calculation of the pentagons shown in Fig. 2. The virtual EW correction to the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process can be expressed as

$$\begin{aligned} \Delta \sigma _{v}= \frac{1}{2!}\frac{(2\pi )^4}{2s}\int \mathrm{d}\Phi _{3} \overline{\sum _\mathrm{spin}} 2Re\left\{ \mathcal{M}_{v}\mathcal{M}_\mathrm{LO}^*\right\} , \end{aligned}$$
(2.4)

where \(\mathcal{M}_{v}\) is the amplitude of all the virtual EW correction Feynman diagrams.

Fig. 2
figure 2

The pentagon diagrams for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process which are calculated by using the codes with quadruple precision. The diagrams with exchanging the final two photons are not drawn

The one-loop Feynman diagrams with possible Higgs and \(Z\)-boson on-shell effects for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process are shown in Fig. 3. Due to the Landau–Yang theorem [32, 33], the contribution from Fig. 3(2) is vanished. The interference between the amplitude of Fig. 3(1) and the LO amplitude leads to a propagator factor of \(\frac{1}{(M_{\gamma \gamma }^2-M_{H}^2)}\), which is divergent in the vicinity of \(M_{\gamma \gamma }^2\sim M_{H}^2\). We regulate it by making the replacement of \(\frac{1}{(M_{\gamma \gamma }^2-M_{H}^2)} \rightarrow \frac{1}{(M_{\gamma \gamma }^2-M_{H}^2 +iM_H\Gamma _H)}\). We find that the contribution of this interference term is so tiny that it can be ignored in the total NLO EW correction.

Fig. 3
figure 3

The one-loop diagrams with possible on-shell internal Higgs or \(Z\)-boson for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process

The amplitude for all the one-loop Feynman diagrams contains both the ultraviolet (UV) and the infrared (IR) singularities. We adopt the dimensional regularization scheme [34], in which the dimensions of spinor and space-time manifolds are extended to \(D=4-2\epsilon \) to regularize the UV divergences in loop integrals, and the IR singularities are regulated by adopting infinitesimal fictitious photon mass as it commonly is applied to photon radiation in EW processes. The relevant fields are renormalized by adopting the on-mass-shell (OMS) renormalization scheme and the explicit expressions for the renormalization constants are detailed in Refs. [35, 36]. As we expect, the UV divergence contained in the loop virtual amplitude can exactly be canceled by that in the counterterm amplitude.

In order to get an IR-finite cross section for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process at the EW NLO, we consider the real photon emission process \(e^+(p_1)+e^-(p_2)\rightarrow Z(p_3)+\gamma (p_4)+\gamma (p_5)+\gamma (p_6)\). The contribution of the real photon emission process has the form as

$$\begin{aligned} \Delta \sigma _\mathrm{real}=\frac{1}{3!}\frac{(2\pi )^4}{2s}\int \mathrm{d}\Phi _{4} \overline{\sum _\mathrm{spin}} |\mathcal{M}_\mathrm{real}|^2, \end{aligned}$$
(2.5)

where \(\frac{1}{3!}\) is due to the final three identical photons. The phase space element of the four particles is defined as

$$\begin{aligned} {\mathrm{d}\Phi _4}=\delta ^{(4)} \left( p_1+p_2-\sum _{i=3}^6p_i \right) \prod _{j=3}^6 \frac{d^3 \vec {p}_j}{(2 \pi )^3 2 E_j}. \end{aligned}$$
(2.6)

By employing the dipole subtraction method we extract the IR singularities from the real photon emission correction and combine them with the virtual contribution. In this method the IR finite real correction is obtained by subtracting an auxiliary function from the squared amplitude of the real photon emission process before integrating over phase space due to the subtraction function having the same singular structure as the squared amplitude pointwise in phase space. 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 partons [3739] and subsequently was generalized to photon radiation off charged particles with arbitrary mass by Dittmaier [40]. In our calculations, we follow the approach of Ref. [40], and we check the independence on the parameter \(\alpha \in (0,1]\), which essentially controls the region of dipole phase space, such as \(\alpha =1\) means the full dipole subtraction being considered [4143]. Then the cancelation of IR singularities is verified and the result shows that the NLO EW corrected cross section for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process is independent on the IR regulator \(m_{\gamma }\) in our calculation.

To analyze the origin of the NLO EW corrections clearly, we calculate the photonic (QED) and the genuine weak corrections separately. The QED correction includes two parts: the QED virtual correction \(\Delta \sigma _{v}^\mathrm{QED}\), which comes from the diagrams with virtual photon exchange loop and the corresponding QED parts of the counterterms, and the real photon emission correction \(\Delta \sigma _\mathrm{real}\). The rest of the virtual electroweak correction part is called the weak correction \(\Delta \sigma _{v}^{W}\). Therefore, the full NLO EW corrected cross section can be expressed as

$$\begin{aligned} \sigma _\mathrm{NLO} \!&= \! \sigma _\mathrm{LO} \!+\! \Delta \sigma _{v}\!+\!\Delta \sigma _\mathrm{real} \!=\! \sigma _\mathrm{LO} \!+\! \Delta \sigma _{v}^\mathrm{QED}\!+\!\Delta \sigma _{v}^{W}\!+\!\Delta \sigma _\mathrm{real} \nonumber \\&= \sigma _\mathrm{LO} + \Delta \sigma ^\mathrm{QED}+\Delta \sigma _{v}^{W} =\sigma _\mathrm{LO}(1+\delta ^\mathrm{QED}+\delta ^{W}) \nonumber \\&= \sigma _\mathrm{LO}(1+\delta ^\mathrm{EW}), \end{aligned}$$
(2.7)

where the \(\delta ^\mathrm{QED}, \delta ^{W}\) and \(\delta ^\mathrm{EW}\) are the pure QED, genuine weak and full EW relative corrections, respectively.

3 Numerical results and discussions

3.1 Input parameters and kinematic cuts

For the numerical evaluation we adopt the \(\alpha \)-scheme and take the following SM input parameters [44]:

$$\begin{aligned}&M_{W} \!=\! 80.398~\mathrm{GeV}, M_{Z} \!=\! 91.1876~\mathrm{GeV}, \Gamma _Z \!=\! 2.4952~\mathrm{GeV} \nonumber \\&m_e \!=\! 0.510998929~\mathrm{MeV}, m_{\mu }\!=\!105.6583715~\mathrm{MeV},\nonumber \\&m_{\tau }\!=\!1.77682~\mathrm{GeV}, \nonumber \\&m_u = 66~\mathrm{MeV}, m_c=1.275~\mathrm{GeV},m_t=173.5~\mathrm{GeV}, \nonumber \\&m_d = 66~\mathrm{MeV}, m_s=95~MeV,m_b=4.65~\mathrm{GeV}. \end{aligned}$$
(3.1)

We take the fine structure constant \(\alpha (0)=1/137.035999074\) defined in the Thomson limit. The current masses for light quarks (\(m_u\) and \(m_d\)) can reproduce the hadronic contribution to the shift in the fine structure constant \(\alpha (M_Z)\) [45]. We take the Higgs boson mass as \(M_H=126~\)GeV, and its decay width is estimated by using the HDECAY program [46]. The CKM matrix, whose matrix elements appear only in the loop contribution, is set to be the unity matrix.

We apply the Cambridge/Aachen (C/A) jet algorithm [47] to photon candidates. For a three-photon event originating from the real emission correction, if the two final photons with the smallest separation \(R\) satisfy the constraint of \(R = \sqrt{\Delta y^2 + \Delta \phi ^2} < 0.4\), where \(\Delta y\) and \(\Delta \phi \) are the differences of rapidity and azimuthal angle between the two photons, we combine this pair of photons as one new photon track and this event is called a ‘two-photon’ event including the merged photon with four-momentum \(p_{ij,\mu }=p_{i,\mu }+p_{j,\mu }\), and contrariwise, it is called a ‘three-photon’ event. In our calculation we consider only the ‘two-photon’ and ‘three-photon’ events with all the final photons satisfying the constraints as

$$\begin{aligned} p^{\gamma }_T \ge 15~\mathrm{GeV},~~~|y_{\gamma }|\le 2.5,~~~ R_{\gamma \gamma }\ge 0.4. \end{aligned}$$
(3.2)

Thereby we can exclude the inevitably infrared (IR) singularity at the tree level. We name the photon in one event with the largest photon transverse energy as the leading photon, while the photon with the next-to-largest photon transverse energy is named as the subleading photon. In the ‘inclusive’ event selection scheme we collect all the ‘two-photon’ and ‘three-photon’ events with the limitations on photons shown in Eq. (3.2). In the ‘exclusive’ event selection scheme, we include only the so-called ‘two-photon’ events satisfying the constraints as shown in Eq. (3.2). In following discussion we adopt the ‘inclusive’ scheme for event selection as default unless otherwise stated.

3.2 Total cross section

The dependence of the LO integrated cross section for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process in the SM on the colliding energy was presented in Fig.1 of Ref. [10]. When we take the same input parameters as in that reference, the coincident numerical results can be obtained. In Fig. 4a, we plot the LO, NLO EW and pure NLO QED corrected integrated cross sections as the functions of the colliding energy \(\sqrt{s}\) in the ‘inclusive’ event selection scheme, and in Fig. 4b we show the corresponding NLO EW and pure NLO QED relative corrections, \(\delta ^\mathrm{EW}\equiv \frac{\sigma _\mathrm{NLO}-\sigma _\mathrm{LO}}{\sigma _\mathrm{LO}}\) and \(\delta ^\mathrm{QED}\equiv \frac{\sigma _\mathrm{NLO}^\mathrm{QED}-\sigma _\mathrm{LO}}{\sigma _\mathrm{LO}}\). Some representative numerical results read out from Fig. 4a and b are listed in Table 1. From these figures we find all the curves for the cross sections decrease quickly with the increment of \(\sqrt{s}\), and the LO cross section is always enhanced by the NLO EW correction in the whole \(\sqrt{s}\) range plotted. When \(\sqrt{s}\) goes up from \(250\) GeV to \(1\) TeV, the NLO EW relative correction \(\delta ^\mathrm{EW}\) varies from 2.32 to \(9.61\,\%\). We also see that the pure NLO QED correction part always increases the LO cross section when \(\sqrt{s}> 270\) GeV and the pure NLO QED relative correction becomes more and more notable with the increment of \(\sqrt{s}\). In order to make a comparison of the results in different event selection schemes, we also present corresponding numerical results by adopting the ‘exclusive’ event selection scheme in Table 2. We can see that with the same \(\sqrt{s}\) the NLO EW and pure NLO QED corrected cross sections in Table 2 are less than the corresponding ones by adopting the ‘inclusive’ event selection scheme, due to the fact that all the so-called ‘three-photon’ events are abandoned in the ‘exclusive’ event selection scheme.

Fig. 4
figure 4

a The LO, NLO EW and the pure NLO QED corrected cross sections (\(\sigma _\mathrm{LO}, \sigma _\mathrm{NLO}, \sigma _\mathrm{NLO}^\mathrm{QED}\)) for the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process as the functions of the colliding energy \(\sqrt{s}\) in the ‘inclusive’ event selection scheme at the ILC. b The corresponding NLO EW and pure NLO QED relative corrections (\(\delta ^\mathrm{EW}, \delta ^\mathrm{QED}\))

Table 1 The total LO, NLO EW, pure NLO QED corrected integrated cross sections (\(\sigma _\mathrm{LO}, \sigma _\mathrm{NLO}\) and \(\sigma ^\mathrm{QED}_\mathrm{NLO}\)), and the corresponding EW and QED relative corrections (\(\delta ^\mathrm{EW}\) and \(\delta ^\mathrm{QED}\)) for the \(e^+ e^- \rightarrow Z \gamma \gamma \) process in the ‘inclusive’ event selection scheme
Table 2 The total LO cross section (\(\sigma _\mathrm{LO}\)), NLO EW, pure NLO QED corrected integrated cross sections (\(\sigma _\mathrm{NLO}\) and \(\sigma ^\mathrm{QED}_\mathrm{NLO}\)), and the corresponding EW and QED relative corrections (\(\delta ^\mathrm{EW}\) and \(\delta ^\mathrm{QED}\)) for the \(e^+ e^- \rightarrow Z \gamma \gamma \) process in the ’exclusive’ event selection scheme

3.3 Kinematic distributions

We present the LO and NLO EW corrected transverse momentum and rapidity distributions of the final \(Z\)-boson in Fig. 5a and a, respectively. The corresponding EW relative corrections \(\delta ^\mathrm{EW}\) are also plotted in Fig. 5b and b, separately. There the results are obtained by taking \(\sqrt{s} = 500\) GeV and applying the ’inclusive’ event selection scheme. From Fig. 5a,b we can see that the NLO EW correction enhances the LO differential cross section \(\mathrm{d}\sigma _\mathrm{LO}/\mathrm{d}p_T^Z\) in the low \(p_T^Z\) region. The NLO relative correction always goes down with the increment of \(p_T^Z\), and changes from being positive to negative when \(p_T^Z\) arrives at the position about \(145\) GeV. In Fig. 6a,b, we find that the LO rapidity distribution is strengthened obviously by the NLO EW correction in the central rapidity region of \(Z\)-boson at the ILC, while weakened by the quantum correction in the regions of \(|y_Z|\ge 1.4\).

Fig. 5
figure 5

a The LO and NLO EW corrected transverse momentum distributions of \(Z\)-boson with \(\sqrt{s} =500\) GeV in the ‘inclusive’ event selection scheme. b The corresponding NLO EW relative corrections

Fig. 6
figure 6

a The LO and NLO EW corrected rapidity distributions of \(Z\)-boson with \(\sqrt{s} =500\) GeV in the ‘inclusive’ event selection scheme. b The corresponding NLO EW relative corrections

The transverse momentum distributions of the leading photon (labeled by \(\gamma _1\)) and the subleading photon (labeled by \(\gamma _2\)) are plotted in Fig. 7a and b, respectively. The rapidity distributions of the leading and subleading photons are presented in Fig. 8a and b, separately. In these four figures we adopt the ’inclusive’ event selection scheme, and we take \(\sqrt{s} =500\) GeV, the cuts on photons being declared in Eq. (3.2). It can be seen from both Fig. 7a and b that the LO \(p_T\) distributions for the leading and subleading photons are enhanced in the lower \(p_T\) region (i.e., \(p_T^{\gamma _1} < 145\) GeV and \(p_T^{\gamma _2} < 75\) GeV, separately), but they are suppressed in the rest of \(p_T\) regions by the NLO corrections. The LO and NLO EW corrected transverse momentum distributions for the leading photon reach their maxima at the position about \(50\) GeV. While for the subleading photon, the LO and NLO EW corrected transverse momentum distributions always decrease with the increment of \(p_T\). Figure 8a,b show that the rapidity distributions of the leading and subleading photons are both reinforced by the NLO EW corrections in the whole plotted rapidity region. From the figures we see that both the LO and NLO corrected rapidity distributions for the leading photon have two peaks, which are located at the positions of \(|y|\sim 1\), in contrast the subleading photon rapidity distributions reach their maxima in the central rapidity region.

Fig. 7
figure 7

The LO and NLO EW corrected transverse momentum distributions of the final photons with \(\sqrt{s} =500\) GeV in the ‘inclusive’ event selection scheme. a For the leading photon. b For the subleading photon

Fig. 8
figure 8

The LO and NLO EW corrected rapidity distributions of the final photons with \(\sqrt{s} =500\) GeV in the ‘inclusive’ event selection scheme. a For the leading photon. b For the subleading photon

The LO and NLO EW corrected distributions of the separation \(R_{\gamma \gamma }\) between the final leading and subleading photons are plotted in Fig. 9a. It shows that at both the LO and the NLO the preferred kinematical configuration of the leading and subleading photons has wide separation in the rapidity–azimuthal-angle plane, and the LO and NLO \(R_{\gamma \gamma }\) distributions reach their maxima at the location of \(R_{\gamma \gamma } \sim 3\). In Fig. 9b, we depict the LO and NLO EW corrected distributions of the invariant mass of the leading and subleading photons (denoted as \(M_{\gamma \gamma }\)). It demonstrates that both the LO and NLO EW corrected \(M_{\gamma \gamma }\) distributions reach their maxima in the vicinity of \(M_{\gamma \gamma }\sim 100\) GeV, and the NLO EW correction enhances the LO differential cross section in the range of \(M_{\gamma \gamma }\le 270\) GeV.

Fig. 9
figure 9

a The LO and NLO EW corrected distributions of the separation \(R_{\gamma \gamma }\) between the final leading and subleading photons. b The LO and NLO EW corrected invariant mass \(M_{\gamma \gamma }\) distributions

From Figs. 5, 6, 7, 8, 9, we can see that the phase space dependence of the NLO EW correction is nontrivial and sizable, and the NLO EW correction does not observably change the LO distribution line shape in the case of taking the ’inclusive’ event selection scheme.

As we know that one of the most important reactions at the ILC for Higgs boson precision study is the \(e^+e^- \rightarrow ZH\) process followed by \(H \rightarrow \gamma \gamma \) decay, while this signal process is accompanied by a serious background process \(e^+e^- \rightarrow Z \gamma \gamma ~\). The one-loop radiative corrections to this signal process \(e^+e^- \rightarrow ZH\) within the SM were calculated by Denner et al. [48]. Here we follow the strategy used in Ref. [48] for the calculation of the \(e^+e^- \rightarrow ZH\) process up to the EW NLO within the SM, and we adopt the input parameters presented in our work (see Section III.1) to calculate the LO and NLO EW corrected results for the \(e^+e^- \rightarrow ZH \rightarrow Z\gamma \gamma \) signal process. The decay width of the SM Higgs is obtained by using the program HDECAY [46]. Since the kinematics of the signal events is distinctively different from that of background events. This difference can be used to suppress the background and enhance the ratio of signal to background (S/B). Taking advantage of the kinematic difference, we expect that we can impose the optimal cuts to extract the signal \(e^+e^-\rightarrow ZH\rightarrow Z\gamma \gamma \) from the SM background \(e^+e^- \rightarrow Z \gamma \gamma ~\)efficiently. To illustrate the distribution differences between the signal and the background, we present the normalized LO and NLO EW corrected distributions of various kinematic observables of the final particles for the signal process \(e^+e^-\rightarrow ZH\rightarrow Z\gamma \gamma \) and the background process \(e^+e^- \rightarrow Z \gamma \gamma ~\)in Fig. 10a–f. All the results are presented in conditions of \(\sqrt{s} = 500~\mathrm{GeV}, p^{\gamma }_T \ge 15~\mathrm{GeV}, |y_{\gamma }|\le 2.5, R_{\gamma \gamma }\ge 0.4\), and the ’inclusive’ event selection scheme. We show the \(p_T\) distributions of the final \(Z\)-boson, leading photon and subleading photon in Fig. 10a–c, separately. We see that, compared to the background, the typical feature of the signal is that the final state particles are distributed in the large transverse momentum regions, especially for the final \(Z\)-boson and the leading photon. The normalized rapidity distributions of the final leading and subleading photons are demonstrated in Fig. 10d and e, respectively. It shows that the leading and subleading photons from Higgs boson decay, mainly appear in the central rapidity region, while the corresponding distributions for the background process are rather flat. We can see also from Fig. 10f that the leading and subleading photons produced in the background are more dramatically separated than in the signal process \(e^+e^- \rightarrow ZH \rightarrow Z\gamma \gamma \). From all of Fig. 10a–f, we can conclude that if we take some proper cuts on kinematic variables of the final \(Z\)-boson and photons, the background from the \(e^+e^- \rightarrow Z \gamma \gamma ~\)process can be significantly suppressed.

Fig. 10
figure 10

The normalized kinematic distributions for the signal process \(e^+e^- \rightarrow ZH \rightarrow Z\gamma \gamma \) and the background \(e^+e^- \rightarrow Z \gamma \gamma ~\)process at ILC at \(\sqrt{s}=500\) GeV. All curves in the six figures are normalized by their total cross sections. a Transverse momentum distributions of the final \(Z\) boson. b Transverse momentum distributions of the leading photon. c Transverse momentum distributions of the subleading photon. d Rapidity distributions of the leading photon. e Rapidity distributions of the subleading photon. f Distributions of the separation \(R_{\gamma \gamma }\) between the final leading and subleading photons

4 Summary

The \(e^+e^- \rightarrow Z \gamma \gamma ~\) process is very important for understanding the nature of the Higgs boson and searching for new physics beyond the SM. In this work we report on our calculation of the full NLO EW contributions to the \(e^+e^- \rightarrow Z \gamma \gamma \) process in the SM, and we analyze the EW quantum effects on the total cross section and the kinematic distributions of the final particles. We study the dependence of the \(Z\gamma \gamma \) production rate on the event selection scheme and provide distributions of some important observables. We find that the full NLO EW corrections can enhance the LO total cross sections quantitatively from 2.32 to \(9.61\,\%\) when colliding energy goes up from \(250\) GeV to \(1\) TeV, and the size of the NLO correction exhibits a strong dependence on the observable and on phase space. We conclude that in studying the signal process \(e^+e^- \rightarrow ZH \rightarrow Z \gamma \gamma ~\), the background events of \(e^+e^- \rightarrow Z \gamma \gamma ~\) process can be suppressed significantly if we take appropriate kinematic cuts on the final products.