Next-to-next-to-leading order QCD ${\otimes }$ EW corrections to Z-boson pair production at electron-positron colliders

The discovery of the Higgs boson [1, 2] at the CERN Large Hadron Collider (LHC) represents a landmark validation of the Standard Model (SM). Nevertheless, achieving a thorough and precise understanding of the SM is essential for unraveling the mechanism of electroweak (EW) symmetry breaking and advancing our search for new physics beyond the SM. In this pursuit, the study of neutral gauge boson pair production assumes a pivotal role. It not only sheds light on the study of the EW gauge structure and exploration of triple gauge couplings (TGCs) but also serves as a crucial avenue for uncovering traces of new physics. Additionally, Z-boson pairs production is an irreducible background to the production of a Higgs boson decaying into two neutral vector bosons. The production of a Z-boson pair has been comprehensively examined in $ e^+e^- $ collisions at the LEP [3−8], in $ p\bar{p} $ collisions at the Tevatron [9−11], and in $ pp $ collisions at the LHC [12−19]. Proposals for several future lepton colliders, including the International Linear Collider (ILC) [20−22], the Circular Electron-Positron Collider (CEPC) [23, 24] and the Future Circular Collider (FCC-ee) [25, 26], hold immense potential for advancing our comprehension of the EW interactions within the SM and probing potential candidates for new physics beyond the SM. These lepton colliders, benefiting from cleaner experimental environments when compared to hadron colliders, are anticipated to achieve measurements with permille precision [21, 24, 26]. Therefore, it is imperative that the theoretical predictions match or exceed the level of precision of hadron colliders. Extensive research has been devoted to studying the leading order (LO) predictions and next-to-leading order (NLO) corrections for the $ e^+e^- \rightarrow ZZ $ process over the past few decades [27−34]. However, the next-to-next-to-leading order (NNLO) corrections to $ e^+e^- \rightarrow ZZ $ still remain elusive. These NNLO corrections comprise two components: the pure EW $ {\cal{O}}(\alpha^2) $ corrections and mixed $ \text{QCD} \otimes \text{EW} $ $ {\cal{O}}(\alpha \alpha_s) $ corrections. The calculation of the pure EW corrections poses a significant challenge due to the tremendous number of two-loop Feynman diagrams with multiple scales. Conversely, the mixed $ \text{QCD} \otimes \text{EW} $ corrections, though intricate, are more tractable, and may hold greater numerical significance due to the relatively large QCD coupling constant. In light of these considerations, in this study, we focus on the analytic calculation of NNLO $ \text{QCD} \otimes \text{EW} $ corrections to $ e^+e^- \rightarrow ZZ $.

All mixed $ \text{QCD} \otimes \text{EW} $ two-loop Feynman diagrams for $ e^+e^- \rightarrow ZZ $ are factorizable and originate from the corrections to TGCs. After implementing the tensor reduction procedure, the two-loop scattering amplitude is expressed as a linear combination of scalar Feynman integrals, which can be categorized into several families. Within each family, the scalar integrals are interrelated and can be systematically reduced to a set of master integrals (MIs) using integration-by-parts (IBP) identities [35, 36]. We conduct analytic calculations of these two-loop MIs by utilizing the canonical differential equation method [37, 38]. Leveraging leading singularities [37−41], dlog integrals [42, 43], Magnus exponential [41, 44−46], and various other techniques, we can establish an appropriate basis comprising a set of uniform transcendental (UT) integrals that satisfy the canonical differential equations. Furthermore, we present Chen's iterated integral [47] solution of the canonical differential system in terms of Goncharov polylogarithms (GPLs) [48], following the rationalization of all square roots appearing in the system. Ultimately, we apply our analytic expressions of the canonical MIs to the computation of the NNLO $ \text{QCD} \otimes \text{EW} $ corrections to $ e^+e^- \rightarrow ZZ $ and provide detailed numerical predictions for integrated and differential cross sections.

Before delving into further analyses, it is essential to provide an overview of the current state-of-the-art in analytic studies of two-loop Feynman integrals, stemming from the NNLO $ \text{QCD} \otimes \text{EW} $ corrections. The two-loop two-point MIs involved in the $ \text{QCD} \otimes \text{EW} $ two-loop vector-boson self-energies have been thoroughly examined in Refs. [49−54]. The analytic results of the two-loop triangle MIs for the QCD corrections to $ H \rightarrow V^{\ast}V^{\ast} $ ($V=W,\; Z$) with massless internal particles and arbitrary external momenta have been effectively expressed in terms of GPLs in Refs. [55−57]. These results can be employed in the calculation of the massless two-loop three-point MIs involved in $ {\cal{O}}(\alpha\alpha_s) $corrections to the $ V^{\ast}ZZ $ vertex. Additionally, all two-loop three-point MIs necessary for the QCD corrections to $ H \rightarrow \gamma \gamma $ and $ H \rightarrow Z \gamma $ with a massive top-quark loop have been expressed in terms of harmonic polylogarithms [58−60] and GPLs [61, 62], respectively. Moreover, all two-loop triangle MIs for QCD corrections to the interaction vertex of an off-shell neutral boson with a pair of W bosons are well documented in terms of GPLs [63−65]. A canonical set of two-loop triangle MIs, contributing to $ {\cal{O}}(\alpha\alpha_s) $ corrections to the $ HZV^{\ast} $ ($V=Z, \;\gamma$) vertex and $ H \rightarrow ZZ^{\ast} $ decay, has been examined in Refs. [66, 67]. The MIs for $ H \rightarrow ZZ^{\ast} $, characterized by a massive internal loop and three distinct external momentum squares, are expressed in terms of Chen's iterated integrals with dlog one-forms that contain a residual non-rationalizable square root [67]. The two-loop triangle MIs with two on-shell Z-boson legs and a massive internal loop, arising from the NNLO $ \text{QCD} \otimes \text{EW} $ corrections to $ e^+e^- \rightarrow ZZ $, can be considered a special case of the two-loop MIs in the $ H \rightarrow ZZ^{\ast} $ decay, where the two Z-boson external legs are forced on their mass-shell. We analytically calculate these two-loop three-point MIs associated with two on-shell Z-boson legs and present them in a more user-friendly format, namely GPLs. This representation offers a faster, more accurate, and more convenient numerical evaluation when compared to the iterated integrals with dlog one-forms.

The rest of this paper is structured as follows. In Sec. II, we introduce our notations and conventions, present the framework for our calculation, delineate the workflow of our computation, and engage a discussion concerning $ \gamma^5 $schemes. In Sec. III, we dedicate to the elaboration on the analytic calculation for the two-loop three-point MIs originating from $ {\cal{O}}(\alpha \alpha_s) $ corrections to $ e^+e^- \rightarrow ZZ $. In Sec. IV, we utilize analytic MIs derived in Sec.III to compute NNLO $ \text{QCD} \otimes \text{EW} $ corrected integrated and differential cross sections for Z-boson pair production at electron-positron colliders. Finally, we conclude with a concise summary in Sec. V.

In this study, we examine NNLO $ \text{QCD} \otimes \text{EW} $ corrections to the scattering process

All the external particles are on their mass-shell, i.e., $ p_1^2 = p_2^2 = m_e^2 $ and $ p_3^2 = p_4^2 = m_{ Z}^2 $. The Mandelstam invariants for this process are defined by

which satisfy $ s + t + u = 2\, (m_e^2 + m_{ Z}^2) $ due to energy-momentum conservation. The four-momentum of the final-state Z boson can be parameterized as

where β denotes the magnitude of momentum normalized by $ \sqrt{s}/2 $, θ denotes the polar angle with respect to the direction of the incident positron beam, and ϕ denotes the azimuthal angle around the positron beam axis. In the center-of-mass frame, β can be denoted as

Leptons and light quarks $ (u, d, c, s) $ are treated as massless particles throughout our calculation unless explicitly stated otherwise. Then, the Mandelstam variables t and u can be expressed as

The lowest-order unpolarized differential cross section for this $ 2 \rightarrow 2 $ scattering process is simply provided by

where $ {\cal{M}}_{0} $ denotes the tree-level Feynman amplitude. Throughout our calculation, we adopt the 't Hooft-Feynman gauge.

The on-shell (OS) renormalization scheme [68, 69] is employed to remove the ultraviolet (UV) divergences, which are regularized by using dimensional regularization (DR) [70, 71] in both one-loop and two-loop calculations. An infinitesimal fictitious photon mass is introduced to regularize the infrared (IR) divergences induced by the virtual photon in loops, which can be canceled exactly by those from the real photon emission. We adopt the phase space slicing method [72−74] to extract the IR singularities from the real photon emission. Additionally, the initial-state QED corrections are implemented in the leading-logarithmic approximation with the structure-function method [73, 75]. Furthermore, the $ G_{\mu} $ scheme is adopted for the renormalization of electric charge, wherein the fine structure constant can be derived from the Fermi constant $ G_{\mu} $ via

It is important to note that the corresponding charge renormalization constant in the $ G_{\mu} $ scheme should be modified to $ \delta Z_{e, {G_\mu}} = \delta Z_{e,{\alpha(0)}} - \Delta r/2 $, where the subtraction term $ \Delta r $ comprises the higher-order corrections to muon decay [69, 76]. The explicit expression for the renormalization constant in the $ \alpha(0) $ scheme, $ \delta Z_{e,{\alpha(0)}} $, has been detailed discussed in Refs. [68, 69]. For the computation of the LO cross section and NLO EW corrections, we employ the modified $ \mathrm{FormCalc}$ and $\mathrm{LoopTools} $ packages [77, 78]. To validate the accuracy and reliability of our NLO results, we use our developed computational tools to conduct rigorous numerical cross-checks against the benchmarks provided in Ref. [33], confirming excellent concordance. For a comprehensive elucidation of the NLO EW corrections to $ e^+e^- \rightarrow ZZ $, please refer to Refs. [31, 33, 34].

It is well known that the cyclicity of Dirac trace conflicts with the anticommutativity relation $ \{\gamma^5,\, \gamma^{\mu}\} = 0 $ in the framework of DR. Consequently, two prominent strategies exist for the practical treatment of $ \gamma^5 $ in DR: either preserving or disregarding the anticommutativity relation. To circumvent the challenging nontrivial renormalizations and computational complexities associated with the 't Hooft-Veltman scheme [70] and its subsequent iterations [79−83], wherein the anticommutativity relation is sacrificed, we adopt the Körner-Kreimer-Schilcher (KKS) scheme [84−86] and its reformulated version [87, 88]. Within the KKS scheme, the cyclicity of the Dirac trace with odd numbers of $ \gamma^5 $ is deliberately abandoned, prioritizing the preservation of the anticommutativity property of $ \gamma^5 $. To address the ambiguity arising from the dependence of the Dirac trace on the reading start point in a fermion chain, a specific reading prescription is introduced. In this prescription, the final expression of the non-cyclic $ \gamma^5 $ trace is defined as the average of all possible reading points, starting from the head and tail of an axial-vector three-point fermion-gauge vertex that contains the maximal one-particle-irreducible non-singlet-type loop correction [86, 88]. Detailed technical descriptions for the KKS scheme and its reformulated version can be found in Refs. [86, 88]. For additional discussions on $ \gamma^5 $ schemes in $ d = 4- 2\, \epsilon $ dimensions, we refer readers to Ref. [89], which provides deeper insights and comprehensive considerations.

The NNLO $ \text{QCD} \otimes \text{EW} $ $ {\cal{O}}(\alpha \alpha_s) $ corrections to $ e^+ e^- \rightarrow ZZ $ are categorized into the following three distinct contributions:

1. Two-loop vertex corrections induced by attaching one gluon to each one-loop quark line in every possible way.

2. One-loop vertex corrections with insertions of $ {\cal{O}}(\alpha_s) $ quark mass counterterm ¹.

3. $ {\cal{O}}(\alpha\alpha_s) $ vertex counterterms.

Some representative Feynman diagrams for these $ {\cal{O}}(\alpha\alpha_s) $ corrections are depicted in Fig. 1. Furthermore, the $ {\cal{O}}(\alpha_s) $ quark mass renormalization constant $ \delta m_q $ in the OS scheme is provided by [90]

Figure 1.  (color online) Representative Feynman diagrams for NNLO $ \text{QCD} \otimes \text{EW} $ corrections to $ e^+ e^- \rightarrow ZZ $. The black crosses symbolize the heavy quark mass counterterm at $ {\cal{O}}(\alpha_s) $, and the red circled crosses signify the vertex counterterm at $ {\cal{O}}(\alpha\alpha_s) $.

where $ C(\epsilon) = (4\pi)^{\epsilon}\, \Gamma(1+\epsilon) $ and $ C_F = 4/3 $. The explicit expressions for the $ {\cal{O}}(\alpha \alpha_s) $ renormalization constants can be derived by substituting the one-loop gauge-boson self-energies in the NLO EW renormalization constants [68] with their corresponding two-loop $ {\cal{O}}(\alpha\alpha_s) $ counterparts [91].

We employ a comprehensive set of tools and methodologies to calculate the NNLO $ \text{QCD} \otimes \text{EW} $ corrections to $ e^+e^- \rightarrow ZZ $. Initially, $\mathrm{FeynArts} $ [92] is used to generate the Feynman diagrams and amplitudes for the Z-pair production process. Then, we utilize an in-house $\mathrm{Mathematica} $ code based on $ \mathrm{FeynCalc}$ [93, 94] to manipulate the Dirac algebra and express the two-loop amplitudes as linear combinations of numerous scalar Feynman integrals, which are not independent. To reduce these extensive Feynman integrals to a small set of irreducible MIs, we use the publicly available package $\mathrm{Kira} $ [95, 96], which implements the Laporta algorithm [97] for solving IBP identities. Subsequently, we analytically calculate all MIs appearing in the NNLO $ \text{QCD} \otimes \text{EW} $ corrections using the method of canonical differential equations and express them in terms of GPLs. Additionally, several $\mathrm{Mathematica} $ ancillary files containing the analytic expressions of the MIs are provided as supplementary materials accompanying the arXiv submission of this paper. To validate our analytical results, numerical checks for the MIs are conducted using the $\mathrm{AMFlow} $ [98, 99] package, showing perfect agreement between our analytical and numerical results. Further details on the analytic calculation of MIs is presented in Sec.III.

In this section, we concentrate on the analytic calculation of the massive two-loop triangle MIs for the NNLO $ \text{QCD} \otimes \text{EW} $ corrections to the $ e^+ e^- \rightarrow ZZ $ process.

The integral family $ {\cal{F}} $ of the two-loop three-point MIs with two on-shell Z boson legs and a massive internal quark loop, as depicted in the panel (a) of Fig. 1, is defined by

where $ d=4-2\,\epsilon $ is the spacetime dimension, and the integration measure is defined as

The integral family $ {\cal{F}} $ is identified by the following set of propagators:

with $ p_3^2=p_4^2=m_{ Z}^2 $. The three topologies illustrated in Fig. 2 correspond to the following three top-sectors of the integral family $ {\cal{F}} $:

Figure 2.  (color online) Three top-level topologies of the integral family $ {\cal{F}} $. The red lines represent massive propagators, while the thin black lines denote massless propagators. The purple lines indicate an external off-shell leg with momentum squared s, whereas the black external lines signify the on-shell Z bosons.

In this section, we elaborate on the construction of the canonical basis for the integral set $ {\cal{S}} $ induced by the three top-sectors,

where the subscript $ \digamma $ is defined by

For brevity and clarity, the three top-sectors are denoted as $ {\cal{I}}_{134567} $, $ {\cal{I}}_{124567} $ , and $ {\cal{I}}_{123467} $, with the sub-sectors of these top-sectors adhering to the same naming convention. Following the IBP reduction, we derive a set of $ 26 $ MIs for the integral set $ {\cal{S}} $. In Sec. III.B, we focus on the analytic calculation of these MIs using the method of canonical differential equations.

First, we establish a system of standard differential equations,

for a pre-canonical basis ${\mathbf{I}}(\boldsymbol{x},\epsilon)$ that comprises $ 26 $ MIs derived using the Laporta algorithm. This differential system is generated by taking partial derivatives of the MIs with respect to the kinematic variables ($ s, m_{ Z}^2, m_t^2 $) and employing IBP identities for scalar reduction, facilitated by $\mathrm{LiteRed} $ [100, 101] and $\mathrm{Kira} $. As suggested in Refs. [37, 38], we can transform the pre-canonical basis ${\mathbf{I}}(\boldsymbol{x},\epsilon)$ into a canonical basis ${\mathbf{g}}(\boldsymbol{x},\epsilon)$, satisfying the following canonical differential equations:

where $ \mathbb{M}_i $ are constant matrices, and the symbol letters $ \eta_i $ are algebraic functions of $\boldsymbol{x}$. The general solution to the canonical differential equations (16) can be expressed in terms of Chen's iterated integrals [47], formulated as the following path-ordered exponential:

where $ {\cal{P}} $ stands for the path ordering along the integration path γ from the start point $\boldsymbol{x}_0$ to the point of interest $\boldsymbol{x}$, and $\mathbf{g}(\boldsymbol{x}_0,\epsilon)$ is the vector of boundary constants at $\boldsymbol{x}_0$. The path-ordered exponential is a shorthand notation for the following series:

The n-th term in this expansion represents an n-fold iterated integral,

where $\mathbb{K}(t)\, \mathrm{d}t$ is the pullback of the matrix-valued differential one-form $\mathrm{d} \mathbb{A}$ to the unit interval $ [0, 1] $. It is crucial to note that these iterated integrals do not depend on the choice of the integration path γ, unless it encounters the singularities or branch cuts of $\mathrm{d}\mathbb{A}$.

We adopt a "bottom-up" sector-wise methodology to construct the canonical differential equations by arranging all MIs sequentially from lower to higher sectors and then systematically transforming them into the canonical basis, sector by sector. We employ Baikov leading singularity analysis techniques [39−41] to construct the UT integrals or candidate UT integrals for a given sector, aided by the publicly available package $ \mathrm{DlogBasis}$ [43]. Furthermore, if required, the contributions arising from sub-sectors are accounted for to transform the off-diagonal blocks into canonical forms. Along with the application of the Magnus exponential method [44−46], we derive the following canonical basis:

where $ \lambda_{1,2,3} $ are three square roots related to kinematic variables,

The canonical basis $\mathbf{g}(\boldsymbol{x},\epsilon)$ defined in Eq. (20) is dimensionless and can be expressed as a Taylor series in $ \epsilon $. For the two-loop three-point integral family with three distinct off-shell legs and a massive internal loop, four square roots that cannot be simultaneously rationalized emerge [66, 67]. However, in our specific case, these four square roots degenerate into the three square roots defined in Eq.(21). Advantageously, the square roots $ \lambda_{1,2,3} $ stemming from the NNLO $ \text{QCD} \otimes \text{EW} $ corrections to $ e^+e^- \rightarrow ZZ $ can be simultaneously rationalized as

by the following change of variables:

Consequently, the canonical differential equations for $ \mathbf{g}(x,y,\epsilon) $ can be cast into the following $ d \log $ form,

where the $ 18 $ rational symbol letters are given as

and the explicit expressions of the constant matrices $ \mathbb{M}_i $ are available in the supplementary file "dlog-fom_Matrix.m." In the Euclidean region delineated by

all letters are real and positive, and thus, all MIs are real functions of the dimensionless variables x and y. Furthermore, it is evident that the differential one-form $ d\mathbb{A} $ is composed exclusively of rational functions, as all symbol letters in Eq. (25) are rational functions of kinematic variables. Therefore, the iterated integral solution (18) of this canonical differential system can be expressed in terms of GPLs. The GPLs are recursively defined by [48]

with

where $(a_1,\;a_2,\cdots, a_n)$ is referred to as the weight vector of $G(a_1,\;a_2, \cdots, a_n; z)$, and the length of the weight vector is called the weight of the GPL.

Within the framework of the differential system constructed using the canonical basis as outlined in Eq. (20), we determine the integration constants by an independent and simpler calculation of the canonical basis at the boundary point $\boldsymbol{x}_0 = (0, \,1)$, which corresponds to the vanishing external momenta $ s = m_{ Z}^2 = 0 $. At this specific kinematic point, all MIs in the pre-canonical basis ${\mathbf{I}}(\boldsymbol{x},\epsilon)$ degenerate into equal-mass vacuum integrals. It is noteworthy that all canonical MIs $g_i\, (i=1,\;2,\cdots,26)$ are regular at the boundary point $\boldsymbol{x}_0$. Through a straightforward analysis of the canonical transformation in Eq. (20), we conclude that all canonical MIs vanish at $\boldsymbol{x}_0$ except for $ g_1 $. Under the normalization of the integration measure specified in Eq. (10), the boundary conditions are formulated as

Leveraging these boundary conditions, we are able to derive the analytic solution of the canonical differential system. All canonical MIs are expressed as Taylor series in $ \epsilon $, with coefficients articulated in terms of GPLs. The analytic expressions of these canonical MIs, $g_i\; (i=1, 2,\cdots,$ 26), are presented in Appendix A up to the order of $ \epsilon^2 $. For further convenience and accessibility, the complete analytic expressions for all MIs up to $ {\cal{O}}(\epsilon^4) $ are available in the supplementary file "analytic_MIs.m," which accompanies the arXiv submission of this paper.

We utilize the $\mathrm{Mathematica} $ package $\mathrm{PolyLogTools} $ [102−104] and the $\mathrm{C++} $ library $\mathrm{GiNaC} $ [105, 106] for the symbolic computation and numerical evaluation of GPLs. To ensure the accuracy and reliability of our analytic results, we also perform numerical checks on the $ 26 $ MIs defined by Eq. (20) in comparison with the numerical results from the publicly available package $\mathrm{AMFlow} $ [98, 99]. We find excellent agreement between our analytical and numerical results, with an accuracy exceeding $ 10^{-99} $, i.e.,

where $ g_i^{{ \mathrm{(GPL)}}} $ and $ g_i^{{ \mathrm{(AMF)}}} $ denote the numerical results obtained from our analytic expressions and those obtained using the $\mathrm{AMFlow} $ package, respectively. In Table 1, we present the numerical values for two representative MIs, $ g_{23} $ and $ g_{26} $, evaluated at the Euclidean point $ (x, y)= $ (0.5, 0.2) up to $ {\cal{O}}(\epsilon^4) $, accurate to $ 90 $ decimal places. The complete numerical results for all $g_i\; (i = 1, \;2,\cdots, 26)$ at this validation point with 100-digit precision are provided in supplementary files "numMIs_GPL.m" and "numMIs_AMF.m," respectively.

MI	Order	Numerical value
$ g_{23} $	$ \epsilon^0 $	$ 0 $
	$ \epsilon^1 $	$ 0 $
	$ \epsilon^2 $	$ -\,7.26101564334485215986593566256166202846742579354269938897134626597079498884306857349066237 $
	$ \epsilon^3 $	$ \,11.86463522232032158447131561240614580110409186179734960547391102577769656727839159997486514 $
	$ \epsilon^4 $	$ -\,9.28510545026443974881656592262564288473103903999133306919325410029473865878053069212244617 $
$ g_{26} $	$ \epsilon^0 $	$ 0 $
	$ \epsilon^1 $	$ 0 $
	$ \epsilon^2 $	$ \quad7.26101564334485215986593566256166202846742579354269938897134626597079498884306857349066237 $
	$ \epsilon^3 $	$ -\,9.68816618676192114865212264256976565734577230132356179420044708343871148651833620502830076 $
	$ \epsilon^4 $	$ \quad3.91185381115700898808123227066867633509446778857243575613151816779834900833108054340221357 $

Table 1.  Numerical values for two representative MIs from sectors $ {\cal{I}}_{13467} $ and $ {\cal{I}}_{12457} $, evaluated at $ (x, y) = (0.5, 0.2) $ up to $ {\cal{O}}(\epsilon^4) $, with a precision of $ 90 $ decimal places.

In this section, we calculate the integrated cross section and various kinematic distributions for $ e^+ e^- \rightarrow ZZ $ using the analytic expressions of MIs derived in Sec.III, and we present a detailed phenomenological analysis based on our calculations. We neglect the masses of all fermions except for the top and bottom quarks. All relevant SM input parameters are taken as follows [107]:

The renormalization scale for the strong coupling constant is selected as $ \mu = m_{ Z} $. For comparison purposes, the numerical computations are performed in the $ \alpha(0) $ and $ G_{\mu} $ schemes.

To facilitate the presentation of the integrated cross section for $ e^+ e^- \rightarrow ZZ $, incorporating NLO EW corrections and NNLO mixed $ \text{QCD} \otimes \text{EW} $ corrections, we express the cross section as

where $ \delta_{ \mathrm{EW}} $ and $ \delta_{ \mathrm{QCD \otimes EW}} $ denote the relative corrections for NLO EW and NNLO mixed $ \text{QCD} \otimes \text{EW} $ contributions, respectively, and they are defined as follows:

The LO and NNLO corrected integrated cross sections for the process $ e^+ e^- \rightarrow ZZ $, as functions of $ e^+e^- $ colliding energy $ \sqrt{s} $, are depicted in Fig. 3 for the $ \alpha(0) $ scheme (left) and $ G_{\mu} $ scheme (right). The corresponding EW and $ \text{QCD} \otimes \text{EW} $ relative corrections are presented in the lower panels. As illustrated in this figure, the production cross sections exhibit similar dependencies on $ \sqrt{s} $ in both schemes. The LO integrated cross sections are sensitive to the colliding energy, increasing sharply near the Z-pair production threshold as $ \sqrt{s} $ increases, peaking around $ \sqrt{s}\sim 207\; \mathrm{GeV} $, and then decreasing smoothly at higher energies. As shown in the lower panels of Fig. 3, the EW correction significantly suppresses the LO cross section near the threshold region, which can be attributed to the Coulomb singularity effect [68, 75, 108, 109]. The NLO EW relative correction increases from approximately −45% to around 2.7% in the $ \alpha(0) $ scheme and to approximately 6.3% in the $ G_{\mu} $ scheme as $ \sqrt{s} $ increases from the threshold to $ 1000\; \mathrm{GeV} $. However, the $ \text{QCD} \otimes \text{EW} $ correction enhances the LO cross section across the entire plotted $ \sqrt{s} $ region. The relative correction exceeds 1.1% in the $ \alpha(0) $ scheme and amounts to approximately 0.3% in the $ G_{\mu} $ scheme. Compared to the $ \text{QCD} \otimes \text{EW} $ correction in the $ \alpha(0) $scheme, the $ \text{QCD} \otimes \text{EW} $ correction in the $ G_{\mu} $ scheme is reduced by approximately an order of magnitude. This effect is attributed to the incorporation of certain significant higher-order corrections into the LO cross section [68, 69, 75]. In the $ \text{QCD} \otimes \text{EW} $ relative correction, a resonance peak appears at $ \sqrt{s}=2\, m_{t} $ due to the resonance effect caused by top-quark loop integrals at the pseudo-threshold of top-quark pair production. To present the numerical results more precisely, we summarize in Table 2 the LO, NLO, and NNLO integrated cross sections, as well as the corresponding EW and $ \text{QCD} \otimes \text{EW} $ relative corrections, for some representative colliding energies in the $ \alpha(0) $ and $ G_{\mu} $ schemes.

Figure 3.  (color online) LO and NNLO corrected integrated cross sections, along with the corresponding EW and $ \text{QCD} \otimes \text{EW} $ relative corrections, for $ e^+e^- \rightarrow ZZ $as functions of $ e^+e^- $colliding energy in the $ \alpha(0) $ (left) and $ G_{\mu} $ (right) schemes.

$\sqrt{s} \text{/GeV}$	Scheme	$\sigma_{ {\text{LO} } } \text{/pb}$	$\sigma_{ {\text{NLO} } } /\text{pb}$	$ \delta_{ {\text{EW}}} $(%)	$\sigma_{ {\text{NNLO} } } \text{/pb}$	$ \delta_{ {\text{QCD} \otimes \text{EW}}} $(%)
240	$ \alpha(0) $	1.077615	0.94682	-12.138	0.95924	1.1529
	$ G_{\mu} $	1.158461	1.05435	-8.987	1.05759	0.2792
250	$ \alpha(0) $	1.019355	0.90740	-10.983	0.91915	1.1528
	$ G_{\mu} $	1.095830	1.01048	-7.789	1.01354	0.2791
350	$ \alpha(0) $	0.632563	0.59909	-5.291	0.60635	1.1477
	$ G_{\mu} $	0.680020	0.66712	-1.897	0.66898	0.2738
500	$ \alpha(0) $	0.385130	0.37745	-1.994	0.38191	1.1581
	$ G_{\mu} $	0.414023	0.42026	1.507	0.42144	0.2848
1000	$ \alpha(0) $	0.140378	0.14419	2.721	0.14582	1.1562
	$ G_{\mu} $	0.150909	0.16051	6.368	0.16094	0.2830

Table 2.  LO, NLO, and NNLO corrected integrated cross sections, as well as the corresponding EW and $ \text{QCD} \otimes \text{EW} $ relative corrections, for $ e^+e^- \rightarrow ZZ $ at some representative colliding energies in $ \alpha(0) $ and $ G_{\mu} $ schemes.

In this subsection, we analyze the kinematic distributions of the final-state Z boson for $ e^+e^- \rightarrow ZZ $ at $ \sqrt{s} = 240\; \mathrm{GeV} $. The EW and $ \text{QCD} \otimes \text{EW} $ differential relative corrections with respect to the kinematic variable x are defined as

Given the indistinguishability of the two Z bosons in the final state, the kinematic distributions of the final-state Z boson are defined as the average of the distributions for the two identical Z bosons.

Due to $ {\cal{CP}} $ conservation and Bose symmetry, the scattering angle distribution of the final-state Z boson exhibits forward-backward symmetry, i.e., ${\mathrm{d}\sigma}/{\mathrm{d}\cos\theta}$ is symmetric with respect to $ \cos\theta = 0 $. In Fig. 4, we illustrate LO and NNLO corrected scattering angle distributions of final-state Z boson, as well as the corresponding EW and $ \text{QCD} \otimes \text{EW} $ relative corrections, for $ e^+e^- \rightarrow ZZ $ at $ \sqrt{s} = 240\; \mathrm{GeV} $ in both $ \alpha(0) $ and $ G_{\mu} $ schemes. As shown in the figure, the scattering angle distribution exhibits strong peaks in the forward and backward directions, indicating that Z-boson pairs are predominantly produced along the directions of the electron and positron beams. The EW correction suppresses the LO differential cross section across the entire range of $ \cos\theta $. The corresponding relative correction is approximately –10% in the $ \alpha(0) $ scheme and approximately –7% in the $ G_{\mu} $ scheme at $ \cos\theta = 0 $, decreasing by approximately 5% as $ |\cos\theta| $ increases to $ 1 $. Conversely, the $ \text{QCD} \otimes \text{EW} $ correction slightly enhances the LO differential cross section, with specific increases of approximately 1.15% in the $ \alpha(0) $ scheme and approximately 0.28% in the $ G_{\mu} $ scheme, particularly in the forward and backward directions. Notably, the forward-backward symmetry in the scattering angle distribution of the final-state Z boson is clearly shown in Fig. 4.

Figure 4.  (color online) LO and NNLO corrected scattering angle distributions of the final-state Z boson, along with the corresponding EW and $ \text{QCD} \otimes \text{EW} $ relative corrections, for $ e^+ e^- \rightarrow ZZ $ at $ \sqrt{s } = 240\; \mathrm{GeV} $ in the $ \alpha(0) $ (left) and $ G_{\mu} $ (right) schemes.

The LO and NNLO corrected transverse momentum distributions of the final-state Z boson, as well as the corresponding EW and $ \text{QCD} \otimes \text{EW} $ relative corrections, at $ \sqrt{s} = 240\; \mathrm{GeV} $ in the $ \alpha(0) $ and $ G_{\mu} $ schemes, are presented in Fig. 5. The LO differential cross section increases smoothly with increasing $ p_T $ in the low $ p_T $ region, while it rises sharply with $ p_T $ in the high $ p_T $ region. The EW correction provides a moderate enhancement to the LO differential cross section in the low $ p_T $ region, while suppressing it in the high $ p_T $ region, where the EW relative correction decreases rapidly as $ p_T $ increases. The substantial magnitude of the EW relative correction at extremely high $ p_T $ can be attributed to the Sudakov effect. Conversely, the $ \text{QCD} \otimes \text{EW} $ correction results in a slight enhancement of the LO differential cross section across the entire $ p_T $ region. The relative correction decreases very slowly, remaining steady at approximately 1.15% in the $ \alpha(0) $ scheme and approximately 0.28% in the $ G_{\mu} $ scheme. The Sudakov effect is not observed in the NNLO $ \text{QCD} \otimes \text{EW} $ correction because there are no real emission corrections at the $ \alpha\alpha_s $ order.

Figure 5.  (color online) Same as Fig. 4 but for transverse momentum distributions of the final-state Z boson.

Based on the discussion of the scattering angle and transverse momentum distributions, we conclude that the NNLO $ \text{QCD} \otimes \text{EW} $ corrections will have a significant impact on the interpretation of future high-precision experimental measurements. In the $ \alpha(0) $ scheme, the $ \text{QCD} \otimes \text{EW} $ relative correction exceeds 1.15%, while in the $ G_{\mu} $ scheme, it amounts to approximately 0.28% across most of the phase space. Near the top-pair production threshold, the $ \text{QCD} \otimes \text{EW} $ relative correction is marginally reduced. However, it still exceeds approximately 1.12% and 0.25% in the $ \alpha(0) $ and $ G_{\mu} $ schemes, respectively.

The production of Z-boson pair serves as a crucial avenue for comprehending the EW gauge structure, investigating anomalous TGCs, and searching for new physics. Inspired by the anticipated high-precision measurements and cleaner hadronic environment of future lepton colliders, we present a detailed calculation of NNLO $ \text{QCD} \otimes \text{EW} $ corrections to Z-boson pair production at electron-positron colliders and provide a more precise theoretical prediction in this paper. Using the canonical differential equation method, we perform a comprehensive analytic calculation of all two-loop MIs arising from the $ {\cal{O}}(\alpha\alpha_s) $ corrections. Specifically, we derive the canonical MIs for the massive two-loop triangle MIs and express them in terms of GPLs up to the order of $ \epsilon^4 $. The analytic expressions offer a faster and more accurate numerical evaluation for any kinematic configuration when compared to purely numerical methods, such as sector decomposition, which is essential for precise theoretical predictions and phenomenological analysis.

We apply our analytic expressions for these two-loop canonical MIs to compute the NNLO $ \text{QCD} \otimes \text{EW} $ corrected integrated cross sections for $ e^+e^- \rightarrow ZZ $, as well as the differential distributions with respect to the scattering angle and transverse momentum of the final-state Z boson. Our findings indicate that the NNLO $ \text{QCD} \otimes \text{EW} $ corrections are non-negligible, leading to an enhancement of the LO cross section across the entire phase space. In the $ \alpha(0) $ scheme, the $ \text{QCD} \otimes \text{EW} $ relative correction reaches approximately 1.15%, whereas in the $ G_{\mu} $ scheme, it is reduced to approximately 0.28%. This appears to indicate that the perturbative convergence in the $ G_{\mu} $ scheme is superior to that in the $ \alpha(0) $ scheme. To align with the anticipated permille level of accuracy in experimental measurements at future lepton colliders, these NNLO $ \text{QCD} \otimes \text{EW} $ corrections should be incorporated into the theoretical predictions.

Zhe Li is grateful to Ming-Ming Long for collaboration in the early stages of the project.

In this appendix, we present the explicit analytic expressions of $g_{i}\; (i = 1, \; 2,\dots, 26)$ in terms of GPLs up to $ {\cal{O}}(\epsilon^2) $ as follows: