$Z_2$ breaking effects in 2-loop RG evolution of 2HDM

We investigate the effects of a $Z_2$ symmetry in the CP-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 $Z_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 $Z_2$ symmetric 2HDM is very constrained, a soft breaking of the $Z_2$ symmetry extends the valid parameter space regions. The effects of a hard $Z_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 $Z_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 $Z_2$ breaking in the Yukawa sector at most gives moderate $Z_2$ breaking in the potential; whereas the FCNCs can become quite sizable far away from the $Z_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; 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 ref. [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) analyses 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. 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 [7,8]. Although the flavor alignment solves the problem of FCNCs at a specific energy scale, the ansatz is not stable during RG evolution [9][10][11]; which also induces Z 2 breaking parameters in the scalar potential.
At 2-loop order, the RG analyses becomes more complicated since there is a stronger coupling between the scalar and Yukawa sector; 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 ref. [18][19][20], the authors employ 2-loop RGEs in their analyses, but on a CP-conserving 2HDM, with a 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 [22].
The structure of this paper is as follows: In section 2, we review the necessary 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. [23][24][25] 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, v = 2m W /g 2 ≈ 246 GeV 2 .
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 [24]; 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 [23,26], 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, 3 Defined by which one is the lighter one; m h ≤ mH .
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 [24] Thus from inspection of the Higgs potential in eq. (2.7), it follows that Y 1 , Y 2 , Z 1 , Z 2 , Z 3 , Z 4 are invariant, while 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 (2.15) where κ F are the diagonal mass matrices, 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 [7], 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 analyses to real coefficients. As is well known, the alignment ansatz is not stable under RG evolution [9][10][11], 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 [27] 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 [12] 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 "sizeable" 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. [24] for details. 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 .
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 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 [23,[28][29][30] 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 [31] Im λ * 5 (m 2 12 ) 2 = 0, and 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 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. [31] and is summarized in table 2, which we will make use of in upcoming sections.  Table 2. Hybrid basis for a softly broken Z 2 symmetric CP-invariant 2HDM from ref. [31]. Imposing an exact Z 2 symmetry fixes Z 7 according to eq. (2.28), withm 2 = 0.
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 [31] 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 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 analyses is thus useful when looking for potential fine-tuning, 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 Nloop matching [19]. In this work, we are however 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.
Even though we will not impose any precision observable constraints on the model, we will impose theoretical constraints, namely perturbativity, unitarity and stability, that are needed to arrive at a consistent model.

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 analyses useless. The couplings can 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 and thus are the ones affected by the unitarity limit. In ref. [32], 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 ref. [33,34], 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.

Derivation of 2-loop RGEs
Although there are some works using 2-loop RGEs for analyzing 2HDMs [18][19][20], 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 [39][40][41]. These were then supplemented with the 2-loop RGEs of massive parameters in ref. [42], 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 ref. [39][40][41][42] 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. [21,38] 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. [43] 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. [44]. 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 [22].

Technical details of numerical code
Evolving the 2HDM in energy corresponds to solving a system of 129 7 coupled ordinary differential equations. We have the developed the C++ code 2HDME [22] to perform this task. It uses GSL [45] to solve the RGEs and the library Eigen [46] for linear algebra operations.
To match the 2HDM to the µ dependent parameters of the SM at the top quark pole mass scale, M t = 173.4 GeV [47], we used the following input as boundary conditions: 6 Refs. [19,20] are using 2-loop RGEs derived by SARAH [35,36] and ref. [18] is using PyR@TE [37]. As have been pointed out in ref. [21,38], 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. [38] 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.
• The MS fermion masses are used to fix the Yukawa matrix elements in the fermion mass eigenbasis. We use the ones from ref. [48]: The numerical values λ = 0.22453, are extracted from the PDG [50].
To calculate the 1-loop pole masses of the Higgs bosons, we made use of the fortran code SPheno [51,52], together with model files generated by SARAH [35,36].

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 fixed to be proportional to the κ F matrices, either by a Z 2 symmetry or flavor ansatz. One can then go to the fermion weak eigenbasis in the scalar generic basis with Here, we made the arbitrary choice to put the CKM matrix in the down sector by setting V D L = V † CKM and V U L = I. This choice does not affect the results in any way.
• The parameters in the generic basis, i.e. the set {g i , v a , η F,0 a , m 2 ij , λ i }, are evolved according the RGEs. This includes the VEVs, v a ≡v a v, which obey the same RGE as the scalar fields, i.e.
with the consequence 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. We define Λ as the lowest energy scale where either perturbativity, unitarity or stability is violated.

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

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 the RG evolution algorithm 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 9 . 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 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. 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; 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: 9 An exact Z2 symmetric 2HDM also implies CP-invariance. 10 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.
• 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 independently in the ranges a U ∈ [0.01, 10] and a D ∈ [−100, 100].
• 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]. For this scenario, we evolve up to µ = 10 10 GeV and only save models where neither perturbativity, unitarity or stability is violated.    Before we investigate each respective scenario in more detail, we study which of the constraints that are violated and consequently gives rise to the breakdown energy Λ. As can be seen in figure 3, the most constraining requirement of the 2HDM during the RG running is most often unitarity; however, perturbativity violation usually follows soon after 11 . While unitarity is broken in the evolution of more than 99% of the parameter points, stability is only broken before perturbativity in ∼ 0.02% (0.2%) of the cases in scenario I (scenario II). Also visible in figure 3 is the summed quartic couplings, i |λ i | as a function of the perturbativity and stability breakdown in scenario II; from which it is clear that to have a viable 2HDM above 10 8 GeV one needs i |λ i | 5.
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 ref. [19,20]. It should be mentioned though that the difference between loop corrected and tree-level masses decreases, when considering models which are valid to increasing energy scales; the 1-loop contributions are significantly smaller when working with models that are valid up to energies above ∼ 10 10 GeV.

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. Regions with large ( 500 GeV) BSM Higgs masses break down already at TeV scales, see figure 5. Smaller masses ( 250 GeV) are favored because of the corresponding small quartic couplings. That larger masses are directly correlated to larger quartic couplings can be seen in figure 6. There, we show the Z 3 coupling, that is essentially fixed and determined by all other quartic couplings, as a function of m H ± . 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, one concludes that imposing an exact Z 2 symmetry on the 2HDM forces large quartic couplings to get 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.

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 and 6 to figure 7, which shows the larger Higgs masses that are viable 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. Furthermore, the region of valid parameter points in tan β is larger, as shown in figure 9. 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 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, section 4.2, 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, 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. The breakdown scale of perturbativity as functions of λ 6 and λ 7 is shown in figure 10. From that plot, we deduce that having λ 6,7 ≥ 0.2 is problematic if one wants viable models above ∼ 10 8 GeV.
To investigate the FCNCs that are induced, we plot the λ D 23 element, which is the largest one in more than 99% of the cases, as a function of λ 6,7 at the perturbativity breakdown energy in figure 10. If we take the limit of "sizeable" FCNCs to be ∼ 0.1 [12], we see that the generated Yukawa element is not problematic in this scenario. However, as previously mentioned, having λ 6,7 ≥ 0.2 gives no valid points above ∼ 10 8 GeV and the model thus breaks down before any significant FCNCs are being 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  Figure 11. Evolution of the parameter point in eq. (4.1); which is an example point from scenario III (hard broken Z 2 in the potential).

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 independently vary each a F with the scalar potential in eq. noted that, stability and 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. That the maximum occurs when a D ∝ 1/a U is a general phenomena that is more clearly visible in figure 15.
The evolution of an example parameter point from scenario IV.a with the alignment coefficients a U = 0.8, a D = a L = cot β ≈ 10 −1 is shown in figure 13. Stability breaks down for this parameter point already at 1 TeV.
In scenario IV.b, we perform a random parameter scan of a softly broken Z 2 2HDM, just as in scenario II, but fix tan β = 2. We keep only parameter points that satisfies all theoretical constraints during the entire RG running up to 10 10 GeV. Furthermore, here, we also vary the alignment parameters in the quark sector, a U,D . Thus we get a quantitative estimate of the generated Z 2 breaking parameters.
The generated λ 6,7 are shown for a scan of 10 4 surviving parameter points in figure 14. It is clear that the generated Z 2 breaking parameters are more sensitive of 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 15. 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 [10]. Even though these relations imply that there is a Z 2 Yukawa sector in one particular basis, at 2-loop order the scalar and Yukawa sector mix and induces 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 ≈ 8/a U , where the maximum FCNCs are generated; similar to what was seen when varying a D independently in figure 12.

Conclusion
We have derived the complete set of 2-loop RGEs for a general, potentially complex, 2HDM and implemented them in the C++ code 2HDME [22]. Using this software, we performed a RG analyses of the CP-conserving 2HDM and investigated the parameter space with different levels of Z 2 symmetry (exact, softly broken or explicitly broken). The case of a softly broken Z 2 symmetry is the most unconstrained type of a 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 explicitly 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 explicitly by having non-zero λ 6,7 or made a flavor Yukawa ansatz at the EW scale.
We have found that having λ 6,7 ≥ 0.2 at the EW scale is problematic if one wants a model that is theoretically valid up to energy scales above 10 8 GeV. 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 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 quartic couplings in the random parameter scan of scenario II with a softly broken Z 2 symmetry, described in section 4.2, as functions of breakdown energy in RG evolution are shown in the generic basis in figure 16 and in the Higgs basis in figure 17. The color coding is the maximum breakdown-energy from perturbativity, unitarity or stability in each bin. . Figure 16. Breakdown energy as a function of quartic couplings in the generic basis in the parameter scan of scenario II.

B Tree-level unitarity conditions
The tree-level unitarity conditions for a general 2HDM have been worked out in ref. [32]. There, they work out the following scattering matrices:

C Stability
Here, we give the conditions for the scalar potential to be bounded from below, as worked out in ref. [33,34]. When working out these conditions, ref. [33,34] 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 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,