refname refsname disposition  EFI 22-8 FERMILAB-PUB-22-740-T

# New tools for dissecting the general 2HDM

Henning<EMAIL_ADDRESS>Department of Physics and Enrico Fermi
Institute, University of Chicago, 5720 South Ellis Avenue, Chicago, IL 60637
USA  Marcela<EMAIL_ADDRESS>Department of Physics and Enrico Fermi
Institute, University of Chicago, 5720 South Ellis Avenue, Chicago, IL 60637
USA  Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, Illinois,
60510, USA  Kavli Institute for Cosmological Physics, 5640 South Ellis Avenue,
University of Chicago, Chicago, IL 60637  Nina M<EMAIL_ADDRESS>Department of Physics and Enrico Fermi Institute, University of Chicago, 5720
South Ellis Avenue, Chicago, IL 60637 USA
Aurora<EMAIL_ADDRESS>Department of Physics and Enrico Fermi
Institute, University of Chicago, 5720 South Ellis Avenue, Chicago, IL 60637
USA  Carlos E.M<EMAIL_ADDRESS>Department of Physics and
Enrico Fermi Institute, University of Chicago, 5720 South Ellis Avenue,
Chicago, IL 60637 USA  Kavli Institute for Cosmological Physics, 5640 South
Ellis Avenue, University of Chicago, Chicago, IL 60637  HEP Division, Argonne
National Laboratory, 9700 Cass Ave., Argonne, IL 60439, USA

###### Abstract

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 ${\cal CP}$ violation.

###### Contents

  1. 1 Introduction
  2. 2 The general 2HDM
  3. 3 Methods for bounding matrix eigenvalues
    1. 3.1 Frobenius norm
    2. 3.2 Gershgorin disk theorem
    3. 3.3 Principal minors
  4. 4 Perturbative unitarity
    1. 4.1 Numerical bound
    2. 4.2 A necessary condition for perturbative unitarity
    3. 4.3 Sufficient conditions for perturbative unitarity
    4. 4.4 Numerical comparison
  5. 5 Boundedness from below
    1. 5.1 Necessary conditions for boundedness from below
    2. 5.2 Sufficient conditions for boundedness from below
    3. 5.3 Numerical analysis
  6. 6 Vacuum stability
    1. 6.1 Sufficient conditions for stability
      1. 6.1.1 Gershgorin bounds
      2. 6.1.2 Frobenius bounds
      3. 6.1.3 Principal minors
    2. 6.2 Numerical comparison
    3. 6.3 Vacuum stability in the Higgs basis
  7. 7 ${\cal CP}$ Violation in the general 2HDM
  8. 8 Conclusions
  9. A Higgs basis conversion
  10. B Connection between $\mathbf{\zeta}$ and $M_{H^{\pm}}^{2}$

##  1 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
$\mathbb{Z}_{2}$ symmetry. The different possible charge assignments for this
$\mathbb{Z}_{2}$ symmetry then fix the Higgs–fermion coupling choices and
define different types of 2HDMs.

This $\mathbb{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 [17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 5, 29] focus either on the
$\mathbb{Z}_{2}$-symmetric case or only provide a numerical procedure to test
these constraints in the general 2HDM (see Song:2022gsz 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 ${\cal 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.

##  2 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_{11}^{2}$ and $m_{22}^{2}$) are real while the
third, $m_{12}^{2}$, 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, $\lambda_{1}$ and
$\lambda_{2}$, must be real and, due to vacuum stability, positive. There are
two couplings associated with Hermitian combinations of the Higgs fields,
$\lambda_{3}$ and $\lambda_{4}$, which must be real, though not necessarily
positive. The coupling $\lambda_{5}$ is associated with the square of the
gauge invariant bilinear of both Higgs fields, and it may therefore be
complex. The couplings $\lambda_{6}$ and $\lambda_{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 $\lambda_{5}$, they
may be complex. The most general scalar potential for a complex 2HDM is
therefore:

$\begin{split}V&=m_{11}^{2}\Phi_{1}^{\dagger}\Phi_{1}+m_{22}^{2}\Phi_{2}^{\dagger}\Phi_{2}-(m_{12}^{2}\Phi_{1}^{\dagger}\Phi_{2}+h.c.)\\\
&+\frac{\lambda_{1}}{2}(\Phi_{1}^{\dagger}\Phi_{1})^{2}+\frac{\lambda_{2}}{2}(\Phi_{2}^{\dagger}\Phi_{2})^{2}+\lambda_{3}(\Phi_{1}^{\dagger}\Phi_{1})(\Phi_{2}^{\dagger}\Phi_{2})+\lambda_{4}(\Phi_{1}^{\dagger}\Phi_{2})(\Phi_{2}^{\dagger}\Phi_{1})\\\
&+\left[\frac{\lambda_{5}}{2}(\Phi_{1}^{\dagger}\Phi_{2})^{2}+\lambda_{6}(\Phi_{1}^{\dagger}\Phi_{1})(\Phi_{1}^{\dagger}\Phi_{2})+\lambda_{7}(\Phi_{2}^{\dagger}\Phi_{2})(\Phi_{1}^{\dagger}\Phi_{2})+h.c.\right]\,,\end{split}$
(1)

with $\Phi_{1,2}=(\Phi_{1,2}^{+},\Phi_{1,2}^{0})^{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 $\mathbb{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 $\lambda_{6}$
and $\lambda_{7}$ would violate parity symmetry and hence should vanish. The
mass parameter $m_{12}^{2}$ 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 $\lambda_{6}$ and $\lambda_{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 $\lambda_{6}$ and $\lambda_{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
$\mathbb{Z}_{2}$-symmetric 2HDM with vanishing $\lambda_{6}$ and $\lambda_{7}$
to avoid Higgs induced flavor violation in the fermion sector.111Note that
non-vanishing $\lambda_{6,7}$ induces flavor violation via Higgs mixing. This
effect is, however, loop suppressed. 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 (${\cal 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
$\mathbb{Z}_{2}$-symmetric 2HDM, however, ${\cal CP}$ is necessarily
conserved, as will be explained in detail in Section 7. In the full 2HDM, on
the other hand, one can have ${\cal CP}$ violation whilst remaining in exact
alignment thanks to the presence of $\lambda_{6}$ and $\lambda_{7}$ terms.
This ${\cal CP}$ violation could manifest in the neutral scalar mass
eigenstates as well as bosonic couplings (see also Haber:2022gsn), providing
many potential experimental signatures. With this motivation in mind, we keep
$\lambda_{6}$ and $\lambda_{7}$ non-zero throughout this work.

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

###  3.1 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: $|e_{i}|\leq||A||$.

where a matrix norm is defined as:

Definition: Given two $m\times n$ matrices $A$ and $B$, the matrix norm
$||A||$ satisfies the following properties:

  * •

$||A||\geq 0$,

  * •

$||A||=0\Leftrightarrow A=0_{m,n}$,

  * •

$||\alpha A||=|\alpha|||A||$,

  * •

$||A+B||\leq||A||+||B||$.

The above theorem holds for any choice of matrix norm, and thus one may employ
the Frobenius norm [42], $||A||=\sqrt{\text{Tr}(A^{\dagger}A)}$, to find the
following result:

$|e_{i}|\leq\sqrt{\text{Tr}(A^{\dagger}A)}$ (2)

This bound on the eigenvalues will be used to derive sufficient bounds in the
following sections.

###  3.2 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.222One can also construct
the radius by summing the magnitudes of the $n-1$ column entries. 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\times n$ matrix with entries $a_{ij}$, and
let $R_{i}$ be the sum of the magnitudes of the non-diagonal entries of the
$i^{\text{th}}$ row, $R_{i}=\sum_{j\neq 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}^{\pm}$ will be given by the sum
and difference, respectively, of the center and the radius,

$x_{i}^{\pm}\equiv a_{ii}\pm R_{i}\,,\,\,\,\text{with}\,\,\,R_{i}=\sum_{j\neq
i}|a_{ij}|\,.$ (3)

We then identify which $x_{i}^{\pm}$ 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,

$\text{min}(x_{i}^{-})\leq e_{k}\leq\text{max}(x_{i}^{+})\,.$ (4)

This may be rephrased into an upper bound on $|e_{k}|$ as:

$\displaystyle|e_{k}|\leq\mathrm{max}_{i}\left(\sum_{j}|a_{ij}|\right),$ (5)

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
$\sum_{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||=\text{max}_{i}(\sum_{j}|a_{ij}|)$.

###  3.3 Principal minors

In order to derive necessary conditions, one may employ Sylvester’s criterion
in a clever way, as proposed in Bento:2022vsb. Sylvester’s criterion involves
the principal minors $D_{k}$ of a matrix, where $D_{k}$ is the determinant of
the upper-left $k\times k$ sub-matrix. The statement of Sylvester’s criterion
is the following:

Theorem: Let M be a Hermitian $n\times n$ matrix. M is positive definite if
and only if all of the principal minors $D_{k}(M)$ are positive.

We further need the following result about Hermitian matrices333This theorem
is in fact used in the proof of Sylvester’s criterion:

Theorem: Let M be a Hermitian matrix. Then M is positive definite if and only
if all of its eigenvalues are positive.

One can apply this to derive an upper bound on the eigenvalues of a
diagonalizable matrix in the following way:

Theorem: Let M be an $n\times n$ diagonalizable, Hermitian (symmetric) matrix
and let $c$ be a positive real number. The eigenvalues $e_{i}$ of M are
bounded as $|e_{i}|<c$ if and only if all principal minors
$D_{k}(c\mathbb{1}-M)$ and $D_{k}(c\mathbb{1}+M)$ are positive for all
$k=1...n$.

To see this, consider applying a unitary transformation which diagonalizes $M$
to the matrix $c\mathbb{1}\pm\mathrm{diag}(M)$. Then for symmetric or
Hermitian matrices, the statement that $c\mathbb{1}\pm 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 $c\mathbb{1}\pm 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 $c\mathbb{1}-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, lower-right, and
central $2\times 2$ sub-matrices of a $4\times 4$ matrix. We will do so in
later analyses to strengthen the lower-$k$ bounds.

This use of sub-determinants has been proposed in Bento:2022vsb 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 Bento:2022vsb 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.

##  4 Perturbative unitarity

Tree-level constraints for perturbative unitarity in the most general 2HDM
have already been investigated in the literature [22, 45, 29]. However, for a
non-zero $\lambda_{6}$ and $\lambda_{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 $\lambda_{6}$ and
$\lambda_{7}$.

###  4.1 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\pi$. 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\rightarrow CD$, where the fields
$A...D$ represent any combination of the physical444Technically since we are
working in the high energy limit, the equivalence theorem allows us to replace
$W_{L}^{\pm}$ and $Z_{L}$ by their corresponding Goldstone bosons. scalars
$(H_{1},H_{2},H_{3},H^{\pm},W_{L}^{\pm},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}^{\pm},h_{i},z_{i})$, appearing as

$\Phi_{i}=\begin{pmatrix}w_{i}^{+}\\\
\frac{1}{\sqrt{2}}(v_{i}+h_{i}+iz_{i})\end{pmatrix}\,,$ (6)

with $v=\sqrt{\sum_{i}v_{i}^{2}}$ = 246 GeV.
Out of these fields, we can construct 14 neutral two-body states:
$\ket{w_{i}^{+}w_{i}^{-}}$, $\frac{1}{\sqrt{2}}\ket{z_{i}z_{i}}$,
$\frac{1}{\sqrt{2}}\ket{h_{i}h_{i}}$, $\ket{h_{i}z_{i}}$,
$\ket{w_{1}^{+}w_{2}^{-}}$, $\ket{w_{2}^{+}w_{1}^{-}}$, $\ket{z_{1}z_{2}}$,
$\ket{h_{1}h_{2}}$, $\ket{z_{1}h_{2}}$, and $\ket{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
Ginzburg:2005dt,Kanemura:2015ska,Jurciukonis:2018skr). For the neutral
scattering channels, this is:

$a_{0}^{(0)}=\frac{1}{16\pi}\begin{pmatrix}X_{(0,0)}&&&\\\ &X_{(0,1)}&&\\\
&&X_{(1,1)}&\\\ &&&X_{(1,1)}\end{pmatrix}\,,$ (7)

where the subscript of each submatrix denotes the quantum numbers $(Y,I)$ of
the corresponding states. The entries are:

$\displaystyle X_{(0,0)}$
$\displaystyle=\begin{pmatrix}3\lambda_{1}&2\lambda_{3}+\lambda_{4}&3\lambda_{6}&3\lambda_{6}^{*}\\\
2\lambda_{3}+\lambda_{4}&3\lambda_{2}&3\lambda_{7}&3\lambda_{7}^{*}\\\
3\lambda_{6}^{*}&3\lambda_{7}^{*}&\lambda_{3}+2\lambda_{4}&3\lambda_{5}^{*}\\\
3\lambda_{6}&3\lambda_{7}&3\lambda_{5}&\lambda_{3}+2\lambda_{4}\end{pmatrix}\,,$
(8a) $\displaystyle X_{(0,1)}$
$\displaystyle=\begin{pmatrix}\lambda_{1}&\lambda_{4}&\lambda_{6}&\lambda_{6}^{*}\\\
\lambda_{4}&\lambda_{2}&\lambda_{7}&\lambda_{7}^{*}\\\
\lambda_{6}^{*}&\lambda_{7}^{*}&\lambda_{3}&\lambda_{5}^{*}\\\
\lambda_{6}&\lambda_{7}&\lambda_{5}&\lambda_{3}\end{pmatrix}\,,$ (8b)
$\displaystyle X_{(1,1)}$
$\displaystyle=\begin{pmatrix}\lambda_{1}&\lambda_{5}&\sqrt{2}\lambda_{6}\\\
\lambda_{5}^{*}&\lambda_{2}&\sqrt{2}\lambda_{7}^{*}\\\
\sqrt{2}\lambda_{6}^{*}&\sqrt{2}\lambda_{7}&\lambda_{3}+\lambda_{4}\end{pmatrix}\,.$
(8c)

For the 8 singly-charged two-body states $\ket{w_{i}^{+}z_{i}}$,
$\ket{w_{i}^{+}h_{i}}$, $\ket{w_{1}^{+}z_{2}}$, $\ket{w_{1}^{+}h_{2}}$,
$\ket{w_{2}^{+}z_{1}}$, $\ket{w_{2}^{+}h_{1}}$, the block diagonal 8$\times$8
singly-charged S-matrix is given by:

$a_{0}^{(+)}=\frac{1}{16\pi}\begin{pmatrix}X_{(0,1)}&&\\\ &X_{(1,0)}&\\\
&&X_{(1,1)}\end{pmatrix}\,,$ (9)

where the new entry $X_{(1,0)}$ is just the one-dimensional eigenvalue:

$X_{(1,0)}=\lambda_{3}-\lambda_{4}\,.$ (10)

Finally, the 3$\times$3 S-matrix for the three doubly-charged 2-body states
$\ket{w_{i}^{+}w_{i}^{+}}$, $\ket{w_{1}^{+}w_{2}^{+}}$ is given by:

$a_{0}^{(++)}=\frac{1}{16\pi}X_{(1,1)}\,.$ (11)

We impose perturbative unitarity by demanding that the eigenvalues of the
scattering matrix are smaller than $8\pi$ implying that $|a_{0}|<\frac{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

$|e_{i}|<8\pi\,.$ (12)

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 $\lambda_{1...7}$ to be equal, the strongest constraint arises
from the 4$\times$4 matrix $X_{00}$. If we set all $\lambda_{i}\equiv\lambda$
and solve for the eigenvalues, we find the bound:

$\lambda<\frac{2\pi}{3}\,.$ (13)

This value is an order of magnitude smaller than $8\pi$, implying that if all
quartic couplings are sizable (i.e. of ${\cal O}(1)$), perturbative unitarity
may be lost even at values of the couplings much smaller than $4\pi$, which is
a bound often encountered in the literature to ensure perturbativity.

Figure 1: Plot showing the fraction of points that pass the unitarity bound
$|e_{i}|<8\pi$ for different choices of $\lambda_{max}$, in units of multiples
of $\pi$. The values of the $\lambda_{i}$ are each chosen randomly such that
$|\lambda_{i}|<\lambda_{max}$. We test 20,000 random sets of $\lambda_{i}$ for
each $\lambda_{max}$.

We investigate the validity of such an upper bound further in Fig. 1. For this
figure, we randomly choose each $\lambda_{i}$ within the range
$|\lambda_{i}|<\lambda_{\text{max}}$ and then show the fraction of test points
which pass the numerical unitarity constraint, as a function of
$\lambda_{\text{max}}$. For $\lambda_{\text{max}}\lesssim\pi$ almost all
points survive the perturbative unitarity constraint. For larger
$\lambda_{\text{max}}$ values the survival rate quickly drops to almost zero
for $\lambda_{\text{max}}\gtrsim 4\pi$. This highlights again that the
simplified perturbativity bound of $|\lambda_{i}|<4\pi$, which is often
encountered in the literature, is too loose. Based on the results in Fig. 1, a
better choice of bound might be $|\lambda_{i}|\lesssim\pi$ or
$|\lambda_{i}|\lesssim 3\pi/2$.

###  4.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\pi$, as required by perturbative unitarity,
amounts to requiring:

$D_{k}(8\pi\mathbb{1}+X)>0\,\,\text{and}\,\,D_{k}(8\pi\mathbb{1}-X)>0\,,$ (14)

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\times 2$ submatrices. Taking the determinant, we
have:

$D_{2}^{L}(8\pi\mathbb{1}+X_{(0,0)})>0\,\,\,\Rightarrow\,\,\,64\pi^{2}+24\pi(\lambda_{1}+\lambda_{2})+9\lambda_{1}\lambda_{2}-(2\lambda_{3}+\lambda_{4})^{2}>0\,,$
(15a)
$D_{2}^{L}(8\pi\mathbb{1}-X_{(0,0)})>0\,\,\,\Rightarrow\,\,\,64\pi^{2}-24\pi(\lambda_{1}+\lambda_{2})+9\lambda_{1}\lambda_{2}-(2\lambda_{3}+\lambda_{4})^{2}>0\,.$
(15b)

Clearly the latter constraint coming from $D_{2}^{L}(8\pi\mathbb{1}-X)$ will
be the stronger of the two, since $\lambda_{1},\lambda_{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\times 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,

$\displaystyle D_{2}^{R}(8\pi\mathbb{1}-X_{(0,0)})$
$\displaystyle=64\pi^{2}+(\lambda_{3}-16\pi)\lambda_{3}+(4(\lambda_{3}+\lambda_{4})-32\pi)\lambda_{4}-9|\lambda_{5}|^{2}$
$\displaystyle>0\,,$ (16a) $\displaystyle D_{2}^{C}(8\pi\mathbb{1}-X_{(0,0)})$
$\displaystyle=64\pi^{2}-8\pi(\lambda_{3}+3\lambda_{2}+2\lambda_{4})+3\lambda_{3}\lambda_{2}+6\lambda_{2}\lambda_{4}-9|\lambda_{7}|^{2}$
$\displaystyle>0\,.$ (16b)

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_{2}^{L}(8\pi\mathbb{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_{2}^{C,R}(8\pi\mathbb{1}+X)$ constraints in our analysis of the $k=2$ bound
as well:

$\displaystyle D_{2}^{R}(8\pi\mathbb{1}+X_{(0,0)})$
$\displaystyle=64\pi^{2}+(\lambda_{3}+16\pi)\lambda_{3}+(4(\lambda_{3}+\lambda_{4})+32\pi)\lambda_{4}-9|\lambda_{5}|^{2}$
$\displaystyle>0\,,$ (17a) $\displaystyle D_{2}^{C}(8\pi\mathbb{1}+X_{(0,0)})$
$\displaystyle=64\pi^{2}+8\pi(\lambda_{3}+3\lambda_{2}+2\lambda_{4})+3\lambda_{3}\lambda_{2}+6\lambda_{2}\lambda_{4}-9|\lambda_{7}|^{2}$
$\displaystyle>0\,.$ (17b)

Next, we look at the upper left $3\times 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\pi\mathbb{1}-X_{(0,0)}$ matrix. We additionally consider the lower-right
$3\times 3$ sub-matrix. These give the following bounds:

$\begin{split}D_{3}^{L}(8\pi\mathbb{1}-X_{(0,0)})=&(8\pi-\lambda_{3}-2\lambda_{4})((8\pi-3\lambda_{1})(8\pi-3\lambda_{2})-(2\lambda_{3}+\lambda_{4})^{2})\\\
&-9(8\pi-3\lambda_{2})|\lambda_{6}|^{2}-9(8\pi-3\lambda_{1})|\lambda_{7}|^{2}\\\
&-9(2\lambda_{3}+\lambda_{4})(\lambda_{6}\lambda_{7}^{*}+\lambda_{6}^{*}\lambda_{7})>0\,,\end{split}$
(18a)
$\begin{split}D_{3}^{R}(8\pi\mathbb{1}-X_{(0,0)})=&(8\pi-\lambda_{3}-2\lambda_{4})^{2}(8\pi-3\lambda_{2})-9|\lambda_{5}|^{2}(8\pi-3\lambda_{2})\\\
&-27(\lambda_{5}^{*}\lambda_{7}^{2}+\lambda_{5}(\lambda_{7}^{*})^{2})-18(8\pi-\lambda_{3}-2\lambda_{4})|\lambda_{7}|^{2}>0\,.\end{split}$
(18b)

Meanwhile, considering $8\pi\mathbb{1}+X_{(0,0)}$ gives

$\begin{split}D_{3}^{L}(8\pi\mathbb{1}+X_{(0,0)})=&(8\pi+\lambda_{3}+2\lambda_{4})((8\pi+3\lambda_{1})(8\pi+3\lambda_{2})-(2\lambda_{3}+\lambda_{4})^{2})\\\
&-9(8\pi+3\lambda_{2})|\lambda_{6}|^{2}-9(8\pi+3\lambda_{1})|\lambda_{7}|^{2}\\\
&+9(2\lambda_{3}+\lambda_{4})(\lambda_{6}\lambda_{7}^{*}+\lambda_{6}^{*}\lambda_{7})>0\,,\end{split}$
(19a)
$\begin{split}D_{3}^{R}(8\pi\mathbb{1}+X_{(0,0)})=&(8\pi+\lambda_{3}+2\lambda_{4})^{2}(8\pi+3\lambda_{2})-9|\lambda_{5}|^{2}(8\pi+3\lambda_{2})\\\
&+27(\lambda_{5}^{*}\lambda_{7}^{2}+\lambda_{5}(\lambda_{7}^{*})^{2})-18(8\pi+\lambda_{3}+2\lambda_{4})|\lambda_{7}|^{2}>0\,.\end{split}$
(19b)

While the $D_{3}(8\pi\mathbb{1}-X_{(0,0)})$ provide the strongest constraints
for values of $|\lambda_{i}|\lesssim 4\pi$, the inclusion of the
$D_{3}(8\pi\mathbb{1}+X_{(0,0)})$ and $D_{2}$ bounds improves the performance
of the $k=3$ bounds at higher $|\lambda_{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).

###  4.3 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\pi$, we obtain a sufficient condition
for perturbative unitarity.

Figure 2: Plot comparing the number of points that pass the exact numerical
bound $|e_{i}|<8\pi$ (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\pi\mathbb{1}\pm
X_{(0,0)})>0$ (dark red), and the necessary condition
$D_{2,3}(8\pi\mathbb{1}\pm X_{(0,0)})>0$ (light red). The $\lambda_{i}$ values
are randomly chosen from the range of values satisfying
$|\lambda_{i}|<\lambda_{max}$, where $\lambda_{max}$ is given by the x-axis in
units of multiples of $\pi$. The minimal bounded from below condition
$\lambda_{1,2}\geq 0$ is enforced. The $\lambda_{5,6,7}$ values are allowed to
be complex. The total number of points checked for each $\lambda_{max}$ is
10,000.

We first construct the intervals $x_{i}^{(Y,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\pi$. Thus we arrive at the sufficient condition:

$\text{max}(x_{i}^{(Y,I)})<8\pi\,.$ (20)

For each of the $X$ matrices, we can work out the $x_{i}^{(Y,I)}$ explicitly.
For $X_{(0,0)}$, we obtain:

$\displaystyle x_{1}^{(0,0)}$
$\displaystyle=3|\lambda_{1}|+(|2\lambda_{3}+\lambda_{4}|+6|\lambda_{6}|)\,,$
(21a) $\displaystyle x_{2}^{(0,0)}$
$\displaystyle=3|\lambda_{2}|+(|2\lambda_{3}+\lambda_{4}|+6|\lambda_{7}|)\,,$
(21b) $\displaystyle x_{3}^{(0,0)}$
$\displaystyle=x_{4}^{(0,0)}=|\lambda_{3}+2\lambda_{4}|+3(|\lambda_{5}|+|\lambda_{6}|+|\lambda_{7}|)\,.$
(21c)

For $X_{(0,1)}$, they are:

$\displaystyle x_{1}^{(0,1)}$
$\displaystyle=|\lambda_{1}|+(|\lambda_{4}|+2|\lambda_{6}|)\,,$ (22a)
$\displaystyle x_{2}^{(0,1)}$
$\displaystyle=|\lambda_{2}|+(|\lambda_{4}|+2|\lambda_{7}|)\,,$ (22b)
$\displaystyle x_{3}^{(0,1)}$
$\displaystyle=x_{4}^{(0,1)}=|\lambda_{3}|+(|\lambda_{5}|+|\lambda_{6}|+|\lambda_{7}|)\,.$
(22c)

For $X_{(1,1)}$, we have:

$\displaystyle x_{1}^{(1,1)}$
$\displaystyle=|\lambda_{1}|+(|\lambda_{5}|+\sqrt{2}|\lambda_{6}|)\,,$ (23a)
$\displaystyle x_{2}^{(1,1)}$
$\displaystyle=|\lambda_{2}|+(|\lambda_{5}|+\sqrt{2}|\lambda_{7}|)\,,$ (23b)
$\displaystyle x_{3}^{(1,1)}$
$\displaystyle=|\lambda_{3}+\lambda_{4}|+\sqrt{2}(|\lambda_{6}|+|\lambda_{7}|)\,.$
(23c)

Finally, $X_{(1,0)}$, we have:

$\displaystyle x_{1}^{(1,0)}$ $\displaystyle=|\lambda_{3}-\lambda_{4}|\,.$
(24)

In examining these conditions, the leading coefficient of 3 in the first set
suggests that the $x_{i}^{(0,0)}$ 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

$\text{max}\big{(}x_{i}^{(0,0)}\big{)}<8\pi\,.$ (25)

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:

$\displaystyle\sqrt{9(\lambda_{1}^{2}+\lambda_{2}^{2})+10(\lambda_{3}^{2}+\lambda_{4}^{2})+16\lambda_{3}\lambda_{4}+18(|\lambda_{5}|^{2}+2|\lambda_{6}|^{2}+2|\lambda_{7}|^{2})}\leq
8\pi$ (26)

The dependence here on the signs of the $\lambda_{i}$ is similar to the
dependence seen in the Gershgorin disk conditions: the bound is not sensitive
to the signs of any $\lambda_{i}$ except for the relative sign between
$\lambda_{3}$ and $\lambda_{4}$.

###  4.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 $\lambda_{max}$. For each $\lambda_{max}$,
we consider 10,000 randomly-drawn values for the $\lambda_{i}$ within the
range $|\lambda_{i}|<\lambda_{max}$. For the necessary conditions, the results
are derived from the combination of all possible $2\times 2$ ($3\times 3$)
sub-matrices along the diagonal for the $D_{2(3)}$ bound. We enforce the
minimal bounded from below condition $\lambda_{1,2}>0$, which has been derived
e.g. in Branco:2011iw. We find that the necessary
$D_{3}(8\pi\mathbb{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 $\lambda_{max}\lesssim\pi$ all tested points
satisfy perturbative unitarity.

##  5 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 Song:2022gsz,
though for the case of explicit ${\cal CP}$ conservation. We will make no such
assumption. There are also previous analyses of the BFB condition using the
eigenvalues of a 4$\times$4 matrix (Ivanov:2006yq); however, these analyses do
not lead to analytical expressions, and we will follow an alternative
approach.

We begin by reparameterizing the potential via
$\Phi_{1}^{\dagger}\Phi_{1}=\frac{1}{2}h_{1}^{2}$,
$\Phi_{2}^{\dagger}\Phi_{2}=\frac{1}{2}h_{2}^{2}$,
$\Phi_{1}^{\dagger}\Phi_{2}=\frac{1}{2}h_{1}h_{2}\rho e^{i\eta}$, with
$\rho\in[0,1]$. Moreover, we decompose the complex couplings $\lambda_{5}$,
$\lambda_{6}$, $\lambda_{7}$ into real and imaginary parts as:
$\lambda_{i}e^{i\eta}+\lambda_{i}^{*}e^{-i\eta}=2\text{Re}[\lambda_{i}]\cos\eta-2\text{Im}[\lambda_{i}]\sin\eta$.
The quartic part of the potential then becomes:

$\begin{split}V_{\text{quartic}}={}\frac{1}{4}\bigg{\\{}&\frac{\lambda_{1}}{2}h_{1}^{4}+\frac{\lambda_{2}}{2}h_{2}^{4}+\big{[}\lambda_{3}+\big{(}\lambda_{4}+\text{Re}[\lambda_{5}]\cos
2\eta-\text{Im}[\lambda_{5}]\sin
2\eta\big{)}\rho^{2}\big{]}h_{1}^{2}h_{2}^{2}\\\
&+2\big{(}\text{Re}[\lambda_{6}]\cos\eta-\text{Im}[\lambda_{6}]\sin\eta\big{)}\rho\,h_{1}^{3}h_{2}\\\
&+2\big{(}\text{Re}[\lambda_{7}]\cos\eta-\text{Im}[\lambda_{7}]\sin\eta\big{)}\rho\,h_{1}h_{2}^{3}\bigg{\\}}\,.\end{split}$
(27)

We can then cast $V_{\text{quartic}}$ into the form

$V_{\text{quartic}}=\frac{1}{4}h_{2}^{4}\left[a\left(\frac{h_{1}}{h_{2}}\right)^{4}+b\left(\frac{h_{1}}{h_{2}}\right)^{3}+c\left(\frac{h_{1}}{h_{2}}\right)^{2}+d\left(\frac{h_{1}}{h_{2}}\right)+e\right]\,,$
(28)

with

$\displaystyle a$ $\displaystyle=\frac{\lambda_{1}}{2},\hskip 8.5359pt\
e=\frac{\lambda_{2}}{2}\,,$ (29a) $\displaystyle b$
$\displaystyle=2\big{(}\text{Re}[\lambda_{6}]\cos\eta-\text{Im}[\lambda_{6}]\sin\eta\big{)}\rho\,,$
(29b) $\displaystyle c$
$\displaystyle=\left[\lambda_{3}+\big{(}\lambda_{4}+\text{Re}[\lambda_{5}]\cos
2\eta-\text{Im}[\lambda_{5}]\sin 2\eta\big{)}\rho^{2}\right]\,,$ (29c)
$\displaystyle d$
$\displaystyle=2\big{(}\text{Re}[\lambda_{7}]\cos\eta-\text{Im}[\lambda_{7}]\sin\eta\big{)}\rho\,.$
(29d)

Clearly $a>0$ and $e>0$ has to be fulfilled as a minimum condition for BFB. We
can then divide out by $e\,h_{2}^{4}$ and define the simplified polynomial:

$f(x)=x^{4}+\alpha x^{3}+\beta x^{2}+\gamma x+1\,,$ (30)

with $x=\frac{a^{1/4}}{e^{1/4}}\frac{h_{1}}{h_{2}}$ and:

$\alpha=ba^{-3/4}e^{-1/4},\hskip 14.22636pt\beta=ca^{-1/2}e^{-1/2},\hskip
14.22636pt\gamma=da^{-1/4}e^{-3/4}\,.$ (31)

We then define the following quantities:

$\displaystyle\Delta$
$\displaystyle=4[\beta^{2}-3\alpha\gamma+12]^{3}-[72\beta+9\alpha\beta\gamma-2\beta^{3}-27\alpha^{2}-27\gamma^{2}]^{2}\,,$
(32a) $\displaystyle\chi_{1}$
$\displaystyle=(\alpha-\gamma)^{2}-16(\alpha+\beta+\gamma+2)\,,\,\,\,$ (32b)
$\displaystyle\chi_{2}$
$\displaystyle=(\alpha-\gamma)^{2}-\frac{4(\beta+2)}{\sqrt{\beta-2}}\left(\alpha+\gamma+4\sqrt{\beta-2}\right)\,.$
(32c)

The positivity of $V_{\text{quartic}}$ is ensured if and only if one of the
following conditions holds [47]:

$\begin{split}&(1)\,\,\,\beta<-2\,\,\text{and}\,\,\Delta\leq
0\,\,\text{and}\,\,\alpha+\gamma>0\,,\\\ &(2)\,\,\,-2\leq\beta\leq
6\,\,\text{and}\,\,\Delta\leq 0\,\,\text{and}\,\,\alpha+\gamma>0\,,\\\
&(3)\,\,\,-2\leq\beta\leq 6\,\,\text{and}\,\,\Delta\geq
0\,\,\text{and}\,\,\chi_{1}\leq 0\,,\\\
&(4)\,\,\,\beta>6\,\,\text{and}\,\,\Delta\leq
0\,\,\text{and}\,\,\alpha+\gamma>0\,,\\\
&(5)\,\,\,\beta>6\,\,\text{and}\,\,\alpha>0\,\,\text{and}\,\,\gamma>0\,,\\\
&(6)\,\,\,\beta>6\,\,\text{and}\,\,\Delta\geq 0\,\,\text{and}\,\,\chi_{2}\leq
0\,.\end{split}$ (33)

If any of these conditions is true for a given set of input parameters
$\lambda_{1...7}$ and for all possible values of $\rho\in[0,1]$,
$\eta\in[0,2\pi)$, then the potential is BFB.

Note that under the transformation $\eta\rightarrow\eta+\pi$, both $\alpha$
and $\gamma$ are anti-symmetric ($\alpha\rightarrow-\alpha$ and
$\gamma\rightarrow-\gamma$). This in turn implies that conditions (1), (2),
(4), and (5) are always violated for some value of $\eta$, and therefore can
never guarantee the positivity of $V_{\text{quartic}}$. Consequently, we are
left with only two conditions under which the potential is BFB:

$\text{The potential is BFB if and only if}:\,\,\Delta\geq
0\,\,\text{and}\,\,\begin{cases}-2\leq\beta\leq
6\,\,\text{and}\,\,\chi_{1}\leq 0\,,\text{or}\\\
\beta>6\,\,\text{and}\,\,\chi_{2}\leq 0\,.\end{cases}$ (34)

Note that upon setting $\lambda_{6}=\lambda_{7}=0$, and after extremizing with
respect to $\eta$, Eq. 34 becomes

$\begin{split}\beta+2\geq
0\;\Rightarrow\;\lambda_{3}+\rho^{2}(\lambda_{4}-|\lambda_{5}|)\geq-\sqrt{\lambda_{1}\lambda_{2}}\,.\end{split}$
(35)

This is a monotonic function of $\rho$, and hence the strongest constraints
are derived for either $\rho=1$ or $\rho=0$, namely

$\displaystyle\lambda_{3}+\lambda_{4}-|\lambda_{5}|$
$\displaystyle\geq-\sqrt{\lambda_{1}\lambda_{2}}\,,\;\;\textrm{and}$ (36)
$\displaystyle\lambda_{3}$
$\displaystyle\geq-\sqrt{\lambda_{1}\lambda_{2}}\,,$ (37)

These conditions reproduce the well-known conditions for BFB in the
$\mathbb{Z}_{2}$-symmetric 2HDM [5]. Let us also stress that, for $\rho=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.

###  5.1 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 $\rho$ and $\eta$, 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 $\rho=1$ and $\eta=\frac{n\pi}{4}$, with
$n=\\{0,...,7\\}$, in Eq. 27 and Eq. 30. Applying this procedure to Eq. 27
leads to the following conditions:

$\displaystyle\frac{\lambda_{1}+\lambda_{2}}{2}+\lambda_{3}+\lambda_{4}+\lambda_{5}^{R}-2|\lambda_{6}^{R}+\lambda_{7}^{R}|>0\,,$
(38a)
$\displaystyle\frac{\lambda_{1}+\lambda_{2}}{2}+\lambda_{3}+\lambda_{4}-\lambda_{5}^{R}-2|\lambda_{6}^{I}+\lambda_{7}^{I}|>0\,,$
(38b)
$\displaystyle\frac{\lambda_{1}+\lambda_{2}}{2}+\lambda_{3}+\lambda_{4}+\lambda_{5}^{I}-\sqrt{2}\;\Big{|}(\lambda_{6}^{R}+\lambda_{7}^{R})+(\lambda_{6}^{I}+\lambda_{7}^{I})\Big{|}>0\,,$
(38c)
$\displaystyle\frac{\lambda_{1}+\lambda_{2}}{2}+\lambda_{3}+\lambda_{4}-\lambda_{5}^{I}-\sqrt{2}\;\Big{|}(\lambda_{6}^{R}+\lambda_{7}^{R})-(\lambda_{6}^{I}+\lambda_{7}^{I})\Big{|}>0\,,$
(38d)

while applying the same procedure to Eq. 30 leads to the conditions:

$\displaystyle\sqrt{\lambda_{1}\lambda_{2}}+\lambda_{3}+\lambda_{4}+\lambda_{5}^{R}-2\left|\tilde{\lambda}_{6}^{R}+\tilde{\lambda}_{7}^{R}\right|>0\,,$
(39a)
$\displaystyle\sqrt{\lambda_{1}\lambda_{2}}+\lambda_{3}+\lambda_{4}-\lambda_{5}^{R}-2\left|\tilde{\lambda}_{6}^{I}+\tilde{\lambda}_{7}^{I}\right|>0\,,$
(39b)
$\displaystyle\sqrt{\lambda_{1}\lambda_{2}}+\lambda_{3}+\lambda_{4}+\lambda_{5}^{I}-\sqrt{2}\left|(\tilde{\lambda}_{6}^{R}+\tilde{\lambda}_{7}^{R})+(\tilde{\lambda}_{6}^{I}+\tilde{\lambda}_{7}^{I})\right|>0\,,$
(39c)
$\displaystyle\sqrt{\lambda_{1}\lambda_{2}}+\lambda_{3}+\lambda_{4}-\lambda_{5}^{I}-\sqrt{2}\left|(\tilde{\lambda}_{6}^{R}+\tilde{\lambda}_{7}^{R})-(\tilde{\lambda}_{6}^{I}+\tilde{\lambda}_{7}^{I})\right|>0\,.$
(39d)

Note that we have combined the $\eta,\eta+\pi$ 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
$\chi_{1,2}<0$. Under the transformation $\eta\rightarrow\eta+\pi$,
$\chi_{1}\leq 0$ produces two conditions that must be satisfied
simultaneously: $(\alpha-\gamma)^{2}-16(\alpha+\beta+\gamma+2)\leq 0$ and
$(\alpha-\gamma)^{2}-16(-\alpha+\beta-\gamma+2)\leq 0$. We can add these
together to obtain the simplified condition: $(\alpha-\gamma)^{2}\leq
16(\beta+2)$. Similarly, demanding that $\chi_{2}\leq 0$ for both $\eta$ and
$\eta+\pi$ gives us the simplified condition $(\alpha-\gamma)^{2}\leq
16(\beta+2)$. So, we see that demanding $\chi_{1}\leq 0$ and $\chi_{2}\leq 0$
are equivalent, and both translate to the constraint:

$\chi_{1}\leq 0\,,\chi_{2}\leq 0\,\,\Rightarrow\,\,(\alpha-\gamma)^{2}\leq
16(\beta+2)\,.$ (40)

In this way, the condition for the potential to be BFB can be reduced to the
form:

$\Delta\geq
0\,\,\text{and}\,\,\beta\geq-2\,\,\text{and}\,\,(\alpha-\gamma)^{2}\leq
16(\beta+2)\,.$ (41)

Note that both $\beta\geq-2$ and $(\alpha-\gamma)^{2}\leq 16(\beta+2)$
restrict $\beta$, but that the latter will always be a stronger condition
since $(\alpha-\gamma)^{2}\geq 0$. Then this necessary but not sufficient BFB
condition simplifies further to

$\Delta\geq 0\,\,\text{and}\,\,\beta\geq\frac{1}{16}(\alpha-\gamma)^{2}-2\,.$
(42)

This condition still depends on $\rho$ and $\eta$. Without loss of generality,
we set555This gives us one necessary condition. We could obtain others by
choosing $\rho<1$, but these tend to be less constraining in most cases.
$\rho=1$. As for $\eta$, we need to find the value which extremizes the
expression for each condition. Take for instance the latter condition of Eq.
42 and define

$f(\eta)\equiv\beta-\frac{1}{16}(\alpha-\gamma)^{2}+2\geq 0\,.$ (43)

Using the definitions of Eqs. 29 and 32, we can recast everything in terms of
$\cos 2\eta$ and $\sin 2\eta$ such that $f(\eta)$ only depends on these
quantities. We can then easily determine the extremal value of
$\eta_{\text{min}}$ which gives the minimal $f_{\text{min}}$. After some
algebra, the positivity condition $f_{\text{min}}\geq 0$ reads:

$2(\lambda_{1}\lambda_{2}+\sqrt{\lambda_{1}\lambda_{2}}(\lambda_{3}+\lambda_{4}))-\frac{1}{2}\big{|}\tilde{\lambda}_{6}-\tilde{\lambda}_{7}\big{|}^{2}-\big{|}2\sqrt{\lambda_{1}\lambda_{2}}\,\lambda_{5}-\frac{1}{2}(\tilde{\lambda}_{6}-\tilde{\lambda}_{7})^{2}\big{|}\geq
0\,,$ (44)

where we have defined the rescaled couplings:

$\tilde{\lambda}_{6}\equiv\left(\frac{\lambda_{2}}{\lambda_{1}}\right)^{1/4}\lambda_{6}\,,\,\,\,\tilde{\lambda}_{7}\equiv\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{1/4}\lambda_{7}\,.$
(45)

Eqs. (38), (39), and (44) present simplified necessary conditions for BFB and
are the main result of this section.

###  5.2 Sufficient conditions for boundedness from below

It is also useful to have simplified sufficient conditions which allow one to
quickly determine if the potential is BFB for a given parameter set. Consider
the top branch of Eq. 34 — i.e. $\Delta\geq 0$ and $-2\leq\beta\leq 6$ and
$\chi_{1}\leq 0$. One can show666Note that $\alpha=-\frac{1}{2}(\beta+2)$ and
$\gamma=-\frac{1}{2}(\beta+2)$ are the directions along which $\Delta=0$. In
order to satisfy the $\chi_{1}\leq 0$ condition, relevant for $\beta\leq 6$,
we must have $\alpha+\frac{1}{2}(\beta+2)\geq 0$ and
$\gamma+\frac{1}{2}(\beta+2)\geq 0$. So long as $\beta\geq-2$, these
conditions combined will always yield positive $\Delta$. Then Eq. 46 is a
sufficient condition that follows from the top branch of Eq. 34. that a
stronger condition (which will lead to a sufficient condition) is $\beta\leq
6$ and $\alpha+\frac{1}{2}(\beta+2)>0$ and $\gamma+\frac{1}{2}(\beta+2)>0$. In
terms of the $\lambda$’s, this translates to the sufficient condition:

$\begin{split}&3\sqrt{\lambda_{1}\lambda_{2}}-(\lambda_{3}+|\lambda_{4}|+|\lambda_{5}|)\geq
0\,,\\\
&\text{and}\,\,\,\sqrt{\lambda_{1}\lambda_{2}}+\lambda_{3}-(|\lambda_{4}|+|\lambda_{5}|+4|\tilde{\lambda}_{6}|)>0\,,\\\
&\text{and}\,\,\,\sqrt{\lambda_{1}\lambda_{2}}+\lambda_{3}-(|\lambda_{4}|+|\lambda_{5}|+4|\tilde{\lambda}_{7}|)>0\,.\end{split}$
(46)

Now consider the bottom branch of Eq. 34 — i.e. $\Delta\geq 0$ and $\beta>6$
and $\chi_{2}\leq 0$. To arrive at an analytic sufficient condition, consider
the stronger bound $\beta>6$ and $\alpha+2\sqrt{\beta-2}>0$ and
$\gamma+2\sqrt{\beta-2}>0$. In terms of the potential parameters, this
condition reads:

$\begin{split}&\lambda_{3}-(3\sqrt{\lambda_{1}\lambda_{2}}+|\lambda_{4}|+|\lambda_{5}|)\geq
0\,,\\\
&\text{and}\,\,\,\sqrt{\lambda_{3}-(\sqrt{\lambda_{1}\lambda_{2}}+|\lambda_{4}|+|\lambda_{5}|)}-\frac{\sqrt{2}}{(\lambda_{1}\lambda_{2})^{1/4}}|\tilde{\lambda}_{6}|>0\,,\\\
&\text{and}\,\,\,\sqrt{\lambda_{3}-(\sqrt{\lambda_{1}\lambda_{2}}+|\lambda_{4}|+|\lambda_{5}|)}-\frac{\sqrt{2}}{(\lambda_{1}\lambda_{2})^{1/4}}|\tilde{\lambda}_{7}|>0\,.\end{split}$
(47)

Eqs. (46) and (47) are the main result of this section.

###  5.3 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 $\\{\lambda_{1},...\lambda_{7}\\}$. We take
as allowed ranges $\lambda_{1,2}\in[0,2\pi]$ and
$|\lambda_{3,4,5,6,7}|\leq\frac{\pi}{2}$, with $\lambda_{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.

Figure 3: The white circle represents 10,000 randomly chosen points in the
7-dimensional parameter space of couplings $\\{\lambda_{1},...\lambda_{7}\\}$.
We take as priors $\lambda_{1,2}\in[0,2\pi]$ and
$|\lambda_{3,4,5,6,7}|\leq\frac{\pi}{2}$, with $\lambda_{5,6,7}$ allowed to be
complex. The red circle encompasses the points which pass our analytic
necessary conditions of Eqs. (38), (39), and (44); the black circle contains
the points which pass the necessary and sufficient BFB condition of Eq. 34;
and the innermost blue circle contains the points which pass our sufficient
condition of Eq. 46.

An examination of the analytical expressions indicates that the quantity
$\sqrt{\lambda_{1}\lambda_{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
$\lambda_{i}$, we plot a histogram in $\sqrt{\lambda_{1}\lambda_{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 $\lambda_{1,2}\in[0,2\pi]$ and
$|\lambda_{3,4,5,6,7}|\leq\pi/2$ with $\lambda_{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 $\sqrt{\lambda_{1}\lambda_{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 $\sqrt{\lambda_{1}\lambda_{2}}$.

Figure 4: Histogram of $\sqrt{\lambda_{1}\lambda_{2}}$ displaying the fraction
of tested points per bin which pass the necessary conditions Eqs. (44, 38, 39)
(red), numerical test Eq. 34 (black), and sufficient condition Eq. 46 (blue).
As in Fig. 3, we take as priors $\lambda_{1,2}\in[0,2\pi]$ and
$|\lambda_{3,4,5,6,7}|\in[0,\frac{\pi}{2}]$, with $\lambda_{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 $\lambda_{6,7}\neq
0$) do exist. For example, the authors of Branco:2011iw find as a necessary
condition:

$\displaystyle\frac{1}{2}(\lambda_{1}+\lambda_{2})+\lambda_{3}+\lambda_{4}+\lambda_{5}-2|\lambda_{6}+\lambda_{7}|>0\,.$
(48)

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 $\lambda_{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.

##  6 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 Barroso:2013awa,Ivanov:2015nea 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:

$V=-M_{\mu}\mathsf{r}^{\mu}+\frac{1}{2}\Lambda_{\mu\nu}\mathsf{r}^{\mu}\mathsf{r}^{\nu}-\frac{1}{2}\zeta\mathsf{r}^{\mu}\mathsf{r}_{\mu}\,,$
(49)

where $M_{\mu}$ encodes the mass terms:

$M_{\mu}=\begin{pmatrix}-\frac{1}{2}(m_{11}^{2}+m_{22}^{2}),&(m_{12}^{2})^{R},&-(m_{12}^{2})^{I},&-\frac{1}{2}(m_{11}^{2}-m_{22}^{2})\end{pmatrix}\,,$
(50)

$\mathsf{r}^{\mu}$ is a vector of field bilinears:

$\mathsf{r}^{\mu}=\begin{pmatrix}|\Phi_{1}|^{2}+|\Phi_{2}|^{2},&2\text{Re}[\Phi_{1}^{\dagger}\Phi_{2}],&2\text{Im}[\Phi_{1}^{\dagger}\Phi_{2}],&|\Phi_{1}|^{2}-|\Phi_{2}|^{2}\end{pmatrix}\,,$
(51)

and $\Lambda_{\mu\nu}$ encodes the quartic terms:

$\Lambda_{\mu\nu}=\frac{1}{2}\begin{pmatrix}\frac{1}{2}(\lambda_{1}+\lambda_{2})+\lambda_{3}&\lambda_{6}^{R}+\lambda_{7}^{R}&-(\lambda_{6}^{I}+\lambda_{7}^{I})&\frac{1}{2}(\lambda_{1}-\lambda_{2})\\\
(\lambda_{6}^{R}+\lambda_{7}^{R})&(\lambda_{4}+\lambda_{5}^{R})&-\lambda_{5}^{I}&\lambda_{6}^{R}-\lambda_{7}^{R}\\\
-(\lambda_{6}^{I}+\lambda_{7}^{I})&-\lambda_{5}^{I}&\lambda_{4}-\lambda_{5}^{R}&-(\lambda_{6}^{I}-\lambda_{7}^{I})\\\
\frac{1}{2}(\lambda_{1}-\lambda_{2})&\lambda_{6}^{R}-\lambda_{7}^{R}&-(\lambda_{6}^{I}-\lambda_{7}^{I})&\frac{1}{2}(\lambda_{1}+\lambda_{2})-\lambda_{3}\end{pmatrix}\,.$
(52)

The last term in Eq. 49 is a Lagrange multiplier we have introduced to enforce
the condition $\mathsf{r}^{\mu}\mathsf{r}_{\mu}=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 Barroso:2013awa,Ferreira:2004yd
for more details). In the above equations, indices are raised and lowered
using a Minkowski metric.

Provided the matrix $\Lambda_{\mu\nu}$, 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:

$\Lambda_{\mu\nu}^{\text{diag}}=\begin{pmatrix}\Lambda_{0}&0&0&0\\\
0&-\Lambda_{1}&0&0\\\ 0&0&-\Lambda_{2}&0\\\
0&0&0&-\Lambda_{3}\end{pmatrix}\,,$ (53)

with $\Lambda_{0}$ the “timelike” eigenvalue and $\Lambda_{i}$ “spacelike”.
Let us define the “signature matrix” $S$ as $S\equiv\Lambda_{\mu\nu}-\zeta
g_{\mu\nu}$. In diagonal form, it looks like:

$S=\begin{pmatrix}\Lambda_{0}-\zeta&0&0&0\\\ 0&\zeta-\Lambda_{1}&0&0\\\
0&0&\zeta-\Lambda_{2}&0\\\ 0&0&0&\zeta-\Lambda_{3}\end{pmatrix}\,.$ (54)

The discriminant is generically given by the determinant of the signature
matrix:

$D=\det S\,.$ (55)

By using the diagonal form above, we can write this as:

$D=(\Lambda_{0}-\zeta)(\zeta-\Lambda_{1})(\zeta-\Lambda_{2})(\zeta-\Lambda_{3})\,.$
(56)

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 $\Lambda_{0}$, and Lagrange multiplier $\zeta$.

$\text{We are in a global minimum if and only
if}:\begin{cases}D>0\,,\,\text{or}\\\
D<0\,\,\,\text{and}\,\,\,\zeta>\Lambda_{0}\,.\end{cases}$ (57)

For our purposes, it is more useful to work with the “Euclideanized" version
of $\Lambda_{\mu\nu}$ obtained by lowering one of the indices with the
Minkowski metric, $\Lambda_{E}\equiv\Lambda^{\mu}_{\nu}$. Explicitly:

$\Lambda_{E}=\frac{1}{2}\begin{pmatrix}\frac{1}{2}(\lambda_{1}+\lambda_{2})+\lambda_{3}&\lambda_{6}^{R}+\lambda_{7}^{R}&-(\lambda_{6}^{I}+\lambda_{7}^{I})&\frac{1}{2}(\lambda_{1}-\lambda_{2})\\\
-(\lambda_{6}^{R}+\lambda_{7}^{R})&-(\lambda_{4}+\lambda_{5}^{R})&\lambda_{5}^{I}&-(\lambda_{6}^{R}-\lambda_{7}^{R})\\\
\lambda_{6}^{I}+\lambda_{7}^{I}&\lambda_{5}^{I}&-(\lambda_{4}-\lambda_{5}^{R})&\lambda_{6}^{I}-\lambda_{7}^{I}\\\
-\frac{1}{2}(\lambda_{1}-\lambda_{2})&-(\lambda_{6}^{R}-\lambda_{7}^{R})&\lambda_{6}^{I}-\lambda_{7}^{I}&-\frac{1}{2}(\lambda_{1}+\lambda_{2})+\lambda_{3}\end{pmatrix}\,.$
(58)

In terms of $\Lambda_{E}$, the discriminant is:

$D=-\det[\Lambda_{E}-\mathbb{1}\zeta]\,.$ (59)

The other quantity necessary for formulating the discriminant is the Lagrange
multiplier $\zeta$. This may be obtained by looking at any component of the
minimization condition:

$\Lambda^{\mu}_{\nu}r^{\nu}-M^{\mu}=\zeta r^{\mu}\,.$ (60)

We parameterize the vacuum expectation values (vevs) of the doublets as:

$\expectationvalue{\Phi_{1}}=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\
v_{1}\end{pmatrix}\,,\,\,\,\expectationvalue{\Phi_{2}}=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\
v_{2}e^{i\eta}\end{pmatrix}\,.$ (61)

Then the expectation value of field bilinears
$r^{\mu}\equiv\expectationvalue{\mathsf{r}^{\mu}}$ is:

$r^{\mu}=\begin{pmatrix}\frac{1}{2}(v_{1}^{2}+v_{2}^{2}),&v_{1}v_{2}\cos\eta,&v_{1}v_{2}\sin\eta,&\frac{1}{2}(v_{1}^{2}-v_{2}^{2})\end{pmatrix}\,.$
(62)

The expression for $\zeta$ is particularly simple if we choose the “1”
component. In particular if we take $\eta=0$, then:

$\zeta=\frac{(m_{12}^{2})^{R}}{v_{1}v_{2}}-\frac{1}{2}\left(\frac{v_{1}}{v_{2}}\lambda_{6}^{R}+\frac{v_{2}}{v_{1}}\lambda_{7}^{R}+(\lambda_{4}+\lambda_{5}^{R})\right)\,.$
(63)

Note that this has the interpretation of the charged Higgs mass over the vev
$v$ squared,

$\zeta=\frac{M_{H^{\pm}}^{2}}{v^{2}}\,,$ (64)

as first demonstrated in Ivanov:2006yq,Maniatis:2006fs. 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
$\Lambda_{0}$ with $\zeta$: we are in a global minimum if $\zeta>\Lambda_{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, $\Lambda^{\mu\nu}$ is positive definite
and $\Lambda_{0}>\Lambda_{1,2,3}$. Then from Eq. (56), $D>0$ necessarily
implies that we have the ordering $\Lambda_{0}>\zeta>\Lambda_{1,2,3}$.
Similarly $D<0$ and $\zeta>\Lambda_{0}$ necessarily implies the ordering
$\zeta>\Lambda_{0}>\Lambda_{1,2,3}$. So we see that the relative ordering of
$\Lambda_{0}$ and $\zeta$ does not actually matter$-$all that matters for a
potential which has been verified to be BFB is that $\zeta$ be larger than the
spatial eigenvalues, $\zeta>\Lambda_{1,2,3}$.

###  6.1 Sufficient conditions for stability

####  6.1.1 Gershgorin bounds

As in Section 4.3, we can bound the eigenvalues of $\Lambda_{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 $\Lambda_{E}$ and define the endpoint of each interval as
$\Gamma_{i}\equiv a_{ii}+R_{i}$, with $R_{i}=\sum_{j\neq i}|a_{ij}|$:

$\displaystyle\Gamma_{1}$
$\displaystyle=\frac{1}{4}(\lambda_{1}+\lambda_{2})+\frac{\lambda_{3}}{2}+\frac{1}{2}\left(|\lambda_{6}^{R}+\lambda_{7}^{R}|+|\lambda_{6}^{I}+\lambda_{7}^{I}|+\frac{1}{2}|\lambda_{1}-\lambda_{2}|\right)\,,$
(65a) $\displaystyle\Gamma_{2}$
$\displaystyle=-\frac{1}{4}(\lambda_{1}+\lambda_{2})+\frac{\lambda_{3}}{2}+\frac{1}{2}\left(|\lambda_{6}^{R}-\lambda_{7}^{R}|+|\lambda_{6}^{I}-\lambda_{7}^{I}|+\frac{1}{2}|\lambda_{1}-\lambda_{2}|\right)\,,$
(65b) $\displaystyle\Gamma_{3}$
$\displaystyle=-\frac{1}{2}(\lambda_{4}+\lambda_{5}^{R})+\frac{1}{2}\bigg{(}|\lambda_{5}^{I}|+|\lambda_{6}^{R}+\lambda_{7}^{R}|+|\lambda_{6}^{R}-\lambda_{7}^{R}|\bigg{)}\,,$
(65c) $\displaystyle\Gamma_{4}$
$\displaystyle=-\frac{1}{2}(\lambda_{4}-\lambda_{5}^{R})+\frac{1}{2}\bigg{(}|\lambda_{5}^{I}|+|\lambda_{6}^{I}+\lambda_{7}^{I}|+|\lambda_{6}^{I}-\lambda_{7}^{I}|\bigg{)}\,.$
(65d)

We know that all eigenvalues must be less than the endpoint of the interval
extending the furthest in the $+\hat{x}$ direction,

$\text{max}[\Gamma_{i}]\geq\Lambda_{0,1,2,3}\,.$ (66)

Meanwhile, an extremum will be the global minimum if $\zeta>\Lambda_{1,2,3}$.
Thus, it is sufficient to demand:

$\zeta>\text{max}[\Gamma_{i}]\,.$ (67)

####  6.1.2 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
$\zeta$, which gives the constraint:

$\displaystyle\zeta>\frac{1}{2}\sqrt{\lambda_{1}^{2}+\lambda_{2}^{2}+2(\lambda_{3}^{2}+\lambda_{4}^{2}+|\lambda_{5}|^{2})+4(|\lambda_{6}|^{2}+|\lambda_{7}|^{2})}$
(68)

Note that in this case, the Frobenius bound is insensitive to the signs of the
$\lambda_{i}$, while the Gershgorin condition is sensitive to the signs of
$\lambda_{3}$, $\lambda_{4}$, and $\lambda_{5}^{R}$. We thus expect the
Gershgorin bound to be the stronger of the two.

####  6.1.3 Principal minors

In the case of a non-symmetric matrix such as $\Lambda_{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 $\Lambda_{E}$ to obtain a
sufficient condition.

The symmetric part of $\Lambda_{E}$ is given by

$\displaystyle\Lambda_{E}^{S}$
$\displaystyle=\frac{1}{2}(\Lambda_{E}^{T}+\Lambda_{E})$
$\displaystyle=\frac{1}{2}\begin{pmatrix}\frac{1}{2}(\lambda_{1}+\lambda_{2})+\lambda_{3}&0&0&0\\\
0&-(\lambda_{4}+\lambda_{5}^{R})&\lambda_{5}^{I}&-(\lambda_{6}^{R}-\lambda_{7}^{R})\\\
0&\lambda_{5}^{I}&-(\lambda_{4}-\lambda_{5}^{R})&\lambda_{6}^{I}-\lambda_{7}^{I}\\\
0&-(\lambda_{6}^{R}-\lambda_{7}^{R})&\lambda_{6}^{I}-\lambda_{7}^{I}&-\frac{1}{2}(\lambda_{1}+\lambda_{2})+\lambda_{3}\end{pmatrix}$
(69)

We require the matrix $\zeta\mathbb{1}-\Lambda_{E}^{S}$ to be positive-
definite. Since the lower-right 3$\times$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$\times$3 submatrix. This gives the
following set of conditions:

$\zeta-\frac{1}{4}(\lambda_{1}+\lambda_{2})-\frac{1}{2}\lambda_{3}>0\,,$ (70a)
$\zeta+\frac{1}{2}\lambda_{4}+\frac{1}{2}\lambda_{5}^{R}>0\,,$ (70b)
$\big{(}\zeta+\frac{1}{2}\lambda_{4}\big{)}^{2}-\frac{1}{4}|\lambda_{5}|^{2}>0\,,$
(70c)
$\begin{split}(4\zeta+\lambda_{1}+\lambda_{2}&-2\lambda_{3})((2\zeta+\lambda_{4})^{2}-|\lambda_{5}|^{2})\\\
&+\lambda_{5}(\lambda_{6}^{*}-\lambda_{7}^{*})^{2}+\lambda_{5}^{*}(\lambda_{6}-\lambda_{7})^{2}-2(2\zeta+\lambda_{4})|\lambda_{6}-\lambda_{7}|^{2}>0\,.\end{split}$
(70d)

Taken together, the Eqs. (70) provide a sufficient condition for vacuum
stability.

###  6.2 Numerical comparison

Figure 5: Comparison of the fraction of points that pass the exact stability
conditions, Eq. (57) (black dots), with respect to the fraction passing the
three sufficient conditions for vacuum stability: principal minors Eq. 70
(blue), Gershgorin disk theorem Eq. 67 (light blue), and Frobenius norm Eq. 68
(purple). We plot the fraction of points that pass each condition as a
function of $M_{H^{\pm}}$. We take $\lambda_{1,2}\in[0,2\pi]$ and
$|\lambda_{3,4,5,6,7}|\leq\frac{\pi}{2}$, with $\lambda_{5,6,7}$ allowed to be
complex, and restrict to examining points which are BFB.

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^{\pm}}$. We compare these results with the fraction that pass the
exact stability condition, Eq. 57. As in previous sections, we choose the
$\lambda_{i}$ randomly with $\lambda_{1,2}\in[0,2\pi]$ and
$|\lambda_{3,4,5,6,7}|\leq\frac{\pi}{2}$, with $\lambda_{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 $\lambda_{i}$ are small compared to the ratio
$M_{H^{\pm}}^{2}/v^{2}$.

###  6.3 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 ${\cal CP}$-even Higgs bosons $hhH$, etc.
(see e.g. Carena:2015moc 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 $\eta=0$), the scalar
which obtains a vev is denoted by $\phi_{1}^{0}$. The mass matrix for the
neutral scalars $\vec{\phi}=(\phi_{1}^{0},\phi_{2}^{0},a_{0})^{\text{T}}$
reads:

$\mathcal{M}^{2}=v^{2}\begin{pmatrix}Z_{1}&Z_{6}^{R}&-Z_{6}^{I}\\\
Z_{6}^{R}&\frac{M_{H^{\pm}}^{2}}{v^{2}}+\frac{1}{2}(Z_{4}+Z_{5}^{R})&-\frac{1}{2}Z_{5}^{I}\\\
-Z_{6}^{I}&-\frac{1}{2}Z_{5}^{I}&\frac{M_{H^{\pm}}^{2}}{v^{2}}+\frac{1}{2}(Z_{4}-Z_{5}^{R})\end{pmatrix}\,,$
(71)

where $M_{H^{\pm}}^{2}$ is the charged Higgs mass:

$M_{H^{\pm}}^{2}=M_{22}^{2}+\frac{1}{2}Z_{3}v^{2}\,.$ (72)

We will restrict ourselves to the alignment limit, which is the limit in which
$\phi_{1}^{0}$ 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 taking $M_{H^{\pm}}^{2}+\frac{1}{2}(Z_{4}\pm Z_{5}^{R})v^{2}\gg
Z_{1}v^{2}$. Under this limit, the heavy mass eigenstates $h_{2}$ and $h_{3}$
and the heavy charged Higgs $H^{\pm}$ decouple from the light mass eigenstate,
leaving $h_{1}$ aligned with $\phi_{1}^{0}$. 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}|\ll 1$, for which the mixing between
$\phi_{1}^{0}$ and the other neutral scalars vanishes, leading to the
identification of $\phi_{1}^{0}$ with the mass eigenstate $h_{1}$. For the
following discussion we will take $|Z_{6}|\ll 1$ and work in the alignment
without decoupling limit.

We define $h_{1}\equiv h$ to be the SM-like Higgs boson, which has a mass
given by

$M_{h}^{2}=Z_{1}v^{2}\,.$ (73)

To obtain a physical Higgs mass close to the experimental value of 125 GeV, it
is required that we fix $Z_{1}\approx 0.25$. The remaining 2$\times$2 mass
matrix can be diagonalized to obtain the masses of the remaining scalars
$h_{2}$ and $h_{3}$:

$M_{h_{3},h_{2}}^{2}=M_{H^{\pm}}^{2}+\frac{1}{2}(Z_{4}\pm|Z_{5}|)v^{2}\,.$
(74)

There are two possibilities for the ${\cal CP}$ properties of these states. So
long as $Z_{5}^{I}\neq 0$, $h_{2}$ and $h_{3}$ have mixed ${\cal CP}$
properties. In the limit of $Z_{5}^{I}=0$, meanwhile, the non-standard Higgs
mass matrix becomes diagonal, and we obtain mass eigenstates $H$ and $A$ with
definite ${\cal CP}$ character,

$M_{H,A}^{2}=M_{H^{\pm}}^{2}+\frac{1}{2}(Z_{4}\pm Z_{5}^{R})v^{2}\,.$ (75)

The masses of the general mass eigenstates and the states of definite ${\cal
CP}$ character can be related by

$\displaystyle M_{H}^{2}$ $\displaystyle=$ $\displaystyle
M_{h_{3},h_{2}}^{2}+\frac{1}{2}(Z_{5}^{R}\mp|Z_{5}|)v^{2}$ $\displaystyle
M_{A}^{2}$ $\displaystyle=$ $\displaystyle
M_{h_{3},h_{2}}^{2}-\frac{1}{2}(Z_{5}^{R}\pm|Z_{5}|)v^{2}.$ (76)

In the following analysis we will make no assumptions about the ${\cal 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 $\Gamma$’s in terms of physical
masses, a sufficient condition for vacuum stability becomes:

$\begin{split}&M_{h_{2}}^{2}>\frac{1}{2}(|Z_{5}^{I}|-Z_{5}^{R}-|Z_{5}|)v^{2}+\frac{1}{2}|Z_{7}^{R}|v^{2}\,\,\,\text{and}\\\
&M_{h_{3}}^{2}>\frac{1}{2}(|Z_{5}^{I}|+Z_{5}^{R}+|Z_{5}|)v^{2}+\frac{1}{2}|Z_{7}^{I}|v^{2}\,\,\,\,\,\text{and}\\\
&M_{H^{\pm}}^{2}>\frac{1}{2}\,\text{max}[M_{h}^{2},\,Z_{2}v^{2}]+\frac{1}{2}(Z_{3}+|Z_{7}^{R}|+|Z_{7}^{I}|)v^{2}\,.\end{split}$
(77)

Next, we can recast the Frobenius sufficient condition, Eq. 68, in terms of
the physical masses; doing so results in the following condition:

$2M_{H^{\pm}}^{2}\big{(}M_{h_{2}}^{2}+M_{h_{3}}^{2}\big{)}-\big{(}M_{H^{\pm}}^{4}+M_{h_{2}}^{4}+M_{h_{3}}^{4}+\frac{1}{4}M_{h}^{4}\big{)}>\frac{1}{4}(Z_{2}^{2}+2Z_{3}^{2}+4|Z_{7}|^{2})v^{4}\,.$
(78)

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:

$\begin{split}&4M_{H^{\pm}}^{2}-M_{h}^{2}>(Z_{2}+2Z_{3})v^{2},\hskip
21.33955pt\\\ &M_{h_{3}}^{2}>\frac{1}{2}(|Z_{5}|-Z_{5}^{R}),\hskip
113.81102pt\\\ &M_{h_{2}}^{2}M_{h_{3}}^{2}>0,\hskip 102.43008pt\\\
4M_{h_{2}}^{2}M_{h_{3}}^{2}(4M_{H^{\pm}}^{2}+M_{h}^{2}+Z_{2}v^{2}&-2Z_{3}v^{2})-2(M_{h_{2}}^{2}+M_{h_{3}}^{2})|Z_{7}|^{2}+Z_{5}{Z_{7}^{*}}^{2}+Z_{5}^{*}Z_{7}^{2}>0\,.\end{split}$
(79)

##  7 ${\cal 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 ${\cal CP}$ violation permitted in the
alignment limit. This possibility has largely been neglected in the many
previous studies which restrict themselves to the $\mathbb{Z}_{2}$-symmetric
2HDM. This is understandable since exact alignment implies ${\cal CP}$
conservation in the $\mathbb{Z}_{2}$-symmetric case. When working in the fully
general 2HDM, however, it is possible to have ${\cal 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_{12}^{2},Z_{5},Z_{6},Z_{7}\\}$. One of these is fixed by the
minimization condition $M_{12}^{2}=-\frac{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
${\cal 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 $J_{1}\sim\text{Im}[Z_{5}^{*}Z_{6}^{2}]$,
$J_{2}\sim\text{Im}[Z_{5}^{*}Z_{7}^{2}]$,
$J_{3}\sim\text{Im}[Z_{6}^{*}Z_{7}]$. It is then clear that the condition for
the Higgs sector to be ${\cal CP}$ invariant is:

$\text{Im}[Z_{5}^{*}Z_{6}^{2}]=\text{Im}[Z_{5}^{*}Z_{7}^{2}]=\text{Im}[Z_{6}^{*}Z_{7}]=0\,.$
(80)

There are two ways in which this can be satisfied [52]: either
$Z_{5}^{I}=Z_{6}^{R}=Z_{7}^{R}=0$ or $Z_{5}^{I}=Z_{6}^{I}=Z_{7}^{I}=0$. Note
that in the limit of exact alignment we have $Z_{6}=0$, so this reduces to
demanding either $Z_{5}^{I}=Z_{7}^{R}=0$ or $Z_{5}^{I}=Z_{7}^{I}=0$.

Meanwhile in the 2HDM with a (softly broken) $\mathbb{Z}_{2}$ symmetry, the
fact that $\lambda_{6}=\lambda_{7}=0$ implies the following two relations
between the parameters in the Higgs basis (see Eq. (88) in Appendix A) [53,
54]

$\displaystyle Z_{6}+Z_{7}$ $\displaystyle=\frac{1}{2}\tan
2\beta(Z_{2}-Z_{1}),$ (81) $\displaystyle Z_{6}-Z_{7}^{*}$
$\displaystyle=\frac{1}{\tan 2\beta}(Z_{1}+2Z_{6}\cot 2\beta-
Z_{3}-Z_{4}-Z_{5}).$ (82)

It immediately follows that:

$\displaystyle Z_{7}^{I}$ $\displaystyle=-Z_{6}^{I},$ (83) $\displaystyle
Z_{5}^{I}$ $\displaystyle=2\;\frac{1-\tan^{2}2\beta}{\tan 2\beta}Z_{6}^{I}.$
(84)

These conditions imply that in the exact alignment limit (i.e. $Z_{6}=0$), it
will necessarily be the case that $Z_{5}^{I}=Z_{7}^{I}=0$. Thus, exact
alignment directly leads to ${\cal CP}$ conservation in the
$\mathbb{Z}_{2}$-symmetric or softly broken $\mathbb{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}|\gtrsim
0$), |$Z_{5}^{I}$|, which controls the mixing between the $H$ and $A$ bosons
(see Eq. 71), can still be large for large $\tan\beta$.

Figure 6: Upper left: Parameter scan of the 2HDM with a softly broken
$\mathbb{Z}_{2}$ symmetry in the $(M_{h_{3}}-M_{h_{2}},\tan\beta)$ parameter
plane with parameter range $\lambda_{1,2}\in[0,2\pi]$ and
$|\lambda_{3,4,5}|\in[0,\pi]$. The conditions $|\sin\theta_{12}|<0.1$ and
$|\sin\theta_{13}|<0.01$ are imposed. The colour indicates the minimal value
of $|\theta_{23}/\pi-1/4|$ in each hexagonal patch. Upper right: Same as upper
left, but constraints from perturbative unitary, BFB, and vacuum stability are
applied in addition. Lower left: Same as upper left, but the scan is performed
in the general 2HDM without a (softly broken) $\mathbb{Z}_{2}$ symmetry
($|\lambda_{6,7}|\in[0,\pi]$). Lower right: Same as lower left, but
constraints from perturbative unitary, BFB, and vacuum stability are applied
in addition.

The physical consequences of the difference in the ${\cal CP}$ properties in
the alignment limit between the (softly broken) $\mathbb{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. ACME:2018yjb,Altmannshofer:2020shb), 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\theta_{12}|<0.1$ and
$|\sin\theta_{13}|<0.01$, where $\theta_{12}$ and $\theta_{13}$ are the
respective mixing angles.777The stated bounds on $\theta_{12}$ and
$\theta_{13}$ are just illustrative. We checked that our conclusions do not
depend strongly on the values inserted. The colour code in the figure
indicates the minimum value of $|\theta_{23}/\pi-1/4|$ in each hexagonal
patch, where $\theta_{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
${\cal 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 $\mathbb{Z}_{2}$ symmetry in the
$(M_{h_{3}}-M_{h_{2}},\tan\beta)$ parameter plane. We observe that for large
$\tan\beta$, a large mixing between the neutral BSM Higgs bosons can be
realized if their mass difference is below $\sim 70\;\text{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 $\frac{1}{2}Z_{5}^{I}v^{2}$ term. For lower $\tan\beta$, the
condition $|\sin\theta_{13}|<0.01$ directly implies that $Z_{5}^{I}$ 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
$\mathbb{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)
$\mathbb{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}}\gtrsim 100\;\text{GeV}$.

Finally, we want to remark that ${\cal CP}$ violation can become manifest not
only in the neutral mass matrix but also in the bosonic couplings. This occurs
if either $\tilde{Z}_{7}^{R}\neq 0$ or $\tilde{Z}_{7}^{I}\neq 0$, since these
couplings enter into couplings like $g_{h_{1}h_{2}h_{3}}$ [52]. Exotic decays
like $h_{3}\rightarrow h_{1}h_{2}\rightarrow 3h_{1}$ would then be indicative
of ${\cal CP}$ violation in the bosonic sector (see e.g. [41]).

##  8 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
$\mathbb{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
flavor-changing 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.

Conditions | Perturbative unitarity | Bounded from below | Vacuum stability
---|---|---|---
Exact | Eq. 12 | Eq. 34 | Eq. 57
Necessary |  Eqs. 15 and 18 | Eq. 39 | —
Sufficient | Eq. 26 | Eq. 46 | Eq. 70
Table 1: Overview of the primary results of this paper; further constraints
and their analysis may be found in the main text.

Based on this motivation, in this work we present a step towards a systematic
exploration of the non-$\mathbb{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
${\cal CP}$-violating mixing between the BSM Higgs bosons can be realized in
the general 2HDM in comparison to a 2HDM with a (softly broken)
$\mathbb{Z}_{2}$ symmetry. While we found that large ${\cal CP}$-violating
mixing can only be realized for large $\tan\beta$ in the 2HDM with a softly
broken $\mathbb{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 $\mathbb{Z}_{2}$ symmetry for future work.

## Acknowledgments

C.W. would like to thank H. Haber and P. Ferreira for useful discussions and
comments. Fermilab is operated by Fermi Research Alliance, LLC under contract
number DE-AC02-07CH11359 with the United States Department of Energy. M.C. and
C.W. would like to thank the Aspen Center for Physics, which is supported by
National Science Foundation grant No. PHY-1607611, where part of this work has
been done. C.W. has been partially supported by the U.S. Department of Energy
under contracts No. DEAC02- 06CH11357 at Argonne National Laboratory. The work
of C.W. and N.C. at the University of Chicago has also been supported by the
DOE grant DE-SC0013642. H.B. acknowledges support by the Alexander von
Humboldt foundation. A.I. is supported by the U.S. Department of Energy,
Office of Science, Office of Workforce Development for Teachers and
Scientists, Office of Science Graduate Student Research (SCGSR) program, which
is administered by the Oak Ridge Institute for Science and Education (ORISE)
for the DOE. ORISE is managed by ORAU under contract number DE-SC0014664.

##  Appendix A Higgs basis conversion

The phenomenological properties of the Higgs sector are more easily analyzed
in the Higgs basis, in which only one of the doublets possesses a vev888This
is technically not enough to uniquely define the Higgs basis. The $U(1)$
diagonal subgroup of the $SU(2)$ symmetry in Higgs flavor space remains intact
following SSB. This corresponds to transformations $\Phi_{1}\rightarrow
e^{i\chi}\Phi_{1}$, $\Phi_{2}\rightarrow e^{-i\chi}\Phi_{2}$. As a result, we
have a one-dimensional family of Higgs bases parameterized by $\chi$:
$\\{e^{-i\chi}H_{1},e^{i\chi}H_{2}\\}$.. We parameterize the doublets as:

$H_{1}=\begin{pmatrix}G^{+}\\\
\frac{1}{\sqrt{2}}(v+\phi_{1}^{0}+iG^{0})\end{pmatrix}\,,\,\,\,H_{2}=\begin{pmatrix}H^{+}\\\
\frac{1}{\sqrt{2}}(\phi_{2}^{0}+ia^{0})\end{pmatrix}\,,$ (85)

where $G^{\pm}$ and $G^{0}$ are the Goldstones that become the longitudinal
components of $W^{\pm}$ and $Z$, $H^{\pm}$ is the physical singly charged
scalar state, and $(\phi_{1}^{0},\phi_{2}^{0},a^{0})$ are the neutral scalars.
The potential in the Higgs basis reads:

$\begin{split}V&=M_{11}^{2}H_{1}^{\dagger}H_{1}+M_{22}^{2}H_{2}^{\dagger}H_{2}-(M_{12}^{2}H_{1}^{\dagger}H_{2}+h.c.)\\\
&+\frac{1}{2}Z_{1}(H_{1}^{\dagger}H_{1})^{2}+\frac{1}{2}Z_{2}(H_{2}^{\dagger}H_{2})^{2}+Z_{3}(H_{1}^{\dagger}H_{1})(H_{2}^{\dagger}H_{2})+Z_{4}(H_{1}^{\dagger}H_{2})(H_{2}^{\dagger}H_{1})\\\
&+\left[\frac{1}{2}Z_{5}(H_{1}^{\dagger}H_{2})^{2}+Z_{6}(H_{1}^{\dagger}H_{1})(H_{1}^{\dagger}H_{2})+Z_{7}(H_{2}^{\dagger}H_{2})(H_{1}^{\dagger}H_{2})+h.c.\right]\,.\end{split}$
(86)

The conversion between the potential parameters in the weak eigenstate basis
and those in the Higgs basis have been worked out in Davidson:2005cw; to be
self-contained, we reproduce them here. They are obtained by a rotation by an
angle $\beta$ in field space of the original two Higgs doublets. The mass
terms in the two bases are related as:

$\displaystyle m_{11}^{2}$
$\displaystyle=M_{11}^{2}c_{\beta}^{2}+M_{22}^{2}s_{\beta}^{2}+\text{Re}[M_{12}^{2}e^{i\eta}]s_{2\beta}\,,$
(87a) $\displaystyle m_{22}^{2}$
$\displaystyle=M_{11}^{2}s_{\beta}^{2}+M_{22}^{2}c_{\beta}^{2}-\text{Re}[M_{12}^{2}e^{i\eta}]s_{2\beta}\,,$
(87b) $\displaystyle m_{12}^{2}e^{i\eta}$
$\displaystyle=\frac{1}{2}(M_{22}^{2}-M_{11}^{2})s_{2\beta}+\text{Re}[M_{12}^{2}e^{i\eta}]c_{2\beta}+i\,\text{Im}[M_{12}^{2}e^{i\eta}]\,,$
(87c)

where $\tan\beta=v_{2}/v_{1}$ with range $0\leq\beta\leq\frac{\pi}{2}$, and
$\eta$ is the phase accompanying $v_{2}$ in the general basis parameterization
of the doublets in Eq. (61). The relations between the quartic couplings are:

$\displaystyle\lambda_{1}$
$\displaystyle=Z_{1}c_{\beta}^{4}+Z_{2}s_{\beta}^{4}+\frac{1}{2}Z_{345}s_{2\beta}^{2}-2s_{2\beta}\left(\text{Re}[Z_{6}e^{i\eta}]c_{\beta}^{2}+\text{Re}[Z_{7}e^{i\eta}]s_{\beta}^{2}\right)\,,$
(88a) $\displaystyle\lambda_{2}$
$\displaystyle=Z_{1}s_{\beta}^{4}+Z_{2}c_{\beta}^{4}+\frac{1}{2}Z_{345}s_{2\beta}^{2}+2s_{2\beta}\left(\text{Re}[Z_{6}e^{i\eta}]s_{\beta}^{2}+\text{Re}[Z_{7}e^{i\eta}]c_{\beta}^{2}\right)\,,$
(88b) $\displaystyle\lambda_{3}$
$\displaystyle=\frac{1}{4}\left(Z_{1}+Z_{2}-2Z_{345}\right)s_{2\beta}^{2}+Z_{3}+\text{Re}[(Z_{6}-Z_{7})e^{i\eta}]s_{2\beta}c_{2\beta}\,,$
(88c) $\displaystyle\lambda_{4}$
$\displaystyle=\frac{1}{4}\left(Z_{1}+Z_{2}-2Z_{345}\right)s_{2\beta}^{2}+Z_{4}+\text{Re}[(Z_{6}-Z_{7})e^{i\eta}]s_{2\beta}c_{2\beta}\,,$
(88d) $\displaystyle\lambda_{5}e^{2i\eta}$
$\displaystyle=\frac{1}{4}(Z_{1}+Z_{2}-2Z_{345})s_{2\beta}^{2}+\text{Re}[Z_{5}e^{2i\eta}]+i\,\text{Im}[Z_{5}e^{2i\eta}]c_{2\beta}$
(88e) $\displaystyle\hskip
28.45274pt+\text{Re}[(Z_{6}-Z_{7})e^{i\eta}]s_{2\beta}c_{2\beta}+i\,\text{Im}[(Z_{6}-Z_{7})e^{i\eta}]s_{2\beta}\,,$
$\displaystyle\lambda_{6}e^{i\eta}$
$\displaystyle=\frac{1}{2}(Z_{1}c_{\beta}^{2}-Z_{2}s_{\beta}^{2}-Z_{345}c_{2\beta}-i\,\text{Im}[Z_{5}e^{2i\eta}])s_{2\beta}$
(88f) $\displaystyle\hskip
28.45274pt+\text{Re}[Z_{6}e^{i\eta}]c_{\beta}c_{3\beta}+i\,\text{Im}[Z_{6}e^{i\eta}]c_{\beta}^{2}+\text{Re}[Z_{7}e^{i\eta}]s_{\beta}s_{3\beta}+i\,\text{Im}[Z_{7}e^{i\eta}]s_{\beta}^{2}\,,$
$\displaystyle\lambda_{7}e^{i\eta}$
$\displaystyle=\frac{1}{2}(Z_{1}s_{\beta}^{2}-Z_{2}c_{\beta}^{2}+Z_{345}c_{2\beta}+i\,\text{Im}[Z_{5}e^{2i\eta}])s_{2\beta}$
(88g) $\displaystyle\hskip
28.45274pt+\text{Re}[Z_{6}e^{i\eta}]s_{\beta}s_{3\beta}+i\,\text{Im}[Z_{6}e^{i\eta}]s_{\beta}^{2}+\text{Re}[Z_{7}e^{i\eta}]c_{\beta}c_{3\beta}+i\,\text{Im}[Z_{7}e^{i\eta}]c_{\beta}^{2}\,,$

where we have defined $Z_{345}\equiv(Z_{3}+Z_{4}+\text{Re}[Z_{5}e^{2i\eta}])$.
For the reverse conversion from the Higgs basis to the general basis, one can
perform the same series of identifications, but substituting
$\lambda_{i}\leftrightarrow Z_{i}$ and $\beta\leftrightarrow-\beta$.

##  Appendix B Connection between $\mathbf{\zeta}$ and $M_{H^{\pm}}^{2}$

In Ivanov:2006yq,Maniatis:2006fs, it was mathematically demonstrated that the
charged Higgs mass squared $m_{H^{\pm}}^{2}$ is related to the Lagrange
multiplier $\zeta$ as:

$\zeta=\frac{M_{H^{\pm}}^{2}}{v^{2}}\,.$ (89)

We explain here why this is the case in the general 2HDM. Let the constrained
function we want to extremize be:
$V=-M_{\mu}r^{\mu}+\frac{1}{2}\Lambda_{\mu\nu}r^{\mu}r^{\nu}-\frac{\zeta}{2}(r^{\mu}r_{\mu}-2c)$,
with $M^{\mu}$, $r^{\mu}$, and $\Lambda_{\mu\nu}$ defined in Eqs. (50), (51)
and (52), respectively. In the case presented in the main text, we have $c=0$,
corresponding to the fact that the constraint is $r^{\mu}r_{\mu}=0$ (i.e. no
charge breaking minima are allowed). However, we will for the moment retain
non-zero $c$ to see what happens when we allow the constraint to relax
slightly.

The claim is that $\zeta$ should be interpreted as the rate of change of the
extremal value999In this appendix stars are used to denote quantities
evaluated at an extremum. For convenience they are dropped in the body of the
text. $\overline{V}_{*}$ with respect to the constraint
$c=\frac{1}{2}r^{\mu}r_{\mu}$, $\zeta=\frac{d\overline{V}_{*}}{dc}$; that is,
it is the price one would pay to relax the constraint $r^{\mu}r_{\mu}=0$. To
see this, let $(r_{*}^{\mu},\zeta_{*})$ be the field configuration and
Lagrange multiplier value which extremizes $\overline{V}$:

$\nabla\overline{V}\big{|}_{r_{\star}^{\mu},\zeta_{\star}}=0\,,$ (90)

where $\nabla=\left(\frac{\partial}{\partial
r^{\mu}},\frac{\partial}{\partial\zeta}\right)$. Now if we elevate $c$ to a
variable, then the extremal values $(r_{*}^{\mu}(c),\zeta_{*}(c))$ and minimum
$\overline{V}_{*}(c)$ will vary with it. The total derivative of
$V_{*}(r_{*}^{\mu},\zeta_{*},c)$ with respect to $c$ is:

$\frac{d\overline{V}_{*}}{dc}=\frac{\partial\overline{V}_{*}}{\partial
r_{*}^{\mu}}\frac{dr_{*}^{\mu}}{dc}+\frac{\partial\overline{V}_{*}}{\partial\zeta_{*}}\frac{d\zeta_{*}}{dc}+\frac{\partial\overline{V}_{*}}{\partial
c}=\frac{\partial\bar{V}_{*}}{\partial c}=\zeta_{*}\,,$ (91)

where we have used the fact that at the extremal point
$\frac{\partial\overline{V}_{*}}{\partial r_{*}^{\mu}}=0$ and
$\frac{\partial\overline{V}_{*}}{\partial\zeta_{*}}=0$. Thus it is indeed the
case that the Lagrange multiplier tells us how the extremum (here the minimum
value of the potential) changes when we relax the constraint:

$\zeta_{*}=\frac{d\overline{V}_{*}}{dc}\,.$ (92)

In our case when we allow for $r^{\mu}r_{\mu}>0$, the vacuum location shifts
and becomes charge breaking, and the potential depth changes as well. How does
the charged Higgs mass evaluated in the neutral minimum compare with the
cost101010We know this is a cost rather than a gain because it was shown in
Ferreira:2004yd that the potential value in a charge breaking minimum is
always higher than that in the corresponding neutral minimum,
$V_{CB}-V_{N}>0$. of increasing $r^{\mu}r_{\mu}$? To see this, we will work
out the expression for $c$ in the case of a charge breaking minimum.

We parameterize111111The minimization conditions forbid both $w$ and the
complex phase $\eta$ from being non-zero (i.e. we cannot have a vacuum which
breaks both charge and CP), so we set $\eta=0$. the expectation values of the
doublets as:

$\expectationvalue{\Phi_{1}}=\frac{1}{\sqrt{2}}\begin{pmatrix}0\\\
v_{1}\end{pmatrix}\,,\,\,\,\expectationvalue{\Phi_{2}}=\frac{1}{\sqrt{2}}\begin{pmatrix}w\\\
v_{2}\end{pmatrix}\,.$ (93)

If we now expand $\frac{1}{2}r_{*}^{\mu}r_{*\,\mu}$ about this minimum and
examine the fluctuations in the charged Higgs fields $\phi_{1}^{\pm}$ and
$\phi_{2}^{\pm}$, after a bit of algebra one can show the quadratic part is:

$\begin{split}c\equiv\frac{1}{2}r^{\mu}r_{\mu}&=2|\Phi_{1}^{\dagger}\Phi_{1}||\Phi_{2}^{\dagger}\Phi_{2}|-2|\Phi_{1}^{\dagger}\Phi_{2}|^{2}\\\
&=(w^{2}+v_{2})^{2}\phi_{1}^{+}\phi_{1}^{-}+v_{1}^{2}\phi_{2}^{+}\phi_{2}^{-}-w^{2}\phi_{1}^{+}\Phi_{1}^{-}-v_{1}v_{2}(\phi_{1}^{+}\phi_{2}^{-}+\phi_{2}^{+}\phi_{2}^{-})\\\
&=\begin{pmatrix}\phi_{1}^{-},&\phi_{2}^{-}\end{pmatrix}\begin{pmatrix}v_{2}^{2}&-v_{1}v_{2}\\\
-v_{1}v_{2}&v_{1}^{2}\end{pmatrix}\begin{pmatrix}\phi_{1}^{+}\\\
\phi_{2}^{+}\end{pmatrix}\,.\end{split}$ (94)

Similarly, in expanding about the minimum121212Note that
$\overline{V}_{*}=V_{*}$, since by design the constraint vanishes on the
solution $c-\frac{1}{2}r_{*}^{\mu}r_{*\mu}=0$. of the potential $V_{*}$, the
quadratic part for the charged Higgs fields is:

$V_{*}\supset\frac{m_{H^{\pm}}^{2}}{v^{2}}\begin{pmatrix}\phi_{1}^{-},&\phi_{2}^{-}\end{pmatrix}\begin{pmatrix}v_{2}^{2}&-v_{1}v_{2}\\\
-v_{1}v_{2}&v_{1}^{2}\end{pmatrix}\begin{pmatrix}\phi_{1}^{+}\\\
\phi_{2}^{+}\end{pmatrix}\,,$ (95)

with $H^{\pm}$ related to $\phi^{\pm}$ as:

$H^{\pm}=\frac{v_{2}}{v}\phi_{1}^{\pm}-\frac{v_{1}}{v}\phi_{2}^{\pm}\,.$ (96)

Thus we immediately see that

$\frac{dV_{*}}{dc}=\frac{M_{H^{\pm}}^{2}}{v^{2}}\,.$ (97)

Since $\zeta_{*}=\frac{dV_{*}}{dc}$, we verify the claim131313Technically the
$\zeta$ of Eq. (89) is $\zeta_{*}$, since it is evaluated in a vacuum field
configuration (recall that the star notation denoted a field configuration
extremizing the potential). In the rest of the paper we drop the star
notation. of Eq. (89). Thus, it will generically be the case that the Lagrange
multiplier accompanying the constraint $r^{\mu}r_{\mu}=0$ is proportional to
the charged Higgs mass squared in the neutral vacuum.

## References

  * [1] M. Tanabashi  “Review of Particle Physics”  In _Phys. Rev. D_ 98.3, 2018, pp. 030001  DOI: 10.1103/PhysRevD.98.030001
  * [2] F. Englert and R. Brout  “Broken Symmetry and the Mass of Gauge Vector Mesons”  In _Phys. Rev. Lett._ 13, 1964, pp. 321–323  DOI: 10.1103/PhysRevLett.13.321
  * [3] Peter W. Higgs  “Broken Symmetries and the Masses of Gauge Bosons”  In _Phys. Rev. Lett._ 13, 1964, pp. 508–509  DOI: 10.1103/PhysRevLett.13.508
  * [4] G.. Guralnik, C.. Hagen and T… Kibble  “Global Conservation Laws and Massless Particles”  In _Phys. Rev. Lett._ 13, 1964, pp. 585–587  DOI: 10.1103/PhysRevLett.13.585
  * [5] G.. Branco et al.  “Theory and phenomenology of two-Higgs-doublet models”  In _Phys. Rept._ 516, 2012, pp. 1–102  DOI: 10.1016/j.physrep.2012.02.002
  * [6] Howard E. Haber and Ralf Hempfling  “The Renormalization group improved Higgs sector of the minimal supersymmetric model”  In _Phys. Rev. D_ 48, 1993, pp. 4280–4309  DOI: 10.1103/PhysRevD.48.4280
  * [7] Apostolos Pilaftsis and Carlos E.. Wagner  “Higgs bosons in the minimal supersymmetric standard model with explicit CP violation”  In _Nucl. Phys. B_ 553, 1999, pp. 3–42  DOI: 10.1016/S0550-3213(99)00261-8
  * [8] Marcela Carena, J.. Espinosa, M. Quiros and C… Wagner  “Analytical expressions for radiatively corrected Higgs masses and couplings in the MSSM”  In _Phys. Lett. B_ 355, 1995, pp. 209–221  DOI: 10.1016/0370-2693(95)00694-G
  * [9] Gabriel Lee and Carlos E.. Wagner  “Higgs bosons in heavy supersymmetry with an intermediate mA”  In _Phys. Rev. D_ 92.7, 2015, pp. 075032  DOI: 10.1103/PhysRevD.92.075032
  * [10] M. Carena et al.  “CP Violation in Heavy MSSM Higgs Scenarios”  In _JHEP_ 02, 2016, pp. 123  DOI: 10.1007/JHEP02(2016)123
  * [11] Henning Bahl and Wolfgang Hollik  “Precise prediction of the MSSM Higgs boson masses for low MA”  In _JHEP_ 07, 2018, pp. 182  DOI: 10.1007/JHEP07(2018)182
  * [12] Henning Bahl, Stefan Liebler and Tim Stefaniak  “MSSM Higgs benchmark scenarios for Run 2 and beyond: the low $\tan\beta$ region”  In _Eur. Phys. J. C_ 79.3, 2019, pp. 279  DOI: 10.1140/epjc/s10052-019-6770-z
  * [13] Nick Murphy and Heidi Rzehak  “Higgs-Boson Masses and Mixings in the MSSM with CP Violation and Heavy SUSY Particles”, 2019  arXiv:1909.00726 [hep-ph]
  * [14] Henning Bahl, Nick Murphy and Heidi Rzehak  “Hybrid calculation of the MSSM Higgs boson masses using the complex THDM as EFT”  In _Eur. Phys. J. C_ 81.2, 2021, pp. 128  DOI: 10.1140/epjc/s10052-021-08939-7
  * [15] Henning Bahl and Ivan Sobolev  “Two-loop matching of renormalizable operators: general considerations and applications”  In _JHEP_ 03, 2021, pp. 286  DOI: 10.1007/JHEP03(2021)286
  * [16] Ian Low, Witold Skiba and David Tucker-Smith  “Little Higgses from an antisymmetric condensate”  In _Phys. Rev. D_ 66, 2002, pp. 072001  DOI: 10.1103/PhysRevD.66.072001
  * [17] K.. Klimenko  “On Necessary and Sufficient Conditions for Some Higgs Potentials to Be Bounded From Below”  In _Theor. Math. Phys._ 62, 1985, pp. 58–65  DOI: 10.1007/BF01034825
  * [18] Shinya Kanemura, Takahiro Kubota and Eiichi Takasugi  “Lee-Quigg-Thacker bounds for Higgs boson masses in a two doublet model”  In _Phys. Lett. B_ 313, 1993, pp. 155–160  DOI: 10.1016/0370-2693(93)91205-2
  * [19] Andrew G. Akeroyd, Abdesslam Arhrib and El-Mokhtar Naimi  “Note on tree level unitarity in the general two Higgs doublet model”  In _Phys. Lett. B_ 490, 2000, pp. 119–124  DOI: 10.1016/S0370-2693(00)00962-X
  * [20] P.. Ferreira, R. Santos and A. Barroso  “Stability of the tree-level vacuum in two Higgs doublet models against charge or CP spontaneous violation” [Erratum: Phys.Lett.B 629, 114–114 (2005)]  In _Phys. Lett. B_ 603, 2004, pp. 219–229  DOI: 10.1016/j.physletb.2004.10.022
  * [21] J. Horejsi and M. Kladiva  “Tree-unitarity bounds for THDM Higgs masses revisited”  In _Eur. Phys. J. C_ 46, 2006, pp. 81–91  DOI: 10.1140/epjc/s2006-02472-3
  * [22] I.. Ginzburg and I.. Ivanov  “Tree-level unitarity constraints in the most general 2HDM”  In _Phys. Rev. D_ 72, 2005, pp. 115010  DOI: 10.1103/PhysRevD.72.115010
  * [23] A. Barroso, P.. Ferreira and R. Santos  “Charge and CP symmetry breaking in two Higgs doublet models”  In _Phys. Lett. B_ 632, 2006, pp. 684–687  DOI: 10.1016/j.physletb.2005.11.031
  * [24] I.. Ivanov  “Minkowski space structure of the Higgs potential in 2HDM” [Erratum: Phys.Rev.D 76, 039902 (2007)]  In _Phys. Rev. D_ 75, 2007, pp. 035001  DOI: 10.1103/PhysRevD.75.035001
  * [25] M. Maniatis, A. Manteuffel, O. Nachtmann and F. Nagel  “Stability and symmetry breaking in the general two-Higgs-doublet model”  In _Eur. Phys. J. C_ 48, 2006, pp. 805–823  DOI: 10.1140/epjc/s10052-006-0016-6
  * [26] A. Barroso, P.. Ferreira and R. Santos  “Neutral minima in two-Higgs doublet models”  In _Phys. Lett. B_ 652, 2007, pp. 181–193  DOI: 10.1016/j.physletb.2007.07.010
  * [27] Igor P. Ivanov  “Minkowski space structure of the Higgs potential in 2HDM. II. Minima, symmetries, and topology”  In _Phys. Rev. D_ 77, 2008, pp. 015017  DOI: 10.1103/PhysRevD.77.015017
  * [28] P.. Ferreira and D… Jones  “Bounds on scalar masses in two Higgs doublet models”  In _JHEP_ 08, 2009, pp. 069  DOI: 10.1088/1126-6708/2009/08/069
  * [29] D. Jurčiukonis and L. Lavoura  “The three- and four-Higgs couplings in the general two-Higgs-doublet model”  In _JHEP_ 12, 2018, pp. 004  DOI: 10.1007/JHEP12(2018)004
  * [30] Yisheng Song  “Co-positivity of tensors and Stability conditions of two Higgs potentials”, 2022  arXiv:2203.11462 [hep-ph]
  * [31] Antonio Pich and Paula Tuzon  “Yukawa Alignment in the Two-Higgs-Doublet Model”  In _Phys. Rev. D_ 80, 2009, pp. 091702  DOI: 10.1103/PhysRevD.80.091702
  * [32] Nina M. Coyle, Bing Li and Carlos E.. Wagner  “Wrong sign bottom Yukawa coupling in low energy supersymmetry”  In _Phys. Rev. D_ 97.11, 2018, pp. 115028  DOI: 10.1103/PhysRevD.97.115028
  * [33] Albert M Sirunyan  “Combined measurements of Higgs boson couplings in proton–proton collisions at $\sqrt{s}=13\,\text{Te}\text{V}$”  In _Eur. Phys. J. C_ 79.5, 2019, pp. 421  DOI: 10.1140/epjc/s10052-019-6909-y
  * [34] Georges Aad  “Combined measurements of Higgs boson production and decay using up to $80$ fb-1 of proton-proton collision data at $\sqrt{s}=$ 13 TeV collected with the ATLAS experiment”  In _Phys. Rev. D_ 101.1, 2020, pp. 012002  DOI: 10.1103/PhysRevD.101.012002
  * [35] John F. Gunion and Howard E. Haber  “The CP conserving two Higgs doublet model: The Approach to the decoupling limit”  In _Phys. Rev. D_ 67, 2003, pp. 075019  DOI: 10.1103/PhysRevD.67.075019
  * [36] Marcela Carena, Ian Low, Nausheen R. Shah and Carlos E.. Wagner  “Impersonating the Standard Model Higgs Boson: Alignment without Decoupling”  In _JHEP_ 04, 2014, pp. 015  DOI: 10.1007/JHEP04(2014)015
  * [37] P.. Bhupal Dev and Apostolos Pilaftsis  “Maximally Symmetric Two Higgs Doublet Model with Natural Standard Model Alignment” [Erratum: JHEP 11, 147 (2015)]  In _JHEP_ 12, 2014, pp. 024  DOI: 10.1007/JHEP12(2014)024
  * [38] Marcela Carena et al.  “Complementarity between Nonstandard Higgs Boson Searches and Precision Higgs Boson Measurements in the MSSM”  In _Phys. Rev. D_ 91.3, 2015, pp. 035003  DOI: 10.1103/PhysRevD.91.035003
  * [39] Marcela Carena et al.  “Alignment limit of the NMSSM Higgs sector”  In _Phys. Rev. D_ 93.3, 2016, pp. 035013  DOI: 10.1103/PhysRevD.93.035013
  * [40] Jérémy Bernon et al.  “Scrutinizing the alignment limit in two-Higgs-doublet models: mh=125 GeV”  In _Phys. Rev. D_ 92.7, 2015, pp. 075004  DOI: 10.1103/PhysRevD.92.075004
  * [41] Howard E. Haber, Venus Keus and Rui Santos  “P-even, CP-violating Signals in Scalar-Mediated Processes”, 2022  arXiv:2206.09643 [hep-ph]
  * [42] Kenneth Ross Garren  “Bounds for the Eigenvalues of a Matrix”  In _Dissertations, Theses, and Masters Projects_ , 1965  DOI: 10.21220/s2-p2tn-7046
  * [43] S. Gershgorin  “Über die Abgrenzung der Eigenwerte einer Matrix”  In _Bulletin de l’Academie des Sciences de l’URSS_ 6.3, 1931, pp. 749–754
  * [44] Miguel P. Bento, Jorge C. Romão and João P. Silva  “Unitarity bounds for all symmetry-constrained 3HDMs”, 2022  arXiv:2204.13130 [hep-ph]
  * [45] Shinya Kanemura and Kei Yagyu  “Unitarity bound in the most general two Higgs doublet model”  In _Phys. Lett. B_ 751, 2015, pp. 289–296  DOI: 10.1016/j.physletb.2015.10.047
  * [46] G.C. Branco et al.  “Theory and phenomenology of two-Higgs-doublet models”  In _Physics Reports_ 516, 2011, pp. 1–102  DOI: 10.1016/j.physrep.2012.02.002
  * [47] Gary Ulrich and Layne T. Watson  “Positivity Conditions for Quartic Polynomials”  In _SIAM Journal on Scientific Computing_ 15.3, 1994, pp. 528–544  DOI: 10.1137/0915035
  * [48] A. Barroso, P.. Ferreira, I.. Ivanov and Rui Santos  “Metastability bounds on the two Higgs doublet model”  In _JHEP_ 06, 2013, pp. 045  DOI: 10.1007/JHEP06(2013)045
  * [49] I.. Ivanov and Joao P. Silva  “Tree-level metastability bounds for the most general two Higgs doublet model”  In _Phys. Rev. D_ 92.5, 2015, pp. 055017  DOI: 10.1103/PhysRevD.92.055017
  * [50] L. Lavoura and Joao P. Silva  “Fundamental CP violating quantities in a SU(2) x U(1) model with many Higgs doublets”  In _Phys. Rev. D_ 50, 1994, pp. 4619–4624  DOI: 10.1103/PhysRevD.50.4619
  * [51] Sacha Davidson and Howard E. Haber  “Basis-independent methods for the two-Higgs-doublet model” [Erratum: Phys.Rev.D 72, 099902 (2005)]  In _Phys. Rev. D_ 72, 2005, pp. 035004  DOI: 10.1103/PhysRevD.72.099902
  * [52] Ian Low, Nausheen R. Shah and Xiao-Ping Wang  “Higgs alignment and novel CP-violating observables in two-Higgs-doublet models”  In _Phys. Rev. D_ 105.3, 2022, pp. 035009  DOI: 10.1103/PhysRevD.105.035009
  * [53] Howard E. Haber and Oscar Stål  “New LHC benchmarks for the $\mathcal{CP}$ -conserving two-Higgs-doublet model” [Erratum: Eur.Phys.J.C 76, 312 (2016)]  In _Eur. Phys. J. C_ 75.10, 2015, pp. 491  DOI: 10.1140/epjc/s10052-015-3697-x
  * [54] Rafael Boto et al.  “Basis-independent treatment of the complex 2HDM”  In _Phys. Rev. D_ 101.5, 2020, pp. 055023  DOI: 10.1103/PhysRevD.101.055023
  * [55] V. Andreev  “Improved limit on the electric dipole moment of the electron”  In _Nature_ 562.7727, 2018, pp. 355–360  DOI: 10.1038/s41586-018-0599-8
  * [56] Wolfgang Altmannshofer, Stefania Gori, Nick Hamer and Hiren H. Patel  “Electron EDM in the complex two-Higgs doublet model”  In _Phys. Rev. D_ 102.11, 2020, pp. 115042  DOI: 10.1103/PhysRevD.102.115042
  * [57] Sacha Davidson and Howard E. Haber  “Basis-independent methods for the two-Higgs-doublet model”  In _Phys. Rev. D_ 72.3, 2005, pp. 035004  DOI: 10.1103/PhysRevD.72.035004