# Di-Higgs production as a probe of flavor changing neutral Yukawa couplings

Shi-Ping He<EMAIL_ADDRESS>Center for Future High Energy Physics and
Theoretical Physics Division, Institute of High Energy Physics, Chinese
Academy of Sciences, Beijing 100049, China

###### Abstract

Top partners are well motivated in many new physics models. Usually, vector
like quarks $T_{L,R}$ are introduced to avoid the quantum anomaly. It is
crucial to probe their interactions with the standard model particles.
However, flavor changing neutral couplings are always difficult to detect
directly in the current and future experiments. In this paper, we will show
how to constrain the flavor changing neutral Yukawa coupling $Tth$ through the
di-Higgs production indirectly. We consider the simplified model including a
pair of gauge singlet $T_{L,R}$. Under the perturbative unitarity and
experimental constraints, we choose $m_{T}=400~{}\mathrm{GeV},s_{L}=0.2$ and
$m_{T}=800~{}\mathrm{GeV},s_{L}=0.1$ as benchmark points. After the analysis
of amplitude and evaluation of the numerical cross sections, we find that the
present constraints from di-Higgs production have already surpassed the
unitarity bound because of the $(y_{L,R}^{tT})^{4}$ behavior. For the case of
$m_{T}=400~{}\mathrm{GeV}$ and $s_{L}=0.2$, $\mathrm{Re}y_{L,R}^{tT}$ and
$\mathrm{Im}y_{L,R}^{tT}$ can be bounded optimally in the range $(-0.4,0.4)$
at HL-LHC with $2\sigma$ CL. For the case of $m_{T}=800~{}\mathrm{GeV}$ and
$s_{L}=0.1$, $\mathrm{Re}y_{L,R}^{tT}$ and $\mathrm{Im}y_{L,R}^{tT}$ can be
bounded optimally in the range $(-0.5,0.5)$ at HL-LHC with $2\sigma$ CL. The
anomalous triple Higgs coupling $\delta_{hhh}$ can also affect the constraints
on $y_{L,R}^{tT}$. Finally, we find that the top quark electric dipole moment
can give stronger bounds of $y_{L,R}^{tT}$ in the off-axis regions for some
scenarios.

###### Contents

  1. I Introduction
  2. II Framework of flavor changing neutral couplings
    1. II.1 Minimal singlet vector-like quark model
    2. II.2 Simplified model
  3. III Constraints on the simplified model
  4. IV The analysis of double Higgs production
    1. IV.1 New physics results of the amplitude
    2. IV.2 Heavy quark expansion
    3. IV.3 The cross section analysis
  5. V The numerical results and constraint prospects
    1. V.1 The benchmark point $m_{T}=$ 400 GeV and $s_{L}=0.2$
    2. V.2 The benchmark point $m_{T}=$ 800 GeV and $s_{L}=0.1$
    3. V.3 Comments on the doublet and triplet vector-like quarks
  6. VI Summary and conclusions
  7. A Asymptotic behaviors of the loop functions
    1. A.1 The shorthand notations of $C_{0}$ and $D_{0}$ functions
    2. A.2 Heavy quark expansion of $C_{0}$ function
    3. A.3 Heavy quark expansion of $D_{0}$ function

##  I Introduction

The standard model (SM) of elementary particle physics has been proposed for
more than fifty years [1, *Weinberg:1967tq, *Salam:1968rm], and it is proved
to be a very effective description of this field [4]. The electro-weak
symmetry breaking (EWSB) [5, *Higgs:1964ia, *Guralnik:1964eu, *Kibble:1967sv]
mechanism predicts the existence of a physical Higgs boson, which is observed
at the Large Hadron Collider (LHC) in 2012 [9, *Chatrchyan:2012xdj]. Although
the SM goes such strong, there are still some clouds in the sky of particle
physics. The typical problems are Higgs mass naturalness, gauge coupling
unification, fermion mass hierarchy, electro-weak vacuum stability, dark
matter, matter anti-matter asymmetry, and so on. Many new physics beyond the
SM (BSM) are aimed at solving or partially solving these problems. In some of
these BSM models, vector-like quarks (VLQs)[11, 12] are introduced to avoid
the quantum anomaly.

For the up-type VLQ, the $T$ quark can interact with the SM particles via the
$TbW,TtZ$, and $Tth$ interactions. Constraints on these couplings are crucial,
because they may help us unveil the nature of the EWSB. For the strong
interaction mediated pair production of VLQs, we can only constrain the
partial decay branching ratios of the $T$ quark. To bound these couplings,
single production of VLQ needs to be considered. Unfortunately, it is
difficult to detect the flavor changing neutral (FCN) interactions $TtZ,Tth$
directly because of the suppression of the single $T$ production from $tZ,th$
fusion. Here, we will focus on the FCN Yukawa (FCNY) interaction $Tth$. After
the Higgs boson discovery, we can get more and more information on new physics
from the Higgs precision measurements [13, *Dittmaier:2012vm,
*Heinemeyer:2013tqa, *deFlorian:2016spz]. The FCNY interaction can show up in
loop induced processes, for example $h\rightarrow\gamma Z$ and $gg\rightarrow
hh$. In our previous work [17], we show that it is possible to constrain the
FCNY coupling through the $h\rightarrow\gamma Z$ decay mode indirectly. The
double Higgs production is also an appealing channel for unravelling the FCNY
interaction, which is free of electro-weak gauge interactions. The constraints
from $h\rightarrow\gamma Z$ and $gg\rightarrow hh$ are independent of exotic
decay modes and total width of the $T$ quark.

In this paper, we build the framework of FCN couplings in Sec. II first. Sec.
III is devoted to the theoretical and experimental constraints on the
simplified model. In Sec. IV, we compute the new physics contributions to the
parton level cross section of $gg\rightarrow hh$. Then we perform the
numerical constraints on the FCNY interactions in Sec. V. Finally, we give the
summary and conclusions in Sec. VI.

##  II Framework of flavor changing neutral couplings

###  II.1 Minimal singlet vector-like quark model

The SM gauge group is $SU_{C}(3)\otimes SU_{L}(2)\otimes U_{Y}(1)$, under
which the singlet up-type VLQs have the representation (3, 0, 2/3). Let us
start with the minimal extension of SM by adding a pair of singlet
$T_{L},T_{R}$ [11, 12], which is named as the VLQT model. Note that the mass
mixing term $\bar{T}_{L}t_{R}$ can be eliminated with field redefinition [18,
19]. Then, the Lagrangian can be written as [12]

$\displaystyle\mathcal{L}=\mathcal{L}_{SM}+\mathcal{L}_{T}^{Yukawa}+\mathcal{L}_{T}^{gauge},$
$\displaystyle\mathcal{L}_{T}^{Yukawa}=-\Gamma_{T}^{i}\bar{Q}_{L}^{i}\widetilde{\Phi}T_{R}-M_{T}\bar{T}_{L}T_{R}+\mathrm{h.c.},\quad\mathcal{L}_{T}^{gauge}=\bar{T}_{L}i/\\!\\!\\!\\!DT_{L}+\bar{T}_{R}i/\\!\\!\\!\\!DT_{R},$
(1)

where $\widetilde{\Phi}=i\sigma_{2}\Phi^{\ast}$ and the covariant derivative
is defined as $D_{\mu}=\partial_{\mu}-ig^{\prime}Y_{T}B_{\mu}$. $Y_{T}$ and
$Q_{T}$ are the $U_{Y}(1)$ and electric charge of the $T$ quark, respectively.
The Higgs doublet is parametrized as
$\Phi^{T}=[\phi^{+},~{}\frac{v+h+i\chi}{\sqrt{2}}]$. It is reasonable to
neglect the mixings between heavy particles and the first two generations
because of mass hierarchy and the bounds from flavor physics [20, 21, 22].
Here, we only consider the mixings between the third generation and heavy
quarks for simplicity.

To diagonalize the $t$ and $T$ quark mass terms, we can perform the
transformations

$\displaystyle\left[\begin{array}[]{c}t_{L}\\\
T_{L}\end{array}\right]\rightarrow\left[\begin{array}[]{cc}\cos\theta_{L}&\sin\theta_{L}\\\
-\sin\theta_{L}&\cos\theta_{L}\end{array}\right]\left[\begin{array}[]{c}t_{L}\\\
T_{L}\end{array}\right],\quad\left[\begin{array}[]{c}t_{R}\\\
T_{R}\end{array}\right]\rightarrow\left[\begin{array}[]{cc}\cos\theta_{R}&\sin\theta_{R}\\\
-\sin\theta_{R}&\cos\theta_{R}\end{array}\right]\left[\begin{array}[]{c}t_{R}\\\
T_{R}\end{array}\right].$ (14)

In fact, we have the relation $m_{T}\tan\theta_{R}=m_{t}\tan\theta_{L}$. In
the following, we will take $s_{L},c_{L},s_{R},c_{R}$ as shorthands for
$\sin\theta_{L},\cos\theta_{L},\sin\theta_{R},\cos\theta_{R}$, respectively.
Then, we can obtain the mass eigenstate Yukawa interactions

$\displaystyle\mathcal{L}_{Yukawa}\supset-
m_{t}\bar{t}t-m_{T}\bar{T}T-\frac{m_{t}}{v}c_{L}^{2}h\bar{t}t-\frac{m_{T}}{v}s_{L}^{2}h\bar{T}T-\frac{m_{T}}{v}s_{L}c_{L}h(\bar{t}_{L}T_{R}+\bar{T}_{R}t_{L})-\frac{m_{t}}{v}s_{L}c_{L}h(\bar{T}_{L}t_{R}+\bar{t}_{R}T_{L}),$
(15)

and gauge interactions

$\displaystyle\mathcal{L}_{gauge}\supset\frac{g}{c_{W}}Z_{\mu}[(\frac{1}{2}c_{L}^{2}-\frac{2}{3}s_{W}^{2})\bar{t}_{L}\gamma^{\mu}t_{L}+(\frac{1}{2}s_{L}^{2}-\frac{2}{3}s_{W}^{2})\bar{T}_{L}\gamma^{\mu}T_{L}+\frac{1}{2}s_{L}c_{L}(\bar{t}_{L}\gamma^{\mu}T_{L}+\bar{T}_{L}\gamma^{\mu}t_{L})$
$\displaystyle-\frac{2}{3}s_{W}^{2}\bar{t}_{R}\gamma^{\mu}t_{R}-\frac{2}{3}s_{W}^{2}\bar{T}_{R}\gamma^{\mu}T_{R}]+\frac{gc_{L}}{\sqrt{2}}(W_{\mu}^{+}\bar{t}_{L}\gamma^{\mu}b_{L}+W_{\mu}^{-}\bar{b}_{L}\gamma^{\mu}t_{L})+\frac{gs_{L}}{\sqrt{2}}(W_{\mu}^{+}\bar{T}_{L}\gamma^{\mu}b_{L}+W_{\mu}^{-}\bar{b}_{L}\gamma^{\mu}T_{L}).$
(16)

Here, we have two independent extra parameters $m_{T}$ and $\theta_{L}$. For
more details, please refer to our previous work [17].

###  II.2 Simplified model

Besides the singlet VLQs, the scalar sector can also be enlarged in the non-
minimally extended models. For example, we can also introduce a real gauge
singlet scalar [23, 24], a Higgs doublet [25], and even both the singlet-
doublet scalars at the same time [26, 25]. In these models, the $T$ quark can
exist other decay channels [27, 28]. Here, we will adopt a general framework
[29]. Then, the simplified related mass eigenstate interactions can be read as
[17]

$\displaystyle\mathcal{L}\supset-
m_{t}\bar{t}t-m_{T}\bar{T}T-eA_{\mu}\sum_{f=t,T}Q_{f}\bar{f}\gamma^{\mu}f+eZ_{\mu}[\bar{t}\gamma^{\mu}(g_{L}^{t}\omega_{-}+g_{R}^{t}\omega_{+})t+\bar{T}\gamma^{\mu}(g_{L}^{T}\omega_{-}+g_{R}^{T}\omega_{+})T$
$\displaystyle+\bar{t}\gamma^{\mu}(g_{L}^{tT}\omega_{-}+g_{R}^{tT}\omega_{+})T+\bar{T}\gamma^{\mu}(g_{L}^{tT}\omega_{-}+g_{R}^{tT}\omega_{+})t]-\frac{m_{t}}{v}h\bar{t}(\kappa_{t}+i\gamma^{5}\widetilde{\kappa}_{t})t+h\bar{T}(y_{T}+i\gamma^{5}\widetilde{y}_{T})T$
$\displaystyle+h\bar{t}(y_{L}^{tT}\omega_{-}+y_{R}^{tT}\omega_{+})T+h\bar{T}((y_{L}^{tT})^{*}\omega_{+}+(y_{R}^{tT})^{*}\omega_{-})t+\frac{gc_{L}}{\sqrt{2}}(W_{\mu}^{+}\bar{t}_{L}\gamma^{\mu}b_{L}+W_{\mu}^{-}\bar{b}_{L}\gamma^{\mu}t_{L})$
$\displaystyle+\frac{gs_{L}}{\sqrt{2}}(W_{\mu}^{+}\bar{T}_{L}\gamma^{\mu}b_{L}+W_{\mu}^{-}\bar{b}_{L}\gamma^{\mu}T_{L})-\lambda_{hhh}h^{3},$
(17)

where $\omega_{\pm}$ are the chirality projection operators
$(1\pm\gamma^{5})/2$ and the gauge couplings are listed as

$\displaystyle
g_{L}^{t}=\frac{1}{s_{W}c_{W}}(\frac{1}{2}c_{L}^{2}-\frac{2}{3}s_{W}^{2}),~{}g_{L}^{T}=\frac{1}{s_{W}c_{W}}(\frac{1}{2}s_{L}^{2}-\frac{2}{3}s_{W}^{2}),~{}g_{L}^{tT}=\frac{s_{L}c_{L}}{2s_{W}c_{W}},$
$\displaystyle\qquad\qquad~{}~{}~{}g_{R}^{t}=-\frac{2s_{W}}{c_{W}},~{}g_{R}^{T}=-\frac{2s_{W}}{c_{W}},~{}g_{R}^{tT}=0.$
(18)

The triple Higgs coupling $\lambda_{hhh}$ can deviate from the SM value
$\lambda_{hhh}^{SM}=\frac{m_{h}^{2}}{2v}$ in many new physics models [30, 31,
32, 33, 34, 35, 36, 37, 38, 39]. Here,
$m_{T},\theta_{L},\kappa_{t},\widetilde{\kappa}_{t},y_{T},\widetilde{y}_{T}$
are all real parameters, while $y_{L}^{tT},y_{R}^{tT}$ can be complex. From
now on, we will turn off the parameters $\widetilde{\kappa}_{t}$ and
$\widetilde{y}_{T}$ for simplicity.

In the following context, we will show how to constrain the FCNY couplings
through the $gg\rightarrow hh$ channel. Although the FCNY couplings
$y_{L,R}^{tT}$ are not free parameters in the above VLQT model, they can be
free in more complex models. If we can extend the SM by the singlet
$T_{L},T_{R}$ and many new scalars, there can be enough degrees of freedom.
Thus, we can take them as free parameters to make a general analysis.

##  III Constraints on the simplified model

In this section, we will review the theoretical and experimental constraints
on the simplified model. Specific details are already given in our previous
study [17]. $S$-wave unitarity will lead to the bound [17]

$\displaystyle\sqrt{(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})^{2}+12|y_{L}^{tT}|^{2}|y_{R}^{tT}|^{2}}+|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2}\leq
16\pi.$ (19)

In Fig. 1, we show the unitarity allowed region in the plane of
$|y_{L}^{tT}|-|y_{R}^{tT}|$.

Figure 1: The region allowed by the perturbative unitarity.

Higgs signal strength and top quark physics give quite loose constraints.
Direct search can bound the VLQ mass as light as 400 GeV without specific
assumptions [28]. The strongest constraints on $m_{T}$ and $s_{L}$ come from
the electro-weak precision measurements. Here, we consider the $S$ and $T$
parameters [40, 41, 42, 22, 43]. In Fig. 2, we show the allowed parameter
space regions from the global fits at $1\sigma$ and $2\sigma$ confidence level
(CL).

Figure 2: The constraints on $m_{T},s_{L}$ from the $S$ and $T$ parameters.
Here, the green and red label the allowed regions at $1\sigma$ and $2\sigma$
CL, respectively.

In this paper, the input parameters are chosen as
$m_{Z}=91.1876~{}\mathrm{GeV}$, $m_{W}=80.387~{}\mathrm{GeV}$,
$m_{h}=125.09~{}\mathrm{GeV}$, $m_{t}=172.74~{}\mathrm{GeV}$,
$G_{F}=1.1664\times 10^{-5}~{}\mathrm{GeV}^{-2}$, and $c_{W}=m_{W}/m_{Z}$ [4].

Then, we turn to the constraints from the top quark electric dipole moment
(EDM) [44, 45, 46, 47, 48]. If there exists CP violation in the FCNY
interactions, it will contribute to the EDM type interaction
$-\frac{i}{2}d_{t}^{EDM}\bar{t}\sigma^{\mu\nu}\gamma^{5}tF_{\mu\nu}$. The
$d_{t}^{EDM}$ is computed as

$\displaystyle
d_{t}^{EDM}=\frac{eQ_{T}m_{T}[y_{R}^{tT}(y_{L}^{tT})^{*}-y_{L}^{tT}(y_{R}^{tT})^{*}]}{16\pi^{2}}C_{1},$
(20)

with $C_{1}$ defined as

$\displaystyle
C_{1}=\frac{1}{4m_{t}^{2}}[B_{0}(m_{t}^{2},m_{T}^{2},m_{h}^{2})-B_{0}(0,m_{T}^{2},m_{T}^{2})+(m_{T}^{2}-m_{t}^{2}-m_{h}^{2})C_{0}(m_{t}^{2},0,m_{t}^{2},m_{h}^{2},m_{T}^{2},m_{T}^{2})].$

As can be seen from the identity
$[y_{R}^{tT}(y_{L}^{tT})^{*}-y_{L}^{tT}(y_{R}^{tT})^{*}]=2i(\mathrm{Re}y_{L}^{tT}\mathrm{Im}y_{R}^{tT}-\mathrm{Re}y_{R}^{tT}\mathrm{Im}y_{L}^{tT})$,
$d_{t}^{EDM}$ will vanish if the imaginary parts of $y_{L,R}^{tT}$ are both
turned off. If we take $m_{T}=400~{}\mathrm{GeV}$, top EDM sets the upper
limit of $|y_{R}^{tT}(y_{L}^{tT})^{*}-y_{L}^{tT}(y_{R}^{tT})^{*}|$ to be 0.12
at $90\%$ CL. If we take $m_{T}=800~{}\mathrm{GeV}$, the corresponding upper
limit of $|y_{R}^{tT}(y_{L}^{tT})^{*}-y_{L}^{tT}(y_{R}^{tT})^{*}|$ is 0.24 at
$90\%$ CL. The larger $m_{T}$ is, the looser will the constraints be.

##  IV The analysis of double Higgs production

###  IV.1 New physics results of the amplitude

The double Higgs production is a hot topic in the field of Higgs physics. The
di-Higgs production cross section has been calculated in SM for many years
[49, 50]. The new physics effects have also drawn much attention of this
community [51, 52, 53, 54, 55]. Some works on di-Higgs production are based on
the SM effective field theory (EFT) [56, 57, 58, 59, 60, 61] and non-linearly
realized EFT [62, 63, 64, 65]. There are also many studies considered in
specific models, for example, Higgs singlet model [66, 67, 68], two Higgs
doublet model [69, 70, 71, 72], VLQ models [73, 74, 75], composite Higgs
models [76, 77], minimal supersymmetric standard model (MSSM) [78, 79, 80],
next-to-MSSM [81, 82], and many other new physics models.

For VLQ models, there are additional fermion contributions: the pure new quark
loops and the loops with both SM and new quarks. In Fig. 3 and Fig. 4, we show
the Feynman diagrams from the pure quark loops and mixed quark loops,
respectively 111The diagrams are drawn by JaxoDraw [83].. The latter will
induced by the FCNY interactions. The FCNY contributions are less considered
in most of the studies, because they are small compared to the same flavor
terms. As a second thought, this channel can be sensitive to large FCNY
couplings.

Figure 3: Typical Feynman diagrams contributing to the
$g(k_{1},\mu)g(k_{2},\nu)\rightarrow h(p_{1})h(p_{2})$ production with pure
top (also $T$) quarks running in the loops, where the counter-clockwise
diagrams should be included. Figure 4: Typical Feynman diagrams contributing
to the $g(k_{1},\mu)g(k_{2},\nu)\rightarrow h(p_{1})h(p_{2})$ production with
both top and $T$ quarks running in the loops, where the counter-clockwise
diagrams should be included.

Starting from the Lagrangian in Eq. (II.2), the amplitude of $gg\rightarrow
hh$ can be parametrized as

$\displaystyle
i\mathcal{M}(\hat{s})=-i\frac{g_{s}^{2}\hat{s}}{16\pi^{2}v^{2}}\epsilon_{\mu}^{a,r_{1}}(k_{1})\epsilon_{\nu}^{a,r_{2}}(k_{2})(A^{\mu\nu}f_{A}+B^{\mu\nu}f_{B}+C^{\mu\nu}f_{C}),$
(21)

where $a$ and $r_{1,2}$ are the color and spin indices, and the tensor
structures $A^{\mu\nu},~{}B^{\mu\nu},~{}C^{\mu\nu}$ are given by

$\displaystyle A^{\mu\nu}\equiv
g^{\mu\nu}-\frac{k_{2}^{\mu}k_{1}^{\nu}}{k_{1}\cdot k_{2}},\qquad\qquad\qquad
C^{\mu\nu}\equiv\frac{k_{1\rho}k_{2\sigma}\epsilon^{\mu\nu\rho\sigma}}{k_{1}\cdot
k_{2}},$ $\displaystyle B^{\mu\nu}\equiv
g^{\mu\nu}+\frac{m_{h}^{2}k_{2}^{\mu}k_{1}^{\nu}}{p_{T}^{2}(k_{1}\cdot
k_{2})}-\frac{2(k_{1}\cdot p_{1})k_{2}^{\mu}p_{1}^{\nu}}{p_{T}^{2}(k_{1}\cdot
k_{2})}-\frac{2(k_{2}\cdot p_{1})p_{1}^{\mu}k_{1}^{\nu}}{p_{T}^{2}(k_{1}\cdot
k_{2})}+\frac{2p_{1}^{\mu}p_{1}^{\nu}}{p_{T}^{2}}\quad(p_{T}^{2}\equiv\frac{\hat{t}\hat{u}-m_{h}^{4}}{\hat{s}}).$
(22)

They have the orthonormal relations
$A^{\mu\nu}A_{\mu\nu}=B^{\mu\nu}B_{\mu\nu}=C^{\mu\nu}C_{\mu\nu}=2$ and
$A^{\mu\nu}B_{\mu\nu}=A^{\mu\nu}C_{\mu\nu}=B^{\mu\nu}C_{\mu\nu}=0$. The
coefficients $f_{A,B,C}$ receive contributions from the $t$ and $T$ quark
loops. They can be written as $f_{A,B,C}\equiv
f_{A,B,C}^{t}+f_{A,B,C}^{T}+f_{A,B,C}^{tT}$. Here $f_{A,B,C}^{t(T)}$ mark the
contributions from pure $t(T)$ quark loops, while $f_{A,B,C}^{tT}$ labels the
contributions from mixed $t$ and $T$ quark loops. When we set
$\kappa_{t}^{2}=1$ and turn off the couplings $y_{L,R}^{tT}$, they will go to
the SM result. After some lengthy calculations, we can obtain their explicit
expressions.

The pure top quark contribution to $f_{A}$ is given by222During the
calculations, we have used the FeynCalc to simplify the results [84, 85].
$f_{A}^{t}=\kappa_{t}f_{A}^{t,\triangle}+\kappa_{t}^{2}f_{A}^{t,\Box 1}$.
Here, $f_{A}^{t,\triangle}$ and $f_{A}^{t,\Box 1}$ are defined as

$\displaystyle
f_{A}^{t,\triangle}=\frac{6m_{h}^{2}m_{t}^{2}}{\hat{s}(\hat{s}-m_{h}^{2})}[2+(4m_{t}^{2}-\hat{s})C_{0}^{t}(\hat{s})],$
$\displaystyle f_{A}^{t,\Box
1}=\frac{2m_{t}^{2}}{\hat{s}}\Big{\\{}4m_{t}^{2}C_{0}^{t}(\hat{s})+\frac{2(m_{h}^{2}-4m_{t}^{2})}{\hat{s}}[(\hat{t}-m_{h}^{2})C_{0}^{t}(\hat{t})+(\hat{u}-m_{h}^{2})C_{0}^{t}(\hat{u})]$
$\displaystyle+m_{t}^{2}(8m_{t}^{2}-\hat{s}-2m_{h}^{2})[D_{0}^{t}(\hat{t},\hat{s})+D_{0}^{t}(\hat{u},\hat{s})+D_{0}^{t}(\hat{t},\hat{u})]+2+\frac{\hat{t}\hat{u}-m_{h}^{4}}{\hat{s}}(4m_{t}^{2}-m_{h}^{2})D_{0}^{t}(\hat{t},\hat{u})\Big{\\}}.$
(23)

The pure $T$ quark contribution to $f_{A}$ is given by
$f_{A}^{T}=(-\frac{vy_{T}}{m_{T}})f_{A}^{T,\triangle}+(\frac{vy_{T}}{m_{T}})^{2}f_{A}^{T,\Box
1}$. Here, $f_{A}^{T,\triangle}$ and $f_{A}^{T,\Box 1}$ are defined as

$\displaystyle
f_{A}^{T,\triangle}=\frac{6m_{h}^{2}m_{T}^{2}}{\hat{s}(\hat{s}-m_{h}^{2})}[2+(4m_{T}^{2}-\hat{s})C_{0}^{T}(\hat{s})],$
$\displaystyle f_{A}^{T,\Box
1}=\frac{2m_{T}^{2}}{\hat{s}}\Big{\\{}4m_{T}^{2}C_{0}^{T}(\hat{s})+\frac{2(m_{h}^{2}-4m_{T}^{2})}{\hat{s}}[(\hat{t}-m_{h}^{2})C_{0}^{T}(\hat{t})+(\hat{u}-m_{h}^{2})C_{0}^{T}(\hat{u})]$
$\displaystyle+m_{T}^{2}(8m_{T}^{2}-\hat{s}-2m_{h}^{2})[D_{0}^{T}(\hat{t},\hat{s})+D_{0}^{T}(\hat{u},\hat{s})+D_{0}^{T}(\hat{t},\hat{u})]+2+\frac{\hat{t}\hat{u}-m_{h}^{4}}{\hat{s}}(4m_{T}^{2}-m_{h}^{2})D_{0}^{T}(\hat{t},\hat{u})\Big{\\}}.$
(24)

The top and $T$ quark mixed contribution to $f_{A}$ is given by
$f_{A}^{tT}=(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})f_{A}^{tT,\Box
1}+[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]f_{A}^{tT,\Box 2}$.
Here, $f_{A}^{tT,\Box 1}$ and $f_{A}^{tT,\Box 2}$ are defined as

$\displaystyle f_{A}^{tT,\Box 1}$ $\displaystyle=$
$\displaystyle\frac{2v^{2}}{\hat{s}}\Big{\\{}2m_{t}^{2}C_{0}^{t}(\hat{s})+2m_{T}^{2}C_{0}^{T}(\hat{s})+\frac{m_{h}^{2}-m_{t}^{2}-m_{T}^{2}}{\hat{s}}[(\hat{t}-m_{h}^{2})(C_{0}^{tT}(\hat{t})+C_{0}^{Tt}(\hat{t}))+(\hat{u}-m_{h}^{2})(C_{0}^{tT}(\hat{u})+C_{0}^{Tt}(\hat{u}))]$
$\displaystyle+(m_{t}^{2}+m_{T}^{2}-m_{h}^{2})[m_{t}^{2}(D_{0}^{tT}(\hat{t},\hat{s})+D_{0}^{tT}(\hat{u},\hat{s})+D_{0}^{tT}(\hat{t},\hat{u}))+m_{T}^{2}(D_{0}^{Tt}(\hat{t},\hat{s})+D_{0}^{Tt}(\hat{u},\hat{s})+D_{0}^{Tt}(\hat{t},\hat{u}))]$
$\displaystyle+2+\frac{\hat{t}\hat{u}-m_{h}^{4}}{\hat{s}}(m_{t}^{2}+m_{T}^{2}-m_{h}^{2})D_{0}^{tT}(\hat{t},\hat{u})\Big{\\}},$
$\displaystyle f_{A}^{tT,\Box
2}=\frac{m_{t}m_{T}v^{2}}{\hat{s}^{2}}\Big{\\{}4(m_{h}^{2}-\hat{t})[C_{0}^{tT}(\hat{t})+C_{0}^{Tt}(\hat{t}))]+4(m_{h}^{2}-\hat{u})[C_{0}^{tT}(\hat{u})+C_{0}^{Tt}(\hat{u})]$
$\displaystyle+\hat{s}(4m_{t}^{2}-\hat{s})[D_{0}^{tT}(\hat{t},\hat{s})+D_{0}^{tT}(\hat{u},\hat{s})+D_{0}^{tT}(\hat{t},\hat{u})]+\hat{s}(4m_{T}^{2}-\hat{s})[D_{0}^{Tt}(\hat{t},\hat{s})+D_{0}^{Tt}(\hat{u},\hat{s})+D_{0}^{Tt}(\hat{t},\hat{u})]$
$\displaystyle+4(\hat{t}\hat{u}-m_{h}^{4})D_{0}^{tT}(\hat{t},\hat{u})\Big{\\}}.$
(25)

The pure top quark contribution to $f_{B}$ is given by
$f_{B}^{t}=\kappa_{t}^{2}f_{B}^{t,\Box 1}$. Here, $f_{B}^{t,\Box 1}$ is
defined as

$\displaystyle f_{B}^{t,\Box
1}=\frac{m_{t}^{2}}{\hat{s}}\Big{\\{}-2\hat{s}C_{0}^{t}(\hat{s})+2(m_{h}^{2}-\hat{t})C_{0}^{t}(\hat{t})+2(m_{h}^{2}-\hat{u})C_{0}^{t}(\hat{u})-2(8m_{t}^{2}+\hat{s}-2m_{h}^{2})C_{0}^{t}(m_{h}^{2})$
$\displaystyle+2m_{t}^{2}(8m_{t}^{2}+\hat{s}-2m_{h}^{2})[D_{0}^{t}(\hat{t},\hat{s})+D_{0}^{t}(\hat{u},\hat{s})+D_{0}^{t}(\hat{t},\hat{u})]$
$\displaystyle+\frac{1}{\hat{t}\hat{u}-m_{h}^{4}}\Big{[}\hat{s}\hat{t}(8m_{t}^{2}\hat{t}-\hat{t}^{2}-m_{h}^{4})D_{0}^{t}(\hat{t},\hat{s})+\hat{s}\hat{u}(8m_{t}^{2}\hat{u}-\hat{u}^{2}-m_{h}^{4})D_{0}^{t}(\hat{u},\hat{s})$
$\displaystyle+(8m_{t}^{2}+\hat{s}-2m_{h}^{2})\Big{(}\hat{s}(\hat{s}-2m_{h}^{2})C_{0}^{t}(\hat{s})+\hat{s}(\hat{s}-4m_{h}^{2})C_{0}^{t}(m_{h}^{2})+2\hat{t}(m_{h}^{2}-\hat{t})C_{0}^{t}(\hat{t})+2\hat{u}(m_{h}^{2}-\hat{u})C_{0}^{t}(\hat{u})\Big{)}\Big{]}\Big{\\}}.$
(26)

The pure $T$ quark contribution to $f_{B}$ is given by
$f_{B}^{T}=(\frac{vy_{T}}{m_{T}})^{2}f_{B}^{T,\Box 1}$. Here, $f_{B}^{T,\Box
1}$ is defined as

$\displaystyle f_{B}^{T,\Box
1}=\frac{m_{T}^{2}}{\hat{s}}\Big{\\{}-2\hat{s}C_{0}^{T}(\hat{s})+2(m_{h}^{2}-\hat{t})C_{0}^{T}(\hat{t})+2(m_{h}^{2}-\hat{u})C_{0}^{T}(\hat{u})-2(8m_{T}^{2}+\hat{s}-2m_{h}^{2})C_{0}^{T}(m_{h}^{2})$
$\displaystyle+2m_{T}^{2}(8m_{T}^{2}+\hat{s}-2m_{h}^{2})[D_{0}^{T}(\hat{t},\hat{s})+D_{0}^{T}(\hat{u},\hat{s})+D_{0}^{T}(\hat{t},\hat{u})]$
$\displaystyle+\frac{1}{\hat{t}\hat{u}-m_{h}^{4}}\Big{[}\hat{s}\hat{t}(8m_{T}^{2}\hat{t}-\hat{t}^{2}-m_{h}^{4})D_{0}^{T}(\hat{t},\hat{s})+\hat{s}\hat{u}(8m_{T}^{2}\hat{u}-\hat{u}^{2}-m_{h}^{4})D_{0}^{T}(\hat{u},\hat{s})$
$\displaystyle+(8m_{T}^{2}+\hat{s}-2m_{h}^{2})\Big{(}\hat{s}(\hat{s}-2m_{h}^{2})C_{0}^{T}(\hat{s})+\hat{s}(\hat{s}-4m_{h}^{2})C_{0}^{T}(m_{h}^{2})+2\hat{t}(m_{h}^{2}-\hat{t})C_{0}^{T}(\hat{t})+2\hat{u}(m_{h}^{2}-\hat{u})C_{0}^{T}(\hat{u})\Big{)}\Big{]}\Big{\\}}.$
(27)

The top and $T$ quark mixed contribution to $f_{B}$ is given by
$f_{B}^{tT}=(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})f_{B}^{tT,\Box
1}+[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]f_{B}^{tT,\Box 2}$.
Here, $f_{B}^{tT,\Box 1}$ is defined as

$\displaystyle f_{B}^{tT,\Box 1}$
$\displaystyle=\frac{2v^{2}}{\hat{s}}\Big{\\{}-\frac{\hat{s}}{2}[C_{0}^{t}(\hat{s})+C_{0}^{T}(\hat{s})]-\frac{1}{2}(\hat{s}-2m_{h}^{2}+2m_{t}^{2}+2m_{T}^{2})[C_{0}^{tT}(m_{h}^{2})+C_{0}^{Tt}(m_{h}^{2})]+\frac{m_{h}^{2}-\hat{t}}{2}[C_{0}^{tT}(t)+C_{0}^{Tt}(t)]$
$\displaystyle+\frac{m_{h}^{2}-\hat{u}}{2}[C_{0}^{tT}(u)+C_{0}^{Tt}(u)]+\frac{1}{2}(m_{t}^{2}-m_{T}^{2})(\hat{s}+m_{t}^{2}+m_{T}^{2}-m_{h}^{2})[D_{0}^{tT}(\hat{t},\hat{s})+D_{0}^{tT}(\hat{u},\hat{s})-D_{0}^{Tt}(\hat{t},\hat{s})-D_{0}^{Tt}(\hat{u},\hat{s})]$
$\displaystyle+\frac{(m_{t}^{2}+m_{T}^{2})}{4}(\hat{s}-2m_{h}^{2}+2m_{t}^{2}+2m_{T}^{2})[D_{0}^{tT}(\hat{t},\hat{s})+D_{0}^{tT}(\hat{u},\hat{s})+D_{0}^{tT}(\hat{t},\hat{u})+D_{0}^{Tt}(\hat{t},\hat{s})+D_{0}^{Tt}(\hat{u},\hat{s})+D_{0}^{Tt}(\hat{t},\hat{u})]$
$\displaystyle+\frac{1}{4(\hat{t}\hat{u}-m_{h}^{4})}\Big{[}\hat{s}(\hat{s}-2m_{h}^{2})(\hat{s}-2m_{h}^{2}+2m_{t}^{2}+2m_{T}^{2})\Big{(}C_{0}^{t}(\hat{s})+C_{0}^{T}(\hat{s})\Big{)}$
$\displaystyle+2\hat{s}(m_{T}^{2}-m_{t}^{2})(\hat{s}-2m_{h}^{2}+2m_{T}^{2}+2m_{t}^{2})\Big{(}C_{0}^{t}(\hat{s})-C_{0}^{T}(\hat{s})\Big{)}+\hat{s}(\hat{s}-4m_{h}^{2})(\hat{s}-2m_{h}^{2}+2m_{t}^{2}+2m_{T}^{2})\Big{(}C_{0}^{tT}(m_{h}^{2})+C_{0}^{Tt}(m_{h}^{2})\Big{)}$
$\displaystyle+2\hat{t}(m_{h}^{2}-\hat{t})(\hat{s}-2m_{h}^{2}+2m_{t}^{2}+2m_{T}^{2})\Big{(}C_{0}^{tT}(\hat{t})+C_{0}^{Tt}(\hat{t})\Big{)}+2\hat{u}(m_{h}^{2}-\hat{u})(\hat{s}-2m_{h}^{2}+2m_{t}^{2}+2m_{T}^{2})\Big{(}C_{0}^{tT}(\hat{u})+C_{0}^{Tt}(\hat{u})\Big{)}$
$\displaystyle+\hat{s}(m_{t}^{2}-m_{T}^{2})^{2}(\hat{s}-2m_{h}^{2}+2m_{t}^{2}+2m_{T}^{2})\Big{(}D_{0}^{tT}(\hat{t},\hat{s})+D_{0}^{tT}(\hat{u},\hat{s})+D_{0}^{tT}(\hat{t},\hat{u})+D_{0}^{Tt}(\hat{t},\hat{s})+D_{0}^{Tt}(\hat{u},\hat{s})+D_{0}^{Tt}(\hat{t},\hat{u})\Big{)}$
$\displaystyle+2\hat{s}\hat{t}(m_{t}^{2}-m_{T}^{2})(\hat{s}-2m_{h}^{2}+2m_{t}^{2}+2m_{T}^{2})\Big{(}D_{0}^{tT}(\hat{t},\hat{s})-D_{0}^{Tt}(\hat{t},\hat{s})\Big{)}$
$\displaystyle+2\hat{s}\hat{u}(m_{t}^{2}-m_{T}^{2})(\hat{s}-2m_{h}^{2}+2m_{t}^{2}+2m_{T}^{2})\Big{(}D_{0}^{tT}(\hat{u},\hat{s})-D_{0}^{Tt}(\hat{u},\hat{s})\Big{)}$
$\displaystyle-\hat{s}\hat{t}\Big{(}\hat{t}^{2}+m_{h}^{4}-2\hat{t}(m_{t}^{2}+m_{T}^{2})\Big{)}\Big{(}D_{0}^{tT}(\hat{t},\hat{s})+D_{0}^{tT}(\hat{u},\hat{s})\Big{)}-\hat{s}\hat{u}\Big{(}\hat{u}^{2}+m_{h}^{4}-2\hat{u}(m_{t}^{2}+m_{T}^{2})\Big{)}\Big{(}D_{0}^{Tt}(\hat{t},\hat{s})+D_{0}^{Tt}(\hat{u},\hat{s})\Big{)}\Big{]}\Big{\\}},$

and $f_{B}^{tT,\Box 2}$ is defined as

$\displaystyle f_{B}^{tT,\Box
2}=\frac{2m_{t}m_{T}v^{2}}{\hat{s}}\Big{\\{}-2[C_{0}^{tT}(m_{h}^{2})+C_{0}^{Tt}(m_{h}^{2})]+2m_{t}^{2}[D_{0}^{tT}(\hat{t},\hat{s})+D_{0}^{tT}(\hat{u},\hat{s})+D_{0}^{tT}(\hat{t},\hat{u})]$
$\displaystyle+2m_{T}^{2}[D_{0}^{Tt}(\hat{t},\hat{s})+D_{0}^{Tt}(\hat{u},\hat{s})+D_{0}^{Tt}(\hat{t},\hat{u})]+\frac{1}{\hat{t}\hat{u}-m_{h}^{4}}\Big{[}\hat{s}(\hat{s}-2m_{h}^{2})\Big{(}C_{0}^{t}(\hat{s})+C_{0}^{T}(\hat{s})\Big{)}-2\hat{s}(m_{t}^{2}-m_{T}^{2})\Big{(}C_{0}^{t}(\hat{s})-C_{0}^{T}(\hat{s})\Big{)}$
$\displaystyle+\hat{s}(\hat{s}-4m_{h}^{2})\Big{(}C_{0}^{tT}(m_{h}^{2})+C_{0}^{Tt}(m_{h}^{2})\Big{)}+2\hat{t}(m_{h}^{2}-\hat{t})\Big{(}C_{0}^{tT}(\hat{t})+C_{0}^{Tt}(\hat{t})\Big{)}+2\hat{u}(m_{h}^{2}-\hat{u})\Big{(}C_{0}^{tT}(\hat{u})+C_{0}^{Tt}(\hat{u})\Big{)}$
$\displaystyle+\hat{s}\Big{(}\hat{t}^{2}+(m_{t}^{2}-m_{T}^{2})^{2}\Big{)}\Big{(}D_{0}^{tT}(\hat{t},\hat{s})+D_{0}^{Tt}(\hat{t},\hat{s})\Big{)}+\hat{s}\Big{(}\hat{u}^{2}+(m_{t}^{2}-m_{T}^{2})^{2}\Big{)}\Big{(}D_{0}^{tT}(\hat{u},\hat{s})+D_{0}^{Tt}(\hat{u},\hat{s})\Big{)}$
$\displaystyle+2\hat{s}\hat{t}(m_{t}^{2}-m_{T}^{2})\Big{(}D_{0}^{tT}(\hat{t},\hat{s})-D_{0}^{Tt}(\hat{t},\hat{s})\Big{)}+2\hat{s}\hat{u}(m_{t}^{2}-m_{T}^{2})\Big{(}D_{0}^{tT}(\hat{u},\hat{s})-D_{0}^{Tt}(\hat{u},\hat{s})\Big{)}+2\hat{s}(m_{t}^{2}-m_{T}^{2})^{2}D_{0}^{tT}(\hat{t},\hat{u})\Big{]}\Big{\\}}.$
(29)

The pure top quark contribution to $f_{C}$ is proportional to
$\widetilde{\kappa}_{t}$, so we set $f_{C}^{t}$ to be zero. Similarly, the
pure $T$ quark contribution to $f_{C}$ is also turned off. The top and $T$
quark mixed contribution to $f_{C}$ is given by
$f_{C}^{tT}=-i[y_{L}^{tT}(y_{R}^{tT})^{*}-y_{R}^{tT}(y_{L}^{tT})^{*}]f_{C}^{tT,\Box}$.
Here, $f_{C}^{tT,\Box}$ is defined as

$\displaystyle
f_{C}^{tT,\Box}=m_{t}m_{T}v^{2}[D_{0}^{tT}(\hat{t},\hat{s})+D_{0}^{tT}(\hat{u},\hat{s})+D_{0}^{tT}(\hat{t},\hat{u})+D_{0}^{Tt}(\hat{t},\hat{s})+D_{0}^{Tt}(\hat{u},\hat{s})+D_{0}^{Tt}(\hat{t},\hat{u})].$
(30)

###  IV.2 Heavy quark expansion

In the limit of $\frac{m_{h}^{2},\hat{s},\hat{t},\hat{u}}{m_{t}^{2}}\ll 1$,
the coefficients of pure top quark loops $f_{A}^{t},~{}f_{B}^{t},~{}f_{C}^{t}$
can be expanded as

$\displaystyle
f_{A}^{t,\triangle}=\frac{2m_{h}^{2}}{\hat{s}-m_{h}^{2}}\left[1+\frac{7\hat{s}}{120m_{t}^{2}}+\mathcal{O}(\frac{1}{m_{t}^{4}})\right],\quad
f_{A}^{t,\Box
1}=-\frac{2}{3}\left[1+\frac{7m_{h}^{2}}{20m_{t}^{2}}+\mathcal{O}(\frac{1}{m_{t}^{4}})\right],$
$\displaystyle f_{B}^{t,\Box
1}=\frac{11(m_{h}^{4}-\hat{t}\hat{u})}{90m_{t}^{2}\hat{s}}+\mathcal{O}(\frac{1}{m_{t}^{4}}),\quad
f_{B}^{t,\Box
2}=\frac{m_{h}^{4}-\hat{t}\hat{u}}{6m_{t}^{2}\hat{s}}+\mathcal{O}(\frac{1}{m_{t}^{4}}),$
$\displaystyle
f_{C}^{t,\triangle}=\frac{3m_{h}^{2}}{\hat{s}-m_{h}^{2}}\left[1+\frac{\hat{s}}{12m_{t}^{2}}+\mathcal{O}(\frac{1}{m_{t}^{4}})\right],\quad
f_{C}^{t,\Box}=-2\left[1+\frac{m_{h}^{2}}{3m_{t}^{2}}+\mathcal{O}(\frac{1}{m_{t}^{4}})\right].$
(31)

For the coefficients of pure $T$ quark loops, they are just the ones with
$m_{t}$ replaced by $m_{T}$.

For the case of $t$ and $T$ quark mixed loops, things are more complicated. In
the limit of $\frac{m_{h}^{2},\hat{s},\hat{t},\hat{u}}{m_{t,T}^{2}}\ll 1$, the
coefficient $f_{A}^{tT},~{}f_{B}^{tT},~{}f_{C}^{tT}$ can be expanded as

$\displaystyle f_{A}^{tT,\Box 1}=\mathcal{O}(\frac{1}{m_{t,T}^{4}}),$
$\displaystyle f_{B}^{tT,\Box
1}=\frac{v^{2}(\hat{t}^{2}-\hat{u}^{2})}{m_{T}^{2}(\hat{t}\hat{u}-m_{h}^{4})}\cdot\frac{(1+r_{tT}^{2})(1+2r_{tT}^{2}\log
r_{tT}^{2}-r_{tT}^{4})}{2r_{tT}^{2}(1-r_{tT}^{2})^{2}}+\mathcal{O}(\frac{1}{m_{t,T}^{4}}),$
$\displaystyle
f_{C}^{tT,\Box}=\frac{v^{2}}{m_{t}m_{T}}+\mathcal{O}(\frac{1}{m_{t,T}^{4}}).$
(32)

###  IV.3 The cross section analysis

When averaging the initial spin and color degrees of freedom, we can get the
partonic cross section of $gg\rightarrow hh$ at leading order (LO) as follows

$\displaystyle\hat{\sigma}_{LO}(gg\rightarrow
hh;\hat{s})=\frac{\alpha_{S}^{2}G_{F}^{2}\sqrt{\hat{s}(\hat{s}-4m_{h}^{2})}}{128(4\pi)^{3}}\int_{-1}^{1}d\cos\theta~{}(|f_{A}|^{2}+|f_{B}|^{2}+|f_{C}|^{2})$
$\displaystyle=\frac{\alpha_{S}^{2}G_{F}^{2}}{64(4\pi)^{3}}\int_{\hat{t}_{min}}^{\hat{t}_{max}}d\hat{t}~{}(|f_{A}|^{2}+|f_{B}|^{2}+|f_{C}|^{2})(\hat{t}_{min}=-\frac{1}{4}(\sqrt{\hat{s}}+\sqrt{\hat{s}-4m_{h}^{2}})^{2},~{}\hat{t}_{max}=-\frac{1}{4}(\sqrt{\hat{s}}-\sqrt{\hat{s}-4m_{h}^{2}})^{2}),$
(33)

where $f_{A},f_{B},f_{C}$ are calculated as

$\displaystyle
f_{A}=f_{A}^{t}+f_{A}^{T}+f_{A}^{tT}=\kappa_{t}f_{A}^{t,\triangle}+\kappa_{t}^{2}f_{A}^{t,\Box
1}+(-\frac{vy_{T}}{m_{T}})f_{A}^{T,\triangle}+(\frac{vy_{T}}{m_{T}})^{2}f_{A}^{T,\Box
1}$ $\displaystyle+(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})f_{A}^{tT,\Box
1}+[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]f_{A}^{tT,\Box 2},$
$\displaystyle
f_{B}=f_{B}^{t}+f_{B}^{T}+f_{B}^{tT}=\kappa_{t}^{2}f_{B}^{t,\Box
1}+(\frac{vy_{T}}{m_{T}})^{2}f_{B}^{T,\Box
1}+(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})f_{B}^{tT,\Box
1}+[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]f_{B}^{tT,\Box 2},$
$\displaystyle
f_{C}=f_{C}^{t}+f_{C}^{T}+f_{C}^{tT}=-i[y_{L}^{tT}(y_{R}^{tT})^{*}-y_{R}^{tT}(y_{L}^{tT})^{*}]f_{C}^{tT,\Box}.$
(34)

Note that there is a $\frac{1}{2}$ factor in the partonic cross section
because of the identical final states. In general, the anomalous triple Higgs
coupling $\lambda_{hhh}$ will also alter the di-Higgs production cross
section. Its effects can be captured with
$f_{A}^{f,\triangle},f_{C}^{f,\triangle}(f=t,T)$ multiplied by the factor
$1+\delta_{hhh}\equiv\lambda_{hhh}/\lambda_{hhh}^{SM}$.

After folding the partonic cross section with the gluon luminosity, we can get
the hadron level cross section

$\displaystyle\sigma_{LO}(pp\rightarrow
hh)=\int_{\frac{4m_{h}^{2}}{s}}^{1}d\tau\int_{\tau}^{1}\frac{dx}{x}f(x,\mu_{F}^{2})f(\frac{\tau}{x},\mu_{F}^{2})\hat{\sigma}_{LO}(gg\rightarrow
hh;\hat{s}=\tau s),$ (35)

where $f$ represents the gluon parton distribution function (PDF) and
$\mu_{F}$ is the factorization scale.

##  V The numerical results and constraint prospects

Just similar to the VLQT model, we take
$\kappa_{t}=c_{L}^{2},y_{T}=-\frac{m_{T}}{v}s_{L}^{2}$ for simplicity, but let
$\mathrm{Re}(y_{L}^{tT}),\mathrm{Re}(y_{R}^{tT}),\mathrm{Im}(y_{L}^{tT}),\mathrm{Im}(y_{R}^{tT})$
to be free. Then we can choose several benchmark scenarios and estimate the
constraints on the magnitude and sign of the FCNY couplings. Now, we need to
normalize the cross section to the SM ones numerically for fixed $m_{T}$ and
$s_{L}$, which is defined as

$\displaystyle\mu_{hh}\equiv\frac{\sigma_{LO}(pp\rightarrow
hh)}{\sigma_{LO}^{SM}(pp\rightarrow hh)}.$ (36)

Up to LO level, $\mu_{hh}$ can be parametrized as

$\displaystyle\mu_{hh}=1+A_{1}+A_{0}^{hhh}\delta_{hhh}+A_{1}^{hhh}\delta_{hhh}^{2}+(A_{2}+A_{2}^{hhh}\delta_{hhh})(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})$
$\displaystyle+(A_{3}+A_{3}^{hhh}\delta_{hhh})[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]+A_{4}(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})^{2}+A_{5}[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]^{2}$
$\displaystyle+A_{6}(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]-A_{7}[y_{L}^{tT}(y_{R}^{tT})^{*}-y_{R}^{tT}(y_{L}^{tT})^{*}]^{2},$
(37)

From the observation of Eq. (IV.3) and Eq. (IV.3), we can find that
$A_{1},~{}A_{2},~{}A_{3}$,
$A_{0}^{hhh},~{}A_{1}^{hhh},~{}A_{2}^{hhh},~{}A_{3}^{hhh}$ depend on the
choices of both $m_{T}$ and $s_{L}$, while $A_{4},~{}A_{5},~{}A_{6},~{}A_{7}$
only depend on $m_{T}$. Moreover, $A_{1}^{hhh}$, $A_{4},~{}A_{5},~{}A_{7}$ are
always non-negative and $A_{1}$ vanishes as $s_{L}$ goes to zero.

Although $\sigma^{SM}(pp\rightarrow hh)$ has been calculated with high
precision [86, 87, 88, 89, 90, 91, 92, 93, 94], we will not do that hard work
here. We only keep the LO results because a large part of the QCD corrections
can be cancelled in the ratio [57, 95, 96, 65]. To get the numerical results
of cross sections, we write a model file through FeynRules [97, 98], FeynArts
[99] and NLOCT [100]. Then it is linked to MadGraph [101]. Before the
numerical calculations, we take the following default settings:
$\bullet$ Proton contains $b,\bar{b}$, that is, we use the 5FS (5 flavor
scheme).
$\bullet$ We adopt the PDF choice of "MSTW2008lo68cl" (LHAPDF ID 21000).
$\bullet$ The default $dynamical_{\\_}scale_{\\_}choice$ is set to be 3 (see
[102]).
$\bullet$ The input parameters are choose as
$m_{h}=125.09~{}\mathrm{GeV},~{}G_{F}=1.1664\times
10^{-5}~{}\mathrm{GeV}^{-2}$, $m_{t}=172.74~{}\mathrm{GeV}$, and
$\alpha_{s}(m_{Z})=0.1184$. Thus we have $v=246.221~{}\mathrm{GeV}$.

Currently, the Higgs pair production is bounded to be $|\mu_{hh}|\leq 6.9$ at
$95\%$ confidence level (CL) [103, 104]. At the high luminosity LHC (HL-LHC),
di-Higgs production measurement is accessible. The expected signal strength is
$\mu_{hh}=1.00_{-0.39}^{+0.41}$ with $1\sigma$ uncertainty [105]. We take the
benchmark points as $m_{T}=400~{}\mathrm{GeV},~{}s_{L}=0.2$ and
$m_{T}=800~{}\mathrm{GeV},~{}s_{L}=0.1$, and all the following discussions are
based on the two benchmark points. Now, we should determine the specific
values of $A_{1},~{}A_{2},~{}A_{3},~{}A_{4},~{}A_{5},~{}A_{6},~{}A_{7}$ and
$A_{0}^{hhh},~{}A_{1}^{hhh},~{}A_{2}^{hhh},~{}A_{3}^{hhh}$. First of all, we
have $\sigma_{LO}^{SM}(pp\rightarrow hh)=$ 24.7 fb. When setting different
values of $\delta_{hhh},~{}y_{L}^{tT},~{}y_{R}^{tT}$, we can obtain different
normalized cross sections (see Tab. 1 and Tab. 2). Then the numerical values
of $A_{1},...,A_{7}$ and
$A_{0}^{hhh},~{}A_{1}^{hhh},~{}A_{2}^{hhh},~{}A_{3}^{hhh}$ can be solved from
the first seven and last four equations individually. Their results are given
in Tab. 3.

$\delta_{hhh}$ | $(~{}y_{L}^{tT},~{}y_{R}^{tT})$ | expressions of $\mu_{hh}$ | numerical values of $\mu_{hh}$
---|---|---|---
0 | $(0,~{}0)$ | $1+A_{1}$ | 0.8081
$(0,~{}1)$ | $1+A_{1}+A_{2}+A_{4}$ | 1.254
$(0,~{}\frac{1}{2})$ | $1+A_{1}+\frac{1}{4}A_{2}+\frac{1}{16}A_{4}$ | 0.9057
$(1,~{}1)$ | $1+A_{1}+2A_{2}+2A_{3}+4A_{4}+4A_{5}+4A_{6}$ | 10.92
$(1,~{}-1)$ | $1+A_{1}+2A_{2}-2A_{3}+4A_{4}+4A_{5}-4A_{6}$ | 1.695
$(1,~{}i)$ | $1+A_{1}+2A_{2}+4A_{4}+4A_{7}$ | 14.13
$(\frac{1}{2},~{}\frac{1}{2})$ | $1+A_{1}+\frac{1}{2}A_{2}+\frac{1}{2}A_{3}+\frac{1}{4}A_{4}+\frac{1}{4}A_{5}+\frac{1}{4}A_{6}$ | 2.206
1 | $(0,~{}0)$ |  $1+A_{1}+A_{0}^{hhh}+A_{1}^{hhh}$ | 0.3877
$(0,~{}1)$ | $1+A_{1}+A_{0}^{hhh}+A_{1}^{hhh}+A_{2}+A_{2}^{hhh}+A_{4}$ | 0.6996
$(1,~{}1)$ |  $1+A_{1}+A_{0}^{hhh}+A_{1}^{hhh}+2(A_{2}+A_{2}^{hhh})+2(A_{3}+A_{3}^{hhh})$ $+4A_{4}+4A_{5}+4A_{6}$ | 8.235
$-1$ | $(0,~{}0)$ | $1+A_{1}-A_{0}^{hhh}+A_{1}^{hhh}$ | 1.779
Table 1: The normalized cross sections for different $\delta_{hhh},~{}y_{L}^{tT},~{}y_{R}^{tT}$ values with $m_{T}=$ 400 GeV and $s_{L}=0.2$ at $\sqrt{s}$=14TeV. $\delta_{hhh}$ | $(~{}y_{L}^{tT},~{}y_{R}^{tT})$ | expressions of $\mu_{hh}$ | numerical values of $\mu_{hh}$
---|---|---|---
0 | $(0,~{}0)$ | $1+A_{1}$ | 0.9506
$(0,~{}1)$ | $1+A_{1}+A_{2}+A_{4}$ | 1.098
$(0,~{}\frac{1}{2})$ | $1+A_{1}+\frac{1}{4}A_{2}+\frac{1}{16}A_{4}$ | 0.9838
$(1,~{}1)$ | $1+A_{1}+2A_{2}+2A_{3}+4A_{4}+4A_{5}+4A_{6}$ | 5.255
$(1,~{}-1)$ | $1+A_{1}+2A_{2}-2A_{3}+4A_{4}+4A_{5}-4A_{6}$ | 0.3376
$(1,~{}i)$ | $1+A_{1}+2A_{2}+4A_{4}+4A_{7}$ | 5.247
$(\frac{1}{2},~{}\frac{1}{2})$ | $1+A_{1}+\frac{1}{2}A_{2}+\frac{1}{2}A_{3}+\frac{1}{4}A_{4}+\frac{1}{4}A_{5}+\frac{1}{4}A_{6}$ | 1.675
1 | $(0,~{}0)$ |  $1+A_{1}+A_{0}^{hhh}+A_{1}^{hhh}$ | 0.4502
$(0,~{}1)$ | $1+A_{1}+A_{0}^{hhh}+A_{1}^{hhh}+A_{2}+A_{2}^{hhh}+A_{4}$ | 0.5526
$(1,~{}1)$ |  $1+A_{1}+A_{0}^{hhh}+A_{1}^{hhh}+2(A_{2}+A_{2}^{hhh})+2(A_{3}+A_{3}^{hhh})$ $+4A_{4}+4A_{5}+4A_{6}$ | 3.519
$-1$ | $(0,~{}0)$ | $1+A_{1}-A_{0}^{hhh}+A_{1}^{hhh}$ | 2.013
Table 2: The normalized cross sections for different $\delta_{hhh},~{}y_{L}^{tT},~{}y_{R}^{tT}$ values with $m_{T}=$ 800 GeV and $s_{L}=0.1$ at $\sqrt{s}$=14TeV. $\sqrt{s}$ (TeV) | ($m_{T}$/GeV, $s_{L}$) | $A_{1}$ | $A_{2}$ | $A_{3}$ | $A_{4}$ | $A_{5}$ | $A_{6}$ | $A_{7}$
---|---|---|---|---|---|---|---|---
14 | (400, 0.2) | $-0.1919$ | 0.3717 | 1.672 | 0.07449 | 1.114 | 0.3166 | 3.071
(800, 0.1) | $-0.04939$ | 0.1279 | 1.087 | 0.01943 | 0.378 | 0.0711 | 0.9907
$\sqrt{s}$ (TeV) | ($m_{T}$/GeV, $s_{L}$) | $A_{0}^{hhh}$ | $A_{1}^{hhh}$ | $A_{2}^{hhh}$ | $A_{3}^{hhh}$ |  |  |
14 | (400, 0.2) | $-0.6958$ | 0.2754 | $-0.1343$ | $-0.9956$ |  |  |
(800, 0.1) | $-0.7814$ | 0.281 | $-0.04494$ | $-0.5731$ |  |  |
Table 3: The coefficients in Eq. (V) solved through the signal strength values
in Tab. 1 and Tab. 2.

###  V.1 The benchmark point $m_{T}=$ 400 GeV and $s_{L}=0.2$

For the case of $m_{T}=$ 400 GeV and $s_{L}=0.2$, the numerical results of
$\mu_{hh}$ are evaluated as

$\displaystyle\mu_{hh}=1-0.1919-0.6958~{}\delta_{hhh}+0.2754~{}\delta_{hhh}^{2}+(0.3717-0.1343~{}\delta_{hhh})(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})$
$\displaystyle+(1.672-0.9956~{}\delta_{hhh})[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]+0.07449(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})^{2}+1.114[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]^{2}$
$\displaystyle+0.3166(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]-3.071[y_{L}^{tT}(y_{R}^{tT})^{*}-y_{R}^{tT}(y_{L}^{tT})^{*}]^{2}.$
(38)

In this case, the present di-Higgs production experiments give the constraints
$\delta_{hhh}\in(-3.61,6.13)$ and
$\mathrm{Re}y_{L}^{tT},\mathrm{Im}y_{L}^{tT},\mathrm{Re}y_{R}^{tT},\mathrm{Im}y_{R}^{tT}\in(-2.62,2.62)$
at $95\%$ CL by setting one parameter at a time. In Tab. 4, we give the
expected constraints on the parameters
$\delta_{hhh},\mathrm{Re}y_{L}^{tT},\mathrm{Im}y_{L}^{tT},\mathrm{Re}y_{R}^{tT},\mathrm{Im}y_{R}^{tT}$
at HL-LHC. It can be seen that both of the current and expected constraints at
HL-LHC are stronger than the unitarity bound, the reason is that the highest
power in di-Higgs production cross section is proportional to
$(y_{L,R}^{tT})^{4}$.

method parameters | $\delta_{hhh}$ | $\mathrm{Re}y_{L}^{tT}$ | $\mathrm{Im}y_{L}^{tT}$ | $\mathrm{Re}y_{R}^{tT}$ | $\mathrm{Im}y_{R}^{tT}$
---|---|---|---|---|---
individual | $1~{}\sigma$ | (-0.681, 0.327)$\cup$(2.20, 3.21) | (-1.13, 1.13) | (-1.13, 1.13) | (-1.13, 1.13) | (-1.13, 1.13)
$2~{}\sigma$ | (-1.03, 3.56) | (-1.40, 1.40) | (-1.40, 1.40) | (-1.40, 1.40) | (-1.40, 1.40)
marginalized | $1~{}\sigma$ | (-3.76, 3.99) | (-1.91, 1.91) | (-1.91, 1.91) | (-1.91, 1.91) | (-1.91, 1.91)
$2~{}\sigma$ | (-4.48, 4.66) | (-2.09, 2.09) | (-2.09, 2.09) | (-2.09, 2.09) | (-2.09, 2.09)
Table 4: The expected $1\sigma$ and $2\sigma$ bounds at HL-LHC for the
parameters
$\delta_{hhh},\mathrm{Re}y_{L}^{tT},\mathrm{Im}y_{L}^{tT},\mathrm{Re}y_{R}^{tT},\mathrm{Im}y_{R}^{tT}$
under the benchmark point $m_{T}=$ 400 GeV and $s_{L}=0.2$. Here we adopt two
different methods: (1) turn on one parameter at a time, namely the individual
method; (2) turn on all the five parameters, namely the marginalized method.

As mentioned above, there are four interesting parameters
$\mathrm{Re}(y_{L}^{tT}),\mathrm{Re}(y_{R}^{tT}),\mathrm{Im}(y_{L}^{tT}),\mathrm{Im}(y_{R}^{tT})$.
Then we can plot the reached two-dimensional parameter space by setting two of
them to be zero or imposing two conditions. Here we choose six scenarios: ①
$y_{L,R}^{tT}$ are both real (it is similar to the both imaginary number
case); ② $y_{R}^{tT}$ is real and $y_{L}^{tT}$ is imaginary (it is similar to
the real $y_{L}^{tT}$ and imaginary $y_{R}^{tT}$ case); ③ $y_{R}^{tT}=0$
(similar to the $y_{L}^{tT}=0$ case); ④ $y_{L}^{tT}=y_{R}^{tT}$; ⑤
$y_{L}^{tT}=-y_{R}^{tT}$; ⑥ $y_{L}^{tT}=(y_{R}^{tT})^{*}$.

In Fig. 5 and Fig. 6, we show the plots with $\delta_{hhh}=0$ and
$\delta_{hhh}=0.5$, respectively. From these plots, we find that
$\mathrm{Re}(y_{L,R}^{tT})$ and $\mathrm{Im}(y_{L,R}^{tT})$ are constrained to
be in the range $(-2,~{}2)$ roughly at $2~{}\sigma$ CL. In some of these
scenarios, the $2~{}\sigma$ interval can be tight as $(-0.4,~{}0.4)$. The
reach regions of $1~{}\sigma$ and $2~{}\sigma$ are quite different. The value
of $\delta_{hhh}$ has significant effects on the extraction of $y_{L,R}^{tT}$.
For the first (upper left) plot, the first and third quadrants are more
constrained. This can be understood from the Eq. (IV.3), because there is
constructive interference between the positive
$[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]$ term and other box
diagram induced terms. While it is destructive interference if
$[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]$ is negative. The
last five plots are symmetric with respect to the horizontal and vertical
axes. For the second (upper central) and sixth (lower right) plots, they
receive the contributions from the
$[y_{L}^{tT}(y_{R}^{tT})^{*}-y_{R}^{tT}(y_{L}^{tT})^{*}]$ term. Because
$|f_{C}|^{2}$ in Eq. (IV.3) is always positive, the constraints are stronger.
For the fourth (lower left) plot, positive
$[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]$ induces the
constructive interference, thus the bounds are also stronger. For the third
(upper right) plot, both
$[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]$ and
$[y_{L}^{tT}(y_{R}^{tT})^{*}-y_{R}^{tT}(y_{L}^{tT})^{*}]$ vanish, thus it is
less constrained compared to other plots. Besides, there can be more
cancellation between the triangle and box diagrams for larger $\delta_{hhh}$.
Thus the constraints are usually looser than the zero $\delta_{hhh}$ ones.

In fact, we find that the di-Higgs production at HL-LHC can give stronger
constraints than those from perturbative unitarity and $h\rightarrow\gamma Z$
decay. When we take into account the top quark EDM bound, the two scenarios
$\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ and
$y_{L}^{tT}=(y_{R}^{tT})^{*}$ can also be constrained. For the other scenarios
$y_{R}^{tT}=0,y_{L}^{tT}=\pm
y_{R}^{tT},\mathrm{Im}y_{L}^{tT}=\mathrm{Im}y_{R}^{tT}=0$, they are
insensitive to the top quark EDM. Because we have the relation
$d_{t}^{EDM}\sim
y_{R}^{tT}(y_{L}^{tT})^{*}-y_{L}^{tT}(y_{R}^{tT})^{*}=2i(\mathrm{Re}y_{L}^{tT}\mathrm{Im}y_{R}^{tT}-\mathrm{Re}y_{R}^{tT}\mathrm{Im}y_{L}^{tT})$.
For the scenarios $\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ and
$y_{L}^{tT}=(y_{R}^{tT})^{*}$, we compare the bounds from di-Higgs production
and top quark EDM for $\delta_{hhh}=0$ (Fig. 5) and $\delta_{hhh}=0.5$ (Fig.
6), respectively. From these plots, we can find that the off-axis regions can
be strongly bounded by the top EDM, while it will lose the constraining power
in the near axis regions.

Figure 5: The reach regions of $y_{L}^{tT},y_{R}^{tT}$ at HL-LHC with
$\delta_{hhh}=0$ for the case of $m_{T}$ = 400 GeV and $s_{L}=0.2$. In the
above plots, we take $\mathrm{Im}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$
(upper left), $\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ (upper
central), $y_{R}^{tT}=0$ (upper right), $y_{L}^{tT}=y_{R}^{tT}$ (lower left),
$y_{L}^{tT}=-y_{R}^{tT}$ (lower central), and $y_{L}^{tT}=(y_{R}^{tT})^{*}$
(lower right) respectively. We also take into account the top quark EDM
constraint for the scenarios
$\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ (upper central) and
$y_{L}^{tT}=(y_{R}^{tT})^{*}$ (lower right), where the reach regions of
$y_{L}^{tT},y_{R}^{tT}$ are shown in red at $90\%$ CL.

Figure 6: The reach regions of $y_{L}^{tT},y_{R}^{tT}$ at HL-LHC with
$\delta_{hhh}=0.5$ for the case of $m_{T}$ = 400 GeV and $s_{L}=0.2$. In the
above plots, we take $\mathrm{Im}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$
(upper left), $\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ (upper
central), $y_{R}^{tT}=0$ (upper right), $y_{L}^{tT}=y_{R}^{tT}$ (lower left),
$y_{L}^{tT}=-y_{R}^{tT}$ (lower central), and $y_{L}^{tT}=(y_{R}^{tT})^{*}$
(lower right) respectively. We also take into account the top quark EDM
constraint for the scenarios
$\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ (upper central) and
$y_{L}^{tT}=(y_{R}^{tT})^{*}$ (lower right), where the reach regions of
$y_{L}^{tT},y_{R}^{tT}$ are shown in red at $90\%$ CL.

###  V.2 The benchmark point $m_{T}=$ 800 GeV and $s_{L}=0.1$

For the case of $m_{T}=$ 800 GeV and $s_{L}=0.1$, the numerical results of
$\mu_{hh}$ are evaluated as

$\displaystyle\mu_{hh}=1-0.04939-0.7814~{}\delta_{hhh}+0.281~{}\delta_{hhh}^{2}+(0.1279-0.04494~{}\delta_{hhh})(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})$
$\displaystyle+(1.087-0.5731~{}\delta_{hhh})[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]+0.01943(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})^{2}+0.378[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]^{2}$
$\displaystyle+0.0711(|y_{L}^{tT}|^{2}+|y_{R}^{tT}|^{2})[y_{L}^{tT}(y_{R}^{tT})^{*}+y_{R}^{tT}(y_{L}^{tT})^{*}]-0.9907[y_{L}^{tT}(y_{R}^{tT})^{*}-y_{R}^{tT}(y_{L}^{tT})^{*}]^{2}.$
(39)

In this case, the present di-Higgs production experiments give the constraints
$\delta_{hhh}\in(-3.42,6.20)$ and
$\mathrm{Re}y_{L}^{tT},\mathrm{Im}y_{L}^{tT},\mathrm{Re}y_{R}^{tT},\mathrm{Im}y_{R}^{tT}\in(-3.81,3.81)$
at $95\%$ CL by setting one parameter at a time. In Tab. 5, we give the
expected constraints on the parameters
$\delta_{hhh},\mathrm{Re}y_{L}^{tT},\mathrm{Im}y_{L}^{tT},\mathrm{Re}y_{R}^{tT},\mathrm{Im}y_{R}^{tT}$
at HL-LHC.

method parameters | $\delta_{hhh}$ | $\mathrm{Re}y_{L}^{tT}$ | $\mathrm{Im}y_{L}^{tT}$ | $\mathrm{Re}y_{R}^{tT}$ | $\mathrm{Im}y_{R}^{tT}$
---|---|---|---|---|---
individual | $1~{}\sigma$ | (-0.499, 0.541)$\cup$(2.24, 3.28) | (-1.61, 1.61) | (-1.61, 1.61) | (-1.61, 1.61) | (-1.61, 1.61)
$2~{}\sigma$ | (-0.852, 3.63) | (-2.04, 2.04) | (-2.04, 2.04) | (-2.04, 2.04) | (-2.04, 2.04)
marginalized | $1~{}\sigma$ | (-4.12, 3.89) | (-2.82, 2.82) | (-2.82, 2.82) | (-2.82, 2.82) | (-2.82, 2.82)
$2~{}\sigma$ | (-4.80, 4.51) | (-3.05, 3.05) | (-3.05, 3.05) | (-3.05, 3.05) | (-3.05, 3.05)
Table 5: The expected $1\sigma$ and $2\sigma$ bounds at HL-LHC for the
parameters
$\delta_{hhh},\mathrm{Re}y_{L}^{tT},\mathrm{Im}y_{L}^{tT},\mathrm{Re}y_{R}^{tT},\mathrm{Im}y_{R}^{tT}$
under the benchmark point $m_{T}=$ 800 GeV and $s_{L}=0.1$. Here we adopt two
different methods: (1) turn on one parameter at a time, namely the individual
method; (2) turn on all the five parameters, namely the marginalized method.

We will also plot the reached two-dimensional parameter space by setting two
of them to be zero or imposing two conditions. In Fig. 7 and Fig. 8, similar
plots are presented for the six scenarios with $\delta_{hhh}=0$ and
$\delta_{hhh}=0.5$, respectively. From these plots, we find that
$\mathrm{Re}(y_{L,R}^{tT})$ and $\mathrm{Im}(y_{L,R}^{tT})$ are constrained to
be in the range $(-3,~{}3)$ roughly at $2~{}\sigma$ CL. In some of these
scenarios, the $2~{}\sigma$ interval can be tight as $(-0.5,~{}0.5)$. The
reach regions are similar to those in the $m_{T}=$ 400 GeV and $s_{L}=0.2$
case. For the scenarios $\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$
and $y_{L}^{tT}=(y_{R}^{tT})^{*}$, we also compare the bounds from di-Higgs
production and top quark EDM for $\delta_{hhh}=0$ (Fig. 7) and
$\delta_{hhh}=0.5$ (Fig. 8), respectively. When $m_{T}$ becomes larger,
$y_{L,R}^{tT}$ are constrained more loosely. When $s_{L}$ becomes very small,
the pure top quark contributions are SM-like and the pure $T$ quark
contributions are highly suppressed. Thus, the main deviation of $\mu_{hh}$ is
from the FCNY interactions.

Figure 7: The reach regions of $y_{L}^{tT},y_{R}^{tT}$ at HL-LHC with
$\delta_{hhh}=0$ for the case of $m_{T}$ = 800 GeV and $s_{L}=0.1$. In the
above plots, we take $\mathrm{Im}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$
(upper left), $\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ (upper
central), $y_{R}^{tT}=0$ (upper right), $y_{L}^{tT}=y_{R}^{tT}$ (lower left),
$y_{L}^{tT}=-y_{R}^{tT}$ (lower central), and $y_{L}^{tT}=(y_{R}^{tT})^{*}$
(lower right) respectively. We also take into account the top quark EDM
constraint for the scenarios
$\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ (upper central) and
$y_{L}^{tT}=(y_{R}^{tT})^{*}$ (lower right), where the reach regions of
$y_{L}^{tT},y_{R}^{tT}$ are shown in red at $90\%$ CL.

Figure 8: The reach regions of $y_{L}^{tT},y_{R}^{tT}$ at HL-LHC with
$\delta_{hhh}=0.5$ for the case of $m_{T}$ = 800 GeV and $s_{L}=0.1$. In the
above plots, we take $\mathrm{Im}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$
(upper left), $\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ (upper
central), $y_{R}^{tT}=0$ (upper right), $y_{L}^{tT}=y_{R}^{tT}$ (lower left),
$y_{L}^{tT}=-y_{R}^{tT}$ (lower central), and $y_{L}^{tT}=(y_{R}^{tT})^{*}$
(lower right) respectively. We also take into account the top quark EDM
constraint for the scenarios
$\mathrm{Re}(y_{L}^{tT})=\mathrm{Im}(y_{R}^{tT})=0$ (upper central) and
$y_{L}^{tT}=(y_{R}^{tT})^{*}$ (lower right), where the reach regions of
$y_{L}^{tT},y_{R}^{tT}$ are shown in red at $90\%$ CL.

By the way, the FCNY coupling may be probed through other processes too. For
example, we can probe the FCNY coupling $Tth$ through direct production
processes $pp\rightarrow T\bar{t}h,T\bar{t},ThW,Thj$. But they suffer from low
event rate, the detailed analyses in these channels are beyond the scope of
this work.

###  V.3 Comments on the doublet and triplet vector-like quarks

We have assumed $T_{L,R}$ to be singlets throughout this work, while they can
be components of the doublet or triplet VLQs. There are two doublet and two
triplet VLQs containing the $T$ quark:
$(X,T)_{L,R},(T,B)_{L,R},(X,T,B)_{L,R},(T,B,Y)_{L,R}$. Here $X,B,Y$ carry
$\frac{4}{3},-\frac{1}{3},-\frac{4}{3}$ electric charges, respectively. For
the doublet $(X,T)_{L,R}$, the Higgs particle only interact with the
$T_{L,R}$. For the doublet $(T,B)_{L,R}$ and triplets
$(X,T,B)_{L,R},(T,B,Y)_{L,R}$, the $B$ quark can mix with the SM bottom quark.
Thus, the Higgs particle will interact with both the $T$ and $B$ quarks. Let
us denote left (right) up-type and down-type quark mixing angles as
$\theta_{L}^{t}$ ($\theta_{R}^{t}$) and $\theta_{L}^{b}$ ($\theta_{R}^{b}$).
They can be related with each other [12].

  * •

For the triplet $(X,T,B)_{L,R}$, we have the relations
$\tan\theta_{R}^{t}=\frac{m_{t}}{m_{T}}\tan\theta_{L}^{t}$ and
$\tan\theta_{R}^{b}=\frac{m_{b}}{m_{B}}\tan\theta_{L}^{b}$. $\theta_{L}^{t}$
and $\theta_{L}^{b}$ can be related through the identity $\sin
2\theta_{L}^{b}=\sqrt{2}\frac{m_{T}^{2}-m_{t}^{2}}{m_{B}^{2}-m_{b}^{2}}\sin
2\theta_{L}^{t}$. Thus, there is only one independent mixing angle
$\theta_{L}^{t}$.

  * •

For the triplet $(T,B,Y)_{L,R}$, we have the relations
$\tan\theta_{R}^{t}=\frac{m_{t}}{m_{T}}\tan\theta_{L}^{t}$ and
$\tan\theta_{R}^{b}=\frac{m_{b}}{m_{B}}\tan\theta_{L}^{b}$. $\theta_{L}^{t}$
and $\theta_{L}^{b}$ can be related through the identity $\sin
2\theta_{L}^{b}=\frac{1}{\sqrt{2}}\frac{m_{T}^{2}-m_{t}^{2}}{m_{B}^{2}-m_{b}^{2}}\sin
2\theta_{L}^{t}$. Thus, there is only one independent mixing angle
$\theta_{L}^{t}$.

  * •

For the doublet $(X,T)_{L,R}$, we have the relation
$\tan\theta_{L}^{t}=\frac{m_{t}}{m_{T}}\tan\theta_{R}^{t}$. Thus, there is
only one independent mixing angle $\theta_{R}^{t}$.

  * •

For the doublet $(T,B)_{L,R}$, we have the relations
$\tan\theta_{L}^{t}=\frac{m_{t}}{m_{T}}\tan\theta_{R}^{t}$ and
$\tan\theta_{L}^{b}=\frac{m_{b}}{m_{B}}\tan\theta_{R}^{b}$. Thus, there are
two independent mixing angles $\theta_{R}^{t}$ and $\theta_{R}^{b}$.

For the doublet $(X,T)_{L,R}$, the constraints on FCNY couplings from di-Higgs
production are similar to those in the singlet $T_{L,R}$ case. Compared to the
singlet $T_{L,R}$, there are extra $BBh,Bbh$ type Yukawa interactions for the
doublet $(T,B)_{L,R}$ and triplets $(X,T,B)_{L,R},(T,B,Y)_{L,R}$. Thus, it is
expected that the constraints on FCNY interactions are looser.

##  VI Summary and conclusions

Top partners are well motivated in many new physics models and FCNY
interactions can appear between top quark and the new heavy quark. To unveil
the nature of flavor structure and EWSB, it is important to probe the FCNY
interactions. However, it is challenging to constrain the $Tth$ coupling at
both current and future experiments directly.

In this paper, we have introduced a simplified model and summarized the main
constraints from theoretical and experimental viewpoints first. Then we
calculate the amplitude of di-Higgs production. After choosing
$m_{T}=400~{}\mathrm{GeV},s_{L}=0.2$ and $m_{T}=800~{}\mathrm{GeV},s_{L}=0.1$
as benchmark points, we evaluate the numerical cross sections. It is found
that the present constraints from di-Higgs production have already surpassed
the unitarity bound because of the $(y_{L,R}^{tT})^{4}$ behavior in di-Higgs
production cross section. For the case of $m_{T}=400~{}\mathrm{GeV}$ and
$s_{L}=0.2$, $\mathrm{Re}y_{L,R}^{tT}$ and $\mathrm{Im}y_{L,R}^{tT}$ are
expected to be bounded in the range $(-2,2)$, and even $(-0.4,0.4)$ in some
scenarios at HL-LHC with $2\sigma$ CL roughly. For the case of
$m_{T}=800~{}\mathrm{GeV}$ and $s_{L}=0.1$, $\mathrm{Re}y_{L,R}^{tT}$ and
$\mathrm{Im}y_{L,R}^{tT}$ are expected to be bounded in the range $(-3,3)$,
and even $(-0.5,0.5)$ in some scenarios at HL-LHC with $2\sigma$ CL roughly.
The value of $\delta_{hhh}$ can have significant effects on the constraints of
$y_{L,R}^{tT}$. Simply speaking, larger $\delta_{hhh}$ leads to looser
constraints on $y_{L,R}^{tT}$, because there can be more cancellation between
the triangle and box diagrams. Finally, we find that the top quark EDM can
give stronger bounds of $y_{L,R}^{tT}$ in the off-axis regions for some
scenarios.

## Acknowledgements

We would like to thank Gang Li, Zhao Li, Ying-nan Mao, and Hao Zhang for
helpful discussions. We also thank Olivier Mattelaer for MadGraph program
discussions through the launchpad platform.

## References

  * [1] S. L. Glashow.  Partial Symmetries of Weak Interactions.  Nucl. Phys., 22:579–588, 1961.
  * [2] Steven Weinberg.  A Model of Leptons.  Phys. Rev. Lett., 19:1264–1266, 1967.
  * [3] Abdus Salam.  Weak and Electromagnetic Interactions.  Conf. Proc., C680519:367–377, 1968.
  * [4] M. Tanabashi et al.  Review of Particle Physics.  Phys. Rev., D98(3):030001, 2018.
  * [5] F. Englert and R. Brout.  Broken Symmetry and the Mass of Gauge Vector Mesons.  Phys. Rev. Lett., 13:321–323, 1964.
  * [6] Peter W. Higgs.  Broken symmetries, massless particles and gauge fields.  Phys. Lett., 12:132–133, 1964.
  * [7] G. S. Guralnik, C. R. Hagen, and T. W. B. Kibble.  Global Conservation Laws and Massless Particles.  Phys. Rev. Lett., 13:585–587, 1964.
  * [8] T. W. B. Kibble.  Symmetry breaking in nonAbelian gauge theories.  Phys. Rev., 155:1554–1561, 1967.
  * [9] Georges Aad et al.  Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC.  Phys. Lett., B716:1–29, 2012.
  * [10] Serguei Chatrchyan et al.  Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC.  Phys. Lett., B716:30–61, 2012.
  * [11] J. A. Aguilar-Saavedra.  Identifying top partners at LHC.  JHEP, 11:030, 2009.
  * [12] J. A. Aguilar-Saavedra, R. Benbrik, S. Heinemeyer, and M. Pérez-Victoria.  Handbook of vectorlike quarks: Mixing and single production.  Phys. Rev., D88(9):094010, 2013.
  * [13] S. Dittmaier et al.  Handbook of LHC Higgs Cross Sections: 1. Inclusive Observables.  2011\.
  * [14] S. Dittmaier et al.  Handbook of LHC Higgs Cross Sections: 2. Differential Distributions.  2012\.
  * [15] J R Andersen et al.  Handbook of LHC Higgs Cross Sections: 3. Higgs Properties.  2013\.
  * [16] D. de Florian et al.  Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector.  2016\.
  * [17] Shi-Ping He.  Higgs boson to $\gamma Z$ decay as a probe of flavor-changing neutral Yukawa couplings.  Phys. Rev. D, 102(7):075035, 2020.
  * [18] S. Dawson and E. Furlan.  A Higgs Conundrum with Vector Fermions.  Phys. Rev., D86:015021, 2012.
  * [19] Andrea De Simone, Oleksii Matsedonskyi, Riccardo Rattazzi, and Andrea Wulzer.  A First Top Partner Hunter’s Guide.  JHEP, 04:004, 2013.
  * [20] F. del Aguila, M. Perez-Victoria, and Jose Santiago.  Effective description of quark mixing.  Phys. Lett., B492:98–106, 2000.
  * [21] F. del Aguila, M. Perez-Victoria, and Jose Santiago.  Observable contributions of new exotic quarks to quark mixing.  JHEP, 09:011, 2000.
  * [22] J. A. Aguilar-Saavedra.  Effects of mixing with quark singlets.  Phys. Rev., D67:035003, 2003.  [Erratum: Phys. Rev.D69,099901(2004)].
  * [23] Matthew J. Dolan, J. L. Hewett, M. Krämer, and T. G. Rizzo.  Simplified Models for Higgs Physics: Singlet Scalar and Vector-like Quark Phenomenology.  JHEP, 07:039, 2016.
  * [24] Jeong Han Kim and Ian M. Lewis.  Loop Induced Single Top Partner Production and Decay at the LHC.  JHEP, 05:095, 2018.
  * [25] J. A. Aguilar-Saavedra, D. E. López-Fogliani, and C. Muñoz.  Novel signatures for vector-like quarks.  JHEP, 06:095, 2017.
  * [26] Margarete Muhlleitner, Marco O. P. Sampaio, Rui Santos, and Jonas Wittbrodt.  The N2HDM under Theoretical and Experimental Scrutiny.  JHEP, 03:094, 2017.
  * [27] Kingman Cheung, Shi-Ping He, Ying-nan Mao, Po-Yan Tseng, and Chen Zhang.  Phenomenology of a little Higgs pseudoaxion.  Phys. Rev., D98(7):075023, 2018.
  * [28] Giacomo Cacciapaglia, Thomas Flacke, Myeonghun Park, and Mengchao Zhang.  Exotic decays of top partners: mind the search gap.  Phys. Lett. B, 798:135015, 2019.
  * [29] Mathieu Buchkremer, Giacomo Cacciapaglia, Aldo Deandrea, and Luca Panizzi.  Model Independent Framework for Searches of Top Partners.  Nucl. Phys., B876:376–417, 2013.
  * [30] Shinya Kanemura, Yasuhiro Okada, Eibun Senaha, and C.-P. Yuan.  Higgs coupling constants as a probe of new physics.  Phys. Rev. D, 70:115002, 2004.
  * [31] Shi-Ping He and Shou-hua Zhu.  One-Loop Radiative Correction to the Triple Higgs Coupling in the Higgs Singlet Model.  Phys. Lett., B764:31–37, 2017.
  * [32] Shinya Kanemura, Mariko Kikuchi, and Kei Yagyu.  One-loop corrections to the Higgs self-couplings in the singlet extension.  Nucl. Phys., B917:154–177, 2017.
  * [33] Abdesslam Arhrib, Rachid Benbrik, Jaouad El Falaki, and Adil Jueid.  Radiative corrections to the Triple Higgs Coupling in the Inert Higgs Doublet Model.  JHEP, 12:007, 2015.
  * [34] Shinya Kanemura, Mariko Kikuchi, Kodai Sakurai, and Kei Yagyu.  H-COUP: a program for one-loop corrected Higgs boson couplings in non-minimal Higgs sectors.  Comput. Phys. Commun., 233:134–144, 2018.
  * [35] Cheng-Wei Chiang, An-Li Kuo, and Kei Yagyu.  One-loop renormalized Higgs boson vertices in the Georgi-Machacek model.  Phys. Rev. D, 98(1):013008, 2018.
  * [36] Johannes Braathen and Shinya Kanemura.  On two-loop corrections to the Higgs trilinear coupling in models with extended scalar sectors.  Phys. Lett. B, 796:38–46, 2019.
  * [37] Christoph Englert and Joerg Jaeckel.  Probing the Symmetric Higgs Portal with Di-Higgs Boson Production.  Phys. Rev. D, 100(9):095017, 2019.
  * [38] Shinya Kanemura, Mariko Kikuchi, Kentarou Mawatari, Kodai Sakurai, and Kei Yagyu.  H-COUP Version 2: a program for one-loop corrected Higgs boson decays in non-minimal Higgs sectors.  Comput. Phys. Commun., 257:107512, 2020.
  * [39] Christoph Englert, Joerg Jaeckel, Michael Spannowsky, and Panagiotis Stylianou.  Power meets Precision to explore the Symmetric Higgs Portal.  Phys. Lett. B, 806:135526, 2020.
  * [40] Michael E. Peskin and Tatsu Takeuchi.  A New constraint on a strongly interacting Higgs sector.  Phys. Rev. Lett., 65:964–967, 1990.
  * [41] Michael E. Peskin and Tatsu Takeuchi.  Estimation of oblique electroweak corrections.  Phys. Rev., D46:381–409, 1992.
  * [42] L. Lavoura and Joao P. Silva.  The Oblique corrections from vector - like singlet and doublet quarks.  Phys. Rev., D47:2046–2057, 1993.
  * [43] Chien-Yi Chen, S. Dawson, and Elisabetta Furlan.  Vectorlike fermions and Higgs effective field theory revisited.  Phys. Rev., D96(1):015006, 2017.
  * [44] Jacob Baron et al.  Order of Magnitude Smaller Limit on the Electric Dipole Moment of the Electron.  Science, 343:269–272, 2014.
  * [45] V. Andreev et al.  Improved limit on the electric dipole moment of the electron.  Nature, 562(7727):355–360, 2018.
  * [46] Jernej F. Kamenik, Michele Papucci, and Andreas Weiler.  Constraining the dipole moments of the top quark.  Phys. Rev. D, 85:071501, 2012.  [Erratum: Phys.Rev.D 88, 039903 (2013)].
  * [47] V. Cirigliano, W. Dekens, J. de Vries, and E. Mereghetti.  Is there room for CP violation in the top-Higgs sector?  Phys. Rev. D, 94(1):016002, 2016.
  * [48] V. Cirigliano, W. Dekens, J. de Vries, and E. Mereghetti.  Constraining the top-Higgs sector of the Standard Model Effective Field Theory.  Phys. Rev. D, 94(3):034031, 2016.
  * [49] E.W.Nigel Glover and J.J. van der Bij.  HIGGS BOSON PAIR PRODUCTION VIA GLUON FUSION.  Nucl. Phys. B, 309:282–294, 1988.
  * [50] Abdelhak Djouadi.  The Anatomy of electro-weak symmetry breaking. I: The Higgs boson in the standard model.  Phys. Rept., 457:1–216, 2008.
  * [51] Eri Asakawa, Daisuke Harada, Shinya Kanemura, Yasuhiro Okada, and Koji Tsumura.  Higgs boson pair production in new physics models at hadron, lepton, and photon colliders.  Phys. Rev. D, 82:115002, 2010.
  * [52] Matthew J. Dolan, Christoph Englert, and Michael Spannowsky.  New Physics in LHC Higgs boson pair production.  Phys. Rev. D, 87(5):055002, 2013.
  * [53] S. Dawson, A. Ismail, and Ian Low.  What’s in the loop? The anatomy of double Higgs production.  Phys. Rev. D, 91(11):115008, 2015.
  * [54] Hong-Jian He, Jing Ren, and Weiming Yao.  Probing new physics of cubic Higgs boson interaction via Higgs pair production at hadron colliders.  Phys. Rev. D, 93(1):015003, 2016.
  * [55] J. Alison et al.  Higgs Boson Pair Production at Colliders: Status and Perspectives.  In B. Di Micco, M. Gouzevitch, J. Mazzitelli, and C. Vernieri, editors, Double Higgs Production at Colliders, 9 2019.
  * [56] Florian Goertz, Andreas Papaefstathiou, Li Lin Yang, and José Zurita.  Higgs boson pair production in the D=6 extension of the SM.  JHEP, 04:167, 2015.
  * [57] Aleksandr Azatov, Roberto Contino, Giuliano Panico, and Minho Son.  Effective field theory analysis of double Higgs boson production via gluon fusion.  Phys. Rev. D, 92(3):035001, 2015.
  * [58] Chih-Ting Lu, Jung Chang, Kingman Cheung, and Jae Sik Lee.  An exploratory study of Higgs-boson pair production.  JHEP, 08:133, 2015.
  * [59] Qing-Hong Cao, Bin Yan, Dong-Ming Zhang, and Hao Zhang.  Resolving the Degeneracy in Single Higgs Production with Higgs Pair Production.  Phys. Lett. B, 752:285–290, 2016.
  * [60] Qing-Hong Cao, Gang Li, Bin Yan, Dong-Ming Zhang, and Hao Zhang.  Double Higgs production at the 14 TeV LHC and a 100 TeV $pp$ collider.  Phys. Rev. D, 96(9):095031, 2017.
  * [61] Gang Li, Ling-Xiao Xu, Bin Yan, and C.-P. Yuan.  Resolving the degeneracy in top quark Yukawa coupling with Higgs pair production.  Phys. Lett. B, 800:135070, 2020.
  * [62] Roberto Contino, Christophe Grojean, Mauro Moretti, Fulvio Piccinini, and Riccardo Rattazzi.  Strong Double Higgs Production at the LHC.  JHEP, 05:089, 2010.
  * [63] Roberto Contino, Margherita Ghezzi, Mauro Moretti, Giuliano Panico, Fulvio Piccinini, and Andrea Wulzer.  Anomalous Couplings in Double Higgs Production.  JHEP, 08:154, 2012.
  * [64] R. Grober, M. Muhlleitner, and M. Spira.  Higgs Pair Production at NLO QCD for CP-violating Higgs Sectors.  Nucl. Phys. B, 925:1–27, 2017.
  * [65] G. Buchalla, M. Capozi, A. Celis, G. Heinrich, and L. Scyboz.  Higgs boson pair production in non-linear Effective Field Theory with full $m_{t}$-dependence at NLO QCD.  JHEP, 09:057, 2018.
  * [66] Chien-Yi Chen, S. Dawson, and I.M. Lewis.  Exploring resonant di-Higgs boson production in the Higgs singlet model.  Phys. Rev. D, 91(3):035015, 2015.
  * [67] S. Dawson and I.M. Lewis.  NLO corrections to double Higgs boson production in the Higgs singlet model.  Phys. Rev. D, 92(9):094023, 2015.
  * [68] Ian M. Lewis and Matthew Sullivan.  Benchmarks for Double Higgs Production in the Singlet Extended Standard Model at the LHC.  Phys. Rev. D, 96(3):035037, 2017.
  * [69] Lan-Chun Lü, Chun Du, Yaquan Fang, Hong-Jian He, and Huijun Zhang.  Searching heavier Higgs boson via di-Higgs production at LHC Run-2.  Phys. Lett. B, 755:509–522, 2016.
  * [70] Stefania De Curtis, Stefano Moretti, Kei Yagyu, and Emine Yildirim.  Single and double SM-like Higgs boson production at future electron-positron colliders in composite 2HDMs.  Phys. Rev. D, 95(9):095026, 2017.
  * [71] Tadashi Kon, Takuto Nagura, Takahiro Ueda, and Kei Yagyu.  Double Higgs boson production at $e^{+}e^{-}$ colliders in the two-Higgs-doublet model.  Phys. Rev. D, 99(9):095027, 2019.
  * [72] Jing Ren, Rui-Qing Xiao, Maosen Zhou, Yaquan Fang, Hong-Jian He, and Weiming Yao.  LHC Search of New Higgs Boson via Resonant Di-Higgs Production with Decays into 4W.  JHEP, 06:090, 2018.
  * [73] Sally Dawson, Elisabetta Furlan, and Ian Lewis.  Unravelling an extended quark sector through multiple Higgs production?  Phys. Rev. D, 87(1):014007, 2013.
  * [74] Giacomo Cacciapaglia, Haiying Cai, Alexandra Carvalho, Aldo Deandrea, Thomas Flacke, Benjamin Fuks, Devdatta Majumder, and Hua-Sheng Shao.  Probing vector-like quark models with Higgs-boson pair production.  JHEP, 07:005, 2017.
  * [75] Kingman Cheung, Adil Jueid, Chih-Ting Lu, Jeonghyeon Song, and Yeo Woong Yoon.  Disentangling new physics effects on non-resonant Higgs boson pair production from gluon fusion.  3 2020.
  * [76] M. Gillioz, R. Grober, C. Grojean, M. Muhlleitner, and E. Salvioni.  Higgs Low-Energy Theorem (and its corrections) in Composite Models.  JHEP, 10:004, 2012.
  * [77] Ramona Grober, Margarete Muhlleitner, and Michael Spira.  Signs of Composite Higgs Pair Production at Next-to-Leading Order.  JHEP, 06:080, 2016.
  * [78] T. Plehn, M. Spira, and P.M. Zerwas.  Pair production of neutral Higgs particles in gluon-gluon collisions.  Nucl. Phys. B, 479:46–64, 1996.  [Erratum: Nucl.Phys.B 531, 655–655 (1998)].
  * [79] S. Dawson, S. Dittmaier, and M. Spira.  Neutral Higgs boson pair production at hadron colliders: QCD corrections.  Phys. Rev. D, 58:115012, 1998.
  * [80] Abdelhak Djouadi.  The Anatomy of electro-weak symmetry breaking. II. The Higgs bosons in the minimal supersymmetric model.  Phys. Rept., 459:1–241, 2008.
  * [81] Junjie Cao, Zhaoxia Heng, Liangliang Shang, Peihua Wan, and Jin Min Yang.  Pair Production of a 125 GeV Higgs Boson in MSSM and NMSSM at the LHC.  JHEP, 04:134, 2013.
  * [82] Philipp Basler, Sally Dawson, Christoph Englert, and Margarete Mühlleitner.  Showcasing HH production: Benchmarks for the LHC and HL-LHC.  Phys. Rev. D, 99(5):055048, 2019.
  * [83] D. Binosi, J. Collins, C. Kaufhold, and L. Theussl.  JaxoDraw: A Graphical user interface for drawing Feynman diagrams. Version 2.0 release notes.  Comput. Phys. Commun., 180:1709–1715, 2009.
  * [84] R. Mertig, M. Bohm, and Ansgar Denner.  FEYN CALC: Computer algebraic calculation of Feynman amplitudes.  Comput. Phys. Commun., 64:345–359, 1991.
  * [85] Vladyslav Shtabovenko, Rolf Mertig, and Frederik Orellana.  New Developments in FeynCalc 9.0.  Comput. Phys. Commun., 207:432–444, 2016.
  * [86] Ding Yu Shao, Chong Sheng Li, Hai Tao Li, and Jian Wang.  Threshold resummation effects in Higgs boson pair production at the LHC.  JHEP, 07:169, 2013.
  * [87] Daniel de Florian and Javier Mazzitelli.  Higgs Boson Pair Production at Next-to-Next-to-Leading Order in QCD.  Phys. Rev. Lett., 111:201801, 2013.
  * [88] Daniel de Florian and Javier Mazzitelli.  Higgs pair production at next-to-next-to-leading logarithmic accuracy at the LHC.  JHEP, 09:053, 2015.
  * [89] Giuseppe Degrassi, Pier Paolo Giardino, and Ramona Gröber.  On the two-loop virtual QCD corrections to Higgs boson pair production in the Standard Model.  Eur. Phys. J. C, 76(7):411, 2016.
  * [90] S. Borowka, N. Greiner, G. Heinrich, S.P. Jones, M. Kerner, J. Schlenk, U. Schubert, and T. Zirke.  Higgs Boson Pair Production in Gluon Fusion at Next-to-Leading Order with Full Top-Quark Mass Dependence.  Phys. Rev. Lett., 117(1):012001, 2016.  [Erratum: Phys.Rev.Lett. 117, 079901 (2016)].
  * [91] Massimiliano Grazzini, Gudrun Heinrich, Stephen Jones, Stefan Kallweit, Matthias Kerner, Jonas M. Lindert, and Javier Mazzitelli.  Higgs boson pair production at NNLO with top quark mass effects.  JHEP, 05:059, 2018.
  * [92] Julien Baglio, Francisco Campanario, Seraina Glaus, Margarete Mühlleitner, Michael Spira, and Juraj Streicher.  Gluon fusion into Higgs pairs at NLO QCD and the top mass scheme.  Eur. Phys. J. C, 79(6):459, 2019.
  * [93] Long-Bin Chen, Hai Tao Li, Hua-Sheng Shao, and Jian Wang.  Higgs boson pair production via gluon fusion at N3LO in QCD.  Phys. Lett. B, 803:135292, 2020.
  * [94] Long-Bin Chen, Hai Tao Li, Hua-Sheng Shao, and Jian Wang.  The gluon-fusion production of Higgs boson pair: N3LO QCD corrections and top-quark mass effects.  JHEP, 03:072, 2020.
  * [95] Alexandra Carvalho, Martino Dall’Osso, Tommaso Dorigo, Florian Goertz, Carlo A. Gottardo, and Mia Tosi.  Higgs Pair Production: Choosing Benchmarks With Cluster Analysis.  JHEP, 04:126, 2016.
  * [96] Alexandra Carvalho, Martino Dall’Osso, Pablo De Castro Manzano, Tommaso Dorigo, Florian Goertz, Maxime Gouzevich, and Mia Tosi.  Analytical parametrization and shape classification of anomalous HH production in the EFT approach.  7 2016.
  * [97] Adam Alloul, Neil D. Christensen, Céline Degrande, Claude Duhr, and Benjamin Fuks.  FeynRules 2.0 - A complete toolbox for tree-level phenomenology.  Comput. Phys. Commun., 185:2250–2300, 2014.
  * [98] Celine Degrande, Claude Duhr, Benjamin Fuks, David Grellscheid, Olivier Mattelaer, and Thomas Reiter.  UFO - The Universal FeynRules Output.  Comput. Phys. Commun., 183:1201–1214, 2012.
  * [99] Thomas Hahn.  Generating Feynman diagrams and amplitudes with FeynArts 3.  Comput. Phys. Commun., 140:418–431, 2001.
  * [100] Celine Degrande.  Automatic evaluation of UV and R2 terms for beyond the Standard Model Lagrangians: a proof-of-principle.  Comput. Phys. Commun., 197:239–262, 2015.
  * [101] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro.  The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations.  JHEP, 07:079, 2014.
  * [102] Valentin Hirschi and Olivier Mattelaer.  Automated event generation for loop-induced processes.  JHEP, 10:146, 2015.
  * [103] Albert M Sirunyan et al.  Combination of searches for Higgs boson pair production in proton-proton collisions at $\sqrt{s}=$ 13 TeV.  Phys. Rev. Lett., 122(12):121803, 2019.
  * [104] Georges Aad et al.  Combination of searches for Higgs boson pairs in $pp$ collisions at $\sqrt{s}=$13 TeV with the ATLAS detector.  Phys. Lett. B, 800:135103, 2020.
  * [105] M. Cepeda et al.  Report from Working Group 2: Higgs Physics at the HL-LHC and HE-LHC, volume 7, pages 221–584.  12 2019.
  * [106] T. Hahn and M. Perez-Victoria.  Automatized one loop calculations in four-dimensions and D-dimensions.  Comput. Phys. Commun., 118:153–165, 1999.

## Appendix

##  Appendix A Asymptotic behaviors of the loop functions

###  A.1 The shorthand notations of $C_{0}$ and $D_{0}$ functions

The definitions of $C_{0}$ and $D_{0}$ function related with pure top quark
loops are given as

$\displaystyle C_{0}^{t}(\hat{s})\equiv
C_{0}(0,0,\hat{s},m_{t}^{2},m_{t}^{2},m_{t}^{2}),$ $\displaystyle
C_{0}^{t}(m_{h}^{2})\equiv
C_{0}(m_{h}^{2},m_{h}^{2},\hat{s},m_{t}^{2},m_{t}^{2},m_{t}^{2}),$
$\displaystyle C_{0}^{t}(\hat{t})\equiv
C_{0}(0,m_{h}^{2},\hat{t},m_{t}^{2},m_{t}^{2},m_{t}^{2}),$ $\displaystyle
C_{0}^{t}(\hat{u})\equiv
C_{0}(0,m_{h}^{2},\hat{u},m_{t}^{2},m_{t}^{2},m_{t}^{2}),$ $\displaystyle
D_{0}^{t}(\hat{t},\hat{s})\equiv
D_{0}(m_{h}^{2},0,0,m_{h}^{2},\hat{t},\hat{s},m_{t}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2}),$
$\displaystyle D_{0}^{t}(\hat{u},\hat{s})\equiv
D_{0}(m_{h}^{2},0,0,m_{h}^{2},\hat{u},\hat{s},m_{t}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2}),$
$\displaystyle D_{0}^{t}(\hat{t},\hat{u})\equiv
D_{0}(m_{h}^{2},0,m_{h}^{2},0,\hat{t},\hat{u},m_{t}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2}).$
(40)

The definitions of $C_{0}$ and $D_{0}$ function related with pure $T$ quark
loops are given as

$\displaystyle C_{0}^{T}(\hat{s})\equiv
C_{0}(0,0,\hat{s},m_{T}^{2},m_{T}^{2},m_{T}^{2}),$ $\displaystyle
C_{0}^{T}(m_{h}^{2})\equiv
C_{0}(m_{h}^{2},m_{h}^{2},\hat{s},m_{T}^{2},m_{T}^{2},m_{T}^{2}),$
$\displaystyle C_{0}^{T}(\hat{t})\equiv
C_{0}(0,m_{h}^{2},\hat{t},m_{T}^{2},m_{T}^{2},m_{T}^{2}),$ $\displaystyle
C_{0}^{T}(\hat{u})\equiv
C_{0}(0,m_{h}^{2},\hat{u},m_{T}^{2},m_{T}^{2},m_{T}^{2}),$ $\displaystyle
D_{0}^{T}(\hat{t},\hat{s})\equiv
D_{0}(m_{h}^{2},0,0,m_{h}^{2},\hat{t},\hat{s},m_{T}^{2},m_{T}^{2},m_{T}^{2},m_{T}^{2}),$
$\displaystyle D_{0}^{T}(\hat{u},\hat{s})\equiv
D_{0}(m_{h}^{2},0,0,m_{h}^{2},\hat{u},\hat{s},m_{T}^{2},m_{T}^{2},m_{T}^{2},m_{T}^{2}),$
$\displaystyle D_{0}^{T}(\hat{t},\hat{u})\equiv
D_{0}(m_{h}^{2},0,m_{h}^{2},0,\hat{t},\hat{u},m_{T}^{2},m_{T}^{2},m_{T}^{2},m_{T}^{2}).$
(41)

The definitions of $C_{0}$ and $D_{0}$ function related with mixed $t$ and $T$
quark loops are given as

$\displaystyle C_{0}^{tT}(m_{h}^{2})\equiv
C_{0}(m_{h}^{2},m_{h}^{2},\hat{s},m_{t}^{2},m_{T}^{2},m_{t}^{2}),$
$\displaystyle C_{0}^{tT}(\hat{t})\equiv
C_{0}(0,m_{h}^{2},\hat{t},m_{t}^{2},m_{t}^{2},m_{T}^{2}),$ $\displaystyle
C_{0}^{tT}(\hat{u})\equiv
C_{0}(0,m_{h}^{2},\hat{u},m_{t}^{2},m_{t}^{2},m_{T}^{2}),$ $\displaystyle
D_{0}^{tT}(\hat{t},\hat{s})\equiv
D_{0}(m_{h}^{2},0,0,m_{h}^{2},\hat{t},\hat{s},m_{T}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2}),$
$\displaystyle D_{0}^{tT}(\hat{u},\hat{s})\equiv
D_{0}(m_{h}^{2},0,0,m_{h}^{2},\hat{u},\hat{s},m_{T}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2}),$
$\displaystyle D_{0}^{tT}(\hat{t},\hat{u})\equiv
D_{0}(m_{h}^{2},0,m_{h}^{2},0,\hat{t},\hat{u},m_{t}^{2},m_{T}^{2},m_{T}^{2},m_{t}^{2}),$
(42)

and

$\displaystyle C_{0}^{Tt}(m_{h}^{2})\equiv
C_{0}(m_{h}^{2},m_{h}^{2},\hat{s},m_{T}^{2},m_{t}^{2},m_{T}^{2}),$
$\displaystyle C_{0}^{Tt}(\hat{t})\equiv
C_{0}(0,m_{h}^{2},\hat{t},m_{T}^{2},m_{T}^{2},m_{t}^{2}),$ $\displaystyle
C_{0}^{Tt}(\hat{u})\equiv
C_{0}(0,m_{h}^{2},\hat{u},m_{T}^{2},m_{T}^{2},m_{t}^{2}),$ $\displaystyle
D_{0}^{Tt}(\hat{t},\hat{s})\equiv
D_{0}(m_{h}^{2},0,0,m_{h}^{2},\hat{t},\hat{s},m_{t}^{2},m_{T}^{2},m_{T}^{2},m_{T}^{2}),$
$\displaystyle D_{0}^{Tt}(\hat{u},\hat{s})\equiv
D_{0}(m_{h}^{2},0,0,m_{h}^{2},\hat{u},\hat{s},m_{t}^{2},m_{T}^{2},m_{T}^{2},m_{T}^{2}),$
$\displaystyle D_{0}^{Tt}(\hat{t},\hat{u})\equiv
D_{0}(m_{h}^{2},0,m_{h}^{2},0,\hat{t},\hat{u},m_{T}^{2},m_{t}^{2},m_{t}^{2},m_{T}^{2}).$
(43)

As a matter of fact, we have the relation
$D_{0}^{tT}(\hat{t},\hat{u})=D_{0}^{Tt}(\hat{t},\hat{u})$.

###  A.2 Heavy quark expansion of $C_{0}$ function

$C_{0}$ function is defined as

$\displaystyle
C_{0}(k_{1}^{2},k_{12}^{2},k_{2}^{2},m_{0}^{2},m_{1}^{2},m_{2}^{2})$
$\displaystyle\equiv$ $\displaystyle\frac{(2\pi\mu)^{4-D}}{i\pi^{2}}\int
d^{D}q\frac{1}{(q^{2}-m_{0}^{2})[(q+k_{1})^{2}-m_{1}^{2}][(q+k_{2})^{2}-m_{2}^{2}]}$
$\displaystyle=$
$\displaystyle-\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydz\frac{\delta(x+y+z-1)}{xm_{0}^{2}+ym_{1}^{2}+zm_{2}^{2}-xyk_{1}^{2}-xzk_{2}^{2}-yzk_{12}^{2}},$
(44)

where $k_{12}\equiv k_{1}-k_{2}$ and $D$ is the dimension of space time. When
the three internal masses are all equal, $C_{0}$ function can be expanded as
[86]

$\displaystyle
C_{0}(k_{1}^{2},k_{12}^{2},k_{2}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2})$
$\displaystyle=$
$\displaystyle-\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydz\frac{\delta(x+y+z-1)}{m_{t}^{2}-xyk_{1}^{2}-xzk_{2}^{2}-yzk_{12}^{2}}$
$\displaystyle=$
$\displaystyle-\frac{1}{2m_{t}^{2}}-\frac{k_{1}^{2}+k_{2}^{2}+k_{12}^{2}}{24m_{t}^{4}}-\frac{k_{1}^{4}+k_{2}^{4}+k_{12}^{4}+k_{1}^{2}k_{2}^{2}+k_{1}^{2}k_{12}^{2}+k_{2}^{2}k_{12}^{2}}{180m_{t}^{6}}+\mathcal{O}(\frac{k^{6}}{m_{t}^{8}}).$
(45)

Especially, we have the following results:

$\displaystyle
C_{0}^{t}(\hat{s})\approx-\frac{1}{2m_{t}^{2}}(1+\frac{\hat{s}}{12m_{t}^{2}}+\frac{\hat{s}^{2}}{90m_{t}^{4}}),$
$\displaystyle
C_{0}^{t}(m_{h}^{2})\approx-\frac{1}{2m_{t}^{2}}(1+\frac{2m_{h}^{2}+\hat{s}}{12m_{t}^{2}}+\frac{3m_{h}^{4}+2m_{h}^{2}\hat{s}+\hat{s}^{2}}{90m_{t}^{4}}),$
$\displaystyle
C_{0}^{t}(\hat{t})\approx-\frac{1}{2m_{t}^{2}}(1+\frac{m_{h}^{2}+\hat{t}}{12m_{t}^{2}}+\frac{m_{h}^{4}+m_{h}^{2}\hat{t}+\hat{t}^{2}}{90m_{t}^{4}}).$
(46)

The expansion of $C_{0}^{t}(\hat{u})$ functions can be obtained when replacing
the $\hat{t}$ in $C_{0}^{t}(\hat{t})$ by $\hat{u}$. The expansion of
$C_{0}^{T}$ functions can be obtained when replacing the $m_{t}$ in
$C_{0}^{t}$ by $m_{T}$.

When the first two internal masses are equal, $C_{0}$ function can be expanded
as

$\displaystyle
C_{0}(k_{1}^{2},k_{12}^{2},k_{2}^{2},m_{t}^{2},m_{t}^{2},m_{T}^{2})$
$\displaystyle=$
$\displaystyle-\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydz\frac{\delta(x+y+z-1)}{(x+y)m_{t}^{2}+zm_{T}^{2}-xyk_{1}^{2}-xzk_{2}^{2}-yzk_{12}^{2}}$
$\displaystyle=$
$\displaystyle-\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydz\frac{\delta(x+y+z-1)}{(x+y)m_{t}^{2}+zm_{T}^{2}}[1+\frac{xyk_{1}^{2}+xzk_{2}^{2}+yzk_{12}^{2}}{(x+y)m_{t}^{2}+zm_{T}^{2}}+\frac{(xyk_{1}^{2}+xzk_{2}^{2}+yzk_{12}^{2})^{2}}{((x+y)m_{t}^{2}+zm_{T}^{2})^{2}}]+\mathcal{O}(\frac{k^{6}}{m_{t,T}^{8}})$
$\displaystyle\approx$ $\displaystyle\frac{1+\log
r_{tT}^{2}-r_{tT}^{2}}{m_{T}^{2}(1-r_{tT}^{2})^{2}}-\frac{2+6r_{tT}^{2}\log
r_{tT}^{2}+3r_{tT}^{2}-6r_{tT}^{4}+r_{tT}^{6}}{12m_{T}^{4}r_{tT}^{2}(1-r_{tT}^{2})^{4}}k_{1}^{2}$
$\displaystyle+\frac{5+2(1+2r_{tT}^{2})\log
r_{tT}^{2}-4r_{tT}^{2}-r_{tT}^{4}}{4m_{T}^{4}(1-r_{tT}^{2})^{4}}(k_{2}^{2}+k_{12}^{2})$
$\displaystyle-\frac{3-30r_{tT}^{2}-20r_{tT}^{4}(1+3\log
r_{tT}^{2})+60r_{tT}^{6}-15r_{tT}^{8}+2r_{tT}^{10}}{180m_{T}^{6}r_{tT}^{4}(1-r_{tT}^{2})^{6}}k_{1}^{4}$
$\displaystyle+\frac{10+9r_{tT}^{2}+3(1+6r_{tT}^{2}+3r_{tT}^{4})\log
r_{tT}^{2}-18r_{tT}^{4}-r_{tT}^{6}}{9m_{T}^{6}(1-r_{tT}^{2})^{6}}(k_{12}^{4}+k_{2}^{2}k_{12}^{2}+k_{2}^{4})$
$\displaystyle-\frac{3+44r_{tT}^{2}+12r_{tT}^{2}(2+3r_{tT}^{2})\log
r_{tT}^{2}-36r_{tT}^{4}-12r_{tT}^{6}+r_{tT}^{8}}{36m_{T}^{6}r_{tT}^{2}(1-r_{tT}^{2})^{6}}k_{1}^{2}(k_{2}^{2}+k_{12}^{2}).$
(47)

When the first and third internal masses are equal, $C_{0}$ function can be
correlated with the first two mass equal cases through the following relations

$\displaystyle
C_{0}(k_{1}^{2},k_{12}^{2},k_{2}^{2},m_{t}^{2},m_{T}^{2},m_{t}^{2})=C_{0}(k_{2}^{2},k_{12}^{2},k_{1}^{2},m_{t}^{2},m_{t}^{2},m_{T}^{2}),$
$\displaystyle
C_{0}(k_{1}^{2},k_{12}^{2},k_{2}^{2},m_{T}^{2},m_{t}^{2},m_{T}^{2})=C_{0}(k_{2}^{2},k_{12}^{2},k_{1}^{2},m_{T}^{2},m_{T}^{2},m_{t}^{2}).$
(48)

Especially, we have the following results

$\displaystyle C_{0}^{tT}(\hat{t})\approx\frac{1}{m_{T}^{2}}\cdot\frac{1+\log
r_{tT}^{2}-r_{tT}^{2}}{(1-r_{tT}^{2})^{2}}+\frac{m_{h}^{2}+\hat{t}}{m_{T}^{4}}\cdot\frac{5+2(1+2r_{tT}^{2})\log
r_{tT}^{2}-4r_{tT}^{2}-r_{tT}^{4}}{4(1-r_{tT}^{2})^{4}}$
$\displaystyle+\frac{m_{h}^{4}+m_{h}^{2}\hat{t}+\hat{t}^{2}}{m_{T}^{6}}\cdot\frac{10+3(1+6r_{tT}^{2}+3r_{tT}^{4})\log
r_{tT}^{2}+9r_{tT}^{2}-18r_{tT}^{4}-r_{tT}^{6}}{9(1-r_{tT}^{2})^{6}},$
$\displaystyle
C_{0}^{tT}(m_{h}^{2})\approx\frac{1}{m_{T}^{2}}\cdot\frac{1+\log
r_{tT}^{2}-r_{tT}^{2}}{(1-r_{tT}^{2})^{2}}-\frac{\hat{s}}{m_{T}^{4}}\cdot\frac{2+6r_{tT}^{2}\log
r_{tT}^{2}+3r_{tT}^{2}-6r_{tT}^{4}+r_{tT}^{6}}{12r_{tT}^{2}(1-r_{tT}^{2})^{4}}$
$\displaystyle+\frac{m_{h}^{2}}{m_{T}^{4}}\cdot\frac{5+2(1+2r_{tT}^{2})\log
r_{tT}^{2}-4r_{tT}^{2}-r_{tT}^{4}}{2(1-r_{tT}^{2})^{4}}+\frac{m_{h}^{4}}{m_{T}^{6}}\cdot\frac{10+3(1+6r_{tT}^{2}+3r_{tT}^{4})\log
r_{tT}^{2}+9r_{tT}^{2}-18r_{tT}^{4}-r_{tT}^{6}}{3(1-r_{tT}^{2})^{6}}$
$\displaystyle-\frac{m_{h}^{2}\hat{s}}{m_{T}^{6}}\cdot\frac{3+12r_{tT}^{2}(2+3r_{tT}^{2})\log
r_{tT}^{2}+44r_{tT}^{2}-36r_{tT}^{4}-12r_{tT}^{6}+r_{tT}^{8}}{18r_{tT}^{2}(1-r_{tT}^{2})^{6}}$
$\displaystyle-\frac{\hat{s}^{2}}{m_{T}^{6}}\cdot\frac{3-30r_{tT}^{2}-60r_{tT}^{4}\log
r_{tT}^{2}-20r_{tT}^{4}+60r_{tT}^{6}-15r_{tT}^{8}+2r_{tT}^{10}}{180r_{tT}^{4}(1-r_{tT}^{2})^{6}}.$
(49)

Keeping the terms up to $\mathcal{O}(\frac{1}{m_{T}^{4}})$ and considering the
$\log r_{tT}^{2}$ enhanced terms, they can be simplified as

$\displaystyle C_{0}^{tT}(\hat{t})\approx\frac{1}{m_{T}^{2}}\big{[}1+\log
r_{tT}^{2}+r_{tT}^{2}(1+2\log
r_{tT}^{2})+\frac{m_{h}^{2}+\hat{t}}{4m_{T}^{2}}(5+2\log r_{tT}^{2})\big{]},$
$\displaystyle C_{0}^{tT}(m_{h}^{2})\approx\frac{1}{m_{T}^{2}}\big{[}1+\log
r_{tT}^{2}+r_{tT}^{2}(1+2\log
r_{tT}^{2})-\frac{\hat{s}}{12m_{T}^{2}}(\frac{2}{r_{tT}^{2}}+11+6\log
r_{tT}^{2})+\frac{m_{h}^{2}}{2m_{T}^{2}}(5+2\log r_{tT}^{2})\big{]}.$ (50)

Similarly, we can get the following results

$\displaystyle
C_{0}^{Tt}(\hat{t})\approx-\frac{1}{m_{T}^{2}}\cdot\frac{1+r_{tT}^{2}\log
r_{tT}^{2}-r_{tT}^{2}}{(1-r_{tT}^{2})^{2}}-\frac{m_{h}^{2}+\hat{t}}{m_{T}^{4}}\cdot\frac{1+2r_{tT}^{2}(2+r_{tT}^{2})\log
r_{tT}^{2}+4r_{tT}^{2}-5r_{tT}^{4}}{4(1-r_{tT}^{2})^{4}}$
$\displaystyle-\frac{m_{h}^{4}+m_{h}^{2}\hat{t}+\hat{t}^{2}}{m_{T}^{6}}\cdot\frac{1+3r_{tT}^{2}(3+6r_{tT}^{2}+r_{tT}^{4})\log
r_{tT}^{2}+18r_{tT}^{2}-9r_{tT}^{4}-10r_{tT}^{6}}{9(1-r_{tT}^{2})^{6}},$
$\displaystyle
C_{0}^{Tt}(m_{h}^{2})\approx-\frac{1}{m_{T}^{2}}\cdot\frac{1+r_{tT}^{2}\log
r_{tT}^{2}-r_{tT}^{2}}{(1-r_{tT}^{2})^{2}}-\frac{\hat{s}}{m_{T}^{4}}\cdot\frac{1-6r_{tT}^{2}-6r_{tT}^{4}\log
r_{tT}^{2}+3r_{tT}^{4}+2r_{tT}^{6}}{12(1-r_{tT}^{2})^{4}}$
$\displaystyle-\frac{m_{h}^{2}}{2m_{T}^{4}}\cdot\frac{1+2r_{tT}^{2}(2+r_{tT}^{2})\log
r_{tT}^{2}+4r_{tT}^{2}-5r_{tT}^{4}}{(1-r_{tT}^{2})^{4}}-\frac{m_{h}^{4}}{m_{T}^{6}}\cdot\frac{1+3r_{tT}^{2}(3+6r_{tT}^{2}+r_{tT}^{4})\log
r_{tT}^{2}+18r_{tT}^{2}-9r_{tT}^{4}-10r_{tT}^{6}}{3(1-r_{tT}^{2})^{6}}$
$\displaystyle-\frac{m_{h}^{2}\hat{s}}{m_{T}^{6}}\cdot\frac{1-12r_{tT}^{2}-12r_{tT}^{4}(3+2r_{tT}^{2})\log
r_{tT}^{2}-36r_{tT}^{4}+44r_{tT}^{6}+3r_{tT}^{8}}{18(1-r_{tT}^{2})^{6}}$
$\displaystyle-\frac{\hat{s}^{2}}{m_{T}^{6}}\cdot\frac{2-15r_{tT}^{2}+60r_{tT}^{4}+60r_{tT}^{6}\log
r_{tT}^{2}-20r_{tT}^{6}-30r_{tT}^{8}+3r_{tT}^{10}}{180(1-r_{tT}^{2})^{6}}.$
(51)

Keeping the terms up to $\mathcal{O}(\frac{1}{m_{T}^{4}})$ and considering the
$\log r_{tT}^{2}$ enhanced terms, they can be simplified as

$\displaystyle
C_{0}^{Tt}(\hat{t})\approx-\frac{1}{m_{T}^{2}}\big{[}1+r_{tT}^{2}(1+\log
r_{tT}^{2})+\frac{m_{h}^{2}+\hat{t}}{4m_{T}^{2}}\big{]},$ $\displaystyle
C_{0}^{Tt}(m_{h}^{2})\approx-\frac{1}{m_{T}^{2}}\big{[}1+r_{tT}^{2}(1+\log
r_{tT}^{2})+\frac{\hat{s}}{12m_{T}^{2}}+\frac{m_{h}^{2}}{2m_{T}^{2}}\big{]}.$
(52)

###  A.3 Heavy quark expansion of $D_{0}$ function

$D_{0}$ function is defined as:

$\displaystyle
D_{0}(k_{1}^{2},k_{12}^{2},k_{23}^{2},k_{3}^{2},k_{2}^{2},k_{13}^{2},m_{0}^{2},m_{1}^{2},m_{2}^{2},m_{3}^{2})$
$\displaystyle\equiv$ $\displaystyle\frac{(2\pi\mu)^{4-D}}{i\pi^{2}}\int
d^{D}q\frac{1}{(q^{2}-m_{0}^{2})[(q+k_{1})^{2}-m_{1}^{2}][(q+k_{2})^{2}-m_{2}^{2}][(q+k_{3})^{2}-m_{3}^{2}]}$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydzdw\frac{\delta(x+y+z+w-1)}{[xm_{0}^{2}+ym_{1}^{2}+zm_{2}^{2}+wm_{3}^{2}-xyk_{1}^{2}-xzk_{2}^{2}-xwk_{3}^{2}-yzk_{12}^{2}-ywk_{13}^{2}-zwk_{23}^{2}]^{2}},$
(53)

where we have $k_{12}\equiv k_{1}-k_{2}$, $k_{23}\equiv k_{2}-k_{3}$, and
$k_{13}\equiv k_{1}-k_{3}$. When the four internal masses are all equal,
$D_{0}$ function can be expanded as [86]

$\displaystyle
D_{0}(k_{1}^{2},k_{12}^{2},k_{23}^{2},k_{3}^{2},k_{2}^{2},k_{13}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2})$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydzdw\frac{\delta(x+y+z+w-1)}{[m_{t}^{2}-xyk_{1}^{2}-xzk_{2}^{2}-xwk_{3}^{2}-yzk_{12}^{2}-ywk_{13}^{2}-zwk_{23}^{2}]^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{6m_{t}^{4}}\Big{[}1+\frac{k_{1}^{2}+k_{12}^{2}+k_{23}^{2}+k_{3}^{2}+k_{2}^{2}+k_{13}^{2}}{10m_{t}^{2}}+\frac{1}{140m_{t}^{4}}\Big{(}2(k_{1}^{4}+k_{2}^{4}+k_{3}^{4}+k_{1}^{2}k_{2}^{2}+k_{1}^{2}k_{3}^{2}+k_{2}^{2}k_{3}^{2})$
$\displaystyle+2(k_{12}^{4}+k_{13}^{4}+k_{23}^{4}+k_{12}^{2}k_{13}^{2}+k_{12}^{2}k_{23}^{2}+k_{13}^{2}k_{23}^{2})+2k_{1}^{2}(k_{12}^{2}+k_{13}^{2})+2k_{2}^{2}(k_{12}^{2}+k_{23}^{2})+2k_{3}^{2}(k_{13}^{2}+k_{23}^{2})$
$\displaystyle+(k_{1}^{2}k_{23}^{2}+k_{2}^{2}k_{13}^{2}+k_{3}^{2}k_{12}^{2})\Big{)}+\mathcal{O}(\frac{k^{6}}{m_{t}^{6}})\Big{]}.$
(54)

Especially, we have the following results

$\displaystyle
D_{0}^{t}(\hat{t},\hat{s})\approx\frac{1}{6m_{t}^{4}}[1+\frac{2m_{h}^{2}+\hat{s}+\hat{t}}{10m_{t}^{2}}+\frac{6m_{h}^{4}+4m_{h}^{2}(\hat{s}+\hat{t})+2\hat{s}^{2}+2\hat{t}^{2}+\hat{s}\hat{t}}{140m_{t}^{4}}],$
$\displaystyle
D_{0}^{t}(\hat{t},\hat{u})\approx\frac{1}{6m_{t}^{4}}[1+\frac{2m_{h}^{2}+\hat{t}+\hat{u}}{10m_{t}^{2}}+\frac{5m_{h}^{4}+4m_{h}^{2}(\hat{t}+\hat{u})+2\hat{t}^{2}+2\hat{u}^{2}+\hat{t}\hat{u}}{140m_{t}^{4}}].$
(55)

The expansion of $D_{0}^{t}(\hat{u},\hat{s})$ can be obtained when replacing
the $\hat{t}$ in $D_{0}^{t}(\hat{t},\hat{s})$ by $\hat{u}$. The expansion of
$D_{0}^{T}$ functions can be obtained when replacing the $m_{t}$ in
$D_{0}^{t}$ by $m_{T}$.

When three internal masses are equal, $D_{0}$ function can be expanded as

$\displaystyle
D_{0}(k_{1}^{2},k_{12}^{2},k_{23}^{2},k_{3}^{2},k_{2}^{2},k_{13}^{2},m_{T}^{2},m_{t}^{2},m_{t}^{2},m_{t}^{2})$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydzdw\frac{\delta(x+y+z+w-1)}{[xm_{T}^{2}+(y+z+w)m_{t}^{2}-xyk_{1}^{2}-xzk_{2}^{2}-xwk_{3}^{2}-yzk_{12}^{2}-ywk_{13}^{2}-zwk_{23}^{2}]^{2}}$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\\!\int_{0}^{1}\\!\int_{0}^{1}\\!\int_{0}^{1}\\!dxdydzdw\frac{\delta(x+y+z+w-1)}{[xm_{T}^{2}+(y+z+w)m_{t}^{2}]^{2}}[1\\!+\\!\frac{2(xyk_{1}^{2}+xzk_{2}^{2}+xwk_{3}^{2}+yzk_{12}^{2}+ywk_{13}^{2}+zwk_{23}^{2})}{xm_{T}^{2}+(y+z+w)m_{t}^{2}}]\\!+\\!\mathcal{O}(\frac{k^{4}}{m_{t,T}^{8}})$
$\displaystyle\approx$
$\displaystyle\frac{1}{m_{T}^{4}}\cdot\frac{1+2r_{tT}^{2}\log
r_{tT}^{2}-r_{tT}^{4}}{2r_{tT}^{2}(1-r_{tT}^{2})^{3}}+\frac{(k_{1}^{2}+k_{2}^{2}+k_{3}^{2})}{m_{T}^{6}}\cdot\frac{1+6r_{tT}^{2}(1+r_{tT}^{2})\log
r_{tT}^{2}+9r_{tT}^{2}-9r_{tT}^{4}-r_{tT}^{6}}{6r_{tT}^{2}(1-r_{tT}^{2})^{5}}$
$\displaystyle+\frac{(k_{12}^{2}+k_{23}^{2}+k_{13}^{2})}{m_{T}^{6}}\cdot\frac{1-8r_{tT}^{2}-12r_{tT}^{4}\log
r_{tT}^{2}+8r_{tT}^{6}-r_{tT}^{8}}{24r_{tT}^{4}(1-r_{tT}^{2})^{5}}.$ (56)

We only expand it up to $\mathcal{O}(\frac{k^{2}}{m_{t,T}^{6}})$, because the
general results will be quite lengthy. For the integral
$D_{0}^{tT}(\hat{t},\hat{s})$, we obtain the expression up to
$\mathcal{O}(\frac{k^{4}}{m_{t,T}^{8}})$

$\displaystyle
D_{0}^{tT}(\hat{t},\hat{s})=\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydzdw\frac{\delta(x+y+z+w-1)}{[xm_{T}^{2}+(y+z+w)m_{t}^{2}-x(y+w)m_{h}^{2}-xzt-
yws]^{2}}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydzdw\frac{\delta(x+y+z+w-1)}{[xm_{T}^{2}+(y+z+w)m_{t}^{2}]^{2}}\cdot$
$\displaystyle\left(1+\frac{2[x(y+w)m_{h}^{2}+xzt+yws]}{xm_{T}^{2}+(y+z+w)m_{t}^{2}}+\frac{3[x(y+w)m_{h}^{2}+xzt+yws]^{2}}{[xm_{T}^{2}+(y+z+w)m_{t}^{2}]^{2}}\right)+\mathcal{O}(\frac{1}{m_{t,T}^{10}})$
$\displaystyle\approx$
$\displaystyle\frac{1}{m_{T}^{4}}\cdot\frac{1+2r_{tT}^{2}\log
r_{tT}^{2}-r_{tT}^{4}}{2r_{tT}^{2}(1-r_{tT}^{2})^{3}}+\frac{(2m_{h}^{2}+\hat{t})}{m_{T}^{6}}\cdot\frac{1+6r_{tT}^{2}(1+r_{tT}^{2})\log
r_{tT}^{2}+9r_{tT}^{2}-9r_{tT}^{4}-r_{tT}^{6}}{6r_{tT}^{2}(1-r_{tT}^{2})^{5}}$
$\displaystyle+\frac{\hat{s}}{m_{T}^{6}}\cdot\frac{1-8r_{tT}^{2}-12r_{tT}^{4}\log
r_{tT}^{2}+8r_{tT}^{6}-r_{tT}^{8}}{24r_{tT}^{4}(1-r_{tT}^{2})^{5}}+\frac{\hat{s}^{2}}{m_{T}^{8}}\cdot\frac{1-9r_{tT}^{2}+45r_{tT}^{4}-45r_{tT}^{8}+9r_{tT}^{10}-r_{tT}^{12}+60r_{tT}^{6}\log
r_{tT}^{2}}{180r_{tT}^{6}(1-r_{tT}^{2})^{7}}$
$\displaystyle+\frac{\hat{s}(\hat{t}+4m_{h}^{2})}{m_{T}^{8}}\cdot\frac{1-15r_{tT}^{2}-80r_{tT}^{4}+80r_{tT}^{6}+15r_{tT}^{8}-r_{tT}^{10}-60r_{tT}^{4}(1+r_{tT}^{2})\log
r_{tT}^{2}}{120r_{tT}^{4}(1-r_{tT}^{2})^{7}}$
$\displaystyle+\frac{\hat{t}^{2}+2m_{h}^{2}\hat{t}+3m_{h}^{4}}{m_{T}^{8}}\cdot\frac{1+28r_{tT}^{2}-28r_{tT}^{6}-r_{tT}^{8}+12r_{tT}^{2}(1+3r_{tT}^{2}+r_{tT}^{4})\log
r_{tT}^{2}}{12r_{tT}^{2}(1-r_{tT}^{2})^{7}}.$ (57)

Keeping the terms up to $\mathcal{O}(\frac{1}{m_{T}^{6}})$ and considering the
$\log r_{tT}^{2}$ enhanced terms, they can be simplified as

$\displaystyle
D_{0}^{tT}(\hat{t},\hat{s})\approx\frac{1}{2m_{T}^{4}}\big{[}\frac{1}{r_{tT}^{2}}+3+2\log
r_{tT}^{2}+r_{tT}^{2}(5+6\log
r_{tT}^{2})\big{]}+\frac{(2m_{h}^{2}+\hat{t})}{6m_{T}^{6}}(\frac{1}{r_{tT}^{2}}+6\log
r_{tT}^{2}+14)$
$\displaystyle+\frac{\hat{s}}{24m_{T}^{6}}(\frac{1}{r_{tT}^{4}}-\frac{3}{r_{tT}^{2}}-12\log
r_{tT}^{2}-25).$ (58)

Similarly, we can get the following results

$\displaystyle
D_{0}^{Tt}(\hat{t},\hat{s})\approx\frac{1}{m_{T}^{4}}\cdot\frac{1+2r_{tT}^{2}\log
r_{tT}^{2}-r_{tT}^{4}}{2(1-r_{tT}^{2})^{3}}+\frac{(2m_{h}^{2}+\hat{t})}{m_{T}^{6}}\cdot\frac{1+9r_{tT}^{2}+6r_{tT}^{2}(1+r_{tT}^{2})\log
r_{tT}^{2}-9r_{tT}^{4}-r_{tT}^{6}}{6(1-r_{tT}^{2})^{5}}$
$\displaystyle+\frac{\hat{s}}{m_{T}^{6}}\cdot\frac{1-8r_{tT}^{2}-12r_{tT}^{4}\log
r_{tT}^{2}+8r_{tT}^{6}-r_{tT}^{8}}{24(1-r_{tT}^{2})^{5}}+\frac{\hat{s}^{2}}{m_{T}^{8}}\cdot\frac{1-9r_{tT}^{2}+45r_{tT}^{4}-45r_{tT}^{8}+9r_{tT}^{10}-r_{tT}^{12}+60r_{tT}^{6}\log
r_{tT}^{2}}{180(1-r_{tT}^{2})^{7}}$
$\displaystyle+\frac{\hat{s}(\hat{t}+4m_{h}^{2})}{m_{T}^{8}}\cdot\frac{1-15r_{tT}^{2}-80r_{tT}^{4}+80r_{tT}^{6}+15r_{tT}^{8}-r_{tT}^{10}-60r_{tT}^{4}(1+r_{tT}^{2})\log
r_{tT}^{2}}{120(1-r_{tT}^{2})^{7}}$
$\displaystyle+\frac{\hat{t}^{2}+2m_{h}^{2}\hat{t}+3m_{h}^{4}}{m_{T}^{8}}\cdot\frac{1+28r_{tT}^{2}-28r_{tT}^{6}-r_{tT}^{8}+12r_{tT}^{2}(1+3r_{tT}^{2}+r_{tT}^{4})\log
r_{tT}^{2}}{12(1-r_{tT}^{2})^{7}}.$ (59)

Keeping the terms up to $\mathcal{O}(\frac{1}{m_{T}^{6}})$ and considering the
$\log r_{tT}^{2}$ enhanced terms, they can be simplified as

$\displaystyle
D_{0}^{Tt}(\hat{t},\hat{s})\approx\frac{1}{2m_{T}^{4}}(1+3r_{tT}^{2}+2r_{tT}^{2}\log
r_{tT}^{2})+\frac{2m_{h}^{2}+\hat{t}}{6m_{T}^{6}}+\frac{\hat{s}}{24m_{T}^{6}}.$
(60)

When two internal masses are equal individually, we have the following
relations

$\displaystyle
D_{0}(k_{1}^{2},k_{12}^{2},k_{23}^{2},k_{3}^{2},k_{2}^{2},k_{13}^{2},m_{t}^{2},m_{T}^{2},m_{T}^{2},m_{t}^{2})=D_{0}(k_{1}^{2},k_{3}^{2},k_{23}^{2},k_{12}^{2},k_{13}^{2},k_{2}^{2},m_{T}^{2},m_{t}^{2},m_{t}^{2},m_{T}^{2})$
$\displaystyle=D_{0}(k_{23}^{2},k_{3}^{2},k_{1}^{2},k_{12}^{2},k_{2}^{2},k_{13}^{2},m_{T}^{2},m_{t}^{2},m_{t}^{2},m_{T}^{2}).$
(61)

$D_{0}$ function can be expanded as

$\displaystyle
D_{0}(k_{1}^{2},k_{12}^{2},k_{23}^{2},k_{3}^{2},k_{2}^{2},k_{13}^{2},m_{t}^{2},m_{T}^{2},m_{T}^{2},m_{t}^{2})$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydzdw\frac{\delta(x+y+z+w-1)}{[(y+z)m_{T}^{2}+(x+w)m_{t}^{2}-xyk_{1}^{2}-xzk_{2}^{2}-xwk_{3}^{2}-yzk_{12}^{2}-ywk_{13}^{2}-zwk_{23}^{2}]^{2}}$
$\displaystyle=$
$\displaystyle\int_{0}^{1}\\!\int_{0}^{1}\\!\int_{0}^{1}\\!\int_{0}^{1}\\!dxdydzdw\frac{\delta(x+y+z+w-1)}{[(y+z)m_{T}^{2}+(x+w)m_{t}^{2}]^{2}}[1\\!+\\!\frac{2(xyk_{1}^{2}+xzk_{2}^{2}+xwk_{3}^{2}+yzk_{12}^{2}+ywk_{13}^{2}+zwk_{23}^{2})}{(y+z)m_{T}^{2}+(x+w)m_{t}^{2}}]\\!+\\!\mathcal{O}(\frac{k^{4}}{m_{t,T}^{8}})$
$\displaystyle\approx$
$\displaystyle-\frac{1}{m_{T}^{4}}\cdot\frac{2+(1+r_{tT}^{2})\log
r_{tT}^{2}-2r_{tT}^{2}}{(1-r_{tT}^{2})^{3}}-\frac{k_{1}^{2}+k_{2}^{2}+k_{13}^{2}+k_{23}^{2}}{m_{T}^{6}}\cdot\frac{3+(1+4r_{tT}^{2}+r_{tT}^{4})\log
r_{tT}^{2}-3r_{tT}^{4}}{2(1-r_{tT}^{2})^{5}}$
$\displaystyle+\frac{k_{3}^{2}}{m_{T}^{6}}\cdot\frac{1+6r_{tT}^{2}(1+r_{tT}^{2})\log
r_{tT}^{2}+9r_{tT}^{2}-9r_{tT}^{4}-r_{tT}^{6}}{6r_{tT}^{2}(1-r_{tT}^{2})^{5}}+\frac{k_{12}^{2}}{m_{T}^{6}}\cdot\frac{1+6r_{tT}^{2}(1+r_{tT}^{2})\log
r_{tT}^{2}+9r_{tT}^{2}-9r_{tT}^{4}-r_{tT}^{6}}{6(1-r_{tT}^{2})^{5}}.$ (62)

We only expand it up to $\mathcal{O}(\frac{k^{2}}{m_{t,T}^{6}})$, because the
general results will be quite lengthy. For the integral
$D_{0}^{tT}(\hat{t},\hat{u})$, we obtain the expression up to
$\mathcal{O}(\frac{k^{4}}{m_{t,T}^{8}})$

$\displaystyle
D_{0}^{tT}(\hat{t},\hat{u})=\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydzdw\frac{\delta(x+y+z+w-1)}{[(x+w)m_{t}^{2}+(y+z)m_{T}^{2}-(xy+zw)m_{h}^{2}-xzt-
ywu]^{2}}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}dxdydzdw\frac{\delta(x+y+z+w-1)}{[(x+w)m_{t}^{2}+(y+z)m_{T}^{2}]^{2}}\cdot$
$\displaystyle\left(1+\frac{2[(xy+zw)m_{h}^{2}+xzt+ywu]}{(x+w)m_{t}^{2}+(y+z)m_{T}^{2}}+\frac{3[(xy+zw)m_{h}^{2}+xzt+ywu]^{2}}{[(x+w)m_{t}^{2}+(y+z)m_{T}^{2}]^{2}}\right)+\mathcal{O}(\frac{1}{m_{t,T}^{10}})$
$\displaystyle\approx$
$\displaystyle-\frac{1}{m_{T}^{4}}\cdot\frac{2+(1+r_{tT}^{2})\log
r_{tT}^{2}-2r_{tT}^{2}}{(1-r_{tT}^{2})^{3}}-\frac{2m_{h}^{2}+\hat{t}+\hat{u}}{m_{T}^{6}}\cdot\frac{3+(1+4r_{tT}^{2}+r_{tT}^{4})\log
r_{tT}^{2}-3r_{tT}^{4}}{2(1-r_{tT}^{2})^{5}}$ $\displaystyle-$
$\displaystyle\frac{5m_{h}^{4}+4m_{h}^{2}(\hat{t}+\hat{u})+2\hat{t}^{2}+\hat{t}\hat{u}+2\hat{u}^{2}}{m_{T}^{8}}\cdot\frac{11+27r_{tT}^{2}-27r_{tT}^{4}-11r_{tT}^{6}+3(1+9r_{tT}^{2}+9r_{tT}^{4}+r_{tT}^{6})\log
r_{tT}^{2}}{18(1-r_{tT}^{2})^{7}}.$ (63)

Keeping the terms up to $\mathcal{O}(\frac{1}{m_{T}^{4}})$ and considering the
$\log r_{tT}^{2}$ enhanced terms, they can be simplified as

$\displaystyle
D_{0}^{tT}(\hat{t},\hat{u})\approx-\frac{1}{m_{T}^{4}}\big{[}2+\log
r_{tT}^{2}+4r_{tT}^{2}(1+\log
r_{tT}^{2})\big{]}-\frac{4m_{h}^{2}-\hat{s}}{2m_{T}^{6}}(3+\log r_{tT}^{2}).$
(64)

In the above calculations, the $t$ and $T$ quark mixed $C_{0}$ and $D_{0}$
integrals will agree with the pure top quark integrals in the limit of
$m_{t}=m_{T}$ (or $r_{tT}\rightarrow 1$). Besides, these expansion results
have been checked by the LoopTools numerically [106].