Probing CP violation in $e^{+}e^{-}$ production of the Higgs boson and toponia

We study the CP violation in the Higgs boson and toponia production process at the ILC where the toponia are produced near the threshold. With the approximation that the production vertex of the Higgs boson and toponia is contact, and neglecting the P-wave toponia, we analytically calculated the density matrix for the production and decay of the toponia. Under these assumptions, the production spectrum of the toponia is solely determined by the spin quantum number, therefore the toponia can be either singlet or triplet. We find that the production rate of the singlet toponium is highly suppressed, and behaves just like the production of a P-wave toponia. In the case of the triplet toponium, three completely independent CP observables, namely azimuthal angles of lepton and anti-lepton in the toponium rest-frame as well as their sum, are predicted based on our analytical results, and checked by using the tree-level event generator. The non-trivial correlations come from the longitudinal-transverse interferences for the azimuthal angles of leptons, and the transverse-transverse interference for their sum. These three observables are well defined at the ILC, where the rest frame of the toponium can be reconstructed directly. Furthermore, the QCD-strong corrections, which are important near the threshold region, are also studied with the approximation of spin-independent QCD-Coulomb potential. While the total cross section is enhanced, the spin correlations predicted in this paper are not affected.

Precise measurements of various physical properties of the observed Higgs particle h(125) [1,2] are the most important and urgent tasks in the elementary particle physics. Of particular interest is the property under the charge conjugation and the parity transformation, which is called the CP property. In general, the mass eigenstate h(125) can be either CP eigenstate or a mixture of CP-even and CP-odd scalar particles. While only one CP-even scalar particle is predicted in the standard model (SM), many of its extensions not only modify the Higgs couplings to gauge bosons and fermions, but also predict additional scalar and pseudo-scalar particles. Therefore decisive measurement of the CP property of h(125) can tell directly whether the observed boson is the Higgs boson in the SM, or it is described by the model beyond the SM. The CP property of h(125) has been investigated experimentally by both ATLAS and CMS collaborations [3][4][5] through the decays into vector boson pairs, and the experimental results disfavor the pure CP-odd hypothesis by nearly 3σ. However a large CP mixing has not been excluded yet . The reasons are twofold: first, the CP-even coupling of Higgs to the Z boson pair appears at the tree level while the CP-odd coupling appears only at the loop level; second, the branching ratio to the ZZ is small.
Theoretically, the CP property of h(125) can also be measured by studying the spin correlations of the two jets in the pp → hjj process [20], and the spin correlation in the h → τ τ channel [21][22][23][24][25][26][27][28][29]. For the process pp → hjj, the QCD backgrounds can significantly reduce the signal significance [20], even through the jet-matching technique can be useful to select out the signals [30]. On the other hand, the CP properties of the Higgs boson can be investigated by using the Higgs coupling to top quarks which is the largest Yukawa coupling in the SM, y t ∼ O (1). The ATLAS group have studied the tth production with an integrated luminosity of 20.3fb −1 , and set a 95%C.L. limit on the cross section σ tth < 4.1σ SM tth . However, the information on the CP properties of the top-Higgs Yukawa coupling is still lacking. In Ref. [6], the authors shown that the CP mixing parameter is limited in the range ξ htt < 0.6π. In Refs. [7][8][9] constraints on the CP-odd htt coupling is studied by using the LHC run-I data through the hgg and hγγ couplings. These constraints are not strong, and still allowing a wide range of the CP-mixing angles. In Ref. [10], a strong constraint on the CP-odd htt coupling is derived by using the constraints on electric dipole moments for several nucleus. However, this constraint is obtained under the assumption that the CP-odd htt coupling is only the source of CP violation, which means there is no contribution from heavier Higgs bosons, sparticle, electron-Higgs CP-odd couplings, etc. . If there are other sources of CP violation and there is a cancelation between them, the constraint can be weakened.
As well, there have also appeared many papers devoted to find optimized CP observables at hadron colliders [31][32][33][34][35][36][37][38][39] and lepton colliders [40][41][42][43][44]. The simplest one requires the reconstruction of the top-and anti-top-quark momenta from their decay products which is difficult to be measured accurately even at lepton colliders. In principle, one can construct CP-odd observables by replacing the top-and anti-top-quark momenta by the momenta of the b andb-jets from the t andt decays, respectively. However, the sensitivity to the CP violating effects gets diluted in this partial reconstruction. It has also been pointed out that the different phase-space distributions for scalar and pseudo-scalar Higgs boson production rates can be used to determine the CP properties of the tth coupling. In Ref. [40], the authors have demonstrated that the CP properties of Higgs can be assessed by measuring just the total cross section and the top-quark polarization. However, these two observables are CP-even, hence only proportional to the square of the CP-odd coupling. Furthermore, the ratio of the production rates for pseudo-scalar and for scalar is very small unless √ s 1 TeV where the chiral limit is recovered. Therefore, the experimental sensitivities of these observables are not as good as enough to probe small CP-odd coupling. To pin down the CP property of the Higgs boson, true CP-odd observables, which is linearly proportional to the CP-odd coupling are really required. The up-down asymmetry of the momentum direction of the anti-top quark with respect to the top-quark-electron plane is an example of such an observable [41,42]. However, the asymmetry is due to the interferences between the amplitudes involving the tth vertex and those involving the hZZ vertex. It has been shown that the latter contribution is very small, amounting to only a few percent for √ s ≤ 1 TeV [40]. Therefore only about 5% asymmetry can be observed at the largest [41,42].
In this paper, we study the density matrix for the e + e − production of the Higgs boson and toponia analytically, and propose new CP-odd observables for the measurement of the CP property of the Higgs boson at the ILC with √ s = 500 GeV [43][44][45][46][47]. In this energy region, the strong-interaction Coulomb force is known to be important to calculate the total production rate. Because the P-wave tt 1 production is heavily suppressed, we focus on the S-wave toponia production. The contents of this paper is organized as follows. In Sec. II we discuss the effective tt production vertex and the spectrum of the toponia in the e + e − → tth process. In Sec. III we present the helicity amplitudes for the S-wave toponia productions and their decays. In Sec. IV we study the QCD bound-state effects for the tt system. In Sec. V we give the numerical results based on the tree-level event generator, and discuss the CP asymmetries from leptonic observables. Finally, discussions and conclusions are given in Sec. VI.

II. EFFECTIVE t −t − h VERTEX
In this section we study how the tth interactions affect the tt system production near the threshold. We assume that the observed Higgs particle h(125) is a mixture of CP-even (H) and CP-odd (A) particles, where ξ is the Higgs mixing angle which is assumed to be real. For simplicity, we further assume that the Yukawa interactions are CP conserving, such that the source of CP violation is only in the Higgs mixing angle ξ in Eq. (1). The interactions between the mass eigenstate h(125) and the fermion anti-fermion pair are then described by where It is worth noting that the effective strengths of the CP-violating hf f couplings can be different for each fermion, even if the origin of CP-violation is only in the mixing parameter ξ. In this paper, we focus on the htt coupling, and for convenience we use the symbol g h to denote the overall coupling constant g htt , i.e. g h = g htt . The assumption of the real mixing parameter is valid when CP violation in the Higgs sector is mediated mainly by the interactions with new heavy particles. For the s-channel production of tt associated with h(125), h(125) can be emitted from either a very virtual top-quark or anti-top quark as shown in Fig. 1. Even through the Higgs boson can also be produced through the hBB vertexes (B = Z, γ), the contributions are negligible (a few percent for √ s ≤ 1 TeV [40]) because of the far off-shell propagator of the vector bosons. In principle, CP violation can also appear in these vetrices. However CP violating operators are induced at the one-loop level, and hence hugely suppressed compared to the CP-even operators. Therefore, we do not consider them to simplify the vertex function in this section.
Near the production threshold at √ s thr. = 2m t + m h 471 GeV, the tth system is nonrelativistic. According to the uncertainty principle, the virtual top and anti-top quarks can propagate only in a distance ∼ 1/( √ s−m t ), which is considerably shorter than the Coulomb radius r C ∼ 1/(α s m t ), at which the QCD interactions bound top and anti-top quarks to form the bound states toponia. Therefore treating the whole production vertex as a local interaction should be a good approximation near the threshold. By denoting the vertex of tt production from a virtual vector boson B (B = γ, Z) as Γ µ B = g Btt V γ µ + g Btt A γ µ γ 5 , the leading order effective Higgs radiation vertex is given as where Γ h is the abbreviation of the tth vertex which is Γ h = g h for the pure scalar case and Γ h = g h tan ξ htt γ 5 for the pure pseudo-scalar case, and the kinematical variables are defined as in Fig. 2 with Q = k 1 + k 2 = p 1 + p 2 + k. Because both tt and h(125) are non-relativistic, the 3-momenta p 1,2 could be neglected in the denominators i.e. p µ 1,2 ≈ (m t , 0). Then the two radiation channels can be combined into a compact form. For convenience, we expand the spinor structure of this vertex by using the Clifford algebra as follows: where we have used Q 2 = s. The expansion coefficients can be calculated easily, as shown in Table I. The production dynamics are described completely by the vertex function V µ (p 1 , p 2 ) in Eq. (7). Note that the coefficients of the (CP-even) hBB vertexes are not included in Table I  Definitions of the kinematical variables in the e + e − rest frame specified by the axes xy-z, and the tt rest-frame specified by the axes x -y -z . In the e + e − rest-frame, the electron momentum is chosen along the z-axis and the tt momentum lies in the x-z plane with positive x-component. In the tt rest-frame, the h momentum direction is chosen as the opposite of the z -axis, and the y -axis is taken as the same direction as the y-axis. and c µν A . Furthermore, the spin correlation which can be used to measure the CP violation effects does not depend on the coefficients of these operators. The magnitudes of these contributions are discussed in the numerical simulation part in Sec. V.
After the electroweak production of tth, the strong interaction between tt becomes important. In the threshold region, infinite number of Feynman diagrams whose effects are proportional to the powers of α s /β t ∼ O(1) contribute, and their resummation is needed; see Fig. 3.
After the resummation, the vertex function satisfies an integral equation, the Salpeter-Bethe equation [48], which describes the formation of bound states in this region. We will discuss it carefully in Sec. IV. Here we would like to classify the possible bound states that can be produced. Table II lists the possible bound states up to P-wave in the spectrum notation for various bi-spinor combinations of spinors ψ and ϕ (see App. VII A for our conventions of the spinor wave functions in the Dirac representation), and the corresponding spinor vertex structures in the non-relativistic limit. The spin-singlet state can be produced only by the pseudoscalar operator O P and the time component of the axial-vector operator O µ=0 V . All the other operators can generate the spin-triplet state but with different angular momentum.
It should be noted that, all those quantum numbers listed in Table II are also affected by The htt vertex is denoted as Γ h = g h + ig h tan ξ htt γ 5 . The momentum q µ = p µ 1 − p µ 2 is the relative momentum between the top and anti-top quarks. Note that the coefficients of the (CP-even) hBB vertexes are not included for the clarity and compactness of the table. the corresponding expansion coefficients, which are tabled in Table I. In Table III, we show the possible bound states by combining the coefficients and operators. For the scalar operator O S , both the coefficient and bi-spinor are of P-wave for the scalar Higgs boson. Therefore the tt system is D-wave which can be ignored completely. In the case of the pseudo-scalar Higgs boson, the tt system is P-wave because the coefficient is S-wave. However, it is still negligible near the threshold region. For the pseudo-scalar operator O P , a singlet toponium can be produced. The coefficient is S-wave for the scalar Higgs boson, while P-wave for the pseudo-scalar Higgs boson. For the vector and axial-vector operators, O V and O A , only vertexes for scalar Higgs boson production exist. The operator O V can generate the S-wave triplet toponium, while O A generates the P-wave triplet toponium. In addition, the axialvector operator can also generate the S-wave singlet toponium via its time-component. This TABLE II. Quantum numbers of the bi-spinors of top and anti-top quarks in the non-relativistic limit in the rest frame of tt.
Operators Non-relativistic limit Quantum state contribution turns out to be very important, because it is destructive with the contribution of the pseudo-scalar vertex O P , and then makes the total production rate of the singlet toponium highly suppressed. Of particular interest is the production by the tensor operator O T , in which both the bi-spinor and the coefficient contain S-wave and P-wave tt. Here we discuss only the S-wave contributions. For both scalar and pseudo-scalar Higgs bosons, it is the "electric component" of the tensor operator ∝ σ 0i generating the S-wave toponium.

III. HELICITY AMPLITUDES
In this section we give a formula for the full helicity amplitudes in terms of the toponium angular momentum. Near the threshold the QCD-strong interactions become important. Here we assume the QCD corrections can be completely factorized out, i.e. the strong force is spin-independent; see Sec. IV. In this approximation the full physics could be modeled by using pure electroweak htt production and their decays. Then, the toponium helicity is obtained by the spin projection. The spin projection becomes simple when the relative momentum q µ between the top and anti-top quarks is neglected. Furthermore neglecting the relative momentum does not lose essential physics as the top and anti-top quarks have large decay width. Therefore, while we calculate the density matrix without the assumption of |q µ | ≈ 0, some important results can be discussed under this simplification. In subsection III A we give our formalism on the factorization of the QCD correction, as well as that for the spin projection. In subsection III B and III C we give the helicity amplitudes for the production and decays of toponia. The total helicity amplitude and the CP-odd observables are discussed in subsection III D.
TABLE III. Quantum states of the tt and tth systems. The spin and angular momenta are summed first by combining the top and anti-top-quarks system, and then by combining the toponium (ψ t ) and Higgs system.

A. Factorization and projection of the helicity amplitudes
The total amplitude for the process e − +e + → h+( ν¯ b )+(¯ ν b) can be written in general as follows: We focus on the CP violation effects due to the anomalous interaction between the toponia and Higgs boson. This is done by inserting a complete basis of the tt resonance states ψ t with quantum number (J ψt , λ ψt ), then the total helicity amplitude can be written as the product of the production and decay amplitudes of the toponia, 2 However this amplitude cannot be calculated directly in perturbation theory because the tt resonances ψ t are composite states. We therefore expand the helicity amplitudes by using the fundamental fields t andt, and the amplitudes take the following form: where the production, decay and resonance amplitudes are, respectively, Here both the production and decay processes are electroweak, and the QCD corrections are accounted for in the resonance amplitudes. In order to make our discussions more simple and clear, we use the free tt resonance states ψ t (J ψ , λ ψt ) to separate out the spin degrees of freedom. Then the amplitude for the toponium formation from the top-and anti-top-quarks can be written as where we have introduced a pure kinematical operator O J ψ λ ψ t to account for the spin correlations of tt to ψ t . In general the quantum numbers (J ψt , λ ψt ) can be different from (J ψt , λ ψt ) by QCD corrections, for instance when we include the spin-orbital interactions, etc. Here we neglect those spin-dependent corrections, i.e. we take (J ψt , λ ψt ) = (J ψt , λ ψt ). Then the resonance amplitudes can be written as where the factor K J ψ t ,λ ψ t is defined as the squared renormalization factor which gives the pure QCD corrections, and the spin projection operator P The QCD corrections are discussed in Sec. IV. Let us focus on the spin projection first.
In general J ψt can be any integer. However the production rates of toponium states with higher angular momentum L are suppressed by β L t where β t is a velocity of top and anti-top quarks in the toponium rest-frame. Therefore, we discuss only the S-wave resonance. Then ψ t can be either spin-singlet or spin-triplet, i.e. J ψt = 0 or 1. The corresponding projection operators are defined as follows: where √ s ψt , m t andm t are the invariant mass of the toponium, top and anti-top quarks respectively. The normalization factor is chosen such that the spin projection operators are dimensionless (the overall normalization of M ψt is fixed by the total QCD correction). With the help of the spin projection operators the total helicity amplitude can be expressed in terms of the toponium production and decay helicity amplitudes as follows: where the projected production and decay helicity amplitudes are In the next two subsections, we study these two helicity amplitudes.

B. Production helicity amplitudes
In this subsection we give the helicity amplitudes for the production process of toponia in associated with the Higgs boson. The kinematical variables are defined as (see also the The fermion helicities are σ i = ±1/2 for i = e,ē, t,t. For the spin-singlet toponium J ψt = 0, λ ψt = 0, and for the spin-triplet toponium J ψt = 1, λ ψt = 0, ±1. In the rest frame of e + e − the particle momenta are given by Here we use √ s to denote the total collision energy, and √ s ψt to denote the invariant mass of the toponium. β is a velocity of the Higgs boson and toponium in this frame, which is given as In this frame the leptonic current is give by where ε µ ( Q = 0, λ V ) are given in Eq. (95) and Eq. (96) in the Appendix VII B by setting θ = 0 and φ = 0, λ V = σ e − σē = ±1 is the helicity of the virtual vector particle B that can be either photon (B = γ) or Z (B = Z); the helicity-dependent form-factor G e λ V is defined as where the first term stands for the photon pole and the second term stands for the Z pole. The momenta of the toponium, t andt in the rest frame of the toponium are given by where Let us first calculate the projection operators. Because we discuss only the S-wave toponium, there are only two kinds of projection operators: the spin-singlet and spin-triplet projection operators which correspond to the matrix element of operators O P and O V . In the rest frame of the toponium we get where the helicities m = σ t − σt andm = σ t + σt are defined along the top-quark momentum direction, and they are related by the Wigner rotation to the helicity states of the toponium along its moving direction. The function f (m, m) is defined as follows: Here we use D to denote the complex-conjugate-transpose of the Wigner-D functions; see Appendix. As we have worked in the non-relativistic approximation, the relative momentum between top and anti-top quarks is negligible, so the kinematical factor β t in the spin-triplet projection operator can be neglected. The helicity amplitudes of tth production are decomposed by the type of production vertexes. Here we use the notation M P (X; σ t , σt) with X = S, P, A, V, T to denote their contributions, and use subscripts of X to distinguish the contributions of scalar and pseudoscalar components of the Higgs boson. The operators that can generate the toponium in S-wave are listed in Table IV.
For the scalar operator, both the scalar and pseudo-scalar components of the Higgs boson start to contribute at P-wave, so there is no relevant contributions. For the pseudo-scalar operator, only the scalar component of the Higgs boson contributes, and the helicity amplitude is where the kinematical factor which is consistent with our previous explanation in Sec. II that the pseudo-scalar operator can only generate the P-wave state of the singlet toponium and Higgs boson. This is also true for the axial-vector operator. The helicity amplitude is similar with M P (P ; σ t , σt), where the kinematical factor is The important thing is that the contributions of pseudo-scalar and axial vector operators are destructive. Because only the pseudo-scalar and axial vector operators generate the singlet toponium, therefore the total helicity amplitude for the singlet toponium production is just the sum of these two contributions. It is proportional to 1 − 4m 2 t /s ψt , and thus negligible near the threshold.
The triplet toponium can be produced through the vector and tensor operators. The helicity amplitude for the vector operator is Here the helicity λ ψt is quantized along the moving direction of the toponium in the e + e − rest-frame, and related to λ ψt by the Wigner rotations after spin projection. The kinematical factor is This is a S-wave production, and can be represented by an effective operator hψ µ t B µ , where B = γ, Z. The contribution from the tenser operator is also of S-wave production, and can be represented by an effective operator hF ψtµν F µν B , where F µν B = ∂ µ B ν − ∂ ν B µ is the field strength tensor. The corresponding helicity amplitude is where the kinematical factor is In the above calculations we have neglected a contribution of the D-wave production which is proportional to β 2 . Apart from the kinematical factor, the rest is completely the same as the contribution of the vector operator M P (V ; σ t , σt). These two contributions are constructive, and hence make the triplet production rate dominant. Furthermore, the pseudo-scalar component of the Higgs boson also contributes in the S-wave toponium production via the tensor operator. However, the overall production is of P-wave. The corresponding effective operator can be written as h F ψtµν F µν B , where F µν B = 1/2 µναβ F αβ B is the dual strength tensor of the field B. The helicity amplitude for the tensor operator is, where the kinematical factor is and h = g h tan ξ htt . Now we can obtain the projected helicity amplitudes. For the pseudo-scalar and axial vector operators the projected helicity amplitudes are similar with each other; Because only these two operators contribute to the singlet toponium production, the total helicity amplitude for the singlet toponium production is given as As expected this is the usual production helicity amplitude of two scalar particles in P-wave. Because it is strongly suppressed by the kinematical factor 1 − 4m 2 t /s ψt which vanishes near the threshold. Therefore we will neglect the singlet toponium in the following study of spin correlations.
For the triplet toponium production, the vector and tensor operators contributes. Apart from the kinematical factors, the projected helicity amplitudes have also the same structure, and proportional to the Wigner-D function as follows: Here we have used a relation σt,σt This is a usual production helicity amplitude of a vector particle. Because the structures of the helicity amplitudes for these three operators are the same, we can add them up directly. After the summation, the helicity amplitudes are given by where the kinematically suppressed CP phaseξ htt and the suppression factor are defined as follows: The production density matrix is defined as Inserting the helicity amplitudes we get
As we have mentioned the helicity amplitudes of the toponium decay are obtained by using the spin projection of the helicity amplitudes of the tt decay. The helicity amplitudes of the tt decay can be separated into the t andt decay amplitudes as The helicity amplitudes of the t andt decays have been known for long time. In the rest frame of the toponium the kinematical variables are defined as follows: k µ 2 = E (1, sin θ cos φ , sin θ sin φ , cos θ ) .
Here and after we neglect the lepton mass. By using the Fierzt transformation, the kinematical variables of (bν)/(bν ) can be factorized out completely. Thus the anti-lepton and lepton carry all the spin informations of t andt, respectively. Then the helicity amplitudes of t andt decays can be written as, where A t and At which are Lorentz invariant stand for the rest of the helicity amplitude. The tt decay helicity amplitudes M D (σ t , σt) can be obtained by using Eq. (52). In terms of m andm, M D (σ t , σt) can be written as follows: where the functions f m,m and k m,m are defined as follows: g m,m (θ ¯ , θ ) = 1 + m cos θ ¯ 1 + m cos θ .
The projected helicity amplitudes can be obtained by using the projection operators in Eq. (30a) and Eq. (30b). As we have explained above, the production rate of the singlet toponium is highly suppressed near the threshold region. Therefore we give only the decay amplitudes for the triplet toponium. By using Eq. (22), the projected decay amplitudes for the triplet toponium are explicitely written as The decay density matrix is defined as where we have integrated out the phase spaces of the bottom quarks and neutrinos. The corresponding matrix elements are The spin correlations occur if the imaginary part of the decay density matrix is non-zero. The above results indicate that the spin correlations can appear in both the transverse-transverse and transverse-longitudinal interferences.

D. Total helicity amplitudes and CP-odd observables
In this subsection we discuss the interferences among the different helicity states of the triplet toponium. The CP-odd observables are obtained by studying the spin correlations due to the interferences. As we have mentioned there are two kinds of interference: transversetransverse (TT) and longitudinal-transverse (LT) interferences, which are predicted by the total density matrix, where for convenience we have defined an intermediate density matrix as follows: For the TT interference we have (for convenience we useξ to denote the variableξ htt for abbreviation), Therefore the production rate has a following non-trivial distribution with respect to the observable φ ¯ + φ , where σ 0 is the total cross section, and C T T is the coefficient for the T T correlation. For the LT interference we have We can see that the azimuthal-angle distributions of lepton and anti-lepton have differentξ dependence. The lepton momentum has a following non-trivial distribution, where C LT is the coefficient of the LT correlation. For the anti-lepton momentum, we have We can see that the correlations are different for the lepton and anti-lepton. For the lepton, the correlation is positive, while negative for the anti-lepton. On the other hand, the sign and the size of the phase shift is the same for both the lepton and anti-lepton. These two distributions are related with each other by the CP transformation. In the case ofξ = 0, i.e. CP is conserved, these two correlations are symmetric under the CP transformation φ ¯ → π − φ and likewise for φ . However, ifξ = 0, the distributions are asymmetric bỹ ξ → −ξ, and therefore indicates the violation of the CP symmetry.

IV. RADIATIVE CORRECTIONS NEAR THE THRESHOLD REGION
As we have explained in Sec. II, the virtual top or anti-top quark is hugely off-shell. According to the uncertainty principle, it can propagate only a distance ∼ 1/( √ s − m t ) which is considerably shorter than the Coulomb radius r C ∼ 1/(α s m t ) for the tt boundstate. Therefore, near threshold production can be treated by a local source δ 4 (y t − yt)j µ (Q 2 )e −iQ·yt . In this approximation, the Higgs field decouples from the exact vertex function Th(z )t i (y t )t j (yt)V µ (z) by modifying the ttV vertex function which has been examined in Sec. II. The modified production vertexes are then in turn to affect the quantum numbers of the generated toponia, which have been discussed in Sec. II. Here we examine how these vertexes are affected by the QCD radiative corrections. The corrections are described by the relativistic Salpeter-Bethe (SB) equation in general [48]. For a general production vertex Γ µ C (the subscript "C" is used to distinguish possible different Dirac matrix), the SB equation is where U αβ (q − k) is the QCD potential in momentum space, and S F is the Feynman propagator for fermions. This integral equation sums over all the contributions from the relevant ladder diagrams; see Fig. 3. Here we consider only the instantaneous Coulomb-like potential, the contributions from the transverse and rest gluons are suppressed by powers of β t . In the rest frame of tt, the dominant contributions come from the region where | k| m t , and the fermionic propagators are approximated by where γ ± = (1 ± γ 0 )/2 are the non-relativistic projection operators for fermion and antifermion. Observing that the vertex function is independent of the energy q 0 , the variable k 0 can be integrated out and we get In our case the toponium system can be boosted by the recoil of the Higgs boson, therefore we express all the quantities in a Lorentz-invariant manner as follows: The integral equation can be rewritten in a covariant form as We define the dressed non-relativistic projection operators for fermions and anti-fermions as follows: The second terms in both γ + ( q) and γ − ( q) involve the small component of the Dirac spinor which are of P-wave, and therefore suppressed by an additional factor of β t . Therefore in the following calculations we can neglect them. In this approximation, a useful relation can be derived as follows: Multiplying γ + ( q) on the left-hand side and γ − ( q) on the right-hand side of Eq. (75), we get Introducing the non-relativistic reduced vertex function the integral equation Eq. (78) reduces to This is a formal Lippmann-Schwinger (LS) equation [49]. Here we study only the corrections to the production vertex up to terms linear in q. Expanding the vertex Γ µ C (E, q) by q we have, where are the S-and P-wave components, respectively. The corrected vertex function V µ C (E, q) can be expanded in the same way, and we get The expansion coefficients satisfy following integral equations These two integral equations are related to the Green function G( r x , r x ) which satisfies the LS equation in the momentum space as follows: As we have mentioned, the local interaction approximation is excellent in the production vertex, therefore the vertex functions are approximated by the condition r y = 0. Expanding the Green function G(E; k, r y ) by r y as and the plane wave factor e i p· ry ≈ 1 + i p · r y we obtain the following integral equations, The solutions of the above equations are formally written as where G 0 (E; p) is the Green function of a free toponium, These Green functions are related to the correction factors K S and K P as follows: In this study, we employ the method give in Ref. [50] to numerically solve the integral equation. Fig. 5 shows the S-and P-wave Green functions for the binding energy E ≡ √ s ψt − 2m t = [−2, 0, 2, 4] GeV. We can see that at the ground state (E −2 GeV), the Pwave contribution is suppressed. However, the corrections on S-and P-wave are comparable for other states. Fig. 6 shows the counter lines of the absolute values of the Green functions in the plane of the binding energy E and the relative momentum | q|.

V. NUMERICAL RESULTS
Our numerical results are obtained by using MadGraph5 [51] with the HC model [52] at the tree level, and then weighted by the QCD correction factor K S/P (E, q) at the LO. For a smooth connection to the large M tt region, we follow the prescription described in Ref. [53]. Left and right panels in Fig. 7 show the production cross sections at √ s = 500 GeV with respect to the invariant mass of the tt system for the pure scalar and pure pseudo-scalar cases, respectively. We can see that the distribution of total production cross sections have a peak below the threshold energy. At this collision energy, the LO cross section is calculated to be σ LO = 0.29 fb for the pure scalar case (we assume that the electron and position beams are not polarized). Note that when we include the diagram with a hZZ vertex, the total cross section is enhanced by about 1.7% due to the interference with the diagrams which contain the top-Yukawa vertex. The QCD-Coulomb corrections give an enhancement factor of about 2.6. However, it has been pointed out that the NLO corrections are important particularly in the large tt invariant mass region [54,55], which gives an additional overall correction factor of about K = 0.84 [56]. In total, our prediction to the total cross section is about 0.63 fb. We use this total cross section for the overall normalization. On the other hand, the NLO effects are almost uniform in the whole phase space, therefore our LO estimation can be safely used for studying the spin correlations.
With the approximation of the S-wave dominance, we have calculated the azimuthal angle correlations of leptons from the decays of top and anti-top quarks. We have shown that there are three independent CP-odd observables. The first one is the sum of the azimuthal angles of leptons in the toponium rest-frame, which is due to the interference among the transverse components of the triplet toponium. The correlation function has been given in Eq. (66). Fig. 8 (a) shows the correlations for pure scalar Higgs (blacksolid line) and for pure pseudo-scalar Higgs (red-dashed line). Both are symmetric under the sign reflection of φ ¯ + φ , because of the CP conservation separately. However the shapes are completely opposite. In the case of scalar Higgs boson, the interference are constructive when the sum of azimuthal angles is either π or −π. However it is constructive when the sum is 0 for a pseudo-scalar Higgs boson. Therefore the CP violation effect is sensitive to the sign of the mixing angle ξ htt . Fig. 8 (b) show three different CP-mixed cases: tan ξ htt = 0 (black-solid), tan ξ htt = 5 (red-dashed) and tan ξ htt = −5 (blue-dotted).
Here in order to show the differences clearly we have chosen | tan ξ htt | = 5 which means an effectively maximum mixing because of the kinematical suppression factor κ ≈ 0.2, see Eq. (48a) and Eq. (48b). Measuring the CP violation effects from transverse-transverse interference requires the reconstructions of both lepton and anti-lepton in the topponium rest-frame. The branching ratio of top quark to leptons (e, µ) is Br(t → X) = 19%. If we use the h → bb channel to reconstruct the Higgs boson, which has a branching ratio 56.9%, there are 52 signal events with 100% reconstruction efficiency for the projected integrated luminosity 4 ab −1 at √ s = 500 GeV [45]. Simple estimation on the experimental sensitivity is δξ htt = 1.34. However, because of ξ htt =ξ htt /κ with κ 0.2 for √ s = 500 GeV, the accuracy of constraining the non-zero CP-phase may be limited at the ILC with the nominal luminosity. This low sensitivity comes from 1) the low total production rate, and 2) a large kinematical suppression factor κ in Eq. (48b).
Apart from the interference among the transverse polarizations, there are also interferences between the longitudinally and transversely polarized toponia which result in nontrivial azimuthal angle distributions of the lepton and anti-lepton individually in the toponium rest-frame. The correlation functions are given in Eq. (69) and Eq. (70). It is constructive at the origin (φ ,¯ = 0) for the lepton, while destructive for the anti-lepton. For the pure scalar case, this feature is shown in Fig. 9 (a). For the pure pseudo-scalar case, because only the transversely polarized toponium can be produced, there are no interference between longitudinally and transversely polarized states. Therefore the azimuthal angle distribution is flat, which is shown in Fig. 9 (b). Fig. 9 (c) and (d) show the interferences in the three cases:ξ htt = 0 (black-solid),ξ htt = π/4 (red-dashed) andξ htt = −π/4 (blue-dotted) for the lepton and anti-lepton, respectively. We can see that both the lepton and anti-lepton azimuthal distributions are sensitive to the sign of the mixing angleξ htt . Most importantly, measuring CP violation effects through the transverse-longitudinal interferences requires only either of the lepton or anti-lepton momentum being reconstructed. For the h → bb decay channel, there expects about 275 signal events for the projected integrated luminosity of 4 ab −1 at √ s = 500 GeV [45] (for either lepton or anti-lepton) assuming the 100% reconstruction efficiency. Combining the lepton and anti-lepton channels we expect about 550 signal events in total. For this situation, the experimental sensitivity of determining ξ htt is estimated to be δξ htt = 0.4. Taking into account the kinematical suppression factor of κ 0.2, the accuracy of determining ξ htt is estimated to be δξ htt 1.1 for ξ htt 0.

VI. DISCUSSION AND CONCLUSION
In summary, we have studied the CP violation effects in the toponia productions in association with a Higgs boson at the ILC with √ s = 500 GeV. The Higgs boson can be produced by the emissions of the top or anti-top quarks via the Yukawa interaction, or through the gauge interactions between Higgs and vector bosons, Z or γ. The CP violation effects can appear both in the Yukawa and gauge interactions. However observing the effects induced by the gauge interactions is difficult, because the CP-odd hZZ interaction is induced at the loop level while the CP-even interaction appears at the tree level. Hence the CP asymmetry induced by the hZZ coupling is suppressed by a factor of α W /(4π). In addition, for the e + e − production at √ s = 500 GeV, the dominant contributions stem from the Higgs emissions from top quarks, but the contributions from the gauge interactions can reach only up to a few percent. Therefore, CP violation effects which emerge by the top-Yukawa couplings should be thoroughly explored.
For √ s = 500 GeV, the produced toponia are non-relativistic, therefore the production of the P-wave toponium is negligible. The eligibility of this assumption is confirmed by the numerical calculation based on the tree-level event-generator; see Fig. 7. Furthermore, the htt production vertex from a virtual vector boson Z or γ can be modeled by a contact vertex operator. By assuming that the spins of the top and anti-top quarks are not altered by the QCD potential, i.e. the QCD potential is spin-independent, the produced toponia spectrum are studied carefully. In this approximation, the relevant toponia are the spin-singlet 1 S 0 and the spin-triplet 3 S 1 states. However, the production rate of the singlet toponium is found to be highly suppressed, because it is P-wave in the toponium and Higgs boson system. This observation has been checked by using the tree-level event-generator.
Based on the careful analysis for the helicity amplitudes of the production and decay of toponia, we propose three CP-odd observables, namely the phase-shift in the azimuthal angle of the lepton and anti-lepton as well as their sum. These observables are induced by the non-trivial correlations in the longitudinal-transverse interferences in the azimuthal angle distributions of leptons, and in the transverse-transverse interference in their sum. Compared to the up-down asymmetry examined in Refs. [40][41][42] which requires the reconstruction of either the top-or anti-top-qaurk momenta, as well as the small contribution from the diagram which contains hZZ interactions (a few percent for √ s ≤ 1 TeV [40]), our observables do not require the reconstruction of the momentum of the top or anti-top quark individually, and are caused purely by the dominant htt interactions. Furthermore, all the three observables have maximum asymmetries of about 32%, which are more than 6 times larger than the maximum asymmetry (5%) in Refs. [41,42]. Because the CP-odd observables for the longitudinal-transverse interferences can be reconstructed by using only the one lepton momentum, the number of signal events can be increased. The experimental sensitivities for these observables are estimated for an integrated luminosity of L = 4 ab −1 , and found to be δξ htt 1.1 for ξ htt = 0. Since the sensitivity is limited mainly due to the statistical fluctuation, it can be improved by increasing the luminosity as projected in Ref. [45].
Compared to the current constraints on ξ htt by the LHC measurement, which has set ξ htt < 0.6π [6], and further improvements by future LHC measurements, the sensitivities of our observables may be relatively low, δξ htt 1.1 for ξ htt = 0. However, our observables can be used to directly measure the CP phase, rather than to measure the overall rates. Particularly, our observables φ and φ¯ require either top or anti-top decaying to leptons, and therefore the efficiency would be enhanced dramatically.

A. Spinor wave functions in the Dirac representation
For completeness we give our conventions for the spinor wave functions in the Dirac representation. In the Dirac representation, Dirac matrixes are given as follows: The free solutions of the Dirac equation in the Dirac representation are where ξ s and η r are eigenstates of the helicity operators σ · p 1 /| p 1 | and σ · p 2 /| p 2 |, respectively. For completeness we also give the helicity eigenstates as follows: ξ + = cos(θ/2) e iφ sin(θ/2) , ξ − = −e −iφ sin(θ/2) cos(θ/2) .
The spinor wave functions and the Dirac gamma matrices in the Dirac representation are related to the ones in the chiral representation by the following unitary transformation: