Next-to-leading-order corrections to the Higgs strahlung process from electron-positron collisions in extended Higgs models

We present the cross section for $e^{+}e^{-}\to hZ$ with arbitrary sets of electron and $Z$ boson polarizations at the full next-to-leading order in various extended Higgs models, such as the Higgs singlet model (HSM), the inert doublet model (IDM) and the two Higgs doublet model (2HDM). We systematically perform complete one-loop calculations to the helicity amplitudes in the on-shell renormalization scheme, and present the full analytic results as well as numerical evaluations. The deviation $\Delta R^{hZ}$ in the total cross section from its standard model (SM) prediction is comprehensively analyzed, and the differences among these models are discussed in details. We find that new physics effects appearing in the renormalized $hZZ$ vertex almost govern the behavior of $\Delta R^{hZ}$, and it takes a negative value in most cases. The possible size of $\Delta R^{hZ}$ reaches several percent under the theoretical and experimental bounds. We also analyze the deviation $\Delta R^{hZ}_{XY}$ in the total cross section times decay branching ratios of the discovered Higgs boson by utilizing the $\texttt{H-COUP}$ program. It is found that the four types of 2HDMs can be discriminated by analyzing the correlation between $\Delta R^{hZ}_{\tau\tau}$ and $\Delta R^{hZ}_{bb}$ and those between $\Delta R^{hZ}_{\tau\tau}$ and $\Delta R^{hZ}_{cc}$. Furthermore, the HSM and the IDM can be discriminated from the 2HDMs by measuring $\Delta R^{hZ}_{WW}$. These signatures can be tested by precision measurements at future Higgs factories such as the International Linear Collider.


Introduction
Since the discovery of the new particle with the mass of 125 GeV at the LHC in 2012 [1,2], it has turned out that its properties are in agreement with those of the Higgs boson in the standard model (SM) within theoretical and experimental uncertainties [3,4]. While no signal for new physics (NP) beyond the SM has been observed at the LHC up to now, there are phenomena that cannot be explained within the SM such as dark matter, baryon asymmetry of the universe and tiny neutrino masses. In addition to these phenomenological problems, there are conceptual problems in the SM such as the hierarchy problem, no unified description for the gauge group and the flavor structure and so on. Therefore, the SM must be replaced by a more fundamental theory.
While the Higgs boson was found, the structure of the Higgs sector remains unknown. There is no theoretical principle to insist on the minimal structure of the Higgs sector as introduced in the SM, and the possibility that the Higgs sector takes a non-minimal form is not excluded experimentally. Furthermore, such non-minimal Higgs sectors are often introduced in various new physics models, where the abovementioned problems are tried to be solved. Therefore, unraveling the structure of the Higgs sector is one of the central interests of current and future high-energy physics, and the direction of new physics can be determined by reconstructing the Higgs sector experimentally.
The discovery of additional scalar bosons is a clear evidence of extended Higgs sectors, and enormous efforts have been devoted to discover such new particles in a wide variety of the search channels [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. However, no observation of such new particles has been reported, leading to constraints on parameters of extended Higgs models such as masses and coupling constants. Direct searches of new particles are one of the key programs at the LHC as well as at the highluminosity LHC (HL-LHC) [22].
In addition to the direct searches, extended Higgs sectors can be indirectly explored by measuring various properties of the discovered Higgs boson such as cross sections, the width and decay branching ratios because mixings of Higgs bosons and/or radiative corrections of additional Higgs bosons would modify them from their SM values. If deviations from the SM are detected, the magnitude of deviations tells us upper limits on the mass scale of the second Higgs boson by taking into account theoretical consistencies [23][24][25]. In addition, the pattern of deviations gives us information on the structure of the Higgs sector such as the representation of the weak isospin, the number of Higgs fields and the structure of Yukawa interactions [23]. Therefore, the discovered Higgs boson is a probe of new physics.
Precision measurements of the properties of the discovered Higgs boson are also one of the main targets at current and future collider experiments. At the LHC, the Higgs boson couplings have been measured with typically order ten percent accuracy. Most extended Higgs models can accommodate this SM-like situation at the lowest order of the perturbation. Therefore, measurements with a few percent accuracies become more important where quantum corrections play an essential role, and it enables us to extract the signature of new physics. Accuracy of the measurement of Higgs boson couplings is expected to be improved at the HL-LHC and further significantly at future lepton colliders; e.g., the International Linear Collider (ILC) [26][27][28][29], the Future Circular Collider (FCC-ee) [30] and the Circular Electron Positron Collider (CEPC) [31].
In order to compare theoretical predictions with future precision measurements, theoretical calculations comparable with expected experimental accuracy are inevitable. Radia-tive corrections to the SM-like Higgs boson vertices have been studied in various Higgs sectors such as the model with a real isospin singlet Higgs field (HSM) [32][33][34][35], two Higgs doublet models (2HDMs) [36][37][38][39][40][41][42][43][44][45][46][47][48][49], the inert doublet model (IDM) [50,51] and so on. In order to see differences in the prediction among these models, it is quite important to calculate the renormalized SM-like Higgs boson vertices in a consistent and systematic way. The H-COUP program [52,53] enables us to evaluate the decay rates including higher-order corrections for the SM-like Higgs boson in the HSM, the IDM and the 2HDM with four types of Yukawa interactions classified under the softly-broken Z 2 symmetry [54,55]. Also, other numerical tools to evaluate the decay of the Higgs boson with radiative corrections are available; e.g., 2HDECAY [56] and Prophecy4f [57].
In this paper, we present the cross section for e + e − → h Z with arbitrary sets of electron and Z boson polarizations at the full next-to-leading order (NLO) in the HSM, the IDM and the 2HDMs. At the future lepton colliders, not only the decay properties of the discovered Higgs boson but also its production cross section can be precisely measured. Especially, e + e − → h Z is the dominant production process at the center-of-mass (CM) energy of 240-250 GeV, 1 and precise calculations of its cross section in various extended Higgs models are quite important. The NLO calculation has been performed in the SM [61][62][63], the IDM [64] and the 2HDM [65,66]. 2 We systematically perform complete NLO calculations to the helicity amplitudes in each extended Higgs model based on the on-shell renormalization scheme [37,43,71,72], and the full analytic results, as well as numerical evaluations, are presented. We comprehensively analyze the deviations in the cross section from the SM prediction in each model under the constraints of perturbative unitarity [34,[73][74][75][76][77] and vacuum stability [78][79][80][81][82][83], conditions to avoid wrong vacua [84][85][86][87][88] and experimental constraints. We discuss the differences in the predictions of the cross section among these models in detail. We also show correlations of the deviation in the cross section times the decay branching ratios of the SM-like Higgs boson from the SM predictions and discuss the discrimination of the extended Higgs models. This paper is organized as follows. In Sect. 2 we briefly introduce the HSM, the 2HDMs and the IDM. In Sect. 3 we present the helicity amplitude of e + e − → h Z including the electroweak (EW) radiative corrections. In Sect. 4 we show numerical results of deviations in the cross section from the SM prediction in each model. In addition, we show the correlations of the deviation in the cross section times the decay branching ratios of the SM-like Higgs boson from the SM predictions. Conclusions are given in Sect. 5. In Appendix, the input parameters and explicit formulae for the NLO calculations are presented.

Models with non-minimal Higgs sectors
In this section, we briefly review the HSM, the 2HDM and the IDM. Before moving on to the discussion on each extended Higgs model, we review concepts of general constraints on parameter spaces that are independent of models. We then define the extended Higgs models in order.

Constraints on extended Higgs models
First of all, the size of Higgs quartic couplings is constrained by the perturbative unitarity bound, which was originally introduced to obtain the upper limit on the mass of the Higgs boson in the SM [89,90]. Using the equivalence theorem [91], this bound requires that the magnitude of partial wave amplitudes for the elastic scatterings of two-body to two-body scalar boson processes, including the Nambu-Goldstone (NG) bosons, does not exceed a certain value. Each eigenvalue of the s-wave amplitude a i 0 should satisfy where ξ = 1 [89,90] or 1/2. We take ξ = 1/2 in this paper. Next, the vacuum stability bound provides an independent constraint on scalar quartic couplings. This bound requires that the Higgs potential is bounded from below in any direction with large field values. This condition is trivially satisfied in the SM by taking the Higgs quartic coupling to be positive. However, this bound requires a set of inequalities in terms of Higgs quartic couplings in extended Higgs models [79].
Furthermore, in extended Higgs models, wrong local vacua can appear in addition to the true vacuum giving the correct value of the Fermi constant G F . We have to avoid parameter regions where the depth of such wrong vacua becomes deeper than that of the true one. The condition to avoid the wrong vacua can be written by combinations of dimensionful and dimensionless parameters in the Higgs potential, and it provides an independent constraint from the above two constraints.
Apart from these theoretical constraints, we need to take into account bounds from experimental data. At the LEP/SLC experiments, various EW observables have been precisely measured such as the masses and widths of the weak gauge bosons. These precision measurements can be used to constrain the size of new physics effects which can enter into the two-point functions for weak gauge bosons. Such indirect effects, so-called oblique corrections, are conveniently parameterized by the S, T and U parameters [92,93], which are expressed in terms of two-point functions of the weak bosons. From the global fit of EW parameters [94], new physics effects on the S and T parameters under U = 0 are constrained by S = 0.05 ± 0.09, T = 0.08 ± 0.07, (2.2) with the correlation factor of +0.91 and the reference values of the masses of the SM Higgs boson and the top quark being m ref h = 126 GeV and m ref t = 173 GeV, respectively. Flavor experiments also provide important constraints on the parameter space in extended Higgs models, particularly in multi-doublet models. We will discuss these constraints in more detail in Sect. 2.3 about the 2HDM. Additional scalars have been directly searched at the LHC [5-21], and constraints are obtained for parameters of extended Higgs models. In addition, Higgs coupling measurements also give constraints especially on the mixing parameters in the HSM and the 2HDM [3,4]. The application of these constraints to each extended Higgs sector will be discussed in the following subsections.

Higgs singlet model
In the HSM, we have one isospin doublet scalar field with the hypercharge Y = 1/2 and one real singlet field S with Y = 0. We parametrize these scalar fields as where v is the vacuum expectation value (VEV) of the doublet field which is related to the Fermi constant by v = ( √ 2G F ) −1/2 246 GeV, while v S is the VEV of the singlet field. The component fields G ± and G 0 in the doublet field correspond to the NG bosons.
The most general Higgs potential is given by where all the parameters are real. We can take any value of v S without changing physical results [85], and we fix v S = 0 in the following discussion.
In the HSM, we have two physical neutral Higgs bosons. Their mass eigenstates are defined by introducing the mixing angle α as with the shorthand notation for the trigonometric functions as s θ ≡ sin θ and c θ ≡ cos θ . We define the domain of α as −π/2 ≤ α ≤ π/2. We identify h as the discovered Higgs boson with a mass of 125 GeV. After solving the tadpole conditions, the squared masses of neutral Higgs bosons are expressed as where the squared mass matrix elements M 2 i j (i, j = 1, 2) in (s, φ) T basis are given by with M 2 ≡ 2m 2 S . The parameters m 2 and t S are eliminated by using the stationary conditions for φ and s. We can replace the parameters λ, m 2 S and μ S with m 2 H , m 2 h and α by using Eqs. (2.6)-(2.8). We choose the following five parameters to be the free input parameters in the HSM: (2.10) and the two parameters m h and v are fixed by experiments. If the Higgs potential respects an exact discrete Z 2 symmetry, the t S , μ S and μ S terms are forbidden. This corresponds to the case with α → 0 and μ S → 0. The kinetic terms of scalar fields are given by where D μ is the covariant derivative for the Higgs doublet. The gauge-gauge-scalar type interaction terms are given by where g is the weak gauge coupling and g Z = g/c W with θ W being the weak mixing angle. The Yukawa interaction terms are the same as those in the SM and given by (2.13) where = iσ 2 * . In the above equation, Q L and L L are left-handed quark and lepton doublets, respectively, while u R , d R and e R are right-handed up-type quark, down-type quark and charged lepton singlets, respectively. The interaction terms for h and H with fermions are given by (2.14) We note that the SM-like Higgs bosons couplings with the SM particles are universally suppressed by c α as compared to those SM values. The parameters in the Higgs potential are constrained by perturbative unitarity, vacuum stability and the conditions to avoid wrong vacua. For the perturbative unitarity bound, there are four independent eigenvalues given in Refs. [34,73]. The necessary and sufficient conditions to satisfy vacuum stability are given by [78] For the conditions to avoid wrong vacua can be found in Refs. [84][85][86]. The one-loop corrected two-point functions for weak bosons are found in Ref. [95]. Imposing the constraint from the S and T parameters, we can obtain the upper limit on m H depending on the value of c α . Constraints on the mass of the additional Higgs boson and the mixing angle from the LHC data have been studied in Refs. [96,97].

Two Higgs doublet model
In the 2HDM, we have two isospin doublet scalar fields i with the hypercharge Y = 1/2. We parametrize these doublets as where v 1 and v 2 are the VEVs of two doublets with v = v 2 1 + v 2 2 . In the most general 2HDM, flavor-changing neutral currents (FCNCs) appear at tree level, and it is severely constrained by experiments. In order to avoid such FCNCs, we introduce a discrete Z 2 symmetry, where two doublets transform as 1 → 1 and 2 → − 2 [98,99]. One can introduce the soft breaking term of the Z 2 symmetry in the Higgs potential without spoiling the desirable property of the flavor sector.
The most general Higgs potential under the softly-broken Z 2 symmetry is given by Although m 2 3 and λ 5 are generally complex, we take them to be real and consider the CP-conserving case for simplicity. The mass eigenstates of the Higgs fields are defined as where tan β = v 2 /v 1 , and H ± and A are the charged and CP-odd Higgs bosons respectively, while H and h are the CP-even Higgs bosons. We define the domain of β to be 0 ≤ β ≤ π/2. We identify h as the discovered Higgs boson with a mass of 125 GeV. After solving two tadpole conditions for h 1 and h 2 , the squared masses of the charged and CP-odd Higgs bosons are given by where M 2 = m 2 3 /(s β c β ) which describes the softly-breaking scale of the Z 2 symmetry. The squared masses of the neutral Higgs bosons and the mixing angle β − α are given by where M 2 i j (i, j = 1, 2) are the squared mass matrix elements for the CP-even scalar states in the Higgs basis [100] (h 1 , h 2 )R(β): with λ 345 ≡ λ 3 +λ 4 +λ 5 . We define the domain of β −α to be 0 ≤ β − α ≤ π so that s β−α is always positive and c β−α has the opposite sign from M 2 12 [101]. The eight parameters in the Higgs potential are expressed by the following six input parameters: (2.26) and the two parameters m h and v are fixed by experiments.
In addition, we have a degree of freedom of the sign of c β−α . The kinetic terms of the Higgs doublets are given by In the mass eigenbasis of the Higgs bosons, the gauge-gaugescalar type interaction terms are given by The Yukawa interaction terms under the Z 2 symmetry are given by − Y e L L e e R + h.c., (2.29) where u,d,e are either 1 or 2 . As in Table 1, there are four types of Yukawa interactions according to the Z 2 charge assignment [102,103]. The interaction terms for the physical Higgs bosons with the fermions are given by and V ud is the Cabibbo-Kobayashi-Maskawa matrix element. The parameters in the Higgs potential are constrained by perturbative unitarity, vacuum stability and the condition to avoid wrong vacua. For the perturbative unitarity bound, there are twelve independent eigenvalues of the swave amplitude matrix [74][75][76][77]. The vacuum stability bound is sufficiently and necessarily satisfied by imposing the following conditions [79][80][81][82][83] Type-X (lepton specific) Type-Y (flipped) In addition, the wrong vacua can be avoided by taking M 2 ≥ 0 [87]. We thus only take the positive value of M 2 in the following discussion. 3 The expressions of the two-point functions for the weak bosons in the 2HDM are found in Refs. [106][107][108][109][110]. Imposing the constraint of the S and T parameters, we can find that the charged Higgs bosons and one of the neutral Higgs bosons should be approximately degenerate in their mass. This results from the constraint of the T parameter, and it can be satisfied if the Higgs potential respects the custodial symmetry [111,112]. Constraints on the parameters in 2HDMs from the LHC data have been discussed in Refs. [25,101,[113][114][115][116].
In the 2HDM, constraints from flavor experiments are important to be taken into account. These bounds particularly provide the lower limit on the mass of the charged Higgs boson m H ± depending on the type of Yukawa interaction and tan β. For example, from the B → X s γ data, m H ± has to be greater than about 800 GeV at 95% confidence level (CL) in the Type-II and Type-Y 2HDMs with tan β 2 [117], while O(100) GeV of m H ± is allowed in the Type-I and TypeX 2HDMs with tan β 2 [118]. Constraints on m H ± and tan β from various flavor observables are also shown in Ref. [119] in the four types of 2HDMs.

Inert doublet model
The contents of the Higgs bosons in the IDM are the same as those in the 2HDM. In the IDM, we assume an exact Z 2 symmetry and prohibit the m 2 3 term in the Higgs potential which softly breaks the Z 2 symmetry in the 2HDM. We also assume that the second Higgs doublet 2 does not develop the VEV to avoid the spontaneous breaking down of the Z 2 symmetry.
We parametrize the doublets as (2.32) The squared masses of the Higgs bosons are given by where M 2 ≡ m 2 2 . We note that in addition to the absence of the m 2 3 term, there is no tadpole condition for H . Therefore, the mass formulae for the scalar bosons are different from those in the 2HDM. We choose the following five parameters to be the free input parameters in the IDM: 37) and the two parameters m h and v are fixed by experiments. The parameters in the Higgs potential are constrained by perturbative unitarity, vacuum stability and the condition to guarantee the inert vacuum. The same conditions for perturbative unitarity and vacuum stability in the 2HDM can be applied to the IDM, because these bounds are given in terms of the scalar quartic couplings. In addition, there is the condition to guarantee the inert vacuum with ( 0 1 , Since the tadpole condition makes m 2 1 negative, and the vacuum stability condition constraints λ 1 and λ 2 to be positive, the condition given in Eq. (2.38) is satisfied by taking M 2 > 0. We refer to this condition as the one to avoid wrong vacua, according to the other two models discussed above.
For the constraints of the S and T parameters, we can use the same expressions as those in the 2HDM with s β−α = 1. As similar to the case in the 2HDMs, the charged Higgs bosons and one of the neutral Higgs bosons should be approximately degenerate in their mass in order to satisfy the constraint from the T parameter. In the IDM, constraints on the masses of the additional Higgs bosons from collider experiments are relatively weak since the additional Higgs bosons do not couple to the SM fermions. Constraints from the LEP and the LHC have been studied in Ref. [120] and Refs. [121,122], respectively. Dark matter constraints from relic density and direct detection also limit the parameter space; see, e.g., Refs. [122,123] for details. Finally, we summarize the scaling factors for the SMlike Higgs boson couplings to the weak bosons κ V and the fermions κ f in Table 2.

Electroweak corrections to the process e + e − → hZ
In this section, we define the notation for the process e + e − → h Z and discuss the helicity amplitudes based on the formfactor decomposition. We list relevant renormalized quantities for this process and give the formulae of form factors including the one-loop corrections. The differential cross section with arbitrary sets of electron and Z boson polarization is also presented. For the numerical evaluation, we use the SM input parameters given in Appendix A.

Helicity amplitudes and cross section
The process is depicted in Fig. 1. The momenta and helicities of the incoming electron and positron are denoted by ( p e , σ e ) and ( pē, σē), respectively. Correspondingly, (k Z , λ) is used for the outgoing Z boson, and k h is the momentum of the outgoing Higgs boson. The signs '+' and '−' of the variables σ e and σē refer to helicities +1/2 and −1/2, respectively. The helicity λ takes '±' or '0'. In the following discussion, we neglect the mass of the electron whenever it is possible.
The Mandelstam variables are denoted by The process e + e − → h Z with momentum and helicity assignments. The momenta p e and pē are incoming, while k h and k Z are outgoing and they satisfy s In the CM frame of the e + e − collision, the momenta of each external particle are with the beam energy E = √ s/2. We use the scattering angle θ between e − and Z boson;p e ·k Z = cos θ where the hat indicates the unit vector. The scattering angle θ is related to t and u via The helicity amplitudes for e + e − → h Z vanish for σ e = σē in the limit m e → 0 due to the chirality conservation. Therefore, we use σ = σ e = −σē for the nonvanishing amplitudes. The helicity amplitude M σ λ (s, t) can be decomposed into a set of basic matrix elements M i,σ λ and corresponding form factors (3.12) The basic matrix elements are given by where ε * μ (k Z , λ) is the polarization vector for Z boson, and j μ σ ( p e , pē) is the fermion current of the initial electron and positron, with the chirality projection operator In the CM frame, the first elements of basic matrix M 1,σ λ become The third elements M 3,σ λ are The six physical helicity amplitudes are given in terms of the form factors F i,σ by We denote the tree-and one-loop contributions to the helicity amplitude as The helicity-dependent differential cross section at NLO in EW is given by The helicity-dependent cross section σ (σ, λ; s) can be obtained by integrating Eq. (3.28) over the solid angle. In a realistic setup, one needs to introduce the degree of polarization of initial electron P e and positron Pē. We use the convention where a purely left-handed (right-handed) electron corresponds to P e = −1 (+1). The polarized differential cross section is given by The unpolarized cross section σ (λ; s) corresponds to The polarized cross section can be rewritten in terms of the helicity-dependent cross section as [124] σ (P e , Pē, λ; s) = 1 4 (1 − P e )(1 + Pē)σ (−, λ; s) where the effective polarization P eff and the left-right asymmetry A LR are defined as

31)
By using Eq. (3.30), one can easily evaluate the effect of the beam polarization from the helicity-dependent cross sections. Therefore, in the following discussion, we focus on the unpolarized and helicity-dependent cross section to exhibit analytical behaviors.

Tree-level contribution to the helicity amplitudes
At LO, only one diagram of Fig. 2 is relevant since the electron-Higgs coupling is proportional to m e , and it is negligible. The contribution of the tree-level diagram to the form factors is expressed as with the tree level h Z Z coupling The scaling factor κ Z is given in Table 2. If there is no mixing between CP-even scalars, the cross section in the extended Higgs models is the same as that in the SM at LO. The couplings g ± are defined by The coefficients f The lowest order differential cross section is given by In the left panel of Fig. 3, we show the helicity-dependent cross sections at LO as a function of the CM energy. In the numerical evaluation, we take G F as an input given in Appendix A and use the tree-level relation v = ( √ 2G F ) −1/2 . The solid and dashed lines correspond to the transversely (λ = ±) and the longitudinally (λ = 0) polarized Z bosons, respectively. The cross sections peak just above the threshold √ s = m Z + m h and monotonically decrease at higher energies. For energies well above the threshold, the longitudinally polarized Z boson dominates the cross section. This is due to the factor s/m 2 Z originated from the longitudinal polarization vector defined in Eq. (3.15). The cross section for the left-handed electron is larger than that for the right-handed electron because the left-handed electron more strongly couples to the Z boson than the right-handed electron (g − /g + ) 2 1.8. At √ s = 250 GeV, the unpolarized cross section is 242 fb, and the polarized cross section is 379 fb with (P e , Pē) = (−0.8, 0.3). Therefore, the beam polarization significantly changes the size of the cross section. We note that the predicted value of the LO cross section highly depends on the input schemes.
The angular distributions of the LO cross section at √ s = 250 GeV are given in the right panel of Fig. 3. They are determined by d 1 σ,λ (θ ). The cross section for the transversely polarized Z boson is proportional to 1 + cos 2 θ , and it takes maximal value in the forward-backward direction. On the other hand, that for the longitudinally polarized Z boson is proportional to 1 − cos 2 θ , and it vanishes in the forwardbackward direction and takes maximal value at cos θ = 0.

One-loop contributions to the form factors
As shown in Fig. 4, the one-loop contributions to the form factors F (1) i,σ consist of (a) the Z boson self-energy and the Z γ mixing, (b) the Zeē vertex correction, (c) the h Z Z and h Zγ vertex corrections, (d-e) the heē vertex correction and (f) the box diagrams. In addition, the renormalization factors of the weak gauge bosons are not to be unity in our renormalization scheme [71,72]. Therefore, we have the term − Re Z Z (m 2 Z ) /2 from the wave function renormalization of the on-shell Z boson. Furthermore, the EW correction to the Fermi decay constant r appears when one replaces the VEV in the tree-level amplitude with G F , since the tree-level relation between these two parameters is no longer valid at the one-loop level. This replacement corresponds to the resummation of universal higher-order leading corrections such as large logarithms from light fermion masses [62]. The oneloop contributions to the form factors F (1) i,σ are given by where the terms in the first line correspond to the contributions from the diagrams in Fig. 4, while the terms in the second line come from the renormalization procedure. For the computation of these EW corrections, we adopt the modified on-shell renormalization scheme defined in Ref. [43].
In the on-shell renormalization scheme, all the counterterms in the amplitude of e + e − → h Z are determined in terms of the one-particle irreducible (1PI) diagrams for one-and twopoint functions of Higgs bosons, gauge bosons and fermions by imposing a set of the renormalization conditions. Adding these counterterms, one can obtain the ultra-violet (UV) finite one-loop corrected vertices.
In the wide range of extended Higgs models, there are mixings among Higgs bosons, and the gauge dependence appears in the renormalization of these mixing angles. We apply the pinch technique to remove the gauge dependence in the renormalized vertex functions [32,41,43].
Apart from the UV divergences, there are infrared (IR) divergences when we calculate virtual photon loop contributions. In the calculation of individual photon loop con-tributions, we regularize them with a finite photon mass μ. The photon mass dependences in the one-loop calculation are exactly canceled by adding contributions of real photon emissions. The analytic expression of the real photon contribution with the soft-photon approximation is given by [61][62][63], with the photon energy cutoff E. The dependence of E vanishes in the inclusive cross section where one also includes the contribution of hard photon emissions [62]. The inclusive cross section still depends on ln m 2 e /s , and this logarithmic term potentially takes a large value. This dependence can also be eliminated by introducing the electron structure functions as discussed in Ref. [66]. However, the treatment of hard photon emission highly depends on the experimental setup. The hard photon changes the kinematics of the process, and these effects would be eliminated by applying appropriate experimental cuts. In addition, if one considers the scenario with κ Z 1, these effects in extended Higgs models are almost the same as in the SM. Therefore, we do not consider the electromagnetic effects when we focus on the difference between the predictions in the extended models and those in the SM.
In the one-loop calculation, we choose the fine structure constant α em , the Fermi constant G F and the Z boson mass m Z as the input EW parameters. In addition to these EW parameters, we also use the shift of the fine structure constant α em , the strong coupling constant α s and the masses of the fermions and the discovered Higgs boson as the input parameters. The values of these SM input parameters are given in Appendix A. We use the input parameters given in Eqs. (2.10), (2.26) and (2.37) in the HSM, the 2HDMs and the IDM, respectively.

Renormalized vertices
In the following calculation, the Z ff , h Z Z, h Zγ and h ff vertices are relevant, where h Zγ vertex is one-loop induced. Each of these vertices can be decomposed into several form factors depending on their Lorentz structure.
The renormalized Z ff vertex can be decomposed in the massless limit of external fermions as where p f ( pf ) is the incoming four-momentum of the fermion (anti-fermion), and p Z is the outgoing fourmomentum of the Z boson. We can further decompose these vertices into the tree-level and one-loop level contributions The tree-level contribution is given by (3.41) In the massless limit of external fermions, expressions of these vertices in the HSM, the 2HDMs and the IDM are the same as those in the SM. Analytic expressions for the 1PI diagrams and counterterms are presented in Appendix B in Ref. [55].
The renormalized h Z Z vertex can be decomposed as where p 1 and p 2 are incoming four-momenta of the Z bosons, and p h is the outgoing four-momentum of the SM-like Higgs boson. We can further decompose these vertices into the treelevel and one-loop level contributions The tree-level contribution is given by The form factor 3 h Z Z is non-zero only when the SM-like Higgs boson is a CP-mixed state. Therefore this form factor vanishes in the model with CP conservation in the Higgs sector. Analytic expressions for the 1PI diagrams and counterterms are presented in Ref. [34] for the HSM, Ref. [39] for the 2HDMs and Ref. [51] for the IDM.
Similarly, we define the loop-induced h Zγ vertex as The form factor 3 h Zγ vanishes in the model with CP conservation in the Higgs sector. Analytic expressions for the 1PI diagrams and counterterms are presented in Refs. [55,125].
The renormalized h ff vertex can be decomposed as where p f ( pf ) is the incoming four-momentum of the fermion (anti-fermion), and p h is the outgoing fourmomentum of the SM-like Higgs boson. Analytic expressions for the renormalized vertices are presented in Ref. [34] for the HSM, Ref. [39] for the 2HDMs and Ref. [51] for the IDM.

Expression of form factors including one-loop corrections
We list the one-loop contributions to the form factors in terms of the renormalized quantities. The one-loop propagator corrections appear in the sum in Eq. (3.37) as the term with the renormalized self-energies T V V (s) of the neutral vector bosons. The renormalized Zeē corrections appear as The renormalized h Z V (V = Z , γ ) corrections appear as The renormalized heē corrections appear as The 1/t and 1/u terms originated from the fermion propagator are canceled by the vertex corrections [63]. However, the renormalized vertices depend on t and u, and they cause non-trivial cos θ dependence.
There are the five W boson mediated and one Z boson mediated box diagrams in the massless limit of the electron. The amplitudes of the W boson mediated diagrams can be written as The amplitude of the Z boson mediated diagram has a different structure from others, and it can be written as with F 6 1 (s, t) = F 6 (s, t), (3.60) The expressions of C k and F k i (s, t) are given in Appendix B. We define B i,σ (s, t) by (3.63) Finally, the form factors at one-loop level are given by

Numerical results
In this section, we begin with an analysis on the behavior of the NLO weak corrections to the cross section of the e + e − → h Z process in the SM. We then analyze the deviations from the SM values at NLO in the HSM, the 2HDM and the IDM. We evaluate the form factors F (1) i,σ by using the H-COUP program [52,53], where G F is taken as an input. In order to compare our results in the SM with the previous works [63,[126][127][128], we extend the H-COUP program and take m W as an input instead of G F . With this extension, we have confirmed that our results are in agreement with the previous results. In the following, we show the results obtained in the scheme where α em (0), G F and m Z are input parameters.
In order to study the theoretical behavior of the one-loop corrections, we parametrize the differential cross section as dσ = dσ LO (1 + δ weak + δ em ), (4.1) where δ weak and δ em denote the relative weak and electromagnetic corrections, respectively. As we will see below, the NP effects mainly come from 1 h Z Z vertex, and that appear independently of σ and λ. Therefore, we show the results of the unpolarized cross section where the polarization of the Z boson is also summed. In order to analyze the NP effects in each renormalized quantity, we introduce EW X , with EW X = σ X /σ LO . We evaluate σ X by substituting σ λ (s, t) in Eq. (3.28), where F X i,σ is defined in Eq. (3.37). We also evaluate the ratios of the total cross sections to exhibit the deviations from the predictions in the SM, (4.4)

Standard model
In the left panel of Fig. 5, we show the weak one-loop corrections to the helicity-dependent cross sections as a function of the CM energy. The weak corrections to the cross section for the right-handed electron are positive, and they increase the cross section by about 10%. On the other hand, those for the left-handed electron are negative, and the size of these corrections strongly depends on the CM energy. The reason for this difference comes from negative contributions from the box diagrams. Among the six box diagrams, the five W boson mediated diagrams only contribute to the helicity amplitudes for the left-handed electron. They give negative contributions to the helicity-dependent cross sections for the left-handed electron. In addition, their effects become relevant at higher energies and give large negative corrections. The peak around √ s 350 GeV corresponds to the threshold at 2m t in the top-loop contributions. We note that the NNLO electroweak-QCD corrections have been estimated in Refs. [127,128], and the magnitude is about a percent.
In the right panel of Fig. 5, we show the weak one-loop corrections to the differential cross sections as a function of the CM energy. From Eq. (3.63), we can see that only the heē vertex and box corrections cause different cos θ dependence from those at LO. At √ s = 250 GeV, this effect is not so large, and the angular distribution of the Z boson is almost determined by the d 1 σ,λ (θ ) functions. At higher energies, the angular distribution of the Z boson is significantly modified through the box contributions [63]. However, the size of the cross sections decreases in such a higher energy region.

Higgs singlet model
First, we consider the Z 2 symmetric scenario in the HSM where c α = 1 and μ S = 0. There remain three input parameters, m H , λ S and λ S . In the following analysis, we impose perturbative unitarity, vacuum stability, avoiding wrong vacua, and the constraints of the EW S and T parameters. In order to analyze the theoretical behavior, we here do not impose the constraints from the direct searches of the additional Higgs boson and the Higgs coupling measurements. In addition, we impose M 2 > 0 as in the case of the IDM and the 2HDMs. (4.5) As the NP contributions, there are two H propagated diagrams in 1PI hh ( p 2 ). One of them is proportional to λ hh H H , while the other is proportional to λ 2 h H H . However, the former does not contribute to δ Z h because the loop function A(m H ) does not depend on the external momentum p 2 . Therefore, only the latter contributes to the helicity amplitude as the NP effect, with λ h H H = −λ S v. We note that this difference does not directly depend on λ S , but it indirectly determines the possible size of the NP effects through the perturbative unitarity and vacuum stability bounds. The magnitude of EW h Z Z becomes larger when the mass of the extra Higgs boson is taken to be larger up to around 700 GeV. This peak corresponds to the point where the minimum value of M 2 = m 2 H −λ S v 2 changes from zero to nonzero due to the perturbative unitarity bound. At this point, λ S takes the maximal value, and it triggers a sizable effect. 4 In the case of larger values of m H , the magnitude of EW h Z Z monotonically decreases because perturbative unitarity constrains the size of λ S . In such a large mass region, m H is approximately equal to M, and the additional Higgs boson almost decouples following the decoupling theorem [131].
We can also see the relatively large NP effects when the mass of the additional Higgs boson is below 200 GeV. In this region, λ S takes a negative value satisfying the vacuum stability bound thanks to the sizable λ S . While the sign of λ S is flipped, only |λ S | 2 appears in Eq. (4.6). Therefore, EW h Z Z is negative independently of the sign of λ S . In the right panel of Fig. 6, we show the predictions of R h Z in Eq. (4.4) as a function of the mass of the additional Higgs boson in the Z 2 symmetric HSM at √ s = 250 GeV. We here take λ S to 0.1, 1 and 2 and scan λ S for |λ S | < 4π . In the Z 2 symmetric HSM, the Z 2 symmetry prohibits the mixing of the CP-even states, and R h Z = 0 at LO. The behavior of R h Z is only determined by δ Z h , and almost the same as that of We here mention the scenario where μ S softly breaks the Z 2 symmetry. We note that R h Z does not directly depend on μ S , and the behavior of R h Z is the same as that in the Z 2 symmetric HSM. However, μ S indirectly affects the possible size of the NP effects through the conditions for avoiding wrong vacua. For example, the region where m H 300 GeV is excluded if μ S = 100 GeV.
Next, we consider the scenario with the mixing of the CP-even states. In the left panel of Fig. 7, we show EW X as a function of the mass of the additional Higgs boson with c α = 0.99 and μ S = 0 at √ s = 250 GeV. We here take λ S = 1 and scan λ S for |λ S | < 4π . We note that not only the renormalized h Z Z vertex but also the other renormalized quantities differ from the SM values unlike the case in the Z 2 symmetric HSM. However, the magnitude of EW h Z Z is larger than that of the others. For m H 700 GeV, EW h Z Z takes a negative value, while it takes a positive value for m H > 700 GeV. In order to realize the finite mixing of the CP-even states with a large mass of the additional Higgs boson, the Higgs quartic couplings should take large values, and it triggers a so-called non-decoupling effect. In the right panel of Fig. 7, we show the predictions of R h Z as a function of the mass of the additional Higgs boson with c α = 0.99 and μ S = 0 at √ s = 250 GeV. We here take λ S to 0.1, 1 and 2 and scan λ S for |λ S | < 4π . If there is the mixing of the CP-even states, the LO cross section decreases from its SM value. When c α = 0. 99

Inert doublet model
In the IDM, the Z 2 symmetry prohibits the mixing of the CPeven states, and R h Z = 0 at LO, similarly to the case in the Z 2 symmetric HSM. In the following analysis, we assume that the additional Higgs bosons are degenerate in their mass, m ≡ m H = m H ± = m A . There remain two input parameters, λ 2 and M 2 . We impose perturbative unitarity, vacuum stability, avoiding wrong vacua, and the constraints of the EW S and T parameters. In order to analyze the theoretical behavior, we here do not impose the constraints from the direct searches of the additional Higgs boson, the Higgs coupling measurements and the dark matter experiments.
In the left panel of Fig. 8, we show EW X as a function of the mass of the additional Higgs bosons in the IDM at √ s = 250 GeV. We here take λ 2 = 1 and scan M 2 for 0 ≤ M 2 ≤ (3 TeV) 2 . We note that not only the renormalized h Z Z vertex but also the other renormalized quantities differ from the SM values unlike in the Z 2 symmetric HSM. This is because the additional Higgs bosons are charged under the SM gauge group. However, the magnitude of EW h Z Z is larger than that of the others in most cases. In the IDM, there are two main contributions to EW h Z Z . The first one is λ 2 h terms originated from δ Z h , similarly to the case in the Z 2 symmetric HSM. In addition, there are 1PI diagram contributions proportional to  h Z Z , and its behavior is almost the same as that in the Z 2 symmetric HSM. The maximal deviation in the IDM is larger than that in the Z 2 symmetric HSM because we have more than one additional Higgs boson running in the loop in the IDM. It monotonically decreases when m > 600 GeV following the decoupling theorem.
We can also see the sizable NP effects when m is below 200 GeV, similarly to the case in the HSM. In addition,

Two Higgs doublet model
First, we consider the alignment limit with s β−α = 1, where R h Z = 0 at LO. We also assume that the additional Higgs bosons are degenerate in mass, m ≡ m H = m H ± = m A . There remain two input parameters, tan β and M 2 . We impose perturbative unitarity, vacuum stability, avoiding wrong vacua, and the constraints of the EW S and T parameters. In order to analyze the theoretical behavior, we here do not impose the constraints from the direct searches of the additional Higgs boson, the Higgs coupling measurements and the flavor measurements.
We analyze all the four types of 2HDMs, and it turns out that predictions for EW X and R h Z are almost the same. This is because differences among the four types of 2HDMs appear through the down-type quark and lepton Yukawa interactions with the SM-like Higgs boson. As we will see later, the magnitude of EW h Z Z is larger than that of the others in the 2HDMs, similarly to the case in the HSM and the IDM. In EW h Z Z , the top-quark contributions dominate fermionic contributions, and there is no sizable difference among the four types of 2HDMs. Therefore, we show the predictions for R h Z in the Type-I 2HDM as a representative in the following.
In the left panel of Fig. 9, we show EW X as a function of the mass of the additional Higgs bosons in the Type-I 2HDM with s β−α = 1 at √ s = 250 GeV. We here take tan β = 1.5 and scan M 2 for 0 ≤ M 2 ≤ (3 TeV) 2 . We note that not only the renormalized h Z Z vertex but also the other renormalized quantities differ from the SM values because the additional Higgs bosons interact with the gauge bosons, the quarks and the leptons. Qualitative behaviors of EW X are almost the same as those in the Z 2 symmetric HSM and the IDM except for m ≤ 200 GeV. The magnitude of EW h Z Z is larger than that of the others in most of the parameter space. It monotonically decreases when m > 450 GeV following the decoupling theorem.
There is no sizable negative correction below 200 GeV unlike in the Z 2 symmetric HSM and the IDM. This is because the Higgs quartic couplings in the 2HDMs are more constrained by vacuum stability than λ S in the HSM and λ 2 in the IDM. If m is lighter than about 150 GeV, EW X can be positive due to the 1PI diagram contributions, similarly to the case in the IDM.
In the right panel of Fig. 9, we show the predictions of R h Z as a function of the mass of the additional Higgs bosons in the Type-I 2HDM at √ s = 250 GeV. We here take tan β to 1.5, 3 and 5 and scan M 2 for 0 ≤ M 2 ≤ (3 TeV) 2 . The behavior of R h Z is almost the same as that of EW h Z Z . The possible magnitude of R h Z decreases as tan β becomes large due to the perturbative unitarity and vacuum stability bounds.
Next, we consider the scenario with the mixing of the CP-even states. In the top (bottom) left panel of Fig. 10, we show EW X as a function of the mass of the additional Higgs bosons in the Type-I 2HDM with s β−α = 0.99 and c β−α < 0 (c β−α > 0) at √ s = 250 GeV. We here take tan β = 1.5 and scan M 2 for 0 ≤ M 2 ≤ (3 TeV) 2 . The magnitude of EW h Z Z is larger than that of the others independently of the sign of c β−α . In addition, EW h Z Z takes a negative value except for the m H 300 GeV with c β−α < 0 unlike in the HSM. We can see the non-decoupling effect in a large mass region of the additional Higgs bosons because they cannot decouple while Finally, we mention the corrections to the angular distribution of the Z boson. As we have mentioned, the heē vertex and box contributions cause non-trivial cos θ dependence. In the limit of the massless electron, only mixing of CP-even states modifies the heē vertex and box contributions. However, as we can see from Figs. 7 and 10, these effects are rather small at √ s = 250 GeV. Therefore, the angular distribution of Z bosons is almost the same as the predictions in the SM.

Correlation in the cross section times decay branching ratios
In this subsection, we analyze the correlation in the cross section times decay branching ratios of the SM-like Higgs boson in the HSM, the IDM and the four types of 2HDMs.
At future collider experiments such as the ILC, the cross section of e + e − → h Z can be measured without depending on the decay of the SM-like Higgs boson by utilizing the recoil mass technique [26,132]. This makes it possible to measure the decay branching ratio of the SM-like Higgs boson independently of the cross section. However, the cross section times decay branching ratios of the SM-like Higgs boson can be measured more precisely. In Table 3, we summarize the expected accuracy of the cross section times decay branching ratios of the SM-like Higgs boson at the ILC at √ s = 250 GeV with 250 fb −1 for (P e , Pē) = (−0.8, +0.3). The values in Table 3 are taken from Table VI in [26].
In the following, we analyze the predictions for σ (e + e − → h Z) × BR(h → XY ) at one-loop level, where X and Y denote decay products of the SM-like Higgs boson. In order to discuss deviations from predictions in the SM, we evaluate the ratio of the total cross section times the decay branching ratios where we assume the beam polarization (P e , Pē) = (−0.8, +0.3).
In the evaluation of the decay branching ratios of the SMlike Higgs boson with the one-loop EW and QCD corrections, we use the H-COUP program [52,53]. Although the magnitude of R h Z XY depends on the treatment of the QED corrections, we do not consider these corrections and discuss the pattern of the deviations in the correlation of σ (e + e − → h Z) × BR(h → XY ). The QED corrections in the cross section universally change the magnitude of σ (e + e − → h Z) × BR(h → XY ). Therefore, the pattern of the deviations is not changed even if we include the QED corrections following a realistic experimental setup.
We scan the input parameters in each model in the following way. In the HSM, there are five input parameters as given in Eq. We here take μ s = 0 and λ S = 0.1 for simplicity.
In the 2HDMs, we have six input parameters given in Eq. (2.26). We assume that the additional Higgs bosons are degenerate in mass as in the previous subsection. In this scenario, the constraint of the EW T parameter is satisfied independently of the type of the 2HDMs. The degenerate mass m (= m H ± = m H = m A ) is scanned as 400 GeV ≤ m < 2000 GeV for the Type-I and X 2HDMs, (4.10) 800 GeV ≤ m < 2000 GeV for the Type-II and Y 2HDMs. (4.11) The lower bound of m in the Type-I and Type-X 2HDMs comes from the direct search for A → ττ at the LHC [25]. In the Type-I 2HDM, the parameter regions with tan β 2 are not excluded. However, we take the above parameter regions for simplicity. In the Type-II and Type-Y 2HDMs, the lower bound comes from the flavor experiments, especially from B s → X s γ [117]. In addition, we scan the other parameters as 0.98 ≤ s β−α < 1, 2 ≤ tan β < 10, (4.12) The lower bound of tan β comes from the flavor experiments. We analyze both the positive and negative signs of c β−α .
In the IDM, we have five input parameters given in Eq. (2.37). We take m H = 63 GeV which is favored by dark matter constraints. In order to satisfy the constraint from the EW T parameter, we assume that H ± and A are degenerate in mass. The degenerate mass m H ± (= m A ), M 2 and λ 2 are scanned as 100 GeV ≤ m H ± < 1000 GeV, (4.13) 0 ≤ M 2 < (m + 250 GeV) 2 , (4.14) 0 < λ 2 < 4π. (4.15) Over the above parameter spaces, we impose the constraints discussed in Sect. 2 such as perturbative unitarity, vacuum stability, avoiding wrong vacua, and the constraints from the EW S and T parameters. In addition, we take into account the current data of the signal strengths of the discovered Higgs boson at the LHC. We evaluate the decay rates of the SM-like Higgs boson with higher-order corrections by using the H-COUP program [52,53]. We define the scaling factor κ XY = NP h→XY / SM h→XY at the one-loop level, and remove the parameter points, where κ XY deviates from the observed data at 95% CL. In Table 4, we summarize the current measurements of κ XY factors at 1σ accuracy. The values in Table 4 are taken from Table XI in Ref. [3]. We assume that there is no decay mode where the SM-like Higgs boson decays into additional Higgs bosons.
In the Type-II, X and Y 2HDMs, we have parameter points where Yukawa coupling constants take the negative sign with a large value of tan β and c β−α > 0. These parameter points show the sizable deviation both in the Higgs branching ratio and cross section. However, we simply omit such particular The order of loop expansion of R Zh R XY is O(h 2 ), and it is sub-leading. Therefore, the qualitative behavior of R h Z XY can be understood by independently analyzing R Zh and R XY . The behavior of R XY has been studied in Ref. [55] by using the H-COUP program [52,53]. For later convenience, we briefly summarize the behavior of R XY in the HSM, the IDM and the four types of 2HDMs. First, in the HSM, the decay branching ratios of h are almost the same as those in the SM predictions, because the partial decay widths are universally suppressed by the radiative corrections and the mixing of the CP-even states. In our study, both R τ τ and R bb are at most 0.5%. The same argument has been claimed for the IDM in Ref. [55]. However, we find that the parameter regions where both R bb and R τ τ take a few percent deviations. This difference comes from the large value of λ 2 . In Ref. [55], λ 2 has been fixed to 0.1. However, the magnitude of λ 2 indirectly controls the possible size of other Higgs quartic couplings especially through the vacuum stability bound given in Eq. (2.31). We have obtained the almost same results as those in Ref. [55] when we impose λ 2 ≤ 0.1.
In the 2HDMs, the predictions to R τ τ and R bb spread out into different directions according to the type of the Yukawa interactions and the sign of c β−α . The possible magnitudes of the deviations in the Type-II, X and Y 2HDMs are rather large compared to the Type-I 2HDM, the HSM and the IDM. They can reach several tens of percent, and especially R τ τ reaches a hundred percent in the Type-X 2HDM. In Fig. 11, we show the correlations between R h Z τ τ and R h Z bb in the HSM, the IDM and the four types of 2HDMs. We take the color codes where orange, grey, red, blue, green and purple correspond to the HSM, the IDM and the Type-I, II, X, Y 2HDMs, respectively. The lighter color corresponds to the lighter mass scale of the additional Higgs bosons, m ≥ 400, 800, 1200 and 1600 GeV in order. In the left (right) panel, we show the results with c β−α < 0 (c β−α > 0). The results in the HSM and the IDM are the same both in the left and right panels.
As discussed in Eq. (4.16), R h Z XY can be rewritten as the sum of R Zh and R XY , and R Zh takes a negative value in most cases. In the Type-II, X and Y 2HDMs, the typical size of R XY is larger than R Zh . Therefore, the pattern of the deviation is mainly determined by R XY , and it is consistent with previous analysis in Ref. [55]. In these models, the possible sizes of the deviations are large enough to be detected at the ILC if m is about 1 TeV or less.
In the HSM, the Type-I 2HDM and the IDM, we can see sizable deviations both in R h Z τ τ and R h Z bb , and they reach about 10%. Both R h Z τ τ and R h Z bb take larger values than those in R τ τ and R bb . This is because the typical size of R Zh is larger than R XY in these models, and R XY also takes a negative value in the HSM, the Type-I 2HDM with c β−α < 0 and the IDM. In the Type-I 2HDM with c β−α > 0, R XY takes a positive value. However, the typical size of R XY is smaller than R h Z , and R Zh XY takes a negative value in most of the parameter regions. In the Type-I 2HDM, R Zh XY quickly decouples. This is because a non-zero c β−α realizes the maximal deviation in R h Z , especially at LO. If m is large, the possible value of c β−α is constrained mainly by perturbative unitarity. On the contrary, in the other types of 2HDMs with c β−α > 0, the decoupling behavior is not clearly seen. This is because taking the inner parameter tan β large keeps the magnitude of the deviation to be large even in the case of large m . On the other hand, constraints from the Higgs signal strength in the HSM are weaker than those in the 2HDMs, and the deviation in c α realizes the sizable R Zh XY even if m is larger than 1 TeV. In the IDM, the decoupling limit cannot be applied because we fix m H = 63 GeV. Therefore, we have a sizable deviation although there is no mixing between the CP-even states.
In Fig. 12, we show the correlations between R h Z τ τ and R h Z cc in the HSM, the IDM and the four types of 2HDMs. The color codes and gradations are the same as those in Fig. 11. In the left (right) panel, we show the results with c β−α < 0 (c β−α > 0). The results in the HSM and the IDM are the same both in the left and right panels.
Qualitative behavior of the deviations in each model is the same as in Fig. 11. In the Type-II, X and Y 2HDMs, the typical size of R XY is larger than R Zh , and the pattern of the deviation is mainly determined by R XY . On the other hand, we also have sizable deviations in the HSM, the Type-I 2HDM and the IDM, and they reach about 10%. In the Type-I 2HDM with c β−α > 0, R XY takes a positive value. However, the typical size of R XY is smaller than R h Z , and R Zh XY takes a negative value in most of the parameter regions.
If In order to discriminate the Type-I 2HDM from the HSM and the IDM, we can use the correlation between R h Z τ τ and R h Z W W . In Fig. 13, we show the correlations between R h Z τ τ and R h Z W W in the HSM, the IDM and the four types of 2HDMs. The color codes and gradations are the same as those in Fig. 11. In the left (right) panel, we show the results with c β−α < 0 (c β−α > 0). The results in the HSM and the IDM are the same both in the left and right panels. Especially in the case of c β−α < 0, R h Z W W takes a positive value in the Type-I 2HDM, while it takes a negative value in the HSM and the IDM. Even in the case of c β−α > 0, there is a stronger correlation between R h Z τ τ and R h Z W W in the HSM and the IDM than those in the Type-I 2HDM.
Finally, we discuss the discrimination between the HSM and the IDM. The deviation R h Z γ γ might be useful to discriminate these models because it is mainly affected by the contribution of the charged Higgs bosons. The behaviors of R h Z τ τ and R h Z γ γ show a different correlation between the HSM and the IDM. However, the possible size of R h Z γ γ is at most 20%, and it is rather challenging to discriminate them with 95% CL at the ILC. The large uncertainty in R h Z γ γ at the ILC mainly comes from the low statistics, and this would be improved by performing the combined study with the measurements at the ILC and the HL-LHC.

Conclusion
We have calculated the cross section for e + e − → h Z with arbitrary sets of electron and Z boson polarization at the full next-to-leading order in the HSM, the IDM and the four types of 2HDMs. We have systematically performed complete one-loop calculations to the helicity amplitudes in the on-shell renormalization scheme and present the full analytic results as well as numerical evaluations. The deviation R h Z in the total cross section from its SM prediction has been comprehensively analyzed, and the differences among these models have been discussed in detail. We have found that new physics effects appearing in the renormalized h Z Z vertex almost govern the behavior of R h Z . We have also shown that the predictions for the deviations in the total cross sec-tion of e + e − → h Z times the branching ratios of h → XY . It has been found that we can discriminate the four types of 2HDMs by analyzing the correlation between R h Z bb and R h Z τ τ and those between R h Z cc and R h Z τ τ . Furthermore, the Type-I 2HDM might be specified from the HSM and the IDM by measuring the deviation in R h Z W W . These signatures can be tested by precision measurements at future Higgs factories such as the ILC. On the other hand, the discrimination between the HSM and the IDM is rather challenging only by the measurement at the ILC. However, this problem might be solved by taking into account the deviation in h → γ γ signals at the LHC and the HL-LHC.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: There are no external data associated with the manuscript.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Input parameters
We work in the scheme where α em (0), G F and m Z are input parameters following the H-COUP program. In Table 5, we list the SM input parameters. The values of input parameters are taken from Ref. [94]. Other parameters can be evaluated in terms of the above inputs by using tree-level relations.

Appendix B: Analytic formulae for the box diagrams
We here give the analytic expressions for the box diagram contributions denoted by F k i (s, t) in Eqs. (3.55) and (3.59) in terms of the Passarino-Veltman functions defined in Ref. [133]. We calculate the box diagrams in the 't Hooft-Feynman gauge (Fig. 14).
The analytic expressions are given as follows