Probing the CP properties of top Yukawa coupling at an $e^+e^-$ collider

We study consequences of CP violation in the $ht\bar{t}$ Yukawa coupling through the process $e^+e^- \to$ $h(125)t\bar{t}$. The helicity amplitudes are calculated in the $t\bar{t}$ rest frame, where the initial $e^+e^-$ current and the final Higgs boson have the same three-momentum. CP-violating asymmetries appear not only in the azimuthal angle between the $e^+e^-$ plane and the $t\bar{t}$ plane about the Higgs momentum direction, but also in the correlated decay angular distributions of $t$ and $\bar{t}$. Complete description of the production and decay angular distributions are obtained analytically, including both leptonic and hadronic decays of $t$ and $\bar{t}$. We study the ultimate sensitivity to the CP-violating $ht\bar{t}$ coupling at a few center-of-mass energies. Our analysis shows that the possibility of discovering CP violating $ht\bar{t}$ coupling improves significantly by studying $t\bar{t}$ decay angular correlations, and more importantly, by increasing its energy upgrade target from $\sqrt{s} = 500$ GeV to 550 GeV.


I. INTRODUCTION
All the current measurements of the Higgs boson h(125) are consistent with predictions of the Standard Model (SM) of particle physics [1,2]. Detailed studies of all its properties should be of at most importance in our search for the physics of the symmetry breakdown, or the origin of the SM.
In this paper, we show how the top quark Yukawa coupling, the htt coupling, can be studied in the process e + e − → htt at linear colliders such as ILC and CLIC. We develop techniques that allow us to perform the optimal measurements of the htt coupling in a clean e + e − collider environment. In order to quantify the impacts of our proposal to measure full production and decay angular correlations of both semi-leptonic and hadronic decays of top and anti-top quarks, we adopt the following simple effective Lagrangian for the htt coupling L = −g htt ht [cos ξ htt + i sin ξ htt γ 5 ] t (1) with two real parameters, g htt and ξ htt . The effective Lagrangian Eq. (1) reduces to the SM htt coupling when with the top quark mass m t and the vacuum expectation value (VEV) v ≃ 246 GeV of the SM Higgs field. Nonzero values of sin ξ htt term in Eq. (1) implies CP violation, while the magnitude of the g htt is measured with respect to its SM value The sign of cos ξ htt term, or that of g htt cos ξ htt , is measured with respect to the sign of the hZZ coupling, because the amplitudes with the htt couplings and those with the hZZ coupling interfere in the process e + e − → htt. Therefore, the amplitudes depend on the three parameters, κ htt , κ hZZ > 0, − π ≤ ξ htt < π.
It is worth noting that by the time an e + e − collider starts studying the htt production process, we should have constraints on (κ htt , ξ htt , κ hZZ ) from the LHC: κ hZZ can be measured from Γ(h → ZZ * ) and from the weak boson fusion production cross section; (κ htt , ξ htt ) are measured from htt production process [3][4][5][6][7][8]; and (κ htt , ξ htt , κ hZZ ) can be constrained in the single top plus Higgs production process [4,[9][10][11]. Some CP violating asymmetries including top decay lepton distributions are studied in [4]. The couplings are constrained by the perturbative unitarity [12,13], and affect loop induced amplitudes for Γ(h → gg) and Γ(h → γγ). The role of the e + e − collider experiments should be their refinements and possible discovery of non-SM physics in the htt coupling, such as CP violation. Differential cross sections and top polarizations are studied in [14,15], including CP violating observables [14]. Ref. [16] studied asymmetries in decay lepton angular correlations in Higgs plus topponium production at √ s = 500 GeV. At present, only the hZZ coupling strength κ hZZ has been measured to be κ hZZ 0.85 (6) more or less free from detailed model assumptions [2]. We expect the measurement to improve significantly by the time of an e + e − collider experiment, and we set κ hZZ = 1 (7) throughout our analysis. All the results are insensitive to this assumption, which does not change significantly by varying κ hZZ in the 0.85 < κ hZZ ≤ 1 range, because the amplitudes with the hZZ coupling are subleading in e + e − → htt process. We allow the two real parameters (κ htt , ξ htt ) to vary freely in our analysis, with the understanding that they should be constrained significantly by the LHC experiments when the e + e − collider experiments are performed. Before starting our studies, we find it instructive to examine a very specific limit of CP violation in the Higgs couplings where the sole origin of CP violation in the Higgs sector is in the Higgs potential or the Higgs self interactions, while all the other couplings including Yukawa couplings are CP conserving. It is like a milli-weak theory of CP violation for the neutral K system [17], where CP violation is confined to the K 0 −K 0 mixing. The scenario can be realized in any multi-doublet models, where the observed Higgs boson, h(125), is a linear combination where H and H ′ are CP-even while A is a CP-odd neutral components of the Higgs boson. If all the Yukawa interactions of the current states, H, H ′ and A are CP conserving, the sole origin of all the CP-violating couplings of h(125) is the mixing matrix element O hA . We can choose H to be the fluctuation of the full SM Higgs VEV that gives the weak boson masses, and the orthogonal states H ′ and A are fixed uniquely in two Higgs doublet models (2HDM). Specific form of the htt couplings is model dependent. For instance in type II 2HDM, the hZZ and htt couplings are where tan β is the ratio of the two VEV's. It is clear that in this scenario, κ hZZ = 1 implies O hH ′ = O hA = 0 from the orthogonality of the mixing matrix. Even in this scenario, significant CP violation is possible for the htt coupling for small tan β, if κ hZZ ∼ 0.9. The paper is organized as follows. In Section II, we calculate the helicity amplitudes of the e + e − → htt process in the tt rest frame. In Section III, we show the numerical results of production cross section, invariant mass distribution of the htt production process with next-to-leading-order (NLO) QCD corrections including Coulomb resummation for topponium formation. In Section IV, we introduce the density matrix formalism to express full kinematical distributions including t andt decay angular correlations, including both semi-leptonic (t → blν) and hadronic (t → bdu) decays. By accounting for uncertainties in quark jet-flavor identification, we study how the sensitivity to CP violation increases by measuring full distributions, at √ s = 500 GeV, 550 GeV and 1000 GeV. Summary and conclusions are given in Section V. In Appendix A, we show t andt decay density matrices for both semi-leptonic and hadronic decays. In Appendix B, we review HELAS phase convention that affects our production density matrix elements. As shown in Fig. 1, three Feynman diagrams contribute to the process eē → htt. The first two are with the htt coupling in Eq. (1), and the third is with the hZZ coupling for which we assume the SM value as in Eq. (7) in this report.
We calculate the helicity amplitudes of the process in the tt rest frame, where the initial e + e − current and the produced Higgs boson have the same three-momentum. We denote the helicity amplitudes as where α/2 = ±1/2 gives the electron helicity, σ/2 andσ/2 give the helicity of t andt respectively in the tt rest frame. In order to fix our reference frame unambiguously, we start from the e + e − collision (laboratory) frame, where the four momenta are parameterized as Here m tt = Q 2 is the invariant mass of the tt system, are the Lorentz boost factors between the eē and the tt rest frames. The t andt momenta are parameterized in the tt rest frame which is obtained from the laboratory frame by a rotation of −θ about the y-axis and then by a Lorentz boost along the tt momentum direction, which give eē rest f rame tt rest f rame The z-axis is chosen along the negative of the common momentum direction of the Higgs boson and the e + e − system; see Fig. 2. We find that the helicity amplitudes are most easily and systematically calculated in the frame obtained from the above by a rotation ofφ about z-axis, such that the top momentum has vanishing y-component: the h and e + e − four momenta are unchanged from Eq. (14c) and (14d), whereas the initial e andē momenta are now, We calculate the helicity amplitudes in this specific tt rest frame where both the Higgs and the e + e − system have common momenta in the negative z-axis as in Eq. (14c) and (14d), and the azimuthal angleφ between the tt plane and the eē plane are given to the e andē momenta. We first note that the helicity amplitudes in Eq. (10) can be factorized as where q µ = p µ e + p μ e , √ s is the center-of-mass energy, is the Breit-Wigner propagator factor for γ, Z, and also W for later use. In Eq. (18), the γ and Z polarization factor −g µν is replaced by the sum with along the eē four momentum q µ = p µ eē = √ sγ(1, 0, 0, −β), see Eq. (14d), because of the eē current conservation, q µ L V µ (eē; α) = 0 (21) in the m e = 0 limit. The leptonic currents are simply where u α and v α are two-component chiral spinors with α = −1 for L (α = +1 for R), and σ µ ± = (1, ± σ) are the chiral four-vectors of σ matrices. The gauge couplings are with g z = g/ cos θ W . We find where the polarization direction n is along the e three-momentum in the e + e − rest frame, which is obtained from Eq. (19) by a Lorentz boost with β and γ given in Eq. (13), n = (− sin θ cosφ, sin θ sinφ, cos θ).
The leptonic amplitudes are hence expressed explicitly as in terms of Wigner's D functions in the eē rest frame, where α denotes helicities along n in Eq. (25), while λ denotes those along q (−z) direction. Summing up we find where the three D functions are shown explicitly. In this expression Eq. (27), the leptonic amplitudes L V λα (eē) gives the kinematical dependence on the production scattering angle θ and the azimuthal angleφ in terms of Wigner's D functions, whereas the htt production amplitudeŝ M V ασσ depend only on the tt invariant mass m tt , and the polar angle of the top quark momentumθ. After using the equations of motion for t andt, we find a compact expression forM V ασσ , where are the t and Z Breit-Wigner propagator factors. For λ = ±1, the matrix elements are proportional to the Wigner's D functions in the tt rest frame: for (σ −σ)/2 = ±1, and for (σ −σ)/2 = 0. Here we introduce compact notation, Here, E h = √ s− E tt = √ s− m tt γ and p h = p tt = m tt γβ are the Higgs energy and the momentum in the eē rest frame. We note here that the terms proportional to sin ξ in all the 12 amplitudes Eqs. (30,31,34,35) are either proportional top =Êβ, β, or D 1 t − D 2 t , which are all suppressed strongly near the htt production threshold. It is instructive to note here that the above amplitudes satisfy the CP transformation properties which is illustrated in Fig. 3. The upper plot (a) shows the three momenta and the helicities of e andē in the eē rest frame, as well as those of t andt in the tt rest frame. The common z ′ -axis is chosen along the tt momentum direction in the eē rest frame. The electron momentum is in the x ′ -z ′ plane. The vertical arrows (⊙ and ⊗) show that the t-momentum has positive y-component, or 0 <θ,φ < π/2. Shaded arrows show all the fermion spin directions for α = σ =σ = +1. The middle plot (b) is obtained from (a) by the CP transformation, which reverses the sign of all three momenta, but keeps the spin (hence the helicity is reserved), and exchange particles and antiparticles. Note that the initial eē state is invariant under CP transformation for each helicity. The momentum configuration of (b) is expressed in our reference frame as in the bottom plot (c), in which p tt is along the z ′ axis and the y axis is along p e × p tt . The frame (c) is obtained from (b) by rotations, which do not affect the helicity amplitudes. This gives the identity Eq. (36), which show explicitly how CP-violating term proportional to sin ξ htt should behave under exchange of angular variables: The helicity amplitudesM V ασσ contain all the information about the htt (and the hZZ) coupling that we can probe in the process e + e − → htt. Being complex numbers, however, they are not directly observable in experiments. For instance, the differential cross section dσ = 1 2s with the 3-body phase space withβ in Eq. (12) andβ in Eq. (15), measures only the squared sum of all the helicity amplitudes. With e − (and possibly e + ) beam polarization, the sum of e LēR annihilation (α = −1) and that of e RēL annihilation (α = +1) can be resolved. When we study t andt decay distributions and their correlations, we can measure 15 more combinations of the helicity amplitudes (16 each, including the absolute value squared of the 4 amplitudes for α = −1 and α = +1 with beam polarization). In this section, we illustrate this when both t andt decay semi-leptonically, The helicity amplitude for the process can be expressed (for e αē−α collisions) as with the Breit-Wigner propagator factors and the decay amplitudes The differential cross section for the process Eq. (41) is hence (for unpolarized beams) dσ = 1 2s where the 7-body phase space can be decomposed as with p 2 t and p 2 t as the invariant mass squared of the (blνl) and (bℓ ′νl ′ ) systems, respectively. In the narrow width limit of the top quark, holds and the differential cross section in Eq. (45) can be expressed as The above expression can be expressed as The differential decay density matrices in Eqs. (50) are calculated in Appendix A, and take particularly simple form for semi-leptonic decay 33 is the semi-leptonic branching fractions summed over ℓ = e, µ, τ .
Here,θ * andφ * (θ * and φ * ) are the polar and azimuthal angles ofl (ℓ) in the t (t) rest frame where the polar axis is chosen along the t momentum direction in the tt rest frame. More details are explained in Appendix A. By inserting Eq. (51) into Eq. (49), we find It is now clear that all the real and imaginary parts of the product of the helicity amplitudes M ασσ and its complex conjugates M * ασ ′σ′ including σ ′ = σ andσ ′ =σ can be measured by studying the correlated decays t → blν and t →bℓν at all htt phase space point (m tt , cos θ, cosθ,φ). There are four α |M ασσ | 2 terms, six α Re(M ασ ′σ′ M * ασ ′σ′ ) terms and six α Im(M ασσ M * ασ ′σ′ ) terms. With polarized e beams, α = −1 and α = +1 combinations can be resolved.
In Fig. 4, we show theφ distribution of all the 16 combinations of M ασσ M * ασ ′σ′ for α = −1 case (e LēR annihilation) at √ s = 500 GeV (left) and at √ s = 1000 GeV (right). We set cos θ = 0, cosθ = 0.5 and m tt = 350 GeV at both energies where the ξ htt -dependences are found to be significant. We compare ξ htt = π 4 (solid lines) and ξ htt = − π 4 (dashed lines) for √ s = 500 GeV, and ξ htt = ±0.2 for √ s = 1000 GeV. The top panels show the four absolute value squared |M −σσ | 2 and their sum. CP violation appears as a phase shift in theφ distribution whose sign and magnitude are proportional to ξ htt [16]. The difference is reduced significantly when only the total sum of all squared amplitudes, i.e. the htt distributions are observed. The middle panels show the real part of the six off-diagonal (σ ′ = σ orσ ′ =σ) products. In the absence of CP violation, these are even functions ofφ. The CP-violating asymmetries appear again as a phase shift or asymmetries betweenφ > 0 andφ < 0. From Eq. (52), we learn that the 6 real terms are measured as coefficients of cos φ * , cosφ * , cos(φ * + φ * ) or cos(φ * − φ * ).
In the bottom panel we show the corresponding imaginary parts of the six off-diagonal products. They are measured as coefficients of sin φ * , sinφ * , sin(φ * + φ * ) or sin(φ * − φ * ) according to Eq. (52). It is worth noting that these distributions are odd functions ofφ if there is no CP violation (ξ htt = 0). CP violation induces a phase shift in these distributions whose sign and magnitude are proportional to ξ htt .
We also study all the distributions for α = +1 (e RēL annihilation), but they are found to be very similar to the α = −1 case shown in Fig. 4, with significantly smaller magnitudes. Although the results shown in Fig. 4 are for a particular kinematical configuration of htt (m tt = 350 GeV, cos θ = 0, cosθ = 0.5), and for α = −1 (e LēR annihilation), their dependence on the sign of ξ htt shows the possible improvement in the CP-violation discovery potential in e + e − collision experiments, by making use of all the information given by t andt decay angular correlations.

III. CROSS SECTIONS AND QCD CORRECTIONS
In this section we study the energy dependence of the total cross sections and the QCD higher-order corrections, perturbative NLO corrections and resummation of Coulombic corrections that account for topponium formation below and around the threshold, m tt ∼ 2m t . Those studies are made in the two CP-conserving limits, at ξ htt = 0 (h = H, the SM limit), and at ξ htt = π 2 (h = A, the pseudo scalar limit). The results are used to normalize our statistical analysis in the next sections, which are based on the leading-order (LO) matrix elements.  A. Leading-order production cross section We show in Fig. 5 the leading-order total cross section of the e + e − → htt process in the two limits, the pure CP-even (h = H) limit when ξ htt = 0, κ htt = κ hZZ = 1, and the pure CP-odd (h = A) limit when ξ htt = π 2 , κ htt = 1, κ hZZ =0. σ(Htt) is about 0.28 fb at √ s = 500 GeV, reaching 2 fb at √ s ∼ 600 GeV and stays above 2 fb until √ s ∼ 1 TeV. On the other hand, σ(Att) is about 0.0045 fb (below the scale of Fig. 5) at √ s = 500 GeV, rising quickly with energy, reaching 0.43 fb at √ s = 1000 GeV. Because the CP-violating asymmetries appear as interference effects between CP-even and CP-odd amplitudes, we show also the ratio of the two cross sections, R = σ(Att)/σ(Htt) in Fig. 5. We can very roughly expect that the CP asymmetry is proportional to √ R. The scale of R is given along the right-hand vertical axis. It grows rapidly from 0.016 at √ s = 500 GeV to 0.047 at √ s = 550 GeV, reaching 0.2 at √ s = 1000 GeV.

B. QCD corrections
Shown in Fig. 6 are the differential cross sections dσ/dm tt v.s. m tt for h = H (black solid curves) and for h = A (red dashed curves) at √ s = 500 GeV (left) and at √ s = 1000 GeV (right) calculated in the leading order. The h = A cross section at √ s = 500 GeV is multiplied by 10, in order to show the shape difference between the CP-even (h = H) and CP-odd (h = A) limits. It is clear from the two cases shown in the figure that the ratio of the CP-odd and CP-even amplitudes squared is large at low m tt , at all energies. This suggests that the sensitivity to CP asymmetry is high at low m tt , and hence the corrections including topponium formation can have significant impacts on our study of CP violation in the top Yukawa coupling.
We show in Fig. 7 the differential cross section dσ/dm tt for the SM Higgs boson (h = H) with QCD corrections. NLO QCD corrections to the process are evaluated by using MadGraph5 aMC@NLO [18]. In addition, we also consider the corrections by Coulomb resummation [16,[19][20][21]. According to Refs. [21,22], we estimate the Coulomb resummation corrections as follows. First, we evaluate the Born-level helicity amplitudes by including the decays of top-quarks, e + e − → htt → hbW +b W − , in which the top-quarks can be off-shell and m tt can be below 2m t . Then, the amplitudes are corrected by multiplying with the S-wave Green function for non-relativistic tt with the energy E = m tt − 2m t . This prescription is justified by the fact that near the tt threshold where the Coulomb resummation is important the amplitudes are dominated by the S-wave component. To evaluate the Green function we employ the NLO QCD potential with the renormalization scale µ = 40 GeV. Finally, the squared amplitudes are corrected by a hard correction factor K ≃ 0.7 which is numerically extracted by matching with the NLO cross section at an intermediate m tt .
In Fig. 7, we show the differential cross section for the pure scalar case (h = H) at √ s = 500 GeV (left) and 550 GeV (right). LO and NLO cross sections for the on-shell tt limit are plotted by black solid and red solid lines, respectively. Born-level and the Coulomb-resummed cross sections with top-quark off-shellness are plotted by black dotted and blue solid lines, respectively. At √ s = 500 GeV, the NLO corrections enhance the total cross section by around 60%, while the Coulomb resummation (including the NLO effects) enhances the cross section by about 120%. At √ s = 550 GeV, the total cross section is enhanced by 20% and 40%, by the NLO and Coulomb resummation corrections, respectively.
In Fig. 8, we plot the ratio of the NLO cross section to the LO cross section as a function of m tt for the pure scalar (h = H) and pseudoscalar (h = A) cases at √ s = 500 GeV, 550 GeV, and 1000 GeV. We find the K-factors for the m tt distribution are almost the same for the scalar and pseudoscalar processes. QCD corrections are large near the threshold, m tt ≃ 2m t , because of the Coulomb singularity, but become almost flat and small at large m tt .
In Table I   at 1000 GeV. The enhancement factors are larger for CP-odd case (h = A), which is consistent with the softer m tt distribution of the Att process as shown in Fig. 6, where the K-factor is large at small m tt at all energies; see Fig. 8.
The NLO K-factor is smaller than unity for h = H at √ s = 1000 GeV, because the m tt dependent K factor becomes smaller than unity at m tt ≥ 650 GeV above which the m tt distribution is large; see Fig. 6.
Finally in Table II, we show the total cross sections for the CP-even (h = H) and CP-odd (h = A) limits at LO, NLO with stable top quarks, and after Coulomb resummation including the off-shell top quark effects. Despite the factor of 3 enhancements (see Table I), σ(Att) remains tiny, 0.013 fb at 500 GeV. It grows by a factor of 6 to 0.075 fb at 550 GeV, which has a significant impact on our CP-violation search at e + e − colliders.
Although the NLO and topponium corrections are quite large, the NNLO effects to the total cross sections are expected to be marginal, since potentially large higher order corrections in the tt threshold regions are incorporated by Coulomb summation. QCD corrections to the angular distributions or correlations of top quark decay products, which may affect determination on ξ htt , are not studied in this paper. Detailed analysis is desired in the future.  the measurement accuracy depends on the total number of produced events that determines the statistical errors, we first estimate the total cross section as a function of the parameters of our effective Lagrangian (κ htt , ξ htt , κ hZZ ) in subsection IV A. In the next subsection IV B, we explain in detail how we can calculate the full differential distribution including both semi-leptonic and hadronic decays of t andt, that are observable by a perfect detector but not as perfect as capable of distinguishingd (d) jets from u (ū) jets. In subsection IV C, we introduce a very simple estimator χ 2 function that measures all differences in the observable distributions between experiment (κ ex htt , ξ ex htt ) and theory (κ th htt , ξ th htt ), and study the potential of rejecting ξ th htt = −ξ ex htt (observation of CP violation) as a function of |ξ ex htt |.

A. Total cross section
The total cross section of the process e + e − → htt depends not only on the parameters of (κ htt , ξ htt , κ hZZ ) but also on the center-of-mass energy √ s and the beam polarization. Since the NLO and topponium formation corrections are obtained only for the SM limit (κ htt , ξ htt , κ hZZ ) = (1,0,1) and for the purely CP-odd limit (κ htt , ξ htt , κ hZZ ) = (1, ±π/2, 0) for unpolarized beams at √ s = 500 GeV, 550 GeV, 1000 GeV, we make the following simple parameterization Here σ NLO H and σ NLO A are obtained from Table I, whereas we quote the difference between σ NLO and the total cross section after taking account of Coulomb resummation as the 'topponium' cross section, Those cross sections values are tabulated in Table III. All the coefficients of our parameterization Eq. (53), which are tabulated in Table IV, are obtained by using the LO matrix elements as follows. We calculate the total cross sections σ L and σ R , respectively, for purely left-handed (e L ) and right-handed (e R ) beam in the LO for several sets of (κ htt , ξ htt , κ hZZ ) parameters, and obtain the parametrisation Eq. (53)  , respectively. We note here that this is not accurate, and the NLO corrections should be made separately for the (κ htt cos ξ htt ) 2 and the (κ hZZ ) 2 term as well as the interference term proportional to (κ htt cos ξ htt κ hZZ ) in the future. For the topponium coefficients (A ′ α , B ′ α , C α ′ , D ′ α ), we calculate the LO cross sections at m tt = 2m t + 0.1 GeV, just above the threshold, and fix all the coefficients, normalizing the total cross section to the topponium cross sections obtained in Section III and tabulated in Table III. Accordingly, the coefficients in Table IV are normalized as  We believe that this is a good approximation to the S-wave topponium formation, which dominates the coulomb resummation corrections. The cross section for partially polarized beams (|P e | < 1 for e − , |Pē| < 1 for e + ) is obtained from Eq. (53) as and hence the unpolarized cross sections are Here, we give a few comments on the impacts of the beam polarization on the htt coupling measurements. We first note from Table IV that the ratio D α /A α and D ′ α /A ′ α are almost the same for e L and e R , and hence the beam polarization has little impacts (besides the total number of events which can be increased by choosing P e < 0 < Pē) in the measurements of the sign of ξ htt that arise from the interference between the amplitudes whose squares give A α and D α , respectively. The most significant effect of the beam polarization is in the value of B α , the coefficients of the interference between two CP-even amplitudes, those with the htt coupling proportional to κ htt cos ξ htt and those with the hZZ coupling. Their signs are opposite at √ s = 500 GeV and 550 GeV, while the magnitude is different by a factor of six at √ s = 1000 GeV. Accordingly, we find significant improvement in resolving ξ htt = 0 (cos ξ htt = 1) and ξ htt = ±π (cos ξ htt = −1) by using the beam polarization. However, because the beam polarization has little impacts on the CP-violation measurements, the following studies are performed for unpolarized beams (P e = Pē = 0). We would like to study the impacts of the e + e − beam polarizations in a separate report, in which the sensitivity should be compared with those from single top and h production at the LHC [11].
Throughout this study, we adopt as a nominal integrated luminosity at each colliding e + e − energy √ s = 500 GeV, 550 GeV, and 1000 GeV so that we can compare the impacts of increasing the beam energy. For unpolarized beams, we expect the following number of htt events with unpolarized beams for the SM (h = H): which measures the total cross section with statistical errors of 4.0%, 2.6%, 2.2%, respectively. If we assume that these number of events are observed, then at each energy, the value of (κ htt , ξ htt ) are constrained to lie on the curves shown in Fig. 9. We set κ hZZ = 1 throughout our studies. The statistical error on κ htt is about half the cross section error, ∼ 2.0%,∼ 1.3%, ∼ 1.1%, respectively, for √ s = 500 GeV, 550 GeV, 1000 GeV, which are too small to show in the plot. It is clear that the measurement of the total cross section at two different energies can constrain a region of κ htt and ξ htt , and more importantly, the LHC experiments on htt and ht, ht productions will also give such constraints. We therefore assume that the SM Higgs cross section is observed at a given e + e − collision energy, and study the capability of distribution studies to resolve the sign of ξ htt as a function of |ξ htt |. Therefore, one should be careful when comparing ξ htt = π 4 against ξ htt = − π 4 at √ s = 500 GeV, since it is not κ htt = 1 but κ htt = 1.4 which gives the same number of events as in the SM.

B. Differential cross sections
The differential cross sections including t andt decay distributions are obtained as in Eqs. (49) by using the density matrix formalism. In this subsection we explain how we treat the hadronic decays of t andt, QCD corrections to the differential cross sections which strongly depend on m tt , possible beam polarizations, and how we take account of the probability that t andt cannot be distinguished uniquely when both of them decay hadronically.
First, we introduce the following generic form of the t andt decay density matrices respectively, wherex are normalized 4-momenta ofl (d ors) and ℓ (d or s), respectively, in the t andt rest frame. They are expressed in terms of the b (b) angles θ * b , φ * b (θb,φ * b ) in the t (t) rest frame, thel,d,s anglesθ * * ,φ * * in the W + rest frame, and the ℓ, d, s angles θ * * , φ * * in the W − rest frame. Details are given in Appendix A. The density matrices in Eq. (60) are normalized as We use the density matrices Eq. (60) for all leptonic decays. Indeed, integration of Eq. (60) over the b (b) angles for fixed angles ofl (ℓ) in the t (t) rest frame reproduce the well known formula Eq. (51) for B ℓ = 1. For hadronic decays, we introduce the probability Pd u with which we identifyd from u (d fromū), ors from c (s fromc) in W − (W + ) decays where we assume that b andb momenta are uniquely identified. With the probability (1 + Pd u )/2, we identifyd ors (d or s) correctly in W + (W − ) decays, whereas with the probability (1 − Pd u )/2, we misidentify u or s (ū ors) ford ors (d or s). In our analysis, we set By using ρ andρ in Eq. (60) for leptonic decays and ρ h andρ h in Eq. (63) for hadronic decays, the differential cross sections are expressed as follows with the reduced phase space In this report, we do not distinguish lepton or quark flavours, and adopt Production density matrices are expressed in terms of the helicity amplitudes as where −1 < P e , Pē < 1 denote e andē longitudinal beam polarizations. There is one subtlety when both t andt decay hadronically. In this case, identification of t andt can be ambiguous, and we introduce the probability P tt with which t andt can be identified correctly. In this report, we adopt which is approximately twice the semi-leptonic decay branching fraction of B mesons. Charge discrimination of hadronic jets from W + and W − decays may also help. With this probability, the distribution proportional to In the last term with the probability (1 − P tt )/2, t andt are misidentified, not only the t angles (θ,φ) in the tt rest frame are replaced by those oft angles (π −θ, π +φ), but also all the 4-angles of t andt decays are exchanged. Although this looks complicated, it is straightforward to implement it in a numerical program, and we can keep the maximum surviving information of the matrix elements. Finally, we implement the perturbative NLO corrections and the topponium contribution as follows. For the NLO corrections, we find in Section III that the m tt -dependent K-factors are almost identical to the CP-even and CP-odd cases, see Fig. 8. We therefore ignore their small differences, and multiply the m tt -dependent K-factor of the CP-even (h = H) cross section to our differential cross section in Eq. (65) We confirm that this simple prescription reproduces all the NLO cross sections listed in Table III. For the topponium contribution at √ s = 500 GeV and 550 GeV, we evaluate the differential cross section Eq. (65) at m tt = 2m t + 0.1 GeV, and normalize the total integral such that it agrees with the production cross section of topponium + Higgs listed in Table II,   .
With this prescription, we obtain all the correlated decays of both t andt decays, reproducing the results of Ref. [16]. Sophisticated simulation program with topponium formation and decays is needed to evaluate the cross sections and distributions for the SM and its extensions. We believe, however, that the model dependences are approximated reasonably well by using the leading-order matrix elements as outlined in this subsection.

C. Results
As an estimator for the sensitivity of our CP-violation measurement in e + e − → htt process, we introduce the following χ 2 function Here dσ ex /dΦ represents the observed experimental cross section calculated by assuming that a set of parameters (κ ex htt , ξ ex htt ) are true and dσ th /dΦ is calculated for an arbitrary set of (κ htt , ξ htt ) values at each energy. L is the integrated luminosity of the process for unpolarized beams. Throughout our study, we set L = 1000 fb −1 , P e = Pē = 0, at all energies. The χ 2 function in Eq. (73) accounts for all possible difference between data dσ ex /dΦ and theory dσ th /dΦ at all kinematical configuration Φ, with the corresponding statistical error proportional to (Ldσ th /dΦ) −1/2 . It maybe regarded as an ultimate possible sensitivity for a perfect detector with infinite resolution and no errors. Because we calculate the data dσ/dΦ by using our theoretical formula without fluctuations, it is obvious that When we fix the test value of ξ htt , the minimum becomes where κ htt (ξ htt ) is the value which gives the same total cross section with the SM. The trajectory of κ htt (ξ htt ) has been given in Fig. 9 for √ s = 500 GeV, 550 GeV, and 1000 GeV.
In the left panel of Fig. 10, we show χ 2 min (ξ htt ) as a function of ξ htt when the SM distribution (κ htt , ξ htt ) = (1, 0) is assumed for dσ ex /dΦ. The blue dashed curve is obtained when only the htt distributions are measured. In our differential cross section formula in Eq. (71) with K(m tt ) for NLO correction, these limits are obtained simply by replacing all our decay density matrices by a unit matrix: As compared to the simple calculation of e + e − → htt differential cross section, we find slightly smaller χ 2 min because of finite tt discrimination probability P tt = 0.4 when both t andt decay hadronically. The red dash-dotted curve   shows our results when leptonic decay angular correlations are assumed. This limit is obtained in our dσ/dΦ formula by setting only the hadronic decay density matrix to be a unit matrix.
Only small improvements are found over the htt distributions only case. This is mainly because of the smallness of the branching fraction factor of the dilepton case B 2 ℓ ∼ 0.11 even including the ℓ = τ and ℓ =τ modes, even though they give full t andt decay angular correlations. It is also because the single (t ort) decay angular distributions with the probability of 2B ℓ (1 − B ℓ ) ∼ 0.44 have relatively small sensitivity to the CP phase, ξ htt . The green solid curve is obtained by including both leptonic and hadronic decay angular correlations by using our full differential cross sections dσ/dΦ in Eq. (71). The significant improvement over the red dash-dotted curve shows that the angular correlation between semi-leptonic and hadronic decays keep the resolving power even with the assumption of Pd u = Ps c = 0 in Eq. (64). Finally, the red solid curve is obtained by adding the topponium contributions. Although the full decay angular distribution studies based on our formalism improve the measurement significantly, it is clear from the plot that the SM data at √ s = 500 GeV is consistent with |ξ htt | = π/4 at 1σ, even with 1000 fb −1 .
In the right panel of Fig. 10, we show the results when the data follow the prediction of the maximum CP phase at (κ ex htt , ξ ex htt ) = (1.4, π/4), where the total cross section is the same as the SM. Here the χ 2 min value at ξ htt = −ξ ex htt = −π/4 gives the sensitivity to CP violation. Apparently, the htt distribution has no sensitivity (see the blue dashed curve), as may be expected from Fig. 4 in Section II, where the difference between ξ htt = π/4 and −π/4 is tiny for the sum of squared of all the amplitudes. The inclusion of leptonic decays improves (red dash-dotted line), and the impact of including hadronic decays (green solid line) is quite significant. Adding the topponium contributions, the χ 2 value reaches 4, or we might find a 2σ hint of CP violation.
We find the results discouraging, because this result is for the largest possible CP phase, with 1000 fb −1 , and a prefect detector with no fluctuations in data are assumed in the analysis. This leads us to re-examine the target energy of a linear collider. As explained in Section II , the disappointing results at √ s = 500 GeV are probably unavoidable, because the small ratio R = σ(Att)/σ(Htt) = 0.016 at √ s = 500 GeV shown in Fig. 5 implies that the CP-odd amplitudes are much smaller than the CP-even amplitudes. Our analytic form of the helicity amplitudes obtained in Section II in Eq. (27) and Eqs. (30-35) confirms that all the amplitudes proportional to sin ξ htt are either proportional toβ, the t velocity in the tt rest frame, β, the velocity of the tt system in the eē rest frame, or D 1 t − D 2 t , the difference between the two top quark propagator factors, which are all strongly suppressed at energies near htt threshold, m h + 2m t = 471 GeV.
We therefore examine the possibility of upgrading the target energy by 10% to 550 GeV, where the ratio R = σ(Att)/σ(Htt) increases to 0.047 ∼ (0.22) 2 . The results are shown in Fig. 11. In the left panel, we show χ 2 (ξ htt ) when the SM (κ ex htt , ξ ex htt ) = (1, 0) is assumed for the data. Again, the htt distributions (modified for the probability  h (1 − P tt ) = 0.27 that t andt cannot be resolved) shown by blue dashed curve give little sensitivity to ξ htt when |ξ htt | ≤ π/4. With the inclusion of t ort semi-leptonic decay angular correlations, the sensitivity shown by the red dash-dotted curve shoots above 1σ at |ξ htt | = π/4. Moreover with inclusion of the hadronic decay angular correlations, as shown by green solid curve χ 2 grows above 5σ. The topponium contribution is not very significant as may be expected from the relatively small topponium formation cross section at √ s = 550 GeV. It may be worth noting that ξ htt = 0 and π can be resolved at 2σ level by using decay angular corrections even without beam polarization.
In the right panel, we show the results for (κ ex htt , ξ ex htt ) = (1.38, π/4). As in the √ s = 500 GeV case, htt distribution has almost no power in resolving CP violation (blue dashed curve). χ 2 min value at ξ htt = −ξ ex htt = −π/4 becomes 4 (2σ) if we study lepton angular correlations shown by red dash-dotted curve. It jumps to above 4σ once we include hadronic decay angular corrections. We can now hope for a meaningful measurement of CP violation in the top Yukawa coupling at √ s = 550 GeV with dedicated efforts. In Fig. 12, we show the difference as a function of ξ ex htt , which gives a measure of the sensitivity of the e + e − → htt experiments to CP violation. We believe that by the time the proposed experiment can be done at the ILC we should have relatively strong constraints on the magnitude of ξ htt from the htt and single top+h production cross sections and their distributions at the LHC experiments. The role of ILC experiments should hence be to test that there is indeed CP violation in the Higgs Yukawa couplings, when nonzero magnitude of ξ htt is favored by those data. This can only be done by observing CP violating asymmetries which determine the sign of ξ htt , whose significance is proportional to the difference in Eq. (79).
In the figure, the blue dashed curve shows that the htt distribution has almost no sensitivity to CP violation. Inclusion of leptonic decay angular correlation (red dash-dotted line) improves, but the most significant improvement is found by including hadronic decay angular correlations as shown by green solid curve. We find that the contributions from the modes where one of t andt decays leptonically, and the other decays hadronically, whose branching fractions account for 2B ℓ B h ≃ 0.44, increase the χ 2 min value most significantly. In this mode, t andt are uniquely identified, and the leptonic decay angular distribution is exact. Even though the hadronic decay angular correlations suffer fromd v.s. u (ors v.s. c) misidentification, significant correlation between the t andt decays remains to resolve CP violation. We may tell from Fig. 12 that a 2σ hint of CP violation can be found if ξ htt > 0.12π ∼ 0.38, while a 3σ evidence of CP violation can be found if ξ htt > 0.18π ∼ 0.55 at √ s = 550 GeV with L = 1000 fb −1 .
We note in passing that our χ 2 function Eq. (79) adopts the experimental distribution dσ ex /dΦ calculated analytically for ξ htt = ξ ex htt without fluctuation for realistic binned data. The significance which can be read off from  to obtain realistic estimates with statistical fluctuation of experimental data. 1 We close our discussions by repeating the study at √ s = 1000 GeV with L = 1000 fb −1 . The left panel of Fig. 13 shows the constraints on ξ htt , when the SM is assumed for the data. At √ s = 1000 GeV, the htt distribution (blue dashed curve) has significant sensitivity, and it alone can reject |ξ htt | = π/4 at 5σ level. The inclusion of both leptonic and hadronic decay angular correlations (green solid curve) doubles the χ 2 min value. In the right panel, we show the constraint when (κ ex htt , ξ ex htt ) = (1.2, π/4) is assumed for data. χ 2 min value at ξ htt = −ξ ex htt = −π/4 tells that CP violation can be seen at 5σ even with htt distribution only (blue dashed curve), while it exceeds 10σ once we study full decay angular correlations.
The discovery potential of CP violation at √ s = 1000 GeV with 1000 fb −1 is summarized in Fig. 14. From the green solid curve that is obtained from our full correlation studies, we may hope to find a 3σ evidence for CP violation if ξ ex htt ≥ 0.074π ∼ 0.23, and a 5σ discovery if ξ ex htt ≥ 0.12π ∼ 0.38.

V. SUMMARY AND DISCUSSIONS
In this report, we study the potential of future e + e − collider experiments to discover CP violation in the top Yukawa coupling, whose magnitude is largest among all the couplings in the SM. We examine e + e − → htt production with all t andt decay modes, including both semi-leptonic (t → blν,t →bℓν) and hadronic decays. The full kinematical distributions including the decay angular correlations are obtained by using the density matrix formalism. In addition to the well known t andt decay density matrices for semi-leptonic decays, where the full information on the t andt polarization is given by thel and ℓ decay angular distributions in the t andt rest frame respectively, we introduce novel decay density matrices for the hadronic decay modes which keep significant t andt polarization resolving power even when one cannot distinguish betweend v.s. u ors v.s. c jets (d v.s.ū or s v.s.c jets) in W + (W − ) decays.
We approximate the full differential cross section by using the decay angular correlations of the leading-order matrix elements with CP-violating htt coupling phase ξ htt , corrected for the m tt dependent NLO correction factor and the topponium formation at m tt ∼ 2m t . QCD corrections to the angular distributions or correlations of top-quark decay products have not been studied in this paper. Detailed analysis will be required in future work. The sensitivity to CP violation has been estimated by comparing the full distributions between ξ htt and −ξ htt in the range 0 ≤ |ξ htt | < π/4, by adjusting the magnitude of the htt coupling (κ htt ) to reproduce the total cross section of the SM Higgs. We study the case at √ s = 500 GeV with L = 1000 fb −1 , and find that even though the semi-leptonic and hadronic decay angular correlations (as well as topponium contributions) give sensitivity to the sign of ξ htt , it reaches only up to the 2σ level after including all the contributions for the maximum CP phase, |ξ htt | = π/4. We therefore examine the possibility of increasing the beam energy by 10% to √ s = 550 GeV. The sensitivity grows significantly, giving us the possibility of a 2σ hint for |ξ htt | ∼ 0.12π ∼ 0.38 and that of a 3σ evidence for |ξ htt | ∼ 0.18π ∼ 0.55. We therefore propose that the next target energy of a linear e + e − collider should be increased from the original 500 GeV to 550 GeV. Impacts of such a 10% increase in the beam energy on other key measurements (such as the hhh coupling) should be studied.
At √ s = 1000 GeV, we expect a 3σ evidence at |ξ htt | ∼ 0.074π ∼ 0.23, and a 5σ discovery at |ξ htt | ∼ 0.12π ∼ 0.38. All these numbers should be regarded as possible maximum sensitivity by a perfect detector that can resolve all information given by the matrix elements.
We believe that the t andt decay density matrices we obtain for both semi-leptonic and hadronic decay will be a powerful tool in all processes with t and/ort in the final states. Their distributions are explained in detail in Appendix A. The use of numerical helicity amplitudes evaluated by amplitude calculators can also have wide applications, especially for all the processes whose amplitudes are so complicated that analytic expressions are of little use. We therefore give in Appendix B a review of the fermion wave function phase convention in Refs. [24,25] that is adopted in MadGraph [26].
in the limit of neglecting all the final fermion masses including the b-quark.
The decay matrix element for the semileptonic decay (A1) can be written as with the W propagator factor, D W (q) = 1/(q 2 − m 2 W + im W Γ W ), and σ µ ± = (1, ± σ) are the chiral four vectors of the Pauli matrices. By using the Fierz identity we can express the amplitudes as where the bi-fermionic products in each big bracket are separately Lorentz invariant 2 . The spinor products are calculated easily as in the b + ν rest frame, while in the t rest frame we find for the top helicity σ/2. The amplitudes (A7) are hence simply We normalize the top quark decay density matrix distribution as so that the total integral gives Note that the trace of (A12) gives twice the branching fraction. Inserting the amplitudes (A10) into (A11), we find In the case of semi-leptonic decays (A1), it is most convenient to parameterize the invariant 3-body phase space as The invariance of each bi-spinors is manifest in the spinor rotation where undotted lower indices give left-handed spinors and the dotted indices are for their complex conjugates. ǫ kℓ = ǫkl are antisymmetric sign factors with ǫ 12 = −ǫ 12 = 1.
where θ * * b denotes the polar angle of the b quark in the b + ν rest frame measured from the − pl direction. In the zero width limit of the W , we can integrate out cos θ * * b as We integrate outx in the region m 2 W /m 2 t <x < 1, and find In the matrix form reproducing the normalization Eq. (A12).
In the case of the hadronic decay (A2) of the top quark, matrix elements are exactly the same as those of the semileptonic decays Eq. (A10), where (x,θ * ,φ * ) are now the normalized energy and the angles of thed quark (or s in the case of W →sc decay); see Eq. (A4). We need, however, a careful treatment of the decay phase space integral, because it is difficult to identifyd from u (s from c) in collider experiments. We introduce a parameter Pd u (0 ≤ Pd u ≤ 1), which measures the probability thatd v.s. u (s v.s. c) are distinguished correctly, in the decay density matrix distributions. Although all the results presented in this report are for Pd u = Ps c = 0 (absolutely no distinction between the two decaying jets of the W ), we hope that our formalism encourages further efforts fors v.s. c discrimination (or evend v.s. u).
We start with the differential density matrices Eq. (A13), multiplied by a factor of 3 for the color, where (x,θ * ,φ * ) are parameterizingd (s) quark momentum in the top quark rest frame. Becaused and u (s and c) are difficult to distinguish, we parameterize the t → bdu phase space by using the W →du rest frame angles: where θ * b and φ * b give the b quark momentum in the top quark rest frame, whileθ * * andφ * * give thed momentum in the W →du rest frame which is obtained from the top quark rest frame by rotations (−φ * b about the z-axis, and then −θ * b about the y-axis) and a Lorentz boost along the z-axis by The integration over m 2 du is done in the zero W width limit to obtain and we find Here,x µ = (x,x x ,x y ,x z ) is the normalizedd quark four momentum in the top quark rest framē where γ * = mt With the above parameterization, we can express the top quark decay density matrix distribution under realistic experimental conditions with finite (or zero) probability Pd u ford v.s. u (s v.s. c) discrimination as follows The meaning of Eq. (A24) is that, ifd v.s. u discrimination is perfect, i.e. Pd u = 1, then the density matrix distribution is exactly the same (up to the branching fraction factors) as that for the semi-leptonic decays (A17). When P dū = 0, which we assume in our numerical studies, thed quark (atθ * * andφ * * ) is not distinguished from the u quark (at π −θ * * and π +φ * * ) in the W rest frame, and hence their average measures the top quark polarization.
With the probability P dū = Pd u for d v.s.ū discrimination, thet decay density matrix for hadronic decays becomes dρ(P dū ) = 6B(t →bdu) It is worth noting here that the t → blν andt →bℓν amplitudes as shown in Eq. (A10) and Eq. (A29), respectively, are expressed as squared roots of complex numbers because space rotations on spinors give half angles. These elegant expressions cannot be obtained if one adopts a phase convention that depends on the azimuthal angle, such as the one in Refs. [24,25] adopted by a Feynman amplitude calculator MadGraph [26]. The amplitudes Eq. (A10) and Eq. (A29) can be recovered by supplying the phase factor in Eq. (B7), as will be explained in detail in Appendix B. Because the phase factor associated with thel and ℓ spinors, respectively, for t andt decay amplitudes are common for both helicity amplitudes, the decay density matrices in Eqs. (A17, A22, A30, A32) are independent of thel and ℓ spinor phase convention. and with the charge conjugation v L (p, σ) = iσ 2 u R (p, σ) * , v R (p, σ) = −iσ 2 u L (p, σ) * , we find that straightforward calculation of the Lorentz transformation, from the t and (t) rest frame p µ = (m, 0, 0, 0) to the frame in which the t(t) four momentum becomes Eq. (B2), The fermion wave functions in Ref. [24], adopted in HELAS [25] are obtained by chopping off the phase factors which depend on the azimuthal angle and the helicity of t andt, u(p, σ) = u(p, σ) HZ e −i σφ Hence, for all tt production processes, the helicity amplitudes in the HZ phase convention differ from the naive ones as when t andt momenta have non-zero azimuthal angles, φ andφ, respectively. In our study, we calculate the helicity amplitudes in the tt rest frame wherē holds. 3 Therefore the helicity amplitudes in HZ convention (the HELAS amplitudes) differ from the conventional amplitudes by the following phase factors They introduce unphysical φ-dependences in real and imaginary parts of the density matrices, when t andt momentum has non-zero azimuthal angles. Note that the phase does not appear in the chirality conserving amplitudes in the massless limit, and hence the e + e − currents are free from the azimuthal angle phases. Such convention dependent φ-dependence does not survive in the physical distributions. For instance, the phases in Eq. (B11) are precisely canceled by the corresponding phases in the t andt decay density matrices, if we use the same wave functions for production and decay. In order to recover the rotational invariance of the conventional helicity amplitudes, we may supply the phase factor in Eq. (B8) or simply stick to the frame where t andt momenta have no y-components. We adopt in Section II the latter approach by choosing the frame in which t andt momenta are in the x-z plane, whereas the azimuthal angles are given to e andē momenta.