No strong $CP$ violation up to the one-loop level in a two-Higgs-doublet model

We put forward a two-Higgs-doublet model, furnished with a $Z_3$ symmetry, wherein $CP$ is conserved in the dimension-four terms of the Lagrangian and is softly broken in the scalar potential. The new particles of our model are one neutral scalar $H$, one neutral pseudoscalar $A$, and two charged scalars $H^\pm$. In our model the only locus of $CP$ violation is the CKM matrix. Strong $CP$ violation is absent both at the tree and one-loop levels. We work out the phenomenological constraints on our model, which features flavour-changing neutral Yukawa interactions, showing that the new scalar particles may in some cases be lighter than 500\,GeV.


Introduction
Non-perturbative effects in Quantum Chromodynamics (QCD) may lead to P and CP violation, characterized by a parameter θ, in hadronic processes. The experimental upper bound on the electric dipole moment of the neutron necessitates θ 10 −9 . 1 The presence in the Lagrangian of this unnaturally small parameter is known as the 'strong CP problem'.
The angle θ is the sum of two terms, θ QCD and θ QFD . Here, θ QCD is the value of a Pand CP -violating angle in the QCD vacuum, and θ QFD originates in the chiral rotation of the quark fields needed to render the quark masses real and positive. Let p and n denote the three up-type quarks and the three down-type quarks, respectively, in a weak basis. Let the mass terms of those quarks be given by where M p and M n are 3 × 3 matrices in flavour space. Then, θ QFD = arg det (M p M n ). There are two general approaches to solving the strong CP problem. In the first approach it is claimed that θ has no significance or physical consequences; theories with different values of θ are equivalent and one may set θ to zero without loss of generality. This may happen either because one of the quarks is massless 2 or because of the presence in the theory of a Peccei-Quinn symmetry [3]; there are also claims that QCD dynamics itself cures the strong CP problem. 3 The second approach, which we shall follow, acknowledges the strong CP problem and tries to find some symmetry that naturally leads to the smallness of θ. One firstly assumes the dimension-four part of the Lagrangian to be either CP -symmetric or P -symmetric; this assumption sets θ QCD to zero. The CP or P symmetry must be either softly or spontaneously broken; one performs this breaking in such a way that θ QFD turns out to be zero at the tree level, because of some peculiar form of M p and M n . Still, it is difficult to avoid loop contributions to θ QFD arising from the quark self-energies Σ; they add to the tree-level mass matrices M and then arg det (M + Σ) ≈ Im tr M −1 Σ is in general nonzero. Artful models are able to obtain Im tr M −1 p Σ p + Im [tr (M −1 n Σ n )] equal to zero at the one-loop level and sometimes even at the two-loop level.
There are various ways to achieve quark mass matrices displaying arg det (M p M n ) = 0. Most of those ways, collectively known as Barr-Nelson-type models, employ extra quarks. There are also many models for solving the strong CP problem that use extra gauge symmetries, especially the left-right symmetry SU (2) L × SU (2) R × U (1) B−L . In this paper we propose a simple extension of the Standard Model (SM), with gauge group SU (2) L × U (1) and without any extra fermions, that partially solves the strong CP problem. Our model is a two-Higgs-doublet model (2HDM) [5].
In a 2HDM the quark Yukawa Lagrangian is Lj Φ a (Γ a ) jk n Rk +Φ a (∆ a ) jk p Rk + H.c.
where Φ a = φ + a , φ 0 a T andΦ a = φ 0 a * , −φ − a T for a = 1, 2 are scalar doublets of SU (2) L . Furthermore,Q Lj = p Lj ,n Lj , and the Γ a and ∆ a are four 3 × 3 matrices in flavour space containing the Yukawa coupling constants. We expand the scalar doublets as where v a exp (iℵ a ) √ 2 = 0 |φ 0 a | 0 and the v a are non-negative real by definition. We define v = v 2 1 + v 2 2 = 2m W /g = 246 GeV and tan β ≡ v 2 /v 1 ; then, where the angle β is in the first quadrant. There is one physical pseudoscalar A and one unphysical (Goldstone boson) pseudoscalar G 0 : (From now on, s ξ ≡ sin ξ and c ξ ≡ cos ξ for any needed angle ξ.) There is a pair of physical charged scalars H ± and a pair of unphysical (Goldstone bosons) charged scalars G ± : There are two physical neutral scalars h and H: The neutral scalar h is chosen to coincide with the LHC-observed particle with mass 125 GeV.
In our specific 2HDM there are softly-broken CP and Z 3 symmetries such that h and H do not mix with A. The interaction Lagrangian between a scalar and a pair of gauge bosons is where θ w is Weinberg's angle. Because of the LHC data we now know that |s β−α | ≈ 1.
The quark mass matrices are Let the unitary matrices U n,p L,R bi-diagonalize M n and M p as The CKM matrix is We define and Then, the Yukawa interactions in the physical basis are given by where P L = (1 − γ 5 )/ 2 and P R = (1 + γ 5 )/ 2 are the projectors of chirality. Also, u and d are column vectors subsuming the fields of the physical up-type and down-type quarks, respectively. In equation (15) we have omitted the Yukawa interactions of the Goldstone bosons; they have the same Lagrangian as in the SM. In this paper we put forward a 2HDM where CP is conserved in the dimension-four terms of the Lagrangian, hence θ QCD = 0, and det (M p M n ) is real both because of CP and because of a Z 3 symmetry. (Another symmetry-furnished 2HDM that also purported to alleviate the strong CP problem was proposed long time ago [6].) Remarkably, our model provides for the absence of strong CP violation even at the one-loop level.
In section 2 we explain our model. In section 3 we demonstrate that strong CP violation vanishes at the one-loop level in our model. Section 4 is devoted to the phenomenological constraints on the model. Section 5 contains the main conclusions of our work. Appendices A, B, and C present some formulas used in the analysis of section 4; appendix D gives a benchmark point for the parameters of the model.

The model
Our model is a 2HDM supplemented by the standard CP symmetry and by a Z 3 symmetry. Let ω = exp (2iπ/3), then the Z 3 symmetry reads This represents just a slight change from the Z 3 symmetry of the 2HDM of ref. [7]. Both CP and Z 3 are softly broken by terms in the quadratic part of the scalar potential.
Soft breaking of a symmetry consists in that symmetry holding in all the Lagrangian terms of dimension higher than some value, but not holding for the Lagrangian terms of dimension smaller than, or equal to, that value. In our case, both CP and Z 3 hold for terms of dimension four but are broken by terms of dimension two, viz. by the terms with coefficient µ 3 in line (17a). In principle, a model with a softly broken symmetry should eventually be justified through an ultraviolet completion, viz. a more complete model, with extra fields active at higher energies, which effectively mimics at low energy scales the model with the softly-broken symmetry. Unfortunately, such a ultraviolet completion is often quite difficult to construct explicitly -we attempted such a construction by adding singlet scalar fields to the theory, who would develop vevs at some high scale and then be "integrated out", leaving the desired low energy potential with only two doublets. However, we were unable to build such extensions that left both the CP and Z 3 symmetries intact at low energies -which of course does not mean such an UV completion does not exist. In the absence of any such explicit construction, a softly broken symmetry constitutes a strong, non-trivial assumption. This is, certainly, a weakness of the model in this paper.
The softly broken CP and Z 3 scalar potential is written as where µ 3 is real and positive by definition. The terms with coefficient µ 3 break the symmetry Z 3 softly. The phase ℵ breaks CP softly.
The vacuum expectation values 0 |φ 0 Thus, there is one gauge-invariant vacuum phase that offsets one phase in the potential, with the consequence that the potential of the physical scalar fields is CP -invariant, in particular there is no mixing between the scalars h and H and the pseudoscalar A. The stationarity equations for the vacuum are Referring to equation (4) and defining for a = 1, 2, the potential is then The potential (20) is invariant under the CP transformation CP : for a = 1, 2, where x = (t, r) andx = (t, − r). The phase λ in the CP transformation (21) is arbitrary.
The physical potential contains seven parameters µ 1,2,3 and λ 1,2,3,4 , since the phase ℵ in line (17a) is cancelled out by the vacuum phase. Instead of those seven parameters we will use as input v = 246 GeV, m h = 125 GeV, the angles α and β, and the masses m H of H, m A of A, and m H + of H ± . Then [8], In order for the potential to be bounded from below, one must impose the conditions [5] In order to avoid the situation of 'panic vacuum' [8] one must enforce the condition [9] The conditions in order for tree-level unitarity not to be violated are Because of the Z 3 symmetry (16), the matrices Γ a and ∆ a are The dimensionless numbers b a , d a , f a , p a , q a , and r a (a = 1, 2) are real because of the CP symmetry. Clearly, is real, hence θ QFD = 0. Because of the assumed CP invariance of the quartic part of the Lagrangian, θ QCD = 0 too. Thus, θ = θ QCD + θ QFD = 0, i.e. there is no strong CP violation at the tree level. We now define for a = 1, 2; and All the phases in equations (29) and (30) are either 0 or π because b a , d a , . . . , r a are real. We define the diagonal matrices where ℵ = ℵ 2 − ℵ 1 . We then have The matrices M n and M p are real, therefore they may be bi-diagonalized through real orthogonal matrices O Ln , O Rn , O Lp , and O Rp as Therefore, in the notation of equations (11), The CKM matrix is then One sees that the CKM matrix is complex because of the presence of the phase 3ℵ. When we compute the matrices N p and N n defined in equations (13), we find that and then, from equations (14), The matrices N d and N u are real.
Thus, in our model 1. The CKM matrix is complex.
2. The matrices N u and N d are real.
3. There is one pseudoscalar A that does not mix with the scalars h and H.

4.
There is no CP violation in the cubic and quartic interactions of the scalars.
In our model CP violation is located solely in the CKM matrix and originates entirely in the phase 3ℵ. This is the same that happened in the model of ref. [7]; however, in that model there was strong CP violation, while in the present model strong CP violation is absent at the tree level.

No strong CP violation at the one-loop level
At one-loop level the diagonal and real quark mass matrices M q (where q may be either u or d) get corrected by self energy diagrams: M q → M q + Σ q . If the diagonal elements of Σ q are complex, then Im tr M −1 q Σ q may be nonzero and strong CP violation may arise. In our model there are no complex phases except in the CKM matrix. Since the matrices N q are real, and since the scalars h and H do not mix with the pseudoscalar A, the Σ q generated through the emission and reabsorption (E&R) by the quarks of either h or H or A are real, hence innocuous. The same happens with the Σ q generated through the E&R of Z 0 gauge bosons. On the other hand, diagrams with the E&R of W ± gauge bosons do not generate mass renormalization (they just produce wavefunction renormalization), since the coupling of W ± to the quarks is purely left-handed. Therefore, the only diagrams where the complex matrix V arises, and might produce complex Σ q , are the ones with E&R of charged scalars H ± .
The Yukawa interactions of the charged scalars are given by lines (15g) and (15h). They contain two complex matrices, The one-loop self-energy of an up-type quark u α caused by the E&R of H + and a down-type quark d j is and we perform the computation in a space-time of dimension d. Thus, the only potentially complex part of the self-energy is where γ is Euler-Mascheroni's constant. One must sum the expression in the right-hand side of equation (40) over the flavour j of the quark d j .
The one-loop value of the strong-CP parameter θ is Im tr One easily finds that, in our model, both matrices M n M † n and M n N † n M −1 p † N † p , and all of their products too, have a structure of phases of the form where p † N † p , and hence its trace, are real, no matter what the function f is.
In this way we have demonstrated that tr (M −1 u Σ u ) is real. In a similar way one may show that tr M −1 d Σ d is also real, hence strong CP violation vanishes at the one-loop level in our model. 4 Phenomenological analysis of the model 4

.1 Constraints
We proceed to analyse how our model conforms to the experimental results. The model has tree-level flavour-changing neutral currents (FCNC) coupling to the scalars, so there is a wealth of flavour-physics observables that need to be taken into account whilst performing a fit of the model to the experimental data. Our procedure involves a global fit of the model's parameters, simultaneously requiring compliance with the theoretical and experimental bounds from the gauge, scalar, and fermionic sectors.
One may rotate the right-handed quarks n R1 and n R2 between themselves in such a way that the entry f 2 of the Yukawa-coupling matrix Γ 2 becomes zero. Similarly, one may rotate p R1 and p R2 so that q 2 becomes zero. 4,5 We use as input the ten entries b 1 , d 1 , f 1 , b 2 , d 2 , p 1 , q 1 , r 1 , p 2 , and r 2 (f 2 and q 2 are set to zero) of the Yukawa-coupling matrices (26), allowing those entries to be either positive or negative. We further input the CP -violating phase ℵ. We fit these eleven parameters in order to reproduce the quark masses [11] and the CKM-matrix observables [11] |V us | = 0.2243 ± 0.0005, (45a) We furthermore input tan β = v 2 /v 1 . We compute the matrices N n and N p through equations (36), and the matrices N d and N u through equations (37). Finally, we input α and get to know L physical in equation (15). 7 In the scalar potential (17) there are seven independent parameters µ 1,2,3 and λ 1,2,3,4 . We input instead the seven observables GeV, which produces the correct masses for the electroweak gauge bosons W ± and Z 0 ; 2. the lightest CP -even-scalar mass m h = 125 GeV, corresponding to the Higgs boson observed at the LHC; 4 Notice that n R1 and n R2 transform in the same way under the Z 3 symmetry (16), and p R1 and p R2 also transform in the same way under that symmetry. 5 With f 2 = q 2 = 0, the phases ∆ϕ and ∆χ in equations (30) become meaningless. That has no impact on our reasonings, in particular the matrix V in equation (35) does not depend on those phases. 6 We have doubled the uncertainty intervals quoted in ref. [11] for the masses of the light quarks u, d, and s; we have done this because of the large theoretical indefinition, due to QCD considerations, as to what exactly should be interpreted as the value of those masses. 7 The angle α may be restricted to lie either in the first quadrant or in the fourth quadrant [5].
3. the angle β = arctan (v 2 /v 1 ); 4. the angle α; 5. the remaining scalar masses-m H of the second CP-even scalar, m A of the pseudoscalar, and m H + of the charged scalar.
The last five parameters must be found through the fitting procedure. We have constrained m H + to be larger than 100 GeV and m H and m A to be larger than 130 GeV. We have furthermore assumed all three masses to be smaller than 1.2 TeV; values of the masses larger than 1.2 TeV would certainly be allowed by the fitting procedure.
The quartic couplings of the model are determined via equations (22). We check that the scalar potential is bounded from below, that it does not have a panic vacuum, and that it satisfies unitarity, viz. we check conditions (23, 24, 25). The constraints from the electroweak oblique parameters S and T are also imposed, by using the expressions for the 2HDM in refs. [12,13].
We will now go into detail about the further constraints that we have imposed.
• We implement the b → sγ bound described in appendix A, including the contributions from both the neutral and the charged scalars.
• The most relevant bounds on the off-diagonal entries of the matrices N d and N u come from flavour-physics observables, specifically the K, B d , B s , and D neutral-meson mass differences, and the CP -violating parameter K . We detail the computation of those quantities, and the requirements on them that we use in our fit, in appendix B.
• The Z → bb constraints described in appendix C are also taken into account. In this case we use only the charged-scalar contributions; the neutral-scalar ones should be negligible.
• For the regions of parameter space where m t > m q + m H + , q being a down-type quark, or where m t > m q +m S , q being either c or u and S being either h or H or A, we require that the branching ratio for each of the kinematically viable t → light quark + scalar decays be smaller than 5 × 10 −3 , in accordance with the current results on FCNC top decays and on the total top-quark width [11].
• In order that the scalar h of our model complies with the observational data from the LHC-it should be SM-like in its behaviour-we require that its couplings to the electroweak gauge bosons and to the top and bottom quarks do not deviate significantly from the SM expectations. We achieve this by focusing on the coupling modifiers κ X defined as cf. equations (15c) and (15d). In the first stage of the fit we constrain these couplings to obey 0.8 ≤ κ V ≤ 1, 0.8 ≤ κ t ≤ 1.2, and 0.8 ≤ |κ b | ≤ 1.2, 8,9 in order to roughly reproduce the LHC results. A second stage of the analysis further constrains these couplings, as detailed below.
A numerical scan of the parameter space of the model, in both the scalar and Yukawa sectors, was performed to discover points that obey all the constraints described above. It must be stressed that we have introduced nowhere in our scan a 'no-fine-tuning' assumption: we have tolerated any set of input values that led to the right outputs, even if either the input values or any intermediate computations displayed either 'fine-tunings' or 'unnatural cancelations'. Strong fine-tunings are often required in order to fit the D-meson mass difference constraint whenever m H and m A are not very high; for the other constraints, fine-tunings are at most moderate and do not occur at all for many points in our fit. As a matter of fact, even for the D-meson mass difference, there are many choices of parameters for which no fine-tuning is necessary and one of the scalars has relatively low mass; one such case is presented in appendix D.
With those points we have proceeded to compute the LHC production cross sections of the neutral scalars in the model, using the software SusHI [16,17] to include the NNLO QCD corrections. We have limited ourselves to the gluon-gluon production process, which is the dominant one in the LHC environment. Regarding the vector boson-fusion process, no differences will occur in this model vis a vis the usual 2HDM, as the couplings of the scalars to the gauge bosons are the same in both models.
The results for h are expressed in terms of the ratios where X may be either Z 0 Z 0 , W + W − , bb, ττ , or γγ. The value µ X = 1 indicates exact SM-like behaviour. We require that all the µ X be within 20% of 1, which is a fair description of the current LHC results, taking into account the uncertainties. With this imposition, the ranges of variation of κ V = s β−α , κ t , and κ b become much smaller than initially allowed in the fit: we obtain 0.929 ≤ κ V ≤ 1 and 0.952 ≤ {κ t , κ b } ≤ 1.04. For comparison, we will also present results for the tighter constraint |µ X − 1| < 0.1.
In principle, we should also consider the leptons. The Z 3 symmetry in the quark sector must be extended to the leptonic sector. Since flavour violation with leptons is much more constrained than with quarks, the best choice would be to extend Z 3 to the leptonic sector in a way identical to the flavour-preserving 2HDMs, allowing only one of the two doublets Φ a to couple to the leptons and give them mass. We would then have have two possibilities for the couplings of the scalars to the charged leptons-either Φ 1 couples to the charged leptons or Φ 2 does. The coupling modifier κ τ = g hττ g SM hττ is given by for the first and second choices, respectively. For definiteness, in our fit we have adopted the second option in equation (48), viz. we have imposed to the points in our fit. However, since the extension of our model to the leptonic sector is largely arbitrary, we have refrained from taking into account any other constraints on our model that might arise from processes involving leptons. We point out, though, that flavour-changing constraints from processes like K L → µ + µ − or B s → µ + µ − may pose serious challenges to our model.

General results
The bounds on the scalar sector-unitarity, oblique parameters, and vacuum stabilityproduce the same contraints on the model's parameters than those found in the usual version of the 2HDM. On the other hand, since the symmetry that we are considering affects in a non-trivial way the quark Yukawa matrices, there are major differences relative to other 2HDMs when the flavour-physics bounds are imposed. The flavour constraints from meson observables and from the top-quark FCNC decays and total width, previously described, constrain severely the magnitudes of the off-diagonal elements of the matrices N u and N d .
Our fit achieves to keep FCNC under control even with extra scalars of "low" masses-H and A may have masses below 500 GeV. This is in contrast with the often-made assumption that models with tree-level FCNC imply masses above 1 TeV; this had already been shown not to necessarily apply in the previous version of the current model [7]. In ref. [18] it has been argued that contributions from the scalar and pseudoscalar particles (H and A, respectively) to FCNC meson observables tend to cancel each other; we have explicitly observed that, for many points in our fit (cf. the point given in appendix D), the arguments of ref. [18] apply, and this is the reason why masses of the extra scalars lower than 1 TeV are possible. One consequence of the present model is the fact that the scalar h, which we have taken to be the 125 GeV state observed at the LHC, has tree-level FCNCs, as indicated by its Yukawa interactions in lines (15c) and (15d). Thus, unlike in the (tree-level) SM, h has the possibility of FCNC decays to final states sb, db, ds, uc, and their charge-conjugate states. However, for all the points resulting from our fit, these decays are extremely suppressedthe sum of the branching ratios for all of them being at most 2 × 10 −6 but usually much lower. Therefore, in our model the FCNC decays of the 125 GeV-particle are impossible to observe at the LHC, and almost certainly even at future e + e − colliders such as the ILC; the existence of those decays has no measurable impact on the phenomenology of the scalar h. The FCNC also raise the possibility of alternative production mechanisms for h, such as ds → h or uc → h; such production channels would be favoured by larger proton PDFs relatively to the SM production mechanism bb → h. However, once again these FCNC processes are found to be extremely small in our fit. According to (15c) and (15d), all the FCNC h interactions are suppressed by their proportionality to c β−α , which is required to be quite small by the SM-like behaviour of the 125 GeV scalar h, by the ratio between a light-quark mass and v = 246 GeV, and, sometimes, by the smallish off-diagonal N d and N u matrix elements induced by compliance with meson-physics bounds. In fig. 1 we show the points generated by our fit, displayed on the tan β-m H + plane. (The observed low density of points is merely a consequence of the difficulty in achieving good fits-further searches would yield more points and fill many more regions in the plot; lack of points in some areas has no physical meaning, it is just an artifact of the limited parameter space scan.) A clear conclusion from fig. 1 is that 1/20 < tan β < 20 in our model 10 ; this is mainly a consequence of the b → sγ bounds and, to a lesser degree, of the Z → bb bounds. One also sees in fig. 1 that values of the charged-scalar mass as low as 130 GeV are easily attained; this is in stark contrast with the findings for the type II 2HDM, where a lower bound on m H + of roughly 580 GeV exists. Unlike in the usual type-I and type-II 2HDMs, in our model the quark mass matrices do not emerge from the Yukawa couplings to a single scalar doublet, but rather from the couplings to both Φ 1 and Φ 2 . As such, although we employ the standard definition tan β = v 2 /v 1 , the usual wisdom about the values of this parameter does not apply.
The matrices N u and N d also exist in the usual flavour-preserving 2HDMs, but there they are diagonal and proportional to the quark mass matrices. In fact, in the type-I 2HDM whereas in the type-II 2HDM It appears that, in this model, most regions of parameter space yield either approximate type-I behaviour or approximate type-II behaviour for bottom quarks. Note that, although the blue points appear superimposed on the green and red lines, they are not exactly on them-the type-I and type-II behaviours displayed are approximate and there are deviations from them, which indeed can be large, as we observe in particular for low values of tan β. In fig. 2 (b) we observe the same behaviour for top quarks-most points have either |(N u ) 33 |/ m t ≈ tan β or |(N u ) 33 |/ m t ≈ cot β. From equations (50) it follows that, in the type I 2HDM, is equal to one, whereas in the type II 2HDM, from equations (51), is equal to one. In fig. 3 we display the quantities (52) and (53) plotted against each other. That figure shows that, for most points, our model is more similar to the type-I 2HDM, at least in what concerns giving mass to the third generation, i.e. the top-quark and bottomquark masses originate mostly in the Yukawa couplings to the same scalar doublet. However, there are also many allowed points for which the quantity (52) is not unity; for those points another regularity applies, namely the quantity (53) is very close to zero.
Still, one should remember that in our model there is flavour violation in the Yukawa interactions, and one obtains different results from those in figs. 2 and 3 for both the (1, 1) and (2, 2) entries of both N d and N u . In fact, the deviations from either type I-or type IIlike behaviour for the first and second generations are much more pronounced than what one observes in fig. 2. But, the corresponding Yukawa couplings being much smaller, that has much less importance for the Higgs-boson phenomenology than the third-generation couplings that we have discussed in those figures.

Properties of the extra scalars
We now turn to the extra neutral scalars in the model, H and A. The LHC Collaborations have been looking for neutral scalars other than the 125 GeV boson by investigating the production of W + W − , Z 0 Z 0 , and ττ , among other channels. The non-observation thus far of meaningful excesses in the cross sections, relatively to their SM expectations, imposes bounds on the masses and couplings of new particles. In our model, the imposition of the top-and meson-physics constraints should force the off-diagonal entries of the matrices N d and N u to be smallish, hence we expect that, just as h, the scalars H and A will have reduced flavour-changing interactions. Therefore, our model is expected to behave very much like the flavour-preserving 2HDMs in what concerns the possibility of evading the current experimental non-observation bounds for the extra scalars. As we will now show, there is a vast parameter space still allowed by the experimental constraints.
We firstly consider the limits coming from the search for resonant Z 0 Z 0 pairs by both the ATLAS [19,20,21,22,23] and CMS [24,25,26] Collaborations. This is a good channel to look for the heavy CP-even scalar H, which may decay at tree level as H → Z 0 Z 0 . 11 In fig. 4 we show the points which obey all the constraints described in section 4.1; the points in blue correspond to the requirement that all the µ X are within 20% of 1, and the points in red have all the µ X less than 10% away from 1. The yellow line is the upper 2σ bound from the observed limit from ref. [23]. We observe that most of the allowed parameter space yields a pp → gg → H → Z 0 Z 0 cross section below the experimental upper bound; only a few low-m H points exceed the bound, but even for those low values of m H there are plenty of points which are still allowed. The tighter constraint of 10% on the h production rates does not qualitatively change the picture. There is a simple explanation for why low values of the pp → gg → H → Z 0 Z 0 event rate should be obtained, namely, in any CP -conserving 2HDM (or, indeed, multi-Higgs-doublet model) there is the sum rule for the couplings of the CP -even neutral scalars to gauge-boson pairs. Therefore, if the coupling of h to Z 0 (and W ± ) pairs is very close to its SM value, then the coupling of H to such pairs will be suppressed. Equation (54) is normally expressed through g 2HDM hZZ = s β−α g SM hZZ and g 2HDM HZZ = c β−α g SM hZZ ; SM-like behaviour of h means s β−α 1, which implies c β−α 0.
In fig. 5 (a) we show the gluon-gluon production cross section for a pseudoscalar A, multiplied by its branching ratio to a tt pair, at LHC. The gluon-gluon production and decay to hh of the scalar H versus its mass. All the points displayed obey the constraints described in subsection 4. For the blue points, h has production rates within 20% of its SM-expected values; for the red points those production rates are within 10% of the SM value. The green points are a subset of the red ones, for which the width of the scalar in each plot is larger than 10% its mass. The yellow line is the 2σ upper bound in figure 6 of ref. [29].
2HDM. In fig. 5 (a), the yellow line is the upper 2σ bound in figure 11 of ref.
[28]. 12 The published results only extend down to m A 500 GeV, but it is clear that no exclusion will occur even for A masses lower than that. As before, the red (blue) points indicate a cut of 10% (20%) on the µ X ratios for the Higgs boson h, meant to ensure its SM-like behaviour and compliance with the LHC results. The green points in the same plot are the subset of the red ones for which the width of A is larger than 10% of its mass: Γ A / m A > 0.1. We have thus far been assuming the validity of the narrow-width approximation and neglecting eventual interferences between backgrounds and signal; by marking these large-width points in green, we want to draw attention to the only regions where that approximation might fail. 13 The conclusion to draw from fig. 5 (a) is that the current exclusion bounds from the tt resonance searches are easily evaded by our model. In fig. 5 (b) we investigate the possibility of the heavy CP -even scalar H being observed through its decay to two 125 GeV scalars h. This hh channel is being thoroughly studied at the LHC, considering several possible decay channels for both h particles [29, 30, 31, 32, 33, 34, 35,36,37]; the yellow line in fig. 5 (b) is the 2σ upper bound of figure 6 of ref. [29]. The blue and red points are the same as before; the green points are the subset of the red ones for which Γ H / m H ≥ 0.1. (Therefore, the green points in fig. 5  green points in fig. 5 (a).) Unlike in fig. 5 (a), the green points, corresponding to scalars H with a large width, 14 correspond to smaller values of σ × BR. Just as in the previous figures, we see that virtually all of our parameter space, except a few low-mass points, complies with the existing experimental bounds.
Since the model that we are studying differs from usual versions of the 2HDM through the existence of FCNC, we have considered the possibility of single-top decays of the heavy (pseudo)scalars H and A. Indeed, the non-diagonal Yukawa interactions lead to the possibility of decays like H → tū and A → ct, which might be observed as top-quark + jet events at the LHC; such events should be quite challenging to study in an hadronic machine such as the LHC, but the recent progress in charmed-jet identification algorithms may be a significant contribution for a future analysis. 15 . In fig. 6 we present the expected cross section times branching ratio for both (a) the H and (b) the A. In fig. 6 we have grouped together all the FCNC decays of the scalars with a single top in the final state, viz. the decays to tū,tu, tc, andtc. The top-and meson-physics constraints described in subsection 4.1 usually produce N u matrices with smallish off-diagonal elements, and this yields very small branching ratios for FCNC decays of both H and A-the maximum values that we have obtained were smaller than 5 × 10 −4 , but usual values were much smaller than that. Therefore, the model predicts values for σ × BR usually orders of magnitude below the fentobarn. It is difficult to find experimental bounds on such a search channel, but the search for a W decaying to a single top quark plus a bottom quark [38] can at least give a rough idea of the current sensitivity of the LHC for a top + jet resonance analysis. Though the mass range is different (the analysis of ref. [38] starts at 1 TeV), the bounds shown in that paper for the cross section times the branching ratio are of order 0.1 pb, and therefore much above the predicted σ × BR shown for our model in fig. 6.

Conclusions
In this paper we have presented a two-Higgs-doublet model that attempts a partial solution of the strong CP problem by relegating a possible generation of a nonzero θ to the two-loop level. Our model achieves this by postulating a soft CP violation that transfers itself just to the CKM matrix, with no CP violation anywhere else in the model, especially no CP violation in scalar-pseudoscalar mixing.
We do not claim that our model achieves a full solution of the strong CP problem, because the θ generated at two-loop level might still be too large. However, since in our model CP violation exists only in the CKM matrix, one may expect θ to be proportional to J ∼ 10 −5 , the only CP -violating invariant quantity in that matrix. Adding in a two-loop factor (16π 2 ) −2 ∼ 10 −4 and [39] probable suppression factors m q /m W , where m q is a generic second-generation quark mass, one might well reach a sufficiently small θ.
Of course, a 2HDM where CP violation only occurs in the CKM matrix eschews one of the motivations for multi-Higgs-doublet models, namely, obtaining extra sources of CP violation in order to reach a sufficiently large baryon number of the Universe. We have nothing to say about this insufficiency.
We have investigated the compatibility of our model with the outstanding experimental constraints, in particular on the flavour-changing neutral currents. Our model can easily evade them, at the price of cancelations that might be qualified as fine-tuning. We find that the new scalars in our model may in some cases be little heavier than the observed Higgs particle of mass 125 GeV. That will not necessarily make them easy to discover, though, as we have seen in section 4. A The decayB → X s γ The decays of a bottom-flavoured meson to a strange-flavoured meson and a photon proceed via the quark transition b → sγ. Those decays constitute one of the most relevant constraints on the parameter space of a multi-Higgs-doublet model, because they receive important contributions from loops with charged scalars. This is because the interactions of a charged scalar with down-type quarks may be substantially enhanced by ratios of VEVs. Consequently, in a 2HDM the constraints from b → sγ typically eliminate substantial regions of the m H + -tan β plane.
In the model under discussion in this paper, the occurrence of tree-level FCNC means that the neutral scalars also contribute to b → sγ, unlike what happens in flavour-conserving 2HDMs. We follow the general analysis of ref. [40] to take into account both the charged and the neutral scalars' contributions. We write the Yukawa interactions of our model in the notation with coefficients defined as 16 The Wilson coefficients required for the computation of b → sγ are [40] and The functions I 3,4,5,6 are given in equations (40)-(43) of ref. [40].
To compute the overall branching ratio of b → sγ we follow refs. [41,42]. We use the effective operators described above, defined as being at the Fermi scale µ = m W , and include NLO QCD corrections by choosing m b as the renormalization scale. Let η = α S (m W ) /α S (m b ) = 0.5651 [41] be the ratio of the running strong coupling constant between scales m W and m b . We compute and then [43] BR (b → sγ) = BR (b → sγ) SM (61a) where BR (b → sγ) SM = 3.15 × 10 −4 . We consider our model to be in compliance with the b → sγ data if it yields a branching ratio within twice the experimental error bar, viz. we require 2.4406 × 10 −4 < BR (b → sγ) < 3.8594 × 10 −4 .

B The neutral meson-antimeson observables
The FCNC induced by the off-diagonal entries of N d and N u lead to tree-level contributions to flavour observables such as CP violation through the parameter K and the mass differences in the K 0 , B 0 d , B 0 s , and D 0 meson-antimeson systems. These are sensitive observables and new-physics contributions to them may easily be overwhelming. Thus, we must make sure that the contributions to them from the scalar sector of our model conform to the current data. We use the numbers listed in ref. [11].
B.1 K 0 -K 0 observables Two K 0 meson observables are sensitive to the tree-level FCNC contributions from the scalar sector: the CP -violating parameter K and the mass difference between K S and K L . Both There are also contributions to the mass difference in the meson system D 0 -D 0 . As in the K 0 -K 0 system, there are considerable uncertainties in the calculation of M SM 21 . Therefore, once again, we resort to requiring only the New Physics contribution not to be too large. We have with f D = 0.212 GeV and m D = 1.865 GeV. Conservatively, we require 2 M NP 21 to be smaller than the measured mass difference 6.253 × 10 −15 GeV.

C Z → bb constraints
A potentially very important constraint to two-Higgs-doublet models (2HDMs) stems from the measurement of the decay Z → bb. We follow the treatment of that decay in refs. [45,46,47]. The Lagrangian for the Zbb vertex is written as where the coefficientsḡ L,R b are, at tree level in the SM,ḡ L b = −1/2 + s 2 W /3 andḡ R b = s 2 W /3. In both the SM and in extensions thereof, these coefficients get one-loop contributions. To wit, in the 2HDM the contributions of loops with charged scalars toḡ L b andḡ R b are given by where The contributions of loops with neutral scalars are expected to be small, both for 2HDMs with flavour conservation [45,46,47] or without it [48]; we neglect them. In order to take into account the current experimental results on the observable quantities R b and A b (see refs. [45,11]), we have required that the charged-scalar contribution added to the SM one, viz.ḡ L b = −0.42112 + δḡ L b andḡ R b = 0.07744 + δḡ R b , does not deviate from the SM prediction by more than 2σ, viz. 2 ḡ L b 2 + 2 ḡ R b 2 = 0.36782 ± 0.00143.

D A benchmark point
To illustrate the model, we provide a specific point, which is meant to serve only as an example. The input in the scalar sector is (72c) p 2 = −1.2295, q 2 = 0, r 2 = −9.2531 × 10 −3 . (72d) We also input ℵ 1 = 0, ℵ 2 = 1.33 rad. One thus obtains quark masses and a CKM matrix in agreement with equations (44,45). The matrices that parameterize the FCNC are One sees that some off-diagonal matrix elements of N d are not very small, and some offdiagonal matrix elements of N u -which is almost a triangular matrix-are pretty large.
There is a cancelation of about one part in 47 between the contributions to the neutral-Dmeson mass difference of the neutral scalars h and H on the one hand and of the pseudoscalar A on the other hand; there are analogous, yet milder, cancelations among the contributions to the other neutral-meson mass differences. In general, we have found that the D-meson mass difference constraint requires quite strongly fine-tuned cancelations when the neutral scalars have low masses, while the constraints from all other neutral-meson systems are much easier to satisfy and mostly require no fine-tuning. Still, notice that this benchmark point has one particle (the pseudoscalar A) with relatively low mass. With this benchmark point, the coupling modifiers defined in equations (46) and (48)    [30] ATLAS Collaboration, Search for Higgs boson pair production in the bbγγ final state using pp collision data at √ s = 13 TeV with the ATLAS detector, ATLAS-CONF-2016-004.
[31] CMS Collaboration, Search for Higgs boson pair production in the final state containing two photons and two bottom quarks in proton-proton collisions at √ s = 13 TeV, CMS-PAS-HIG-17-008.  [34] ATLAS Collaboration, Search for pair production of Higgs bosons in the bbbb final state using proton-proton collisions at √ s = 13 TeV with the ATLAS detector, ATLAS-CONF-2016-049.