New Tools for Dissecting the General 2HDM

Two Higgs doublet models (2HDM) provide the low energy effective field theory (EFT) description in many well motivated extensions of the Standard Model. It is therefore relevant to study their properties, as well as the theoretical constraints on these models. In this article we concentrate on three relevant requirements for the validity of the 2HDM framework, namely the perturbative unitarity bounds, the bounded from below constraints, and the vacuum stability constraints. In this study, we concentrate on the most general renormalizable version of the 2HDM -- without imposing any parity symmetry, which may be violated in many UV extensions. We derive novel analytical expressions that generalize those previously obtained in more restrictive scenarios to the most general case. We also discuss the phenomenological implications of these bounds, focusing on CP violation.

Two Higgs doublet models (2HDM) provide the low energy effective theory (EFT) description in many well motivated extensions of the Standard Model. It is therefore relevant to study their properties, as well as the theoretical constraints on these models. In this article we concentrate on three relevant requirements for the validity of the 2HDM framework, namely the perturbative unitarity bounds, the bounded from below constraints, and the vacuum stability constraints. In this study, we concentrate on the most general renormalizable version of the 2HDM -without imposing any parity symmetry, which may be violated in many UV extensions. We derive novel analytical expressions that generalize those previously obtained in more restrictive scenarios to the most general case. We also discuss the phenomenological implications of these bounds, focusing on CP violation.

Introduction
The Standard Model (SM) [1] relies on the introduction of a Higgs doublet, whose vacuum expectation value breaks the electroweak symmetry [2][3][4]. This mechanism generates masses for the weak gauge bosons and charged fermions, as well as potentially the neutrinos (although there may be other mass sources for the latter). The Standard Model Higgs sector is the simplest way of implementing the Higgs mechanism for generating the masses of the known elementary particles. However, it is not the only possibility, and may be easily extended to the case of more than one Higgs doublet without violating any of the important properties of the SM. Moreover, one of the simplest of these extensions -two Higgs doublet models (2HDMs) [5] -appears as a low-energy effective theory of many well motivated extensions of the Standard Model (e.g. those based on supersymmetry [6][7][8][9][10][11][12][13][14][15] or little Higgs [16]). Two Higgs doublet models may differ in the mechanism of generation of fermion masses. If both Higgs doublets couple to fermions of a given charge, their couplings will be associated to two different, complex sets of Yukawa couplings, which would form two different matrices in flavor space. The fermion mass matrices would be the sum of these, each multiplied by the corresponding Higgs vacuum expectation value. So diagonalization of the fermion mass matrices does not lead to the diagonalization of the fermion Yukawa matrices. Such theories are then associated with large flavor violating couplings of the Higgs bosons at low energies -a situation which is experimentally strongly disfavored. Hence, it is usually assumed that each charged fermion species couples only to one of the two Higgs doublets. In most works related to 2HDM, this is accomplished by implementing a suitable Z 2 symmetry. The different possible charge assignments for this Z 2 symmetry then fix the Higgs-fermion coupling choices and define different types of 2HDMs.
This Z 2 symmetry not only fixes the Higgs-fermion couplings but also forbids certain terms in the Higgs potential that are far less problematic with respect to flavor violation. As a starting point for an investigation of the phenomenological implications of these terms, we will in this work discuss the theoretical bounds on the boson sector of the theory (without any need to specify the nature of the Higgs-fermion couplings). We will concentrate on the constraints that come from the perturbative unitarity of the theory, the stability of the physical vacuum, and the requirement that the effective potential is bounded from below. Existing works [5,[17][18][19][20][21][22][23][24][25][26][27][28][29] focus either on the Z 2 -symmetric case or only provide a numerical procedure to test these constraints in the general 2HDM (see Ref. [30] for a recent work on analytic conditions for boundedness-from-below). We will go beyond current studies by deriving analytic bounds that apply to the most general, renormalizable realization of 2HDMs. Our conditions will be given in terms of the mass parameters and dimensionless couplings of the 2HDM tree-level potential. At the quantum level, however, these parameters are scale dependent; although we will refrain from doing so here, one can apply these conditions at arbitrarily high energy scales by using the renormalization group evolution of these parameters.
Our article is organized as follows. In Section 2 we introduce the scalar sector of the most general 2HDM that defines the framework for most of the work presented in this article. Section 3 reviews three theorems from linear algebra which will allow us to derive analytic bounds in the coming sections. In Section 4, we concentrate on the requirement of perturbative unitarity. Section 5 presents the bounds coming from the requirement that the tree-level potential be bounded from below. In Section 6, we discuss the vacuum stability. Finally, we reserve Section 7 for a brief analysis of the phenomenological implications (focusing on CP violation) and Section 8 for our conclusions. A Table listing the most relevant findings of our work may be found at the beginning of the Conclusions.

The general 2HDM
As emphasized above, we focus on the scalar sector of the theory. In general, gauge invariance implies that the potential can only include bilinear and quartic terms. Each of the three bilinear terms has a corresponding mass parameter, of which two (m 2 11 and m 2 22 ) are real while the third, m 2 12 , is associated with a bilinear mixing of both Higgs doublets and may be complex.
Regarding the quartic couplings in the scalar potential, the two associated with self interactions of each of the Higgs fields, λ 1 and λ 2 , must be real and, due to vacuum stability, positive. There are two couplings associated with Hermitian combinations of the Higgs fields, λ 3 and λ 4 , which must be real, though not necessarily positive. The coupling λ 5 is associated with the square of the gauge invariant bilinear of both Higgs fields, and it may therefore be complex. The couplings λ 6 and λ 7 are associated with the product of Hermitian bilinears of each of the Higgs fields with the gauge invariant bilinear of the two Higgs fields, and, as with λ 5 , they may be complex. The most general scalar potential for a complex 2HDM is therefore: with Φ 1,2 = (Φ + 1,2 , Φ 0 1,2 ) T being complex SU (2) doublets with hypercharge +1. One way to prevent Higgs-induced flavor violation in the fermion sector is to introduce a Z 2 parity symmetry under which each charged fermion species transforms as even or odd. The Higgs doublets are assigned opposite parities and couple only to those charged fermions that carry their own parity. In such a scenario, the terms accompanying the couplings λ 6 and λ 7 would violate parity symmetry and hence should vanish. The mass parameter m 2 12 is also odd under the parity symmetry but induces only a soft breaking of this symmetry, which does not affect the ultraviolet properties of the theory. Thus it is usually allowed.
There are alternative ways of suppressing flavor violating couplings of the Higgs to fermions which do not rely on a simple parity symmetry and hence allow for the presence of λ 6 and λ 7 terms. One example is the flavor-aligned 2HDM [31]. Alternatively, one can impose a parity symmetry in the ultraviolet but allow the effective low energy field theory to be affected by operators generated by the decoupling of a sector where this symmetry is broken softly by dimensionful couplings which do not respect the parity symmetry properties. One example of such a theory is the NMSSM in the presence of heavy singlets, as discussed in Ref. [32]. In this case, the presence of the couplings λ 6 and λ 7 is essential to allow for the alignment of the light Higgs boson with a SM-like Higgs, leading to a good agreement with precision Higgs physics even in the case of large Higgs self couplings.
So we see that it is not necessary to restrict to the Z 2 -symmetric 2HDM with vanishing λ 6 and λ 7 to avoid Higgs induced flavor violation in the fermion sector. 1 Further, it is phenomenologically interesting to study the 2HDM in full generality with these terms present. One consequence would be the possibility of having charge-parity (CP) violation in the bosonic sector. Indeed, to keep good agreement with Higgs precision data [33,34], one is normally interested in studying 2HDM in (or close to) the exact alignment limit -the limit in which one of the neutral scalars carries the full vacuum expectation value and has SM-like tree-level couplings [35][36][37][38][39][40]. If one imposes exact alignment in the Z 2 -symmetric 2HDM, however, CP is necessarily conserved, as will be explained in detail in Section 7. In the full 2HDM, on the other hand, one can have CP violation whilst remaining in exact alignment thanks to the presence of λ 6 and λ 7 terms. This CP violation could manifest in the neutral scalar mass eigenstates as well as bosonic couplings (see also Ref. [41]), providing many potential experimental signatures. With this motivation in mind, we keep λ 6 and λ 7 non-zero throughout this work.

Methods for bounding matrix eigenvalues
In this work, much of the analysis of perturbative unitarity and vacuum stability involves placing bounds on matrix eigenvalues. In the most general 2HDM, analytic expressions for these constraints are either very complicated or simply can not be formulated. To obtain some analytic insight, we derive conditions which are either necessary or sufficient. Their derivation is based on three linear algebra theorems which we briefly review here.

Frobenius norm
One may derive a bound on the magnitude of the eigenvalues of a matrix using the matrix norm. The following definition and theorem are needed: Theorem: The magnitude of the eigenvalues e i of a square matrix A are bounded from above by the matrix norm: where a matrix norm is defined as: Definition: Given two m × n matrices A and B, the matrix norm ||A|| satisfies the following properties: • ||A|| ≥ 0, The above theorem holds for any choice of matrix norm, and thus one may employ the Frobenius norm [42], ||A|| = Tr(A † A), to find the following result: This bound on the eigenvalues will be used to derive sufficient bounds in the following sections.

Gershgorin disk theorem
We will use the Gershgorin disk theorem [43] in upcoming sections to derive sufficient conditions for perturbative unitarity and vacuum stability of the 2HDM potential. The theorem is typically used to constrain the spectra of complex square matrices. The basic idea is that one identifies each of the diagonal elements with a point in the complex plane and then constructs a disk around this central point, with the radius given by the sum of the magnitudes of the other n − 1 entries of the corresponding row. 2 The theorem says that all eigenvalues must lie within the union of these disks. Formally, we have the following definition and theorem: Definition: Let A be a complex n × n matrix with entries a ij , and let R i be the sum of the magnitudes of the non-diagonal entries of the i th row, R i = j =i |a ij |. Then the Gershgorin disk D(a ii , R i ) is defined as the closed disk in the complex plane centered on a ii with radius R i .
Theorem: Every eigenvalue of A lies within at least one such Gershgorin disk D(a ii , R i ).
This theorem can be used to derive an upper bound on the magnitude of the eigenvalues of a matrix. We will use this technique below when discussing perturbative unitarity and vacuum stability. Since all the matrices we will consider in the subsequent sections on perturbative unitarity and boundedness from below are Hermitian matrices, each eigenvalue will lie within the intervals formed by the intersection of the Gershgorin disks with the real axis.
We shall proceed in the following manner: We will first construct the intervals containing the eigenvalues of each matrix A. For each interval, the rightmost and leftmost endpoints x ± i will be given by the sum and difference, respectively, of the center and the radius, We then identify which x ± i extends furthest in the positive or negative direction. We know that every eigenvalue e k must lie within the endpoints of the largest possible total interval, This may be rephrased into an upper bound on |e k | as: where the left-hand side of Eq. (5) represents the absolute value of any given eigenvalue e k and the right-hand side represents the maximum value of j |a ij | over all rows i of the matrix A. In fact, this condensed statement of Gershgorin circle theorem is an application of the matrix norm theorem, employing the norm ||A|| = max i ( j |a ij |).

Principal minors
In order to derive necessary conditions, one may employ Sylvester's criterion in a clever way, as proposed in Ref. [44]. Sylvester's criterion involves the principal minors D k of a matrix, where D k is the determinant of the upper-left k × k sub-matrix. The statement of Sylvester's criterion is the following: To see this, consider applying a unitary transformation which diagonalizes M to the matrix c1 ± diag(M ). Then for symmetric or Hermitian matrices, the statement that c1 ± M is positive definite becomes a statement on the relative values of e i and c. In this manner, the application of Sylvester's criterion to these specific matrices allows one to put an upper bound on the magnitude of the eigenvalues without diagonalizing the matrix M . Note that for the absolute value |e i | to be bounded by c, we require the use of both c1 ± M matrices. On the other hand, if one has only an upper bound on e i , i.e. e i < c, as we will have in the case of vacuum stability, then one only requires the principal minors of the matrix c1 − M to be positive.
We note that the use specifically of the upper-left sub-matrices in Sylvester's criterion is an arbitrary choice, and basis-dependent. One could instead consider the lower-right sub-matrices, or any matrices along the diagonal. As such, it is possible to derive further conditions using this criteria by further considering, for example, the upper-left, lowerright, and central 2 × 2 sub-matrices of a 4 × 4 matrix. We will do so in later analyses to strengthen the lower-k bounds.
This use of sub-determinants has been proposed in Ref. [44] as a method to increase the efficiency of parameter scans in models with large scattering matrices. For such theories (e.g. the model with N Higgs doublets, NHDM, considered in Ref. [44] for higher N ), the numerical calculation of the scattering matrix eigenvalues is computationally expensive. We note that the use of the Gershgorin disk theorem proposed in Section 3.2 provides an additional complementary method to speed-up parameter scans.

Perturbative unitarity
Tree-level constraints for perturbative unitarity in the most general 2HDM have already been investigated in the literature [22,29,45]. However, for a non-zero λ 6 and λ 7 , no exact analytic conditions have been obtained yet. Here, we will first review the existing literature and then derive analytic expressions for the case of non-vanishing λ 6 and λ 7 .

Numerical bound
Perturbative unitarity is usually imposed by demanding that the eigenvalues e i of the scalar scattering matrix at high energy be less than the unitarity limit, |e i | < 8π. Thus to derive the constraints on the quartic couplings, one must construct the scattering matrix for all physical scalar states.
We are interested in all processes AB → CD, where the fields A...D represent any combination of the physical 4 scalars (H 1 , H 2 , H 3 , H ± , W ± L , Z L ). The interactions and hence S-matrix take a complicated form in terms of the physical states. However since we are only interested in the eigenvalues of the S-matrix, we may choose any basis related to the physical basis by a unitary transformation. The derivation is simplest in the basis of the original Higgs fields (w ± i , h i , z i ), appearing as with v = i v 2 i = 246 GeV. Out of these fields, we can construct 14 neutral two-body states: and |h 1 z 2 . By constructing states which are linear combinations with definite hypercharge and total weak isospin, denoted by (Y, I), and grouping the ones with the same set of quantum numbers, the matrix of S-wave amplitudes a 0 takes a block diagonal form (for more details see Refs. [22,29,45]). For the neutral scattering channels, this is: where the subscript of each submatrix denotes the quantum numbers (Y, I) of the corresponding states. The entries are: For the 8 singly-charged two-body states w + where the new entry X (1,0) is just the one-dimensional eigenvalue: Finally, the 3×3 S-matrix for the three doubly-charged 2-body states We impose perturbative unitarity by demanding that the eigenvalues of the scattering matrix are smaller than 8π implying that |a 0 | < 1 2 . Indeed, the eigenvalues of the submatrices X (0,0) , X (0,1) , X (1,0) , and X (1,1) , which we denote as e i , must all satisfy Obtaining analytic expressions for the eigenvalues requires solving cubic and quartic equations, and the result is complicated and not very useful. Given a choice of input parameters, however, it is easy to check this condition numerically. Assuming all λ 1...7 to be equal, the strongest constraint arises from the 4×4 matrix X 00 . If we set all λ i ≡ λ and solve for the eigenvalues, we find the bound: This value is an order of magnitude smaller than 8π, implying that if all quartic couplings are sizable (i.e. of O(1)), perturbative unitarity may be lost even at values of the couplings much smaller than 4π, which is a bound often encountered in the literature to ensure perturbativity. We investigate the validity of such an upper bound further in Fig. 1. For this figure, we randomly choose each λ i within the range |λ i | < λ max and then show the fraction of test points which pass the numerical unitarity constraint, as a function of λ max . For λ max π almost all points survive the perturbative unitarity constraint. For larger λ max values the survival rate quickly drops to almost zero for λ max 4π. This highlights again that the simplified perturbativity bound of |λ i | < 4π, which is often encountered in the literature, is too loose. Based on the results in Fig. 1, a better choice of bound might be |λ i | π or |λ i | 3π/2.

A necessary condition for perturbative unitarity
To gain some intuition for the perturbative unitarity constraint, we now turn to derive some simplified analytic conditions which are either necessary or sufficient, though not both. In this section we will focus on the former, which can be used to quickly rule out invalid parameter sets which violate perturbative unitarity. One can derive necessary conditions by invoking the method of principal minors, which is reviewed in Section 3.3 and can be used to give an upper bound on the maximal value of the eigenvalues of a Hermitian matrix. Since the scattering matrices are Hermitian, demanding the eigenvalues to be bounded as |e i | < 8π, as required by perturbative unitarity, amounts to requiring: for k = 1, 2, 3, 4. Since satisfying both criteria for all k = 1, 2, 3, 4 is a necessary and sufficient condition, any single k condition provides a necessary condition.
Since the eigenvalues of X (0,0) are generically the largest and therefore the most constraining, we will focus on bounds coming from this matrix. We begin with the upper left 2 × 2 submatrices. Taking the determinant, we have: Clearly the latter constraint coming from D L 2 (8π1 − X) will be the stronger of the two, since λ 1 , λ 2 > 0 if boundedness from below is imposed. Thus the necessary k = 2 condition reduces to Eq. (15b). We also examine the constraints that arise from using the lower-right and center 2 × 2 sub-matrices, as proposed in Section 3.3. The analytic conditions from the lower-right D R and center D C sub-matrices are, respectively, The combination of these three expressions provides a stronger constraint than just the upper-left minor constraint alone.
While it is immediately clear that D L 2 (8π1 − X) is stronger than the addition-based bound for the upper-left matrix, it is not clear for the center and lower-right matrices; in fact, including these bounds provides a slightly more constraining result. We thus employ the D C,R 2 (8π1 + X) constraints in our analysis of the k = 2 bound as well: (17b) Next, we look at the upper left 3 × 3 submatrices. Unlike the k = 2 case, it is not clear that one of these is generically more constraining than the other. To be consistent with the k = 2 case, we will examine the 8π1 − X (0,0) matrix. We additionally consider the lower-right 3 × 3 sub-matrix. These give the following bounds: While the D 3 (8π1 − X (0,0) ) provide the strongest constraints for values of |λ i | 4π, the inclusion of the D 3 (8π1 + X (0,0) ) and D 2 bounds improves the performance of the k = 3 bounds at higher |λ i |. We omit analytic expressions for the k = 4 case, since they cannot be simplified to a useful form. Moreover, the k = 3 expressions already provide constraints very close to the full numerical bound (see Fig. 2).

Sufficient conditions for perturbative unitarity
Next, we turn to derive sufficient conditions for perturbative unitarity by applying the Gershgorin disk theorem, which is reviewed in Section 3.2 and gives an upper bound on the maximal value of the eigenvalues. By demanding that this upper bound is less than 8π, we obtain a sufficient condition for perturbative unitarity.
We first construct the intervals x (Y,I) i containing the eigenvalues of each of the scattering matrices, X (0,0) , X (0,1) , X (1,0) , and X (1,1) . We know that in order to uphold perturbative unitarity, we must have |e i | < 8π. Thus we arrive at the sufficient condition: For each of the X matrices, we can work out the x explicitly. For X (0,0) , we obtain:  Figure 2.: Plot comparing the number of points that pass the exact numerical bound |e i | < 8π (black), the sufficient bound from the Gershgorin disk theorem Eq. (25) (dark blue), the sufficient bound from the Frobenius norm Eq. (26) (light blue), the necessary condition D 2 (8π1 ± X (0,0) ) > 0 (dark red), and the necessary condition D 2,3 (8π1 ± X (0,0) ) > 0 (light red). The λ i values are randomly chosen from the range of values satisfying |λ i | < λ max , where λ max is given by the x-axis in units of multiples of π. The minimal bounded from below condition λ 1,2 ≥ 0 is enforced. The λ 5,6,7 values are allowed to be complex. The total number of points checked for each λ max is 10,000.
For X (0,1) , they are: For X (1,1) , we have: Finally, X (1,0) , we have: In examining these conditions, the leading coefficient of 3 in the first set suggests that the x (0,0) i corresponding to X (0,0) will generically be larger than those corresponding to X (0,1) , X (1,0) , and X (1,1) ; a numerical check confirms this intuition. Thus the sufficient condition for perturbative unitarity simplifies slightly to One may alternatively employ the bound arising from the Frobenius norm, as discussed in Section 3.1. Taking the dominant X (0,0) matrix, one finds the condition: The dependence here on the signs of the λ i is similar to the dependence seen in the Gershgorin disk conditions: the bound is not sensitive to the signs of any λ i except for the relative sign between λ 3 and λ 4 .

Numerical comparison
In order to compare the various bounds derived in this section, in Fig. 2 we plot the number of points which pass the exact, sufficient, and necessary conditions for different values of λ max . For each λ max , we consider 10,000 randomly-drawn values for the λ i within the range |λ i | < λ max . For the necessary conditions, the results are derived from the combination of all possible 2 × 2 (3 × 3) sub-matrices along the diagonal for the D 2(3) bound. We enforce the minimal bounded from below condition λ 1,2 > 0, which has been derived e.g. in Ref. [5]. We find that the necessary D 3 (8π1 − X (0,0) ) condition lies very close to the exact condition and is effective at ruling out parameter sets which fail perturbative unitarity, while for λ max π all tested points satisfy perturbative unitarity.

Boundedness from below
Next, we seek to determine the conditions on the parameters such that the potential of Eq. (1) is bounded from below (BFB). For this, it is necessary to ensure that the quartic part of the potential does not acquire negative values. If negative values were present, one could easily find indefinite negative values of the potential by rescaling all fields to infinity in the same direction as the one in which the negative value was found. We remark that analytic expressions have been formulated previously in Ref. [30], though for the case of explicit CP conservation. We will make no such assumption. There are also previous analyses of the BFB condition using the eigenvalues of a 4×4 matrix (Ref. [24]); however, these analyses do not lead to analytical expressions, and we will follow an alternative approach.
If any of these conditions is true for a given set of input parameters λ 1...7 and for all possible values of ρ ∈ [0, 1], η ∈ [0, 2π), then the potential is BFB.
(34) Note that upon setting λ 6 = λ 7 = 0, and after extremizing with respect to η, Eq. (34) becomes This is a monotonic function of ρ, and hence the strongest constraints are derived for either ρ = 1 or ρ = 0, namely These conditions reproduce the well-known conditions for BFB in the Z 2 -symmetric 2HDM [5]. Let us also stress that, for ρ = 0, Eq. (34) leads to Eq. (37) independently of the value of the other quartic couplings and hence this equation is a necessary condition for the potential stability even in the generic 2HDM case.

Necessary conditions for boundedness from below
The two options of Eq. (34) present a necessary and sufficient condition for BFB. In order to implement this bound, one should scan over all possible values of ρ and η, which can be computationally expensive for large parameter spaces. Thus we present here two simplified necessary (though not sufficient) conditions which can be used to quickly rule out invalid parameter sets and speed up scans. We can first derive generalized versions of the existing literature bounds [5] by setting x = 1 and taking ρ = 1 and η = nπ 4 , with n = {0, ..., 7}, in Eq. (27) and Eq. (30). Applying this procedure to Eq. (27) leads to the following conditions: while applying the same procedure to Eq. (30) leads to the conditions: Note that we have combined the η, η + π conditions in each set to obtain four conditions instead of eight. Alternatively, we can collapse the two conditions of Eq. (34) into a single necessary condition as follows. Consider the two different branches with χ 1,2 < 0. Under the transformation η → η + π, χ 1 ≤ 0 produces two conditions that must be satisfied simultaneously: (α − γ) 2 − 16(α + β + γ + 2) ≤ 0 and (α − γ) 2 − 16(−α + β − γ + 2) ≤ 0. We can add these together to obtain the simplified condition: (α − γ) 2 ≤ 16(β + 2). Similarly, demanding that χ 2 ≤ 0 for both η and η + π gives us the simplified condition (α − γ) 2 ≤ 16(β + 2). So, we see that demanding χ 1 ≤ 0 and χ 2 ≤ 0 are equivalent, and both translate to the constraint: In this way, the condition for the potential to be BFB can be reduced to the form: Note that both β ≥ −2 and (α − γ) 2 ≤ 16(β + 2) restrict β, but that the latter will always be a stronger condition since (α − γ) 2 ≥ 0. Then this necessary but not sufficient BFB condition simplifies further to This condition still depends on ρ and η. Without loss of generality, we set 5 ρ = 1. As for η, we need to find the value which extremizes the expression for each condition. Take for instance the latter condition of Eq. (42) and define Using the definitions of Eqs. (29) and (32), we can recast everything in terms of cos 2η and sin 2η such that f (η) only depends on these quantities. We can then easily determine the extremal value of η min which gives the minimal f min . After some algebra, the positivity condition f min ≥ 0 reads: where we have defined the rescaled couplings: Eqs. (38), (39), and (44) present simplified necessary conditions for BFB and are the main result of this section.

Numerical analysis
To compare the performance of our analytic conditions with the numerical BFB condition, we perform a scan over 10,000 randomly chosen points in the 7-dimensional parameter space of {λ 1 , ...λ 7 }. We take as allowed ranges λ 1,2 ∈ [0, 2π] and |λ 3,4,5,6,7 | ≤ π 2 , with λ 5,6,7 complex, as this choice yields about half of the points BFB. Fig. 3 shows the number of points which pass the numerical condition Eq. (34) as well as the number which pass the combination of our necessary conditions Eqs. (38), (39), and (44) and the number which pass the combination of our sufficient conditions Eqs. (46) and (47). While we display in the figure the result of combining all necessary conditions derived in Section 5.1, we note that Eq. (39) provides the strongest necessary condition, with the combination of all conditions improving the results by a few percent.
We see that our necessary conditions are very effective at eliminating points which are not BFB, with only approximately 51% of points passing this condition, as compared with the 45% of points which actually satisfy BFB. Meanwhile, our analytic sufficient conditions guarantee approximately 11% of points are BFB. Of these points, essentially all are obtained from Eq. (46), which was derived from the upper branch of Eq. (34). The second condition Eq. (47), derived from the lower branch, is too strong and admits almost no points, which also reflects the fact that most of the points sampled fall within the regime of the first branch.
An examination of the analytical expressions indicates that the quantity √ λ 1 λ 2 may play an important role in the determination of BFB. To examine whether the analytical form of our bounds indeed captures the primary underlying behavior of BFB with respect to the λ i , we plot a histogram in √ λ 1 λ 2 of the fraction of tested points which pass the numerical, necessary, and sufficient bounds of Eqs. (34), (38), (39), (44) and (46), respectively. As in Fig. 3, we choose parameters in the range λ 1,2 ∈ [0, 2π] and |λ 3,4,5,6,7 | ≤ π/2 with λ 5,6,7 allowed to be complex. The resulting figure is shown in Fig. 4. As can be seen from the figure, more points pass the BFB condition for higher √ λ 1 λ 2 , as indicated by the forms of the necessary and sufficient conditions. We find that both the necessary and sufficient bounds follow the same behavior as the exact numerical results, indicating that the analytic bounds do indeed capture the relevant behavior in √ λ 1 λ 2 .   Fig. 3, we take as priors λ 1,2 ∈ [0, 2π] and |λ 3,4,5,6,7 | ∈ [0, π 2 ], with λ 5,6,7 allowed to be complex.
Finally, we note that within the existing literature, some simplified analytic BFB constraints for the most general 2HDM (i.e. involving λ 6,7 = 0) do exist. For example, the authors of Ref. [5] find as a necessary condition: This expression, that agrees with Eq. (38) in the appropriate limit, is derived by assuming that the Higgs doublets are aligned in field space, and is limited to the case that all λ i are taken to be real. Restricting ourselves to this regime, we find that the literature expression excludes approximately 17% of points while ours excludes approximately 49%, making our condition the stronger of the two by a large margin.

Vacuum stability
We can also place constraints on the allowed 2HDM potential parameters by demanding the existence of a stable neutral vacuum. Strictly speaking, this not a necessary requirement: it is only necessary that the vacuum is meta-stable, with a lifetime longer than the age of the Universe. Here, we just derive the conditions for absolute stability, more precisely the absence of deeper minima at scales of the order of the TeV scale. The discriminant D introduced in Refs. [48,49] offers a prescription for distinguishing the nature of a solution obtained by extremizing the potential. We summarize the method here, beginning by writing the potential as: where M µ encodes the mass terms: r µ is a vector of field bilinears: and Λ µν encodes the quartic terms: The last term in Eq. (49) is a Lagrange multiplier we have introduced to enforce the condition r µ r µ = 0, which ensures we are in a charge-neutral minimum; we enforce this condition since charge-breaking and normal minima cannot coexist in the 2HDM (see Refs. [20,48] for more details). In the above equations, indices are raised and lowered using a Minkowski metric.
Provided the matrix Λ µν , which contains the coefficients of the quartic terms in the potential, is positive definite, corresponding to a potential which is BFB, it can be brought into a diagonal form by an SO(1, 3) transformation: with Λ 0 the "timelike" eigenvalue and Λ i "spacelike". Let us define the "signature matrix" S as S ≡ Λ µν − ζg µν . In diagonal form, it looks like: The discriminant is generically given by the determinant of the signature matrix: By using the diagonal form above, we can write this as: We finally come to the vacuum stability condition. Suppose we have already verified that our potential is BFB and calculated the discriminant, time-like eigenvalue Λ 0 , and Lagrange multiplier ζ.
We are in a global minimum if and only if : For our purposes, it is more useful to work with the "Euclideanized" version of Λ µν obtained by lowering one of the indices with the Minkowski metric, Λ E ≡ Λ µ ν . Explicitly: In terms of Λ E , the discriminant is: The other quantity necessary for formulating the discriminant is the Lagrange multiplier ζ. This may be obtained by looking at any component of the minimization condition: We parameterize the vacuum expectation values (vevs) of the doublets as: Then the expectation value of field bilinears r µ ≡ r µ is: The expression for ζ is particularly simple if we choose the "1" component. In particular if we take η = 0, then: Note that this has the interpretation of the charged Higgs mass over the vev v squared, as first demonstrated in Refs. [24,25]. For a discussion of why this is generically true, see Appendix B.
If D > 0, then the physical minimum is the global one, implying absolute stability. If, instead, D < 0, we need to compare the timelike eigenvalue Λ 0 with ζ: we are in a global minimum if ζ > Λ 0 ; otherwise, the minimum is metastable. Provided we have already verified that the potential is BFB, however, there is an even simpler way to assess the nature of the extremum.
As an aside, working at the level of eigenvalues the two options of Eq. (57) for an extremum to be the global minimum can actually be collapsed into one. Recall that when the potential is BFB, Λ µν is positive definite and Λ 0 > Λ 1,2,3 . Then from Eq. (56), D > 0 necessarily implies that we have the ordering Λ 0 > ζ > Λ 1,2,3 . Similarly D < 0 and ζ > Λ 0 necessarily implies the ordering ζ > Λ 0 > Λ 1,2,3 . So we see that the relative ordering of Λ 0 and ζ does not actually matter−all that matters for a potential which has been verified to be BFB is that ζ be larger than the spatial eigenvalues, ζ > Λ 1,2,3 .

Gershgorin bounds
As in Section 4.3, we can bound the eigenvalues of Λ E using the Gershgorin disk theorem in order to derive a sufficient condition for a given vacuum solution to be stable. We first construct the intervals containing the eigenvalues of Λ E and define the endpoint of each interval as Γ i ≡ a ii + R i , with R i = j =i |a ij |: We know that all eigenvalues must be less than the endpoint of the interval extending the furthest in the +x direction, Meanwhile, an extremum will be the global minimum if ζ > Λ 1,2,3 . Thus, it is sufficient to demand:

Frobenius bounds
One may also bound the eigenvalues using the Frobenius norm to obtain a single-equation condition. We require the maximum eigenvalue be less than ζ, which gives the constraint: Note that in this case, the Frobenius bound is insensitive to the signs of the λ i , while the Gershgorin condition is sensitive to the signs of λ 3 , λ 4 , and λ R 5 . We thus expect the Gershgorin bound to be the stronger of the two.

Principal minors
In the case of a non-symmetric matrix such as Λ E , Sylvester's criterion no longer holds and so cannot be applied in a straightforward manner. However, the following statement does hold: if the symmetric part of a matrix M is positive-definite, then the real parts of the eigenvalues of M are positive. This statement does not hold in the other direction, and therefore cannot be used to derive necessary conditions. However, we can apply Sylvester's criterion to the symmetric part of Λ E to obtain a sufficient condition.
The symmetric part of Λ E is given by  We plot the fraction of points that pass each condition as a function of M H ± . We take λ 1,2 ∈ [0, 2π] and |λ 3,4,5,6,7 | ≤ π 2 , with λ 5,6,7 allowed to be complex, and restrict to examining points which are BFB.
We require the matrix ζ1 − Λ S E to be positive-definite. Since the lower-right 3×3 matrix decouples from the "11" element, we can analyze them separately when considering positive-definiteness. We require the "11" element to be positive, and apply Sylvester's criterion to the lower-right 3×3 submatrix. This gives the following set of conditions: Taken together, the Eqs. (70) provide a sufficient condition for vacuum stability.

Numerical comparison
In Fig. 5 we plot the performance of the three sufficient conditions for vacuum stability, Eqs. (67), (68), (70), as a function of the charged Higgs mass M H ± . We compare these results with the fraction that pass the exact stability condition, Eq. (57). As in previous sections, we choose the λ i randomly with λ 1,2 ∈ [0, 2π] and |λ 3,4,5,6,7 | ≤ π 2 , with λ 5,6,7 allowed to be complex. The y-axis shows the fraction of tested points which pass the respective stability condition; we restrict to testing points which are BFB, to ensure the validity of the stability conditions implemented in Fig. 5. We find that the set of conditions arising from the application of Sylvester's criterion capture the most stable points, while all three bounds capture more stable points when the λ i are small compared to the ratio M 2 H ± /v 2 .

Vacuum stability in the Higgs basis
It is particularly interesting to study vacuum stability in the Higgs basis, in which only one of the doublets possesses a vev (see Appendix A for a review of the conversion to the Higgs basis as well as our conventions). One advantage of this basis is that the potential parameters are closely related to physical observables: for example, Z 1 controls the trilinear coupling of three SM-like Higgs bosons hhh, Z 6 controls the trilinear coupling of two SM-like and one non-SM-like CP-even Higgs bosons hhH, etc. (see e.g. Ref. [39] for an exhaustive list of couplings). Since none of the bounds obtained in this article have relied on the choice of basis, they can equally well be applied to Higgs basis parameters. Using the close relationship between the Higgs basis parameters and physical quantities, we here aim at obtaining approximate bounds on the physical observables of the model. In our notation, given in Eq. (85) of Appendix A (choosing η = 0), the scalar which obtains a vev is denoted by φ 0 1 . The mass matrix for the neutral scalars φ = (φ 0 1 , φ 0 2 , a 0 ) T reads: where M 2 H ± is the charged Higgs mass: We will restrict ourselves to the alignment limit, which is the limit in which φ 0 1 is aligned with the 125 GeV mass eigenstate. In this case, the 125 GeV Higgs couples to the electroweak gauge bosons and all fermions with SM strength, and the alignment limit is therefore phenomenologically well-motivated by precision Higgs results from the LHC [33,34].
Examining the above matrix, it appears that there are two ways in which one may obtain alignment. The first option, known as the decoupling limit, corresponds to tak- Under this limit, the heavy mass eigenstates h 2 and h 3 and the heavy charged Higgs H ± decouple from the light mass eigenstate, leaving h 1 aligned with φ 0 1 . More interesting from a phenomenological standpoint is the approximate alignment without decoupling limit, as it leaves the non-standard Higgs states potentially within collider reach. This corresponds to taking |Z 6 | 1, for which the mixing between φ 0 1 and the other neutral scalars vanishes, leading to the identification of φ 0 1 with the mass eigenstate h 1 . For the following discussion we will take |Z 6 | 1 and work in the alignment without decoupling limit.
We define h 1 ≡ h to be the SM-like Higgs boson, which has a mass given by To obtain a physical Higgs mass close to the experimental value of 125 GeV, it is required that we fix Z 1 ≈ 0.25. The remaining 2×2 mass matrix can be diagonalized to obtain the masses of the remaining scalars h 2 and h 3 : There are two possibilities for the CP properties of these states. So long as Z I 5 = 0, h 2 and h 3 have mixed CP properties. In the limit of Z I 5 = 0, meanwhile, the non-standard Higgs mass matrix becomes diagonal, and we obtain mass eigenstates H and A with definite CP character, The masses of the general mass eigenstates and the states of definite CP character can be related by In the following analysis we will make no assumptions about the CP character of the mass eigenstates, and will work with the generic physical masses M h 3 ,h 2 .
With the above definitions, we can rephrase our sufficient vacuum stability conditions into constraints on physical quantities. We start with the Gershgorin condition of Eq. (67). Expressing the Γ's in terms of physical masses, a sufficient condition for vacuum stability becomes: Next, we can recast the Frobenius sufficient condition, Eq. (68), in terms of the physical masses; doing so results in the following condition: Finally, Sylvester's criterion provides an additional set of sufficient conditions. A sample set of sufficient conditions for vacuum stability in the alignment limit based on Eq. (70) is:

CP Violation in the general 2HDM
The bounds we have derived in this work have implications for the allowed values of physical parameters in a given 2HDM. This can be seen in a straightforward way in the previous section, where the conditions for vacuum stability were recast into expressions that restrict the physical masses of the bosonic sector. One particularly interesting question to which our bounds can be applied is that of the amount of CP violation permitted in the alignment limit. This possibility has largely been neglected in the many previous studies which restrict themselves to the Z 2 -symmetric 2HDM. This is understandable since exact alignment implies CP conservation in the Z 2 -symmetric case. When working in the fully general 2HDM, however, it is possible to have CP violation whilst still keeping the SM-like Higgs boson fully aligned.
To justify this claim, recall that there are four complex parameters in the 2HDM: {M 2 12 , Z 5 , Z 6 , Z 7 }. One of these is fixed by the minimization condition M 2 12 = − 1 2 Z 6 v 2 , leaving just three independent parameters, which we take to be the couplings {Z 5 , Z 6 , Z 7 }. These complex parameters enter into the three basic CP violating invariants of the 2HDM scalar sector J 1 , J 2 , and J 3 , which can be thought of as analogous to the Jarlskog invariant J of the SM quark sector. They are worked out explicitly in [50,51]; the important fact is that they scale as It is then clear that the condition for the Higgs sector to be CP invariant is: There are two ways in which this can be satisfied [52]: either Z I 5 = Z R 6 = Z R 7 = 0 or Z I 5 = Z I 6 = Z I 7 = 0. Note that in the limit of exact alignment we have Z 6 = 0, so this reduces to demanding either Z I 5 = Z R 7 = 0 or Z I 5 = Z I 7 = 0. Meanwhile in the 2HDM with a (softly broken) Z 2 symmetry, the fact that λ 6 = λ 7 = 0 implies the following two relations between the parameters in the Higgs basis (see Eq. (88) in Appendix A) [53,54] It immediately follows that: These conditions imply that in the exact alignment limit (i.e. Z 6 = 0), it will necessarily be the case that Z I 5 = Z I 7 = 0. Thus, exact alignment directly leads to CP conservation in the Z 2 -symmetric or softly broken Z 2 -symmetric 2HDM. This need not be the case in the fully general 2HDM, where the above relations no longer hold. If we allow for a small misalignment (i.e. |Z 6 | 0), |Z I 5 |, which controls the mixing between the H and A bosons (see Eq. (71)), can still be large for large tan β.
The physical consequences of the difference in the CP properties in the alignment limit between the (softly broken) Z 2 -symmetric and the general 2HDMs are illustrated in Fig. 6, in which we show the results of several parameter scans. Motivated by Higgs precision [33,34] and electric dipole moment bounds (see e.g. Refs. [55,56]), we demand a small mixing of the h and the H states and an even smaller mixing of the h and the A states; we impose these constraints on the mixing in an approximate manner by demanding | sin θ 12 | < 0.1 and | sin θ 13 | < 0.01, where θ 12 and θ 13 are the respective mixing angles. 7 The colour code in the figure indicates the minimum value of |θ 23 /π−1/4| in each hexagonal patch, where θ 23 is the mixing angle between the H and A states. This variable is chosen such that the maximal mixing case corresponds to a value of zero. Correspondingly, a dark blue color indicates that a large CP-violating mixing between the H and A states can be realized; a bright yellow color instead signals that no large mixing can be realized.
In the upper left panel of Fig. 6 we present a parameter scan of the 2HDM with a softly broken Z 2 symmetry in the (M h 3 − M h 2 , tan β) parameter plane. We observe that for large tan β, a large mixing between the neutral BSM Higgs bosons can be realized if their mass difference is below ∼ 70 GeV. Larger mass differences originate from differences between the diagonal terms of the Higgs mass matrix (see Eq. (71)), suppressing possible mixing effects induced by the off-diagonal 1 2 Z I 5 v 2 term. For lower tan β, the condition | sin θ 13 | < 0.01 directly implies that Z I 5 is small (see Eq. (84)), resulting in substantial mixing only when the mass difference is close to zero.
In the upper right panel of Fig. 6, we again consider a softy broken Z 2 symmetry but additionally impose perturbative unitarity, boundedness-from-below, and vacuum stability constraints following the discussions in the previous Sections. Aside from lowering the maximal possible mass difference between h 3 and h 2 , the region in which large mixing between the BSM Higgs bosons can be realized is also reduced.
If we instead investigate the general 2HDM without a (softly broken) Z 2 symmetry (see lower panels of Fig. 6), large mixing between the H and A states can be realized throughout the shown parameter plane. Applying the bounds derived in the previous sections only excludes the region with M h 3 − M h 2 100 GeV.  Finally, we want to remark that CP violation can become manifest not only in the neutral mass matrix but also in the bosonic couplings. This occurs if eitherZ R 7 = 0 or Z I 7 = 0, since these couplings enter into couplings like g h 1 h 2 h 3 [52]. Exotic decays like h 3 → h 1 h 2 → 3h 1 would then be indicative of CP violation in the bosonic sector (see e.g. [41]).

Conclusions
Two Higgs doublet models (2HDMs) present a natural extension of the Standard Model description. In spite of the simplicity of this SM extension, many new parameters appear in this theory, and it is very important to understand the constraints on these parameters which will impact in a relevant way the 2HDM phenomenology. Most existing studies concentrate on the case in which a Z 2 symmetry is imposed on the 2HDM potential and Yukawa sector. While this symmetry is an easy way to avoid flavor-changing neutral currents, it also forbids certain terms in the Higgs potential which do not induce flavorchanging neutral currents at tree level. In fact, in many scenarios in which the 2HDM is the low-energy effective field theory of a more complete high-scale model, these couplings are predicted to be non-zero.
Based on this motivation, in this work we present a step towards a systematic exploration of the non-Z 2 -symmetric 2HDM. We studied three of the most important theoretical constraints on the scalar potential parameters: perturbative unitarity, boundedness from below, and vacuum stability. In all three cases, we concentrated on the most general renormalizable potential (not restricted by any discrete symmetry) extending previous works by deriving analytic necessary and sufficient conditions for these constraints. For convenience, our main results (i.e., those conditions which approximate the exact conditions the best) are summarized in Table 1.
The derivation of our constraints makes use of several relevant mathematical properties, of which many have not been exploited in the literature before. These properties are not only applicable to the 2HDM but are also useful for the exploration of other models with extended Higgs sectors.
As a first phenomenological application of our bounds, we studied how much CPviolating mixing between the BSM Higgs bosons can be realized in the general 2HDM in comparison to a 2HDM with a (softly broken) Z 2 symmetry. While we found that large CP-violating mixing can only be realized for large tan β in the 2HDM with a softly broken Z 2 symmetry, no such theoretical constraints exist for the general 2HDM.
We leave a comprehensive study of the phenomenological consequences of not imposing a Z 2 symmetry for future work.
with H ± related to φ ± as: Thus we immediately see that Since ζ * = dV * dc , we verify the claim 13 of Eq. (89). Thus, it will generically be the case that the Lagrange multiplier accompanying the constraint r µ r µ = 0 is proportional to the charged Higgs mass squared in the neutral vacuum.