One-Loop Charge-Breaking Minima in the Two-Higgs Doublet Model

We analyze the vacuum structure of the one-loop effective potential in the two Higgs doublet model. We find that electroweak-breaking vacuua can coexist with charge breaking ones, contradicting a theorem valid at tree-level. We perform a numerical analysis of the model and supply explicit parameter values for which charge-breaking vacuua can be the global minimum of the theory, and deeper than charge-preserving ones.


Introduction
The discovery by the LHC collaborations of the Higgs boson [1,2] provided the missing piece of the puzzle for the Standard Model (SM) of particle physics. Since then, measurements of the Higgs' properties [3] have shown that this scalar, with mass around 125 GeV, behaves largely as expected in the minimal SM: thus far, and within the measured precision, no significant deviations from SM-like behavior have been observed. But the SM leaves a great many questions unanswered, such as the origin of the matter-antimatter asymmetry, the nature of dark matter as a particle, the observed fermion mass hierarchy, and the strong CP problem. SM extensions are therefore of interest to attempt to provide answers to these, and other, unsolved problems. Models with extended scalar sectors, in particular, are quite popular and widely studied in the literature. One of the simplest beyond the SM theories is the two Higgs doublet model (2HDM), first proposed by Lee in 1973 [4] to provide an additional source of CP violation stemming from the scalar sector through spontaneous symmetry breaking.
In the 2HDM, the gauge and fermion content are the same as in the SM, but instead of a single SU(2) doublet with hypercharge Y = 1,we now have two, Φ 1 and Φ 2 . This leads to a rich phenomenology (see [5] for a review), boasting a richer scalar spectrum than the SM's, with two CP-even scalars, a pseudo-scalar and a charged scalar. The model can have tree-level flavor-changing neutral currents mediated by scalars, can provide a dark matter candidate and have spontaneous CP violation. The 2HDM easily reproduces all experimental results from the SM, and indeed it has a decoupling limit where the extra scalars are very massive and the model's predictions can be made to be virtually indistinguishable from those of the SM. The 2HDM also has a richer vacuum structure then the SM -whereas in the SM the only possible vacuum is the one which breaks electroweak symmetry, in the 2HDM spontaneous CP breaking is also possible, as well as minima where the electromagnetic symmetry U (1) em is broken. These latter minima are unwanted since they would imply a massive photon. The possibility of charge (and color) breaking already arises in SUSY models, and leads to bounds on some of the parameters of the model [6].
The possibility of reducing the 2HDM parameter space in a similar manner -imposing bounds on the model's parameters to avoid global charge-breaking (CB) minima -is very appealing. In many cases, sufficient conditions to avoid CB minima were considered [4,7,8], but at the time it wasn't known whether such conditions were too restrictive. However, in [9,10] a remarkable result was obtained: the structure of the tree-level 2HDM scalar potential is such that, if an electroweak breaking minimum exists, any CB extremum that might then occur is necessarily a saddle point lying above the minimum. Likewise, it was shown that if a CB minimum exists, any stationary point which would break the normal electroweak symmetries is a saddle point lying above it. Analogous results were also proved for the relationship between electroweak extrema and CP breaking ones. Thus a 2HDM electroweak breaking minimum, if it exists, is guaranteed, at tree level, to be stable against tunneling to deeper CB or CP breaking vacuua, since such deeper minima were shown to not exist. This result was further studied in refs. [11][12][13][14][15]. In particular, using a Minkowski formalism to rewrite the 2HDM scalar potential, Ivanov was able to show [13,14] the stability of the different vacuua through geometric arguments. Other results concerning neutral minima in the 2HDM were also obtained -it was shown [13,14,16] that neutral minima can coexist in the 2HDM scalar potential, provided they break the same symmetries. This had implications for the Inert model [17][18][19], a version of the 2HDM where a discrete Z 2 symmetry is preserved by both the Lagrangian and by the vacuum -the vacuum preserves Z 2 since only one of the doublets acquires a non-zero vacuum expectation value. Two possibilities for minima then arise, depending on which of the doublets has the non-zero VEV, the Inert Minimum (where fermions acquire mass after spontaneous symmetry breaking) and the Inert-Like Minimum (where the fermions remain massless). In [14,16] expressions relating the depth of the potential at each of these minima were found, and it was shown that, for specific regions of parameter space, they could coexist.
However, powerful though the demonstrations of [9,10] and [13,14] were, those works dealt with the tree-level potential. The expressions found there comparing the depth of the potential at different extrema depended heavily on tree-level formulae for the scalar masses; for the minimization conditions determining the vacuum expectation values (VEVs); and for the potential itself. A valid question is therefore whether these results are robust when one considers loop corrections to the potential -will the stability theorems deduced for the tree-level potential still hold at one-loop? The first hint that that may not be the case was obtained in [20], where a one-loop calculation was undertaken to analyze the coexistence of neutral minima in the Inert model. The effective potential formalism was employed and it was shown that, in certain cases, a tree-level local minimum could become a one-loop global one, and vice-versa. This (rare) possibility occurred only for regions of parameter space where the tree-level minima were close to degenerate, hence it did not correspond to a breakdown in perturbation theory, rather it implied loop corrections could change the nature of tree-level vacuua. Further, the one-loop calculation enlarged the region of parameter space for which different neutral minima could coexist.
The purpose of this paper is to investigate, using the one-loop effective potential, whether the conclusions concerning the (non-)coexistence of neutral and charge breaking vacuua in the 2HDM hold when radiative corrections are taken into account. We review the tree-level results for the classical 2HDM potential in section 2 then proceed to review the formalism of the one-loop effective potential in section 3, including a discussion of issues related with gauge fixing. The numerical methods we use to carry the minimization of the one-loop potential are detailed in section 4 where we also present results of numerical scans of the model's parameter space and give a few illuminating examples. We draw our conclusions in section 6.
2 The tree-level vacuum structure of the 2HDM The 2HDM contains two hypercharge 1 SU (2) scalar doublets, and the most general scalar potential one can write has a total of 14 real parameters. Since both doublets are identical, any linear combination of them which preserves the scalar kinetic terms should lead to the same physics. This basis invariance, which corresponds to a redefinition of the fields via a 2 × 2 unitary matrix U 1 , allows one to reduce the number of free parameters to 11 [21]. This most general 2HDM will include flavor changing neutral currents (FCNC), mediated by neutral scalars at tree-level, when one considers the full Lagrangian, including fermions. To prevent this, a discrete Z 2 symmetry is introduced [22,23] which is extended to the Yukawa sector in such a way that each class of same-charge fermions (up and down-type quarks and charged leptons) only couple to one of the doublets. This eliminates tree-level FCNC and, due to the several possibilities of extending Z 2 to the Yukawa sector leads to four types of 2HDMs (type I, type II, lepton specific and flipped [5]). So that the model can possess a decoupling limit [8], a softly Z 2 breaking quadratic term, m 2 12 , is usually introduced, so that the scalar potential is characterized by 8 real independent parameters.

Classical Potential
The 2HDM scalar potential we will be studying possesses a softly broken Z 2 symmetry and is therefore given, at tree-level, by where all the parameters are taken to be real 2 . The scalar doublets, Φ 1 and Φ 2 contain a combined eight real component fields which can be parameterized as follows: It is well know that the 2HDM classical potential exhibits three-different types of extrema (for instance, see section 5.8 of [5] for a demonstration): a U (1) EM and CP-conserving extremum, which we denote by Φ i EW with the vev α = 0 breaking electric charge conservation and consequently giving a mass to the photon and a CP-violating extremum: where δ = 0. In this work, we will be focusing on the EW and CB extrema.

Classical Extrema
To investigate the relative depths of the classical potential evaluated at the EW and CB extrema, it is useful to introduce the following gauge-invariant variables [7,9]: 2 By real we mean the CP symmetry that the unbroken Z2 potential had is left unbroken. Considering a complex coefficient m 2 12 would lead to a model with explicitly broken CP, known as the Complex 2HDM [24][25][26][27][28][29][30][31][32][33]. Further promoting the Z2 symmetry to a continuous U (1) but keeping the complex soft breaking term and allowing for the possibility of flavor violation in the quark sector yields models with interesting phenomenology [34,35], but not the subject of the current paper.
In terms of these variables, then, the classical potential of Eqn. (2.1) is written as with real parameters a i and b ij = b ji . In terms of the original parameters of Eqn. (2.1), the a i and b ij are given by with all the unspecified parameters equal to zero. Collecting the x i 's, a i 's and b ij 's in vectors X, A and symmetric matrix B, respectively, we can rewrite the classical potential as 14) The values of the vector X at the EW and CB extrema are given by We can then see that the non-trivial minimisation conditions of the classical potential at the EW extremum may be expressed in terms of X EW as We can therefore see that and we note that this expression implies the following relation: Combining this expression with ∇ X V = A + BX, we obtain Therefore, the classical potential evaluated at the EW extrema is equal to: In a similar manner, we find what the non-trivial minimisation conditions imply for X CB : From these equations it is clear that from which one obtains the value of the classical potential at the CB extrema, to wit Combining the above results, we obtain: Using this result, Eqn. (2.30), Eqn. (2.29) and the fact that B = B T , we find that We thus find that Now suppose that V EW is a local minimum of the theory. It is possible to show that the mass of the charged Higgs is given by: Using this, we finally obtain that The implications of this expression are clear: if the potential has coexisting EW and CB stationary points, and the EW solution is actually a minimum, then all of its squared scalar masses will necessarily be positive; therefore, since the quantity in square brackets is guaranteed to be positive, an EW minimum implies V CB − V EW > 0 and therefore the EW minimum is deeper than the CB stationary point. Further, it can be shown that under these conditions the CB extremum is a saddle point. Thus, in refs. [9][10][11][12][13][14][15] the following tree-level theorem was established: • If the 2HDM tree-level scalar potential has an electroweak minimum, any charge breaking extremum that eventually exists will necessarily lie above that minimum.
• Further, the charge breaking extremum will necessarily be a saddle point.
We will now investigate these properties of the 2HDM vacuum structure at the loop level.

One-Loop Corrections to the Scalar Potential
The classical scalar potential of a quantum field theory is not the true scalar potential. In an interacting quantum field theory, quantum effects will induce corrections to the scalar potential. The standard way of computing the corrections to the classical scalar potential is to use the path integral and background field method. For clarity, we will explain this formalism in a quantum field theory of a single interacting scalar field. We begin by writing down the so-called generating functional of the theory in terms of a path integral over field configurations: where Φ is the scalar field, j is an external source and S[φ] is the action of the theory. By taking n-functional derivatives of Z[j] with respect to j, one can generate the n-point Green's function consisting of the connected and disconnected Feynman diagrams with nexternal propagators. We redefine the field Φ to consist of a classical component φ cl and a fluctuation field φ: Φ = φ cl + φ(x). Here φ cl is chosen to satisfy the classical equations of motion in the presence of the external source j. We can expand the action of the theory around the classical field using: where the · · · represent higher-order functional derivatives of the action (which aren't of interest to use here.) Using the expansion of the action in Z[j], where again, the · · · represent higher-order terms. Since φ cl satisfies the classical equations of motion in the presence of the source j, the term linear in φ vanishes. The quadratic term can be integrated exactly, yielding: , which is the generating functional for the connected Green's functions. To order , this is: From now on, we will drop the log(N ) term, since it a constant and will not play any role. Lastly, we define the effective action Γ[φ] through the Legendre transform of W [j]: whereφ(x) = δW [j]/δj(x). To order , the fieldφ(x) is given by (using Eqn. (3.5)): Using this relationship, we can write φ cl =φ + φ 1 , where φ 1 is of order . We can now replace φ cl in favor ofφ. Given that δS/δφ φ=φ cl = −j, we can think of j as a functional of φ cl , replacing φ cl in favor ofφ. Writing Γ 1 [φ cl ] = i /2 log det δ 2 /δΦ(x)δΦ(y) φ=φ cl , we find, to order : where we dropped terms that go like O φ 2 1 since they are O 2 . If we takeφ to be spacetime independent, the classical action evaluated atφ is simply −(V T )V 0 (φ) where V T is the space-time volume and V 0 is the tree-level scalar potential. We thus define the effective potential as: It is straightforward to evaluate the log det δ 2 S/δΦδΦ term by using the identity log(det(A)) = tr(log(A)). For a real scalar field, one has δ 2 S/δΦδΦ = + m 2 (φ) (where m 2 (φ) is the field-dependent mass computed by diagonalizing ∂ 2 V 0 (φ)/∂φ 2 ) and hence: The integral can be computed by replacing log −p 2 + m 2 with − lim α→0 ∂ ∂α (−p 2 + m 2 ) −α and using standard one-loop integral tables. The result is divergent and requires the couplings of the theory to be renormalized. Once the infinities are canceled off (using MS), the result is: with µ being the renormalization scale. It is straight forward to add in additional scalar fields, gauge bosons, and fermions. The form of effective potential to O( ) is [36]: where i runs over all the particles of the theory, s i is the spin of the particle, n i is the number of degrees of freedom of the particle, µ is the renormalization scale and M 2 i (φ) is the field dependent squared mass. The value of c i is renormalization-scheme-dependent. For MS [36], c i = 5/6 for gauge fields and 3/2 for all other particles. In principle, one needs to take into account the order correction present inφ when computing V 0 (φ). The terms of order arising from V 0 (φ) play an important role in ensuring that the effective potential is gauge-independent order-by-order in . We will discuss this further in Sec. (3.4). In the remaining subsections, we provide results for the various contributions from the particles involved in the 2HDM and discuss the expansion of the effective potential.

Scalar Contributions
For the scalar fields, the one-loop correction to the scalar potential is (3.14) In the expression above, the values of the squared masses are the eigenvalues of the second derivative of the tree-level potential, M 2 ij (Φ): with {φ i , φ j } any of the real components defined in eq. (2.2). In a general gauge, there are additional gauge dependent pieces which contribute to the scalar squared mass matrix. These gauge contribution have the effect of giving the Goldstones masses which are ξ times the corresponding massive gauge bosons (see below for the gauge masses), where ξ is the gauge-fixing parameter. However, we will chose the Landau gauge ξ = 0, where the additional gauge-dependent pieces do not contribute to the scalar mass matrix. For a general field configuration, this 8 × 8 mass matrix is extremely complicated, preventing us from giving explicit expressions to its eigenvalues. It is, however, possible to compute the scalar masses and their derivatives (which we will need) for the cases where c 2 = c 3 = c 4 = i 1 = i 2 = 0. These expressions are lengthy and we will, therefore, omit the results, having in any way developed a numerical procedure to obtain their values for our calculation.

Gauge Contributions
The field-dependent squared masses of the W and Z bosons and the photon are generated from the kinetic terms of the two Higgs doublets: Plugging in the expectation values for the Higgs doublets, the result is: where g and g are the U (1) Y and SU (2) L gauge couplings and σ a = σ 1 , σ 2 , σ 3 are the Pauli-sigma matrices. In order to compute the squared masses of the gauge fields, it is useful to organize the gauge fields into the following vector: Computing the second derivative of L gauge,mass with respect to the components of G µ generates the following 4 × 4 mass squared matrix for the gauge bosons: where we have defined: t W ≡ tan(θ W ) = g /g and the parameters x, y, z and w as: It is possible to explicitly compute the eigenvalues of the gauge mass squared matrix. In terms of the above parameters, the squared masses of the W, Z and photon are 3 : In general, one also needs to consider the effects of ghosts. Ghost fields add additional contributions to the effective potential with squared mass equal to ξ i m 2 g,i , for each of the massive gauge bosons. However, we will work in the ξ = 0 Landau gauge where the ghosts and Golstone bosons are massless.

Top quark Contribution
For simplicity, the only fermion we consider is the top quark. The contributions to the effective potential from other fermions will be significantly smaller than the top quark contribution since the top mass is almost two orders of magnitude greater than the next heaviest fermion. To compute the field-dependent squared mass of the top quark, we consider the Yukawa interactions between the Higgs doublets and the top quark: where Φ = Φ 1 or Φ 2 and the · · · terms represents Yukawa interactions involving the remaining fermions, which we ignore. Following the usual convention, we take Φ = Φ 2 , for which the top quark mass is given by: The Yukawa coupling y t will depend on the values we choose for the r 2 , i 2 , c 3 and c 4 . We define the Yukawa coupling through the EW VEV, i.e. with r 2 = v 2 and i 2 = c 3 = c 4 = 0.
The resulting Yukawa coupling is given by: Given this Yukawa coupling, the field dependent top quark mass is: Given that we will not consider the (much smaller than the top's) contributions from other quarks or leptons, the results we present here will (within that approximation) therefore be valid for the several Yukawa-types of 2HDM (Type I, II, lepton-specific and flipped [5]).

-Expansion
As is well known, the effective potential is a gauge-dependent quantity. In principle, however, physical quantities calculated from the effective potential should not be gauge dependent. In practice, however, how such physical quantities are calculated determines whether or not the gauge dependence appears. The theoretical backbone for these issues are the so-called Nielsen identities [37], which can be cast as the fact that variations of the effective potential with respect to the gauge parameter ξ are proportional to variations with respect to the field itself, The equation above holds order by order in perturbation theory and, in particular, it implies that the value of V eff at critical points, i.e. where ∂V eff /∂φ = 0, is gauge independent.
The key issue with the "brute force" minimization of the effective potential to compute physical quantities lies with the fact that truncating the perturbative expansion means that incomplete higher-order terms are, implicitly, introducing a spurious gauge dependence. The proposal of Ref. [38], known as the expansion method, consists of casting the effective potential (and its derivatives) as a series in , after "reintroducing" the in the partition function. The minimization is then carried out by an "inversion of series" method [38]. Notice that while the -expansion method was originally developed for the finite-temperature effective potential, its applicability extends (in fact, in a much more straightforward way) to the zero-temperature effective potential we are concerned with here.
The -expansion method is manifestly gauge-independent, and unlike "brute force" minimization, it does not introduce an imaginary part in the broken phase. Also, it is valid at all types of extrema, including maxima and saddle points. In practice, the method's prescription is simply to find the extrema of the tree-level potential, with the perturbative series generating the corrections order by order.
As mentioned in Eqn. (3.13), the effective potential can be expanded in terms of . To be consistent, we must also include the order contributions present in the vacuum configuration: The effects of including the O( ) contribution to the vacuum configuration is to introduce additional terms arising from the tree-level potential that contribute to the effective potential at order . The full scalar potential evaluated at φ vac , expanded to order O( 2 ) is: The extrema conditions for the full effective potential are then given by: From this expression, we can immediately interpret the meaning of φ vac is a vacuum configuration that extremizes the classical scalar potential. Eqn. (3.33) also shows us how to find the extrema of the full effective scalar potential to order . One simply needs to determine all the extrema of the tree-level potential. Then, setting the term of order in Eqn. (3.33) to zero, we can determine the one-loop correction to the classical vacuum configuration, φ (1) vec . Once the classical extrema have been determined, the minimum of the effective scalar potential will be the configuration which gives the smallest value of vac extremizes the classical scalar potential. Ref. [39] and [40] revealed IR divergences in the Landau gauge arising from massless Goldstone bosons, and argued that a resummation is necessary. However, e.g. Ref. [41] argued that the -expansion obviates the need for a resummation because the IR divergences cancel order by order in perturbation theory. Notice that this procedure was extended in Ref. [42] to the small mass limit. In the case of small, non-Goldstone masses, a resummation is necessary. If negative masses are found corresponding to a one-loop minimum, which would always be the case in our theory (because remember, a tree-level EW minimum implies that any coexisting tree-level CB extremum is a saddle point), one would additionally need to perform a resummation of the two-point correlation functions to obtain a more accurate result for the scalar masses. In the -expansion method, attempting to find a counterexample to the tree-level theorem, e.g. simultaneous minima at one-loop, would always result in at least one set of negative squared masses (since either the EW or CB vacuum would be a saddle point at tree-level) and thus one would be left with imaginary one-loop potentials. Therefore, using the -expansion method to find counterexamples to the tree-level theorem would always require a resummation to provide a sensible result. We, therefore, will instead perform a numerical minimization of the effective potential and require that the tree-level potential is convex at one-loop minima, the one-loop effective potential thus becoming free of any imaginary pieces. Any gauge dependence of the results will be residual, stemming from the truncated perturbative expansion and arising, at least at the two-loop level, at order O 2 .

Numerical Methods
In this section, we describe the procedure we employ in finding counterexamples to the tree-level theorem on EW vacuum stability against charge breaking at one-loop order. A counterexample to the tree-level theorem is obtained if we can find a set of parameters for which there exist simultaneous EW and CB minima -this, at tree-level, is impossible. Further, we will show that one-loop EW minima may have deeper CB minima and thus their stability is not guaranteed. In brief, the algorithm we use to find counterexamples is as follows: 1. Generate EW and CB VEVs for Φ 1 and Φ 2 by sampling from a uniform distribution.
2. Generate initial random guesses for all eight of the 2HDM parameters: m 2 11 , m 2 22 , m 2 12 , λ 1 , λ 2 , λ 3 , λ 4 and λ 5 by sampling from uniform distributions. These will be used later as initial "seeds" for a numerical minimization of the potential.
3. Extremize the effective potential at both the EW and CB by solving the following five non-linear root equations: We solve these equations by holding the EW and CB VEVs fixed and varying five of the 2HDM parameters. We randomly chose which of the five 2HDM parameters we use to solve these equations each time.
4. Choose a set of 50 random vacuua and perform minimizations at each to find remaining extrema.
5. Categorize all of the extrema (as minimums, maximums or saddle points) by computing the eigenvalues of the effective potential Hessian.
Below, we describe the algorithm is more detail. Note that all the code was written in Julia and is available on GitHub.
Randomly choosing VEVs: Our starting point is to choose EW and CB vacuua at which we will attempt to extremize the effective potential. We characterize the mass scale of our problem in terms of the renormalization scale µ which we set to be the SM Higgs VEV: µ = 246 GeV. Since the effective potential contains logarithms of the form log M 2 /µ 2 , we choose all of our dimensionful parameters to be of the order of the renormalization scale. We do this to avoid unwanted large logarithms which can ultimately spoil our perturbative expansion. As we did in Sec.
(2), we define the EW and CB vacuua as: In terms of the individual components of the fields Φ 1 and Φ 2 (see Eqn. (2.2)), this means that the following real components will have non-zero VEVs, with all other component fields of Eqn. (2.2)) have expectation values equal to zero. As stated above, we choose the scale of the VEVs to be on the order of the renormalization scale. That is, we set: We set v 2 1 + v 2 2 = µ 2 = (246 GeV) 2 to of course in order obtain a SM-like EW vacuum, with gauge boson and quark masses in accordance wit experiment, but with an arbitrary value of tan β = v 2 /v 1 . However, when we search for other minima by numerically minimizing the effective potential w.r.t. the fields r 2 , r 2 and c 1 (see below), we may find deeper EW minima which no longer satisfy this condition v 2 1 + v 2 2 = (246 GeV) 2 (this is a well known property of the 2HDM, already occurring at tree level). However, this condition allows us to find situations where at least there is a SM-like vacuum.
Initializing the 2HDM Parameters: In the 2HDM we consider, there are a total of 3 dimensionful mass parameters and 5 dimensionless quartic couplings (see Eqn. (2.1)): As with the vacuua, we choose the dimensionful mass parameters to be of the same order as the renormalization scale. That is, we choose: As we stated above, we make this choice to avoid generating large scalar masses which in turn could lead to large logarithms. In choosing the values of the dimensionless couplings, we keep in mind that sufficiently large couplings will result in a breakdown of perturbation theory. In practice, the breakdown occurs when dimensionless expansion parameters exceed 4π (since the perturbative expansion is in powers of (expansion parameter)/4π.) To satisfy perturbative unitarity, we keep all of the quartic couplings to be below 10. In addition to perturbative unitarity, we also wish to have a stable potential for the scalars. The tree-level conditions for stability of the scalar potential are: With these conditions and perturbative unitarity in mind, we choose the quartic couplings such that: − λ 1 λ 2 ≤λ 3 ≤ 10 + λ 1 λ 2 (4.13) −1 ≤λ 4 , λ 5 ≤ 1. (4.14) Even with these choices, it is possible to violate the stability conditions. Thus, we generate parameters according to the above prescriptions and then check if the three-level potential is bounded. If it is, we continue, otherwise, we continue to generate parameters until the potential is stabilized.
Extremize the Effective Potential: Our goal is ultimately to have minima at the EW and CB vacuua we have chosen. As a first step, we simultaneously extremize (not knowing ahead of time whether or not we are at a minimum, maximum or saddle point) the effective potential at the EW and CB vacuua. To do this, we must simultaneously solve the following five root equations: The derivatives of the effective potential are given by: Here M 2 s,i (Φ) are the eigenvalues of the scalar squared-mass matrix, M 2 g,i (Φ) are the eigenvalues of the gauge squared-mass matrix and M 2 top (Φ) is the squared top mass. Note that the factor of 3 on the log of the gauge contribution comes from the polarization of the massive gauge fields (the fact that the W boson is charged is also taken into account on the sum over the four eigenvalues of the gauge boson mass matrix of Eqn (3.18)) and the factor of 12 for the top contribution accounts for the 3 colors, 2 spins, and charge of that particle.
To solve the five root equations, we must have five independent parameters which can vary. Since we wish to fix the EW and CB vacuua, we must resort to varying five of the 2HDM parameters. To ensure that we can sample the entire parameter space, we randomly choose any five 2HDM parameters (i.e., any given five of the quadratic or quartic parameters) to vary each time we solve the extremal equations. We employ the NLsolve.jl Julia library [43] using the Trust Region method. Since we allow five of the 2HDM parameters to vary, we could potentially find solutions which make the scalar potential unstable or spoil the perturbative expansion. We thus reject solutions which for which the stability conditions are violated or solutions which have 2HDM parameter which are too large (m 2 ij > (10µ) 2 or |λ i | > 10.) As explained in Sec. (3), there are no analytical expressions for the squared scalar masses for an arbitrary vacuum configuration. They must, therefore, be computed numerically by calculating the eigenvalues of the scalar squared mass matrix. This makes computing the derivatives of the eigenvalues of the scalar mass matrix extremely difficult. To obtain those (first and second-order) derivatives, then, we employ an algorithm using forward-mode automatic differentiation through the use of dual-numbers, which we explain in App. (A). We use the FowardDiff.jl package [44], which implements a dual-number type in Julia. This allows us to simply pass dual-number types into the effective potential, and we obtain automatic derivatives without ever needing to use Eqn. (4.16) 4 .
Finding Additional Minima: As explained above, we solve the minimization conditions of the one-loop effective potential so as two obtain two different extrema. But beyond those two vacuua -one EW breaking, the other CB -the 2HDM potential may yet have other extrema. For instance, if at tree-level the CB minimum is unique (see [9,10]) other neutral minima may exist ( [13,14,16,45]). To search for any remaining minima of the effective potential, we randomly generate 50 vacuum configurations and perform a numerical minimization starting from these vacuua, verifying whether the potential may assume deeper values than the starting points. The minimization is performed using the Broyden-Fletcher-Goldfarb-Shanno algorithm provided from the Optim.jl library [46]. After performing this procedure, we sometimes find a minimum which is not one of the initial solutions found by solving the extremal equations of Eqn. (4.15). For these deeper neutral minima, it is likely that we will now have v 2 1 + v 2 2 = 246 2 GeV 2 , the so-called "panic vacuua" of refs. [47,48].
Characterizing Extrema: After we have found extrema of the effective potential, we need to determine if they are minima, maxima or saddle points. In general, an extremum of a scalar function can be characterized by computing the eigenvalues of the Hessian matrix. The Hessian matrix is the matrix consisting of all second derivatives of the function, which in our case is an 8 × 8 matrix with components: ∂ 2 V eff /∂φ i ∂φ j . The components of the Hessian matrix of the effective potential are given by: If the Hessian matrix is positive semi-definite (i.e. all the eigenvalues are greater than or equal to zero), then the extremum is a minimum (note that the zero eigenvalues signal a flat direction, which we will explain momentarily.) Similarly, if the eigenvalues of the Hessian are negative semi-definite or neither positive nor negative semi-definite, then the extremum is a maximum or saddle point, respectively. When the two Higgs doublets attain their non-trivial VEVs, the SU(3) c × SU(2) L × U(1) Y gauge group is broken. In the case of the EW VEVs, the gauge group is broken down to SU(3) c × U(1) EM and in the case of CB VEVs, the gauge group is broken down to SU(3) c . In either case, we expect there to be Goldstone bosons corresponding to each broken generator of the gauge-group. The Goldstone bosons will manifest themselves as zero eigenvalues 5 of the Hessian matrix, i.e. flat directions of the effective potential. In Tab. (1), we list the various extrema type corresponding to the eigenvalues of the effective potential Hessian. As with computing first derivatives of the effective potential, to compute the second derivatives, we use automatic differentiation. This again allows us to simply pass dual-numbers into the effective potential (in this case we pass nested dual numbers, i.e. dual numbers consisting of dual-numbers, see App. (A)) and we obtain the second derivatives without ever having to use Eqn. (4.17).
Notice that the second derivatives of the one-loop effective potential provide the one-loop ∂ 2 V eff /∂φ i ∂φ j Eigenvalues Extrema Type 3 zero, 5 positive EW minimum 3 zero, 5 negative EW maximum 3 zero, 5 positive and negative EW saddle 4 zero, 4 positive CB minimum 4 zero, 4 negative CB maximum 4 zero, 4 positive and negative CB saddle Table 1: Characterization of the extrema of the 2HDM effective potential. squared scalar masses computed at zero external momentum. For massive scalars they are therefore an approximation to the exact result, but for the massless Goldstones -which must be computed at precisely zero external momentum -they yield the exact result. Obtaining the correct number of massless Goldstones for either the EW or CB extrema is a powerful check of our calculations.

Results
In this section, we describe the results of running the algorithm described in the previous section to find counter-examples to the tree-level theorem described in Sec. (2). To wit, our purpose is to investigate whether at the one-loop level an EW minimum is guaranteed to be stable against charge breakingi.e., whether still at one-loop there is no deeper CB extremum. Further, we will verify whether at one-loop the existence of an EW minimum also implies that any CB extremum must need be a saddle point. All computations were run on a 2015 Mac Book Pro using 8 threads. We developed all of the code for this algorithm using the Julia language, using various well-developed Julia packages. For example, we use the ForwardDiff.jl [44] package for automatic differentiation, NLsolve.jl [43] for solving the root equations of Eqn. (3.33) and Optim.jl [46] for performing minimizations. All the code developed for this project can be viewed/downloaded on GitHub. For more details, the interested reader may e-mail the authors.
To consider a set of vacuua (EW and CB) and 2HDM parameters to yield a counterexample to the tree-level theorem, we set various requirements. First, to be a counterexample to the tree-level theorem, we must have a minimum of the effective potential at an EW and CB vacuum. We consider a vacuum to be an extremum of the effective potential if the infinity-norm of the gradient is less than 10 −5 (although in many cases we obtain much higher accuracy.) We categorize the extremal type (minimum, maximum or saddle) of a vacuum using the conditions in Tab. (1). In particular, we consider an EW vacuum a minimum if the Hessian of the effective potential evaluated at the vacuum contains three zero masses (Goldstone bosons corresponding to the breaking of SU(2) L ×U(1) Y → U(1) EM ) and five positive masses. In the case of a CB vacuum, we require four zero masses (the additional zero mass due to the explicit breaking of U(1) EM ) and four positive masses. In addition, we require that the tree-level potential is bounded from below (following the conditions of Eqns. (4.9-4.10).) We also require that the values of the 2HDM be constrained to be of natural order: |m 2 ij | < (10µ) 2 , |λ i | < 10 where µ is the renormalization scale. We make these requirements to preserve our perturbative expansion and avoid generating large masses which could result in large logarithms.
After running our algorithms for roughly 24 hours, we found ∼ 3000 sets of parameters which yield simultaneous one-loop EW and CB minimathis is the first demonstration that the tree-level vacuum stability theorem is no longer valid at one-loop. Out of these 3000 sets, for ∼ 1000 of then, the global minimum of the one-loop effective potential was the CB vacuum; the remaining ∼ 2000 had the EW vacuum as the global minimum. To get a sense of how common the sets of parameters yielding counter-examples were, we also recorded those sets of parameters for which there was only an EW minimum (no CB) and for which there was only a CB minimum (no EW). The former yielded ∼ 54000 sets of parameters, while the latter ∼ 17000. Thus, we can see that parameters which yield both a CB and an EW minimum are roughly 5% of those which yield a single minimum -thus even at one-loop, we can expect that the exclusion of regions of parameter space due to CB vacuum instability will be rare. Furthermore, only 4 out of the 1000 points which yielded a deeper CB minimum have positive tree-level masses and 10 for the case where the EW was deeper (the remaining contained at least one negative tree-level mass from either the CB or EW vacuum.) We should stress, however, that our purpose is not to perform a thorough scan of the 2HDM parameter space to find charge breaking bounds of the model, but rather to prove that the tree-level vacuum stability theorem no longer holds. As mentioned in Sec. (3.4), for parameters with negative tree-level masses which result in one-loop minima, one likely needs to perform a resummation to obtain a sensible result (i.e. one that doesn't exhibit an apparent instability -imaginary part of the effective potential), which we have not done. But having performed the scan over the model's parameters such that the tree-level masses at the one-loop minima were always positive, the issue of a complex one-loop effective potential is no longer an issue that should worry us.  Table 2: Parameter values. Quadratic parameters in GeV 2 , VEVs in GeV and potential values in GeV 4 . All values have been rounded for readability.
We provide two sets of parameters yielding counterexamples to the tree-level theorem in Tab. (2): one for the case where the CB vacuum is deeper than the EW one (left column) and one where the EW minimum is deeper than the CB one (right column.) 6 We reiterate that both of these points are convex at tree-level, meaning all of the squared scalar masses at tree-level are positive at the one-loop extrema, yielding no complex contributions to the effective potential. To better visualize the behaviour of the potential at both extrema for both sets of parameters, consider Fig. (1). Of course, it is impossible to visualize the full, 8-dimensional potential at these points, but to give some sense of what it looks like in the vicinity of the one-loop vacuua, we resort to one-and two-dimensional "slices". Thus, in Fig. (1), we display the one-loop and tree-level potential evaluated at each vacuum and along a line linearly interpolating the EW and CB vacuua for the parameters/vacuua given in Tab. (1). This means we are evaluating the potential along values of the fields given by, for each component of the doublets, φ(t) = (1 − t) φ EW + t φ CB . Thus, at t = 0, the potential is being evaluated at the EW vacuum and at t = 1 at the CB vacuum. Fig. (1) show that the potential always has minima at the EW and CB extrema, both at tree and one-loop level. This is, however, deceiving -at tree-level, the vacuum stability theorem states that if there is an EW minimum, any CB extrema will be a saddle point. However, the tree-level potential is not in fact at minima for both the EW and CB one-  Figure 1: One dimensional slices of the effective scalar potential. The horizontal-axis represents vacuum configurations interpolating between the CB and EW vacuua. i.e. we interpolate between φ(t) = (1 − t)φ EW + tφ CB . Hence, at t = 0, φ(t = 0) = φ EW and at . The values of the VEVs and parameters are given in Tab. (2).
loop-vacuua. For both points given, the parameters give no solutions for the tree-level minimization conditions, and the one-loop EW vacuum is near the global tree-level vacuum but the one-loop CB is simply at some convex point at tree-level (but not an extremum). It would be easy to see that along some other direction(s) in field space the seeming tree-level extrema would not be minima at all. What is however clear from Fig. (1) is, as soon as we realize that at one-loop both the EW and CB extra are minima, the tree-level vacuum stability theorem is once again violated at the one-loop levelit is possible, at one-loop, to obtain a potential with an electroweak breaking minimum, which also possesses a deeper charge breaking minimum. Thus the absolute stability of 2HDM EW minima found at treelevel is broken by radiative corrections -the quantum mechanical effects on the effective potential can change the vacuum properties of the model. To further illustrate the behaviour of the 2HDM potential close to these extrema consider Figs. (2) and (3). There we display a two-dimensional slice of the effective and tree-level potentials. In these figures, the horizontal axis is identical to that of the onedimensional plots of Fig. (1) -that is, a line interpolating between both one-loop minima. The vertical axis in Figs. (2) and (3) represents variation along a direction s in r 1 − r 2 − c 1 space which is orthogonal to the line interpolating the EW and CB vacuua. These figures give us a slightly more convincing visualization of the minimization at the EW and CB vacuua, and show the distortion induced upon the tree-level potential by the loop corrections. Having said that, they are nonetheless incomplete images of the full 8-dimensional picture and can not illustrate, for instance, the conversion between tree-level saddle points and one-loop extrema that the violation of the tree-level vacuum theorem implies.
We have fixed the renormalization scale µ to 246 GeV, and the procedure we fol-   Figure 2: A two-dimensional slice of the effective potential (left) and the corresponding tree-level potential (right) for the case where V eff (φ CB ) < V eff (φ EW ). The horizontal axis is identical to that of Fig. (1). The vertical axis is an line in r 1 − r 2 − c 1 space orthogonal to the t-axis, i.e. orthogonal to a line connecting φ EW and φ CB . The scale of s vertical axis is identical to the scale of the horizontal axis, i.e. distance in field space from (t = 0, s = 0) and (t = 1, s = 0) is identical to the distance in field space between (t = 0, s = 0) and (t = 0, s = 1) (both these distances are the distances between φ EW and φ CB .) lowed should, obviously, not depend on that choice. To verify that the results we obtained are indeed not dependent on a particular choice of µ, we took the parameters given in Tab. (2) and evolved them according to their RG equations (see the appendix of Ref. [49] for explicit expressions for the RG equations for the THDM parameters.) We use the DifferentialEquations.jl [50] package to perform the RG evolution of the parameters from µ = 246 GeV to µ = 400 GeV. At all renormalization scales between 246 − 400 GeV, we re-minimize the effective potential starting from both the EW and CB vacuua, determining the new VEVs at each minimum for the new values of the parameters of the potential at the new scales. We then compute the value of the one-loop effective potential at each minimum, which we show as a function of the renormalization scale in Fig. (4), for both sets of parameters given in Tab. (2). As we can see, the separation of the EW and CB vacuua is preserved as a function of the renormalization scale. Also, we can see that the difference of the values of the potentials is nearly a constant, which is a consequence of the fact that the effective potential at the minimum is RG independent. The values of the effective potential These plots demonstrate that our results are independent of the particular choice of renormalization scale that we chose.
only change due to us not including the RG evolution of the field-independent piece of the effective potential (which is the same for both the curves). To better understand this point consider the discussion in [51,52]: the one-loop effective potential depends on a set of parameters λ i , fields φ j and renormalization scale µ, and it may be written generically as where the field-independent term Ω is the same for any extremum of the potential. The crucial insight to understand the behavior shown in Fig. (4) is that, unlike what one usually thinks, the sum V 0 + V 1 of the tree-level and one-loop contributions to the potential, is not RG independent. Rather, the independence of the renormalization scale on the full effective potential, dV /dµ = 0, is accomplished at the one-loop level by "compensating" the µ dependence on V 0 + V 1 with that of the Ω term [53]. But since Ω does not depend on the value of the fields it will not change between the EW and CB minima, and as such it is trivial to obtain d(V 0 + V 1 ) EW /dµ = d(V 0 + V 1 ) CB /dµ -meaning, one expects that by varying the value of the renormalization scale, the value of the potential at the EW minimum evolves "parallel" to that of the CB minimum, and that is exactly the behavior one witnesses in Fig. (4). Thus the conclusion is that indeed our one-loop result is not an artifact of a specious choice of renormalization scale, but rather it is independent of the value of µ. However, we emphasize that these parameter sets exemplify the best case scenario obtained in our numerical calculations, boasting nearly perfect RG evolution. Not all parameter sets we found behave as well. In particular, we find some parameter sets for which the RG curves cross. Crossing of the RG curves signals that 2-loop corrections to the effective potential and RG equations are likely important for those particular sets of parameters. Before concluding, it is worth mentioning the consequences of these results. At treelevel, it was clear that, since it is impossible to have simultaneous EW and CB minima, an EW vacuum would be stable against the possibility of charge breaking, i.e. no tunneling could occur that would spoil the residual U(1) EM gauge symmetry and disastrously give the photon a mass. This is no longer the case when one considers the quantum corrections to the classical potential. That is, simultaneous EW and CB minima can exist at oneloop. This implies that, if we lived in a EW vacuum of the one-loop effective potential in a scenario where there is an additional, deeper CB vacuum, it would be possible to tunnel to the CB vacuum. The decay rate of the EW vacuum would be highly suppressed since the simultaneous minima are only realized at one-loop. In particular, we would then expect the decay rate to be a two-loop effect and of the order O 2 .

Conclusions
In this work, we have analyzed the vacuum structure of the one-loop effective potential of the 2HDM. At tree-level, the 2HDM scalar potential is found to have a remarkable stabilityany minimum which breaks the ordinary electroweak symmetries and thus preserves charge conservation (and furthermore, also preserves CP) is guaranteed to be stable against the possibility of charge breaking vacuua -meaning, any CB extremum that eventually might coexist with that minimum is guaranteed to lie above it, and furthermore to be a saddle point. This theorem was found in 2004 via analytical calculations with the 2HDM potential, along with a series of other remarkable results concerning the model's vacuum structure [9,10,13,14].
The first hint that these vacuum stability theorems might not hold at one-loop was obtained analyzing the coexistence of neutral minima in a version of the 2HDM, the Inert Model. Comparing tree-level minima with one-loop ones, using the formalism of the effective potential, it was possible to show that the loop corrections might indeed change the nature of the vacuum -for certain choices of parameters, a minimum which at tree-level was global would become a local one at one-loop [20]. It then became clear that an analysis of the one-loop potential was required to ascertain whether the stability of EW minima against deeper CB vacuua remained a valid conclusion. This work shows that the theorem does not hold at one loop.
We have indeed obtained, through extensive numerical scans of the parameters of the model, many cases where an EW minimum of the one-loop effective potential can coexist with a CB minimum -this is a first violation of the tree-level theorem, which stated that an EW minimum implied necessarily CB saddle points. We have also determined that one-loop EW minima can coexist with deeper CB minima -and hence the tree-level stability against CB of the 2HDM no longer holds at one-loop. The conclusion one must draw from these results is that quantum corrections to the potential may change the vacuum structure of said potential. Conclusions drawn at tree level for which kind of minimum is the global one, and whether it is stable, may well not survive a higher-order calculation. And this, in fact, perhaps should not surprise us -after all, this is indeed what one already obtained in the case of the Coleman-Weinberg potential [54].
Our calculations were performed using tried-and-true computational algorithms and numerical minimization routines which are widely available, and we offer two examples of parameter sets to be checked by interested readers. Issues of gauge dependence of the effective potential should not affect the validity of the conclusions drawn here since we are fundamentally comparing the value of the effective potential in different minima. Though we only included the contribution of the top quark, clearly the results would not qualitatively mutate with the inclusion of further fermions. And the calculations underwent a rigorous check via the computation, at each EW and CB extremum, of the respective one-loop Hessian matrices. That check had a twofold purpose: to verify the nature of any given extremum, so that we could be certain when claiming to have found minima and to verify whether the correct number of Goldstone bosons was found -three for any EW vacuum, four for a CB one. Further, a verification of the independence of our results from the value of the renormalization scale µ we chose was undertaken -an RG evolution of the parameters of the potential in an interval of values of µ was performed, followed by a re-minimization of the potential to obtain the values of the VEVs at each new scale. The comparison of the values of the potential showed that the relative depth of the minima remained unchanged with the renormalization scale, and thus our conclusions are RG stable.
Should this mean that we are witnessing a breakdown in perturbation theory, wherein higher-order corrections invalidate our calculations? Hardly -the RG evolution performed showed us that perturbation theory is working as one would expect. The results of [20] should further illuminate our conclusions -what was found there was that loop corrections could change tree-level expectations for the nature of the vacuua by "swapping" global and local minima, but that this could only occur if both minima were nearly degenerate. Thus, one concludes, at least for the results of [20], the loop corrections are small and acceptable perturbations that "flip" the system between two states of nearly degenerate energies. Likewise, the interpretation of the results we present in the current work points to perturbation theory still holding: a vast numerical scan of the model's parameters only yields counterexamples to the tree-level theorem for a small subset of the parameter space. Also, the tree-level result was strongly dependent on the specific form of the potential; of its derivatives; of the scalar squared masses. At one-loop something remarkable occursthe vacuum is determined, not only by the scalar sector but by all sectors of the theory, gauge and Yukawa included. It is therefore unsurprising that different statements can be made.
To conclude, the 2HDM electroweak vacuua is not guaranteed to be stable against charge breaking vacuua -there may well be, for certain regions of parameter space, deeper CB minima below an EW one, and it may well happen that the tunneling time to the deeper minimum is smaller than the age of the universe. Though we expect this situation to be rare, this work raises the necessity to perform a wide reassessment of bounds imposed upon the parameters of the 2HDM, by fully analyzing the one-loop vacuum structure of the model. The task is not an easy one, for the one-loop effective potential is very complex and unwieldy, especially at CB vacuua. Finally, two comments-first, we have used the results from [9,10] concerning the simplest form one could take for CB vacuua; but those results stem from an analysis of the tree-level potential, so they too might change when considering a one-loop calculation. Second, in [9,10] the tree-level theorems deduced concerned the stability of EW vacuua against, not only CB, but also minima with spontaneous CP violation. As for the case of CB, the conclusion therein obtained was that any EW minimum cannot have a deeper CP breaking extremum, and any such extremum is found to lie above the EW minimum and be a saddle point. Given that we have shown that at one-loop an EW could coexist with a deeper CB minimum, there are strong reasons to believe that the same will apply to coexistence with CP minima. differences in which the derivative is approximated using (with forward finite differences): This method suffers from many issues. First off, to get a good approximation of the derivative, one would like to make as small as possible. However, due to the finite precision of machine numbers, as becomes sufficiently small, round-off errors will seep into the calculation and the error in the approximation will increase [55]. Thus, there is a given value of for which finite-differences yields the smallest error and one can do no better. Another method for evaluating derivatives is to use the complex-step method [56], in which the derivative is approximated using: This method doesn't suffer from the round-off errors that arise from finite differencing. can be taken arbitrarily small. However, the complex step method requires one to only use real numbers (if one mixes complex derivatives with the complex step method, the results will be non-sense.) A slightly more complicated, but robust method of numerically computing derivatives in forward-mode automatic differentiation [55]. The core idea of forward-mode automatic differentiation is the concept of dual-numbers. A dual-number is defined similarly to infinitesimals: where has the property that 2 = 0. The algebra of dual numbers is defined as follows: If we set d = x + (i.e. set b = 1), then we find f (x + ) = f (x) + f (x). We thus obtain f (x) and f (x) by evaluating f at the dual number x + . Using dual-numbers provides us with a method of computing exact derivatives (up to machine precision.) Dual-numbers can also be used to compute higher-order derivatives by nesting dual-numbers: i.e. have dualnumbers of dual-numbers. For example, if d = a + 1 b with a = a 1 + a 2 2 and b = b 1 + b 2 2 , with 2 1 = 2 2 = 0 and 1 2 = 0, then we find the following: = f (a 1 ) + 2 a 2 f (a 1 ) + 1 (b 1 + 2 b 2 ) · (f (a 1 + 2 f (a 1 ))) (A.8) = f (a 1 ) + b 1 1 f (a 1 ) + a 2 2 f (a 1 ) + b 2 1 2 f (a 1 ) + a 2 b 1 1 2 f (a 1 ) (A.9) If we set a = x + 2 and b = 1 + 0 2 , then we obtain: f ((x + 2 ) + 1 (1 + 0 2 )) = f (x) + 1 f (x) + 2 f (x) + 1 2 f (x) (A.10) Hence, the 1 component of the number gives the first derivative of f and the 1 2 -component gives the second derivative of f . If we continue nesting dual-numbers, we can compute arbitrary derivatives of f (x). Given the power of template meta-programming and multiple-dispatch built into Julia, it is an easy task to implement dual-numbers numerically. Below we provide code snippets of how this is done (note that this is not what we use, instead we use ForwardDiff.jl a well-developed Julia package.) The basic idea of implementing dual-numbers is to define a new type, which we call Dual. We then overload all necessary operations that we need, i.e., addition, subtraction, multiplication, division and any other functions we wish to use with dual-numbers. Our type Dual contains two attributes: val (the real component of the dual-number) and eps (the infinitesimal part): where the second component of the dual number is: ∂ ∂x (xy) = y = 2. Another example would be to take the sine of a dual number: julia> sin(d1) Dual{Float64}(0.8415, 0.5403) which we notice has sin(1) in the first component and cos(0) in the second component. Once basic operations like the above have been defined, one can then chain together very complicated functions and easily obtain their derivatives. Additionally, we can easily take second derivatives as well by nesting the dual numbers. If we define a cos overload, we can then take the second derivative of the sine function: Base.cos(z::Dual{T}) where T<:Real = Dual{T}(cos(z.val), -sin(z.eps)) julia> d3 = Dual{Dual{Float64}} (Dual{Float64}(1., 0.), Dual{Float64}(0., 1.)) julia> sin(d3) Dual{Dual{Float64}}(Dual{Float64}(0.8415, 1.0), Dual{Float64}(1.0, -0.8415)) Here, the second component of the first dual is d/dx sin(x) = cos(1), the first component of the second dual is the same and the second component of the second dual is the second derivative of sin at x = 1.