Signature of 2HDM at Higgs Factories

TL;DR

This paper calculates one-loop electroweak and QED corrections to Higgs production in the 2HDM at future electron-positron colliders, highlighting potential measurable deviations from the Standard Model predictions.

Contribution

It provides the first detailed one-loop correction analysis for $e^+ e^- o Z o ext{Higgs}$ processes in the 2HDM, including benchmark scenarios for future collider tests.

Findings

Corrections are sizeable and potentially measurable at future colliders.

Deviations from the Standard Model can be significant in certain benchmark scenarios.

Electroweak and real emission effects are crucial for precise predictions.

Abstract

The full one-loop corrections, both the weak and QED corrections, to the process $e^+ e^- \to Z \phi $ ($\phi=h^0,H^0$) in the two Higgs doublet model (2HDM) at the Higgs factories are presented. Up to the $O(\alpha_{em})$ level, the virtual corrections are evaluated by using the FeynArts/FormCalc packages. The real emission corrections are computed using the Feynman Diagram Calculation (FDC) package and the collinear divergences are regularized by the structure functions of an electron. Using the FeynArts/FormCalc and the FDC packages, we study the corrections in the Standard Model (SM) and the 2HDM, respectively. Gauge dependence arising in the normalization of mixing angles is removed by using the pinch technique. After taking into account experimental constraints from the current LHC data, we propose four interesting benchmark scenarios for future colliders. By using these benchmark scenarios, we evaluate the deviation of $\Delta\sigma(e^+ e^- \to Z \phi)$ from their SM values. We also examine Higgs boson decays $\phi\to b\bar{b}$ and $\phi\to \tau^+\tau^-$, which can have large electroweak (EW) contributions from triple Higgs couplings which are absent in the SM. It is found that for these benchmark scenarios, both EW and real emission corrections are sizeable and could be measured at a future $e^+ e^-$ colliders such as the ILC, CLIC, and CEPC.

Tables

Table 1. Table 1: Experimental results of the observables combined by the Particle Data Group (PDG) and/or Heavy Flavor Averaging Group (HFAG) Collaborations in Refs.  Olive  et al. ( 2014 ) – Amhis  et al. ( 2014 ) . As for ℬ ¯ ​ ( B q 0 → μ + ​ μ − ) ¯ ℬ → subscript superscript 𝐵 0 𝑞 superscript 𝜇 superscript 𝜇 \overline{\mathcal{B}}(B^{0}_{q}\to\mu^{+}\mu^{-}) , the combined results from the LHCb and CMS collaborations are shown as given in Ref.  Archilli ( 2014 ) . Ref.  Hagiwara  et al. ( 2011 ) is used for constraints from ( g − 2 ) μ subscript 𝑔 2 𝜇 (g-2)_{\mu} data.

Observable	Experimental result	SM contribution	Combined error at 1$\sigma$
$ℬ\left(K\to \mu \nu \right)/ℬ\left(\pi \to \mu \nu \right)$	$0.6357\pm 0.0011$ Olive et al. (2014)	$0.6231\pm 0.0071$	0.0071
$\overline{ℬ}{\left(b\to s\gamma \right)}_{E_{\gamma }>1.6\text{GeV}}$	$\left(3.32\pm 0.16\right)\times {10}^{-4}$ Amhis et al. (2017)	$\left(3.36\pm 0.24\right)\times {10}^{-4}$	$0.29\times {10}^{-4}$
$ℬ\left(B\to \tau \nu \right)$	$\left(1.14\pm 0.22\right)\times {10}^{-4}$ Amhis et al. (2014)	$\left(0.78\pm 0.07\right)\times {10}^{-4}$	$0.23\times {10}^{-4}$
$ℬ\left(D\to \mu \nu \right)$	$\left(3.74\pm 0.17\right)\times {10}^{-4}$ Olive et al. (2014); Amhis et al. (2014)	$\left(3.94\pm 0.13\right)\times {10}^{-4}$	$0.21\times {10}^{-4}$
$ℬ\left(D_s\to \tau \nu \right)$	$\left(5.55\pm 0.24\right)\times {10}^{-2}$ Olive et al. (2014); Amhis et al. (2014)	$\left(5.17\pm 0.11\right)\times {10}^{-2}$	$0.26\times {10}^{-2}$
$ℬ\left(D_s\to \mu \nu \right)$	$\left(5.57\pm 0.24\right)\times {10}^{-3}$ Olive et al. (2014); Amhis et al. (2014)	$\left(5.28\pm 0.11\right)\times {10}^{-3}$	$0.26\times {10}^{-3}$
$\overline{ℬ}\left(B_s^0\to {\mu }^{+}{\mu }^{-}\right)$	$\left(2.8\pm 0.7\right)\times {10}^{-9}$ Archilli (2014)	$\left(3.66\pm 0.28\right)\times {10}^{-9}$	$0.75\times {10}^{-9}$
$\overline{ℬ}\left(B_d^0\to {\mu }^{+}{\mu }^{-}\right)$	$\left(3.9\pm 1.5\right)\times {10}^{-10}$ Archilli (2014)	$\left(1.08\pm 0.13\right)\times {10}^{-10}$	$1.50\times {10}^{-10}$
$\Delta M_s$	$\left(17.757\pm 0.021\right){\text{ ps}}^{-1}$ Amhis et al. (2017)	$\left(18.257\pm 1.505\right){\text{ ps}}^{-1}$	$1.5{\text{ ps}}^{-1}$
$\Delta M_d$	$\left(0.510\pm 0.002\right){\text{ ps}}^{-1}$ Amhis et al. (2017)	$\left(0.548\pm 0.075\right){\text{ ps}}^{-1}$	$0.075{\text{ ps}}^{-1}$
${\Delta }_0\left(B\to K^{\ast }\gamma \right)$	$\left(5.2\pm 2.6\right)\times {10}^{-2}$ Olive et al. (2014)	$\left(5.1\pm 1.5\right)\times {10}^{-2}$	$3.0\times {10}^{-2}$
$\delta a_{\mu }$	$\left(261\pm 80\right)\times {10}^{-11}$ Hagiwara et al. (2011)	$-$	$80\times {10}^{-11}$

Table 2. Table 2: Selected benchmark points using Higgs data at 13 TeV with m h subscript 𝑚 ℎ m_{h} = 125 GeV are presented, for BP1-h and BP2-h and with m H subscript 𝑚 𝐻 m_{H} = 125 GeV, for BP1-H and BP2-H. Obviously, in this notation, BP1 (BP2) type is related to type-I (type-II). Here we adopt four physical masses, two angles, and λ 5 subscript 𝜆 5 \lambda_{5} to describe each point in the parameter space. In each BP, we examine the ratio of total decay width over the mass of a Higgs boson, i.e. Γ ϕ / m ϕ subscript Γ italic-ϕ subscript 𝑚 italic-ϕ \Gamma_{\phi}/m_{\phi} , and find it is smaller than 5 % percent 5 5\% with ϕ = H , A italic-ϕ 𝐻 𝐴 \phi=H,A and H ± . superscript 𝐻 plus-or-minus H^{\pm}.

BPs	$\sin \left(\beta -\alpha \right)$	$\tan \beta$	$m_h$ (GeV)	$m_H$ (GeV)	$m_A$ (GeV)	$m_{H^{\pm }}$ (GeV)	${\lambda }_5$
BP1-h	0.99679	14.300	125.00	212.00	98.20	178.27	0.5819
BP2-h	0.99999	2.012	125.00	594.00	512.00	592.00	0.0000
BP1-H	-0.06000	2.830	95.00	125.00	169.00	170.00	-0.3220
BP2-H	-0.03000	2.160	95.00	125.00	600.00	600.00	-5.7800

Table 3. Table 3: NLO SM results in different schemes at s = 250 𝑠 250 \sqrt{s}=250 GeV are shown (in unit of fb). Three schemes are chosen to demonstrate the scale dependence of the results.

scheme	$1/{\alpha }_{em}\left(\mu \right)$	${\sigma }^0$	${\sigma }^{1,weak}$	${\sigma }^1$	${\sigma }^{NLO}={\sigma }^0+{\sigma }^1$
${\alpha }_{em}\left(0\right)$	137.036	223.12(0)	6.09(0)	7.13(2)	230.25(2)
${\alpha }_{em}\left(m_Z\right)$	128.943	252.00(0)	-24.33(0)	-23.07(2)	228.93(2)
${\alpha }_{em}\left(\sqrt{s}\right)$	127.515	257.68(0)	-30.92(0)	-29.63(2)	228.05(2)

Table 4. Table 4: Check for Δ ​ E Δ 𝐸 \Delta E independence at s = 250 𝑠 250 \sqrt{s}=250 GeV (in unit of 10 − 1 superscript 10 1 10^{-1} pb).

${\delta }_s=2\Delta E/\sqrt{s}$	${\sigma }_{S+V}$	${\sigma }_{H\overline{C}}$	${\sigma }_{HC+CT}^{\ast }$	${\sigma }_{SC}$	${\sigma }^1$
${10}^{-1}$	-0.7127(0)	0.1240(0)	-0.1209(0)	0.4794(0)	-0.2302(0)
${10}^{-2}$	-1.4347(0)	0.5445(0)	-0.5306(1)	1.1903(0)	-0.2305(1)
${10}^{-3}$	-2.1567(0)	0.9788(1)	-0.9540(2)	1.9012(0)	-0.2307(2)
${10}^{-4}$	-2.8787(0)	1.4142(1)	-1.3784(2)	2.6121(0)	-0.2308(2)
${10}^{-5}$	-3.6006(0)	1.8497(2)	-1.8027(4)	3.3230(0)	-0.2306(4)
${10}^{-6}$	-4.3227(0)	2.2853(2)	-2.2271(5)	4.0339(0)	-0.2306(5)
${10}^{-7}$	-5.0446(0)	2.7208(2)	-2.6516(6)	4.7448(0)	-0.2306(6)
${10}^{-8}$	-5.7666(0)	3.1564(3)	-3.0762(7)	5.4558(0)	-0.2306(8)

Table 5. Table 5: Check for Δ ​ θ Δ 𝜃 \Delta\theta independence at s = 250 𝑠 250 \sqrt{s}=250 GeV (in unit of 10 − 1 superscript 10 1 10^{-1} pb).

$\Delta \theta$	${\sigma }_{V+S}$	${\sigma }_{H\overline{C}}$	${\sigma }_{HC+CT}^{\ast }$	${\sigma }_{SC}$	${\sigma }^1$
${10}^{-1}$	-2.1567(0)	0.3856(1)	-0.3608(1)	1.9012(0)	-0.2307(1)
${10}^{-2}$	-2.1567(0)	0.6822(1)	-0.6574(1)	1.9012(0)	-0.2307(1)
${10}^{-3}$	-2.1567(0)	0.9788(1)	-0.9540(2)	1.9012(0)	-0.2307(2)
${10}^{-4}$	-2.1567(0)	1.2751(1)	-1.2506(3)	1.9012(0)	-0.2310(3)
${10}^{-5}$	-2.1567(0)	1.5527(1)	-1.5471(3)	1.9012(0)	-0.2499(3)
${10}^{-6}$	-2.1567(0)	1.6227(2)	-1.8437(4)	1.9012(0)	-0.4765(4)

Table 6. Table 6: Check for m e subscript 𝑚 𝑒 m_{e} independence at s = 250 𝑠 250 \sqrt{s}=250 GeV (in unit of 10 − 1 superscript 10 1 10^{-1} pb).

$k$	${\sigma }_{V+S}$	${\sigma }_{H\overline{C}}$	${\sigma }_{HC+CT}^{\ast }$	${\sigma }_{SC}$	${\sigma }^1$
$2^{-4}$	-2.5815(0)	0.9788(1)	-0.9540(2)	2.3260(0)	-0.2307(2)
$2^{-3}$	-2.4753(0)	0.9788(1)	-0.9540(2)	2.2198(0)	-0.2307(2)
$2^{-2}$	-2.3691(0)	0.9788(1)	-0.9540(2)	2.1136(0)	-0.2307(2)
$2^{-1}$	-2.2629(0)	0.9788(1)	-0.9540(2)	2.0074(0)	-0.2307(2)
$2^0$	-2.1567(0)	0.9788(1)	-0.9540(2)	1.9012(0)	-0.2307(2)
$2^1$	-2.0505(0)	0.9788(1)	-0.9540(2)	1.7950(0)	-0.2307(2)
$2^2$	-1.9443(0)	0.9788(1)	-0.9540(2)	1.6888(0)	-0.2307(2)
$2^3$	-1.8381(0)	0.9788(1)	-0.9540(2)	1.5826(0)	-0.2307(2)
$2^4$	-1.7318(0)	0.9788(1)	-0.9540(2)	1.4763(0)	-0.2307(2)
$2^5$	-1.6256(0)	0.9788(1)	-0.9540(2)	1.3701(0)	-0.2307(2)
$2^6$	-1.5194(0)	0.9788(1)	-0.9540(2)	1.2639(0)	-0.2307(2)
$2^7$	-1.4132(0)	0.9788(1)	-0.9540(2)	1.1577(0)	-0.2307(2)
$2^8$	-1.3070(0)	0.9788(1)	-0.9540(2)	1.0515(0)	-0.2307(2)

Table 7. Table 7: Dimensionless couplings λ ^ i subscript ^ 𝜆 𝑖 \hat{\lambda}_{i} with two typical values of v 𝑣 v , the inputs of each of benchmark points are given in Table 2 . λ ^ i subscript ^ 𝜆 𝑖 \hat{\lambda}_{i} are defined from the THCs via Eq. ( 38 ). The values of v 𝑣 v are taken as explained in the second paragraph of Section IV .

$v$(GeV)	BPs	${\hat{\lambda }}_{h^0{h}^0{h}^0}$	${\hat{\lambda }}_{H^0{h}^0{h}^0}$	${\hat{\lambda }}_{h^0{H}^0{H}^0}$	${\hat{\lambda }}_{h^0{H}^{\pm }{H}^{\pm }}$	${\hat{\lambda }}_{h^0{A}^0{A}^0}$	${\hat{\lambda }}_{H^0{H}^{\pm }{H}^{\pm }}$	${\hat{\lambda }}_{H^0{A}^0{A}^0}$	${\hat{\lambda }}_{H^0{H}^0{H}^0}$	$m_{12}^2$(GeV2)
246.220	BP1-h	-0.763	-0.048	-0.810	-0.375	0.353	-0.030	0.028	-0.195	3126.00
	BP2-h	-0.773	-0.049	-3.257	-3.199	-0.285	-2.305	-2.292	-6.917	104480.72
	BP1-H	0.025	-0.240	-0.002	0.048	0.047	-0.895	-0.884	-0.774	2839.78
	BP2-H	-0.036	-0.235	-0.001	0.335	0.335	-11.807	-11.807	-0.773	3567.16
243.137	BP1-h	-0.763	-0.047	-0.789	-0.388	0.340	-0.238	-0.179	-0.814	3064.92
	BP2-h	-0.773	-0.049	-3.257	-3.199	-0.285	-2.305	-2.292	-6.917	104480.72
	BP1-H	-0.037	-0.228	-0.0003	0.027	0.026	-0.880	-0.869	-0.774	2992.36
	BP2-H	-0.793	0.030	0.016	0.083	0.083	-11.527	-11.527	-0.773	6981.21

Table 8. Table 8: Allowed ranges of λ 5 subscript 𝜆 5 \lambda_{5} , λ 5 ​ v 2 subscript 𝜆 5 superscript 𝑣 2 \lambda_{5}v^{2} and m 12 2 superscript subscript 𝑚 12 2 m_{12}^{2} with two typical values of v 𝑣 v . The constraints from vacuum stability, unitarity and perturbativity are used to determine the ranges. The values of v 𝑣 v are taken as explained in the second paragraph of Section IV .

$v$(GeV)	BPs	BP1-h	BP2-h	BP1-H	BP2-H
246.220	${\lambda }_5$	[0.545,0.583]	[-0.482,0.377]	[-1.330,-0.294]	[-6.500,-5.737]
	${\lambda }_5{v}^2$(GeV2)	[33040, 35368]	[-29197, 22874]	[-80612, -17848]	[-394063, -347807]
	$m_{12}^2$(GeV2)	[2970, 3132]	[92844, 113597]	[-16351, 3365]	[-12986, 4649]
243.137	${\lambda }_5$	[0.560, 0.598]	[-0.441, 0.387]	[-1.337, -0.302]	[-6.601, -5.884]
	${\lambda }_5{v}^2$(GeV2)	[33099, 35367]	[-26094, 22872]	[-79032, -17847]	[-390217, -347807]
	$m_{12}^2$(GeV2)	[2974, 3132]	[94081, 113597]	[-15855, 3366]	[-11520, 4649]

Equations

$$V({\Phi }_1,{\Phi }_2)$$

$${\sigma }^0\left(e^{+}{e}^{-}\to Z\varphi \right)={\sin }^2(\beta -\alpha )[{\cos }^2(\beta -\alpha )]{\sigma }_{SM}^0\left(e^{+}{e}^{-}\to Zh\right)\mathrm{for}\varphi =h[H]$$

$${\sigma }_{SM}^0\left(e^{+}{e}^{-}\to Zh\right)=\frac{{\alpha }_{em}^2\pi }{192s{s}_W^4{c}_W^4}\left(v_e^2+a_e^2\right){\lambda }^{\frac{1}{2}}\frac{\lambda +12m_Z^2/s}{{\left(1-m_Z^2/s\right)}^2}$$

$$\delta Z_e=-\frac{1}{2}\delta Z_{AA}-\frac{s_W}{c_W}\frac{1}{2}\delta Z_{ZA}=\frac{1}{2}\Pi (0)-\frac{s_W}{c_W}\frac{{\sum }_T^{AZ}(0)}{m_Z^2}$$

$$\Pi (s)\equiv \frac{{\sum }_T^{AA}(s)}{s},$$

$$\Pi(0)=\lim_{s\rightarrow 0}\dfrac{\sum^{AA}_{T}(s)}{s}=\dfrac{\partial\sum^{AA}_{T}(s)}{\partial s}\biggr{|}_{s=0}.$$

$$\Pi (0)={\Pi }_{\mathrm{hadron}}^{(5)}(0)+{\Pi }_{\mathrm{lepton}}(0)+{\Pi }_{\mathrm{remaining}}(0)\to \mathrm{Re}{\Pi }_{\mathrm{hadron}}^{(5)}(m_Z^2)+\Delta {\alpha }_{\mathrm{hadron}}^{(5)}(m_Z)+{\Pi }_{\mathrm{lepton}}(0)+{\Pi }_{\mathrm{remaining}}(0),$$

$$\delta Z_e(0)=\frac{1}{2}\mathrm{Re}{\Pi }_{\mathrm{hadron}}^{(5)}(m_Z^2)+\frac{1}{2}\Delta {\alpha }_{\mathrm{hadron}}^{(5)}(m_Z)+\frac{1}{2}{\Pi }_{\mathrm{lepton}}(0)+\frac{1}{2}{\Pi }_{\mathrm{remaining}}(0)-\frac{s_W}{c_W}\frac{{\sum }_T^{AZ}(0)}{m_Z^2}.$$

$$\delta Z_e(\mu )\equiv \delta Z_e(0)-\frac{1}{2}\Delta \alpha (\mu )=\frac{1}{2}\mathrm{Re}{\Pi }_{\mathrm{hadron}}^{(5)}({\mu }^2)+\frac{1}{2}\mathrm{Re}{\Pi }_{\mathrm{lepton}}({\mu }^2)+\frac{1}{2}{\Pi }_{\mathrm{remaining}}(0)-\frac{s_W}{c_W}\frac{{\sum }_T^{AZ}(0)}{m_Z^2},$$

$$\Delta \alpha (\mu )$$

$${\alpha }_{em}(\mu )\equiv \frac{{\alpha }_{em}(0)}{1-\Delta \alpha (\mu )}.$$

$$d{\sigma }^1=d{\sigma }_V(\lambda )+d{\sigma }_S(\lambda ,\Delta E)+d{\sigma }_{HC+CT}(\Delta E,\Delta \theta )+d{\sigma }_{H\overline{C}}(\Delta E,\Delta \theta ),$$

$$d{\sigma }_S$$

$$d{\sigma }_{HC}=$$

$$\stackrel{\Delta {\theta }^2\gg m_e^2/s}{\to }$$

$$d\sigma_{CT}=-\dfrac{\alpha_{em}}{2\pi}\log\dfrac{s}{4m_{e}^{2}}P_{ee}^{+}(z,0)\times\biggl{[}d\sigma^{0}(zk_{1})+d\sigma^{0}(zk_{2})\biggr{]}dz$$

$$d{\sigma }_{HC+CT}$$

$$d{\sigma }_{HC+CT}^{\ast }$$

$$d{\sigma }_{SC}$$

$$\delta {\mathcal{L}}_{h^{0}ZZ}$$

$$\delta {\mathcal{L}}_{H^{0}ZZ}$$

$${\sigma }^{NLO}$$

$${\sigma }^1={\sigma }^{\text{mbox}1,weak}+{\sigma }^{\text{mbox}1,QED}.$$

$${\delta }^1={\delta }^{\text{mbox}1,weak}+{\delta }^{\text{mbox}1,QED},$$

$${\lambda }_5=\frac{m_{12}^2-m_A^2{s}_{\beta }{c}_{\beta }}{v^2{s}_{\beta }{c}_{\beta }}$$

$${\Delta }^{weak}=\frac{{\sigma }_{2HDM}^0+{\sigma }_{2HDM}^{1,weak}}{{\sigma }_{SM}^0+{\sigma }_{SM}^{1,weak}}-1,$$

$${\Gamma }_1(\varphi \to f\bar{f})=\frac{N_C{G}_F{m}_f^2}{4\sqrt{2}\pi }{\beta }_f^3({\xi }_f^{\varphi }{)}^2{m}_{\varphi }{\hat{Z}}_{\varphi }[1-\Delta r+2\mathrm{Re}(\Delta {\mathcal{M}}_1)],$$

$${\Delta }_{ff}(\varphi )$$

$${\sigma }^{LO}(s)={\bar{\sigma }}^{W,LO}(s)R(1-\frac{\Delta E}{E},s),$$

$${\sigma }^{NLO}(s)={\bar{\sigma }}^{W,NLO}(s)R(1-\frac{\Delta E}{E},s)$$

$$f_{ee}(x,s)=\delta (1-x)+\frac{{\alpha }_{em}}{2\pi }\log \frac{s}{4m_e^2}{P}_{ee}^{+}(x,0)$$

$$P_{ee}^{+}(z,0)=\frac{1+z^2}{(1-z{)}_{+}}+\frac{3}{2}\delta (1-z),$$

Full text

I Introduction

A Higgs-like boson ($h$) has been discovered in the first Run of the Large Hadron Collider (LHC) in 2012 Chatrchyan et al. (2012); Aad et al. (2012). Based on the data of Run-1, both ATLAS and CMS collaborations have established a Higgs boson mass with $m_{h}=125.09\pm 0.21(stat.)\pm 0.11(syst.)$ GeV Aad et al. (2015a). ATLAS and CMS also performed several Higgs coupling measurements to a precision around $10\%-20\%$, such as Higgs couplings to di-bosons $VV$ with $V=W^{\pm},Z,\gamma$ through a global fit with input from the measurement of Higgs decaying into two weak bosons, weak vector boson associated processes, and VBF processes Aad et al. (2016). Very recently, ATLAS and CMS have also measured Higgs couplings to the fermions of the third generation over 5$\sigma$, via the process $pp\to Vh(h\to b{\bar{b}})$ Aaboud et al. (2018a); Sirunyan et al. (2018a), via the VBF process $pp\to jjh(h\to\tau^{+}\tau^{-})$ Aaboud et al. (2019); Sirunyan et al. (2018b), and $pp\to tth$ with combined Higgs boson decay final states Sirunyan et al. (2018c); Aaboud et al. (2018b). With Run-2 data, high precision for these measured Higgs boson couplings has been achieved Sirunyan et al. (2019); Aad et al. (2020). These measurements demonstrated that the SM is consistent with the current Higgs data.

At the LHC, due to the large theoretical uncertainties like the Parton Distribution Function (PDF) and large experimental background, the precision for the measured Higgs boson couplings can reach $4\%-6\%$ or so at 300 fb*-1* Collaboration (2018) . In contrast, due to its clean environment, an $e^{+}e^{-}$ collider can perform precision measurements on the production and decay of the observed Higgs boson. For example, at 240 GeV $e^{+}e^{-}$ collider, the Higgs-strahlung process $e^{+}e^{-}\to Zh$ is the dominant Higgs production channel, where the Higgs boson can be reconstructed by using the recoiled $Z$ boson via its leptonic decay Garcia-Abia and Lohmann (2000). When a sufficiently large dataset is accumulated, say 1 ab*-1*, the projected precision for the vertex $ZZh$ can reach to $0.1\%$ or so Dong et al. (2018). For a center-of-mass energy of 240-250 GeV and an integrated luminosity of 250 fb*-1* for each of the two detectors in each year, a total ${\cal O}(10^{6})$ Higgs bosons can be produced in a 10 year running (a total integrated luminosity can be estimated as 5 ab*-1*), which can lead to a precision of those currently measured Higgs couplings around a percent level or less Asner et al. (2013); Peskin (2012); Blondel et al. (2013); Belusevic and Higo (2012); Dong et al. (2018). At the international linear collider (ILC) experiments Baer et al. (2013), the luminosity of Higgs-factory run can be around 250 fb*-1* at $\sqrt{s}=250$ GeV. The production rate of Higgs bosons from the Higgs-strahlung process can be enhanced when polarized beams Moortgat-Pick et al. (2008) are used. However, at $\sqrt{s}=500$ GeV and/or $\sqrt{s}=800$ GeV, the vector-boson fusion process $WW\to h$ can be dominant Moortgat-Pick et al. (2008) and ILC will have the capability to reach a precision around 1% for most of Higgs couplings Arbey et al. (2015). Correspondingly, precise theoretical predictions to the physical observables related to Higgs bosons, like cross sections and branching ratios, are required.

In this paper, we propose some benchmark scenarios of 2HDM after taking into account theoretical constraints as well as experimental restrictions from recent LHC data. By using these benchmark scenarios, we will study the effects of electroweak radiative corrections to the production process $e^{+}e^{-}\to Z\phi$ $\phi=h^{0},H^{0}$ and report full next-to-leading order (NLO) calculations on the cross-section in 2HDM by including both virtual corrections up to one-loop and real emissions of a photon. Our results are consistent with Lopez-Val et al. (2010) within the same set of 2HDM parameters. The complete SM one-loop correction of $\sigma(e^{+}e^{-}\to Zh_{SM})$ was first time calculated in Fleischer and Jegerlehner (1983) then followed by works Kniehl (1992a) and Denner et al. (1992), later on had been calculated using the GRACE system***GRACE is a system of program packages for automatic calculation, see https://minami-home.kek.jp in Belanger et al. (2003) and we reproduced the results for comparison. We find that the EW corrections of new physics in our benchmark scenarios can be of order $-10\%\sim-20\%$, while the contribution of real emission has a positive sign. Four benchmark points are studied in this paper and found that both EW and real emission corrections are sizeable, and could be the physical target at future $e^{+}e^{-}$ colliders such as the ILC, CLIC, and CEPC.

One-loop corrections of the Higgs boson couplings have been investigated in various models. QCD and electroweak corrections of $hf\bar{f}$ were calculated in the SM in Refs. Sakai (1980); Inami and Kubota (1981); Drees and Hikasa (1990); Fleischer and Jegerlehner (1981); Kniehl (1992b). Similar calculations had been done in various models beyond the SM like in the minimal supersymmetric SM (MSSM). Due to its sizeable SUSY-QCD effects, it was observed that $hf\bar{f}$ couplings can receive a large corrections Dabelstein (1995); Haber et al. (2001); Guasch et al. (2001). Moreover, $hf\bar{f}$ couplings had been intensively studied in the 2HDM Kanemura et al. (2014); Arhrib et al. (2004, 2017a). In this work, we also examine the deviations of $hf\bar{f}$ ($f=b,\tau$) in both type-I and type-II scenarios and find these deviations are within the reach of future Higgs factories and are helpful to distinguish new physics models.

The plan of the paper is as follows. In section 2, we introduce the benchmark points of the 2HDM model. In section 3, we outline the framework of our calculation and specifies the renormalization scheme we will use. In section 4, we present the numerical results. We conclude this work in section 5.

II The 2HDM model and benchmark points

2HDM was first introduced by T.D. Lee Lee (1973) in order to have spontaneous T parity (or equally CP parity) breaking, and a more comprehensive review can be found in Refs. Branco et al. (2012) and Davidson and Haber (2005). Before proposing benchmark scenarios, we briefly review the main features of 2HDM model related to this work††† In this work, we follow the notations of the reference Arhrib et al. (2017a) where more details can be found.. The most general renormalizable potential which is invariant under $SU(2)\times U(1)$ can be written as

[TABLE]

where $\Phi_{i}$, $i=1,2$ are complex $SU(2)$ doublets with 4 degrees of freedom each and all $m^{2}_{1}$, $m^{2}_{2}$ and $\lambda_{1-4}$ are real which follows from the hermiticity of the potential. However, the parameters $\lambda_{5-7}$ and $m^{2}_{12}$ can be complex in general. Explicit CP-violation may arise in the Higgs sector by considering the imaginary parts of the above complex parameters. In this work we focus on our study in the CP-conserving case so that we assume all parameters to be real. Furthermore, for the sake of simplicity, we only include the soft breaking term proportional to $m_{12}^{2}$ and omit the terms proportional to $\lambda_{6}$ and $\lambda_{7}$ as these lead to hard $Z_{2}$ violation.

The $Z_{2}$ symmetry is defined as ($\Phi_{1}$, $\Phi_{2})\to(\Phi_{1},-\Phi_{2})$, which was introduced in order to suppress the tree-level flavor changing neutral current (FCNC) processes Glashow and Weinberg (1977); Paschos (1977). The exact $Z_{2}$ symmetry will lead to the absence of $\Phi_{2}$ couplings into fermions. which makes $\Phi_{2}$ possible as a dark matter candidate in the so-called inert doublet model Ma (2006); Barbieri et al. (2006); Cao et al. (2007). However in this work, we consider the case $\lambda_{6}=\lambda_{7}=0$ but $m_{12}^{2}\neq 0$, which can allow FCNC processes at loop level Glashow and Weinberg (1977); Gunion and Haber (2003).

In CP-conserving 2HDM case, four realizations have been considered in the 2HDM to avoid FCNC at the tree-level, known as type I, type II, type III (called Lepton specific ) and type IV (also called Flipped) Barger et al. (1990); Grossman (1994). In the type I model, the $Z_{2}$ symmetry is an exact symmetry and the $\Phi_{1}$ doublet gives mass to all fermions. In the type II model, the $\Phi_{1}$ gives mass to leptons and down quarks while the $\Phi_{2}$ couples to up-type quarks. The MSSM has a type II Higgs sector. In the type III model, all quarks couple to the $\Phi_{2}$ doublet while all charged leptons couple to the $\Phi_{1}$ doublet. Finally, in the type IV model, down-type quarks couple to the $\Phi_{1}$ doublet while the rest of fermions couple to the $\Phi_{2}$ doublet. With such four arrangements in flavor space of 2HDM, the tree-level flavour changing neutral current can be suppressed safely.

From the initial eight degrees of freedom, if the $SU(2)$ symmetry is broken, we end up with two CP-even Higgs states usually denoted by $h^{0}$ and $H^{0}$, one CP-odd $A^{0}$, two charged Higgs boson, $H^{\pm}$ and three Goldstone bosons. After the electroweak symmetry breaking, the potential in Eq. (1) can be expressed with 7 independent parameters namely $m_{h^{0}}$, $m_{H^{0}}$,$m_{A^{0}}$, $m_{H^{\pm}}$, $\tan\beta=v_{2}/v_{1}$, $\sin(\beta-\alpha)$ and $\lambda_{5}$ (or equivalently $m^{2}_{12}$ ). The angle $\beta$ is the rotation angle from the weak gauge eigenstates to the mass eigenstates in the CP-odd and charged Higgs sector. The angle $\alpha$ is the corresponding rotation angle for the CP-even sector. From this potential, the triple Higgs couplings needed for the present work can be derived, which are functions of the physical parameters and are given in Eqs. (30-37) of the Appendix B.

The parameter space of the 2HDM is reduced by the following theoretical and experimental constraints:

Vacuum stability conditions that ensure the potential is bounded from below, where we use the conditions derived in Deshpande and Ma (1978). 2. 2.

Perturbative tree-level unitarity for scattering amplitudes of Higgs bosons and longitudinally components of gauge bosons. 3. 3.

The perturbativity of all quartic coefficients of the scalar potential , i.e $|\lambda_{i}|\leq 8\pi$ ($i=1,...,5$), 4. 4.

EW Precision Observables (EWPOs). Due to the contributions of extra Higgs bosons, the universal parameters $S$, $T$ and $U$ provide additional constraints on the mass splitting between these Higgs bosons. To implement the constraints from EWPOs, we consider the following values Olive et al. (2014); Haller et al. (2018) for $S$, $T$ and $U$: $\Delta S=0.05\pm 0.11$, $\Delta T=0.09\pm 0.13$ and $\Delta U=0.01\pm 0.11$. 5. 5.

Indirect experimental constraints from $B$ physics observables, which have been taken into account by using SuperIso Mahmoudi (2009) public code. Several important experimental values are tabulated in Table 1, which will be used to constrain parameters such as $\tan\beta$ and the charged Higgs boson mass of the 2HDM significantly.

The 2HDMC public code Eriksson et al. (2010) allows us to check all the listed constraints above. Furthermore, the proposed benchmark scenarios in this paper satisfy two more constraints: 1) the limits obtained from various searches for additional Higgs bosons at the LHC and other collider data, and 2) by the requirement that there exists a neutral scalar which should match the measured properties of the Higgs-like boson. We evaluate the first constraint with the public code HiggsBounds-5.3.2 Bechtle et al. (2010, 2011, 2012, 2014a, 2015), and the second constraint with the code HiggsSignals-2.2.3 Bechtle et al. (2014b); Stål and Stefaniak (2013); Bechtle et al. (2014c).

Here we stress that after the Higgs-like particle discovery, several theoretical studies have performed global-fit analysis for the 2HDM to pinpoint the allowed regions of parameter space both for a SM-like Higgs $h^{0}$ Ferreira et al. (2012a); Bernon et al. (2015); Dumont et al. (2014); Cheung et al. (2014); Eberhardt et al. (2013); Coleppa et al. (2014); Arhrib et al. (2016); Chiang and Yagyu (2013) as well as for a SM-like Higgs boson $H^{0}$ Bernon et al. (2016); Ferreira et al. (2012b). Before presenting our results, we would like to mention that we have performed a cross check of the results  Lopez-Val et al. (2010) for the subclass of Yukawa corrections considered there and found a perfect agreement.

To explore the effect of the new physics in the 2HDM in our following study, we propose four benchmark points (BPs) as shown in Table 2 based on the best-fit from the latest results of Higgs data using HiggsBounds and HiggsSignals public codes.

To demonstrate the impact of experimental data as well as theoretical constraints upon the parameter spaces, we perform a scan on $\sin(\beta-\alpha)/\cos(\beta-\alpha)$ and $\tan\beta$ plane. The results are shown in Figure (1) where we project the LHC constraints discussed above onto allowed regions at 95% Confidence Level (CL) (yellow color) as well as at 68% CL (green color) in the ($\sin(\beta-\alpha),\tan\beta$) plane for $h$ being SM-like for type-I (up left panel) and type-II (up right panel) 2HDM. In the lower plots, the heavier Higgs boson $H$ is the SM-like Higgs boson with the same coding color. The selected four benchmark points are given in Table 2, each of which corresponds to the red star in each plot of Figure (1).

About the mass of the charged Higgs boson, here we would like to point out that a light charged Higgs boson less than 200 GeV can be consistent with the current LHC Higgs data. As it is well known that in the framework of 2HDM-II and IV, for example, the measurement of the $b\to s\gamma$ branching ratio requires the mass of charged Higgs boson to be heavier than 580 GeV Misiak and Steinhauser (2017); Misiak et al. (2015) for any value of $\tan\beta\geq 1$. Such a limit is much lower for the other 2HDM types Enomoto and Watanabe (2016). In 2HDM-I and III, as long as $\tan\beta\geq 2$, there are allowed regions in the parameter space with a charged Higgs bosons as light as 100 GeV Enomoto and Watanabe (2016); Arhrib et al. (2017b), which are still consistent with all $B$ physics constraints as well as with LEP and LHC limits, as shown in the literature Aad et al. (2015b); Khachatryan et al. (2015a, b); Aad et al. (2013); Abbiendi et al. (2013); Akeroyd et al. (2017).

Due to the tiny Yukawa coupling of Higgs/Goldstone to electron, the overwhelming and leading order (LO) contribution for the process $e^{+}e^{-}\to Zh$ and $e^{+}e^{-}\to ZH$ in the 2HDM comes from the Feynman diagram given in Figure (2).

Then the total tree-level cross-section in 2HDM $\sigma^{0}$ can be expressed as

[TABLE]

where the tree-level cross section of the SM $\sigma^{0}_{SM}$ is defined as

[TABLE]

where $s$ is the collision energy, $a_{e}=-1,v_{e}=-1+4\sin^{2}\theta_{W},\alpha_{em}={e^{2}}/{4\pi},\lambda=\left(1-{m_{h}^{2}}/{s}-{m_{Z}^{2}}/{s}\right)^{2}-{4m_{h}^{2}m_{Z}^{2}}/{s^{2}}$. Here we have used $s_{W}$ and $c_{W}$ for $\sin(\theta_{W})$ and $\cos(\theta_{W})$, and $\theta_{W}$ is the Weinberg angle. With the relation $\alpha_{em}={\sqrt{2}m_{W}^{2}s_{W}^{2}G_{F}}/{\pi}$, it can be seen that our LO result in the SM is consistent with those in Refs. Kilian et al. (1996); Carena et al. (1996); Djouadi (2008) . The formula of tree-level cross sections for Higgs-strahlung process at $e^{+}e^{-}$ colliders in the 2HDM can be also found in Djouadi et al. (1993), where the MSSM was considered which has a built-in type-II 2HDM in the Higgs sector.

III One-loop renormalization and calculation

Calculations of higher-order corrections in perturbation theory in general lead to ultra-violet (UV) divergences. The standard procedure to eliminate these UV divergences consists in the renormalization of bare Lagrangian by redefining the couplings and fields. In the SM, the on-shell renormalization scheme is well elaborated Bohm et al. (1986); Hollik (1990); Denner (1993).

As shown in Ref. Denner (1993), the renormalization constant for charge $\delta Z_{e}$ is obtained as

[TABLE]

with

[TABLE]

and

[TABLE]

where the definitions for $\delta Z_{e}$, $\delta Z_{AA}$ and $\delta Z_{ZA}$ can be found in Eqs. (40) and (47). This corresponds to the running coupling constant obtained in the Thomson limit, taken as $\alpha_{em}(m_{e})$ or $\alpha_{em}(0)$. The vacuum polarization $\Pi(0)$, the first term in $\delta Z_{e}$, is sensitive to the hadronic contribution. Usually a non-perturbative parameter $\Delta\alpha^{(5)}_{\mathrm{hadron}}(m_{Z})$ (“5” here is the number of quark flavors, which means top quark is not included) is used to absorb the hadronic contribution, i.e. $\Pi(0)$ is modified as (the contributions from leptons are also separated for further discussion):

[TABLE]

with $\Delta\alpha^{(5)}_{\mathrm{hadron}}(m_{Z})=0.02764$ according to PDG data Tanabashi et al. (2018). Then $\delta Z_{e}$ is rewritten as

[TABLE]

In the following text and Table 3, we label the scheme defined in Eq. (8) as $\alpha_{em}(0)$ scheme. Meanwhile, whenever $\Pi(0)$ is mentioned, it refers to the one defined in Eq. (7).

Another two schemes labelled as $\alpha_{em}(m_{Z})$ and $\alpha_{em}(\sqrt{s})$ are defined below where the large logarithmic contributions of leptons are also absorbed into the redefinition of running coupling constant Denner (1993); Sun et al. (2017)

[TABLE]

with

[TABLE]

And the running coupling constant is defined as

[TABLE]

In the $\alpha_{em}(m_{Z})$ scheme, we take $\mu=m_{Z}$, while in the $\alpha_{em}(\sqrt{s})$ scheme, we take $\mu=\sqrt{s}$. Results in these two schemes will be independent of $\log(m_{e})$. In the following, we use the $\alpha_{em}(m_{Z})$ scheme as our default choice. In the procedure of renormalization using dimensional regularization, a scale $\mu_{r}$ is introduced. Usually this will generate a logarithmic term $\log(\mu_{r})$. As we use a complete on-shell scheme in both production processes, such terms which can be considered as an overall factor $(\mu_{r}^{2})^{\epsilon}$ in each individual part (counter terms in the modified minimal subtraction ($\overline{\mathrm{MS}}$) scheme do not have such a factor), will vanish upon the cancellation of UV and infrared (IR) divergences when summing over all the parts. This leads to $\mu_{r}$ independence of our results for the production processes. On the other hand, from Eq. (9) a new scale $\mu$ is introduced, which denotes the scale at which charge is renormalized, and our results are dependent on it. From now on, when we talk about renormalization scale, it denotes $\mu$ which is introduced in the renormalization of charge, not $\mu_{r}$.

The renormalization of 2HDM had been plagued by the issue of gauge-dependence of the mixing parameters, like the $\tan\beta$, as in the MSSM case where type II 2HDM is needed Kanemura et al. (2004, 2015, 2017). A reasonable and convenient scheme should be gauge-independent and process-independent and numerically stable as well Freitas and Stockinger (2002). Based on the renormalization scheme worked out by Fleischer and Jegerlehner in Ref. Fleischer and Jegerlehner (1981), which is now usually called FJ tadpole scheme, two groups find a way to fulfill such conditions. The first one is the $\overline{\mathrm{MS}}$ tadpole scheme (MSTS) Denner et al. (2016, 2017, 2018), where the mixing angles are renormalized using $\overline{\mathrm{MS}}$ subtraction. The other is the pinched tadpole scheme (PTS) Krause et al. (2016, 2017), where the pinch technique (see e.g. Ref. Binosi and Papavassiliou (2009)) is used to define gauge-independent counter terms for the mixing angles. More renormalization schemes have been proposed and examined numerically Altenkamp et al. (2017, 2018) and have been implemented in the HDecay package Krause et al. (2020).

In this work, we adopt the on-shell PTS described in Ref. Krause et al. (2016). On-shell here means that the on-shell scale is chosen in the self-energies during the renormalization of the mixing angles. More details about the renormalization scheme in this work can be found in Appendix C.

In our calculation, a small photon mass $\lambda$ is introduced to regularize the soft divergence. Meanwhile, two cutoffs, $\Delta E$ and $\Delta\theta$, are introduced to deal with the IR singularities in real correction processes. The 3-body phase space of real correction process $e^{+}e^{-}\rightarrow Zh\gamma$ is divided into three parts:

• •

soft ($S$) part: where the energy of photon $E_{\gamma}$ is smaller than $\Delta E$

• •

hard collinear ($HC$) part: where $E_{\gamma}\geq\Delta E$ and the angle between photon and the beam $\theta_{\gamma}$ is smaller than $\Delta\theta$

• •

hard non-collinear ($H\overline{C}$) part: the remaining, which is finite.

Thus next leading order (NLO) corrections can be expressed as

[TABLE]

Here $d\sigma_{V}$ denotes the virtual correction including loop diagrams and counter terms from renormalization.

The soft part is given by

[TABLE]

where $d\sigma^{0}$ denotes the tree-level differential cross section. The hard collinear part is obtained in the collinear limit as

[TABLE]

with $0\leq z\leq 1-\delta_{s}=1-2\Delta E/\sqrt{s}$, which is also related to the tree-level differential cross section $d\sigma^{0}$. $CT$ denotes the “counter term” from electron structure function, originated from the 2nd term in Eq. (28):

[TABLE]

with $0\leq z\leq 1$. Thus the combination of HC and CT parts can be expressed as

[TABLE]

Both the soft and virtual parts are obtained with FormCalc, while the other parts are obtained with the help of FDC Wang (2004). The $\lambda$ dependence has been checked when we combine the soft and virtual parts, more detailed checking in the SM can be found in Appendix A.

In Appendix A, we examine the dependence of the SM cross-section on those unphysical parameters, such as $\delta_{s}$ in Table 4 and $\Delta\theta$ in Table 5. From these results, it is observed that within a reasonable region where $\delta_{s}\ll 1$ and $m_{e}/\sqrt{s}\ll\Delta\theta\ll 1$ are satisfied, our results are indeed independent of those unphysical parameters introduced by the algorithm. Meanwhile, as we adopt the $\alpha_{em}(m_{Z})$ scheme, the cross-section should also be independent of the logarithm term of electron mass $\log(m_{e})$. This independence is shown in Table 6. After these checks, we apply our procedure to the 2HDM, i.e. we calculate the one-loop radiative corrections to the Higgs-strahlung process $e^{+}e^{-}\to Z\phi$.

Here we calculate the radiative corrections to the tree level $e^{+}e^{-}\to Zh^{0},ZH^{0}$ processes in 2HDM in the Feynman-’t Hooft gauge, including all the particles of the model in the loops. Counter-terms are constructed by specific renormalization conditions which allow us to cancel all the UV divergences with one-loop diagrams. Inserting these redefinitions into the Lagrangian, we find the following counter term for $h^{0}ZZ$ and $H^{0}ZZ$:

[TABLE]

The one-loop Feynman diagrams related to the processes $e^{+}e^{-}\to Z\phi_{i}$ are displayed in Figures (3-6), which are conveniently dubbed as the vertex corrections, the box contributions, self-energy corrections and the counter terms, respectively. In these Feynman diagrams, the labels $V,S,F$ denote all insertions of vector, scalar and fermionic states of the 2HDM. As we used the on-shell PTS in our renormalization scheme, we should also include all possible tadpole diagrams, which have been merged to corresponding counter term diagrams in Figure (6) and have not been shown explicitly.

There is one issue on the effects of the width of particles in the s-channel diagrams which need to clarify. For the processes that we consider here, the intermediate $Z$ boson is always far away from its mass shell as the collision energy is assumed to be 250 GeV or even higher. Thus it is plausible to neglect the effect of the width of Z boson in our calculation. Meanwhile, it is observed that contributions from initial state $e^{+}e^{-}\phi$ vertices vanish in the limit $m_{e}\to 0$ since both $e^{+}$ and $e^{-}$ are on-shell, then it is justifiable to neglect the widths of scalar bosons. Obviously, such an argument holds for the rest of scalar boson exchange diagrams in s-channel. As mentioned in the caption of Table 2, the ratio $\Gamma/m$ for all the new Higgs boson particles is less than $5\%$ and the narrow width approximation can hold. Therefore, in this work, for those s-channel diagrams, we simply neglect the widths of Z boson and other scalar bosons.

Finally, the total cross section at NLO, denoted as $\sigma^{NLO}$, is defined as the sum of LO cross section, $\sigma^{0}$, and the NLO corrections $\sigma^{1}$, i.e.

[TABLE]

where $\delta^{1}$ is defined as the ratio which measures the relative strength of next to leading order corrections over the tree-level result.

As described in Section 3.1 of Ref.Beenakker et al. (1991), the NLO electroweak corrections $\sigma^{1}$ can be safely grouped into two gauge-invariant parts:

• •

the ”QED” part, which includes all the diagrams which contain an extra photon attached to the LO diagrams, such as diagrams $v_{23}$, $v_{29}$ and $v_{52}$ in Figure (3) when the vector boson denotes a photon. Meanwhile, the photon’s contribution to the wave-function renormalization of electron is also grouped into this part.

• •

the ”weak” part, which contains all the remaining contributions.

In term of the convention introduced in Ref.Beenakker et al. (1991), we can divide $\sigma^{1}$ as

[TABLE]

Correspondingly, the $\delta$ defined in Eq. (19) can also be divided into two parts

[TABLE]

which represent the relative strength of weak and QED corrections, respectively.

IV Numerical results and discussions

Below we present the numerical results that are obtained from the analysis of the processes $e^{+}e^{-}\rightarrow Z+h/H$ at the one-loop level in both the SM and the 2HDM. We focus on the following two quantities: the ratio of the weak correction $\delta^{\mbox{1,weak}}$ and the ratio of the full one-loop corrections (including real emissions) to the leading order results $\delta^{1}$. ‡‡‡In the QED part, the extra photon can only attach to an initial electron pair in this process either in the SM or in the 2HDM. Thus $\delta^{\mbox{1,QED}}$ should be the same in both SM and 2HDM results.

In our numerical calculation, parameters in Table 2 are used as physical input in the Higgs sector, while in the gauge sector two different physical parameterizations are used:

• •

$\{\alpha_{em},G_{F},m_{Z}\}$: They are used in the scan of the 2HDM parameter space, which finally leads to Table 2. They are also used in determining the allowed range of $\lambda_{5}$ and in the calculation of decay part.

• •

$\{\alpha_{em},m_{W},m_{Z}\}$: They are used in the calculation of production processes, as we have used on-shell condition for the renormalization of $W$ boson self-energy there.

The values of above parameters are taken from PDG Tanabashi et al. (2018) as $\alpha_{em}=1/137.036$, $G_{F}=1.16638\times 10^{-5}$ Gev*-2*, $m_{Z}=91.1876$ GeV and $m_{W}=80.385$ GeV. It should be noted that the Vacuum Expectation Value (VEV) $v$ is different in these two cases. In the former case, it is determined as $v=1/\sqrt{\sqrt{2}G_{F}}\approx 246.220$ GeV. While in the latter case, it is determined as $v=m_{W}\sqrt{1-m_{W}^{2}/m_{Z}^{2}}/\sqrt{\pi\alpha_{em}}$ and varies as the scale of $\alpha_{em}$ changes. In the $\alpha_{em}(m_{Z})$ scheme, it is $v\approx 243.137$ GeV.

The total NLO cross-section for $e^{+}e^{-}\rightarrow Zh_{SM}$ including the real emissions can be found in Table 3. The results in the $\alpha_{em}(0)$ scheme still depend on $\log(m_{e})$ while the other two are independent. At the LO, the cross-section $\sigma^{0}$ in these three schemes reads as 223.12 fb, 252 fb, and 257.68 fb, respectively, where the maximal difference is 34.56 fb. While at NLO, the cross-section $\sigma^{NLO}$ reads as 230.25, 228.93, and 228.05 fb, where the maximal difference is 1.20 fb. From these numbers, we can see that scale dependence has been greatly improved at NLO, which can be exposed more clearly in Figure (7). For the sake of comparison with other literature, our definition of the weak part agrees to the one defined in Refs. Sun et al. (2017); Gong et al. (2017) and we also find that our numerical results in the SM agree with those given in Refs. Sun et al. (2017); Gong et al. (2017).

In Figure (7), the dependence on the renormalization scale $\mu$ and the collision energy of the cross-section in the SM are explicitly shown. In Figure (7a), it is observed that the leading order results can change drastically when the unphysical renormalization scale $\mu$ varies from 0 to $\sqrt{s}$, i.e. the difference between the results of $\alpha_{em}(0)$ and $\alpha_{em}(\sqrt{s})$ can reach up to $15\%$ as given in Table 3. In contrast, the scale dependence of the results is significantly reduced at the next-leading order. For example, the difference between the NLO results of $\alpha_{em}(0)$ and $\alpha_{em}(\sqrt{s})$ is $1\%$ or so. In Figure (7b), the lineshapes of the LO and NLO cross sections varying with the collision energy $\sqrt{s}$ are shown. The largest difference between different renormalization scales occurs near the threshold.

In Figure (8), we compare the results of leading order and those of the full NLO in the SM. In Figure (8a), the distribution of the transverse momentum of the Higgs boson in the process $e^{+}e^{-}\to Zh$ in the SM is shown where the collision energy $\sqrt{s}$ is taken as 250 GeV. When the angle of the outgoing Higgs boson with reference to the incoming electron direction is defined as $\theta$ in the laboratory frame, the transverse momentum of the Higgs boson is defined as $P_{t}(h)=|\overrightarrow{P}(h)|\sin\theta$, where $|\overrightarrow{P}(h)|$ denotes the magnitude of 3-dimensional momentum of Higgs boson in the laboratory frame. The cutoff of $P_{t}(h)$ is determined by the total collision energy. When the collision is set to be 250 GeV, the maximum of $P_{t}(h)$ is $62.12$ GeV which is determined by kinematics. We find that the line-shapes of the LO and full next-to-leading order (NLO) results are similar except a global scaling factor. In Figure (8b), the energy dependence of the cross-section is shown.

As a cross check, we compare our LO cross section of the SM with that computed using Whizard Kilian et al. (2011) (say in the $\alpha_{em}(m_{Z})$ scheme), and a good agreement is found. Moreover, it is observed that although the lineshape of $P_{t}(h)$ given in Figure (8a) is different from Figure (2b) in Chen et al. (2018) where the distribution is a normalized one with the effects of initial state radiation (ISR) included, a good agreement is found when we also normalize the lineshape of $P_{t}(h)$ given in Figure (8a) with the tree-level results of Figure (2b) of Chen et al. (2018) simulated using MadGraphAlwall et al. (2014).

Based on the results of the SM, since the NLO results have less dependence on the renormalization scale $\mu$, in the following study on the 2HDM, we will fix a renormalization scale $\mu=m_{Z}$ to present the results.

Now we start to present the results of the four benchmark points of 2HDM. The parameter of $\lambda_{5}$ is related to the triple Higgs boson couplings, which is given as Gunion and Haber (2003)

[TABLE]

Once $\lambda_{5}$ is fixed, all the triple Higgs boson couplings are fixed, as shown in Eqs. (30-37) of Appendix B. To demonstrate the effects of triple Higgs boson self-couplings, we try to vary it near its value for each benchmark point in Table 2. From now on, we treat each BP as a scenario. In each scenario, we fix the mass spectra, $\alpha$ and $\tan\beta$ as given in Table II, and let $\lambda_{5}$ vary. Due to the constraints from vacuum stability, unitarity and perturbativity to the theoretical parameters, the parameters $\lambda_{5}$ for 4 benchmark scenarios are determined to be in the ranges [0.545,0.583], [-0.48, 0.37], [-1.32,-0.30] and [-6.50,-5.74], respectively. We present our results in 2HDM with three typical values of $\lambda_{5}$, the upper bound, lower bound and the one in the benchmark point.

In Figure (9), the ratios $\delta^{1,weak}$ and $\delta^{1}$, which defined in Eq. (21), are shown, where the results are obtained in the $\alpha_{em}(m_{Z})$ scheme and label ”$h$” denotes that the lighter CP-even Higgs boson is assumed to be the SM-like Higgs boson. $\delta^{1}$ in which both weak and QED corrections are included is only presented for the SM, as we have discussed before that $\delta^{1,QED}$ is exactly the same in both SM and 2HDM. It is found that the signs of weak part and QED part are different in the SM. The QED part tends to increase the total cross-section by $+0.5\sim+1\%$ as the collision energy increases.

Also, in the BP1-h and BP2-h scenarios, the difference among the curves with different values of $\lambda_{5}$ is small and insignificant. For BP1-h scenario, this may be attributed to the fact that the allowed range for $\lambda_{5}$ is too narrow and the triple Higgs boson couplings are small, as shown in Figure (15a). In contrast, for the BP2-h scenario, the allowed range is wide enough and the triple Higgs couplings can be large enough, as shown in Figure (15b), this small difference instead can be attributed to the decoupling effects of heavy Higgs boson. It is found that for BP2-h scenario, this difference between the prediction of the SM is small and typically less than $0.2\%$ at most, which means at NLO new physics plays a less important role. While for BP1-h scenario, it it observed that the difference between the new physics and the SM can reach $0.8\%$ or so from 300 GeV to 500 GeV.

In Figure (10), two benchmark points for an alternative interpretation of the SM-like Higgs boson are considered, in such a scenario the heavier CP even Higgs boson has a mass near 125 GeV. The detailed information on the mass spectra and parameters of these two benchmark points are presented in Table 2 and labeled as BP1-H and BP2-H scenarios. Again, three typical values of $\lambda_{5}$ are used to show the effects of new physics. We observe that the ratio for BP1-H scenario (the pink dashed curve) varies from $-13\%$ to $-23\%$ or so, while for BP2-H scenario it varies from $-13.8\%$ to $-22\%$. Also in this alternative interpretation, the ratio difference for different values of $\lambda_{5}$ (say $\lambda_{5}=-1.32$ and $\lambda_{5}=0.30$ for BP1-H scenario, and $\lambda_{5}=-6.5$ and $\lambda_{5}=-5.74$ for BP2-H scenario) can reach more than $5\%$. It is found that the weak corrections of 2HDM for these two benchmark points and the typical values of $\lambda_{5}$ have the same sign in the SM, which holds for both the type-I and type-II as well as for the process $e^{+}e^{-}\to Zh$ and $e^{+}e^{-}\to ZH$. The last but not the least, all the lines have a bump near the collision energy 350-400 GeV, which can be attributed to the contribution of top quark pair in the loop functions.

In order to expose the contribution of new physics, in Figures (11) and (12), we present the contribution of new physics $\Delta^{weak}_{Zh}$ and $\Delta^{weak}_{ZH}$, respectively with different values $\lambda_{5}$. The quantity $\Delta^{weak}$ is defined as

[TABLE]

which describes the contribution of new physics compared with the results of the SM. We have discussed before that $\delta^{1,QED}$ are exactly same in both the SM and 2HDM, thus $\Delta^{weak}$ is enough to show the differences. Obviously, in this formula, the pure contribution of the SM is subtracted, while the contribution of new physics and the interference terms between the new physics and the SM are counted.

In Figure (11), the quantities $\Delta^{weak}$ of benchmark scenarios BP1-h and BP2-h are demonstrated. From the figure in the left panel, it is observed that the deviations caused by the 2HDM can change from $-0.28\%$ to $-0.29\%$ near the threshold and can decrease to $-0.54\%$ and $-0.52\%$ when the collision energy increases to 1000 GeV. A bump near the $2m_{H^{\pm}}=356$ GeV can be seen for BP1-h scenario. From the figure in the right panel, the deviations are flat in the explored energy region and can change from $-0.03\%$ to $-0.37\%$ at the threshold and can decrease to $+0.02\%$ and $-0.30\%$ at 1 TeV. Meanwhile, the deviations are more sensitive to $\lambda_{5}$ in BP1-h scenario where new Higgs bosons are relatively lighter than those in BP2-h scenario, although the allowed range of $\lambda_{5}$ is narrower. For three cases of BP2-h scenario, the deviations show no obvious bumps within the explored collision energy.

In Figure (12), the quantities $\Delta^{weak}$ of the benchmark scenarios BP1-H and BP2-H are shown. In the figure of left panel for the BP1-H scenario, the deviations are around $-0.58\%$ at the threshold and are similar for both $\lambda_{5}=-0.32$ and $\lambda_{5}=-0.30$. For two cases with $\lambda_{5}=-0.32$ and $\lambda_{5}=-0.30$, when the collision energy increases to 1000 GeV, the deviations change to $-1.5\%$ or so. While the case with $\lambda_{5}=-1.32$, the deviation changes from $-5.2\%$ at the threshold to $-6.6\%$ at 1000 GeV. In the figure of right panel for the BP2-H scenario, the deviations can reach from $-1.6\%$ to $-5.8\%$ which depends upon the values of $\lambda_{5}$. When the collision energy increases to 1000 GeV, the deviations change to $0.1\%$ to $-4.4\%$.

For the case of BP1-H scenario with $\lambda_{5}=-1.32$, a little bump near the energy region $\sqrt{s}=2m_{H^{\pm}}=340$ GeV can be observed. In contrast, for the two cases of BP1-H scenario with $\lambda_{5}=-0.30$ and $\lambda_{5}=-0.32$, the new physics contributions show no clear structure with the increase of collision energy due to smaller triple Higgs couplings. For the BP2-H case with $\lambda_{5}=-6.50$, an obvious dip near the energy region $\sqrt{s}=m_{A^{0}}+m_{h^{0}}=695$ GeV can be attributed to the contribution of $A^{0}$ and $H^{0}$ via the diagrams like $v_{12}$, $v_{14}$, $v_{15}$ and $v_{16}$, in Figure (3).

Except for the differential cross-sections of $e^{+}e^{-}\to Z\phi_{i}$ for four benchmark points in the 2HDM, it might also be interesting to explore the effects of new physics to the branching fractions of Higgs boson decays, which can be measured to a precision of $1\%$ or less at future Higgs factories.

The decay width of neutral Higgs boson in a model including one-loop correction can be formulated as

[TABLE]

where $N_{C}$ is number of color, $G_{F}$ is the Fermion constant, $m_{f}$ denotes the fermion’s mass, $\beta_{f}$ is the speed of fermion in natural units, and $\xi^{\phi}_{f}$ is the ratio of couplings of neutral Higgs boson in the model over those of the SM. The $\widehat{Z}_{\phi}$ refers to the finite wave function renormalization of the external $\phi$ while $\Delta{\cal M}_{1}$ is the amplitude of the one-loop vertex diagrams of $\phi\to f\bar{f}$. The quantity $\Delta r$ can be found in Denner (1993), which is obviously model dependent and renormalization scheme dependent. As this part is a work following Ref. Arhrib et al. (2004), the same renormalization scheme as in Ref. Arhrib et al. (2004) is used. Our results for the decays are renormalization scale dependent since the $\overline{\mathrm{MS}}$ scheme is used in the renormalization of Higgs boson wave functions. We have chosen the scale as $\mu_{r}=m_{\phi}$, the mass of decaying particle. For the study on $\phi\to f\bar{f}$ decay, as it is a work following Ref. Arhrib et al. (2004), the same renormalization scheme as in Ref. Arhrib et al. (2004) is used and will not be discussed here.

We can define a quantity $\Delta_{ff}(\phi)$ from Eq. (24) as a deviation from the predictions of the SM, which has the following form

[TABLE]

More details can be found in Arhrib et al. (2017a).

It is necessary to point out that the value of $\Delta_{ff}(\phi)$ in the SM is [math], and it is a measurement of new physics effects at quantum level. There are several studies which evaluate the BSM effects at the tree-level and one-loop radiative corrections to the decays of Higgs bosons Arhrib et al. (2004) and recently in Kanemura et al. (2020); Krause and Mühlleitner (2020).

In the type-I and type-II models, the couplings of $h^{0}$ to $b\bar{b}$ and $\tau^{+}\tau^{-}$ are proportional to the same factor $\xi^{h^{0}}_{f}$. It is $\cos\alpha/\sin\beta$ in type-I and $-\sin\alpha/\cos\beta$ in type-II. The couplings of $H^{0}$ to $b\bar{b}$ and $\tau^{+}\tau^{-}$ are also proportional to another same factor $\xi^{H^{0}}_{f}$. In type-I it is $\sin\alpha/\sin\beta$ while in type-II, it is $\cos\alpha/\cos\beta$. Obviously, there exist two special cases: 1) the couplings $|\xi^{h^{0}}_{f}|=1$ if $h^{0}$ corresponds to the SM-like Higgs boson; 2) the couplings $|\xi^{H^{0}}_{f}|=1$ if $H^{0}$ corresponds to the SM-like Higgs boson. Apparently, from tree-level expression, it is impossible to separate these cases from the SM. Below we explore whether it is possible to distinguish these cases by using the quantum corrections, i.e. by resorting to the information of $\Delta_{ff}(h^{0})$ and $\Delta_{ff}(H^{0})$.

We use HiggsSignals Bechtle et al. (2014b); Stål and Stefaniak (2013); Bechtle et al. (2014c) which incorporates Simplified Template Cross section (STXS) framework de Florian et al. (2016) to evaluate the $\chi^{2}$ of each point in the parameter space. The $\chi^{2}$ of HiggsSignals we used here includes two parts: 1) a signal strength part $\chi^{2}_{\mu}$ and 2) a mass peak part $\sum_{i}\chi^{2}_{m_{i}}$. In the signal strength part, when N observables for a given Higgs boson are considered, the $\chi^{2}_{\mu}$ (here $\mu$ denotes one of Higgs bosons) equals to $\chi^{2}_{\mu}=\sum_{\alpha=1}^{N}\chi^{2}_{\mu,\alpha}=(\widehat{\mu}-\mu)^{T}C_{\mu}^{-1}(\widehat{\mu}-\mu)$, where $(\widehat{\mu}-\mu)$ should be understood as a N-dimensional vector and $C_{\mu}$ is the $N\times N$ signal strength covariance matrix. Theoretical observables in the $\chi^{2}$ of HiggsSignals are defined as signal strengths $\widehat{\mu}$. When the final states of Higgs decay are specified, a signal strength is defined as $\widehat{\mu}=\sum_{i}(\sigma_{i}^{2HDM}Br^{2HDM})/(\sigma_{i}^{SM}Br^{SM})$, where $\sigma_{i}$ denotes the cross section of a production process and $Br$ denotes the decay branching fraction of the Higgs boson. The mass peak part takes into account all neutral Higgs bosons in a model, which is defined as $\chi^{2}_{m_{i}}=\sum_{\alpha=1}^{N}\chi^{2}_{m_{i},\alpha}=(\widehat{m}-m_{i})^{T}C_{m_{i}}^{-1}(\widehat{m}-m_{i})$, where N neutral Higgs boson is assumed. In our study of 2HDM, the $\chi^{2}$ mainly takes into account $\chi^{2}_{\mu}$ and the $\chi^{2}_{m_{i}}$ when a $m_{i}$ is assigned to be the SM-like Higgs boson. For example, when $H^{0}$ is assigned to be the SM-like Higgs boson in BP1-H and BP2-H scenarios, the corresponding $\chi^{2}_{\mu}$ and $\chi^{2}_{m_{H^{0}}}$ are computed.

In Figures (13) and (14), in order to examine $\Delta_{ff}(\phi)$ in the 2HDM, we use the mass parameters from Table II and $|\xi^{\phi}_{f}|=1$ to fix $\alpha$. Then we allow two parameters $\tan\beta$ and $\lambda_{5}$ to change. These two parameters are allowed to vary within the range $2<\tan\beta\leq 15$ and $-6\leq\lambda_{5}\leq 1$ while only points which satisfy all theoretical and experimental constraints are accepted. In Figure (13), we examine numerically $\Delta_{bb}(h)$ and $\Delta_{\tau\tau}(h)$ for D1-h and D2-h scenarios where $h^{0}$ is assumed to be the SM-like Higgs boson. We use D1-h and D2-h to label these two different scenarios which have the same mass spectra of the BP1-h and BP2-h in Table II, respectively. In this practice, it is obvious that the sign of $\cos(\beta-\alpha)$ is allowed to vary, therefore we assume it to be either positive or negative. In Figure (14), we examine $\Delta_{bb}(H)$ and $\Delta_{\tau\tau}(H)$ for D1-H and D2-H scenarios where $H^{0}$ is assumed to be the SM-like Higgs boson, likewise, D1-H and D2-H scenario share the same mass spectra like BP1-H and BP2-H, respectively. Meanwhile, the sign of $\sin(\beta-\alpha)$ is allowed to be both positive and negative, so we label the scenarios as $\pm$ D1-H and $\pm$ D2-H. Similarly we label other scenarios as $\pm$ D1-h and $\pm$ D2-h. In total, we can have 8 scenarios.

By scanning over the parameter space of 2HDM, we can compute its corresponding $\chi^{2}$ value. Obviously, the $\chi^{2}$ can be determined for each point in the parameter space and then the difference of $\Delta\chi^{2}=\chi^{2}-\chi^{2}_{min}$ can be computed. The total of number of freedom $\nu$ is defined as $\nu=n_{o}-n_{p}$, where $n_{o}$ is the number of observables taken from experiments measurements which is equal to 78 in this work and $n_{p}$ is the number of model parameters in 2HDM which is equal to 7. For a general fit, as shown in Figure (1) and the results provided in Table 2, the number of freedom is $78-7=71$. Now, for the specific fits shown in Figures (13) and (14), due to the constraint of $|\xi^{\phi}_{f}|=1$ and that we have fixed four mass parameters, there are only two free parameters. Therefore, the degree of freedom is $78-2=76$.

In Figures (13) and (14) , we show points that have satisfied the condition $\Delta\chi^{2}\leq 11$ ($<3\sigma$ for a two-parameter fit). For all 8 scenarios (labelled as $\pm$Dh-1, $\pm$Dh-2, $\pm$DH-1, and $\pm$DH-2), the main $\Delta\chi^{2}$ is mainly the contributions of $\chi^{2}_{\mu}$ from the signal strength observables since $\chi^{2}$ of the mass part is tiny and negligible. Plots are shown up to 99% C.L. It is remarkable that our results are in good agreement with previous publication Arhrib et al. (2017a); Kanemura et al. (2018). There are a few comments :

• •

A positive correlation between $\Delta_{\tau\tau}(h)$ and $\Delta_{bb}(h)$ is observed for both $\pm$D1-h and $\pm$D2-h, as demonstrated in Figure (13), where color bar indicates regions close to the minimum $\chi^{2}$. The best-fit point is indicated by a red star in each of these figures.

It turns out that in the D1-h scenarios, the corrections to $\Delta_{bb}(h)$ and $\Delta_{\tau\tau}(h)$ depend on the sign of $\cos(\beta-\alpha)$, i.e. when $\cos(\beta-\alpha)>0$ ($\cos(\beta-\alpha)<0$) the corrections are positive (negative). It can deviate from the prediction of the SM by $+15\%$($-6\%$) for the best fit point. In contrast, it is opposite for the BP2-h scenario and allowed to deviate from the prediction of the SM by less than $-6\%$($+5\%$).

• •

An interesting observation is that due to the positive correlation between $\Delta_{bb}(h)$ and $\Delta_{\tau\tau}(h)$, it is possible to determine the sign of $\cos(\beta-\alpha)$ by using a precise measurements of the branching fractions of $h\to b\bar{b}$ and $h\to\tau\tau$. As pointed out in Refs. Dawson et al. (2013); Zeppenfeld et al. (2000); Gianotti and Pepe-Altarelli (2000); Arhrib et al. (2017a), if the couplings of Higgs to fermions can be determined up to $8\%$, it becomes possible to discriminate the types of Yukawa interactions.

• •

At the CEPC or the ILC Dong et al. (2018); Arbey et al. (2015); Bambade et al. (2019), the precision of Higgs coupling measurements can reach typically ${\cal O}$(1)$\%$, and there is no doubt on their capability to discriminate the types of Yukawa interactions, like 4 scenarios introduced here.

The deviations of H decaying into $b\bar{b}$ and $\tau\tau$ are presented in Figure (14), where the correlations of $\Delta_{bb}(H)$ and $\Delta_{\tau\tau}(H)$ for $\pm$D1-H and $\pm$D2-H scenarios are demonstrated. One remarkable difference from Figure (13) is that the allowed deviation of $\Delta_{bb}(H)$ and $\Delta_{\tau\tau}(H)$ from those of the SM is smaller than $2.5\%$ ($-3\%$) for the best fit point. Such a deviation is even smaller (less than $1\%$) as shown by the BP2-H in the D2-H scenarios. Such deviations are quite distinctive due to Yukawa structure in these four scenarios. In other words, when a deviation is observed, an upper limit can be set on the mixing angles $\cos(\beta-\alpha)$ or $\sin(\beta-\alpha)$.

V Discussions and Conclusion

In this work, we have proposed 4 benchmark scenarios of the 2HDM after taking into account the current Higgs data from the LHC and evaluated the radiative corrections to the process $e^{+}e^{-}\to Z\phi$ in the SM and these 4 benchmark scenarios up to one-loop level.

In the Monte Carlo simulation results Mo et al. (2016) obtained using Whizard, it is found that the ISR can decrease the total LO cross-section of $e^{+}e^{-}\to Zh$ in the SM by more than $10\%$ when the contribution of high order logarithms are resummed. In contrast, the real emission in this work can increase the total LO cross-section by only a factor $+0.5\%$. This simply manifests the difference between the resummed and the fixed order results. Meanwhile, we find that the 1-loop weak correction of the SM can reduce the total NLO cross section by $9.15\%$, as demonstrated in Table 3.

As an estimation, we can express the total cross section with ISR as

[TABLE]

which had been proposed by the reference Kuraev and Fadin (1985) when $\frac{\Delta E}{E}<<1$ ($E$ denotes the energy of electron/positron, and $\Delta E$ is the energy of soft photon), where $R(1-\frac{\Delta E}{E},s)$ denotes the convolution of incoming fluxes of electron and positron, i.e. ISR, and ${\bar{\sigma}^{LO}}(s)$ denotes the cross section without weak corrections. It should be pointed out that ${\bar{\sigma}^{W,LO}}$ has included high order QED corrections but no weak corrections. The numerical results given in Mo et al. (2016) demonstrated that $R(1-\frac{\Delta E}{E},s)$ is roughly equal to $0.90$.

When the one-loop weak correction is taken into account, the total cross section at NLO with ISR can be expressed as

[TABLE]

where the $R(1-\frac{\Delta E}{E},s)$ (i.e. ISR) is the same as in the LO case. Our results demonstrate that weak correction can reduce ${\bar{\sigma}^{W,LO}}$ by a factor $9.15\%$ with $\mu=m_{Z}$. We can use the NLO cross section ${\bar{\sigma}^{W,NLO}}=228$ fb with $\mu=m_{Z}$ to estimate the total cross section with both ISR and weak corrections, which yields $\sigma^{NLO}(s)=228\times 0.9=205.2$ fb, which is $18.6\%$ smaller than the tree level cross section in the $\alpha_{em}(m_{Z})$ scheme given in Table 3.

Here we add one more comment on the ISR. As it is well-known that due to the soft-collinear divergences of QED, the large logarithmic terms in the form of $\alpha_{em}^{n}\log^{n}\frac{s}{m_{e}^{2}}\log^{m}\frac{m_{e}}{\lambda}$ (here $\lambda$ denotes the mass of photon or IR cutoff) from the ISR should be resummed. The effects of ISR at electron-positron colliders can affect cross-sections of physics processes significantly Berends and Kleiss (1985); Jadach et al. (1999). In Whizard Kilian et al. (2011), the effect of ISR to all orders Gribov and Lipatov (1972); Kuraev and Fadin (1985); Skrzypek and Jadach (1991) has been implemented in the structure functions of incoming electron and positron. The ISR has been implemented in other Monte Carlo codes Barberio and Was (1994); Jadach et al. (1999, 2000); Skrzypek et al. (1996); Jadach et al. (2001) and applied to LEP I Barate et al. (2000) and LEP II Schael et al. (2013) experimental analysis. Therefore, it is expected that the results for the SM from Whizard with ISR is more accurate in capturing large logs. While for the study of the 2HDM, currently, ISR effect has not implemented in Whizard and therefore there is no alternative code for us to check it.

We compute the decays $\phi\to b\bar{b}$ and $\phi\to\tau^{+}\tau^{-}$ with $\phi=h^{0}$ ( $H^{0}$ ) in Dh-1 and Dh-2 (DH-1 and DH-2) by including the EW corrections. We have shown that in Dh-2 and DH-2 scenarios, the electroweak radiative corrections in these two decay processes are rather small due to the fact that the heavy states $A^{0}$, $H^{0}$ and $H^{\pm}$ have masses of order of 600 GeV while they could be sizeable for D1-h and D1-H, as shown in Figure (13) and Figure (14). Considering the recent progress in the gauge independent renormalization schemes Krause et al. (2016); Denner et al. (2016); Altenkamp et al. (2017), it will be interesting to evaluate the differences caused by different renormalization schemes.

Our results demonstrate that $e^{+}e^{-}$ colliders (especially the Higgs factories with $\sqrt{s}=250$ GeV) can offer us the potential to distinguish various 2HDM models by looking at quantum effects in Higgs observables. Except for performing precision measurements of the SM-like Higgs boson, linear colliders also have the potential to discover new physics. For example, it can directly produce the light charged Higgs boson pair in BP1-h and BP1-H scenarios via $e^{+}e^{-}\to H^{+}H^{-}$ processes Komamiya (1988); Arhrib and Moultaka (1999). The energy scan of $e^{+}e^{-}$ colliders can also help to detect the mass spectra of Higgs bosons and even triple couplings of Higgs bosons, as shown in Figure (11) and Figure (12). If nature chooses 2HDM as the new physics beyond the SM at TeV region, and it would be hopeful to probe parameters of the Higgs potential sector by using the data of productions and decays. The results of the proposed benchmark points show that precision measurement on the cross-sections at different collision energy is helpful to explore the mass spectra and Higgs couplings of the 2HDM.

Acknowledgments

We thank Abdesslam Arhrib and David Lopez-Val for helpful discussion at the early stage of this work. This work is supported by the Moroccan Ministry of Higher Education and Scientific Research MESRSFC and CNRST: Project PPR/2015/6. RB was supported in part by the Chinese Academy of Sciences (CAS) President, International Fellowship Initiative (PIFI) program (Grant No. 2017VMB0021). Q.S.Yan is supported by the Natural Science Foundation of China under grant No. 11475180 and No. 11875260. B. Gong is supported by the Natural Science Foundation of China No. 11475183.

Appendix A The IR beheavior of $e^{+}e^{-}\to Zh$ in the SM

The dependence of the corrections at NLO, $\sigma^{1}$, on $\Delta E$ and $\Delta\theta$ are shown in Tables 4 and 5. In Table 4, the dependence is seen in a wide range and we choose $\delta_{s}=10^{-3}$ as our default choice. In Table 5, the result becomes cut dependent when $\Delta\theta$ is smaller than $10^{-4}$. It is because the approximation in Eq. (14) demands $\Delta\theta\gg m_{e}/\sqrt{s}\sim 2\times 10^{-6}$. Thus we choose $\Delta\theta=10^{-3}$ in our setting.

One-loop radiation corrections include collinear singularities which can be infinite when $m_{e}$ goes to zero, and they become terms proportional to $\log(m_{e})$ in this limit. Some of these terms are cancelled when virtual and real corrections are summed up, some of them can be absorbed into the redefinition of running coupling constant as mentioned above, but the collinear singularities can not be removed completely, which demands a careful manipulation.

To deal with these remaining collinear singularities, we used following fixed order electron structure function which can be derived from Eq. (11) of Ref. Kuraev and Fadin (1985):

[TABLE]

with

[TABLE]

being the regularized Altarelli-Parisi splitting function. One-loop structure function is used here to ensure the cancellation of collinear singularities, instead of using most-commonly-used resummed ones. The cancellation is shown in Table 6. We vary the mass of electron with a factor of $k$ from $2^{-4}$ to $2^{8}$ and find the result unchanged. Also, we can see that singular terms only appear in $\sigma_{V+S}$ and $\sigma_{SC}$ parts.

Appendix B Triple Higgs couplings in 2HDM

The Higgs potential of 2HDM determines the self-couplings among Higgs bosons. Among them, the triple Higgs couplings (THC) can be parameterized as a function of the 2HDM parameters $m_{h^{0}}$, $m_{H^{0}}$, $m_{A^{0}}$, $m_{H^{\pm}}$, $\tan\beta$, $\alpha$ and $\lambda_{5}$. At tree level, these couplings are independent of the Yukawa types and they are given as follows

[TABLE]

where $v$ is the VEV and $m_{12}^{2}$ can be derived from $m_{A}$, $\beta$, $\lambda_{5}$ and $v$, according to the relation given in Eq. (22). We have used the notation $s_{\theta}$ and $c_{\theta}$ as short-hand notations for $\sin(\theta)$ and $\cos(\theta)$, respectively. The mixing angle $\beta$ is defined by via $\tan\beta=v_{2}/v_{1}$.

As mentioned in the second paragraph of section IV, two different physical parameterizations have been used in this paper, which yields different values of $v$ and hence different THCs. Besides the overall factor $1/v$, the remaining terms of Eqs. (30-37) still depend on $v$, due to the fact that $m_{12}^{2}$ is obtained via Eq. (22) and terms proportional to $\lambda_{5}v^{2}$ will appear. In order to further investigate the dependence, we rewrite the THCs as:

[TABLE]

where the dimensionless couplings $\hat{\lambda}_{i}$ are just the terms in the curly brackets of Eqs. (30-37) divided by (246.220 GeV)2.

In Table 7, all the $\hat{\lambda}_{i}$ in the two cases are tabulated, as well as corresponding values of $m_{12}^{2}$. It can be seen from the table that in the case of BP2-h (BP2-H), some of the THCs are quite large compared with $h^{0}h^{0}h^{0}$ ($H^{0}H^{0}H^{0}$). Although these THCs are sizeable, we have checked that they still respect the conditions of perturbativity and the unitarity of scattering amplitudes of Higgs bosons in the model, as listed on page 3.

On the other hand, although the deviation between the two VEVs are smaller than 1.5%, it is observed that the values of $\hat{\lambda}_{i}$ can be significantly different, such as $\hat{\lambda}_{H^{0}H^{\pm}H^{\pm}}$, $\hat{\lambda}_{H^{0}A^{0}A^{0}}$ and $\hat{\lambda}_{H^{0}H^{0}H^{0}}$ in the case of BP1-h, and $\hat{\lambda}_{h^{0}h^{0}h^{0}}$, $\hat{\lambda}_{H^{0}h^{0}h^{0}}$, $\hat{\lambda}_{h^{0}H^{\pm}H^{\pm}}$ and $\hat{\lambda}_{h^{0}A^{0}A^{0}}$ in the case of BP2-H, which indicates the corresponding THCs are very sensitive to the value of $v$. As $v$ only appears in $\hat{\lambda}_{i}$ in terms of $\lambda_{5}v^{2}$, this also means those THCs are very sensitive to the value of $\lambda_{5}$ (a deviation smaller than 3% will change the corresponding $\hat{\lambda}_{i}$ drastically). Usually the effects caused by the different values of $v$ are regarded as higher order corrections. But from Table 7 we can see that in some cases this will lead to totally different theoretical predictions, e.g. the decay width of $H^{0}\rightarrow A^{0}A^{0}$ in the case of BP1-h. Fortunately in this work, this difference in $\hat{\lambda}_{i}$ does not change our results very much (see also the discussion below).

In Figure (15), we show the dependence of $\hat{\lambda}_{i}$ upon the parameter $\lambda_{5}$ for each of benchmark mass spectra and fixed parameters $\alpha$ and $\tan\beta$ given in Table 2, where $v=243.137$ GeV is used. In the benchmark points BP1-h(H) and BP2-h(H) where $h^{0}$($H^{0}$) is assumed to be the SM-like Higgs boson, the $\hat{\lambda}_{i}$ related to $h^{0}$($H^{0}$) are depicted as solid lines while the others are depicted by dashed lines as they can also contribute via the renormalization of wave functions or mixing angles. We have merged several couples of lines in the figure since they are identical or very close to each other, which can be easily learned from Eqs. (30-37) and the mass spectra in Table 2.

The range of $\lambda_{5}$ between two vertical dashed lines denotes the allowed region for each benchmark point. It is observed from the figure that the gradients of some lines are quite large, thus the corresponding $\hat{\lambda}_{i}$ change drastically in the allowed range of $\lambda_{5}$, such as $\hat{\lambda}_{h^{0}h^{0}h^{0}}$ in the case of BP1-H which changes from $7.5$ to $-0.05$. As we discussed before, the dependence of $\hat{\lambda}_{i}$ on $v$ can be converted into the dependence on $\lambda_{5}$. This explains the discrepancy of some $\hat{\lambda}_{i}$ caused by the different values of $v$ in Table 7. Also, it can be estimated that the effect caused by different values of $v$ (after removing the overall factor) is much smaller than the variation of $\lambda_{5}$ in the allowed range, since the difference between the values of $v$ can be taken as a deviation smaller than $3\%$ on $\lambda_{5}$.

We have discussed the dependence of THCs on $v$. In fact, the allowed ranges of $\lambda_{5}$ also depend on it. In Table 8, we show the allowed ranges of $\lambda_{5}$ for all the BPs. As mentioned in Section IV, we determine the allowed ranges of $\lambda_{5}$ by fixing other parameters and using the constraints from vacuum stability, unitarity and perturbativity. In this procedure we have taken $v=246.220$ GeV, same as the one used in the scan of the 2HDM parameter space. However in the calculation of production processes, $v=243.137$ GeV is used.

It can be observed from Table 8 that both the upper and lower bounds of $\lambda_{5}$ will change when a different $v$ is used. The largest deviation is found in the lower bound of BP2-h, which is about 10%. This can not be simply explained by an overall factor related to $v$ only, and needs further investigation.

The constraints used to determine the ranges of $\lambda_{5}$ can be expressed as several dimensionless inequalities involving $v$ and the 2HDM parameters: $m_{h^{0}}$, $m_{H^{0}}$, $m_{A^{0}}$, $m_{H^{\pm}}$, $\tan\beta$, $\alpha$ and $\lambda_{5}$. Besides an overall factor $1/v^{2}$, the only dependence on $v$ in these inequalities is the term $\lambda_{5}v^{2}$, similar to the case of THCs. This means we can remove part of the effects from $v$ by looking into $\lambda_{5}v^{2}$ or $m_{12}^{2}$. In Table 8, the ranges of $\lambda_{5}v^{2}$ and $m_{12}^{2}$ are also presented. It is found that the upper bounds of $\lambda_{5}v^{2}$ and $m_{12}^{2}$ are almost same even different values of $v$ are used, while the lower bounds are quite different. We have checked the constraints related to each bound and found that all the upper bounds are determined by the constraints from vacuum stability while all the lower bounds are determined by the constraints from unitarity. The constraints from vacuum stability are homogeneous so that they are not affected by the choice of $v$. On the other hand the constraints from unitarity are inhomogeneous, and some of them are even nonlinear. This leads to more complicated $v$ dependence of $\lambda_{5}v^{2}$ and $m_{12}^{2}$. From the table we can see that the effects of different $v$ on the lower bounds can be negligible such as $m_{12}^{2}$ in BP1-h, but they can also be large such as $m_{12}^{2}$ in BP2-H.

From above analysis we can conclude that deviations of the upper bounds of $\lambda_{5}$ can be easily understood as by multiplying the ratio of two different $v^{2}$, while the lower bounds are more complicated. What’s more, together with the situation of THCs, it hints that it is better to choose independent physical parameters with same mass dimension.

Appendix C On-shell pinched tadpole scheme in the 2HDM

In this appendix we present more details about the on-shell PTS Krause et al. (2016, 2017), which is used in the calculation of the production $e^{+}e^{-}\rightarrow Z\phi$. For the study on $\phi\to f\bar{f}$ decay, as it is a work following Ref. Arhrib et al. (2004), the same renormalization scheme as in Ref. Arhrib et al. (2004) is used and will not be discussed here.

C.1 Definition of renormalization constants

In the following we introduce renormalized quantities and renormalization constants in the on-shell PTS, where an index 0 is used to label bare quantities.

• •

Gauge and fermion sectors: things are similar as in the SM for these two sectors. Renormalization of gauge boson and fermion masses are given by

[TABLE]

and the electric charge by

[TABLE]

Meanwhile, renormalization of the gauge fields are defined as

[TABLE]

while for the fermion field, due to its left- and right-handed chirality, it has independent field renormalization constants:

[TABLE]

• •

Higgs sector: the renormalization is performed in physics Higgs basis. Renormalization of Higgs masses are defined as

[TABLE]

while renormalization of Higgs fields are given by

[TABLE]

For the mixing angles, they are defined as

[TABLE]

We have skipped the renormalization of $\lambda_{5}$ as it is not needed in our calculation. In order to expose the difference between the on-shell PTS and the other tadpole schemes, we also provide the renormalization of the tadpoles in the standard tadpole scheme, which are defined as given below:

[TABLE]

In contrast, the on-shell PTS is based on the FJ tadpole scheme Fleischer and Jegerlehner (1981), where the shifts of the VEVs $v_{1,2}\rightarrow v_{1,2}+\Delta v_{1,2}$ are introduced. We can choose the values of $\Delta v_{1,2}$ order by order to make renormalized tadpoles vanishing in the perturbation expansion. In such a procedure, we are able to appropriately allocate the tadpole contributions and define all counter terms in a gauge-independent way.

It is noteworthy that in the on-shell PTS, there is no need to define tadpole counter terms given in Eq. (72) (i.e. $\delta T_{H^{0},h^{0}}$ ) anymore. Instead, their roles are replaced by $\Delta T_{H^{0},h^{0}}$ which are caused by the shift of the VEVs and are given in Krause et al. (2016, 2017) as

[TABLE]

To avoid confusion, we have deliberately changed the convention of $\delta T$ in Krause et al. (2016, 2017) to $\Delta T$ to indicate that the origin and meaning of $\Delta T$ is different from $\delta T$ given in Eq. (72). Such a convention has also been adopted in Denner et al. (2016).

C.2 Renormalization conditions

$\Delta v_{1,2}$/$\Delta T_{H^{0},h^{0}}$ will appear in self-energies and as well as vertices. It has already been revealed in Refs. Krause et al. (2016); Denner et al. (2016) that the effect of $\Delta v_{1,2}$/$\Delta T_{H^{0},h^{0}}$ can be easily included by adding possible tadpole diagrams to generic one-particle irreducible diagrams. In this scheme, well-known on-shell conditions are used everywhere. The differences between well-known on-shell scheme and this scheme are :

Tadpole counter terms $\delta T_{H^{0},h^{0}}$ defined in Eq. (72) vanish in this scheme; 2. 2.

Self-energies $\Sigma(k^{2})$ are replaced with $\Sigma^{\mathrm{tad}}(k^{2})$ by adding possible tadpole contributions.

It should be pointed out that the derivatives of self-energies $\Sigma^{\prime}(k^{2})$ remain unchanged due to the fact that tadpole contributions are independent of external momentum $k$, thus most wave function renormalization constants are same.

In the following, we present more necessary information about this on-shell PTS. For more details, please refer to Ref. Krause et al. (2016).

In gauge and fermion sectors, same on-shell conditions as in the SM are used Denner (1993). Among all the renormalization constants in these sectors, three are different: $\delta m_{W}^{2}$, $\delta m_{Z}^{2}$ and $\delta m_{f}$. Only the first two are used in this work, they are given by:

[TABLE]

In Higgs sector, the renormalized self-energy of the Higgs field $\phi$ is the following finite combination of the unrenormalized self-energy

[TABLE]

while for the mixing one, it is

[TABLE]

with $(\phi_{1},\phi_{2})=(G^{\pm},H^{\pm}),(G^{0},A^{0})$ and $(H^{0},h^{0})$.

The renormalization conditions for the masses and wave functions in Higgs sector are given as:

• •

For the masses, on-shell conditions are used:

[TABLE]

These conditions ensure that the physical masses of Higgs bosons are pole masses of corresponding renormalized propagators, from which we obtain

[TABLE]

• •

For the wave functions, we require the residue of renormalized propagator of $\phi$ to be one:

[TABLE]

and the mixing self-energies vanish when one of the external particle is on shell:

[TABLE]

From the equations we obtain

[TABLE]

The renormalization of the mixing angles might be the most important part of this appendix. As the FJ tadpole scheme is applied, gauge-independent definitions for the counter terms $\delta\alpha$ and $\delta\beta$ are possible. In Ref. Denner et al. (2016), $\overline{\mathrm{MS}}$ subtraction is used, while in this work, we follow the procedure in Ref. Krause et al. (2016) where the pinch technique (see e.g. Ref. Binosi and Papavassiliou (2009)) is applied:

• •

First, the pinched self-energies $\overline{\Sigma}$ are obtained with the help of the pinched technique. It is found out to a sum of two parts:

[TABLE]

where $\xi$ stands for the gauge fixing parameters $\xi_{Z}$, $\xi_{W}$ and $\xi_{\gamma}$ of the $R_{\xi}$ gauge. As shown in Eq. (80), the first part has the same form as the tadpole self-energies evaluated in the Feynman gauge, and the second one is an additional contribution which is explicitly independent of the gauge fixing parameter $\xi$.

• •

The counter terms are then defined through above pinched self-energies with the on-shell scale, namely

[TABLE]

where the definition for $\delta\beta$ from CP-odd Higgs sector has been chosen.

More details about this on-shell PTS can be found in Ref. Krause et al. (2016). For completeness, we list here the additional parts in the pinched self-energies as the end of this appendix:

[TABLE]

where $B_{0}$ is the scalar two-point function ’t Hooft and Veltman (1979).