High Scale Validity of the DFSZ Axion Model with Precision

Satsuki Oda, 2 Yutaro Shoji, and Dai-suke Takahashi 2 Okinawa Institute of Science and Technology Graduate University (OIST), Onna, Okinawa 904-0495, Japan Research Institute, Meio University, Nago, Okinawa 905-8585, Japan Kobayashi-Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Nagoya, Aichi 464-8602, Japan Abstract With the assumption of classical scale invariance at the Planck scale, the DFSZ axion model can generate the Higgs mass terms of the appropriate size through technically natural parameters and may be valid up to the Planck scale. We discuss the high scale validity of the Higgs sector, namely the absence of Landau poles and the vacuum stability. The Higgs sector is identical to that of the type-II two Higgs doublet model with a limited number of the Higgs quartic couplings. We utilize the state-of-the-art method to calculate vacuum decay rates and find that they are enhanced at most by 1010 compared with the tree level evaluation. We also discuss the constraints from flavor observables, perturbative unitarity, oblique parameters and collider searches. We find that the high scale validity tightly constrains the parameter region, but there is still a chance to observe at most 10% deviation of the 125 GeV Higgs couplings to the fermions.


I. INTRODUCTION
An invisible axion [1][2][3][4][5][6][7][8][9][10][11] is one of the plausible solutions to the strong CP problem and is also an excellent dark matter candidate. We focus on the DFSZ axion model [10,11], where the standard model (SM) is extended with a SM singlet complex scalar and an additional Higgs doublet. Since the Higgs doublets have non-zero Peccei-Quinn (PQ) charges, the Higgs couplings are tightly restricted by the PQ symmetry. For example, dangerous flavor changing neutral currents (FCNC) are forbidden at the tree level and the CP is not broken spontaneously in the scalar sector.
In this paper, we discuss the possibility that the DFSZ axion model remains valid up to the Planck scale. In such a scenario, one of the disadvantages is that we need to give up a complete explanation of the hierarchy between the Planck scale and the electroweak (EW) scale. However, if there is a mechanism that realizes classical scale invariance at the Planck scale, the hierarchy problem may be solved without introducing supersymmetry or compositeness [12][13][14]. Since the scale invariance is violated at the quantum level, the PQ breaking scale can appear through the dimensional transmutation. If the PQ sector and the Higgs sector are connected by (technically natural) tiny couplings, the PQ breaking can also generate the Higgs mass terms without causing a hierarchy problem [15,16]. Since the PQ breaking sector decouples from the Higgs sector due to the tiny couplings, the model is well approximated by the type-II two Higgs doublet model (THDM) with a restricted number of coupling constants. Importantly, the additional Higgs bosons should be around the EW scale in this scenario since there is no technically natural parameter that accommodates a hierarchy among the Higgs boson masses.
Another disadvantage is that the model does not explain the neutrino masses, the baryon asymmetry of the Universe, or inflation. However, they can be explained without affecting the Higgs sector. For example, one may consider the see-saw mechanism with right handed neutrinos having a few orders of magnitude smaller masses than the PQ breaking scale [17][18][19]. It can explain the neutrino masses and also the baryon asymmetry of the Universe with tan β 4 [20]. However, it does not cause the hierarchy problem thanks to the tiny Yukawa couplings of the right handed neutrinos. As for inflation, one may attach an inflation sector to the model and assume tiny couplings between the inflation sector and the Higgs sector, which is, at least, technically natural.
For the model to be valid up to the Planck scale, Landau poles should not appear during the renormalization group (RG) evolution and the lifetime of the EW vacuum should be long enough. We refer to these two conditions as the high scale validity. Similar discussions can be found in the context of THDMs [21][22][23][24][25][26][27][28][29][30][31][32][33]. As we will see and as found in the previous studies, these conditions are complementary and become very restrictive if combined. Thus, the model becomes more predictive and it is important to determine the allowed parameter space precisely.
The lifetime of the EW vacuum is estimated by the bubble nucleation rate [34,35], which has a form of where B is the Euclidean action of the so-called the bounce, and A represents quantum corrections to B having mass dimension four. In many papers, A is assumed to lie around the typical scale of the problem, but it has been pointed out [36] that such an estimation leads to theoretical uncertainty of e −B×O(10%) in the nucleation rate. As we will see later, it can become comparable with the uncertainties coming from those of the top mass and the strong coupling. Thus, it is important to calculate both of A and B to get a precise vacuum decay rate.
Recently, the one-loop calculation of A for the SM has been completed [37][38][39] and the analytic expression for A at the one-loop level has been derived [37,39] for an approximately scale invariant theory. Since they are applicable to the case where the bounce is composed of a single field, we extend them to a multi-field case in this paper. Differently from the single-field case, there can be more than one unstable directions and there can appear an additional zero mode due to a global symmetry breaking. In addition, the electromagnetic U (1) symmetry can also be broken spontaneously.
Before the analysis of the high scale validity, we impose the constraints from flavor observables, perturbative unitarity, oblique parameters and collider searches. For the constraints from flavor observables, we obtain the 95% exclusion limit in Appendix A using the recent experimental values.
We determine the allowed parameter space by utilizing the Monte Carlo method. We show how much the high scale validity narrows down the parameter space and discuss the implications on the Higgs couplings and the Higgs mass splittings.
This paper is organized as follows. In Section II, we briefly explain the DFSZ axion model. Section III is devoted to the details of the analysis on the bubble nucleation rate for the multi-field case. Then, in Section IV, we discuss the low energy constraints. In Section V, we execute numerical analysis and discuss the consequence of the high scale validity.
Finally, we summarize in Section VI.

II. DFSZ AXION MODEL
In this section, we briefly review the DFSZ axion model [10,11]. The scalar sector consists of two Higgs doublets, H 2 and H 1 , and a SM singlet complex scalar, Φ. We choose the PQ charges of H 1 , H 2 and Φ to be x 1 , x 2 and (x 2 − x 1 )/2, respectively. Here, we assume x 1 = x 2 so that Φ has a non-zero PQ charge.
The general scalar potential is given by where v 2 Φ ,m 2 i 's, λ i 's,λ Φ andκ i 's are constants. We assume the classical scale invariance and setm 2 1 andm 2 2 to zero at the Planck scale. Then, the Higgs mass terms are assumed to be generated through the PQ symmetry breaking. In order to obtain the EW scale,κ i 's should be very small since Φ has to develop a huge vacuum expectation value (VEV) to avoid the constraints on the axion decay constant, 10 9 GeV f a 10 12 GeV [40,41]. Due to the smallness ofκ i 's, Φ decouples from the Higgs sector and the potential reduces to where PQ Charge Assignment Here, we took m 2 3 to be real and positive by the redefinition of the phase of H 1 . Notice that PQ violating quartic couplings can be generated after the PQ symmetry breaking, but they are suppressed byκ i 's and hence are negligible.
With the PQ charge assignment shown in Table I, the Higgs doublets couple to the SM fermions as Here, iσ 2 is the completely anti-symmetric matrix and Q, L, U , D and E represent the left quark doublets, the left lepton doublets, the up-type quarks, the down-type quarks and the charged leptons in the SM, respectively. The model is thus regarded as the type-II THDM with a limited number of Higgs quartic couplings.
Let us define the mass eigenstates and the mixing angles. We expand the Higgs fields as where v i 's are the VEVs of the Higgs fields, tan β = v 2 /v 1 and Here, h is the 125 GeV Higgs boson, H is the additional CP-even Higgs boson, A is the CP-odd Higgs boson, H + is the charged Higgs boson, and G 0 and G + are the would-be Nambu-Goldstone bosons. The SM-like limit for h is given by β − α → π/2.

III. VACUUM DECAY RATE
Since the quantum corrections to the effective potential depend on the VEVs of the Higgs fields, the shape of the effective potential is non-trivial at large Higgs VEVs. When there is a deeper vacuum or the effective potential is unbounded from below, the EW vacuum is not absolutely stable and decays through quantum tunneling. Even in such a case, we can live in the meta-stable vacuum if it has a much longer lifetime than the age of the Universe.
In this section, we discuss the precise determination of vacuum decay rates for the DFSZ axion model.
Recently, the analytic formulas for the prefactor, A, at the one-loop level have been derived [37,39], which are applicable to the case where the theory is approximately scale invariant and the bounce consists of a single field. In the following, we extend their results to the case where the bounce consists of more than one fields.
Since the field value at the true vacuum is typically much larger than the EW scale 1 , the Higgs potential is approximately given by For the moment, we fix the renormalization scale and will discuss the running effect later.
The bounce is a solution to the Euclidean equations of motion given by with boundary conditions, 1 The another vacuum may be close to the EW vacuum, which happens when the low energy potential already has an instability and the RG running cures it above the EW scale. We will put an IR cut-off on the size of the bounce to avoid such a situation.
and their complex conjugates. Here, r is the radius from the center of the bubble and i = 1, 2 labels the components of the doublet. Without loss of generality, we parameterize the Higgs fields as where σ i 's are the Pauli matrices. Then, the potential is expressed as where λ φ (Θ, Ω) = 1 2 λ 1 cos 4 Ω + λ 2 sin 4 Ω + 2(λ 3 + λ 4 cos 2 Θ) sin 2 Ω cos 2 Ω , We put an ansatz that Θ and Ω do not depend on r since it apparently minimizes the action.
Then, the equations of motion reduce to with boundary conditions dφ dr From Eq. (18), we can see that a minimum of λ φ satisfies cos 2 Θ = 0 for λ 4 ≥ 0, and cos 2 Θ = 1 for λ 4 < 0. Then, from Eq. (21), we get the following solutions; (a) Ω = 0, Notice that (c) exists only when (λ 1 −λ)/(λ 2 −λ) > 0. For each case, λ φ is given by If λ φ < 0, the solution to Eqs. (22) and (23) is given by which gives with R being a free parameter that fixes the radius of the bounce. Notice that B is independent of R, which is due to the (approximate) classical scale invariance.
Since all the possible bounces contribute to the vacuum decay rate, the total vacuum decay rate is expressed as where λ φ is summed over its minima with λ φ < 0. Now, the problem is reduced to the single field case for each λ φ and we can use the results of [37][38][39].
Before an explicit calculation of a vacuum decay rate, let us discuss the convergence of the R integral. From the dimensional analysis and the renormalization scale independence of the vacuum decay rate, the R-dependence of the integrand can be determined as λ φ is the one-loop beta function for λ φ . Thus, if we integrate it over R ∈ (0, ∞), the integration does not converge. However, as discussed in [37][38][39], the result can be convergent if we include higher-loop corrections. Although it is very difficult to calculate them, their R-dependence is completely determined by the beta functions and we can sum them up by taking µ ∼ 1/R for each bounce with radius R.
From [39], the differential vacuum decay rate is expressed as where ln ln Here, [ln A (X) ] MS 's are defined in [39]. The degrees of freedom, n (X) i , and the couplings, κ i , y i and g 2 i , are summarized below for each case. For case (c), the symmetry breaking pattern depends on the sign of λ 4 . Thus, we divide it into two cases; (c.1): λ 4 < 0 and (c.2): λ 4 > 0.
case (a): The scalar contributions: The fermion contributions: The gauge boson contributions: case (b): The scalar contributions: The fermion contributions: The gauge boson contributions: case (c.1): The scalar contributions: The fermion contributions: .
The gauge boson contributions: case (c.2): The scalar contributions: The fermion contributions: 10 The gauge boson contributions: .
Here, g Y and g 2 are the gauge couplings for U (1) Y and SU (2) L , respectively. The group volume is V G = 2π 2 for cases (a), (b) and (c.1), and V G = 4π 3 for case (c.2).
Notice that we can determine whether a solution, φ, is a minimum or not from the sign of κ i − λ φ . Let χ i be a scalar orthogonal to φ. Then, its potential can be written as where δ breaks the rotational symmetry of (φ, χ i ). Since κ i can be read off from the mass term for χ i , we get Thus, for φ to be a minimum of the action, we need κ i − λ φ > 0 for all i.
For case (c.1), we have κ 3 = λ φ and thus there appears a zero mode, which is due to the spontaneous breaking of the PQ symmetry. Its treatment is discussed in Appendix E of [39] and we replace with V σ = 2π.
Next, let us discuss the range of the R-integral. So far, we have assumed R ∈ (0, ∞), but it is not appropriate for our setup. First, we need an IR cut-off because we have ignored the dimensionful couplings 2 . Second, we need a UV cut-off because we do not consider gravitational corrections. Thus, we set the integration region as where M Pl is the reduced Planck scale. We also impose the same limits on the field value of the bounce as when λ φ < 0. In addition, we exclude the region where any of | ln A|'s become larger than 80% of B since the perturbative expansions become unreliable. Such a region appears where λ φ is very close to zero. Since the integrand of the vacuum decay rate is positive definite, these limits always make the vacuum decay rate small. Thus, what we get with these limits is a lower bound on the vacuum decay rate.
The condition for the stability of the EW vacuum is then given by where H 0 67.66 (km/s)/Mpc [42] is the current Hubble constant. Next, we calculate the differential vacuum decay rate, dγ/d(ln R), for case (b). We take µ = 1/R. The result is shown with the solid line in the top right panel of Fig. 1. Integrating it over ln R, we get where the 1st, 2nd, 3rd and 4th errors are those from m t , α s , m h and µ, respectively 3 . We use the SM values and uncertainties given in Table III and α s = 0.1181 (11). We estimate the renormalization scale uncertainty by taking µ = 2/R and µ = 1/(2R). The vacuum decay rate is much larger than log 10 [H 4 0 × Gyr Gpc 3 ] −3 and thus the parameter set is excluded by the vacuum stability.
Let us see the difference between the "tree level" vacuum decay rate and our result. For the tree level vacuum decay rate, we adopt where the maximum value is searched in the same region as the integration region of the one-loop vacuum decay rate. In the top right panel of Fig. 1 with the dashed line. We get Thus, the 1-loop calculation enhances the vacuum decay rate by about 10 6 , which is comparable with the uncertainties from those of the top mass and the strong coupling constant.
We show each quantum contribution in the bottom panel of Fig. 1. The vertical black dashed line corresponds to the maximum of the differential vacuum decay rate. Around the In Fig. 2, we show the binned plot of the vacuum decay rates at the tree level and at the 1-loop level by using the data accumulated for Fig. 4. We observe that the enhancement of the vacuum decay rate is generic for γ H 4 0 and that it is enhanced at most by 10 10 . For γ H 4 0 , the vacuum decay rate can be either suppressed or enhanced.

IV. LOW ENERGY CONSTRAINTS
Before the discussion of the high scale validity, let us discuss the low energy constraints; flavor observables, perturbative unitarity, oblique parameters, and collider searches. In this section, we do not consider the constraints from the signal strengths of the 125 GeV Higgs boson since they are the outputs of our analysis.

A. Flavor Observables
The additional Higgs bosons contribute to flavor observables and the current strongest constraints for the type-II THDM come from the branching ratios of B → τ ν, B s → µµ and b → sγ, and the B s −B s mixing as discussed, for example, in [43][44][45]. We obtain the constraints following the analysis of [46] with the current experimental values. The details are given in Appendix A.
In Fig. 3 4 These parameters affect only BR(B 0 s → µ + µ − ). As we will see later, the high scale validity requires a small cos(β − α) and mass differences. Then, the result is not so much affected as discussed in [46]. 5 Since there are choices of input parameters and of the treatment of theoretical uncertainty, O(10%) difference of the constraints is acceptable.

B. Low Energy Perturbative Unitarity
For the study of the high scale validity, the perturbative unitarity is necessary because otherwise all the calculations, including the matching conditions to the MS couplings, become unreliable.
At the tree level, the Higgs quartic couplings are related to the Higgs masses and mixing as Then, we impose the condition of the s-wave unitarity, which is given by [47,48] at the tree level. Since we will use the same condition to detect Landau poles later, we refer to the perturbative unitarity with the tree level matching conditions as "the low energy (LE) perturbative unitarity".

C. Oblique Parameters
The oblique parameters, especially the S-parameter and the T -parameter, are affected by the additional Higgs doublet. We use the general formulas for multi-Higgs-doublet models [49,50] (for the THDM, see [51] 6 ) to calculate these parameters.
The current constraints are given by [52] S = 0.02 ± 0.07, with the assumption of U = 0. The correlation coefficient is ρ = 0.92. We adopt 95% exclusion limit on S and T , which is given by where S cent and T cent are the central values of S and T , respectively. If it is satisfied, we then check the vacuum stability, where we take µ = 1/R.
We reduce the number of free parameters by choosing three slices of parameter space; (i) m H + = 600 GeV, 1.8 < tan β < 25, (ii) m H + = 900 GeV, 0.8 < tan β < 33, which satisfy the flavor constraints of Fig. 3. Since the flavor constraints do not depend so much on the other Higgs masses or cos(β − α) in the region of interest, we do not further check the flavor constraints to reduce computational complexity. In addition, we assume sin(β − α) > 0 in this analysis.
For each slice, we generate random 2M data points that satisfy all of the other low energy constraints, namely, LE perturbative unitarity, oblique parameters and collider searches.
The scattering range covers all of the parameter space where the LE perturbative unitarity is satisfied. The details of data generation are in Appendix C.
In Fig. 4, we show the binned plots of the allowed data points. All the colored points satisfy the low energy constraints. In the upper panels, the large tan β region is excluded by the H → τ τ channel. For slice (i), the upper and the lower bounds on cos(β − α) are determined by the constraints on the H → V V and the H → 2h → 4b channels, respectively.
For slices (ii) and (iii), the upper and the lower bounds on cos(β − α) are mostly determined by the constraints on the LE perturbative unitarity and the oblique parameters, respectively.
As for the lower panels, the concave shape is due to the constraint on the oblique parameters and the horns have the ends due to the other constraints.
Next, the orange and the green points satisfy the HE perturbative unitarity. The allowed parameter space is reduced especially for slice (i), but the reduction is not so drastic.
Finally, the green points satisfy the vacuum stability condition. As we can see, the parameter space is reduced drastically. It is because of the complementarity of the HE perturbative unitarity and the vacuum stability. It can be understood from the one-loop beta functions of λ 1 and λ 2 , which are given by Since y t , y b , y τ and g 2 are UV free, β λ 1 and β λ 2 generically become positive at a high energy scale. To avoid Landau poles, the quartic couplings should be small enough. In addition, negative λ 1 or λ 2 are preferable since they delay the appearance of Landau poles. Thus, the potential easily becomes unstable and a large part of the parameter space is constrained by the vacuum stability.
A similar condition as the vacuum stability is the bounded-from-below condition, which is given by Here, we regard those couplings as the MS couplings at µ = m t and impose it only at low energy. Notice that the condition is obtained from the discussion of Section III and is equivalent to that in [65]. We expect that the combination of the HE perturbative unitarity and the bounded-from-below condition should give a similar result 7 , which we will see below.
In Fig. 5, we pick up the parameter space defined by the region surrounded by the white dashed lines in Fig. 4 and prepare additional 5M points satisfying all the low energy constraints for each region. The scattering region is taken so that it can cover all the green points. The distribution of the new data points is uniform in the space of tan β, Here, | cos(β − α)| max is the maximum value of cos(β − α) depending on tan β, which is shown in Fig. 4. The red points satisfy the bounded-from-below condition and the HE perturbative unitarity. The lighter and the darker green points correspond to the green points in Fig. 4 and are plotted over the red points. Thus, in the red region appearing in the figure, the potential is stable at low energy, but always becomes unstable at high energy. The darker green points satisfy both the vacuum stability and the bounded-from-below conditions. Thus, in the lighter green region, the potential is always unstable at low energy, but the instability is cured at high energy. Notice that the vacuum decay rates may be affected by the IR cut-off for the R integral in the lighter green region.
As we can see from the figure, the bounded-from-below condition has a similar effect as the vacuum stability condition, but the allowed regions do not overlap completely. In particular, a large part of the region with m H < m A is excluded by the vacuum stability, where λ 2 tends to become negative during the RG evolution. In addition, a negative cos(β − α) is more favored by the vacuum stability.
Let us discuss the implication on the Higgs couplings. At the tree level, the SM-value normalized couplings of the 125 GeV Higgs boson are given by where U, D, L, and V represent the up-type quarks, the down-type quarks, the leptons, and the gauge bosons, respectively. Since | cos(β − α)| 0.04 for all the slices, we have 0.9992 g hV V ≤ 1, which is not possible to be distinguished from unity even with HL-LHC plus 1 TeV ILC [66]. It also means that the model can not be valid up to the Planck scale if we observe larger deviations of g hV V couplings. On the other hand, the other couplings can deviate by more than 1% because of the second term of the above equations.
In Fig. 6, we cast the data points in Fig. 5 into the g hU U vs g hDD = g hLL plane. The colors are the same as in Fig. 5. As we can see, g hU U can be reduced by about 2% − 6% for each m H + , but can not be enhanced so much. On the other hand, g hDD and g hLL can deviate by about 5% − 10% for each m H + , and tend to be enhanced.
The current constraints on these couplings are given by [67] g hZZ = 1.10 ± 0.08, g hW W = 1.05 ± 0.08, with the assumption that there is no new particles in loops and decays. Thus, they have already started to touch the parameter space. Future measurements of the Higgs couplings by, for example, the combination of HL-LHC and ILC will reach the precision of a few percent level [66] and will possibly find deviations from the SM values.
Let us discuss the dependence on the IR cut-off and the UV cut-off, which are introduced in Eqs. (67) and (68). Since the beta functions for λ 1 and λ 2 become generically positive at high energy, a factor change of the UV cut-off rarely affects the vacuum decay rate. As for the IR cut-off, we have checked that Figs. 5 and 6 are not affected even if we use 1 TeV for the IR cut-off, instead of 10 TeV.
Finally, we comment on the effect of A, which we have calculated precisely. Although the vacuum decay rates are enhanced compared with the tree level ones around γ ∼ H 4 0 , we find that the effect is not large enough to change Fig. 6. It is because of the strong dependence of the vacuum decay rates on the Higgs quartic couplings. However, if we find the additional Higgs bosons in future, the vacuum decay rate can be determined precisely from the measurements of the mass differences and the couplings of the Higgs bosons, which will become very important for the scenario.

VI. SUMMARY
In this paper, we analyzed the high scale validity of the DFSZ axion model, namely the HE perturbative unitarity and the vacuum stability. The model has been widely studied since it can explain the strong CP problem and dark matter elegantly. Once we admit a mechanism that forces classical scale invariance at the Planck scale, the Higgs mass terms of the appropriate size can be generated through the technically natural parameters and may be valid up to the Planck scale. In addition, the model can be extended without affecting the Higgs sector to explain the neutrino masses, the baryon asymmetry of the Universe and inflation. Thus, the high scale behavior of the Higgs sector is worth discussing.
We utilized the state-of-the-art method to calculate the vacuum decay rate precisely. We extended the results of [37][38][39] to accommodate bounces that are composed of more than one fields. Then, we showed that A can enhance the vacuum decay rates at most by 10 10 , which can become comparable with the uncertainties from those of the top mass and the strong coupling constant.
We performed the parameter scan and found the parameter space that satisfies the constraints from flavor observables, LE/HE perturbative unitarity, oblique parameters, collider searches, and vacuum stability. Due to the complementarity of the HE perturbative unitarity and the vacuum stability, the allowed parameter space becomes very small. We observe that it still accommodates at most 10% enhancement of the hDD and hLL couplings, and at most 6% suppression of the hU U couplings. These are around the current experimental constraints and will be searched at future experiments such as HL-LHC and ILC. On the other hand, the deviation of the hV V couplings are found to be smaller than 0.1% and the scenario may be excluded if we observe large deviations of these couplings. In this appendix, we follow [46] and obtain the flavor constraints with the current experimental values.

CKM Matrix Elements
We first determine the CKM matrix elements by using observables that are insensitive to the additional Higgs bosons.

Flavor Constraints
For the theoretical evaluation of the flavor observables, we use the formulas given in [46].
In the calculation, we utilize RunDec [63,64] to get the running masses of quarks.
Let us clarify the statistical method that we adopt. For observable X that depends on known parameters {x i ± δx i } and model parameters {y i }, we define  where X exp ± δX exp is the experimental result, X({y i }) is the theoretical result for inputs {x i } and {y i }, and Then, we use the offset corrected χ 2 defined as and the 95% CL exclusion limit is given by ∆χ 2 ({y i }) 3.84. Here, the minimum value is searched over the parameter space of Fig. 3 and the SM limit. The known parameters summarized in Tables II and III, where we ignore uncertainties of O(0.1%) 8 .
The result is shown in Fig. 3.

Appendix B: Matching Conditions
This appendix is devoted to the explanation of the one-loop matching conditions for the dimensionless coupling constants. The matching scale is taken to be µ t = m t . For a detailed discussion of the renormalization scheme, see [74,75]. We assume that tan β and cos(β − α) are renormalized with the MS scheme. In the calculation of one-loop threshold corrections, we utilize the public codes of SARAH [58,59], FeynArts [60], FeynCalc [61, 62].

Gauge Couplings
We first evaluate the electric charge at µ t as where [52] [α SM ( It is then matched to the THDM electric charge as where Next, we calculate the MS masses of the gauge bosons as with V = W, Z. Here,Σ T V (p 2 ) is the self energy for the transverse mode with 1/ε being subtracted. Here, where D is the spacetime dimension and γ E is the Euler number. The superscript, OS, indicates the on-shell mass.
Using these, the Weinberg angle is calculated as Then, the MS gauge couplings are given by

Yukawa Couplings
The MS tau mass is obtained from whereΣ S τ (p 2 ),Σ L τ (p 2 ) andΣ R τ (p 2 ) are the scalar, the left-handed and the right-handed parts of the self energy with 1/ε being subtracted.
As for the MS masses of the top quark and the bottom quark, we include the four-loop QCD corrections by using RunDec [63,64]. Then, we add the non-QCD one-loop threshold corrections to the output of RunDec as Here, m RunDec f (m t ) is the output of RunDec and the subscript g 3 = 0 indicates that the strong coupling is switched off in the calculation.
Then, the MS Yukawa couplings are given by . (B16)

Higgs Quartic Couplings
To adjust the Higgs VEVs order by order in perturbative expansions, we extend the scalar potential with tadpole terms as where T h and T H are zero at the tree level.
The MS values of these couplings are chosen as is the tadpole contributions to the effective action with 1/ε being subtracted.
As for the scalars, the MS masses are given by with X = h, H, A, H + . Here,Σ X (p 2 ) is the self energy with 1/ε being subtracted.
The MS Higgs quartic couplings are then obtained as where all the quantities appearing in the right-hand side are the MS values; we suppressed the renormalization scale, µ t , for visibility.

Appendix C: Generation of Data Points
In our analysis, we need to generate data points that consist of (m H , m A , tan β, cos(β−α)) for a fixed m H + . We first take a random tan β, which is uniformly distributed in the ranges defined in Eqs. (88)-(90). The other variables are generated with the procedure described in this section. The generated data points are then filtered by the perturbative unitarity conditions and are passed to the next analysis. In this appendix, we answer the following questions: (i) what is the appropriate range for m H , m A and cos(β − α) that covers all the points allowed by the perturbative unitarity? (ii) how can we effectively generate data points that are allowed by the perturbative unitarity?
A naive answer to question (i) is that |m 2 H + − m 2 H,A | 8πv 2 and | cos(β − α)| ≤ 1, where v 246 GeV. However, they are too weak to be used for the parameter scan. As we can see from Fig. 4, the allowed mass differences are smaller than about 200 GeV. However, one realizes that √ 8πv 1.2 TeV. It also means that H can be as light as h and the mixing angle can become large, which is why we naively expect no constraint on the mixing angle.
However, the allowed | cos(β − α)| is smaller than about 0.02 and becomes much smaller in the large tan β regime. Thus, if we scattered the data points over this naive range, we could get only a very few points that satisfy the low energy constraints. That is why we have question (ii).
We find that the speed of the data generation is fast enough and 50% − 60% of the generated data points satisfy the perturbative unitarity conditions.

Appendix D: Couplings and Partial Decay Widths of Heavy Higgs Bosons
In this appendix, we summarize couplings and the partial decay widths of the Higgs bosons, which are used for the inputs of HiggsBounds [53][54][55][56][57]. We use the results of [76][77][78]. The tree level couplings of the neutral Higgs bosons are shown in Table V. They are normalized by the corresponding couplings of the SM Higgs boson having the same mass as that of the decaying particle.
In the following, we use the running mass for the quark mass; where X represents the decaying particle. The reference value m q (µ 0 ) is calculated with RunDec [63,64] with µ 0 = 500 GeV. We define the following variables;