Probing the CP properties of top Yukawa coupling at an e+e− collider

We study consequences of CP violation in the htt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ ht\overline{t} $$\end{document} Yukawa coupling through the process e+e−→h125tt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {e}^{+}{e}^{-}\to h(125)t\overline{t} $$\end{document}. The helicity amplitudes are calculated in the tt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ t\overline{t} $$\end{document} 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 tt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ t\overline{t} $$\end{document} plane about the Higgs momentum direction, but also in the correlated decay angular distributions of t and t¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{t} $$\end{document}. Complete description of the production and decay angular distributions are obtained analytically, including both leptonic and hadronic decays of t and t¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{t} $$\end{document}. We study the ultimate sensitivity to the CP-violating htt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ ht\overline{t} $$\end{document} coupling at a few center-of-mass energies. Our analysis shows that the possibility of discovering CP-violating htt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ ht\overline{t} $$\end{document} coupling improves significantly by studying tt¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ t\overline{t} $$\end{document} decay angular correlations, and more importantly, by increasing its energy upgrade target from s=500\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s}=500 $$\end{document} GeV to 550 GeV.


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.1) with two real parameters, g htt and ξ htt . The effective Lagrangian eq. (1.1) reduces to the SM htt coupling when

JHEP02(2018)180
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.1) implies CP violation, while the magnitude of the g htt is measured with respect to its SM value (1. 3) The sign of cos ξ htt term, or that of g htt cos ξ htt , is measured with respect to the sign of the hZZ coupling, (1.4) 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 < π. (1.5) 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 (1.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 (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 JHEP02(2018)180 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 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 2, we calculate the helicity amplitudes of the e + e − → htt process in the tt rest frame. In section 3, 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 4, we introduce the density matrix formalism to express full kinematical distributions including t andt decay angular correlations, including both semi-leptonic (t → b¯ ν) 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 5. 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.

Helicity amplitudes and density matrix
As shown in figure 1, three Feynman diagrams contribute to the process eē → htt. The first two are with the htt coupling in eq. (1.1), and the third is with the hZZ coupling for which we assume the SM value as in eq. (1.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  Figure 2. Kinematics of the eē → htt process in the eē rest frame and in the tt rest frame.
The z-axis is chosen along the negative of the common momentum direction of the Higgs boson and the e + e − system; see figure 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. (2.5c) and (2.5d), whereas the initial e andē momenta are now, p μ e = √ s 2 (γ(1 + β cos θ), sin θ cosφ, − sin θ sinφ, −γ(cos θ + β)). (2.8) 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. (2.5c) and (2.5d), and the azimuthal angleφ between the tt plane and the eē plane are given to the e and e momenta. We first note that the helicity amplitudes in eq. (2.1) can be factorized as

JHEP02(2018)180
where q µ = p µ e + p μ e , √ s is the center-of-mass energy, D V (q) = 1/(q 2 − m 2 V + im V Γ V ) is the Breit-Wigner propagator factor for γ, Z, and also W for later use. In eq. (2.9), the γ and Z polarization factor −g µν is replaced by the sum along the eē four momentum q µ = p µ eē = √ sγ(1, 0, 0, −β), see eq. (2.5d), because of the eē current conservation, 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. (2.10) by a Lorentz boost with β and γ given in eq. (2.4), n = (− sin θ cosφ, sin θ sinφ, cos θ). (2.16) 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. (2.16), while λ denotes those along q (−z) direction. Summing up we find where the three D functions are shown explicitly.

JHEP02(2018)180
In this expression eq. (2.18), 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 amplitudesM 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 ασσ , 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. (2.21), (2.22), (2.25), (2.26) 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 M α,σ,σ (θ,θ,φ; ξ) = M α,−σ,−σ (π − θ, π −θ, −φ; −ξ) (2.27) which is illustrated in figure 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. (2.27), which show explicitly how CP-violating term proportional to sin ξ htt should behave under exchange of angular variables: (θ,θ,φ) ↔ (π − θ, π −θ, −φ), (2.28) and when the t andt helicities are exchanged as 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 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. (2.32) is hence (for unpolarized beams) 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 (b¯ ν¯ ) and (b ν¯ ) systems, respectively. In the narrow width limit of the top quark, holds and the differential cross section in eq. (2.36) can be expressed as The above expression can be expressed as The differential decay density matrices in eqs. (2.41) 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 of¯ ( ) 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. (2.42) into eq. (2.40), 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 → b¯ ν andt →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.

JHEP02(2018)180
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. (2.43), we learn that the 6 real terms are measured as coefficients of cos φ * , 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. (2.43). 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 figure 4, with significantly smaller magnitudes. Although the results shown in figure 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.

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.
Because the CP-violating asymmetries appear as interference effects between CPeven and CP-odd amplitudes, we show also the ratio of the two cross sections, R = σ(Att)/σ(Htt) in figure 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.

QCD corrections
Shown in figure 6    (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 figure 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]    In figure 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 1 (see table 1), σ(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

JHEP02(2018)180
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.
4 Sensitivity of e + e − → htt experiments on ξ htt In this section, we study the potential of e + e − → htt experiments to discover CP violation in the top Yukawa coupling at a future linear e + e − collider by postulating a perfect detector with no systematic uncertainties. Because 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 4.1. In the next subsection 4.2, 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 4.3, 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 |.

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 σ α (κ htt , ξ htt , κ hZZ ) = σ NLO Here σ NLO H and σ NLO A are obtained from table 1, 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 3. All the coefficients of our parameterization eq. (4.1), which are tabulated in table 4, 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.   Table 3. e + e − → htt cross sections at NLO and for the topponium formation, whose sum gives the total cross section after Coulomb resummation in table 1. h = H for the SM Higgs, h = A for the CP-odd Higgs with (κ htt , ξ htt , κ hZZ ) = (1, ±π/2, 0).  , 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 3 and tabulated in table 3. Accordingly, the coefficients in table 4 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. (4.1) 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 4 that the ratio D α /A α and D α /A α are almost the JHEP02(2018)180 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 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 figure 9. We set κ hZZ = 1

JHEP02(2018)180
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.

Differential cross sections
The differential cross sections including t andt decay distributions are obtained as in eqs. (2.40) 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 are normalized 4-momenta of¯ (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, thē ,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. (4.8) are normalized as We use the density matrices eq. (4.8) for all leptonic decays. Indeed, integration of eq. (4.8) over the b (b) angles for fixed angles of¯ ( ) in the t (t) rest frame reproduce the well known formula eq. (2.42) for B = 1.

JHEP02(2018)180
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 Pd u = Ps c = 0 (4.12) for simplicity, but Ps c may be significant. The decay density matrices ρ h andρ h in eqs. (4.11) keep the maximum information given by the matrix elements. If we ignore the top spin sensitivity in thed ors (d or s) distributions by integrating over the W decay angles (θ * * ,φ * * ), only poor resolution power of b angular variables θ remains. By using ρ andρ in eq. (4.8) for leptonic decays and ρ h andρ h in eq. (4.11) for hadronic decays, the differential cross sections are expressed as follows with the reduced phase space (4.14) 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 P tt = 0.4 (4.17)

JHEP02(2018)180
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 B 2 h 0.45 in eq. (4.13) should be 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 3 that the m tt -dependent K-factors are almost identical to the CP-even and CP-odd cases, see figure 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. (4.13) We confirm that this simple prescription reproduces all the NLO cross sections listed in table 3. For the topponium contribution at √ s = 500 GeV and 550 GeV, we evaluate the differential cross section eq. (4.13) 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 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.

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, (4.22) at all energies. The χ 2 function in eq. (4.21) 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 χ 2 min = χ 2 (κ ex htt , ξ ex htt ) = 0 (4.23) 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 figure 9 for √ s = 500 GeV, 550 GeV, and 1000 GeV.
In the left panel of figure 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. (4.19) 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

JHEP02(2018)180
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. (4.19). 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. (4.12). 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 figure 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 figure 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 2 , 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 figure 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. (2.18) and eqs. (2.21)-(2.26) 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 figure 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 B  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 figure 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. (4.27).
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 figure 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. (4.27) 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 figure 12 should hence be regarded only as a first optimistic estimate. Dedicated studies with event generation may be needed 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 figure 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 figure 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.

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 → b¯ ν,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 the¯ 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 JHEP02(2018)180 resolving power even when one cannot distinguish betweend v.s. u ors v.s. c jets (d v.s. u 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][26] that is adopted in MadGraph [27].

A Top and antitop quark decay density matrix
In this appendix, we give top quark decay density matrices for helicity amplitudes are calculated). In the frame, we parameterize the four momentum of b and¯ (d) as (1, sinθ * cosφ * , sinθ * sinφ * , cosθ * ) (A. 4) in the limit of neglecting all the final fermion masses including the b-quark.
The decay matrix element for the semileptonic decay (A.1) 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 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 = k˙ are antisymmetric sign factors with 12 = − 12 = 1.

JHEP02(2018)180
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 → b¯ ν andt →b ν amplitudes as shown in eq. (A.10) and eq. (A.29), 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][26] adopted by a Feynman amplitude calculator MadGraph [27]. The amplitudes eq. (A.10) and eq. (A.29) can be recovered by supplying the phase factor in eq. (B.7), as will be explained in detail in appendix B. Because the phase factor associated with the¯ and spinors, respectively, for t andt decay amplitudes are common for both helicity amplitudes, the decay density matrices in eqs.

B HELAS phase convention for the density matrix
As explained above, we use HELAS amplitudes as obtained by MadGraph to test our analytic calculations, and to generate differential cross sections with full decay correlations. We encounter a subtle frame dependence of the htt production density matrix elements due to the specific phase convention, hereafter called HZ convention [24], which has been adopted by HELAS subroutines [25,26]. Because the use of MadGraph-generated HELAS amplitudes can be a powerful tool for studying massive fermion spin correlations in various experiments, we show in this appendix the origin of the frame dependence and ways to avoid possible inconsistencies.
In ref. [24], massive fermion wave functions are expressed in terms of two helicity spinors where p = | p|, 0 ≤ θ < π so that cos θ 2 , sin θ 2 > 0 and 0 ≤ φ < 2π. Note that χ − ( p) = −iσ 2 χ + ( p) * (B.3a) and with the charge conjugation v L (p, σ) = iσ 2 u R (p, σ) * , v R (p, σ) = −iσ 2 u L (p, σ) * , (B.5) 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. (B.2), The fermion wave functions in ref. [24], adopted in HELAS [25,26]  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. (B.11) 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. (B.8) or simply stick to the frame where t andt momenta have no y-components. We adopt in section 2 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.