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 (SM) 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\beta$) and the non-regular ones (small or large $\tan\beta$) --- and conclude that the tree level expectation fails only for the non-regular type of coexisting minima. For the regular type, the sign of $m^2_{12}$ already precisely indicates which minima is the global one, even at one-loop.

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 m 2 12 already precisely indicates which minima is the global one, even at one-loop.

I. INTRODUCTION
After the discovery of the Higgs boson of 125 GeV mass in 2012 [1], 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 [2].
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 (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. [3] for a review) partly because more fundamental theories, for instance the MSSM [4], require a similar extended scalar sector. The model also allows the possibility of other sources of CP violation [5], a feature that gets even richer when more doublet copies are added [6]. Finally, a more complex scalar sector can generate a strong enough EW first-order phase transition [7,8], a property that is lacking in the SM [9] 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 [10][11][12][13][14]. With additional Higgs doublets, this complexity increases substantially [15]. 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 [11]. 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. [16] 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. [17] for a more general test). Although many of the possible scenarios with coexisting minima are not favored by current LHC data [16], 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 [10] it was found that the potential difference of the coexisting minima can change sign when one-loop corrections are taken into account [18]. 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 Sec. II 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 Sec. III we review the form of the one-loop effective potential for our case while Sec. IV 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 Sec. V. The steps we performed to generate the numerical samples are listed in Sec. VI and the resulting analysis is shown in Sec. VII. Finally, the conclusions can be found in Sec. VIII.

II. COEXISTING NORMAL VACUA AT TREE-LEVEL
The general 2HDM potential at tree-level is 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 [16] 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 [12]. 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 [11]: 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 . We will see that there are two types of coexisting normal vacua depending on the Z 2 symmetric limit. The complete solutions of Eq. (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. [18] 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 12v1v2 , so that the depth difference of the two minima is 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 where with m 2 Hi = m 2 ii + 1 2 λ 345ṽ 2 j , (ij) = (2,1) or (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 2 As long as the solutions for v 2 i give positive solutions.
respectively. The depth difference of the coexisting ZP 1,2 is 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 Fig. 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 [3]. 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 (10): 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. [16] 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 3 For D = 0 but m 2 12 = 0 the discriminant is inconclusive.
automatically taken into account. The discriminating power of D is shown in Fig. 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.

III. 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. [19]. 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 [20].
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 : We reserve the symbols v 1 , v 2 in Eq. (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 Fig. 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.

IV. PARAMETRIZATION AND MINIMIZATION AT ONE-LOOP
We are interested in surveying the cases where the effective potential at one-loop (18) 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 (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 (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. the first reference in [16]. When we add the one-loop contribution, it is clear that the true minimum will be shifted by a small amount from the position (v 1 , v 2 ) at treelevel. 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 (18) 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 [21]. Our scheme, however, avoids the need to deal with infrared divergences coming from the vanishing Goldstone masses [8,22]. This problem is more severe at higher loop orders [23]. Now the minimization equations at one-loop can be separated into a tree-level part which leads to the tree-level equations written in (5), and a one-loop part that defines δm 2 ii by 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 Fig. 3 that the contributions from scalars are 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 (27) must be solved self-consistently. The only remaining task is to write M 2 S (v i ) in terms of the input parameters (23). 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 where c β s β −v 2 λ 5 are the masses squared for the charged Higgs and the pseudoscalar at tree-level; see e.g. [3]. 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 (φ (30), (31) and (32), respectively, from the parametrization where 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 : The explicit forms for δβ +,0 and δα can be seen in appendix D.

V. 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 [24]: It is clear that at one-loop we cannot use this value for the tree-level parameter m h in (23). Instead, we must check that the pole mass corresponds to (37). 6 We follow Ref. [19, b] 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 [19] 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 [26] which, in the notation of Ref. [18], 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. [18]. 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 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 [19] 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.

VI. 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 .
The first four values were taken from Refs. [23,27] 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 [24]. 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 [28]. 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 [29].
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 (23), 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 I, 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 positive and the solutions for the shift δm 2 ii in (27) were real. Then we further selected only the points where the pole massm h calculated as in Sec. V fell in the experimental range of (37). 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 (18) 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 [30]. Given the small number of non-regular 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 I; 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 I for the whole sample but, as we select only 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 [31]. For bounded from below conditions, this is a conservative choice as one-loop corrections may enlarge the possible parameter space [32]. 8 The exact adopted procedure differs slightly with respect to the range of m h : instead of post-selecting only the values of m h for which the pole mass coincided with the experimental range, we randomly selected the input m h in the range [0, 200] GeV and then later varied only this value searching for the correct pole massm h . If a solution were found, we kept that point. This procedure resulted in the approximately homogeneous distribution of m h in the range shown in Table I. This modification speeded up the generation of points.
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.
2. Equality of the pole masses for the CP even higgsses H, h for zero external momentum and the eigenvalues of the explicitly calculated second derivative matrix of the effective potential in the real neutral directions (20) in all three cases of successive addition of the one-loop scalar, vector boson and fermion contributions.

VII. RESULTS
Let us first quantify the shifts δm 2 ii in (27) 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 Fig. 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 Fig. 4 where we show the contribution only from scalars (green points) and all the contributions (red points). The orange curve in Fig. 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. (28) 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 Fig. 3) are much smaller. The contribution from scalars are calculated using the whole sample (G + NR) described in Sec. VI which also includes the points with only one vacuum.  (27). The scalar contributions are positive for δm 2 11 (blue points) and mostly positive for δm 2 22 (green points) while the fermions contribute negatively. The fermion contribution to δm 2 22 (red curve) is practically the same for the type I and II models while the contribution to δm 2 11 (orange curve) applies only to the type II case but vanishes for the type I case. The contributions from the gauge bosons are shown in the purple dashed line.
The deviation of the location of our vacuum when we add δV do V 0 is illustrated in Fig. 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 (25) 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. 2 22 in (27): the green points quantify the scalar contribution only whereas the red points consider all the contribution in the type I or II model. For the case of coexisting vacua, we can see a clear difference in behavior between the regular points and the non-regular ones in Fig. 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 (25). 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 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.

FIG. 4: Different contributions to δm
We can have an idea of the different one-loop contributions in Fig. 7 where we separate the potential difference in Fig. 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 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 nonregular points (right plot), no clear pattern emerges. We can also see in Fig. 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 (17) 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 one-loop potential with two coexisting minima. We can adapt  (1) is considered. This plot should be compared with Fig. 1. Right: The full 1-loop effective potential (25) is considered. The green points represent a change in sign of the tree level prediction for the potential difference. 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 (25). 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 Fig. 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 tree-level. 9 For the regular coexisting minima, Fig. 6 shows that m 2 12 is already a good predictor of the global minimum, but we can test the discriminants in (45). The result of this test is shown in Fig. 9 where the potential difference at one-loop 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 (25), D 1 is the most natural quantity to use and there is no one-loop calculation required. If V 0 + δV is considered as the potential at tree-level, D 2 or D 3 are the natural quantities depending on which minimum is taken.
FIG. 8: One-loop correction to the potential difference relative to the tree level potential difference. The color coding follows Fig. 6.
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. In constrast, the failure of all the discriminants (45) for the non-regular points (of samples G and NR) can be seen in Fig. 10 which shows the potential depth difference as a function of the discriminants, similarly to the regular points in Fig. 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 Fig. 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 Fig. 6 by considering phenomenological constraints implemented in the 2HDMC code [33]. 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 and data from colliders implemented in the HiggsBounds and HiggsSignals packages [34,35], 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 [36] as well as B meson decays [37] 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 [38]. We also checked that all points still respect simple bounded from below and perturbativity constraints [3].

VIII. 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 one-loop, 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. (45)] 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. [16] 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. [18]. 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, is much rarer: only 38% correspond to the non-regular cases and, among then, only 0.3% are coexisting minima.
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 non-regular 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.
Appendix 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 (A1), we get the deviation to first order where U i ≡ ∂U (φ)/∂ϕ i and M ij = ∂ 2 V0(φ) ∂ϕi∂ϕj is the squared-mass matrix aroundφ. The deviation in the value of the minimum can be equally expanded aroundφ as after using (A5). 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.
Appendix B: Finding more than one normal vacua Let us rewrite the minimization equations in Eq. (5) as where and ζ ≡ 2m 2 12 /v 1 v 2 depends on the vevs. To find solutions (v 1 , v 2 ) of (B1) for m 2 12 = 0, we first need an equation for ζ. For that end, we can formally write (B1) 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 m 4 12 ζ 2 → v 2 1 v 2 2 /4. When m 2 12 = 0, we can find the possible values for ζ (extrema) from the quartic equation We know there are at most four real solutions coinciding with the maximal number of extrema in this case [11]. We also can see that the root λ 345 + √ λ 1 λ 2 from the lefthandside is always positive due to bounded below conditions. Possible solutions for (B1) 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 X = A −1 ζ M or, explicitly, After ensuring these expression are positive, we can extract v 1 , v 2 and tan β with the sign convention where v 1 > 0 and sign(v 2 ) = sign(ζm 2 12 ) .
Appendix C: Matrix of second derivatives 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).

Appendix D: Cubic derivatives
The coefficients λ iS in (28) 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 (35).
Explicit calculation leads to where the mixing angles β + , β 0 , α were defined in (36) with The coefficients λ 2S can be obtained from λ 1S for each scalar S above by using the replacement Q 12 in (G4).

Appendix 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. [39] 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. [39]. We also abbreviate e.g. g G 0 G 0 G 0 G 0 as g G 04 . The essential set of quartic couplings is More couplings can be obtained from simple replacements; see table II.

Coupling
Obtainable from Replacement The essential set of cubic couplings is The remaining couplings can be obtained through reparametrization symmetries as shown in table III.
Coupling Obtainable from Replacement Coupling Obtainable from Replacement 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, Appendix 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.
In the CP even sector we have 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 [19]: where N c is the number of colors of the fermion in the loop, the B function was given in (40) 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 whilē Y refers to the transformationΓ = γ 0 Γ † γ 0 in spinor space. These Yukawa couplings are listed in table IV where the coefficients C S f depend on the model used as in table V.
The couplings C x SS and C x S can be read from the following tables.