High scale validity of the DFSZ axion model with precision

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 about 10% deviation of the 125 GeV Higgs couplings to the fermions.


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 JHEP03(2020)011 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.

JHEP03(2020)011
The one-loop calculation of A for the SM was first calculated in [37]. Since the treatment of the gauge zero mode had not been known at that time, the calculation was not complete. Recently, the correct treatment has been found [38] and the one-loop calculation for the SM has been completed [39][40][41]. In addition, the analytic expression for A at the one-loop level has become available [39,41] 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 2, we briefly explain the DFSZ axion model. Section 3 is devoted to the details of the analysis on the bubble nucleation rate for the multi-field case. Then, in section 4, we discuss the low energy constraints. In section 5, we execute numerical analysis and discuss the consequence of the high scale validity. Finally, we summarize in section 6.

DFSZ axion model
In this section, we briefly review the DFSZ axion model [10,11]. The scalar sector consists of two Higgs doublets, H 1 and H 2 , 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λ Φ is moderate so that the VEV of Φ is not affected by those of H 1 and H 2 .
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 JHEP03(2020)011 Table 1. Assignment of the PQ charge in the DFSZ axion model. decay constant, 10 9 GeV f a 10 12 GeV [42,43]. Due to the smallness ofκ i 's, Φ decouples from the Higgs sector and the potential reduces to

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 1, 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.

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.

Formulation
Recently, the analytic formulas for the prefactor, A, at the one-loop level have been derived [39,41], 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 PQ-breaking sector couples to the THDM sector very weakly, the vacuum decay rate can be calculated independently of the PQ-breaking sector, i.e. the decay path, the RG running or the calculation of A is not affected by the PQ-breaking sector. 1 Notice that even when the field value of H 1 or H 2 becomes much larger than the PQ-breaking scale, Φ is almost constant during the tunneling. This is because the typical size of a bounce, i.e.R 1/ |H 1 (0)| 2 + |H 2 (0)| 2 , is too small. Here, H i (0)'s are the field values at the center of the bounce. For example, let us assume that Φ obtains a negative mass squared, m 2 Φ < 0, during the tunneling. Then, the displacement of Φ is roughly estimated as v Φ (e √ |m 2 Φ |R − 1), which is negligible since |m 2 Φ | 1/R 2 . Since the field value at the true vacuum is typically much larger than the EW scale, 2 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 that are given by 1 If the potential of the PQ field itself is unstable, we need to calculate the vacuum decay rate in the PQ sector and add it to that in the THDM sector. In this paper, we assume the stable potential of the PQ field given in eq. (2.1). 2 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.

JHEP03(2020)011
with boundary conditions, (3.4) 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, 3 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 Ω , (3.7) In appendix B, we show that Ω and Θ are constant. 4 Then, the equations of motion reduce to 5 with boundary conditions dφ dr (0) = 0, φ(∞) = 0. (3.12) From eq. (3.7), we can see that a minimum of λ φ satisfies cos 2 Θ = 0 for λ 4 ≥ 0, and cos 2 Θ = 1 for λ 4 < 0. Then, from eq. (3.10), we get the following solutions; JHEP03(2020)011 whereλ = min(λ 3 , λ 3 + λ 4 ). (3.16) Notice that (c) exists only when (λ 1 −λ)/(λ 2 −λ) > 0. If λ φ < 0, the solution to eqs. (3.11) and (3.12) 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 one-loop results of [39][40][41]. The details are in appendix C. 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 at the one-loop level. Here, µ is the renormalization scale and β (1) λ φ is the one-loop beta function for λ φ . Thus, if we integrate it over R ∈ (0, ∞), the integration does not converge. However, as discussed in [39][40][41], 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 up the logarithmic corrections by taking µ ∼ 1/R for each bounce with radius R (for detailed discussion see [41]). If there exists a minimum of the effective action, it dominates the R integral and the result is convergent.
Independently of the convergence of the R integral, we use cut-offs for the R integral for the following reasons. First, we need an IR cut-off because we have ignored the dimensionful couplings. 6 Second, we need a UV cut-off because we do not consider gravitational corrections. Thus, we set the integration region as 7 The effect of the mass term of the bounce field at the false vacuum, m 2 , is discussed in [39] and is shown to be suppressed by R 2 m 2 . We will discuss the cut-off dependence later. 7 We will discuss the cut-off dependence later.

JHEP03(2020)011
where M Pl is the reduced Planck scale. We also impose the same limits on the field value of the bounce as We also exclude the region where the quantum corrections to the action 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 and it always gives a conservative constraint.
The condition for the stability of the EW vacuum is then given by where H 0 67.66 (km/s)/Mpc [44] is the current Hubble constant.

Example
Let us show an example of the calculation. We take We first calculate the MS dimensionless couplings at renormalization scale µ t = m t , where we include the one-loop corrections and the four-loop QCD corrections. The details are in appendix D. Then, we evolve them with the two-loop RG equations. The result is shown in the top left panel of figure 1. In this example, only λ 2 becomes negative and contributes to the vacuum decay rate. 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 figure 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. We use the SM values and uncertainties given in table 3 and α s = 0.1181 (11). We estimate the renormalization scale uncertainty by taking µ = 2/R and µ = 1/(2R). With this parameter set, the vacuum decay rate is close to the upper bound, log 10 [H 4 0 × Gyr Gpc 3 ] −3. Let us see the difference between the "tree level" vacuum decay rate and our result. For the tree level vacuum decay rate, we adopt In figure 2, we show the binned plot of the vacuum decay rates at the tree level and at the one-loop level by using the data accumulated for figure 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 , which is comparable with the uncertainties from those of the top mass and the strong coupling constant. For γ H 4 0 , the vacuum decay rate can be either suppressed or enhanced.

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 JHEP03(2020)011 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.

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 [45][46][47]. We obtain the constraints following the analysis of [48] with the current experimental values. The details are given in appendix A.
In figure 3, we plot the 95% exclusion limits on the (m H + , tan β)-plane with assuming m H = m A = m H + and cos(β −α) = 0. 8 The white region is allowed and the shaded regions are excluded by the observables shown on the regions. As we can see, BR(b → sγ) gives the lower bound of m H + 580 GeV almost independently of tan β. The upper bound and the lower bound on tan β are set by BR(B s → µµ) and ∆M Bs , respectively. Notice that these constraints are stronger than the perturbativity limits of y t , y b √ 4π. The results are consistent with the recent works 9 [45][46][47].

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 [49,50] |λ 1 | < 8π, (4.5) at the tree level. 10 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". JHEP03(2020)011

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 [51,52] (for the THDM, see [53] 11 ) to calculate these parameters. The current constraints are given by [54] S = 0.02 ± 0.07, (4.12) 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.

Collider searches
New scalar particles have been searched extensively at Tevatron, LEP and LHC. We utilize HiggsBounds [55][56][57][58][59] to check the constraints from the collider searches. For simplicity, we consider only on-shell decays for non-SM channels. The couplings and the partial decay widths used in the analysis are summarized in appendix F. For the model to be valid up to the Planck scale, Landau poles should not appear during the RG evolution. We adopt the condition of the tree level perturbative unitarity given in eqs. (4.5)-(4.11) to detect Landau poles and require that they should be satisfied until the Planck scale. We refer to this condition as "the high energy (HE) perturbative unitarity". If it is satisfied, we then check the vacuum stability, where we take µ = 1/R.

High scale validity
We reduce the number of free parameters by choosing three slices of parameter space; (i) m H + = 600 GeV, 1.8 < tan β < 25,

JHEP03(2020)011
which satisfy the flavor constraints of figure 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 two million 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 E.
In figure 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 3 and is equivalent to that in [67]. We expect that the combination of the HE perturbative JHEP03(2020)011  unitarity and the bounded-from-below condition should give a similar result, 12 which we will see below.
In figure 5, we pick up the parameter space defined by the region surrounded by the white dashed lines in figure 4 and prepare additional five million 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 β, cos(β − α)/| cos(β − α)| max , m H − m A and (m H + m A )/2 − m H + . Here, | cos(β − α)| max is the maximum value of cos(β − α) depending on tan β, which is shown in figure 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 figure 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 always becomes unstable at low energy, but the instability is cured at high energy. Notice that the vacuum decay rates can 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 g hU U = sin(β − α) + cot β cos(β − α), (5.7) 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.06 for all the slices, we have 0.9982 g hV V ≤ 1, which is not possible to be distinguished from unity even with HL-LHC plus 1 TeV ILC [68]. It also means that the model cannot 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 figure 6, we cast the data points in figure 5 into the g hU U vs g hDD = g hLL plane. The colors are the same as in figure 5. As we can see, g hU U can be reduced by about 2% − 5% for each slice, but cannot be enhanced so much. On the other hand, g hDD and g hLL can deviate by about 5% − 12% for each slice, and tend to be enhanced.
The current constraints on these couplings are given by [69] g hZZ = 1.10 ± 0.08, (5.11) g hW W = 1.05 ± 0.08, (5.12) 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 [68] 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. (3.21) and (3.22). 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, which we have checked numerically as well. As for the IR cut-off, we have checked that figures 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 figure 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 give an important implication on the scenario.

JHEP03(2020)011 6 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 [39][40][41] 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 12% enhancement of the hDD and hLL couplings, and at most 5% 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.2% and the scenario may be excluded if we observe large deviations of these couplings.

A Flavor
In this appendix, we follow [48] and obtain the flavor constraints with the current experimental values.

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

A.2 Flavor constraints
For the theoretical evaluation of the flavor observables, we use the formulas given in [48].
In the calculation, we utilize RunDec [65,66] 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 JHEP03(2020)011 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 figure 3 and the SM limit. The known parameters {x i ± δx i } are summarized in tables 2 and 3, where we ignore uncertainties of O(0.1%). 13 The result is shown in figure 3.

B Proof of straight bounce
In this appendix, we show that Ω and Θ of the relevant solution do not depend on r. In terms of φ, Ω and Θ, the Euclidean action is expressed as where Ω , Θ and φ are derivatives with respect to r. As introduced in [ Its Euclidean action is given by Let us assume that there exists a minimum, K[φ A ,Ω A ,Θ A ], whereΩ A orΘ A is not constant. Since λ φ (Ω, Θ) is a continuous function, there exist constant Ω B and Θ B satisfying Then, we have The equality holds only whenΩ (r) =Θ (r) = 0 for any r. Notice that when sinΩ A = 0, the field space is not parameterized byΘ A . Then, from eq. (B.7), there exists a bounce with smaller action ifΩ orΘ is not constant. Thus, the bounce with minimum action can only be realized with constant Ω and Θ. 14 13 We have also ignored the uncertainty of the strong coupling constant since its effect is suppressed compared with other uncertainties. 14 Only the bounce with minimum action is relevant for the vacuum decay since the contributions from the others are exponentially suppressed.

JHEP03(2020)011
C One-loop corrections to a vacuum decay rate From [41], the differential vacuum decay rate is expressed as Here, [ln A (X) ] MS 's are defined in [41]. 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.
The scalar contributions: The fermion contributions:

JHEP03(2020)011
The gauge boson contributions: . (C.31) 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 [41] and we replace

D 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 [77,78]. 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 [60,61], FeynArts [62], FeynCalc [63,64].

D.1 Gauge couplings
We first evaluate the electric charge at µ t as It is then matched to the THDM electric charge as Next, we calculate the MS masses of the gauge bosons as 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 , (D.9) . (D.10) Finally, the strong coupling constant is evaluated with RunDec [65,66]. Notice that there are no one-loop threshold corrections from the additional Higgs bosons.
For the later convenience, let us define

D.2 Yukawa couplings
The MS tau mass is obtained from
As for the MS masses of the top quark and the bottom quark, we include the four-loop QCD corrections by using RunDec [65,66]. 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.

D.3 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.

JHEP03(2020)011
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.

E 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. (5.1)-(5.3). 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 figure 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).

E.1 Necessary conditions for perturbative unitarity
Let us first analyze the perturbative unitarity conditions. For arbitrary real numbers A, B and C, the inequality, |A ± B 2 + C 2 | < 1, (E.1)

JHEP03(2020)011
can be reduced to Applying it to the perturbative unitarity constraints, we get constraints on λ 1 and λ 2 as which can be reduced to Here, T = 8π. As for λ 3 and λ 4 , the constraints have the form of Here, a i 's and b i 's are constants and c i 's are functions of λ 1 and λ 2 . The range of λ 3 and λ 4 satisfying eq. (E.10) can be determined by the simplex method of linear programming. We consider simultaneous equations given by where z i 's are the slack variables. Then, we solve them under the constraint of z k = z l = 0 for each pair of (k, l). The solutions satisfying z i ≥ 0 correspond to the corners of the allowed region. We search for such solutions and get T , (E.14)