Global fits of the two-loop renormalized Two-Higgs-Doublet model with soft $Z_2$ breaking

We determine the next-to-leading order renormalization group equations for the Two-Higgs-Doublet model with a softly broken $Z_2$ symmetry and CP conservation in the scalar potential. We use them to identify the parameter regions which are stable up to the Planck scale and find that in this case the quartic couplings of the Higgs potential cannot be larger than 1 in magnitude and that the absolute values of the S-matrix eigenvalues cannot exceed 2.5 at the electroweak symmetry breaking scale. Interpreting the 125 GeV resonance as the light CP-even Higgs eigenstate, we combine stability constraints, electroweak precision and flavour observables with the latest ATLAS and CMS data on Higgs signal strengths and heavy Higgs searches in global parameter fits to all four types of $Z_2$ symmetry. We quantify the maximal deviations from the alignment limit and find that in type II and Y the mass of the heavy CP-even (CP-odd) scalar cannot be smaller than 340 GeV (360 GeV). Also, we pinpoint the physical parameter regions compatible with a stable scalar potential up to the Planck scale. Motivated by the question how natural a Higgs mass of 125 GeV can be in the context of a Two-Higgs-Doublet model, we also address the hierarchy problem and find that the Two-Higgs-Doublet model does not offer a perturbative solution to it beyond 5 TeV.


I. INTRODUCTION
After the discovery of a scalar particle at the LHC [1,2], one of the next questions is whether this is the one Higgs particle predicted by the Standard Model (SM) or whether there are more generations of SU (2) doublets, like it is the case for the fermions. In that sense, the simplest and most straightforward extension of the SM would be the addition of another Higgs doublet, the so-called Two-Higgs-Doublet model (2HDM) [3][4][5]. Furthermore, the measured mass of this new scalar [6] is a peculiar value for the SM: it tells us that the Higgs potential of this model cannot be stable up to very high energy scales [7,8]. However, there is the possibility that the electroweak vacuum may just end up being metastable. So either one has to believe that we live in a metastable universe and then there is no need of new physics beyond the SM, or one has to introduce an additional mechanism to stabilize the Higgs potential. The latter could for instance be achieved by the heavier scalars of the 2HDM. This model might be realized as an intermediate "effective" theory which describes physics at energy scales between the electroweak scale µ ew of order 10 2 GeV and some higher scale µ high . Beyond the latter, a more comprehensive model would be needed to describe "physics beyond the 2HDM". An upper bound on µ high is the Planck scale µ Pl ≈ 10 19 GeV, at which gravitational effects become non-negligible in a quantum field theory framework. Large scale differences between µ ew and µ high bring along hierarchy problems like the finetuning of the 125 GeV Higgs mass, which could be resolved by mechanisms of the "complete" models, but are usually neglected in the effective models. Still one could ask to what extent the 2HDM could possibly mitigate the Higgs mass hierarchy problem and whether it might even be valid up to Planck scale without requiring any other New Physics.
For renormalization group studies, especially the role of Higgs self-couplings is crucial and has been studied in the literature, in the SM (see for instance [7,37]) as well as in the 2HDM [5,[38][39][40][41][42][43], because these quartic couplings tend to destabilize the Higgs potential at some µ high . Since a break-down of stability would mean that our theory would lose validity beyond a certain scale, we want to impose a stable Higgs potential beyond µ ew as a constraint on all couplings. Recently, the impact of stability up to the Planck scale on the parameters of the aligned 2HDM was discussed in [44]. If one wants to solve the Higgs mass fine-tuning problem, one has to guarantee the cancellation of quadratic divergencies of higher order Higgs mass correction terms. The corresponding conditions that need to be fulfilled are called "Veltman conditions" [45] in general, and in the context of the 2HDM also "Newton-Wu conditions" [46]. They have been analyzed at one-loop level [47][48][49] and even leading two-loop contributions have been taken into account in type II [50][51][52]. A recent idea was to only relax the cancellation of the generically large contributions of quadratic divergencies instead of imposing the strict cancellation using Veltman conditions [53].
In this article, we want to improve available results concerning two main aspects: We perform global parameter fits including the most up-to-date ATLAS and CMS results, rather than only using a handful of benchmark scenarios, which might not cover the whole spectrum of interesting features. Secondly, we go beyond leading order precision by employing twoloop renormalization group equations (RGE) in order to analyze vacuum stability of the 2HDM scalar potential. Moreover, we want to make use of the framework of next-to-leading order RGE to find out to what extent Veltman conditions can be fulfilled in the 2HDM.
In the following, we first want to make the reader familiar with the model in Section II, and introduce in Section III its theoretical and experimental constraints and the numerical setup that we use. Then we are ready to compare leading and next-to-leading order renormalization group equations for a benchmark scenario in Section IV A. We go on examining the quartic couplings varying the stability cut-off scale in global fits without experimental inputs in Section IV B. We also address the question of which upper limit to use for the unitarity condition from the perspective of renormalizability. Taking into account experimental data, we analyze the results of global fits to the physical parameters at µ ew and µ Pl in Section IV C. The hierarchy problem is discussed in Section V, before we conclude in Section VI. Explicit expressions for the one-loop and two-loop RGE can be found in the Appendix.

II. MODEL
The Two-Higgs-Doublet model with a softly broken Z 2 symmetry is characterized by the following scalar potential: where Φ 1 and Φ 2 are the two Higgs doublets. In the following, we will use two sets of parameters: the eight potential parameters from Eq. (1), which we assume to be real, and the physical parameters consisting of the vacuum expectation value v, the CP -even Higgs masses m h and m H , the CP -odd Higgs mass m A , the mass of the charged Higgs, m H + , the two diagonalization angles α and β, and the soft Z 2 breaking parameter m 2 12 . The first two physical parameters can be treated as fixed by measurements, assuming that the 125 GeV scalar found at the LHC is the lighter CP -even Higgs. Instead of α and β we will use the combinations β − α and tan β, since they can be directly related to physical observables. The measurements of the light Higgs couplings to fermions and bosons are compatible with the SM, such that the 2HDM is pushed towards the so-called alignment limit [30,54], in which β − α = π/2. In this limit, it has recently been shown that CP violation in the scalar potential of Z 2 symmetric models with a soft breaking term is strongly suppressed [55], which qualifies our above assumption that the potential parameters are real. The masses of the heavy scalars could in general even be lighter than 125 GeV, and are not necessarily in the decoupling limit [4] (which itself is a limiting case of the alignment limit). In the following we will consider them to be in the range between 130 GeV and 10 TeV, that is heavier than the region where the 125 GeV scalar was found, yet still in the TeV range, which will be accessible by future colliders.
Neglecting the first two generations of fermions, the Yukawa part of the 2HDM Lagrangian is In the above Lagrangian, the top quark only couples to Φ 2 by convention; its Yukawa coupling is related to the SM value Y SM t by Y t = Y SM t / sin β. Without breaking the Z 2 symmetry in the Yukawa sector, there are only four possibilities to couple the Higgs fields to the bottom quark and tau lepton at the tree-level. They are called type I, type II, type X or "lepton specific" and type Y or "flipped"; in Table I we show the corresponding Higgs field assignments. Type II is of special interest, as it contains the Higgs part of supersymmetric models. As soon as we consider any one of the above types, only three Yukawa couplings remain as free parameters, and we can speak of Y t , Y b and Y τ without any ambiguity. Type I Type II Type X ("lepton specific") Type Y ("flipped")

III. CONSTRAINTS AND SET-UP
We will apply the following sets of constraints on the parameter space: On the theoretical side, the positivity of the Higgs potential [56] and the unitarity of the eigenvalues of the Φ i Φ j → Φ i Φ j scattering matrix [57] are imposed at all scales and vacuum stability [58] at the electroweak scale. Moreover, we make sure that the quartic couplings λ i (i = 1, 2, 3, 4, 5) and the Yukawa couplings Y i (i = t, b, τ ) do not run into non-perturbative regions. On the experimental side, electroweak precision observables, the branching ratio Br(B → X s γ), the mass difference in the B s system and light and heavy Higgs searches constrain the 2HDM parameters at the electroweak scale. For a detailed description on the various constraints, we refer to [14] and [20] except for the experimental Higgs data: here, we use the most upto-date ATLAS and CMS publications and pre-prints [59][60][61][62][63][64][65][66][67][68][69][70][71][72], applying the narrow width approximation. We do not make use of (semi-)tauonic B decay observables, which would only be relevant in type II [73], because the existing tension between the measurements can only be explained in a 2HDM with explicitly broken Z 2 symmetry [74]. The SM parameters will be fixed to their best fit values [75]; for the SM Yukawa couplings in the MS renormalization scheme at the scale m Z we take Y SM t = 0.961, Y SM b = 0.0172 and Y SM τ = 0.0102. While variations of the strong coupling α s (m Z ) within the 3σ allowed range have no effects on the outcome of our fits, varying the input for m t (m Z ) can have an impact on a specific parameter region like mentioned in [42]. However, we observe that these effects are imperceptible in the results of our global fits.
The two-loop RGE have been obtained with the publicly available package PyR@TE [76]; we neglect all Yukawa couplings except for the top and bottom quarks and the τ lepton. The observables have been calculated with the help of Zfitter [77][78][79], FeynArts [80], FormCalc [81], LoopTools [82], HDECAY [83][84][85], FeynRules [86] and MadGraph5 [87]. The frequentist fits are performed with the CKMfitter package [88]. For the fits involving experimental constraints we use of the naive definition of the p-value (Wilks' theorem) [89]. If not stated differently, exclusion limits are meant to be at the 2σ level, which roughly corresponds to the 95% confidence level.
Since we want to discuss various values for the scale µ in the following, we want to define our notation: The scale range for the running quantities lies between the electroweak scale µ ew = m Z at the lower end and the Planck scale µ Pl =10 19 GeV at the upper end, as mentioned in the introduction. If -starting with a given set of parameters at µ ew and evolving to higher energy scales -one of the theoretical constraints is violated, we denote this breakdown of stability as µ st . When discussing the hierarchy problem, it might be useful to introduce a cut-off scale µ nat , which a priori does not need to be the same as µ st . Furthermore, we want to introduce the parameter t (X) = ln(µ (X) / GeV) as the usual logarithmic scale.

IV. RENORMALIZATION AT NEXT-TO-LEADING ORDER
One can find a plethora of leading order RGE for different realizations of a 2HDM in the literature [5,[90][91][92]; however, we failed to find a complete set for a 2HDM with soft Z 2 breaking, including the mass parameters, so we list the leading order (LO) and next-toleading order (NLO) expressions in the Appendix. Before scanning over the whole parameter space with our fitting set-up, we want to explain some features of the 2HDM RGE looking at a representative example.

A. A benchmark point
In order to compare the LO and NLO RGE, we choose the scenario H-4 from [20] as benchmark scenario, because all quartic couplings are relatively large already at µ ew . It is defined by m H = 600 GeV, m A = 658 GeV, m H + = 591 GeV, β − α = 0.513π, tan β = 4.28, and m 2 12 = 76900 GeV 2 , and compatible with all experimental measurements so far. The cut-off scale, where one of the quartic couplings becomes non-perturbative, is at 19.5 TeV at LO (dashed lines), and at 82 TeV at NLO (solid lines), see the top left panel of Fig. 1. The Landau poles are at 54 TeV and 3.2 · 10 6 TeV, respectively; the former is shown as a vertical dotted line in Figs. 1 and 2. The fact that the higher order contributions "stabilize" the RG evolution, and thus increase both, cut-off and Landau pole scales, holds for all benchmark points from [20] and is a general feature in the 2HDM: All dominant NLO contributions to the RGE of the λ i , which are cubic in the quartic couplings, come with a negative coefficient and thus mitigate the positive LO contribution coming from quadratic λ i terms (see Appendix). Beneath the total values of LO and NLO running we show the relative difference between LO and NLO RGE | with respect to the scale. For this benchmark point, the relative change of λ 1 and λ 3 , is as large as 10% at around 2 TeV and the difference increases even at a faster rate at higher scales. This is a first hint that the effect of the NLO contribution to the RGE in the 2HDM is non-negligible. One can see that r 3 diverges around 35 TeV due to the fact that at this scale λ NLO 3 turns to 0. A better quantitative measure of the NLO vs. LO RGE is the relative distance δ L 12 defined in [93]: For a dimensionless coupling L, we can define the relative distance between the LO and NLO curves L LO (t) and L NLO (t) in the scale range from t 1 to t 2 as | between LO and NLO expressions for the quartic couplings. The Yukawa couplings and the potential mass parameters are shown on the right. All types look the same except for the b and τ Yukawa couplings.
For H-4, δ λ 1 12 is 38%, if we integrate from m Z to the LO cut-off at 19.5 TeV. To quantify the typical size of δ λc 12 , where λ c is the quartic coupling with the lowest perturbativity violating scale, we checked all benchmark points of [20], and found values between 17% and 45%, which indicates that in general the two-loop corrections are not negligible. 1 In Fig. 1 we do not show the running of the gauge couplings g 1 , g 2 and g 3 , since the two-loop corrections are too small to be visible. Also, the running of the Yukawa couplings is not significantly altered going from LO to NLO. However, due to the different assignment of the Higgs fields in the four types, we start with different values at the low scale (see Table I); that is why we denote the Yukawa couplings in the upper right panel of Fig. 1 as introduced in Eq. (2). (Note that only two of them are non-zero, depending on the type of Z 2 .) Among the mass parameters, m 2 12 changes least if we run to higher scales, which we also observe as general feature of all types. or m 2 22 . Furthermore, we have checked the mentioned benchmark scenarios for fixed point behaviour and do not find any below the perturbativity cut-off.
If we switch to the physical parameter basis, we observe that also the RG running of the mixing angles can be sizable, see left side of Fig. 2. The scale at which β − α hits the alignment limit corresponds to vanishing v and m h , which can be seen in the right panel of Fig. 2, where we show the running of all physical mass parameters. We find that the breakdown of the vacuum expectation value at some scale above µ ew is a general feature and occurs for all benchmark scenarios that we have analyzed; this is also observable in the benchmark points of [41].
After scrutinizing one benchmark point, we want to discuss more general features that can be found in comprehensive fits. Especially the general dependence of µ st on the value of λ i (m Z ) will be an interesting question in the following.

B. Fits without experimental data
A parametrization independent way of setting upper limits to the quartic couplings of the Higgs potential is the requirement that the scattering matrix of Φ i Φ j → Φ i Φ j processes is unitary. This corresponds to the condition that its absolute eigenvalues should be smaller than 16π. The tree-level expressions [57,[94][95][96][97] are widely used theoretical constraints for the 2HDM; however, it seems that these bounds are very conservative. Studies involving higher order corrections have shown that the eigenvalues cannot be larger than 2π in the SM [37], and this bound has been adopted for the 2HDM of type II in [20]. Analyzing maximally allowed cut-off scales can shed light on how well this bound is motivated from the RGE perspective.
In this section, we only want to impose the Higgs potential bounds, regardless of experimental constraints, in order to show the impact of the former on the 2HDM parameters. Since the assumption of having a stable potential affects the potential parameters, we express our results in terms of the five quartic couplings and tan β. (The latter modifies the Yukawa couplings as compared to their SM values, see Table I.) Due to the smallness of Y b and Y τ , their influence on the RGE is very weak and differences between the four Z 2 types are not visible in the λ i planes. The red (dark) regions illustrate the allowed regions, if we take 4π instead and force at least one of the eigenvalues to be larger than 2π in magnitude. All types give the same dependence for the λ i . For tan β, we show the possible regions in type I and X as shaded, and the areas below the dashed lines correspond to type II and Y.
In Fig. 3 we show the dependence of the cut-off scale on the values of quartic couplings and tan β at the electroweak scale, for the two cases that either all eigenvalue moduli are smaller than 2π or that at least one of them is larger than 2π. Our fits show that forcing at least one eigenvalue of the S-matrix to have an absolute value larger than 2π reduces the maximal cut-off scale µ st to be at 5 · 10 6 GeV instead of the Planck scale; if we set at least one eigenvalue modulus larger than 4π, the maximal µ st is at a few TeV. If we want to maintain a stable Higgs potential up to µ Pl , the largest eigenvalue can have a magnitude of at most 2.5 (≈ 0.8π). Naturally, a larger upper bound on the eigenvalues allows for larger quartic couplings. But one can also see that cut-off scales larger than 10 TeV are only allowed for a very narrow range of tan β around 0.7 and -only in type II and Yfor an additional narrow range around 80. While the low tan β scenarios are known to be disfavoured for light 2HDM spectra by flavour observables, we will see in the next section that also the large tan β regions are now excluded in type II. So we can conclude for all types but type Y that assuming µ st > 10 TeV and not too heavy new Higgs states all S-matrix eigenvalues need to be smaller than 2π in magnitude. We will use the upper bound of 2π in the following. Fig. 3 also shows the allowed λ i (m Z ) and tan β(m Z ) intervals for µ st at Planck scale. Roughly speaking, stability up to 10 19 GeV requires |λ i (m Z )| < ∼ 1 and tan β(m Z ) > 1.
In this case, we also observe a lower limit on λ 2 (m Z ), which cannot be smaller than 0.15. In type II and Y, tan β(m Z ) is also limited from above and cannot be larger than 60. We list the precise ranges of the parameters in Table II.   TABLE II. Allowed intervals for the quartic couplings and tan β at the electroweak scale, if we assume stability at µ ew (first line) and up to µ Pl (second line). The bounds on λ 5 (m Z ) give us a handle on the question whether the Z 2 symmetry can be exact with stability up to the Planck scale: following [4], we find that the soft breaking parameter can be written as With our choice of m A (m Z ) >130 GeV, we see that even at µ Pl it can cancel with v 2 λ 5 (m Z ), so an exact Z 2 symmetry is compatible with a Higgs potential stable up to µ Pl . This is in contradiction with the result of [44]. The inclusion of experimental bounds has only very little visible impact on the potential parameters, that is why in the following section we switch to the physical basis.

C. Fits with experimental data
In this section we want to show the impact of the experimental results discussed in Section III on the physical parameter space at the electroweak scale, once assuming a stable scalar potential at µ ew and once for stability up to µ Pl . We put special emphasis on the dependence of mass parameters on the relevant angles in order to investigate how large deviations from the alignment limit can still be.
In Fig. 4, we show the tan β-(β − α) plane for type I on the upper left, for type II on the upper right, for type X on the lower left and for type Y on the lower right. For a stable potential at the electroweak scale (orange) we show the 1σ, 2σ and 3σ allowed regions (the 2σ region is shaded, the 1σ and 3σ contours are defined by the dash-dotted and dashed lines, respectively), and at Planck scale (purple shaded) we only present the 2σ region. With stability at µ ew , tan β is not constrained by any observable. For 2HDM masses below h → γγ signal strength values. We have seen in Section IV B that imposing stability up to µ Pl constrains the quartic couplings; at this point, we want to shed light on the effect on the physical parameters. In Fig. 3, we already observed that tan β is constrained from below in type I and type X and additionally from above in type II and type Y. While in type I and X sizable deviations of β − α from π/2 are still possible, stability at µ Pl basically fixes β − α to the alignment limit in the other two types, with the maximal deviation being 0.012π and 0.013π in type II and Y, respectively (see also Table III). Stability up to the Planck scale is also incompatible with the alignment limit for tan β < 2.6 in type II and Y. The reason why this value is smaller than the one found in [44] is that we use NLO RGE. At leading order, we confirm their result that tan β < 3 is excluded in the alignment limit.
In Fig. 5, we show the dependence of the charged Higgs mass bounds on tan β. In type I the strongest constraint for low tan β values comes from the mass difference in the B s system. The other observables have no visible impact on this plane. The same holds for type X, except for m H + < 300 GeV, where direct Higgs searches additionally cut away low tan β values. For type II and type Y, Br(B → X s γ) yields a lower limit on m H + of around 300 GeV; 2 for large masses and low tan β, the bound from the mass difference in the B s 2 During the finishing stage of this work it was shown that due to a reduction of the theoretical uncertainty TABLE III. Allowed intervals for β − α at the electroweak scale (and its sine and cosine) for all types of Z 2 symmetry, if we assume stability at µ ew (first three lines) and up to µ Pl (lines four to six). system is stronger. In case of the type II we also find that a light charged Higgs is excluded for large tan β values; for instance if tan β = 30, we obtain m H + > 700 GeV. This is an effect only visible in a global fit: for large tan β, heavy Higgs searches (mainly the tauonic decays) exclude light m H and m A . Electroweak precision data, however, are not compatible with too large mass splittings between the heavy neutral and the charged Higgs particles, so also the charged Higgs cannot be too light if tan β is large. This also qualifies that we did not use data from (semi-)tauonic B decays, which would give a weaker bound on the same corner of the type II plane. Type Y also features this constraint from the neutral Higgs searches, but it is much weaker and only visible for very large values of tan β because the τ and b couplings to H and A cannot be enhanced simultaneously. Requiring stability up to µ Pl gives almost the same regions as at µ ew , only that tan β gets constrained at the borders to stay within the limits from Table II. While the charged Higgs searches mainly depend on tan β, neutral Higgs signals strongly depend on the deviation from the alignment limit, i.e. the actual value of β − α. Therefore, we show in Fig. 6 the allowed regions in the (β − α)-m H and (β − α)-m A planes. In all types it is visible, how the direct searches for heavy Higgs particles exclude large regions for m H/A < 600 GeV, where sizable deviations from the alignment limit would be allowed by the stability bounds. Like for the charged Higgs mass, there is no lower limit for the neutral masses in type I and X. For type II, however, the fits yield lower 1σ bounds of around 250 GeV for the heavy CP even Higgs mass and around 350 GeV for the CP odd Higgs mass; for type Y, we obtain 1σ limits of m H > 220 GeV and m A > 320 GeV. In the alignment limit, one even gets a 2σ bound of 230 GeV for m A in type II and Y. This lower limit on the pseudoscalar mass excludes an exact Z 2 symmetry if we require stability up to µ Pl in those two types, cf. Eq. (3).
The "lower branches" of Fig. 4 correspond to the regions left of the alignment limit in Fig. 6 and feature m H/A values below 380 GeV in the types II, X and Y. For all types we observe that for neutral masses above 600 GeV the deviation of β − α from π/2 can be 0.05π at most due to the stability bound. If we additionally impose stability up to the Planck scale, this bound holds even for masses down to 240 GeV. For this stability assumption, we also confirm the result of [44] that the heavy masses need to be almost degenerate; in our fits we find an upper limit of 45 GeV on the absolute mass splittings.
of Br(B → X s γ) the lower type II bound on m H + is increased to 480 GeV at 95% C.L. [99].

V. THE HIERARCHY PROBLEM
As we have already mentioned, there is a large scale difference between the Planck scale and the scale at which electroweak symmetry breaking occurs. This gap leads to the hierarchy problem of the Higgs mass: if loop corrections can be of order of µ Pl , why do they cancel each other almost perfectly, such that the Higgs mass is 17 orders of magnitude smaller? The cancellation of these mass corrections to retain a naturally light m h was first proposed by Veltman [45], therefore also referred to as "Veltman conditions", and was first applied at leading order to the 2HDM by Newton and Wu [46]. Unlike in supersymmetry, in the 2HDM there is no mechanism which naturally accounts for these cancellations. Still, there might be a hidden symmetry of which we are not aware, so nevertheless it is interesting to address this question. In the framework of the 2HDM, this hierarchy problem does not only affect m h but in principle also the other scalar masses, if they are not in the decoupling limit. However, since we do not know whether the heavier scalars are decoupled or not, we will only discuss the hierarchy problem of the already discovered 125 GeV scalar in the following. The largest one-loop contributions to the Higgs mass come from terms that are quadratic in µ nat , if we assume that the 2HDM is valid up to a given scale µ nat and use this scale as a cut-off. Leading higher order contributions get an additional factor of [ln(µ nat /µ ew )] n , where n + 1 is the number of loops. If µ nat is large enough, the logarithmic factor might compensate for the loop suppression, and the power series of the higher order corrections no longer converges. So requiring the cancellation of the first order Higgs mass correctionlike often applied in the literature [39,48,49,53] -is not sufficient if we do not know about the higher order terms. Only if we assume perturbativity of the power series, we can make a valid statement about whether the Higgs mass at the electroweak scale can be natural in the 2HDM or at least whether the hierarchy problem can be mitigated. This assumption of perturbativity is analogous to the one applied above on the Yukawa and quartic Higgs couplings. All higher order leading logarithm mass corrections proportional to µ 2 nat are given by As described in [100], especially for low cut-off scales the power series can be perturbative. However, we need to be careful to keep the leading logarithm sufficiently large with respect to the lower powers in the logarithm assuming that the leading logarithm gives the largest contribution. The leading coefficient function can be derived from the one-loop Higgs mass corrections and reads as In order to easily obtain the leading logarithm contributions to higher orders, we use the recursive formula derived by Einhorn and Jones [101], relating the coefficient functions f n+1 to f n and the running of the couplings: This recursive relation is based on the following assumptions: the new theory has only one mass scale (m h ), and the logarithmic factor has to be large enough to suppress the terms with lower powers of logarithms. Two-loop effects on the Veltman condition have already been applied to the 2HDM of type II using this approach [50][51][52]; the authors found that the Higgs mass hierarchy problem can be ameliorated. An obvious choice of µ nat as cut-off would be the breakdown of one of the stability constraints µ st , so we will use it for the moment. In order to analyze whether we can make a statement on the Higgs mass naturalness in a 2HDM which is based on a reliable perturbation series, we want to define k n as the ratio of the n-th correction term of Eq. (4) to the one of order n − 1: The allowed size of the one-loop coefficient f 0 (λ i , Y i , g i ) of the Veltman condition series depends on the naturalness cut-off µ nat . Without experimental constraints (light blue shaded), |f 0 (λ i , Y i , g i )| can be as large as 75. With the inclusion of the measurements (orange shaded), it gets strongly constrained to be smaller than 6. More or less independently of taking into account experimental data, we obtain an upper bound on µ nat of a few TeV. The plane shows the type I fit, but the other Z 2 types look the same.
Apart from the obvious logarithmic dependence on µ nat , the k n depend on the cut-off scale also indirectly: the latter determines which values for the λ i are allowed, see Fig. 3. Now we can re-write Eq. (4) as Only if we impose a small value of the leading order coefficient function and a sufficiently small number for k 1 , we can guarantee a perturbatively stable mitigation of the hierarchy problem of m h , also assuming that the k for > 1 are not too large. Note that if we choose f 0 (λ i , Y i , g i ) to be exactly 0, k 1 diverges. If we constrain the first two factors k 1 and k 2 to be smaller than 1 in magnitude and that |δm 2 h | < m 2 h , we observe negative k 1 and k 2 in most cases, independently of the type of Z 2 symmetry. This indicates that the series is alternating, which in turn means that -except for pathological scenarios -a suppression of the first two k i factors should be sufficient to make the series relatively robust with respect to pertubativity. Cutting the series in Eq. (5) after the second term (i.e. setting k 3 = 0), we find that the maximal µ nat is in the TeV range for all types, depending on the value we choose for f 0 (λ i , Y i , g i ). In Fig. 7, we show this dependence for type I; all other types look very similar. Taking only Higgs potential constraints (as in Section IV B), leaves a lot of freedom for f 0 (λ i , Y i , g i ), which indicate large cancellations between the leading order contribution and higher order terms. This calls into question our assumption that the series can be cut after the second term. The inclusion of the experimental results (cf. Section IV C) limits the choice of f 0 (λ i , Y i , g i ) to be of order 1, and thus presumably stabilizes the perturbative series. In both cases, however, the maximal µ nat is a few TeV for very small values of f 0 (λ i , Y i , g i ).
Finally, one could also impose the perturbativity of the power series in Eq. (4) as constraint and define µ nat as its breakdown scale if it is smaller than µ st . This, however, would not alter the maximal µ nat , nor would it constrain the 2HDM parameters stronger than the conventional constraints. It would leave us with the question of what happens beyond the breakdown of perturbative naturalness already at a few TeV.
To put it in a nutshell: softening the Higgs mass hierarchy problem is very difficult in the context of a perturbative 2HDM and can be achieved only for very low cut-off scales µ nat . Nevertheless, this is an improvement of one order of magnitude as compared to the SM hierarchy problem and might hint at a more complete model beyond the 2HDM at TeV scales.

VI. CONCLUSIONS
We obtain the two-loop renormalization group equations for all four Z 2 symmetric types of the 2HDM using PyR@TE and show that in general, two-loop corrections to the leading order one-loop expressions should not be neglected. We then apply these equations to improve the predictions of renormalization group evolution of the coupling parameters, putting a special emphasis on the quartic couplings λ i , which are usually prone to run into non-perturbative regions. The relative distance between the LO and NLO curves of the λ i can be as large as 45%. The quadratic couplings m 2 11 and m 2 22 from the Higgs potential can vary by an order of magnitude between the electroweak scale and the perturbativity cut-off, while m 2 12 is in general found to be rather stable under RG evolution. We do not observe any fixed point behaviour in the regions with a stable Higgs potential.
Imposing positivity and perturbativity bounds at all scales and stability of the vacuum at the electroweak scale, the magnitudes of the λ i at µ ew which give a stable Higgs potential up to the Planck scale are found to be typically below 1; we also find lower limits of 0.15 for λ 2 and of 1.0 for tan β. We have checked that these results are the same in all types. In type II and type Y we additionally get an upper limit of 60 on tan β with stability up to µ Pl . Moreover, we address the question of which upper limit for the eigenvalues of the tree-level Φ i Φ j → Φ i Φ j scattering matrix is appropriate and show that as soon as at least one of the eigenvalues exceeds 2π in magnitude, the maximal scale up to which the Higgs potential can be stable is 5 · 10 6 GeV. It even reduces to 10 TeV in all types if we assume 1 < tan β < 60. Imposing stability up to µ Pl leads to an upper limit of 2.5 on the magnitude of the eigenvalues.
Including latest results from the LHC as well as all other relevant experimental data, we show the result of our fits for all four types. We observe that deviations from the alignment limit strongly depend on the value of tan β; the maximal deviation of β − α from π/2 is 0.45, 0.38, 0.35 and 0.38 in type I, II, X and Y, respectively. (This corresponds to deviations of sin(β − α) from 1 of at most 0.098, 0.070, 0.059 and 0.071.) Taking the stability constraint up to µ Pl , the bounds on β − α become even stronger and allow for deviations from π/2 of at most 0.40, 0.038, 0.26 and 0.041 in the respective types. The searches for heavy neutral Higgs particles exclude a light charged Higgs boson for large tan β in type II and for very large tan β in type Y. In the m H/A -(β − α) planes it is visible that deviations from the alignment limit by more than 0.05π are possible only for narrow regions featuring relatively light m H and m A in the types II, X and Y and excluded for masses larger than 600 GeV in all types.
We finally discuss whether a reliable statement on the seemingly fine-tuned Higgs mass m h can be made in the context of a 2HDM and whether its hierarchy problem can be solved at least partially. Restricting higher order corrections to the perturbative regime, we observe a maximal naturalness cut-off at 5 TeV. Our conclusion is that within a perturbative framework a natural cancellation of quadratic divergencies cannot be implemented into a Two-Higgs-Doublet model beyond O(TeV) scales.

ACKNOWLEDGMENTS
We thank N. Craig, A. Kundu, and A. de la Puente for fruitful discussions, U. Nierste for helpful advice and proofreading, and the CKMfitter group for allowing us to use their statistical analysis framework. We are very grateful to F. Lyonnet for support concerning the PyR@TE code and to the MadGraph and FeynRules authors for useful input, and we want to thank H. Lacker for providing computational power. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 279972. DC would like to thank the Centro de Ciencias de Benasque Pedro Pascual for its hospitality during the initial stages of this work. He has also received partial support from the Munich Institute for Astro-and Particle Physics (MIAPP) of the DFG cluster of excellence "Origin and Structure of the Universe".

APPENDIX
Here we list the renormalization group equations for the 2HDM with soft Z 2 breaking which we obtained with the PyR@TE code [76]. For any coupling L the complete β functions at NLO can be split into leading and next-toleading order contributions and further divided into bosonic and fermionic parts as follows: Except for the Yukawa RGE the bosonic part does not involve fermionic couplings and is type independent, while the fermionic part in general depends on the type of Z 2 symmetry. If the expressions for the latter differ for the different types, we will replace the index f by the type label I, II, X or Y, respectively.