One-loop considerations for coexisting vacua in the CP conserving 2HDM

The Two-Higgs-Doublet model (2HDM) is a simple and viable extension of the Standard Model with a scalar potential complex enough that two minima may coexist. In this work we investigate if the procedure to identify our vacuum as the global minimum by tree-level formulas carries over to the one-loop corrected potential. In the CP conserving case, we identify two distinct types of coexisting minima — the regular ones (moderate tan β) and the non-regular ones (small or large tan β) — and conclude that the tree level expectation fails only for the non-regular type of coexisting minima. For the regular type, the sign of m122 already precisely indicates which minima is the global one, even at one-loop.


Introduction
After the discovery of the Higgs boson of 125 GeV mass in 2012 [1,2], all the pieces of the SM were firmly established. The Standard Model (SM) however is far from being a complete theory at the smallest physical scales as, to name a few shortcomings, it does not contain a Dark Matter candidate and cannot explain the matter and antimatter asymmetry of the Universe [3].
Among the various proposals, considering an extended Higgs sector is one of the simplest modifications that we can implement in the SM without extending the hidden fundamental forces of nature and also passing the many stringent collider tests at the LHC. In this work we consider one of such extensions namely the Two-Higgs-doublet model JHEP11(2017)106 (2HDM) in which just another scalar doublet is added. This model features five physical Higgs bosons, three neutral and two charged, instead of just one in the SM. It has been extensively studied in the literature (see e.g. ref. [4] for a review) partly because more fundamental theories, for instance the MSSM [5][6][7], require a similar extended scalar sector. The model also allows the possibility of other sources of CP violation [8,9], a feature that gets even richer when more doublet copies are added [10]. Finally, a more complex scalar sector can generate a strong enough EW first-order phase transition [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27], a property that is lacking in the SM [28][29][30][31][32][33] but is necessary to explain the matter-antimatter asymmetry of our universe.
Another feature that a more involved scalar potential encompasses is the possibility of different symmetry breaking patterns and the 2HDM is no exception [34][35][36][37][38][39][40]. With additional Higgs doublets, this complexity increases substantially [41][42][43][44]. In general, for a sufficiently complex potential, there may even exist sets of parameters for which many local minima coexist. In this case, identifying which one is the global minimum might be a nontrivial task that involves solving a system of polynomial equations. Usually, when one is sure that only one minimum is present for a given choice of parameters, such a task is bypassed by trading some of the quadratic parameters in favor of the vevs and ensuring the extremum is a local minimum. In the 2HDM, that assumption does not hold in general and for some parameter ranges it is possible that up to two minima coexist for the same potential [35,36]. So a worrisome possibility arises: we may be living in a metastable vacuum with the possibility to tunnel to the global minimum. That situation was described in ref. [45,46] as our vacuum being a panic vacuum. One way of testing that situation is explicitly calculating the depth of the second minimum and comparing it to the depth of the first. However, finding the location of the minima explicitly may be a difficult or, at least, computationally intensive task. Fortunately, in the same work, the authors developed a method capable of distinguishing if our vacuum is a panic vacuum by calculating a discriminant that depends only on the position of our vacuum (see ref. [47] for a more general test). Although many of the possible scenarios with coexisting minima are not favored by current LHC data [45,46], we are interested here in studying if the simple use of this discriminant can be carried over to the one-loop corrected effective potential. Already for the inert doublet model [34] it was found that the potential difference of the coexisting minima can change sign when one-loop corrections are taken into account [48]. Therefore, the present work aims to verify the validity of such conclusions for a general CP conserving 2HDM with softly broken Z 2 . We focus on the case of two coexisting normal vacua and study the predictive power of tree-level formulas for the depth of the potential when one-loop corrections are considered.
The outline of the paper is as follows: in section 2 we review the properties of the 2HDM with softly broken Z 2 at tree-level focusing on the possibility of two coexisting minima. Some results can not be found in previous literature. In section 3 we review the form of the one-loop effective potential for our case while section 4 explains our procedure for ensuring that our vacuum has the correct vacuum expectation value. The procedure to compute and fix the pole mass of the SM Higgs boson to its experimental value is explained in section 5. The steps we performed to generate the numerical samples are JHEP11(2017)106 listed in section 6 and the resulting analysis is shown in section 7. Finally, the conclusions can be found in section 8.

Coexisting normal vacua at tree-level
The general 2HDM potential at tree-level is (2.1) We will be considering the real softly broken Z 2 symmetric case where m 2 12 and λ 5 are real while λ 6 = λ 7 = 0.
We will also focus on CP conserving vacua where the vacuum expectation values are real and we employ the parametrization These vevs can be further parametrized by modulus and angle as where v = 246 GeV for our vacuum and we use the shorthands c β ≡ cos β, s β ≡ sin β. We call this type of vacuum a normal vacuum and we denote our vacuum by NV [45,46] with v 1 > 0 and v 2 > 0. By ensuring the existence of one normal vacuum (our vacuum), a scalar potential with fixed parameters cannot simultaneously have another minimum of a different type, namely a charge breaking vacuum or a spontaneously CP breaking vacuum [37,38]. Just another coexisting normal vacuum NV with vevs (v 1 , v 2 ) may exist and this is the only case where two minima can coexist in the 2HDM potential at tree level [35,36]: only two minima with the same residual symmetry may coexist. When the coexisting minima exist, we define the potential difference as so that ∆V > 0 indicates that our vacuum is the global minimum. We use this convention for the one-loop potential as well.
To describe the situation of two coexisting normal vacua in more detail, we can write the extremum equations for nonzero v 1 and v 2 : We employ the usual shorthand λ 345 ≡ λ 3 + λ 4 + λ 5 .

JHEP11(2017)106
We will see that there are two types of coexisting normal vacua depending on the Z 2 symmetric limit. The complete solutions of eq. (2.5) for m 2 12 = 0 can be easily found: there are two degenerate extrema that spontaneously break Z 2 -ZB + and ZB − -and two extrema that preserve Z 2 -ZP 1 and ZP 2 ; the latter are often denoted as inert or inert-like vacuum (see ref. [48] and references therein). Only one of the pairs ZP 1,2 or ZB ± may coexist as minima. 1 They are characterized by (v 1 , v 2 ) of the form where we adopt the convention that allv 1 ,v 2 ,ṽ 1 ,ṽ 2 are positive. The specific values of the vevs are given by v 2 for the Z 2 breaking extrema 2 and for the Z 2 preserving minima. The two extrema ZB ± are indeed connected by the spontaneously broken Z 2 symmetry: φ 2 → −φ 2 . We note that simultaneous sign flips of both v 1 , v 2 is a gauge symmetry and do not count as a degeneracy. Hence we adopt the convention that v 1 > 0 while v 2 can attain both signs so that we only analyze the first and fourth quadrant in the (v 1 , v 2 ) plane. As the −m 2 12 v 1 v 2 term is continuously turned on, the Z 2 symmetry is soft but explicitly broken with a negative (positive) contribution in the first (fourth) quadrant when m 2 12 > 0. The opposite is true for negative m 2 12 . The effect of adding the m 2 12 term is different for the two types of coexisting minima which we denote by ZB ± and ZP 1,2 from their m 2 12 → 0 limit. We also denote the ZB ± minima as regular and ZP 1,2 as non-regular simply because it is much more probable to generate models with the former pair than the latter for generic values of tan β and other parameters.
The two degenerate spontaneously breaking minima ZB ± : (v 1 , ±v 2 ) deviate to ZB + : (v 1 , v 2 ) and ZB − : (v 1 , v 2 ), respectively, and the degenerate potential depth, also deviates differently lifting the degeneracy. In first approximation in small m 2 12 and in the deviation of the vevs, the potential depths change respectively by the amount δV ± ≈ ∓m 2 12v 1v2 , so that the depth difference of the two minima is

JHEP11(2017)106
See appendix A for the general formula. As we defined our vacuum to be ZB + , we see that indeed it gets deeper as m 2 12 increases from zero while the non-standard vacuum ZB − is pushed up.
The deviation for the Z 2 preserving minima are different: the first order perturbation to the potential value vanishes. We need the deviations in the locations of the minima ZP 1,2 which, in first order, read (1,2), is the second derivative in the v i direction around ZP j . When m 2 12 > 0, the deviations δv i are positive and the two minima enter the first quadrant. Otherwise they move to the fourth quadrant. The potential depth then changes from the Z 2 limit by respectively. The depth difference of the coexisting ZP 1,2 is (2.16) The first term corresponds to the usual difference between the two inert-like vacua. We conventionally consider ZP 2 to be our vacuum NV corresponding to large tan β. The behaviors described above can be clearly seen in the left panel of figure 1 (blue points) for ZB ± where we show the depth difference calculated exactly against m 2 12 , both normalized by the appropriate power of v = 246 GeV (NV). Only potentials with two minima are selected and the free parameters are taken as with m h = 125 GeV, 1 ≤ tan β ≤ 50, {m H + , m A , m H } ranging from 90 GeV to 1 TeV (m H > m h ), −20,000 GeV 2 ≤ m 2 12 ≤ 6000 GeV 2 and α is constrained near alignment, −0.1 ≤ cos(β − α) ≤ 0.1. Simple bounded from below and perturbativity constraints are also imposed [4]. The blue points end around m 2 12 /v 2 ∼ 0.07 because the nonstandard minimum gets pushed up until the point where it disappears. In contrast, the right panel shows the normalized depth difference with respect to the ratio v /v of the values for NV and NV. We can see for the blue points that the vacuum that lies deeper has a larger vacuum expectation value. The method we employed to calculate the location of NV is described in appendix B.  We confirm the approximation (2.11): for ZB ± the sign of m 2 12 discriminates between our vacuum being the global minimum (m 2 12 > 0) or just a local metastable vacuum (m 2 12 < 0). This behavior, however, does not apply for the points (red) that deviate from the inert-like vacua, ZP 1,2 , where the nonstandard vacuum may lie deeper despite m 2 12 being positive. For the generic values of t β as used above the density of non-regular coexisting minima is very low so that only a handful of coexisting non-regular minima is obtained jointly with the regular points. To generate a sufficient number of non-regular points we further produced another sample (most of the red points) by restricting 20 ≤ t β ≤ 50 and positive m 2 12 . To accurately distinguish among the different cases, ref. [45,46] constructed a very useful discriminant D that ensures that our vacuum is always the global minimum if D is positive. 3 Since that discriminant was derived assuming that v 1 , v 2 are both positive, we cannot apply it to NV when v 2 < 0. So we rederived the discriminant allowing the vevs to be negative with the result where k ≡ (λ 1 /λ 2 ) 1/4 and we have normalized to obtain a dimensionless quantity. This discriminant is useful because it can be obtained by using only the angle β calculated in one vacuum and cases with only one minimum are automatically taken into account. The discriminating power of D is shown in figure 2 where the depth difference is plotted against D calculated using NV. Obviously we could have calculated the discriminant for NV , obtaining a D with sign opposite to D. That implies that the quantity that depends on the vevs, s 2β (t 2 β − k 2 ), must have opposite signs when calculated for NV and NV . Our main goal here is to analyze if the discriminant power of m 2 12 and D carries over to the one-loop effective potential.

Effective potential at one-loop
We can now consider the effective potential with the one-loop contribution The masses M 2 k (ϕ i ) correspond to the scalar-field-dependent eigenvalues of the tree-level mass matrices of all particles of the theory while µ is the renormalization scale. We are already assuming a renormalization scheme with minimal subtraction (MS) and, for the gauge sector, the Landau gauge and dimensional reduction (DRED), following the scheme of ref. [49,50]. The parameters contained in V 0 are thus the renormalized parameters. The integer coefficients |c k | count all the degrees of freedom for each particle k including color, charge and spin, while the sign of c k is determined by its boson/fermion character: positive for bosons and negative for fermions. For example, for the top quark we have c t = −3 × 2 × 2 corresponding to its 3 colors, 2 particle/antiparticle and 2 spin degrees of freedom. We should note that the effective potential is generically a gauge dependent quantity but its value at an extremum is not [51].
As we will focus on normal vacua, we can consider that the effective potential depends only on the two real values ϕ 1 , ϕ 2 in the real neutral directions: 4 4 For a generic field dependence modulo gauge freedom, we would need two more real directions.

JHEP11(2017)106
We reserve the symbols v 1 , v 2 in eq. (2.2) to values at a minimum. So the field-dependent gauge boson masses retain the same functional form as in the SM with v 2 = v 2 1 + v 2 2 : The fermion masses depend on the type of 2HDM we are considering. We only consider models where FCNC are suppressed due to a Z 2 symmetry. We focus more on the type I model but our results are equally valid for all types II, X and Y because of the dominance of the top Yukawa. For the type I we have where y t,b are the Yukawa couplings of the third family quarks normalized to the SM values and the enhancement factor 1/ s β should be considered as the fixed value at the NV minimum at one-loop. We emphasize that information with the bracket . For the type II, we should replace the M b dependence on ϕ 2 and s β by ϕ 1 and c β respectively. We will see that, as usual, the top correction dominates the fermion loops and the difference between type I or type II is negligible for the one-loop corrections except for excessively large tan β which we do not consider. It is also justified that we only consider the effects of the top and bottom quarks; see figure 3 and comments in the text.
For the scalar contribution we need to calculate the eigenvalues of the matrix of second derivatives of V 0 for generic values of ϕ i . These mass matrices are shown in appendix C and their eigenvalues correspond to M 2 S (ϕ i ) of the 8 scalars S ∈ {G ± , H ± , G 0 , A 0 , H 0 , h 0 }. Due to charge and CP conservation, the mass matrices are still separated into three sectors: two charged scalars and its antiparticles, two CP odd scalars and two CP even scalars. We emphasize that e.g. M 2 G 0 (ϕ i ) is nonvanishing at ϕ i away from any tree-level minimum. It is the second derivative of the whole effective potential at one-loop that will vanish in the directions of the Goldstone modes.

Parametrization and minimization at one-loop
We are interested in surveying the cases where the effective potential at one-loop (3.1) continues to have two local minima, one of which should be our vacuum with v = v 2 1 + v 2 2 = 246 GeV. The vevs v 1 , v 2 no longer satisfy the tree-level minimization relations in (2.5) but should now minimize the whole effective potential V 0 + V 1l . We need a convenient parametrization to ensure that one minimum has the appropriate value of v.
To parametrize V 0 , we will use as input the usual 8 quantities where v i satisfy the minimization equations (2.5). It is clear that these quantities define V 0 unambiguously by fixing the 8 parameters {m 2 11 , m 2 22 , m 2 12 , λ 1 , λ 2 , λ 3 , λ 4 , λ 5 }; see e.g. ref. [45]. When we add the one-loop contribution, it is clear that the true minimum will JHEP11(2017)106 be shifted by a small amount from the position (v 1 , v 2 ) at tree-level. Instead of correcting for that shift, we add the finite counterterms to the potential and adjust the values of δm 2 ii so that v 1 , v 2 continue to be a minimum at one-loop. 5 This means that the one-loop effective potential (3.1) is now rewritten as We can see that δm 2 ii → 0 and V 1l → 0 in the limit where we turn off all couplings of the scalars to other particles including self-couplings. For small couplings, it is also expected that the physical masses are close to the masses It is possible to use a different renormalization scheme where all the masses and mixing angles at tree level are maintained at one-loop [52]. Our scheme, however, avoids the need to deal with infrared divergences coming from the vanishing Goldstone masses [22-27, 53, 54]. This problem is more severe at higher loop orders [55]. Now the minimization equations at one-loop can be separated into a tree-level part which leads to the tree-level equations written in (2.5), and a one-loop part that defines We can separate the derivative of V 1l into its contribution from scalars (S), vector bosons (V ) and fermions (F ): for each i = 1, 2. The dimensionless coefficients λ iS are given in appendix D and we use lowercase letters for m W,Z,t,b because they correspond to the actual values in the SM when we use our vacuum NV. For charged particles the contribution of the antiparticle can be taken into account by doubling the contribution of the particle. The fermion part corresponds to the type I model. For the type II model, we must replace s β → c β in the couplings to the b quark. We can see in figure 3 that the contributions from scalars are JHEP11(2017)106 large for δm 2 11 and δm 2 22 while the top contribution is also large and negative for δm 2 22 . The contribution from the bottom quark is negligible for tan β ≤ 50 and there is no appreciable difference between the type I or type II model. Thus for definiteness we consider the type I model. We also note that the scalar masses and coefficients depend on δm 2 11 , δm 2 22 and (4.5) must be solved self-consistently.
The only remaining task is to write M 2 S (v i ) in terms of the input parameters (4.1). We note that M 2 S (v i ) should be computed from the second derivatives of V 0 + δV . But the part coming from V 0 at ϕ i = v i corresponds to the usual masses at tree-level because (v 1 , v 2 ) still corresponds to a minimum of V 0 . Therefore, these matrices will have the generic form Specifically, the mass matrices for the different sectors read c β s β − v 2 λ 5 are the masses squared for the charged Higgs and the pseudoscalar at tree-level; see e.g. [4]. Clearly the first and second matrices contain each a vanishing eigenvalue corresponding to the charged (G ± ) and neutral (G 0 ) Goldstone bosons in the limit where δm 2 ii → 0. We use the basis (φ , (4.9) and (4.10), respectively, from the parametrization (4.12) Using the same notation we can find the eigenvalues shifted by δm 2 ii as where the angles β + , β 0 , α are shifted from β, β, α by a small amount due to δm 2 ii : (4.14) The explicit forms for δβ +,0 and δα can be seen in appendix D.

JHEP11(2017)106 5 Pole masses at one-loop
The previous section showed a way of ensuring that one of the minima of the effective potential at one-loop corresponded to our vacuum with v = 246 GeV and that tan β = v 2 /v 1 could be used as input at one-loop. The following task to describe a realistic 2HDM at one-loop is to ensure that the SM higgs boson mass corresponds to the experimentally measured value [56]: It is clear that at one-loop we cannot use this value for the tree-level parameter m h in (4.1). Instead, we must check that the pole mass corresponds to (5.1). 6 We follow ref. [50] to calculate the pole massesm S of all the scalars S including the SM higgs boson by computing the self-energies of the theory at one-loop. We focus in this section on the CP even sector which will give rise to the pole masses of h and H restricted to the casem H >m h . The self-energy for the other sectors are given in appendix F.
The scalar self-couplings can be extracted from V 0 . Given a set of real scalars S i that interact through the quartic vertex −ig ijkl and cubic vertex −ig ijk , the self-energy Π ij for S i -S j coming from scalars in the loop with incoming momentum p 2 = s is given by [49,50] assuming we are in the basis with diagonalized masses (quadratic part of V 0 + δV ). The A and B functions are the Passarino-Veltman functions [59] which, in the notation of ref. [48], read The A function represents the one-loop graph with one quartic vertex (tadpole) and the B function represents the one-loop graph with two vertices and two internal lines. Hence the factors 1/2 are the symmetry factors that appear in front of the Feynman diagrams for identical fields. These functions appear after renormalizing these diagrams using the MS prescription. Simplifications for vanishing s can be found in ref. [48]. The contributions coming from gauge bosons and fermions in the loop are also shown in appendix F and all the necessary cubic and quartic couplings of the theory are explicitly shown in appendix E. For example, the SM higgs self-energy is given by Another possibility is to use a more physical renormalization condition to ensure this [52]. The renormalization of the entire theory is discussed in ref. [57,58].

JHEP11(2017)106
where g h 2 h 2 = g h 4 and g hh 2 = g h 3 are the quartic and cubic self-couplings for h and all M S refer to M S (v i ). However, the mixing to the other CP even scalar H in Π Hh cannot be neglected. Due to charge and CP conservation, the self-energy for the different scalars decouple into separate pieces for the CP even, CP odd and charged sectors. The self-energy for the CP even sector is given by the matrix The one-loop pole squared massesm 2 H ,m 2 h for the CP even scalars H, h are the solutions s k for [49,50] det where we take only the real part of the self-energy because we are not interested in the decay widths. Them 2 h corresponds to the solution that continually approaches M 2 h in the limit where we turn off the interactions. A similar consideration applies to all pole squared masses.

Numerical survey
We describe here the procedure we used to survey the models. For definiteness, we adopted the relevant fixed parameters of the SM to be g = 0.6483 ; g = 0.3587 ; v = 246.954 GeV ; y t = 0.93697 ; y b = 0.023937 . (6.1) The first four values were taken from refs. [55,60] as the running parameters of the SM at the top pole mass µ = 173.34 GeV and the bottom Yukawa was adapted from its mass value m b = 4.18 GeV [56]. Among these parameters, only the top Yukawa appreciably affects the one-loop effective potential together with the quartic scalar self-couplings. We use a fixed renormalization scale µ = 300 GeV for all calculations and note that the running of y t from the top mass scale only amounts to a small difference. The running for the rest of the parameters are even less relevant. We remark that any choice of the renormalization scale is allowed, since the difference in depth of the potential at two extrema is a renormalization scale-independent quantity [61]. However, from a practical point of view, it is desirable that the logarithms in the effective potential do not become too "large" so as to lead to numerical instabilities [62]. We have checked our calculation for some different values of the renormalization scale and the value of µ = 300 GeV proved to be a stable choice.
Among the input parameters in (4.1), we fixed the standard vev v as above and took the rest of the parameters randomly in the range shown in the first row of table 1, restricted to m H > m h . After checking for simple perturbativity and bounded from below conditions at tree-level, 7 we picked only the points where the shifted masses squared M 2 S (v i ) were 7 Since we work with a fixed renormalization scale where the one-loop corrections are not large, we expect the tree level relations to be valid to a good approximation [64][65][66]. For bounded from below conditions, this is a conservative choice as one-loop corrections may enlarge the possible parameter space [67].

JHEP11(2017)106
Sample Table 1. Input parameter ranges for the general sample (G) and non-regular points (NR). The dash denotes the same range as in the previous row.
positive and the solutions for the shift δm 2 ii in (4.5) were real. Then we further selected only the points where the pole massm h calculated as in section 5 fell in the experimental range of (5.1). 8 In this way, we generated the sample G with 294437 points among which 4525 had two minima at one-loop, 17 of which were non-regular. To find the second minimum, we explicitly minimized the real part of the effective potential (3.1) starting from the non-standard minimum at tree-level and then retained only the points where the value of the potential at that minimum was real [63]. Given the small number of nonregular points, we generated another sample denoted as NR focusing only on non-regular points by imposing large t β as in the second row of table 1; other ranges were kept the same, except for m 2 12 which was chosen positive. Sample NR thus contained 185905 points among which 1563 had two minima at one-loop. Hereafter, unless explicitly specified, we will consider only the joint sample of G and NR. However, we would like to emphasize that, even though samples G and NR have a similar number of points, the parameter space scanned by sample NR represents only a small portion of the parameter space probed by sample G. This justifies our definition of non-regular points.
After the selections described above, the input parameters m H + , m A , m H get roughly confined to the range [90, 500] and the non-standard higgsses acquire pole masses in a similar range with slightly smaller maximal value form H + . In contrast, the distribution for t β is homogeneous in the range of table 1 for the whole sample but, as we select only the points with two minima at one-loop, it gets separated into two ranges, 1 ≤ t β 3.8 for the regular minima, and 24 t β ≤ 50 for the non-regular ones.
To check our numerically implemented formulas, we performed the following consistency checks: 1. Vanishing of the pole masses for the Goldstone bosons G ± , G 0 for zero external momentum for the three cases where we successively add the one-loop scalar, vector boson and fermion contributions.

Results
Let us first quantify the shifts δm 2 ii in (4.5) for the different contributions. The contribution coming from the gauge bosons only depends on the value of v and, for fermions, it depends on v and β. The scalar contribution depends on many parameters coming from the scalar potential. The dependence of the different contributions on tan β are shown in figure 3. We can see that the dominant contribution for δm 2 11 comes from the scalars whereas δm 2 22 also has large positive contributions from scalars but they are partly canceled by the negative contribution from fermions (top). Such a partial cancellation in δm 2 22 can be clearly seen in figure 4 where we show the contribution only from scalars (green points) and all the JHEP11(2017)106 contributions (red points). The orange curve in figure 3, which quantifies the bottom contribution to δm 2 11 in the type II model, shows that it is negligible in the range of tan β we are interested in and our calculation that uses the type I model applies equally well to the type II case. All the different contributions are calculated using eq. (4.6) and the scalar contribution in particular depends on the δm 2 ii themselves and these are taken as the total contributions. The remaining contributions of gauge bosons (purple dashed curve in figure 3) are much smaller. The contribution from scalars are calculated using the whole sample (G + NR) described in section 6 which also includes the points with only one vacuum.
The deviation of the location of our vacuum when we add δV do V 0 is illustrated in figure 5 where we show the ratio of t β for V 0 +δV (t tree+δ β ) to that of V 0 (t β ) against the ratio of the vev for V 0 +δV (v tree+δ ) to that of V 0 (v). We only show the points with two coexisting minima and divide the points between the regular ones (blue) and the non-regular ones (red). We can see that as m 2 ii get shifted by δm 2 ii all points with two minima have their vevs decreased while t β mostly increases for the regular points and mostly decreases for the non-regular points. Note that the non-regular points only consider large t β whereas the regular points only include moderate t β roughly up to 3.8. As the location of (v 1 , v 2 ) for our vacuum is the same for V 0 and the one-loop corrected potential (4.3) we can also interpret this plot as the modification of the vev location of V 0 + δV compared to V 0 + δV + V 1l . If we had considered all the points including the points with only our vacuum, the majority of points would follow the behavior of the regular points in blue.
For the case of coexisting vacua, we can see a clear difference in behavior between the regular points and the non-regular ones in figure 6 where we plot the potential difference with respect to m 2 12 . The left panel shows these quantities using the tree-level potential V 0 while the right panel is the same plot for the one-loop potential (4.3). We can clearly see that the potential difference goes continuously to zero as m 2 12 → 0 for the blue points representing the regular minima. Moreover, for positive m 2 12 our vacuum is guaranteed to be the global minimum in both tree level or one-loop potential. This behavior is opposite JHEP11(2017)106  for negative m 2 12 . The reason is that only m 2 12 breaks explicitly the Z 2 symmetry of the theory and it controls the degeneracy breaking of the spontaneously breaking minima. This behavior is not followed by the non-regular minima that do not have degenerate minima in the m 2 12 → 0 limit. Some points (green) in which our vacuum is not the global one at tree-level even get inverted and become the global minimum as the one-loop corrections are added. The same conclusion is reached if we had compared the one-loop potential to V 0 + δV instead: some cases where our minimum is not the deepest at tree-level becomes the global minimum at one-loop.
We can have an idea of the different one-loop contributions in figure 7 where we separate the potential difference in figure 6 into its different contributions. Regarding the regular points (left plot), it can be seen that the one-loop potential difference is almost entirely due JHEP11(2017)106 to the tree-level contribution (blue) since the contributions from δV (yellow), fermions (red) and scalars (purple) approximately cancel each other while the contribution from gauge bosons (green) is negligible compared to the others. This behavior justifies our choice for the renormalization scale. For the non-regular points (right plot), no clear pattern emerges. We can also see in figure 8 that there are points where the potential difference is raised as well as points where it is lowered by the one-loop corrections for both regular and non-regular points.
Considering that the discriminant (2.18) test is applicable for all cases where at least one normal vacuum is known, we can investigate if it is still a good predictor for the oneloop potential with two coexisting minima. We can adapt the tree-level discriminant to one-loop in the following three ways: The quantity D 1 is the discriminant calculated with V 0 as the whole potential while D 2 is calculated by using V 0 + δV in (4.3). In the latter case, the quadratic parameters m 2 ii get shifted and the location of the minimum, denoted above as (v (0) , β (0) ) or as tree + δ (superscript/subscript) in figure 5, do not coincide with the ones used as input, denoted as (v, β). The last adaptation D 3 considers the shifts in the quadratic parameters but keeps the vevs as (v, β). We will test here if any of these discriminants are capable of distinguishing if our vacuum is the global one at one-loop only using parameters at treelevel. 9 9 Strictly speaking, some calculation at one-loop is required for some of these quantities depending on how the calculation is set up. Using the splitting of the potential in the form (4.3), D1 is the most natural quantity to use and there is no one-loop calculation required. If V0 + δV is considered as the potential at tree-level, D2 or D3 are the natural quantities depending on which minimum is taken.  For the regular coexisting minima, figure 6 shows that m 2 12 is already a good predictor of the global minimum, but we can test the discriminants in (7.1). The result of this test is shown in figure 9 where the potential difference at one-loop is plotted against the three discriminants. We can see in the left and middle plots that D 1 and D 2 correctly predict the global minimum of the one-loop effective potential while the right plot shows that D 3 fails for more than 10% of the points.

JHEP11(2017)106
In constrast, the failure of all the discriminants (7.1) for the non-regular points (of samples G and NR) can be seen in figure 10 which shows the potential depth difference as a function of the discriminants, similarly to the regular points in figure 9; note that the horizontal scale is very different. The green points in the second and fourth quadrants mark the cases where the discriminant D 1 predicts the opposite behavior at one-loop. We can clearly see that all discriminants fail for a significant portion of points, not necessarily for the same ones. From the property of D 1 , we could have seen its wrong prediction in figure 6 as well.
At last, to make sure that the points for which the discriminant test fails include phenomenologically realistic points, we have checked the viability of the red/green points in figure 6 by considering phenomenological constraints implemented in the 2HDMC code [68,69]. We found that in the case of the type I model around half of the green points are allowed by experimental constraints such as the S, T, U precision electroweak parameters JHEP11(2017)106 and data from colliders implemented in the HiggsBounds and HiggsSignals packages [70][71][72][73], while in the type II model they are all excluded. For the type II model we have also included constraints coming from R b measurements [74] as well as B meson decays [75][76][77][78] which prove to be very strong by setting m H + > 480 GeV independently of the value of tan β. For type I, since the red/green points have tan β > 10, these constraints do not impose any further restrictions [79,80]. We also checked that all points still respect simple bounded from below and perturbativity constraints [4].

Conclusions
We have studied the one-loop properties of the real Two-Higgs-doublet model with softly broken Z 2 with respect to the possibility of two coexisting normal vacua. The softly broken nature must remain at one-loop and the case of two coexisting normal minima can be classified into two very distinct types depending on their nature in the vanishing m 2 12 limit: regular minima that spontaneously break the symmetry and the non-regular minima (or minimum) that preserve the symmetry, i.e., they are inert-like in the symmetric limit. Since in the first case the two minima are degenerate in the symmetry limit, even at oneloop, they are connected by the Z 2 symmetry and then they should differ only by the sign of v 2 . After the inclusion of the m 2 12 term, the sign of v 2 continue to be opposite for our vacuum and the non-standard one so that the two regular minima are found in the first and fourth quadrants in the (v 1 , v 2 ) plane. In contrast, the non-regular minima deviate from the inert-like minima and both deviate to the first quadrant when m 2 12 is positive. The vacua that spontaneously break Z 2 in the m 2 12 → 0 limit behave rather regularly and we can distinguish which coexisting minimum is the global one by just examining the sign of m 2 12 : when it is negative our vacuum is only a metastable local one and the opposite is true if it is positive. For this type of coexisting vacua, the discriminant at tree level (D 1,2 in eq. (7.1)) is still a good predictor of the nature of the minimum it is calculated with.
For the non-regular coexisting vacua, m 2 12 is positive in the convention that our vacuum has both vevs positive and it cannot be used as an indicator. At tree level, the discriminant of ref. [45,46] is a very convenient way of testing if our vacuum is the global one because only the location of our minimum is required. However, at one-loop, this discriminant is not a precise indicator of which minima is the global one. We have found realistic cases where our vacuum is not the global minimum at tree-level but it becomes the global one after the addition of the one-loop corrections. Few cases for the opposite behavior were also found. As the discriminant effectively distinguishes the sign of the potential difference between the coexisting minima at tree-level, the latter itself is also a not good indicator for the non-regular minima. This is reminiscent of the exact Z 2 symmetric case (inert model) investigated in ref. [48]. On the other hand, we were unable to find a discriminant that works for both regular and non-regular minima at one-loop. Finding a simple and precise criterion for global minimum does not seem to be an easy task as that was not achieved even in the simpler exact Z 2 limit. We also emphasize that for our parametrization which enforces our vacuum to be a minimum from the start, and for the chosen generic ranges of parameters (sample G), the occurrence of non-regular minima, as suggested by its name,
In summary, the soft-breaking term m 2 12 controls the lifting of the degeneracy of the regular coexisting minima (moderate tan β) even at one-loop and can be used as the sole indicator of which minimum is the global one. That is not true when two coexisting nonregular vacua (large or small tan β) exist: the discriminant that is a precise indicator at tree-level is not reliable at one-loop and explicit calculation of the potential depths must be carried out.

A Deviation of a potential minimum
Take a potential V 0 (ϕ) depending on real scalar fields ϕ i for which we know a minimum (extremum) ϕ i =φ i satisfying We want to quantify how the location of the minimum and the value of the potential deviate when we add a small perturbation U on the potential as We assume V 0 has no flat direction aroundφ. We first quantify the deviation of the location of the minimum as The derivative of the perturbed potential gives Since the first term vanishes due to (A.1), we get the deviation to first order where U i ≡ ∂U (φ)/∂ϕ i and M ij = ∂ 2 V 0 (φ) ∂ϕ i ∂ϕ j is the squared-mass matrix aroundφ. The deviation in the value of the minimum can be equally expanded aroundφ as

JHEP11(2017)106
after using (A.5). Generically the third term contributes to deepen the potential. That is the dominant contribution when the perturbation vanishes on the unperturbed minimum: U (φ) = 0. The latter happens for the m 2 12 term when perturbing the inert-like minima while for the spontaneously broken Z 2 minima the dominant term is the second, linear in the perturbation.
B Finding more than one normal vacua Let us rewrite the minimization equations in eq. (2.5) as and ζ ≡ 2m 2 12 /v 1 v 2 depends on the vevs. To find solutions (v 1 , v 2 ) of (B.1) for m 2 12 = 0, we first need an equation for ζ. For that end, we can formally write (B.1) as X = A −1 ζ M and equate v 2 1 v 2 2 = 4m 4 12 /ζ 2 . We obtain where µ 1 = λ 345 m 2 11 − λ 1 m 2 22 and µ 2 = λ 345 m 2 22 − λ 2 m 2 11 . The m 2 12 → 0 limit is taken as ζ → 0 and We know there are at most four real solutions coinciding with the maximal number of extrema in this case [35,36]. We also can see that the root λ 345 + √ λ 1 λ 2 from the lefthandside is always positive due to bounded below conditions. Possible solutions for (B.1) depends on the sign of ζ and we allow both signs for v 1 v 2 . We can see that for the same ζ, flipping the sign of m 2 12 is equivalent to flipping the sign of v 1 v 2 . Once some solution ζ = ζ 0 is found, we can find the vevs from the relation The mass matrices at tree-level for the charged, CP-odd and CP-even scalar sectors prior to imposing the minimization conditions are given respectively by The term (quadratic) refers to the quadratic contribution given by The derivatives are with respect to the fields around the values in eq. (2.2).

D Cubic derivatives
The coefficients λ iS in (4.6) are defined as where S = {G + , H + , G 0 , A, H, h}, sec = {char, odd, even} refers to the different scalar sectors and U sec diagonalizes M sec . The last subindex SS refers to one of the diagonal entries following the ordering in (4.13). Explicit calculation leads to

JHEP11(2017)106
where the mixing angles β + , β 0 , α were defined in (4.14) with The coefficients λ 2S can be obtained from λ 1S for each scalar S above by using the replacement Q 12 in (G.4).

E Cubic and quartic couplings
Given that the cubic and quartic couplings do not depend on the quadratic parameters, they correspond to the tree-level ones listed, e.g., in ref. [81] if we are in the limit α → α, β 0,+ → β. If we want the couplings with these corrections, we can adapt the rotation angles α → α or β → β + or β → β 0 for the couplings that do not mix the charged sector with the CP-odd sector because in the latter there is an ambiguity in distinguishing β + from β 0 . Another ambiguity arises if the couplings are written in terms of (v, β) instead of (v 1 , v 2 ) because then we need to distinguish β coming from the vevs and the β 0,+ coming from the diagonalization of the shifted mass matrices. We adopt the convention that −ig ABCD is the Feynman rule associated to the vertex ABCD. This is opposite to e.g. ref. [81]. We also abbreviate e.g. g G 0 G 0 G 0 G 0 as g G 0 4 . The essential set of quartic couplings is

JHEP11(2017)106
Coupling Obtainable from Replacement More couplings can be obtained from simple replacements; see table 2.
The essential set of cubic couplings is The remaining couplings can be obtained through reparametrization symmetries as shown in table 3.
Finally we note that, although the reparametrization symmetry already allows a huge simplification in the computation of the cubic and quartic couplings needed for our calculation, their use can be error-prone. In our routines we have adopted a different approach, namely we expanded the tree-level scalar potential in terms of physical fields S = {h, H, A, G 0 , G + , H + , G − , H − } and performed derivatives to obtain the desired couplings. For instance,

F Self-energy for the scalars
We show here the different contributions to the self-energy of scalars due to scalars (S), fermions (F ) and gauge bosons (V ) in the loop. We first list the contributions from scalars in the loop for the different sectors. For the CP odd sector we have Some couplings are absent because CP is conserved and A, G 0 are CP odd.

JHEP11(2017)106
In the CP even sector we have  (F.2f) The self-energy for the charged sector is Note that couplings such as g G + G − A 0 = 0 due to CP conservation. The fermionic corrections to the propagator of a scalar S to a scalar S are given by the general formula [49,50]:  Table 4. Yukawa couplings.
Type I Type II s β+ /s β c β+ /c β We denote the vertex of the scalar S k to the fermionsf and f by −iY k f f and they may contain the γ 5 matrix whileȲ refers to the transformationΓ = γ 0 Γ † γ 0 in spinor space.  G ± s β+−β0 c β+−β0 Table 6. Coefficients for the couplings between two scalars and one gauge boson. Table 7. Coefficients for the couplings between one scalar and two gauge bosons.
where there is an implicit summation on S . The loop functions B SV (m V , m S , s) and B V V (m V , m V , s) are defined as (F.10) The couplings C x SS and C x S can be read from the following tables.

G Reparametrization symmetry
Here we list some useful reparametrization symmetries that allow us to relate different quartic and cubic scalar couplings that are numerous. These relations help us to check different couplings or deduce new ones from a smaller set. Moreover, given simple relations, they minimize errors and speed up numerical implementation. The simplest reparametrization is to exchange fields of the same type, one in φ 1 and the other in φ 2 , together with a 90 • shift in the respective mixing angle: This is what allowed us to relate g HG 0 2 to g hG 0 2 in (E.2a) after α → α − π/2 since g hG 0 2 · hG 0 2 → g hG 0 2 | α →α −π/2 · HG 0 2 by the inverse transformation of Q even . Another example is g HhAG 0 in (E.1k): it is odd by the replacement α → α + π/2 or β 0 → β 0 + π/2. The reparametrization symmetry above arises by noting that the original fields in φ 1,2 are left invariant by the transformations and consequently the potential is also invariant.

JHEP11(2017)106
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.