ℤ2 breaking effects in 2-loop RG evolution of 2HDM

We investigate the effects of a ℤ2 symmetry in the CP\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{C}\mathcal{P} $$\end{document}-conserving Two-Higgs-Doublet-Model (2HDM); which is often imposed to prevent Flavor-Changing-Neutral-Currents (FCNCs) at tree-level. Specifically, we analyze how a breaking of the ℤ2 symmetry spreads during renormalization group evolution; employing general 2-loop renormalization group equations that we have derived. Evolving the model from the electroweak to the Planck scale, we find that while the case of an exact ℤ2 symmetric 2HDM is very constrained, a soft breaking of the ℤ2 symmetry extends the valid parameter space regions. The effects of a hard ℤ2 breaking in the scalar sector as well as the stability of the flavor alignment ansatz are also investigated. We find that while a hard breaking of the ℤ2 symmetry in the potential is problematic, since it speeds up the growth of quartic couplings, the generated FCNCs are heavily suppressed. Conversely, we also find that hard ℤ2 breaking in the Yukawa sector at most gives moderate ℤ2 breaking in the potential; whereas the FCNCs can become quite sizable far away from the ℤ2 symmetric regions.


Introduction
Since the discovery of a 125 GeV scalar particle at the Large Hadron Collider (LHC), by the ATLAS [1] and CMS [2] collaborations, the quest to decipher its true nature has begun. So far, it closely resembles the Standard Model (SM) Higgs boson [3]; making the minimal SM electroweak sector a viable description. However, a thorough experimental investigation is required to determine the exact dynamics of electroweak symmetry breaking.
On the theoretical side, there are still arguments for a non-minimal scalar sector. One of the simplest extensions of the SM is to add another Higgs doublet. There are numerous motivations for doing so, e.g. SUSY models need another Higgs doublet to cancel anomalies; JHEP02(2019)152 realistic models of baryogenesis cannot be constructed out of the SM; and the vacuum metastability of the SM can be rendered stable. It is also interesting in its own right to investigate the effects of having an extended scalar sector, much like there are three generations of fermions in the SM.
Two-Higgs-doublet models (2HDMs) exhibit a rich phenomenology with three new scalar particles and potential sources of CP violation, to name only a few of its features. It has been studied for a long time and we refer to ref. [4] for a recent review.
An immediate problem when introducing an additional Higgs doublet to the SM is that the Yukawa sector, in its general form, gives rise to flavor changing neutral currents (FCNCs) at tree-level, that are severely constrained by experiments. These arise since it is in general not possible to diagonalize the Yukawa matrices for both Higgs doublets simultaneously, when going to the fermion mass eigenbasis. One mechanism, often employed in the literature, that solves this problem is to impose a discrete Z 2 symmetry on the 2HDM as proposed in refs. [5,6].
In this paper, we are interested in the effects of breaking such a Z 2 symmetry. More specifically, we investigate how stable an approximate Z 2 symmetry of the 2HDM is in cases where it is explicitly broken by a small amount. To investigate these effects, we employ a renormalization group (RG) analysis of the parameters of the model. A small breaking of the Z 2 symmetry at one energy scale will, in general, spread during the RG evolution and generate additional Z 2 breaking parameters. Our main aim of this work is to give a quantitative estimate of how large these Z 2 breaking parameters can be when requiring a model that is stable, unitary and perturbative up to some energy scale, where potentially some new physics will be present. Experimental constraints are also taken into consideration by making sure that the model is within the bounds that are implemented in HiggsBounds [7][8][9] and HiggsSignals [10]. We look both at the cases with a softly broken as well as a hard broken Z 2 scalar potential. The softly broken Z 2 scenario is very common in the literature and is also stable; in the sense that no dimensionless Z 2 symmetry breaking parameters are generated during the RG running. A hard breaking in the scalar potential, on the other hand, can have severe effects on the RG flow of the quartic couplings and potentially induce problematically large FCNCs in the Yukawa sector.
As a scenario of Z 2 breaking in the Yukawa sector, we look at the ansatz of flavor alignment in the 2HDM [11,12]. Although the flavor alignment solves the problem of FCNCs at a specific energy scale, the ansatz is not stable during RG evolution [13][14][15]; which also induces Z 2 breaking parameters in the scalar potential.
At 2-loop order, the RG analysis becomes more complicated since there is a stronger coupling between the scalar and Yukawa sectors; the quartic couplings enter the Yukawa beta functions first at 2-loop order. This also makes it more interesting when looking at how a small breaking of the Z 2 symmetry can spread during RG evolution. In refs. [22][23][24], the authors employ 2-loop RGEs in their analysis, but on a CP-conserving 2HDM, with a JHEP02(2019)152 softly broken Z 2 symmetry. We have derived the full set of 2-loop RGEs of the general, potentially complex, 2HDM and implemented them in an open-source C++ program called 2HDME [26].
The structure of this paper is as follows: in section 2, we review the needed theory of the 2HDM. We discuss the technical details of the RG evolution in section 3; how the RGEs are derived and the algorithm to perform the evolution. The RGEs are used in parameter scans that we describe and present the results of in section 4. Finally, we present our conclusions in section 5.

The 2HDM
Here, we will briefly review the content of the most general renormalizable 2HDM, i.e. the standard model, SU(3) c × SU(2) W × U(1) Y , with an additional Higgs doublet. The 2HDM has been studied extensively and for a thorough review see ref. [4].
The 2HDM contains two hypercharge +1 complex scalar SU(2) doublets, Φ 1 and Φ 2 . First of all, since the scalar fields have identical quantum numbers, one can always perform a field redefinition of the scalar fields, i.e. a non-singular complex transformation Φ a → B ab Φ b , 1 where the matrix B depends on 8 real parameters. One uses four of these parameters to transform to canonical diagonal kinetic terms. The Lagrangian of the 2HDM exhibits a U(2) Higgs flavor symmetry, Φ a → U ab Φ b ; since the Lagrangian keeps the same form after such a transformation. We will denote 2HDMs related by such Higgs flavor transformations as different bases of the 2HDM. It is therefore very important, when investigating the general 2HDM, to work with basis-independent quantities. A thorough basis-independent treatment of 2HDM has been developed in refs. [27][28][29] and we will, mostly, follow their notational conventions.

Generic basis
The most general 2HDM gauge invariant renormalizable scalar potential can be written where m 2 12 , λ 5 , λ 6 and λ 7 are potentially complex while all the other parameters are real; resulting in a total of 14 degrees of freedom. Three of these will be removed by the tadpole equations and one can be removed by a re-phasing of the second Higgs doublet. Thus, the most general 2HDM potential exhibits 10 physical degrees of freedom; or 11 if one includes the SM Higgs VEV, 2 v = 2m W /g 2 ≈ 246 GeV. 1 The bar notation keep tracks of complex conjugation. That is, replacing a barred index to an unbarred corresponds to complex conjugation [27][28][29]. 2 We denote the SU(3)c × SU(2)W × U(1)Y gauge couplings as g3, g2, g1 respectively.

JHEP02(2019)152
After electroweak symmetry breaking, SU(2) × U(1) Y → U(1) em , both the scalar fields acquire a VEV, which can be expressed in terms of a unit vector in the Higgs flavor space where the unit vector is normalized tov * av a = 1. By convention, we take 0 ≤ β ≤ π/2 and 0 ≤ ξ ≤ 2π. Here, we have used up all our gauge freedom, when setting the VEV in the lower component of the doublets with a SU(2) transformation and removing any phase in the Φ 1 VEV with a U(1) Y transformation. We also defineŵ b ≡v * a ab , where 12 = − 21 = 1 and 11 = 22 = 0.
The angle β can be defined in terms of the ratio of the Higgs fields, It should be noted that tan β is an unphysical parameter if the two Higgs fields have identical quantum numbers [28]; since there is no preferred basis in that case. In subsequent sections it will however be promoted to a physical quantity when a Z 2 symmetry is imposed on the 2HDM in a particular basis. The physical scalar degrees of freedom, after electroweak symmetry breaking, corresponds to three neutral and one U(1) em charged pair of Higgs bosons, {h 1 , h 2 , h 3 , H ± }. In the CP-conserving case, the neutral mass eigenstates have definite CP properties. Instead of h 1,2,3 , we can then work with the two CP even states, denoted h and H, 3 as well as a CP odd one, A.

Higgs basis
One particularly convenient basis is the Higgs basis [27,30], where only one Higgs field gets a VEV. The Higgs basis fields in terms of the previously defined generic basis fields are which acquire the VEVs The transformation between the bases is The scalar potential in the Higgs basis takes a similar form as in the generic basis,

JHEP02(2019)152
Minimizing this potential leads to the tree-level tadpole equations During a Higgs flavor transformation of the generic basis, Φ a → U ab Φ b , the Higgs fields transform as [28] H Thus from inspection of the Higgs potential in eq. (2.7), it follows that Y 1 , are pseudo-invariants under the Higgs flavor transformation. The Higgs basis is therefore unique up to a rephasing of H 2 , arising from the freedom to perform a Higgs flavor transformation.

Yukawa sector
The Yukawa interactions in the generic basis are 4 where the left-handed fermion fields in the weak eigenbasis are andΦ ≡ iσ 2 Φ * . To go to the fermion mass eigenbasis, we perform the biunitary transformation where F ∈ {U, D, E} is each fermion species. The CKM matrix is composed out of the left-handed transformation matrices, V CKM ≡ V U L V D † L , and the Yukawa matrices in the mass eigenbasis are obtained with biunitary transformations (2.14) Going to the Higgs basis, we get where κ F are the diagonal mass matrices,

JHEP02(2019)152
and ρ F =ŵ * a η F a are arbitrary complex 3 × 3 matrices. If the ρ F matrices are non-diagonal, there are FCNCs present at tree-level. There exist multiple solutions to the problem of FCNCs in the 2HDM. One of them is the idea of alignment [11], which makes the ansatz that η F,0 1 and η F,0 2 are proportional to each other; with the consequence that κ F and ρ F can be diagonalized simultaneously. We will parameterize the alignment ansatz by setting where the alignment parameter a F could be a c-number, but we will restrict our analysis to real coefficients. As is well known, the alignment ansatz is not stable under RG evolution [13][14][15], since there is no symmetry protecting it. Thus, eq. (2.17) is only valid at a specific energy scale and during RG evolution FCNCs will be generated.
One ansatz, that allows for small FCNCs is the Cheng-Sher ansatz [31] which parameterizes the non-diagonal Yukawa matrices as where the λ F ij should be order one. We will use this parameterization when comparing generated FCNCs later on.
One can get limits on the non-diagonal elements from neutral meson mixing at low energy scales, see [16] and references therein. Even though we will discuss FCNCs that are generated at high energy scales, we will use the limit λ F i =j ≤ 0.1 as a general measure of sizable non-diagonal elements. Z 2 symmetry. Another solution to the FCNC problem is the idea of imposing a Z 2 symmetry [5,6], where the Higgs doublets have opposite charge. By only coupling each right-handed fermion to one of the Higgs doublets, one can diagonalize all Yukawa matrices and end up with an alignment like in eq. (2.17), but with the alignment parameters set to be equal to either − tan β or cot β. The Z 2 symmetric Yukawa sector is thus a special case of an aligned Yukawa sector; however, because of the symmetry, a Z 2 symmetric 2HDM is stable under RG evolution, in that no FCNCs are generated.
There are four different variations of this symmetry that are summarized in table 1. It is also worth noting that the previously unphysical parameter tan β = v 2 /v 1 gets promoted to a physical degree of freedom, when imposing a Z 2 symmetry, see ref. [28] for details.
Requiring the generic potential in eq. (2.1) to be symmetric under a Z 2 symmetry fixes m 2 12 = λ 6 = λ 7 = 0. Often in the literature, one relaxes the Z 2 symmetry by letting m 2 12 be non-zero; thus breaking the symmetry softly. It should be noted that the mass parameters of the scalar potential do not enter the beta functions of the quartic or Yukawa couplings, thus a soft symmetry breaking does not alter the RG evolution of the dimensionless parameters.
One can also consider hard breaking of the Z 2 symmetry in the scalar potential, by having small λ 6,7 . This introduces additional terms in the RGEs and also induces FCNCs JHEP02(2019)152 Table 1. Different Z 2 symmetries that can be imposed on the 2HDM. Φ 1 is odd(−1) and Φ 2 is even(+1). For every type of Z 2 symmetry, the ρ F matrices become proportional to the diagonal mass matrices, ρ F = a F κ F . at different energy scales after 2-loop RGE running. It is one of this article's main goal to investigate how severe these effects can be; see section 4.3 for more details.

CP-invariant limit
The scalar potential and vacuum are CP-conserving if and only if [27,[32][33][34] Im Z * 5 Z 2 6 = Im Z * 5 Z 2 7 = Im (Z * 6 Z 7 ) = 0. (2.19) If these conditions are satisfied, it is possible to perform a phase shift of H 2 to end up with a Higgs basis composed of purely real parameters making it a real basis. An exact Z 2 symmetric potential only has one potentially complex parameter, λ 5 , that can be made real with a phase shift. Also the phase of Φ 2 can be removed with a field redefinition; thus one ends up with a 2HDM with all real VEVs and couplings. Thus an exact Z 2 symmetry enforces CP-invariance in the scalar sector. As previously mentioned, one can relax the imposed Z 2 symmetry constraint by allowing for a softly broken symmetry from a non-zero m 2 12 parameter in the generic basis. However, doing so does no longer guarantee CP-invariance. To make sure that the model is still CP-invariant, one can restrict the parameters such that [35] Im λ * 5 (m 2 12 ) 2 = 0, and (2.20) which will guarantee the absence of both explicit and spontaneous CP-violation in the scalar sector. It would be interesting to look at how CP-violation affects the RGE running of the parameters, but this is beyond the scope of this work. We will therefore assume CP-invariance in the scalar sector; making it possible to always work in a real basis. 5

Hybrid basis
When imposing a softly broken Z 2 symmetry on the CP-conserving 2HDM, the number of quartic couplings is reduced to 5 real ones. It is often very convenient to work with a basis where some quartic couplings are substituted for tree-level mass parameters, for physical  Table 2. Hybrid basis for a softly broken Z 2 symmetric CP-invariant 2HDM from ref. [35]. Imposing an exact Z 2 symmetry fixes Z 7 according to eq. (2.28), withm 2 = 0.
clarity. However, when doing numerical parameter scans, it can be computationally expensive to work with bases with many masses. A bad choice of scalar masses easily corresponds to large quartic couplings that break perturbativity and unitarity. It is therefore a good idea to set up a basis with a mixture of scalar masses, angles and quartic couplings. Such a hybrid basis is worked out in ref. [35] and is summarized in table 2, which we will make use of in upcoming sections. The hybrid basis contains the two tree-level masses of the neutral CP even scalars, )v 2 is the mass of the CP-odd scalar and m h ≤ m H . In addition to β, the hybrid basis contains an angle α and we follow the convention that 0 ≤ β − α ≤ π. Here, we simply give the relations [35] which are valid for the case of Z 6 = 0. The three quartic couplings Z 4,5,7 are free parameters in the hybrid basis while There is also the relation For an exact Z 2 symmetry,m 2 = 0 gives that Z 7 is not a free parameter in this basis. Finally, the mass of the charged Higgs boson is given by

Renormalization group evolution
The behavior of the parameters of the 2HDM during RG evolution is the focus of attention in this paper. Even though there are many parameters in the general 2HDM that are essentially free, requiring "good" behavior when evolving the model in energy restricts the parameter space severely. An RGE analysis is thus useful when looking for potential finetuning, stability of the model and valid energy ranges before new physics would need to come in. A careful matching to physical observables should be done to get the most precise determination of the 2HDM's MS-parameters. Including higher order corrections to observables can then be important, since it may shift the parameters by a non-trivial amount. It is especially important when considering RG evolution of the model. A common procedure in the literature is to do the matching at one loop order lower than the loop order of the RGEs; there are however indications of that N-loop running should be combined with N-loop matching [23]. In this work, we are more interested in the general behavior of the MS-parameters during RG evolution and not so much in the exact determination of the physical observables. We will thus restrict ourselves to calculate 1-loop pole masses of the scalar particles as the only quantum corrections. These corrections can indeed be very large in certain parts of the 2HDMs parameter space, as will be shown in subsequent sections.
We will impose theoretical constraints, namely perturbativity, unitarity and stability, that are needed to arrive at a consistent model. To check against experimental data from collider searches, we use HiggsBounds [7][8][9] and HiggsSignals [10].
Perturbativity. Since our analysis of the RG running of the 2HDM uses the RGEs, calculated with perturbation theory, perturbativity must always be satisfied during the RG evolution; breaking it would render the whole analysis meaningless. The couplings can JHEP02(2019)152 therefore not be too large and we will impose the upper limit to specify perturbativity. The RGEs make up a strongly coupled set of ordinary differential equations, which implies that once any coupling gets large others soon follow; unless there is some kind of cancellation occurring among the parameters. Breaking of perturbativity at some energy scale can therefore be seen as an indication of a Landau pole at a scale not far from it.
Unitarity. Unitarity of the model is very important for consistency, since it is needed for a well defined S-matrix, and enforcing it puts constraints on the eigenvalues of the S-matrix for scalar scattering. At high energies, only the quartic couplings contribute to the amplitudes; they are therefore the ones affected by the unitarity limit. In ref. [36], the conditions for tree-level unitarity in a general 2HDM has been worked out. They take the form of upper limits on eigenvalues of scattering matrices, composed of the quartic couplings. For completeness, we give these matrices in appendix B.
Bounded potential. To ensure that the vacuum is a stable minimum of the potential, we require that the potential is bounded from below. This means that the potential should be positive when the field values go to infinity for any direction in field space. The sufficient conditions for a general renormalizable 2HDM have been worked out in refs. [37,38], which we will make use of. Enforcing stability constrains the quartic couplings in the scalar potential and, for completeness, we state the conditions in appendix C.
In this work, we are only enforcing tree-level stability conditions on the potential. Sometimes, this is a too strong constraint since loop corrections can render the potential stable. A more formal test of stability would be to check the stability of the effective potential. Thus, if one instead were to check the stability of the 1-loop effective potential, it would relax the stability constraint as in ref. [24]; however, this is beyond the scope of this work. In some cases, we will omit the stability constraint when investigating the maximal effects of symmetry breaking; in order not to miss any region of parameter space that might exhibit interesting behavior.
Experimental bounds. Constraints coming from LEP, the Tevatron and the LHC are taken into consideration using HiggsBounds [7][8][9]; which uses experimental cross section limits to determine if a certain parameter point is excluded at a 95% C.L. . To ensure that the lightest Higgs boson resembles the 125 GeV Higgs signal observed at the LHC, we use HiggsSignals [10]. It performs a statistical goodness-of-fit test by calculating the χ 2 for each parameter point. The χ 2 is used to calculate a p-value with the number of degrees of freedom set to the number of free parameters in the scalar potential; we exclude points with a p-value less than 0.05. We also use 2HDMC [39] to calculate the decay widths and cross sections needed for HiggsBounds and HiggsSignals.

Derivation of 2-loop RGEs
Although there are some works using 2-loop RGEs for analyzing 2HDMs [22][23][24], they all work with different Z 2 symmetries and the RGEs used are therefore not applicable to our more general scenario. 6 The general expressions for the 2-loop RGEs of massless parameters in any renormalizable gauge theory were first written down in the seminal papers of Machacek and Vaughn [44][45][46]. These were then supplemented with the 2-loop RGEs of massive parameters in ref. [47], which is the source that we have used to derive the 1-and 2-loop RGEs for a general 2HDM. Care should be taken though, when working with quantum field theories with multiple indistinguishable scalar fields, since the formulas in refs. [44][45][46][47] are written for the case of an irreducible representation of the scalar fields. The anomalous dimension of scalar fields with equal quantum numbers is non-diagonal; the fields mix during RG evolution. In 2HDMs, the fields obey the RGEs where γ ij is the 2-dimensional anomalous dimension. Refs. [25,43] discuss this subtlety at great length and we have independently reached the same conclusions.
As cross checks, we compared the SM limit with known SM RGEs as well as to Z 2 symmetric RGEs found in the literature.
In ref. [48] it is argued that one needs to consider scalar kinetic mixing in RG analyses of quantum field theories with multiple scalar fields with equal quantum numbers. We claim that this is not needed, since there is no physical parameter related to scalar kinetic mixing. The issue can be resolved differently depending on the renormalization schemes; which we explain in more detail in ref. [49]. Here, the scalar mixing phenomenon only manifests itself in the previously mentioned 2-dimensional non-diagonal anomalous dimension of the Higgs fields. The anomalous dimension sets the evolution of the Higgs VEVs in our scheme; with the implication that tan β runs in energy. Thus the Higgs flavor transformation matrix in eq. (2.6) is µ-dependent.
In the general case of no Z 2 symmetry, the equations are very long and we will not write them down in this article. They are instead available in the C++ code 2HDME [26].

Technical details of numerical code
Evolving the 2HDM in energy corresponds to solving a system of 129 7 coupled ordinary differential equations. We have developed the C++ code 2HDME [26] to perform this task. It uses GSL [50] to solve the RGEs and the library Eigen [51] for linear algebra operations. 6 Refs. [23,24] are using 2-loop RGEs derived by SARAH [40,41] and ref. [22] is using PyR@TE [42]. As have been pointed out in refs. [25,43], older versions of these software miss non-diagonal anomalous dimension terms in the case of models containing scalar fields with equal quantum numbers. The issue is resolved in ref. [43] and newer versions of SARAH and PyR@TE should give the correct RGEs. 7 3 gauge couplings; 2 complex VEVs; 6 complex Yukawa matrices; 2 real and 1 complex mass parameters; 4 real and 3 complex quartic couplings.

JHEP02(2019)152
To match the 2HDM to the µ dependent parameters of the SM at the top quark pole mass scale, M t = 173.4 GeV [52], we used the following input as boundary conditions: 8 • The MS fermion masses are used to fix the Yukawa matrix elements in the fermion mass eigenbasis. We use the ones from ref. [53]: The numerical values λ = 0.22453, are extracted from the PDG [55].

JHEP02(2019)152
To calculate the loop corrected pole masses of the Higgs bosons, we make use of the Fortran code SPheno [56,57], together with model files generated by SARAH [40,41]. We choose to work only at 1-loop order when calculating the loop corrected masses for computational efficiency. This choice incorporates the largest loop corrections, but also treats all scalars equally.

Evolution algorithm
Here, we briefly describe the algorithm used when scanning the parameter space of the 2HDM: • Initialize the 2HDM with the SM values at µ = M t . The scalar potential is randomly generated in the hybrid basis; however, only potentials with a 1-loop corrected mass of 125±5 GeV for one of the scalar particles are accepted. Note that although we use the hybrid basis as a first input, we transform to and evolve the parameters in the generic basis.
The Yukawa sector is set in the fermion mass eigenbasis. First, κ F is set with the fermion masses in eq. (3.3) as input. The ρ F are then fixed to be proportional to the κ F matrices, either by a Z 2 symmetry or flavor ansatz. We then go to the fermion weak eigenbasis in the scalar generic basis with a =v a κ L +ŵ a ρ L . • The parameters in the generic basis, i.e. the set {g i , v a , η F,0 a , m 2 ij , λ i }, are evolved according to the RGEs. This includes the VEVs, v a ≡v a v, which in our scheme obey the same RGE as the scalar fields, i.e.
where γ ab is defined in the Landau gauge. A consequence of this is that tan β is µ-dependent.
• As previously mentioned, the fact that tan β is µ-dependent implies a µ-dependent transformation of the generic basis to the Higgs basis with the matrix in eq. (2.6). The parameters in the Higgs basis, {Y i , Z i , κ F,0 , ρ F,0 }, are at each step calculated with thisÛ ab (µ). Also, κ F , ρ F and V CKM are calculated with biunitary transformations.
• The RG evolution is stopped when perturbativity is broken or when the (reduced) Planck scale is reached, which we take to be 10 18 GeV. We define the breakdown scale, Λ, as the lowest energy scale where either perturbativity, unitarity or tree-level

JHEP02(2019)152
stability is violated. Note that this means that we continue the evolution even if the unitarity or stability constraints are broken. In this way, we can study each constraint individually and also loosen the stability constraint since this is applied at tree-level.

RG evolution example
As a first example of the RG evolution of a 2HDM and to see the qualitative behavior of the parameters, we here evolve the type-I softly broken Z 2 parameter point 9 tan β = 2.49772,

Parameter space scan
To probe the RG behavior of the 2HDM's parameter space, we construct random scans of the free parameters in scenarios with different levels of Z 2 symmetry. Although CP violation would be very interesting to investigate, we restrict ourselves to CP-invariant scalar potentials for simplicity, where a real basis always exists. The only source of CP-violation is then the δ phase in the CKM matrix, defined in eq. (3.5); however, its effects turn out to be negligible. As explained in section 3.3, we employ the hybrid basis reviewed in section 2.5 as a starting point in the generation of random parameter points. The number of free parameters, however, depends on the nature of the Z 2 symmetry. We will consider three different levels of Z 2 symmetry in the scalar potential: • Scenario I: 2HDM with an exact Z 2 symmetry. 10 This fixes m 2 12 = λ 6 = λ 7 = 0; which makes Z 7 no longer a free parameter since it is determined by eq. (2.28).
• Scenario II: softly broken Z 2 symmetry, with non-zero m 2 12 . In the hybrid basis, this corresponds to having Z 7 as a free parameter.
• Scenario III: hard broken Z 2 symmetry in scalar sector that allows for small real values of λ 6 and λ 7 . In the hybrid basis this means that Z 2 and Z 3 will deviate from eqs. (2.25) and (2.26).
There are 6, 7 and 9 degrees of freedom in each scenario respectively.
The tree-level mass can deviate significantly from the 1-loop corrected mass and we therefore start from a tree-level mass m tree h ∈ [10, 130] GeV. For the heavier CP even Higgs boson, we scan a tree-level mass m H ∈ [150, 1000] GeV. To have m h to be SM-like, c β−α should be close to zero; however, we use the generous range c β−α ∈ [−0.5, 0.5]. We also take tan β ∈ [1.1, 50]. 11 The quartic couplings Z 4,5,7 are in the range -π to π. 10 An exact Z2 symmetric 2HDM also implies CP-invariance. 11 Although we actually generate random β angles in the corresponding range from a flat distribution.
The parameter region tan β = 1 is excluded for convenience, since it is singular in our parametrization of the potential; however, the singular behavior can be removed by re-parametrization.

JHEP02(2019)152
In the hard Z 2 breaking scenario III, we add deviations from the softly broken scenario in the generic basis and scan over the range λ 6,7 ∈ 10 −5 , 2 , although logarithmically distributed.
We have performed numerical scans of these three Z 2 scenarios; producing 10 5 random parameter points for each one, that are perturbative, unitary and stable at the electroweak scale with a 1-loop corrected m h = 125 ± 5 GeV. To take into account constraints from collider searches, the events are also filtered by only keeping points that are allowed by HiggsBounds and HiggsSignals. As mentioned in section 3.3, these models are then evolved until perturbativity is lost or the Planck scale at 10 18 GeV has been reached and we define the breakdown energy of each parameter point, Λ, as the energy scale where either perturbativity, unitarity or stability is violated.
The specific type-choice of Yukawa symmetry has been checked not to matter for our conclusions about Z 2 symmetry breaking effects; the scans of scenario I-III are performed with a type I symmetry.
Lastly, we also investigate the effects of breaking the Z 2 in the Yukawa sector: • Scenario IV.a: starting from the softly broken Z 2 2HDM defined in eq. (3.10), we make an alignment ansatz ρ F = a F κ F and vary each a F ∈ R separately, i.e. one-by-one in the ranges a U ∈ [0.01, 10] and a D ∈ [−100, 100]; for larger values the corresponding Yukawa couplings become too large for perturbation theory to apply.
• Scenario IV.b: we make a similar softly broken Z 2 parameter scan as in scenario II; however, we fix tan β = 2 and draw random a U ∈ [−1, 1] and a D ∈ [−50, 50]; which are limited to the most relevant ranges. For this scenario, we only consider models where neither perturbativity, unitarity nor stability is violated and we want to compare the generated Z 2 breaking parameters at a common scale. The scale choice is motivated by being able to find parameter points that are valid throughout the selected a F ranges. A higher scale would make the parameter scan more costly in terms of computer time, but we have checked that the dependence on the scale used is very small. In the end, we choose the scale to be µ = 10 10 GeV.
We impose the same theoretical and experimental constraints at the electroweak scale as in scenario I-III. All of these scenarios correspond to a bottom-up running of the parameters. Hence, we are only probing how large Z 2 breaking parameters can be at the EW scale and how fast an alignment ansatz breaks down. Formally, there are of course no limits on the Z 2 breaking parameter at high energy scales; which have not been probed experimentally. However, following ref. [16] we assume that their values should not differ by a large amount compared to the ones at the EW scale. A large sensitivity, such that the Z 2 breaking parameters change with many orders of magnitude in the RGE evolution, would indicate that the underlying assumptions of the model are either fine-tuned or not stable. It would also be interesting to look at top-down scenarios, where one starts of with small deviations from a Z 2 symmetric parameter point at a high scale and then see how large the Z 2 breaking effects are at the EW scale. Such a scenario with a more symmetrical theory at higher scales would  Table 3. How many parameter points that pass the constraints from HiggsBounds and HiggsSignals in scenario I to III.  perhaps be deemed a more natural one, but we will leave such an investigation to future works. However, our analysis provides an indication of how large such effects could be. Before we investigate each respective scenario in more detail, we study at which scales the individual constraints are violated and which one gives rise to the breakdown energy Λ. First of all, the statistics of how many parameter points that pass the constraints coming from HiggsBounds and HiggsSignals is shown in table 3 for scenario I to III. For the points that are within these limits, we show the breakdown energies of unitarity, perturbativity and stability in figure 3. As can be seen in the figure, the constraining requirement of the 2HDM during the RG running is most often unitarity; however, perturbativity violation usually follows soon after. 12 While unitarity is broken in the evolution of more than 99% of the parameter points, stability is only broken before perturbativity in ∼ 0.02% (2.5%) of the cases in scenario I (scenario II). In order to illustrate the importance of the magnitude of the quartic couplings at the starting scale, figure 3 also shows the summed quartic couplings, i |λ i | vs. the perturbativity and stability breakdown scales in scenario II. As a quantitative example of the i |λ i | dependence, it is clear that to have a viable 2HDM above 10 8 GeV one needs i |λ i | 5 and 10 18 GeV requires i |λ i | 2.
An additional remark, valid for all scenarios, is that some parts of the parameter space of 2HDMs exhibit very large loop corrections to the scalar masses, see figure 4. This is also pointed out in refs. [23,24]. It should be mentioned though that the difference between JHEP02(2019)152  loop corrected and tree-level masses decreases, when considering models which are valid to increasing energy scales. For example, the 1-loop contributions are significantly smaller when working with models that are valid up to energies above ∼ 10 10 GeV; thereby it is sufficient to only consider 1-loop contributions in our analysis.

Scenario I: exact Z 2 symmetry
The parameter space with an exact Z 2 symmetry is severely constrained and the evolution of most parameter points breaks down already at scales a few orders of magnitude higher than the EW scale. Low Higgs boson masses are heavily favored and although we put an upper limit of a tree-level m H ≤ 1 TeV, essentially no heavy scalars, {H, A, H ± } above ∼ 600 GeV, were found to be valid parameter points already at the EW scale due to perturbativity constraints. In figure 5 we show the allowed masses and the correlations between them in more detail. From the figure we see that regions with large ( 500) BSM Higgs masses break down already at TeV scales and that smaller ( 250 GeV) are favored because of the corresponding small quartic couplings. From the leftmost figure we also see that regions with m A ∼ 150 GeV are disfavored, which is due to experimental constraints, and the same applies to m A ≈ m H and m A ≈ m H ± shown in the rightmost panel. To illustrate that larger masses are directly correlated to larger quartic couplings, we show in the middle panel the correlation between Z 3 and m H ± , where the Z 3 coupling is fixed and essentially determined by all other quartic couplings, as can be seen in eq. (2.26). Choosing any other BSM mass on the x-axis would produce a similar plot. Furthermore, in figure 5 one also sees that all masses lie at roughly the same scale, since large mass differences are also connected to large quartic couplings.
In short, we conclude that imposing an exact Z 2 symmetry on the 2HDM requires large quartic couplings to be able to obtain large scalar masses. This results in rapid RGE running, which makes it essentially impossible to find models that are valid all the way up to the Planck scale without some sort of fine-tuning; the highest energy reached in the scan was ∼ 10 13 GeV. Similar findings have been found in ref. [58].

Scenario II: softly broken Z 2 parameter scan
As seen in the previous section, imposing an exact Z 2 symmetry on the general 2HDM makes it hard, if not impossible, to find parameter points that correspond to a valid model all the way up to the Planck scale; furthermore, demanding large BSM Higgs masses makes it even more difficult. The main effect from breaking the Z 2 softly is a positive one; it opens up regions in parameter space that are much more stable during RG evolution. The introduction of m 2 12 helps to give larger BSM masses without increasing the magnitude of the quartic couplings; thus, the strong correlation of large Higgs boson masses and large quartic couplings is reduced when comparing to the case of an exact Z 2 symmetry. This can be seen by comparing figure 5 to figure 6, which shows the larger Higgs masses that are viable in scenario II. Also, experimental constraints are not as constraining in scenario II. Similarly to scenario I, there is still a strong correlation between the BSM masses for parameter space points that are valid all the way up to the Planck scale. In order to illustrate the dependence on couplings to other particles we show plots of the maximum breakdown energy as a function of cos(β − α) and tan β in figure 7. From the figure, we note that the region of valid parameter points in tan β is larger in scenario II compared to scenario I. Similar plots for all the quartic couplings in the generic basis, as well as the Higgs basis, are collected in appendix A. Although it has been checked that the specific Yukawa Z 2 symmetry type does not matter for our conclusions about Z 2 breaking effects in this case, the plots in figure 7 differs quite a bit if one were to consider a type-II symmetry instead of a type-I. In addition, one should then also include constraints from B-physics such as b → sγ.
Note that there are irregular patterns in the regions where some parameter points are valid all the way to the Planck scale in figures 6, 7 and also 9 below; the scale of nearby bins differ several orders of magnitude. This effect is purely statistical in nature and a larger data sample and bin choice would smooth out the plots.
To illustrate the difference induced by allowing soft Z 2 breaking, a histogram of parameter points sorted by breakdown energy is shown in figure 8 for both scenario I and II. If we compare the different scenarios, we see that breaking the Z 2 softly opens up an entire new region of the 2HDM, where the model is theoretically viable all the way to the Planck scale.
A comparison of the 2-loop and 1-loop breakdown scales is also shown for both scenario I and II in figure 8. In general, using the 2-loop RGEs increases the breakdown energy scale, Λ, with a factor of O(1) to O(10); although there are a few parameter points with a larger difference, usually related to the point being close to stability boundaries. This is a consequence from the fact that RG running of the parameters slows down when going to 2-loop RGEs, which can also be seen in the example scenario in figure 2.

Scenario III: hard scalar Z 2 breaking
In the last section, we studied the scenario of a soft broken Z 2 symmetry. A softly broken Z 2 does not spread in the RG evolution in that no further Z 2 breaking parameters are generated; since the mass parameters do not enter the RGEs for the dimensionless parameters. Now, we break the symmetry hard by having small non-zero λ 6 , λ 7 that can have major implications for the RG evolution. In general, if there is no kind of fine-tuning present, JHEP02(2019)152 non-zero λ 6 , λ 7 speed up the running of the other quartic couplings. At 2-loop order these Z 2 breaking parameters induce Z 2 breaking in the Yukawa sector as well, giving rise to FCNCs at tree-level.
Since we expect the induced FCNCs to be very suppressed, we here ignore if the parameter points break stability to be able to generate as large FCNCs as possible. For simplicity, we also ignore the unitarity constraint; it is usually violated very close to the perturbativity breakdown energy and thus has minor implications for the analysis. However, even with these relaxed conditions, we find that no sizable FCNCs are generated, as we will now show.
To see the dependency on the magnitude of λ 6,7 , the breakdown scale of perturbativity as functions of λ 6 and λ 7 is shown in figure 9. From that plot, we deduce that to have a model that is good all the way up to 10 18 GeV one needs λ 6,7 0.1 and models that have λ 6,7 1 break down already at scales ∼ 10 5 GeV.

JHEP02(2019)152
Example from scenario III  Figure 10. Evolution of the parameter point in eq. (4.1); which is an example point from scenario III (hard broken Z 2 in the potential).

2-loop
To investigate the FCNCs that are induced, we plot the λ D 23 element at the perturbativity breakdown energy, which is the largest element in more than 99% of the cases, as a function of λ 6,7 (m t ) in figure 9. Even though, as already mentioned, there are no strict limits on the FCNCs at energies which have not been probed experimentally, we assume that the values of the Yukawa couplings should not be widely different at disparate scales. Therefore, we will use a generic limit on the non-diagonal Yukawa elements as a measure of the induced FCNCs. A discussion of neutral meson mixings can be found in [16] and we will use the resulting limit of sizable FCNCs to be λ F i =j 0.1. As will be presented below, we see that the generated Yukawa element is not problematic in this scenario; the model breaks down before any significant FCNCs are generated. In a way this is natural since the hard Z 2 breaking in the potential only affects the Yukawa sector at 2-loop level, whereas the other potential parameters are affected already at 1-loop level.
As an example of this scenario, we plot the evolution of the parameter point . For completeness, it should be noted that stability and unitarity are violated at ≈ 5 TeV and ≈ 10 11 GeV respectively for this point, but as argued above we ignore this since we are looking for maximal Z 2 -breaking effects.

Scenario IV: Yukawa alignment
The ansatz of flavor alignment in the Yukawa sector as in eq. (2.17) is not stable under RG running, except for the Z 2 symmetric values of the a F parameters in table 1. In this section we investigate how large the generated FCNCs can become by varying these a F parameters. Since the flavor alignment is an example of a Z 2 breaking, it will also spread to the scalar sector and induce non-zero λ 6,7 , which could be a problem in the evolution of the potential.
In scenario IV.a we separately vary each a F one-by-one, while keeping the other fixed, with the scalar potential in eq. (3.10), see figure 11. The Z 2 symmetric values, a U = cot β and a D = cot β or − tan β, are clearly visible and are displayed as vertical lines. As expected, the up-sector is much more sensitive because of the large top Yukawa coupling. sizable non-diagonal FCNCs can be generated by deviating from the Z 2 symmetric values in both the up and down sector. We especially note that, perturbativity in the scalar potential easily breaks down for large deviations of a U 1 in the up sector.
When varying a D , the maximum generated FCNC occurs for a value a D ≈ 5/a U , where a U = cot β. In this region, λ D 23 can become very large, O 10 2 , without perturbativity breaking down. The phenomenon that the maximum occurs when a D ∝ 1/a U is general and is more clearly visible in figure 14 below.
To illustrate the RG evolution, an example parameter point from scenario IV.a with the alignment coefficients a U = 0.8, a D = a L = cot β ≈ 0.4 is shown in figure 12. Stability breaks down for this parameter space point already at 1 TeV, but as before we ignore that.

JHEP02(2019)152
Example from scenario IV.a  In scenario IV.b, we perform a random parameter scan of a 2HDM with a softly broken Z 2 potential, just as in scenario II, but fix tan β = 2. Furthermore, here we also vary the alignment parameters in the quark sector, a U,D . This way we get a quantitative estimate of the induced Z 2 breaking parameters λ 6,7 and λ F i =j . As discussed earlier, we only keep parameter points that satisfies all theoretical constraints during the entire RG running up to 10 10 GeV.
The induced λ 6,7 are shown for a scan of 10 4 surviving parameter points in figure 13. It is clear that the generated Z 2 breaking parameters are more sensitive to variations in a U , because of the large top Yukawa coupling. Note the vertical area around a U = cot β = 1/2 which gives the smallest λ 6,7 .
As previously mentioned, it is the λ D 23 = ρ D 23 v 2 /(2m s m b ) which is the largest generated non-diagonal element in more than 99% of all cases; therefore we use it as a measure of the induced FCNCs in figure 14. The dark regions, where the smallest FCNCs are generated, arise since the relations a U = a D and a D = −1/a U correspond to an aligned Yukawa sector that is 1-loop stable under RG evolution [14]. Even though these relations imply that there is a Z 2 symmetric Yukawa sector in one particular basis, at 2-loop order the scalar and Yukawa sectors mix and induce FCNCs if the a F do not take on the correct Z 2 symmetric values set by the scalar potential. We also note the region a D ≈ 10/a U , where the maximum FCNCs are generated; similar to what was seen when varying a D separately in figure 11.

Conclusions
We have derived the complete set of 2-loop RGEs for a general, potentially complex, 2HDM and implemented them in the C++ code 2HDME [26]. Using this software, we performed a RG analysis of the CP-conserving 2HDM and investigated the parameter space with different levels of Z 2 symmetry (exact, soft or hard breaking). The case of a softly broken Z 2 symmetry is the most unconstrained type of 2HDM, since it allows for large scalar masses, larger values of tan β and we found regions in the parameter space with models that are valid all the way up to the Planck scale.
Breaking the Z 2 symmetry hard at the EW scale induces dimensionless Z 2 breaking parameters during the RG running which can be potentially dangerous and severely limit the validity of the model. We have looked at two scenarios, where we have either broken the Z 2 symmetry hard in the scalar potential, by having non-zero λ 6,7 , or in the Yukawa sector, by making a flavor Yukawa ansatz at the EW scale.
We have found that for λ 6,7 0.1 at the EW scale, it is possible to find models that are valid all the way up to the Planck scale at 10 18 GeV; while models break down already at ∼ 10 5 GeV if λ 6,7 1. The generated FCNCs, that spread from the scalar sector first at 2-loop order, are however heavily suppressed.
A flavor ansatz can be a viable solution to the FCNC problem, even though it is not protected during RG evolution. Deviations of a F from the Z 2 symmetric limit, a D = a U = tan β, a D = a U = cot β or a D = −1/a U , should however not be too large. The up sector is especially sensitive when it comes to induced λ 6,7 ; which contribute to a rapid growth of the quartic couplings. In the down sector, the induced FCNCs pose the primary limitation. Finally we observe that there is a region in the flavor alignment parameter space, a D ∝ 1/a U , which gives the maximum FCNCs. The origin of this effect is unclear and requires further study. Here, we give the conditions for the tree-level scalar potential to be bounded from below, as worked out in refs. [37,38]. When working out these conditions, refs. [37,38] constructed a Minkowskian formalism of the 2HDM that uses gauge-invariant field bilinears, These can be used to create a four-vector r µ = (r 0 , r); where one can raise and lower the indices as usual with the flat Minkowski metric η µν = diag(1, −1, −1, −1). In this formalism, the scalar potential is conveniently written as −Im (Z 6 + Z 7 ) Im (Z 5 ) −Z 4 + Re (Z 5 ) Im (Z 6 − Z 7 ) The scalar potential is bounded from below if and only if all of the below requirements are fulfilled: • All the eigenvalues of Λ ν µ are real.
• There exist four linearly independent eigenvectors; one V (a) for each eigenvalue Λ a .
• The eigenvector corresponding to the largest eigenvalue is timelike, while the others are spacelike,

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