Renormalization of the C2HDM with FeynMaster 2

We present the one-loop electroweak renormalization of the CP-violating 2-Higgs-Doublet Model with softly broken $\mathbb{Z}_2$ symmetry (C2HDM). The existence of CP violation in the scalar sector of the model leads to a quite unique process of renormalization, since it requires the introduction of several non-physical parameters. The C2HDM will thus have more independent counterterms than independent renormalized parameters. As a consequence, different combinations of counterterms can be taken as independent for the same set of independent renormalized parameters. We compare the behaviour of selected combinations in specific NLO processes, which are assured to be gauge independent via a simple prescription. FeynMaster 2 is used to derive the Feynman rules, counterterms and one-loop processes in a simultaneously automatic and flexible way. This illustrates its use as an ideal tool to renormalize models such as the C2HDM and investigate them at NLO.


Introduction
The discovery of a scalar particle at the Large Hadron Collider (LHC) in 2012 [1,2] marked the beginning of the unveiling of the scalar sector of particle physics. Since then, several properties of the newfound particle have been ascertained and they are so far compatible with those of the Higgs boson predicted by the Standard Model (SM) [3][4][5]. But it is possible that such particle belongs to a theory with an extended scalar sector. Actually, since the SM is unable to explain open problems like dark matter or baryogenesis, New Physics beyond the SM (BSM) is expected. One of the most promising directions of BSM physics is precisely the possibility of extra scalars. In recent years, the LHC has been actively searching for extra Higgs bosons [6,7], while plans for future colliders clearly aim at a better understanding of the scalar sector [8,9]. Among the various classes of models with an extended scalar sector, the simplest one that can provide a new source of CP violation required by the three Sakharov criteria for baryogenesis [10] is the 2-Higgs-Doublet Model (2HDM) [11] (for reviews of the 2HDM, cf. refs. [12][13][14]). One of the simplest versions of 2HDM with a CP-violating scalar sector is the so-called complex 2HDM (C2HDM). This model, characterized by a softly broken Z 2 symmetry and complex parameters in the scalar potential, was originally discussed in ref. [15], having been later developed and used by many authors [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33]. It describes a rich phenomenology, featuring five physical scalars: a charged Higgs pair and three neutral scalars which are in general a mixture of a CP-even and a CP-odd component. The model becomes especially relevant in face of the very recent suggestion [34] that the so-called "real 2HDM"-a widely known variant of 2HDM where CP conservation is enforced in the scalar sector-is theoretically inconsistent. It turns out that the C2HDM is the simplest consistent version of such variant which can account for CP violation.
In the last few months, the C2HDM has aroused particular interest; for example, a 2loop renormalization group evolution of the model was discussed in ref. [35]; ref. [36] used the model to study phase transition dynamics and gravitational wave signals; ref. [37] investigated a basis-independent treatment of the model; studies on electric dipole moments (EDM) of the C2HDM have been put forward in refs. [38,39]; ref. [40] discussed the impact of the discovery of a new scalar particle on the parameter space of the model; ref. [41] considered CP-violating gauge-scalar interactions; ref. [42] investigated Higgs alignment and signatures of CP violation in the context of the model; and ref. [43] derived constraints on the model from the phenomenology of light scalars.
However, precise predictions in the C2HDM have not yet been systematically considered. Precise predictions are the key to probe the SM and BSM models, since they are indispensable for a) a sound interpretation of the observed results, b) a correct determination of the parameter space of models and c) a proper distinction between different BSM models. Precise predictions require the inclusion of one-loop electroweak corrections, which in turn require the one-loop electroweak renormalization of the model. 1 The literature devoted to the renormalization of models with an extended Higgs sector is vast (for the SM, reviews can be found in refs. [44][45][46][47][48]); the "real 2HDM" has been subject to several studies [49][50][51][52][53][54][55][56]; the next-to-minimal 2HDM has been considered in ref. [57], and the scalar sector of a variant of the Minimal Supersymmetric extension of the SM with complex parameters (cMSSM) was explored for instance in refs. [58][59][60][61][62][63][64][65][66][67][68][69].
In this paper, we present the renormalization of the C2HDM. To our knowledge, it is the first time that the renormalization of a model with explicit CP violation in the scalar sector is put forward. This leads to a unique process of renormalization, as one is forced to introduce several parameters that can be rephased away in the context of the renormalized parameters, but must be considered anyway in order to assure the generation of all the necessary counterterms. As a consequence, there will be more independent counterterms than independent renormalized parameters. In particular, different combinations of counterterms can be taken as independent for the same set of independent renormalized parameters.
A crucial element in all of this is FeynMaster [70], a multi-tasking software for particle physics studies. Combining FeynRules [71,72], QGRAF [73] and FeynCalc [74][75][76], FeynMaster is able to perform the totality of the following list of tasks in a consistent, automatic and flexible way: A new version of the program-FeynMaster 2-was recently made publicly available at https://porthos.tecnico.ulisboa.pt/FeynMaster/, containing several improvements over the first version. FeynMaster 2 turns out to be an ideal tool for building and investigating models, and especially for renormalizing them. In this paper, we apply it to the renormalization of the C2HDM and illustrate the advantages of doing so. The paper is organized as follows. In section 2, we describe the C2HDM in such a way that we aim at the one-loop renormalization of the model, presenting different combinations of independent parameters. Section 3 introduces the treatment of the theory when considered at up to one-loop level, clarifying how the theory can contain counterterms for quantities that do not show up in the set of renormalized quantities. Sections 4 and 5 proceed with the aforementioned treatment: the former describes the selection of the true vacuum expectation value (vev), whereas the latter is devoted to the counterterms of the C2HDM. Finally, we present our results in section 6. After the conclusions in section 7, some appendices are included. Appendix A provides a simple description of the on-shell subtraction scheme and investigates the scenarios where the mass of a particle is a dependent parameter. Appendix B is devoted to CP violation in fermionic 2-point functions and its influence on counterterms. Then, in appendix C, we show how counterterms can be fixed through symmetry relations. Finally, we dedicate appendix D to FeynMaster 2 and to its application to the C2HDM.

The C2HDM
The complete Lagrangian of the C2HDM can be written as a sum of the partial Lagrangians for the different sectors of the theory: L C2HDM = L Gauge + L Fermion + L Higgs + L Yukawa + L GF + L Ghost . (2.1) The terms L Gauge and L Fermion are just those of the SM, while L GF and L Ghost can be easily derived from the SM ones (cf. e.g. ref. [77]). Here, we study in detail the Higgs and the Yukawa partial Lagrangians. The former can be split into kinetic terms and potential, We start by studying the potential V in section 2.1, then the scalar kinetic terms in section 2.2, and finally the Yukawa Lagrangian in section 2.3. We shall use a general parameterization, which implies that we will introduce several parameters that are absent in the usual tree-level description of the model (cf. e.g. ref. [31]). 2 As it will become clear, this is required in order to assure the one-loop renormalization of the theory.

The potential
Assuming a softly broken Z 2 symmetry, under which Φ 1 → Φ 1 , Φ 2 → −Φ 2 , the potential can be written as [31] V = m 2 3) The hermiticity of the Lagrangian obliges all parameters to be real except m 2 12 and λ 5 , which are in general complex. In what follows, we define: 12R ≡ Re m 2 12 , m 2 12I ≡ Im m 2 12 , λ 5R ≡ Re (λ 5 ) , λ 5I ≡ Im (λ 5 ) . (2.4) After spontaneous symmetry breaking (SSB), each of the Higgs doublets acquires a vev. These are in general complex, and in general relatively complex when compared to the scalar fields. We thus parameterize the Higgs doublets as: , (2.5) where v i are real parameters, ζ i (real) phases, φ + i complex fields, and ρ i and η i real fields (i = 1, 2). The writing of the vevs with modulus and phase allows us to define the real parameters v and β such that with c β = cos β, s β = sin β.

Charged bilinear terms
The mass matrix of the charged scalars, defined by is hermitian, so that one needs in general a unitary matrix to diagonalize it. We then define the unitary matrix X such that: 4 where ζ a and ζ b are (real) phases, χ a (real) angle, and the fields S c (with S c = (G + H + ) T ) are the charged states in the mass basis, with G + corresponding to the (massless) charged would-be Goldstone boson. By definition, X is such that: (2.11) which in turn implies M 2 c = X † D 2 c X. (2.12) It turns out that there are two non-trivial relations between the elements of M 2 c , namely, which must also be verified for the right-hand side of eq. 2.12, thus leading to two identities: (2.14)

Neutral bilinear terms
The mass matrix of the neutral scalars, defined by is symmetric, which means that one needs an orthogonal matrix to diagonalize it. Hence, we define the orthogonal matrix Q such that: where S n (with S n = (h 1 h 2 h 3 G 0 ) T ) are the neutral states in the mass basis, with G 0 corresponding to the (massless) neutral would-be Goldstone boson. Concerning Q, since it is a 4 × 4 orthogonal matrix, one in general needs 6 angles to parameterize it, which we take to be α 0 , α 1 , α 2 , α 3 , α 4 , α 5 , such that: . By definition, Q is such that: where m i represents the mass of h i , such that m 1 < m 2 < m 3 . In order to find relations in this case, it is convenient to rewrite the doublets as: , (2.20) where the neutral fields with primes are related to the original ones (in eq. 2.5) via: Then, where we defined: There are simple relations between the elements of M 2 n : On the other hand, eqs. 2.19 and 2.22 imply: (2.25) which means that the elements of the right-hand side of this equation must obey the same nine relations of eq. 2.24. We thus have nine relations, involving in general the parameters: It turns out that only four relations are independent, which means that one can only fix four of the parameters of eq. 2.26. Obviously, there are many choices (or combinations) for the possible set of dependent parameters. In this work, we consider the combinations C i (i = 1, 2, 3, 4) described in table 1. Some notes are in order here.
Combination Dependent parameters Table 1: The dependent parameters associated to the different combinations C i . First, m 2 3 is chosen as a dependent parameter in all combinations; this is not accidental, and can be justified as follows. When considering the model up to one-loop level, we shall find that there can be a clear relation between the renormalized parameters, on the one hand, and the parameters of the usual tree-level description, on the other. For the latter, we use ref. [31] as a guiding reference; as it turns out, m 2 3 is taken as a dependent parameter in that reference. Hence, in order for this parameter to show up as dependent in the renormalized parameters of the model considered up to one-loop level, it must necessarily be taken as dependent at tree-level. This also allows to explain the four combinations chosen in Table 1. Indeed, if we wish to obtain the aforementioned clear relation, the parameters taken as independent in ref. [31] must also be taken as independent here; in particular, α 1 , α 2 , α 3 , β, m 2 1 , m 2 2 must always be taken as independent. As a consequence, besides m 2 3 (which is always dependent), the only parameters available to be dependent are α 0 , α 4 , α 5 , ζ 1 and ζ 2 . We decide to take ζ 2 as independent, which we will justify below. Therefore, since there are only four dependence relations, there are only four possible combinations of dependent counterterms-precisely those of Table 1. Finally, the system of equations at stake is non-linear, and a rather complex one. We shall only solve it when we consider the model up to one-loop level. 5 For now, the dependent parameters must be understood simply as functions of the independent parameters.

Combined sectors
We have already derived some dependence conditions that resulted from relations we found within the squared mass matrices. Specifically, we found relations between the elements of M 2 c , as well as relations between the elements of M 2 n . The former lead to the two conditions in eq. 2.14, whereas the latter to the dependence relations implied in table 1. But there are still relations we have not yet used, which come from two equalities: that of eq. 2.12 and that of eq. 2.25. These lead to six independent relations, that we may use to rewrite the six λ i (i = 1, 2, 3, 4, 5R, 5I) in terms of other parameters. 6

Scalar kinetic sector
We now consider the first term of the right-hand side of eq. 2.2, where the covariant derivative is defined by (2.28) Here, g 1 and g 2 are the gauge couplings of the U(1) and SU(2) gauge groups, respectively, with B µ and W a µ (a = 1, 2, 3) being the corresponding gauge fields, and Y and τ a the corresponding group generators, respectively (we are following the conventions of ref. [77] with all η's positive). 7 After SSB, the physical gauge fields W ± µ , A µ and Z µ are obtained from the original fields through the relations: with s w = sin θ w , c w = cos θ w , where θ w is the weak mixing parameter, which is such that: (2.30) Expanding the bilinear terms in the gauge fields in eq. 2.27, one easily finds that the photon A µ is massless, whereas the squared masses of the W + µ and Z µ bosons are, respectively, Finally, defining the electric charge as e = g 1 g 2 we can take v, g 1 and g 2 as dependent parameters, which are then written as: (2.33)

Yukawa sector
The Yukawa Lagrangian in general leads to flavour-changing neutral currents at tree-level. A simple way to avoid this is to assure that each right-handed fermionic singlet couples to only one Higgs doublet. This in turn can be accomplished if the Z 2 symmetry is extended to the fermion fields, such that: Here, a, b, c, d, e are general powers, q L = (p L n L ) T and L L = (ν L l L ) T are the quark and lepton left-handed SU(2) doublets, respectively, and p R , n R and l R are the up-type quark, down-type quark and lepton right-handed singlets, respectively. There are four different combinations of the powers a to e, each combination corresponding to a different type of C2HDM, as can be seen in Table 2. In this article, we restrict ourselves to the Type II in which case we have, in the Type II C2HDM: (2.36) where Y d , Y u and Y l are the Yukawa matrices for the down-type quarks, up-type quarks and leptons, respectively. 9 We can, however, rewrite this equation in a more meaningful way. We start by noting that the quarks are rotated to the mass basis through unitary transformations in such a way that the interaction with the vevs of the Higgs doublets generates the mass terms, that is, In this equation, the SU(2) product is shown explicitly, but the sum over fermion generations is left implicit by the matrix notation. Had we written it explicitly, the quantitiespL,nL, Y d , nR in the first term on the right-hand side of eq. 2.36, for example, would have beenpL i ,nL i , Y ij d , nR j , respectively.
On the other hand, the mass terms must also obey: we find: 10 .
which allows us to write eq. 2.36 as Finally, note that the fermion masses m f are in general complex parameters, since the Yukawa matrices are general complex matrices. 11 Usually, one performs a chiral rotation (that is, a rotation of the Weyl spinors) in order to render the masses real. 12 But the circumstance that the masses are in general complex means that, when the theory is considered up to one-loop level, the counterterms for the masses will also be in general complex. And although this is an irrelevant detail in most models, it becomes most relevant whenever there is CP violation in fermionic 2-point functions, as in the present model at one-loop level. We discuss in detail complex mass counterterms-as well as their relation with CP violation in fermionic 2-point functions-in appendix B. For now, we assume m f to be in general complex.

Parameters
The relations introduced in the previous sections allow us to replace the original sets of parameters in the Yukawa and Higgs sectors, respectively given by by new sets of parameters: 10 We assume there are no right-handed neutrinos, which implies we can take Y l to be diagonal without loss of generality. 11 f represents the physical down-type quarks, up-type quarks and leptons, whose masses are respectively contained in M d , Mu and M l . 12 Cf. e.g. section 29.3.2 of ref. [78].
where we define: We thus see that, in the Higgs sector, there are four possibilities, {p H C i } (i = 1, 2, 3, 4), respectively associated to the combinations C i introduced in section 2.1.3. For a given C i , the parameters in the new total set (including parameters from the Yukawa and Higgs sectors) are all independent, and are more convenient than the ones in eqs. 2.43. One can then express the first four terms of the right-hand side of eq. 2.1 in terms of the parameters of eqs. 2.44.

Going up to one-loop level
So far, we have been discussing the theory at tree-level, i.e. in lowest order (LO) in perturbation theory. When we include the next order-that is, when we consider the theory up to one-loop level, or at the next-to-leading order (NLO)-, we must start by modifying the notation. 13 The theory we have been considering, indeed, is to be identified with a bare theory, and the parameters and fields therein contained as bare parameters and fields. Moreover, we make sure to select the true vev of the theory up to one-loop level; such selection is discussed in section 4. The bare parameters and fields are identified with an index "(0)" and they can renormalized.
To do so, one starts by splitting them into renormalized quantities and their corresponding counterterms; for example, for a generic bare parameter p (0) and a generic bare field ψ (0) , Here, p represents the renormalized parameter and δp the corresponding counterterm, and ψ the renormalized field and δZ ψ the corresponding counterterm. 14 The renormalization is 13 The notion 'NLO' does not refer to one-loop level only, but to both tree-level and one-loop level. For the sake of clarity, we will often use the notion 'up to one-loop' instead. Both notions are equivalent and shall be used indifferently in what follows. What we identify here as up-to-one-loop theory is then part (namely, the part which includes the tree-level and one-loop level) of what is usually identified as the effective theory (cf. e.g. ref. [47]). 14 The renormalized quantities in the context of the theory up to one-loop level (which have no index whatsoever) should not be confused with the bare quantities in the context of the theory at tree-level (which, in that context, did not have the index "(0)"). Note also that δp is dubbed counterterm for the parameter (or parameter counterterm) and δZ ψ is dubbed counterterm for the field (or field counterterm).
completed by the calculation of the counterterms, which is performed in section 5. 15 Now, in order to assure that all S-matrix elements are ultraviolet (UV) finite, one must renormalize the independent parameters. The first step is choosing an independent set of parameters. 16 As independent set of parameters, besides those of eq. 2.44a we will consider the four possibilities of independent sets of eq. 2.44b. All these parameters are taken as bare parameters in the context of the up-to-one-loop theory, and thus read: (3.2b) Afterwards, and as suggested above, these bare parameters are separated into renormalized parameter and parameter counterterm, and the latter must be calculated.
On the other hand, if one intends to obtain finite Green's functions (GFs) and not only finite S-matrix elements, fields must be renormalized as well. This can be done in the symmetric form of the theory [45,81]-i.e. before the SSB of the symmetry-, in which case it is enough to consider one field counterterm for each multiplet [47]. As an alternative, one can renormalize the fields in the mass basis [44,46,47]. Here, counterterms for the mixing of fields (mixed fields counterterms) must in general be considered whenever fields mix at loop-level.
In the present section, we start by identifying the bare quantities of the theory with the sum of renormalized ones and counterterms (we do not do this for gauge parameters nor ghost fields, since it would be irrelevant for the calculation of S-matrix elements [47]). 17 Then, in section 3.2, we review the need to consider a general parameterization, and show that the renormalized parameters can be described by a simpler parameterization.

Expansion of the bare quantities
In eq. 3.1 above, we considered generic expansions. We now need to apply them to the specific content of the C2HDM. Concerning the parameters, we include not only those of eqs. 3.2, but also some dependent parameters (for convenience). Some renormalized 15 In general, the LSZ factors (present in the LSZ reduction formula) must also be calculated in order to obtain correct S-matrix elements [47,79]. However, they become trivial when fields are renormalized in the on-shell subtraction scheme, so that they can be ignored. For details, cf. ref. [80]. 16 Dependent parameters need not be renormalized, although they can be renormalized for convenience; for details, see ref. [80]. 17 Actually, in linear gauges (as the Feynman gauge), it can be shown that one can simply restrain from applying any renormalization transformation whatsoever to the Gauge Fixing Lagrangian [82][83][84]. That is, all the parameters and fields therein are assumed to be already renormalized.
parameters will be set equal to non-obvious quantities, which will be justified below. As for the fields, we choose those in the mass basis. We start with the gauge boson masses and fields. For the masses, we have: while for the fields: In the fermions, the masses obey in such a way that m f,i are real, and we define: 18 The fields obey: Considering now the scalar sector, we have for the masses: We renormalize m 2 3(0) for convenience: The renormalized squared mass in this case explicitly includes an index R (from 'renormalized'). 19 The reason why this is necessary in this case is that, since m 3(0) is a dependent parameter, it is fixed by other parameters. As a consequence, when the bare quantities are identified with a renormalized quantity plus a counterterm, the renormalized squared mass m 2 3R will also be fixed by other parameters, which means that it cannot be set equal to the pole squared mass-which we identify in the following as m 2 3P . Moreover, δm 2 3 will be a dependent counterterm. All this implies that h 3 plays a special role in theory, which is described detail in appendix A. For the scalar fields, The mixing parameters are such that: 20 The electric charge and the CKM matrix obey: As for the remaining parameters that are taken as independent (ζ 1(0) will be independent in the combination C 4 ), we have: For convenience, we also renormalize several dependent parameters: 21 (3.14)

General vs. simple parameterization
When studying the theory at tree-level, we considered the most general parameterization of the Higgs doublets, as well as of the diagonalization of the scalar states-recall section 2.1. This general information was required in order to have general counterterms (which are necessary to absorb the UV divergences at one-loop). But it is not required to describe the renormalized parameters. Indeed, if the only reason we considered such a general parameterization was to account for the generality of the counterterms, there is no need to apply such general description to the renormalized parameters as well.
To better understand the claim we are making, consider for example the phase ζ a . At tree-level, we had the freedom to rephase the fields G + and H + to absorb ζ a . Had we used that freedom, though, we would have run into problems; for when we were to consider the theory up to one-loop level, we would have no counterterm δζ a (as there would have been no bare parameter ζ a(0) to start with). It turns out that such counterterm is necessary 20 Recall that α 0(0) , α 4(0) and α 5(0) will be independent according to the combination Ci chosen. 21 In the combinations Ci for which the parameters α 0(0) , α 4(0) , α 5(0) and ζ 1(0) are dependent, we also renormalize them for convenience (according to eqs. 3.11c and 3.13). We do not renormalize χ, ζ b , g1, g2 or v. This means that, before renormalization, they must be replaced by the corresponding expressions (cf. eqs. 2.14 and 2.33), and the latter must be renormalized. More details can be found in ref. [80].
for the absorption of the divergences of the theory, so that it cannot be discarded. That is to say, we could not use the freedom to rephase ζ a away at tree-level. But once in the up-to-one-loop theory, and given that δζ a was already generated, there is no longer any reason not to use the rephasing freedom. We can thus rephase the renormalized parameter ζ a away in the up-to-one-loop theory; and when we do so, the counterterm δζ a (that was meanwhile generated because we did not rephase ζ a away at tree-level) will not vanish. To see this explicitly, we may consider a simplified version of eq. 2.10, where we take ζ b = χ = 0 (which are irrelevant for the argument): It is obvious that, in this tree-level relation, one is free to rephase ζ a away through G + → e −iζa G + , H + → e iζa H + . Yet, as we want to generate a counterterm δζ a , we do not use that rephasing freedom. So, when considering the theory up to one-loop level, we take the bare version of eq. 3.15; then, expanding the charged scalar states in the mass basis and ζ a(0) into renormalized quantities and counterterms, according to eq. 3.10a and we obtain, to first order, Since δζ a was already generated, the renormalized parameter ζ a can be rephased away This simple example illustrates several relevant points. First, it shows that one must start with a general description including the tree-level parameter ζ a , for otherwise the counterterm δζ a will not be generated. Second, after δζ a was generated, the renormalized parameter ζ a can be rephased away, and in such a way that δζ a does not vanish. Third, instead of starting by splitting ζ a(0) into a non-zero ζ a plus a counterterm (as in eq. 3.16), and afterwards rephasing ζ a away, we can assume beforehand that we will work in a basis for the renormalized parameters where ζ a = 0, which lets us equate ζ a(0) immediately with the counterterm only. This is what we did in eq. 3.11a above. Finally, what we obtained for the renormalized parameters looks just the same as that which we would have obtained at tree-level should we not have been aiming at renormalizing the theory. As a matter of fact, if we were not concerned with the generation of counterterms, but rather in a mere tree-level description, we would not have needed ζ a , so that we could have rephased it away in eq. 3.15. In that case, we would have obtained the same which we obtained in eq. 3.18 for the renormalized terms (i.e. excluding counterterms). 22 This also holds for the other parameters that were not rephased away for the sake of renormalization only. Therefore, the basis for the renormalized parameters (i.e. the basis that we can and shall use for the renormalized parameters) is considerably simpler than the general one used in section 2.1, and is just the same as the one usually employed when considering the theory solely at tree-level-as in ref. [31]. Comparing this reference with the general parameterization used in section 2.1, we conclude that, for the renormalized parameters in the up-to-one-loop theory, all the fermion masses are real, and 23 As a consequence, the relations one obtains for the renormalized parameters in the context of the up-to-one-loop theory are precisely those that were derived when considering the theory solely at tree-level. We checked this explicitly for different expressions: the minimum equations [31], the relations for the λ i [26], as well as the relation for m 2 3R [31]. The latter reads: . (3.20) Here, the matrix R is given by: and is related to the matrix Q in a simple way. Indeed, given eq. 3.19, Q becomes: The original parameterization of Q in eqs. 2.17 and 2.18 thus allows a simple connection between the matrix Q for the renormalized parameters and the matrix R. The matrix X becomes orthogonal, and is simply: Finally, the total independent renormalized parameters become: 22 Care should be taken with the notion 'the same', because in one case what is at stake are renormalized quantities in the context of the up-to-one-theory, whereas in the other case there are only tree-level quantities in a tree-level theory. Hence, they are formally different. Nonetheless, the physical description in both cases is equivalent. 23 The reason for identifying the renormalized phases ζ1R and ζ2R with the index R is given in note 26.
These parameters are precisely those which are usually taken as independent in the theory considered solely at tree-level [31]. Moreover, the set of renormalized parameters of the Higgs sector, {p H reno. }, is the same whichever the set of bare parameters {p H C i (0) } taken as the independent (recall eqs. 3.2b). Actually, not only it is the same, but it contains less parameters than each of the sets {p H C i (0) }, since the renormalized versions of the last three independent bare parameters of each set were rephased away. On the other hand, for each bare independent parameter, there is an independent counterterm. Thus, whichever the combination of independent bare parameters chosen, there will be more independent counterterms than independent renormalized parameters. In fact, identifying {δp H C i } with the set of independent counterterms of the Higgs sector in the combination C i , we have: Something similar also happens in the Yukawa sector, as the phases of the renormalized fermion masses were rephased away. This means that, while the counterterms δm f are complex, the renormalized masses m f are real. In summary, the presence of more independent counterterms than independent renormalized parameters is a consequence of CP violation, which forces the introduction of several parameters; these are required for the one-loop renormalization of the model, but their renormalized versions can be rephased away.

Selection of the true vev
When considering the theory up to one-loop level, it is convenient to assure that scalar fields have zero vacuum expectation value. Such assurance corresponds to the selection of the true vev, and is discussed in detail in ref. [80]. In the present work, we perform that selection using the Fleischer-Jegerlehner tadpole scheme (FJTS) [54,56,85]. 24 We start by considering the bare version of eq. 2.5: The quantities with a bar correspond to true quantities, in the sense that they guarantee that no proper tadpoles (i.e. 1-point GFs) show up in the up-to-one-loop theory. In the FJTS, they are split as follows: 24 Other relevant studies concerning the vev can be found e.g. in refs. [86,87].
Thus, each quantity with a bar is split into two terms: the bare one and an extra quantity. These extra quantities are identified in what follows as the ∆ quantities and are responsible for the elimination of proper tadpoles (they are introduced with this purpose). To see this, it is convenient to write eq. 4.1 in an alternative formulation: Comparing with eqs. 4.1 and 4.2, we obtain: We neglected terms of order ∆ 2 , since the ∆ quantities are of one-loop order (as we are about to show). This also means that the bare quantities can be replaced by their renormalized versions. 25 Hence, using the fact that we choose a basis for the renormalized parameters such that the renormalized phases ζ 1R and ζ 2R vanish (eq. 3.19), we get: 26 (4.6) Now, when the Lagrangian is expanded with eq. 4.3, the ∆v φn of eq. 4.4 must assure that proper tadpoles vanish. In order for this to happen, the terms with ∆v φ n,j that contribute to the 1-point function of a certain neutral scalar field φ n,i must precisely cancel the oneloop 1-point function in that field. That is, if T φ n,i represents the one-loop tadpole for φ n,i (i.e. the one-loop contribution to the 1-point function in φ n,i ), we must have: Using the fact that every ∆v φ n,j always shows up together with the corresponding field φ n,j (cf. eqs. 4.3 and 4.4), we can write: where the mass matrix M 2 n obeys eq. 2.19. Then, going to the diagonal basis,  which implies: where we used the fact that T G 0 = 0, which is a consequence of the Goldstone theorem. Some comments are in order. First, by obeying eq. 4.7, the ∆ quantities-in the physical basis, ∆v h 1 , ∆v h 2 , ∆v h 3 , ∆v G 0assure that no proper tadpoles (1-point GFs) show up in the theory. In other words, they guarantee that the true vev up to one-loop level is selected. But as they are contained inside the scalar doublets, they will contribute to other GFs of theory. For example, although no 1point GF for h 2 remains in the theory (that is, no proper tadpole for h 2 remains), ∆v h 2 will contribute to some 2-point and 3-point functions; and given eq. 4.10, the one-loop tadpole T h 2 will contribute to those GFs; we dub these one-loop tadpole structures contributing to GFs other than the 1-point GFs broad tadpoles. In this way, ∆v h 2 assures that no proper tadpole for h 2 exists, but introduces broad tadpoles in the theory. Broad tadpoles are thus a consequence of the selection of the true vev and must be considered for consistency. Second, all broad tadpoles can be accounted for by considering all possible one-loop tadpole insertions in all possible GFs (as usual in the FJTS [54]). Indeed, as suggested above, each renormalized scalar field will have the corresponding ∆ quantity added to it: ∆v h 1 + h 1 , ∆v h 2 + h 2 , and so on. This means that, for every term in the Lagrangian with (say) ∆v h 2 , there will be a term with h 2 instead. As a consequence, for every n-point function with a ∆v h 2 contribution, there is a (n + 1)-point function with a h 2 field contribution instead; and with that field, a reducible diagram with a one-loop tadpole can be formed. But it is easy to see that such reducible diagram precisely equals the n-point function with ∆v h 2 : one just needs to consider the definition of the one-loop tadpole for h 2 , and the zero-momentum h 2 propagator, and compare with eq. 4.10. 27 The same argument obviously applies to the remaining broad tadpoles. In sum, although broad tadpoles must be considered for consistency, in the FJTS they can be accounted for by including all possible one-loop tadpole insertions in all possible GFs (cf. also refs. [52,57]). Third, the alternative just proposed (calculating reducible diagrams with one-loop tadpoles instead of tracing down the terms with ∆ showing up in the Lagrangian) will be the one adopted in the following. The reason is simply that, with softwares such as FeynMaster, the calculation of reducible diagrams with one-loop tadpoles is trivial; 28 in contrast, the task of expanding the Lagrangian and tracing down all the terms with ∆ is quite cumbersome; and since the two methods are equivalent, we opt for the former. It follows that all GFs that may receive an insertion of a one-loop tadpole must receive such insertion. In general, then, non-renormalized 2-point functions will get contributions from two different types of terms: where the full lines represent any type of particle and the dashed lines represent h 1 , h 2 , h 3 . 29 The same will happen for non-renormalized 3-point functions, which then include: (4.14)

Counterterms
As mentioned in section 3, the renormalization of the theory is completed by fixing the counterterms. Independent counterterms are a priori not fixed, and must be fixed through independent renormalization conditions. Now, counterterms in general contain a divergent part and a finite part; whichever the renormalization condition used to fix a certain counterterm, its divergent part will always be the same. But the same is not true for the finite part, which in general depends on the prescription used. The different prescriptions that can be used to define the finite part of counterterms are called subtraction schemes. In what follows, most counterterms will be fixed either through the modified minimal subtraction (MS) scheme or through the on-shell subtraction (OSS) scheme. The former is simpler: one selects a process that depends on the counterterm to be calculated; then, the renormalization condition simply states that the terms of the renormalized one-loop amplitude proportional to ∆ are zero. Here, ∆ is 28 For details, cf. appendix D. 29 We draw small black circles on the dashed lines to highlight the presence of a tree-level propagator.
where = 4 − d (d being the dimension) and γ E is the Euler-Mascheroni constant. The counterterm at stake will then be proportional to ∆ only. When MS is used to calculate a parameter counterterm, the corresponding renormalized parameter will have no clear physical meaning. In contrast, OSS is defined precisely so that all renormalized parameters have a direct physical meaning [47]; the renormalization conditions in this case are less direct, though, and are discussed in detail in appendix A. In what follows, we adopt the conventions therein defined. In section 5.1, we fix all of the counterterms of the theory, proceeding by sectors. Several counterterms are fixed through usual renormalization conditions (cf. e.g. ref. [47]); in those cases, we omit the derivation, and present just the results. 30 The remaining cases will be considered in detail. Finally, in section 5.2, we discuss the gauge dependence of the counterterms.

Gauge masses and fields
The counterterms for the gauge sector are determined through OSS. The results are as follows: for the charged sector [47,80], while for the neutral sector [47,80], (5.3b)

Fermion masses and fields
The counterterms for fermion masses and fields are also fixed through OSS conditionswhich, in this case, leave some freedom. We show in detail in appendix B how to derive the counterterms. The results are: 30 A detailed derivation of all the counterterms can be found in ref. [80].

Scalar masses and fields
The counterterms relative to the charged scalar sector are easily defined in the OSS scheme; the results are [57,80]: The case of the neutral fields is more subtle, and is considered in detail in appendix A. The independent counterterms are: 31 and the pole squared mass m 2 3P is given by:

Electric charge
The counterterm for the electric charge is also fixed through OSS; the result is [57,80]: In the case Sn = h3, m 2 Sn represents m 2 3R .

CKM matrix
For the CKM matrix, we follow the prescription originally presented in ref. [88]. We fix the expression in the Feynman gauge, which will be justified in section 5.2 below. We then have: where ξ = 1 represents the Feynman gauge.

Mixing parameters
The renormalization of mixing parameters has been recurrently discussed in recent years [50][51][52][53][54][55][56][57], especially in ref. [56]. 32 In the latter, several methods for fixing counterterms for mixing parameters are proposed and analysed. In the present work, we consider one of those methods, which involves fixing the independent counterterms for mixing parameters using symmetry relations. This method is simple and leads to well-behaved results, as shall be discussed below. A description of the procedure, as well as the derivation of the expressions for the counterterms for mixing parameters of the C2HDM, can be found in appendix C.
In this section, we present the results therein derived. For the charged sector, For the neutral sector, the set of independent counterterms for mixing parameters depends on the combination C i (recall eq. 3.26). There are three counterterms which are independent in all the combinations: Then, the combinations C 1 , C 2 and C 3 have δα 5 , δα 4 and δα 0 as independent, respectively, given by: 33 32 Mixing parameters usually reduce to mixing angles. Yet, in the present model (and due to CP violation), there are not only mixing angles (like β), but also mixing phases (like ζa). Hence, we promote the designation 'mixing angles' to 'mixing parameters'. 33 The notation C i = clarifies which combination Ci the expression is applicable to. For example, eq. 5.12a is only valid if C1 is chosen; in other combinations, δα5 will be a dependent counterterm, so that eq. 5.12a is not valid. Finally, note that δα5, δα4 and δα0 are all dependent counterterms in the combination C4.

Remaining parameters
By now, all the field counterterms are determined; as for the independent parameter counterterms, considering eq. 3.26 we see that we still need to fix δµ 2 and δζ 2 (both required for all combinations C i ) and, in the specific case of the combination C 4 , also δζ 1 . We calculate the latter and δµ 2 in the MS scheme. As described above, this requires selecting a process where the counterterm to be fixed intervenes. For δµ 2 , we choose h 3 → H + H − . Here, the total counterterm is such that: Then, by requiring the renormalized 3-point function to be such that the terms proportional to ∆ are zero, 14) we fix δµ 2 . Concerning δζ 1 in C 4 , we choose h 3 →τ τ . Knowing that: where the ellipses represent terms proportional to γ 5 , we fix δζ 1 in C 4 by imposing: As for δζ 2 , we checked that its divergent part is related to that of δζ 1 via: 34 Then, similarly to what can be done with symmetry relations (cf. appendix C), we fix δζ 2 as a whole by requiring that eq. 5.17 also holds for the finite parts. That is, we fix δζ 2 by imposing: 35 Note that this condition applies to all four combinations C i . Accordingly, δζ 2 will have a general finite part in the first three combinations, whereas in C 4 it will be proportional to ∆ only.

Gauge dependence
Physical observables cannot depend on the gauge chosen, which implies that S-matrix elements must be gauge independent. In a renormalized up-to-one-loop process, both the (non-renormalized) one-loop diagrams and the counterterms are in general gauge dependent. Now, it can be shown that, if parameter counterterms are gauge independent, the gauge dependences of the one-loop diagrams and those of the field counterterms precisely cancel, thus rendering the S-matrix elements gauge independent [79]. Conversely, gauge dependent parameter counterterms in general lead to gauge dependent S-matrix elements. Motivated by this, several authors have insisted in recent years on the importance on finding gauge independent parameter counterterms [52,54,57,89]. However, as noted in ref. [56], the circumstance that parameter counterterms are gauge dependent does not by itself ruin the gauge independence of S-matrix elements. In fact, all that is required for S-matrix elements to be gauge independent is that parameter counterterms are considered gauge independent. In order to assure the gauge independence of S-matrix elements, indeed, it does not really matter whether the parameter counterterms are truly gauge independent or not; all what is required is that they are considered gauge independent, i.e. that they do not change when the gauge is changed [56,79]. 36 But this can be easily implemented by calculating a gauge dependent counterterm in a particular gauge. When calculating observables in other gauges, one must stick to the value of the parameter counterterm in the gauge it was calculated in. Then, the gauge dependence of field counterterms cancels the gauge dependence of the one-loop diagrams, and S-matrix elements will be gauge independent. In face of this, there does not seem to be a clear advantage of true gauge independence in parameter counterterms over considered gauge independence. In both scenarios, S-matrix elements end up being gauge-independent, and both lead to correct relations between physical observables. 37 34 The divergent part of δζ2 can be easily calculated by requiring e.g. the process h3 →ūu to be finite. 35 Eq. 5.17 thus provides a simple prescription to fix δζ2. This simplicity is the reason why we chose this counterterm as independent in all the four combinations discussed above. 36 We thank Ansgar Denner for clarifications on this point. 37 It goes without saying that, when using considered gauge independent counterterms, one has to make sure that no inconsistencies are introduced (i.e. one must make sure not to change the gauge in the Now, this procedure is obviously irrelevant when parameter counterterms are truly gauge independent. In the FJTS, this happens for parameter counterterms fixed in either the OSS or the MS schemes. 38 However, the method described in the previous paragraph is very convenient for the remaining independent parameter counterterms, which in this model (and according to our choices) consist of those for the CKM matrix elements, those for the independent mixing parameters, and δζ 2 (sections 5.1.5, 5.1.6 and 5.1.7 above, respectively). We thus applied the method in those cases, following the suggestion proposed in ref. [56] of selecting the Feynman gauge as the gauge in which the counterterms are calculated, which avoids artificially large parameters.

Results
We now study the behaviour of the different combinations C i introduced above. This agenda is quite novel in studies of renormalization, due to the mere existence of combinations. In fact, there is usually no such thing as different combinations of independent parameters, as the number of independent renormalized parameters usually equals the number of independent counterterms. The studies of the renormalization of models thus tend to be focused either on the choice of different subtraction schemes for the counterterms, or on the choice of different renormalization scales. 39 In our case, the fact that there are more independent counterterms than independent renormalized parameters means that, for the same set of the latter, there can be different combinations of the former. Our goal, then, is to compare those combinations, ascertaining how the choice of one or the other combination (i.e. the choice of one or the other set of independent counterterms) affects the predictions of the theory at NLO. 40 To do that, we consider three specific NLO processes: the decays of h 2 to ZZ, h 1 Z and h 1 h 1 .
In the present section, we start by describing both the computational tools and the simulation procedure used to study these decays. Then, we discuss how the decays are influenced by the different combinations C i . Finally, after presenting the expressions for the one-loop decay widths, we show numerical results that allow to compare the different combinations.

Computational tools and simulation procedure
The generation of Feynman rules for both the tree-level and the counterterm interactions, the drawing of the Feynman rules and diagrams, and the calculation of one-loop amplitudes, counterterms and one-loop decay widths were all performed with FeynMaster 2.
gauge dependent parameter counterterms). In this sense, truly gauge independent counterterms may be preferable. It has been claimed that true gauge independence is a desirable property, as it allows a "more physical interpretation" of non-observable parameters [90]. However, a non-observable parameter is always non-physical anyway; the fact that it is truly gauge independent does not render it any more physical (in whatever physically meaningful way). 38 For details, cf. ref. [80]. 39 Should there be counterterms fixed in a subtraction scheme that introduces a renormalization scale dependence. 40 In this paper, we restrict ourselves to the four combinations Ci implied in eq. 3.26 and, for a certain combination, to the subtraction schemes proposed in section 5.
These different tasks are discussed in detail in appendix D. FeynMaster 2 was also used for two additional tasks. The first was a confirmation that the model is renormalized; by considering processes of the different sectors, we numerically checked that all the counterterms calculated in section 5.1 lead to finite results. The second task was a conversion of the results to Fortran, which were then numerically evaluated using LoopTools [91].
For the scatter plots of section 6.4 below, we generated points in the parameter space of the C2HDM, assuming that h 1 corresponds to the SM-like Higgs boson, such that m 1 = 125 GeV. We imposed both theoretical and experimental restrictions on the points. From the theory side, we required them to comply with boundeness from below [92], perturbative unitarity [93,94], vacuum globality [95] and the oblique parameters S, T and U [13]. As for the experimental restrictions, we required points to have a 2σ compatibility with the results coming from different experiments. One of them is B → X s γ [96][97][98][99][100], which forces m H ± > 580 GeV (for the Type II model); the constraints from R b are also included [96,101]; concerning the different signals for the SM-like Higgs boson, we required compatibility with the fits presented in ref. [102], whereas exclusion bounds from additional Higgs searches are taken into account via HiggsBounds5 [103]. We also considered experiments that restrict the amount of CP violation of the C2HDM; we used the most stringent limit on the electron EDM, |d e | < 1.1 × 10 −29 e cm at 90% confidence level, provided by the ACME collaboration [104]. The Higgs boson production cross sections were calculated with SusHiv1.6.0 [105,106]; for the ratio of rates, we followed the procedure described in ref. [31]. Finally, we allowed the not-fixed input parameters from eq. 3.25 to vary according to:

The influence of the different combinations
The renormalized NLO amplitude for the process j (with j = {h 2 → ZZ, h 2 → h 1 Z, h 2 → h 1 h 1 }), which we represent asM j , can generically be written as: where M tree j andM loop j respectively represent the tree-level amplitude and the renormalized one-loop amplitude, the latter of which can in turn be split according to: where the terms on the right-hand side represent the non-renormalized one-loop amplitude and the counterterms, respectively. Here, M CT j includes all the counterterms that, after expanding the bare quantities into renormalized plus counterterm, end up contributing to the process j. Now, the set of counterterms that constitute M CT j is just a consequence of such an expansion, and is the same for all the combinations C i . The differences between the combinations start when the counterterms are to be given an expression, since the latter in general depend on the combination. For example, suppose that δζ 1 contributes to M CT j ; recalling eq. 3.26, it is clear that the expression for this counterterm depends on C i . Indeed, in C 4 , δζ 1 is an independent counterterm, so that it is determined by eq. 5.16. In the remaining combinations, in contrast, δζ 1 is dependent, which means that it is determined by a function (which in general depends on the particular combination at stake) of independent counterterms. As a consequence, M CT j will in general have different values according to the combination C i chosen. In this way, the different C i in general lead to different predictions for the renormalized NLO amplitudes. At this point, it is worth noting that, in a renormalized function such that all its counterterms are fixed through OSS or other momentum-dependent subtraction schemes, the dependence on the renormalization scale µ R drops out. 41 Conversely, if there is at least a counterterm which is fixed in MS, such counterterm will not cancel the dependence of the non-renormalized function on the renormalization scale, so that the renormalized function at stake will depend on µ R . In section 5.1, all independent counterterms were fixed through a momentum-dependent subtraction scheme, except for δµ 2 and (whenever it is independent) δζ 1 . 42 These exceptions will be quite relevant for the impact of the different combinations on the three processes j introduced above. In order to better illustrate such impact, we show in table 3 all the independent counterterms contributing to the different M CT j , for the four combinations C i . 43 Some comments are in order. First, in all three processes, the complete set of independent counterterms always depends on the combination C i (note that the last four columns of the table reflect eq. 3.26). Moreover, processes h 1 → ZZ and h 2 → h 1 Z should have no scale dependence in C 1 , C 2 and C 3 , since they do not involve any counterterm fixed in MS in those combinations. By contrast, the combination C 4 in all processes is expected to lead to scale dependent results, as it includes the MS-fixed δζ 1 . Finally, process h 2 → h 1 h 1 depends on δµ 2 , so that we expect scale dependence in all combinations of that process.

Decay widths
Using FeynMaster 2, the expressions for the renormalized NLO decay width Γ NLO j for the process j can be easily calculated from the renormalized NLO amplitudeM j . We assume 41 A simple way to see this is to note that, in a one-loop integral calculated through dimensional regularization, the UV divergent term ( 2 ) and the renormalization scale dependent term (ln µ 2 R ) always come together ( 2 + ln µ 2 R ). Now, a momentum-dependent subtraction scheme is such that the counterterms are defined as a function of a one-loop integral evaluated at a specific momentum. Thus, since neither 2 nor ln µ 2 R depend on the momentum, the counterterms will contain both. And just as they (i.e. the counterterms) assure that the renormalized function is UV finite (i.e. assure that the 2 of the counterterms cancel against that of the original one-loop integral), they also assure that it does not depend on the renormalization scale (i.e. the ln µ 2 R cancel). A final note: the renormalization scale µR should not be confused with the parameter µ 2 introduced in eq. 2.45. 42 Whenever δζ1 is fixed in MS, δζ2 will also be (cf. eq. 5.18). 43 These counterterms are fixed according to the prescriptions presented in section 5.1. We should mention that δζ2 contributes to all processes, in all combinations; and since it is an independent counterterm, it should in principle be included in table 3. However, we fixed it by reference to δζ1 (recall eq. 5.18), which is only independent in C4. In the remaining combinations, the counterterms which δζ1 depends on are already accounted for in the table. Hence, for the sake of clarity, we do not show δζ2.

Process j
Common to all the C i Table 3: Total set of independent counterterms contributing to M CT j for the combinations C i . For each process j, the counterterms are separated in two groups: those that are common to all the four combinations (second column), and those that are specific to each combination (last four columns).
that all external particles are on-shell (OS). The final expressions are considerably clear, as long as they are written in terms of form factors. 44 We start with the decay h 2 → ZZ. Defining the momenta and Lorentz indices such that h 2 (p 1 ) → Z ν (q 1 )Z σ (q 2 ), we have: where the different f k correspond to form factors. 45 These in general contain contributions from the counterterms, so that they will in general have different values in the different combinations. The form factors f 15 , f 18 and f 21 do not contribute, since the fact that the Z bosons are OS implies ε * ν (q 1 ) q ν 1 = ε * σ (q 2 ) q σ 2 = 0. The decay width is: As for the decay h 2 → h 1 Z, we define the momenta and Lorentz indices such that 44 In what follows, we follow the conventions of FeynMaster 2; we omit the Feynman diagrams, as well as the expressions for the form factors. 45 We identify tree-level form factors with the superscript 'tree'. The form factor f33 corresponds to the CP-violating component predicted in ref. [41]. It turns out that it does not contribute to Γ NLO h 2 →ZZ , as can be seen in eq. 6.5.

M loop
h 2 →h 1 Z = ε * σ (q 2 ) f 10 q σ 1 + f 11 q σ 2 + f 12 p σ 1 , (6.6b) but, as the Z boson is OS, f 11 does not contribute. The decay width is: Finally, the decay h 2 → h 1 h 1 is such that the renormalized NLO amplitude is simply given by:M 8) and the decay width is:

Numerical results
At last, we investigate the behaviour of the different combinations on the decay widths using numerical results. First and foremost, we checked that all combinations in all processes lead to gauge independent decay widths. This is in agreement with the discussion presented in section 5.2, and validates the simple prescription for the renormalization of mixing parameters proposed in section 5.1.6. We compare the different combinations by ascertaining how their results affect the numerical stability of the perturbative expansion. If perturbation theory is to be trusted, higher orders should contribute with smaller corrections. In particular, the NLO results should generate small corrections when compared to the LO ones. This leads us to define the quantity [57,107]: where Γ LO j represents the decay width at LO for the process j. In figure 1, we show ∆Γ h 2 →ZZ in percentage as a function of Γ LO h 2 →ZZ for the combination C 1 . Several aspects should be stressed here. First if all, we checked numerically that the results are invariant under the renormalization scale, as expected from the discussion of section 6.2. Besides, the plot shows that HiggBounds-5 hinders points with large values of Γ LO ; this is consistent with the fact that, had h 2 → ZZ a large value of Γ LO , h 2 would probably already have been discovered by now. The majority of points passing all constraints thus obey Γ LO < 0.7 GeV; points with Γ LO almost up to 0.8 GeV are still allowed, albeit very scarcely (a stringent fine tuning is required). All in all, the points describe a smooth and well-behaved pattern; in particular, there is no region of parameter space generating singularities in ∆Γ. Finally, although there are points with small ∆Γ for all allowed values of Γ LO -indicating that, for those points, perturbation theory is valid-, larger values of ∆Γ are also possible for smaller values of Γ LO . 46 For those points, one 46 From eq. 6.10, larger values of ∆Γ are expected for smaller values of Γ LO . would in principle need to calculate the following order in perturbation theory to check whether or not a perturbative description is possible. Similar to figure 1 are the two panels of figure 2, where we show equivalent results for the combinations C 2 and C 3 . We checked that, as before, there is no scale dependence. The values of ∆Γ h 2 →ZZ are essentially the same as in figure 1; moreover, we find the same smooth pattern in the points. This allows us to attest the quality of our prescription for the determination of the counterterms for the mixing parameters, proposed in section 5.1.6. 47 Indeed, from table 3, the combinations C 1 , C 2 and C 3 for the process h 2 → ZZ, 47 Except for δζa, which we do not investigate in this paper. besides depending all on δβ, δα 1 , δα 2 and δα 3 , depend individually on δα 5 , δα 4 and δα 0 , respectively. Then, by considering figures 1 and 2, we conclude that all these counterterms lead to well-behaved results. Actually, not only are the results all well-behaved, but they are also very similar. This is consistent with the circumstance that δα 5 , δα 4 and δα 0 were all determined via symmetry relations, so that one would expect similar behaviours from the different lines of eqs. 5.12. All in all, then, the combinations C 1 , C 2 and C 3 are essentially equivalent for the process h 2 → ZZ.
The results for C 4 are shown in figure 3. Here, as expected from the discussion of section 0.4 GeV, these turn out to be quite similar to those of the points passing all constraints in the remaining combinations: the ranges of values of ∆Γ are essentially the same, and the same well-behaved pattern is observed. This could hardly be expected from the expressions for the different counterterms; not only is the specific counterterm in C 4 (δζ 1 ) fixed through a method completely different from the one used to fix the specific counterterms of the combinations C 1 , C 2 and C 3 (δα 5 , δα 4 and δα 0 , respectively), but the role of counterterms is also different: δζ 1 is a counterterm for a phase of the potential, whereas δα 5 , δα 4 and δα 0 are counterterms for mixing parameters. For smaller values of Γ LO h 2 →ZZ , however, and whichever the scale chosen, C 4 leads to more unstable results than the remaining combinations. In figure 4, we show the results for the decay h 2 → h 1 Z, for the four combinations. The conclusions are similar to the ones for h 2 → ZZ, which we just discussed. In particular, the . This similarity is expected from table 3, which shows that the only difference in the set of independent counterterms contributing to the two decays concerns field counterterms; but since all field counterterms are fixed through the same method (namely, OSS), it is not surprising that the two processes turn out to lead to similar NLO corrections. In figure 5, we show the results for h 2 → h 1 h 1 , for the four combinations. Here, as predicted in section 6.2 above, all combinations are scale dependent, due to the contribution of the counterterm δµ 2 ; we plot the intermediate scale only, which we checked is again the most stable one. The four combinations lead to similar results; a well-behaved pattern is once again to be observed in all of them. In contrast to what happened previously, HiggBounds-5 does not preclude points with large values of Γ LO h 2 →h 1 h 1 ; we checked that these points lead to small values of Γ LO h 2 →ZZ and Γ LO h 2 →h 1 Z , and that only large values of these would be excluded by the experimental results. Moreover, values of ∆Γ as high as 40% are allowed for all the region of allowed points in all combinations. We found that more stable results can be obtained by requiring smaller values for the different λ i (parameters of the potential), as can be seen in figure 6 for the combination C 1 , where we constrained |λ i | < 2 for all λ i . This can be explained by the circumstance that smaller values for the λ i correspond to a more stringent constraint regarding the tree-level perturbative unitarity. Hence, it is not surprising that such values lead to a more stable perturbative expansion [54,107].

Conclusions
We presented the renormalization of the C2HDM, which is one of the simplest models containing an extended scalar sector with CP violation. We showed in detail how the presence of CP violation in the scalar sector leads to a rather peculiar program of renormalization, due to two main reasons. First, several non-physical parameters must be introduced in order to assure the generation of the necessary counterterms. This happens in such a way that, once the counterterms are obtained, the corresponding renormalized parameters are not needed and can be rephased away. Second, for the same set of independent physical parameters, multiple combinations C i of independent counterterms can be chosen. It thus becomes relevant to study the impact of different combinations on observables. We considered four combinations and applied them to three specific one-loop decays of h 2 ; the first three combinations choose as independent a counterterm for a mixing parameter, and they fix it through symmetry relations, whereas C 4 takes a parameter of the potential as independent, and fixes it through MS-thus yielding a dependence on the renormalization scale µ R . We checked that all combinations in all processes lead to gauge-independent and well-behaved results, which validates the prescriptions used to fix the counterterms. For all processes, we found that the combination C 4 leads (at least for some values of µ R and of the LO decay width) to results very similar to those of the remaining combinations. Among the points passing all constraints in all processes, one finds regions in the four combinations that lead to a numerical stability of the perturbative expansion; points leading to large NLO corrections are also allowed, in which case the results of higher orders should be investigated. We stressed several aspects related to an adequate treatment of the C2HDM when considered up to one-loop level. We used the Fleischner-Jegerlehner tadpole scheme to assure the selection of the true vev; this was the first time that this scheme was applied to a model with CP violation in scalar sector (and thus to a model with generally complex vevs). We discussed the renormalization of a field whose mass is a dependent parameter. We showed that the OSS conditions in the fermion sector leave some freedom, which can be used to define complex counterterms for the fermion masses. We proposed a simple and easily generalizable prescription to calculate counterterms for mixing parameters, which assures that observables are gauge independent. Finally, FeynMaster 2 was shown to be an ideal tool for renormalizing models such as the C2HDM and study them at NLO, since it calculates a multiplicity of crucial elements-Feynman rules, counterterms, one-loop amplitudes, decay widths, etc.-in a simultaneously automatic and flexible way. Both FeynMaster 2 and the model file for the C2HDM can be downloaded at: https://porthos.tecnico.ulisboa.pt/FeynMaster/ Given its pioneering character, this work opens the door to an exploration of precise predictions in the C2HDM, as well as to the renormalization of models with CP violation in the scalar sector. One interesting direction of future investigation concerns processes in the C2HDM with external h 3 fields, whose pole mass is in general different from the renormalized mass.

A OSS and dependent masses
The on-shell subtraction (OSS) scheme, originally proposed in ref. [83], has become one of the most common subtraction schemes in the renormalization of gauge theories, having been recurrently described and used by many authors [44-47, 78, 81, 108]. In this appendix, we present a clear and simple description of the one-loop renormalization conditions that characterize the OSS scheme, explaining how they allow to derive the mass and field counterterms. We consider both the scenario with no mixing of fields, as well as that with mixing; we also study the scenarios where a dependent mass prevents the application of the OSS conditions, and we show how the counterterms are determined in those cases. This appendix can thus be seen as an introduction to sections 5.1.1 to 5.1.3 above. We will assume scalars particles, for simplicity; in the case of fermions and gauge bosons, the 2-point functions must be projected onto (i.e. multiplied by) physical states: the spinor u(p) for fermions and the polarization vector ε µ (k) for gauge bosons.

A.1 No mixing, independent mass
We start with the simplest case: a generic single scalar neutral field of a generic theory. Suppose a real scalar field φ with mass m, such that: We assume that m is an independent parameter, i.e. that it is not written in terms of other parameters. When considering the theory up to one-loop level, we follow the usual procedure, namely: we identify the tree-level quantities as bare ones, and split them into renormalized quantities and counterterms, Note that we explicitely represent the renormalized mass m R with an index R. Moreover, no mixed fields counterterms are introduced, since we are assuming that φ mixes with no other field at one-loop. Identifying the up-to-one-loop 2-point function with Γ and the one-loop 2-point function with Σ, we have: where the caret means 'renormalized'. Then, expanding the bare version of eq. A.1 with eq. A.2, we find: Now, the pole mass m P is defined as the value of the momentum for whichΓ(p 2 ) is zero, that is: 49 This allows us to see how the OSS conditions fix the counterterms δm 2 and δZ φ . In the case we are considering (a single field), there are only two such conditions. The first one simply states that the renormalized mass is equal to the pole mass, thus giving a physical meaning to the renormalized mass. In that case, eq. A.5 implies: which, together with eqs. A.3 and A.4, fixes δm 2 to: Recall that, as seen in eq. 4.13, the non-renormalized one-loop 2-point function Σ(p 2 ) in general includes diagrams with one-loop tadpoles if the FJTS is adopted. 49 We are using the definition of real pole mass (also known as Breit-Wigner mass), where mP is real as a consequence of the fact that we are considering the real part ofΣ(m 2 P ). This is what is done in the traditional OSS scheme [44,47]. We are using the operator Re (and not Re), which neglects the imaginaryalso called absorptive-parts of loop integrals, while keeping the imaginary parts of complex parameters. The use of Re in the renormalization conditions is also common practice in the traditional OSS scheme.
One should not confuse eqs. A.5 and A.7: the former is the definition of the pole mass, valid in every subtraction scheme. In contrast, eq. A.7 is specific of OSS, and results from the application of the OSS condition A.6 to eq. A.5. The second OSS condition takes the (renormalized, real) residue of the renormalized propagator at the pole mass, and sets it to one. Actually, one usually sets the inverse of the residue equal to one, which means: 50 lim (A.10) Then, using eqs. A.3 and A.4, as well as L'Hôpital's rule, the expression for δZ φ follows: (A.11)

A.2 No mixing, dependent mass
If we suppose instead that the bare mass m (0) is a dependent parameter, then both the renormalized mass m R and the counterterm δm 2 are also dependent, i.e. fixed. 51 Now, since m R is fixed, it cannot be set to be equal to the pole mass m P , so that the first OSS condition cannot be used. In that case, then, the two masses will in general be different, and eq. A.5 (which is always valid) assures that their difference is of one-loop order. But this means that, to first order in perturbation theory, m 2 P can be replaced by m 2 R when m 2 P is an argument of a one-loop function. In particular, the relation in valid to first order, so that eqs. A.5 and A.4 imply, to that order, Hence, the pole mass is completely determined. 52 Just as in the case of the independent mass, δZ φ can be fixed by setting the residue equal to one. Only, one must be careful with the fact that m 2 P is different from m 2 R . Inserting eqs. A.3 and A.4 inside eq. A.9, equating the latter to one, using L'Hôpital's rule and the relation (A.14) 50 In this equation, mP and mR can be used indifferently, as they are chosen to be the same in OSS. 51 The case of dependent masses is common in the minimal supersymmetric extension of the SM (MSSM); see e.g. ref. [109]. 52 Note that m 2 P is finite, as it should be. A simple way to see this is to recall that the divergent parts of counterterms must be the same whichever the subtraction scheme chosen. Then, the divergent part of the left-hand side of eq. A.8 is equal to that of the right-hand side, for all subtraction schemes. Therefore, the divergent parts of the right-hand side of eq. A.13 cancel.
which holds to first order, we find: which is precisely the same as the field counterterm fixed through OSS in the case of the independent mass, eq. A.11.

A.3 Mixing, all masses independent
We now consider a generalization of eqs. A.1 to A.4 for two fields: where i, j = {1, 2}, and where both m 1(0) and m 2(0) are assumed to be independent parameters. In OSS, the counterterms are once again calculated by reference to the pole masses and the residues. But, now that there is mixing, what is the definition of the pole masses and the residues? BecauseΓ(p 2 ) is now a matrix (eq. A.19), the pole masses are determined as the values of the momentum for which its eigenvalues are zero. And while that leads to complicated expressions, these can be simplified by fixing the off-diagonal field counterterms in an appropriate way. More specifically, if we choose those counterterms to be such that: then the definition of the pole masses m iP is a trivial generalization of eq. A.5, namely, in which case the residues become: Eq. A.21 assures that OS particles do not mix with each other, and allows to calculate the off-diagonal field counterterms (δZ 12 and δZ 21 ) through eqs. A.19 and A.20. Now, OSS assumes not only eq. A.21, but also that the renormalized masses are equal to the pole masses, m iR OSS = m iP (which we can always do, since the renormalized masses m iR are assumed to be free parameters). In that case, eq. A.21 becomes: while eq. A.22 becomes: which allows to determine the mass counterterms. Besides, and also in a trivial generalization of the case with no mixing, OSS fixes the diagonal field counterterms by requiring the different inverse residues to be equal to one: The expressions for the complete set of mass and field counterterms are, thus,

A.4 Mixing, one dependent mass
Finally, we consider the situation where m 2(0) is a dependent parameter (while m 1(0) is still independent). It follows that m 2 2R and δm 2 2 are dependent. Nonetheless, we can always impose eq. A.21, which not only fixes the counterterms δZ 12 and δZ 21 , but also implies that eqs. A.22 and A.23 are valid. Now, while we have the freedom to choose m 1R OSS = m 1P (so that eqs. A.25 and A.26 hold for i = 1, and these respectively allow to calculate δm 2 1 and δZ 11 ), we do not have the freedom to fix m 2 2R . In this case, though, we have a trivial generalization of section A.2: the expression for m 2 2P follows from eq. A.22, while δZ 22 can be determined by fixing the residue of the OS propagator of φ 2 to one. 53 The expressions then become: where we neglected higher order terms. Hence, eq. A.29 only differs from eq. A.27 in the fact that δm 2 2 is a priori fixed. The generalization to the case with four fields with one dependent mass (as in section 5.1.3 above) is straightforward.

B CP violation in fermionic 2-point functions
In this appendix, we study CP violation in fermionic 2-point functions. As we will show, this is particularly relevant for renormalization purposes-and hence for the renormalization of the C2HDM. However, we do not restrict ourselves to a particular theory; in fact, the appendix can be applied to every model with Dirac fermions. We organize it as follows: in a first part, we show how a complex mass shows up in the Dirac representation of the Lorentz group; then, we consider CP violation in fermionic 2-point functions and prove that a complex mass is intrinsically related with it; after that, we investigate the consequences of a CP-violating one-loop 2-point function on counterterms; finally, we ascertain how to handle CP violation using the OSS scheme.

B.1 Complex mass
Let us consider two fields: a left-handed Weyl spinor f w L -transforming in the 1 2 , 0 representation of the Lorentz group-and a right-handed Weyl spinor f w R -transforming in the 0, 1 2 representation of the Lorentz group. 54 The Lagrangian for the mass term is such that: − Note that the parameter m f is in general complex. Moreover, f w L and f w R are two separate and unrelated fields, which implies in particular that they can be rotated in independent ways. Later on, they (or some rotated versions of them) will be embedded in a Dirac spinor as its left-handed and right-handed components. But until then, f w L and f w R are not connected with each other in any way. Now, the fact that we can redefine them differently means that the phase of m f can be absorbed in the fields. Indeed, defining m f = |m f |e iθm , we can for example rephase f w R as to render m f real. Such rephasing (or another one equivalent to it) is usually performed, since m f is usually identified with the (real) pole mass. Yet, it is illustrative to keep m f complex. In particular, we want to ascertain how a complex mass shows up in the Dirac representation 1 2 , 0 ⊕ 0, 1 2 of the Lorentz group. Thus, we now embed the Weyl spinors in a Dirac spinor f , such that: In the Weyl representation, 54 We consider a single flavour for simplicity.
which means that eq. B.1 can be rewritten as: Separating m f into its real and imaginary parts, we find: Using the properties γ † 5 = γ 5 and {γ 5 , γ µ } = 0, it is straightforward to conclude that so that eq. B.6 becomes: In the Dirac representation, then, a mass m f with non-null imaginary part is equivalent to the presence of a term of the formf γ 5 f . 55 We now show that such term violates CP.

B.2 CP violation in 2-point functions
We start by considering the 2-point function for f . We can parameterize it as [47]: or equivalently as with the correspondence (B.11) Now, eq. B.10 corresponds to the effective Lagrangian In Table 4, we show the action of the operators C, P and CP on the four bilinear structures of eq. B.12, but for general fields ψ i and ψ j . In the particular case ψ i = ψ j = f , there is CP violation in eq. B.12 if and only if Γ f,P = 0. Using eq. B.11, we thus conclude: Comparing eqs. B.8 and B.12, it is clear that the term in eq. B.8 proportional to Im[m f ] violates CP. One might wonder why this conclusion is relevant; after all, and as we showed, one can always avoid such term by a rephasing of the Weyl spinors. But in that case, what happens at loop level? In particular, if we perform such rephasing and if CP violation contributes to fermionic 2-point functions at one-loop, how are we assured that there will be a counterterm for such contribution? It is to these questions that we now turn. 55 This term has no physical meaning whatsoever, as it depends on the basis chosen: given eq. B.2, that term will not show up if, before we embed the Weyl spinors f w L and f w R in the Dirac spinor f , we rotate them to render m f real. Table 4: Results of the application of the symmetries C, P and CP to the terms in the left column.

B.3 Fermionic counterterms and CP violation
When renormalizing the theory, the tree-level parameters are identified as usual with bare parameters. Since the tree-level parameter m f was in general complex, the bare parameter m f (0) is also in general complex. Writing eq. B.5 in terms of bare quantities, we have: such that Here, m f and δm f are in general complex. However, based on what we showed one could do at tree-level, we can take the freedom in the basis choice of the (renormalized) chiral spinors f L and f R to render m f real. After that, δm f is still in general complex. 56 Separating the renormalized terms and the counterterms, eq. B.14 then becomes: Eqs. B.17 clarify the importance of the discussion of the previous section: we see that, due to our basis choice for the chiral spinors, there is no CP violation in the renormalized terms, eq. B.17a. But it is also clear that, even after that basis choice, there is CP violation in the counterterms, eq. B.17b, as long as Im[δM f ] = 0 (recall eqs. B.12 and B.13). Therefore, if CP-violating effects show up at one-loop level, there will in general be counterterms to absorb their divergences.

B.4 CP violation in the OSS scheme
We now focus on the consequences of the OSS renormalization conditions on the different terms of δM f . We start by noting that, motivated by the fact that we make the identification: which implies, in particular, Now, the fermionic renormalized 2-point function reads: and we can use the parameterization of eq. B.9 for renormalized quantities: which, in the particular case of the renormalized one-loop GFs, reads: Expanding the Lagrangian using eqs. 3.5 and 3.7, we find so that, using eqs. B.22, B.24 and B.25, we obtain: The OSS renormalization conditions are (no summation in the indices is implicit): Inserting eq. B.23, we have respectively: Then, inserting eqs. B.26 in the two eqs. B.28a for i = j, we obtain respectively: 57 Now, the hermiticity of the up-to-one-loop action implies: so that, in particular, Hence, inserting eqs. B.26 in eq. B.28b (again for i = j), and subtracting the complex conjugate of eq. B.29a to eq. B.29b while using eqs. B.21 and B.31, we get: There are thus 6 real degrees of freedom (dof)-two for each of the three complex counterterms δm f , δZ f,L and δZ f,R -, but the diagonal OSS conditions only imply the 4 relations of eqs. B.29 and B.32. Hence, there are two free dof, so that many solutions satisfy those conditions. Before considering them, note that, for the specific case of the non-renormalized diagonal one-loop 2-point function, eq. B.13 implies: Now, among the many possible solutions, we restrict ourselves to two. The first chooses Im δZ f,L = Im δZ f,R = 0. In this case, eq. B.33 together with eqs. B.29 and B.32 leads to: This relation shows that, whenever there is no CP violation in one-loop fermionic diagonal 2-point functions, there is no need to introduce δm L f nor δm R f : they would be the same, and both equal to δm f , which is thus preferred. Conversely, if there is CP violation in Re Σf f , the choice Im δZ f,L = Im δZ f,R = 0 forces δm L f and δm R f to be different (or, what is equivalent, forces δm f to be complex). This is the solution we adopted in eqs. 5.4 for the case i = j. 58 Another simple solution uses one of the two free dofs to impose δm L f = δm R f . In this case, by subtracting eq. B.29a to eq. B.29b, taking the imaginary part of the result and using eqs. B.31 and B.33, we get: It is then clear that, in the context of OSS, CP violation in fermionic diagonal 2-point functions does not necessarily force the relation δm L f = δm R f . Only, by choosing δm L f = δm R f , it follows that δZ f,L and δZ f,R cannot be both real. This is the solution usually adopted in the cMSSM (cf. e.g. refs. [67,68]).

C Symmetry relations
To our knowledge, the procedure of fixing counterterms through symmetry relations has been originally proposed in ref. [50], having been later worked out and used by several authors [52,53,56,57]. As mentioned in section 3, the theory can be renormalized in its symmetric form (before SSB); in that case, it is enough to consider one field counterterm for each multiplet in order to cancel UV divergences [47]. On the other hand, the theory can also be renormalized in the broken phase, where counterterms for mixing parameters and mixed fields counterterms in general appear. The two methods are not independent whenever fields in the symmetry basis are related with fields in the physical basis through mixing parameters. In this appendix, we start by showing that one can exploit such dependence to write the divergent parts of the counterterms for mixing parameters in terms of the divergent parts of field counterterms; this we do in section C.1, considering a simple example. Then, in section C.2, we apply the same method to the specific case of the mixing parameters of the scalar sector of the C2HDM, and discuss the limit of CP conservation. The symmetry relations as such only have to be true for the divergent parts, since the common ground between the renormalization in the two forms (symmetric and broken) was the elimination of UV divergences. Yet, as we will see in section C.3, in the case of independent counterterms for mixing parameters (which need to be fixed) one can decide to fix them by requiring that the symmetry relations also hold for the finite parts. This is the reason why (independent) counterterms for mixing parameters can be fixed through symmetry relations. In that final section, we discuss the generality of the method, and consider the consequences of symmetry relations on dependent counterterms.

C.1 A simple example
Let us consider an example involving SM particles only. From the bare version of eq. 2.29, we have: According to the previous paragraphs, from the point of view of the cancellation of divergences, W a µ and B µ can be renormalized with a single counterterm each. Hence, 59 so that, in particular, On the other hand, from eqs. 3.4 and 3.14, (C.4) and the renormalized quantities are defined such that: As a consequence, eq. C.3 becomes: (C.6) 59 We define the counterterm δZ W with a prime to distinguish it from the counterterm for the physical On the other hand, we can apply eq. C.4 to the right-hand side of eq. C.1 and consider the divergent parts only to get: Finally, equating eqs. C.6 and C.7 to first order and noting that δθ w = −δc w /s w , one concludes that: This is the symmetry relation; it fixes the divergent part of the counterterm for the mixing parameter θ w in terms of the divergent parts of the counterterms for the physical fields associated with θ w . This relation is also valid in the SM, as well as in all its extensions that keep its gauge structure. 60 The procedure just used is generalizable to other mixing parameters (or sets of mixing parameters) of different models and can be summarized in four steps. First, one considers the relation involving the mixing parameters at stake (i.e. transforming the fields in the symmetry basis into the mass basis) in its bare form; in the above example, this is simply eq. C.1. Second, the bare fields in the symmetry basis are expanded-aiming at the elimination of divergences, so that only one field counterterm for multiplet is taken and only the divergent parts are considered-and, after that, the renormalized relation between the two bases is used; in the above example, the result of this step is eq. C.6. In the third step, one expands instead the bare mixing parameters and bare physical fields and takes the divergent parts; in the above example, the result is eq. C.7. Finally, one equates the two previous steps to first order, which yields the divergent parts of the counterterms for the mixing parameters in terms of the divergent parts of field counterterms.

C.2 The case of the C2HDM
Applying this procedure to the scalar sector of the C2HDM, we get, for the counterterms for mixing parameters of the charged fields, 61 60 The same is true for another relation that also follows from equating eqs. C.6 and C.7 to first order, namely, (C.9) 61 Besides, we also find the following relations for the field counterterms: Had we started with a parameterization of X that included an overall phase, the symmetry relations would have forced the divergent part of its counterterm to be zero. This means that such phase is really not needed, which explains why we ignored it in the first place (cf. eq. 2.10 and note 4).
and, for those of the neutral fields, Note that the expressions for the divergent parts of δα 0 , δα 4 and δα 5 depend only on the divergent parts of fields counterterms for the mixing between G 0 and other neutral scalar fields. This means that, should there be no such mixing at one-loop, the divergent parts of the corresponding mixed fields counterterms would vanish-which would then imply the vanishing of the divergent parts of δα 0 , δα 4 and δα 5 . In other words, if G 0 did not mix with the fields h 1 , h 2 and h 3 at one-loop, the counterterms δα 0 , δα 4 and δα 5 would not be needed. This is consistent with description of the theory considered solely at tree-level, where the angles α 0 , α 4 and α 5 are not introduced. However, since G 0 does mix with the remaining neutral scalar fields at one-loop, the divergent parts of δα 0 , δα 4 and δα 5 are in general not zero, so that these counterterms must be considered-which in turn means that α 0 α 4 and α 5 have to be included at tree-level when aiming at the one-loop renormalization of the theory. It is also interesting to consider the limit of CP conservation; one possible limit is [17]: in which case [17]: with h and H being CP-even and A 0 CP-odd. Then, given that there would be no one-loop mixing between fields with different CP values, eqs. C.12 would become: This is consistent with the symmetry relations in a 2HDM with CP conservation (cf. e.g. ref. [52]). Moreover, still in the context of CP conservation, Z G + H + would be real, which would imply that δζ a ∆ would vanish; consequently, if CP were a symmetry of the theory, there would also be no need to consider δζ a .

C.3 Dependent and independent mixing parameters
The symmetry relations have interesting consequences both when the counterterm for the mixing parameter at stake is dependend and when it is independent. In the first case, one can use the symmetry relations to establish relations between the divergent parts of different independent counterterms. For example, given that [46,80]: one can use eq. C.8 to conclude that: Note that the same relation is in general not valid for the finite parts. 62 In the case of counterterms for independent mixing parameters, the finite parts are a priori not fixed. In fact, even if we are able to determine their divergent parts, the finite parts are free, and must be fixed by renormalization conditions. But here is the trick: we can choose as renormalization conditions the symmetry relations including finite parts. Said otherwise, we can decide that the way through which we fix the independent counterterms is requiring that the symmetry relations hold not only for the divergent parts, but also for the finite parts [56]. In particular, we can decide that the relations involved in eqs. C.11 and C.12 also hold for the finite parts-whenever the counterterm for mixing parameters at stake is independent (which depends on the combination C i ). This fixes the counterterms involved, thus justifying the expressions presented in section 5.1.6. 63 The prescription we just described-of promoting symmetry relations for independent mixing parameter counterterms to renormalization conditions-is completely general, not being restricted to the C2HDM. In fact, just as the procedure presented in section C.1 is easily generalizable to other models, so is this prescription to fix independent counterterms for mixing parameters. Actually, not only it is easily generalizable, as it is very simple, leading to gauge independent observables (as discussed in section 5.2). 62 This is the reason why the counterterms involved are still independent (even if their divergent parts are related). Again, eq. C.19 is also valid in the SM, as well as in all extensions of this model that do not modify its gauge structure. 63 As discussed in section 5.2, the independent counterterms for the mixing parameters are calculated in the Feynman gauge.

D FeynMaster 2: an ideal tool for renormalizing the C2HDM
In this appendix, we introduce FeynMaster 2 and describe its application to the renormalization of the C2HDM. We start by presenting the main differences regarding the previous version of the software. We then summarize its usefulness in the renormalization of models. Finally, we explain how it can be exploited to renormalize a model such as the C2HDM, providing illustrations of the different steps involved.

D.1 The new version
An overall description of FeynMaster was already presented in the Introduction. The second version of the program, FeynMaster 2, has been very recently made available online in the webpage: https://porthos.tecnico.ulisboa.pt/FeynMaster/.
There, the new software can be downloaded, and a complete and auto-sufficient manual can be found. FeynMaster 2 has significant improvements over the previous version of the program. First, it is much more practical, since a) a model is now specified by a single model file, b) the file controlling the FeynMaster run now contains just the data relevant for the run, and c) the presence of QGRAF has virtually vanished. The last point means that there is no longer such thing as a QGRAF style or convention, which in turn implies that there is a single set of conventions (the FeynRules ones) used throughout the whole program. A second major improvement in FeynMaster 2 is that the calculated results can now be automatically stored. This implies that a certain process can be calculated once and for all, which considerably simplifies the manipulation of results and the use of FeynMaster as a whole (this improvement is especially important for renormalization, as shall be seen). In the third place, the renormalization is significantly faster; for example, the generation and impression of the total set of Feynman rules (tree-level and counterterms) for the SM is now completed in just 5 minutes in normal laptops. Finally, several functions were corrected, improved or simply created.

D.2 An ideal tool for renormalizing models
Having seen this, and given the list of tasks presented in the Introduction, it is easy to conclude that FeynMaster 2 is an extremely useful tool in the renormalization of a model. Such usefulness is manifest, first and foremost, in the fact that the different steps required for the renormalization of a model-definition of the model, generation of Feynman rules, generation of counterterms, one-loop calculations and renormalization conditions-can all be performed with FeynMaster 2. One starts by writing the model, with the desired conventions. This allows to automatically generate the Feynman rules not only for the tree-level interactions, but also for the counterterms. With the former, one-loop processes can be automatically calculated and stored, which allows to fix all the counterterms in a desired subtraction scheme, thus completing the renormalization of the model.
But this is not all. One can then study different one-loop processes (meanwhile renormalized), disposing not only of all the tools of FeynCalc, but also of some specific functions of FeynMaster 2. These can also be combined with the numerical interface with Fortran, thus allowing a more complete analysis. In the end of the day, renormalization and the study of NLO processes has become quite simple. A constant element in the chain of processes just described should be stressed: flexibility. This property begins in the moment when, in the writing of the model, the user defines the conventions. Not only he or she has the freedom to define those conventions at will, but-and what is perhaps more relevant-they are kept throughout the whole process: in Feynman rules, one-loop calculations, counterterms, etc. Indeed, by combining all the required steps in a single software, FeynMaster 2 assures that the conventions do not change when changing from one step to the other. 64 Flexibility is present in many other ways. Intermediate and final results can be manipulated: not only due to the intuitive structure of FeynMaster 2, but also to the user-friendly character of FeynCalc and Mathematica. Notebooks for processes are automatically written, in which all of the results already calculated (Feynman rules, counterterms, oneloop amplitudes, etc.) are immediately available. The writing of Feynman diagrams and analytical expressions in L A T E X is immediate. The treatment of Dirac structures and the calculation of (one-loop) decay widths and cross sections is remarkably simple given the functions available in FeynMaster 2. The numerical interface includes the beginning of a Fortran main file (adapted to the process at stake), as well as a template of a makefile. These examples suffice to make our point: FeynMaster 2 is very convenient for renormalizing a model and to study it at NLO. We now illustrate these features in the particular case of the C2HDM.

D.3 The C2HDM
We organize this section in three parts. In a first one, we briefly describe some particularities concerning the implementation of the C2HDM as a model for FeynMaster 2. Then, we explain how to calculate the counterterms. In the end, we show how to use FeynMaster 2 to study processes at NLO.

D.3.1 The model
The C2HDM as a model for FeynMaster 2 is available online in the webpage of Feyn-Master 2 (in the Type II version, cf. section 2.3). Here, we do not describe in detail all the different components of the model (the manual can be checked for details). Rather, we stress some aspects which may be less obvious. First, the model for the C2HDM is an extension of the model for the SM, in the sense that it was built upon the model file for the SM. This implies, in particular, that it also contains the η's from ref. [77], which are very useful to compare different conventions. Now, the theory (the Lagrangian, with its particles and parameters) is originally written as a bare theory. And since we intend to renormalize it, we need to write it in a general basis. The bare quantities are identified with the sum of a renormalized quantity and its corresponding counterterm, as usual. The renormalized quantities obey eqs. 3.19 to 3.24, which are enforced through the restrictions file. Counterterms such as δλ 1 will show up because the model is written in such a way that most of the dependent parameters are renormalized. 65 Although we did not need to renormalize them (they are dependent, so that one could rewrite them in terms of independent parameters, and renormalize the latter), we decided to renormalize them for convenience. This is indeed very convenient: although it introduces several unnecessary counterterms, it renders the expressions considerably more compact. 66 By running FeynMaster 2 with the variables FRinterLogic and RenoLogic (of the Control.m file) set to True, the Feynman diagrams and rules for the tree-level and counterterm interactions are automatically generated. In figures D.1 and D.2, we present excerpts of the files containing the diagrams and rules written in L A T E X. The corresponding files with the internal rules (i.e. the rules written in a FeynCalc-readable format, to be used in the calculations) are also automatically generated. We show in figure D.3 an excerpt of one of them, CTini.m, which is the one concerned with the counterterms. There, each non-empty line begins with the interaction at stake, containing after it the total set of counterterms contributing to it. Note the presence of dependent counterterms, such as δQ 41 and 65 Recall note 21 and eq. 3.14. The circumstance that we renormalize m 2 11 , m 2 22 and m 2 12I for convenience implies that FeynMaster will generate counterterms for these variables, in such a way that it does not know their dependences-i.e. does not know how they depend on other counterterms. In that case, counterterms for the linear terms in h1, h2, h3 and G0 are apparently generated. Only, they are zero (which can be proven by replacing the bare values of m 2 11 , m 2 22 and m 2 12I for their expressions, and renormalizing the latter). This is a consequence of the fact that we are using the FJTS. For more details, see ref. [80]. 66 Actually, if no unnecessary counterterms were used, the renormalization would probably be unfeasible: as it is, it takes about 50 minutes in our laptops to generate the total set of Feynman rules for the tree-level as well as for the counterterm interactions. It would take much longer if unnecessary counterterms were avoided.  δs w (represented by dQ41 and dsw, respectively). Note also the presence of renormalized parameters such as R21, corresponding to the entry 21 of the matrix R. 67 Finally, we have set the option M$PrMassFL of the FeynMaster model file (defined in the Extra.fr file) to False. This means that FeynMaster does not extract the mass of a certain propagator from the term in the Lagrangian which is bilinear in the corresponding field, but rather writes a simplified version of the propagator. 68 But it also means that the gauge boson propagators are written in the Feynman gauge. In order to study gaugedependences, however, they must be written in an arbitrary R ξ gauge. As the totality of propagators that depend on the gauge exist also in the SM, we can copy them-assuming the user has already generated the Feynman rules for SM-from the file Propagators.m of the SM (which lies inside the folder FeynmanRules, contained in the FeynMaster output 67 We used eqs. 3.22 and 3.23 to replace the elements of matrix Q by those of matrix R. We could also have used eq. 3.21 to rewrite the entries of the matrix R, but the Feynman rules would then become longer. 68 Otherwise, the masses in the propagators of the scalar particles would be written in terms of parameters in the potential. For more details, cf. the manual.

D.3.2 Calculation of counterterms
As we just saw, the running of FeynMaster already generated the totality of Feynman rules for the counterterm interactions. That is, for each process, the total set of counterterms contributing to it is already known. Now, those counterterms need to be calculated, which can be done in different ways, i.e. through different subtraction schemes. Here, we adopt the schemes defined in section 5. We thus need to calculate several GFs: mostly 2-point functions, but also some 3-point functions (due to MS renormalization). This can be done in a single FeynMaster run, as illustrated in figure D.5. Some clarifications are in order, concerning this figure.
First, we only show some of the GFs involved in section 5; the remaining ones should be added in the place of the ellipses "(* ... *)". Second, in the calculation of Σb b , we exclude diagrams with gluons via {avoid,G,1,9}; the same applies for all GFs with external fermions. Third, note the options noPVauto and nosigma. The former must be employed in some GFs to prevent meaningless divergences. 69 The latter assures two things: on the one hand, that no one-loop corrections to the external legs are included (such corrections 69 The conditions of section 5 predict that some counterterms are calculated taking the limit of zero external momentum in some 2-point functions. Now, as explained in the manual of FeynMaster 2, loop integrals are performed with the function OneLoopTID, which employs the FeynCalc function TID. By default, OneLoopTID uses TID with the option PaVeAutoReduce -> True, which simplifies some special cases of Passarino-Veltman functions. But it turns out that, in a general R ξ gauge, that option causes spurious divergences to appear in some amplitudes when the limit of zero external momentum is taken. We verified this phenomenon for the cases Σ AA , Σ AZ , Σ G 0 G 0 and Σ G + G + . This can be solved by making PaVeAutoReduce -> False, which in turn can be obtained by selecting the option noPVauto in the Control.m file. Finally, note that nothing of this happens in the Feynman gauge: all the counterterms are well-behaved in the limit of zero external momentum in that case. are not needed for renormalizing purposes); 70 on the other hand, that reducible diagrams with broad tadpoles are included, as is required by the FJTS (recall section 4). We thus prove what we claimed in that section (4), namely, that the inclusion of broad tadpoles in FeynMaster is trivial. Once all the GFs required in section 5 are calculated, the total set of counterterms can be directly determined in an analytical form through the results derived in that section.

D.3.3 NLO processes
As suggested above, the usefulness of FeynMaster does not end when the counterterms are determined. After that, multiple tasks are enabled by the tools provided in FeynMaster. We now briefly illustrate three of them. A first one consists in verifying the finiteness of a certain one-loop process. Once in the possession of the expressions for the counterterms, this task is almost immediate. One starts by selecting the process at stake in Control.m with the options Comp and SumLogic set to True. By running FeynMaster, the total divergent part of the process is calculated, and a notebook for the process is automatically generated. One then opens the latter, having 70 Actually, this aspect is only relevant in the case of GFs with more than two external particles. thus immediate access not only to the total divergent part (through the variable resDtot), but also to the total counterterms contributing to the process at stake (through the variable CTprocess, where process corresponds to the names of the incoming and the outgoing particles joined together). The expressions for the counterterms can be imported and replaced inside CTprocess. The divergent part of the resulting expression can be obtained with the function GetDiv and compared with resDtot. 71 If the process is finite, the difference must be zero. 72 A second task that becomes quite simple when using FeynMaster is the calculation of a decay width at up-to-one-loop level. This requires the calculation of the tree-level amplitude, as well as of the renormalized one-loop one. The user can start with the former (choosing the process at stake at tree-level in Control.m with Comp set to True). Then, in a similar way, he or she calculates the non-renormalized one-loop amplitude (if the variable Draw in Control.m is set to True, FeynMaster automatically draws the Feynman diagrams, as in figure D.6). Here, as in the previous paragraph, the total counterterm should be added to obtain the (finite) renormalized one-loop amplitude. After summing up the two amplitudes (tree-level and renormalized one-loop), one will in general have a complicated expression. In this case, although the FeynMaster function DecayWidth can be directly applied to it, it will probably take quite some time; it is thus convenient to rewrite the expression in terms of form factors. This can be done automatically using another FeynMaster function, FacToDecay, which yields a two-element list: the first one is precisely the expression rewritten in terms of form factors; the second one is a list of replacements, associating to each form factor the corresponding analytical expression (written using the kinematics of the process). The function DecayWidth can now be easily 71 As a matter of fact, the counterterm should be compared with i times resDtot, since the amplitude calculated corresponds to M, and not to iM. 72 The comparison can be made analytically or numerically. While the former is more robust, it can take some time. The numerical method, on the other hand, requires a software to generate points in the parameter space of the model.
applied to the first element of the output of FacToDecay. 73 In this way, one obtains in a couple of steps a simple formula for the NLO decay width, as well as a list with the different form factors contained in that formula. A final example concerns the numerical interface of FeynMaster 2, enabled by yet another FeynMaster function, FCtoFT. This function converts the expression it was applied to into Fortran format and generates several files to run a Fortran program (cf. the manual for details). FCtoFT is especially convenient in cases where the user has a simple expression written in terms of form factors and a list of replacement rules for each of them-precisely the case described in the previous paragraph. The reason is that FCtoFT can be used with two arguments, respectively: a simple expression in terms of form factors, and the list of replacement rules for them. In summary, the combination of functions FacToDecay, DecayWidth and FCtoFT is very powerful.