Decays A → Zγγ and φ → Zγγ ( φ = h , H ) in two-Higgs doublet models

The one-loop contributions to the decays of the CP -odd and CP -even scalar bosons A→ Zγγ and φ → Zγγ (φ = h,H) are calculated within the framework of CP -conserving THDMs, where they are induced by box and reducible Feynman diagrams. The behavior of the corresponding branching ratios are then analyzed within the type-II THDM in a region of the parameter space around the alignment limit and still consistent with experimental data. It is found that the A→ Zγγ branching ratio is only relevant when mA > mH +mZ , but it is negligible otherwise. For mA > 600 GeV and tβ ' O(1), BR(A → Zγγ) can reach values of the order of 10−5 − 10−4, but it decreases by about one order of magnitude as tβ increases up to 10. A similar behavior is followed by the H → Zγγ decay, which only has a non-negligible branching ratio when mH > mA + mZ and can reach the level of 10−4 − 10−3 for mH > 600 GeV and tβ ' O(1). We also estimated the branching ratios of these rare decays in the type-I THDM, where they can be about one order of magnitude larger than in type-II THDM. As far as the h → Zγγ decay is concerned, since the properties of this scalar boson must be nearly identical to those of the SM Higgs boson, the h→ Zγγ branching ratio does not deviates significantly from the SM prediction, where it is negligibly small, of the order of 10−9. This result is in agreement with previous calculations.


I. INTRODUCTION
The standard model (SM) has provided a successful description of the observed electroweak phenomena at the energy scales explored until now, as confirmed recently with the discovery of the Higgs boson by the ATLAS and CMS experiments at the CERN LHC [1,2].Nonetheless, it is worth to explore whether there is a unique Higgs boson, as predicted by the SM, or the electroweak symmetry breaking (EWSB) mechanism requires additional Higgs bosons.To address some SM flaws, a plethora of extension models have been proposed, several of which contain a scalar sector with more than one Higgs multiplet, thereby predicting more than one physical Higgs boson.If experimental data reveal the existence of any additional Higgs bosons, it will be crucial to test what extension model is consistent with such particles.The simplest of such theories are two-Higgs doublet models (THDMs) [3,4], which are obtained by adding a second complex SU (2) L Higgs doublet to the SM one.These models respect the ρ = 1 relation at the tree-level, contrary to other higher-dimensional Higgs-multiplet models.Also, in spite of its simplicity, THDMs can predict several new phenomena absent in the SM, such as new sources of CP violation, tree-level scalar-mediated flavor changing neutral currents (FCNCs), a dark matter candidate, etc.After EWSB, three of the eight degrees of freedom are removed from the spectrum to provide the longitudinal modes of the W ± and Z gauge bosons.Five physical Higgs bosons remain as remnant: a charged Higgs boson pair H ± and three neutral Higgs bosons h, H, and A. If the scalar sector respects CP invariance, the neutral scalar bosons are CP -eigenstates: h and H are CP -even, whereas A is CP -odd.It is usually assumed that one of the neutral CP -even scalar bosons is the one observed at the LHC.The most general CP -conserving THDMs have tree-level FCNCs [5], which can be removed by imposing a Z 2 discrete symmetry that forbids such interactions at the tree-level [6].In this scenario, there are four THDM types, which are typically known as type-I, type-II, lepton-specific [7] and flipped THDM [8].It turns out that type-II THDM is the most studied in the literature as it has the same Yukawa couplings as the minimal supersymmetric standard model (MSSM), therefore its still-allowed region of parameter space has been considerably studied.
Since the proposal of the Higgs mechanism, the phenomenology of the Higgs bosons has been the focus of considerable attention.As for the dominant tree-level decay modes of a CP -even Higgs boson h → f f and h → V V (V = W, Z), they have been long studied in the literature both in the SM and several of its extensions, along with the one-loop induced decays h → γγ, h → γZ, and h → gg.Although the h → γγ decay has a tiny branching ratio for a 125 GeV Higgs boson, it was very helpful for the detection of the SM Higgs boson.This decay mode has the advantage of a relatively low background, so it was fundamental in the design of the ATLAS and CMS detectors.As for the h → gg decay, it is undetectable but it is fundamental to compute the cross section for Higgs production via gluon fusion.
It is expected that the data collected at the LHC may allow us to search for any other rare decays of the Higgs boson [9], such as lepton flavor changing Higgs decays h → ¯ i j (i = j) or invisible Higgs decays h → E T , which are forbidden in the SM and can shed light on any new physics effect.Even more, with the prospect of a future Higgs boson factory, other exotic decays of the Higgs boson could be at the reach of experimental detection.In particular, the rare decay h → Zγγ is very suppressed in the SM as it arises at the one-loop level via the exchange of charged particles, so it can offer a relatively clean signal of new physics: two energetic photons plus a back to back lepton anti-lepton pair.This process can also provide a test for the couplings of the Higgs boson to the particles running into the loops, which can be SM particles or any new charged particle predicted by other extension models.A similar decay is h → Zgg, which at the leading order can be straightforwardly calculated from the h → Zγγ one.In addition, the study of the hZγγ and hZgg vertices would allow us to obtain the leading order contributions to the cross section of hZ pair production via photon fusion γγ → hZ and gluon fusion gg → hZ [10,11].
On the other hand, a CP -odd scalar boson has fewer decay channels and so it is worth studying some one-loop induced decays of such a particle.At tree level, its dominant decay channels are A → f f , A → Zh(H) and A → W ± H ∓ , when kinematically allowed, whereas at the one-loop level a CP -odd scalar boson can decay as A → gg, A → γγ and A → Zγ [3].These decay channels can have significant branching ratios in some regions of the parameter space.Other one-loop induced decay modes such as A → W W and A → ZZ have already been studied in [12,13], though they are more suppressed than the aforementioned decay channels.
In this work we are interested in studying the A → Zγγ and φ → Zγγ (φ = h, H) decay modes in the context of THDMs, which induce these processes at the one-loop level via box and reducible Feynman diagrams, with contributions from charged fermions, mainly from the top and bottom quarks.The W gauge boson and the charged scalar boson H ± can only contribute through reducible diagrams to the A → Zγγ decay .The respective decay of the SM Higgs boson has already been studied: the decay h → Zγγ was studied in Ref. [14] and the analogue decay h → Zgg was studied in [15,16].To our knowledge, the A → Zγγ decay has not been studied until now.
The organization of this paper is as follows.Section II is devoted to a brief discussion of the general THDM, focusing on the CP -conserving THDMs.In Sec.III we present the details of the calculation of the decays A → Zγγ and φ → Zγγ (φ = h, H) by the Passarino-Veltman reduction scheme.We present the analytical expressions for the invariant amplitudes, the decay widths, as well as the kinematic distributions of the invariant mass of the photons and the energy of the Z gauge boson, which can be useful to disentangle the decay signal from its potential background.The numerical analysis of the branching ratios within type-II THDM is presented in Sec.IV, whereas the conclusions and outlook are presented in Sec.V.The Feynman rules necessary for the calculations and some lengthy formulas are presented in the appendices.

II. TWO-HIGGS DOUBLET MODELS
THDMs have been largely studied in the literature [3].We will present here a brief outline of CP -conserving THDMs, including only those details relevant for our calculation.For the interested reader, a comprehensive review of these models can be found in [4].

A. THDM Lagrangian
In THDMs, two complex SU (2) L Higgs doublets Φ i are introduced in the scalar sector: where v i are the vacuum expectation values (VEVs) of the neutral components, which satisfy where |D µ Φ i | 2 is the kinetic term for the two-Higgs doublets, with D µ the SM covariant derivative, V (Φ 1 , Φ 2 ) is the Higgs potential, L Y uk denotes the Yukawa interactions between Φ i and the SM fermions, and L SM describes the SU (2) L × U (1) Y interactions of fermions and gauge bosons.The most general gauge-invariant renormalizable potential V (Φ 1 , Φ 2 ) for THDMs is a hermitian combination of electroweak invariant combinations.It contains 14 parameters and can give rise to new sources of CP violation [17].However, as long as CP is conserved in the Higgs sector, the scalar potential for the two doublets Φ 1 and Φ 2 with hypercharge +1 can be written in terms of 8 parameters as follows [3,4] After EWSB, three of the eight degrees of freedom of the two Higgs doublets are the Goldstone bosons (G ± , ξ), which are absorbed as longitudinal components of the W ± and Z gauge bosons, whereas the remaining five degrees of freedom become the physical Higgs bosons: there is a pair of charged scalar bosons H ± , two neutral CP -even scalar bosons h and H, where m h < m H by convention, and one neutral CP -odd scalar A. Since all the parameters appearing in the potential are real, there are no bilinear mixing terms, which is why the neutral mass eigenstates are also CP eigenstates.In the neutral sector the following mass term appears with As far as the Yukawa Lagrangian L Y uk is concerned, the scalar-to-fermion couplings are not univocally determined by the gauge structure of the model.The most general Yukawa Lagrangian for THDMs is [4] where Φj = iτ 2 Φ j , Y f are 3 × 3 complex matrices, and the left-and right-handed fermion fields are tree-vectors in flavor space.
To prevent tree-level FCNCs it is usual to introduce a discrete Z 2 symmetry respected by the Φ i doublets and the fermions.Under this symmetry one of the scalar doublets is even Φ 2 → Φ 2 and the other one is odd Φ 1 → −Φ 1 .This gives rise to four types of THDMs, which are usually known as type-I THDM, type-II THDM, lepton-specific THDM and flipped THDM.The way in which each Higgs doublet couples to the fermions in these models is summarized in Table I.On the other hand, if no Z 2 discrete symmetry is imposed, there will be tree-level FCNCs.In such a scenario both doublets couple to the charged leptons and quarks.This model is known as type-III THDM [18,19].In this work, however, we are not interested in this realization of THDMs.There is another version of THDMs, known as type-III THDM, in which both Higgs doublets couple to the leptons and quarks simultaneously, thereby giving rise to tree-level FCNCs [3,4].
Although we will present a rather general calculation within flavor-conserving THDMs, the numerical analysis will be carried out in the context of the type-II THDM, which is by far the most studied THDM since it shares the same Yukawa interactions as the MSSM.The most distinctive difference between the type-II THDM and the MSSM is that the former does not have a strict upper bound on the mass of the lightest Higgs boson, which is an important feature of the latter.In addition, in THDMs the scalar boson self couplings are arbitrary and so is the mixing parameter α, which in the MSSM is given in terms of tan β and the scalar boson masses.
Our calculation is to be performed in the unitary gauge.The Feynman rules for THDMs can be obtained once the Lagrangian is expanded in terms of mass eigenstantes and can be found for instance in Refs.[3,4].We present those Feynman rules required by our calculation in Appendix A.

A. Kinematic conditions
We now turn to present the A → Zγγ and φ → Zγγ (φ = h, H,) decay widths.We first present the kinematics conditions, which are defined according to the following notation for the external 4-momenta The mass-shell conditions thus read p 2 = m 2 φ , k 2 3 = m 2 Z and k 2 1 = k 2 2 = 0. We now introduce the following Lorentz invariant quantities: These variables are not all independent as s 1 + s 2 + s = m 2 φ + m 2 Z by four-momentum conservation.In our calculation, we express all the scalar products between the four-momenta k 1 , k 2 and k 3 in terms of the Lorentz invariant variables s 1 , s 2 and s as well as the scaled variable µ Z = m 2 Z /m 2 φ .In addition, because of the transversality conditions obeyed by the gauge bosons, i.e., k 1 we drop from the invariant amplitudes any terms proportional to k µ 1 , k ν 2 , and k α 3 .All the above kinematic conditions probe useful to simplify the calculation.We now present the invariant amplitudes for the A → Zγγ and φ → Zγγ (φ = h, H) decays, which are induced at the one-loop level at the lowest order in perturbation theory.

B.
A → Zγγ decay invariant amplitude There are two sets of Feynman diagrams that induce this decay: box diagrams and reducible diagrams.Once the invariant amplitude for each Feynman diagram was written down in the unitary gauge, we used the Passarino-Veltman reduction scheme to solve the loop integrals [20], which were reduced down to a combination of two-, three-and fourpoint scalar functions.The algebra was carried out with the aid of the Mathematica package FeynCalc [21].We first present the invariant amplitude arising from the box diagrams.

Box diagram contribution
In Fig. 1 we show the box diagrams that contribute to the A → Zγγ decay.The dynamical content is rather simple in the sense that there is only one kind of particles circulating into the loop, namely, SM charged fermions.Other charged particles do not contribute to this decay at the one-loop level in THDMs: due to CP invariance in the scalar sector, the CP -odd scalar A does not couple to a pair of W gauge bosons or charged scalars H ∓ , though it can couple to a W ± H ∓ pair.However, the V W ∓ H ± vertex (V = γ, Z) is absent at the tree-level and so the A → Zγγ decay cannot proceed via box diagrams with both W ∓ and H ∓ particles.The main contributions of box diagrams are thus expected to arise from the heaviest fermions.For small t β , the top quark contribution would dominate, whereas for large t β the bottom quark contribution would become relevant.This is due to the presence of the factors 1/t β and t β appearing in the Yukawa couplings for the top and bottom quarks, respectively, as will be shown below.
FIG. 1: Box diagrams that contribute to the A → Zγγ decay in the THDM.There are three additional diagrams that are obtained by exchanging the photons.Similar diagrams also contribute to the φ → Zγγ (φ = h, H) decay, after the replacement A → φ.
Once the Passarino-Veltman reduction scheme was applied, we performed several test on our results.First, we verified that the invariant amplitude for all the box diagrams is gauge invariant under U (1) em , i.e. it vanishes when the photon four-momenta are replaced by their polarization vectors.We also verified that Bose symmetry is respected and that ultraviolet divergences cancel out.The invariant amplitude for the A → Zγγ decay can be cast in the following gauge-invariant manifest form with the Lorentz structures given as follows where the form factors F i depend on s 1 , s 2 , s, and µ Z , though we will refrain from writing out such a dependency explicitly.These form factors will receive contributions from both box and reducible diagrams, which means that the latter will not generate additional Lorentz structures.We can thus write , where the notation is self-explanatory.The expressions for the box diagram contributions are too lengthy and they are presented in Appendix B in terms of Passarino-Veltman scalar functions.

Reducible diagram contribution
There are also reducible diagrams in which the A → Zγγ decay proceeds as A → Zφ * → Zγγ (φ = h, H), as depicted in Fig. 2, with the two photons emerging from the intermediate scalar boson via loops carrying charged fermions, the W gauge boson, and the charged scalar boson + + Reducible Feynman diagrams for the A → Zγγ decay in the THDM.For the triangle diagrams there are additional diagrams that are obtained by exchanging the photons.Similar diagrams also contribute to the φ → Zγγ (φ = h, H) decay, except that the intermediate particle is now the CP -odd scalar boson A and there is only contribution from charged fermions in the triangle loop.
As was the case for the box diagram contribution, the reducible diagram contribution is gauge invariant and ultraviolet-finite by its own.It turns out that these diagrams contribute to the gauge-invariant amplitude of Eq. ( 17) only through the form factor F 1 , which includes the contributions of charged fermions, the W gauge boson, and the charged scalar boson with As for the φ → Zγγ (φ = h, H) decay, at the one-loop level it also receives the contributions of the fermion box diagrams of Fig. 1 with A replaced by φ.It is worth noting that although a CP -even scalar boson does couple to charged W ∓ gauge bosons and charged scalar bosons H ∓ , the corresponding box diagram contributions exactly cancel out due to CP invariance.Notice that the amplitude of this vertex must include the Levi-Civita tensor due to CP invariance, but it cannot arise via box diagrams with charged particles other than charged fermions, whose coupling with the Z gauge boson includes a γ 5 matrix.As the invariant amplitude of a fermion loop includes the trace of a chain of Dirac matrices, the term involving the γ 5 matrix would give rise to the required Levi-Civita tensor.The most general Lorentz structure for the φ → Zγγ (φ = h, H) decay can be written in the following gaugeinvariant manifest form where we use the shorthand notation αkpq = αβλρ k β p λ q ρ , etc. Again the form factors G i depend on s, s 1 , s 2 , and µ Z .To arrive to the above equation, we used Schouten's identity.These form factors receive contributions from both box diagrams and reducible diagrams: As far as the contributions from the box diagrams are concerned, they are reported in Appendix B in terms of Passarino-Veltman scalar functions.

Reducible diagram contribution
There are also contributions from reducible diagrams that are analogue to those depicted in Fig. 2, but with the photons emerging from the intermediate CP -odd scalar boson A via loops of charged fermions only.There are also extra reducible diagrams arising from the process φ → ZZ * → Zγγ, as shown in Fig. 3.This diagram involves the well-known triangle anomaly Z * γγ, which receives contributions from charged fermions only.This is due to CP invariance as the amplitude for this vertex must be proportional to the Levi-Civita tensor, which can only arise via the trace of a chain of Dirac matrices including γ 5 , which in turn only is present in a fermion loop.Therefore loops of the charged W gauge boson or the charged scalar boson do not contribute to this vertex.Also, due to the Landau-Yang theorem, the Z * γγ vertex vanishes for real Z, so this diagram does not contribute to the φ → Zγγ decay when the φ scalar boson is kinematically allowed to decay into a pair of real Z gauge bosons.These reducible Feynman diagrams only contribute to the invariant amplitude of the φ → Zγγ (φ = h, H) decay via the form factor G 3 : where G Z 3 and G A 3 are the form factors arising from the diagrams with the vertices Z * γγ and A * γγ, respectively.Explicit expressions in terms of Passarino-Veltman scalar functions are given in Appendix B. Note that we must include the contribution of all fermion families in order to cancel the Z * γγ anomaly.

D. A → Zγγ and φ → Zγγ (φ = h, H) decay widths
There are two scenarios for the φ i → Zγγ (φ = h, H, A) decays, which depend on the value of the mass of the incoming scalar boson φ i as compared to the mass of the exchanged scalar boson, which we denote by φ e : φ e = h, H for φ i = A or φ e = A for φ i = h, H.We will present the expression for the resulting decay width in both scenarios.
In this scenario, the incoming scalar boson φ i will not be heavy enough to produce an on-shell φ e in addition to the on-shell Z gauge boson.Therefore we will have a pure three-body decay induced by both box and reducible diagrams.The corresponding decay width can be written as where we introduced the following scaled variables with ŝ = s/m 2 φi and ŝi = s i /m 2 φi .In the center-of-mass frame of the decaying φ i we have x 1 = 2E Z /m φi , x 2 = 2E γ /m φi , and x 3 = 2E γ /m φi , where E γ (E γ ) stands for the energy of the photon with four-momentum k 1 (k 2 ).From energy conservation, these variables obey The kinematic limits in Eq. ( 21) are as follows The squared average amplitudes for both decays A → Zγγ and φ → Zγγ (φ = h, H) are presented in Appendix C.

m φ i > m φe + mZ
In this scenario the incoming scalar boson φ i is heavy enough to produce an on-shell scalar boson φ e .Therefore the φ i → Zγγ decay proceeds as the pure two-body decay φ i → Zφ e , followed by the decay φ e → γγ.Note that in the case of the decay of a CP -even Higgs boson, although the decay into a pair of real Z gauge bosons will now be kinematically allowed, the Z → γγ decay is forbidden by the Landau-Yang theorem, which means that the contribution of the intermediary Z gauge boson will thus vanish.In this scenario, using the Breit-Wigner propagator for the exchanged scalar boson, Eq. ( 21) can be integrated and the φ i → Zγγ decay width can be written as with the decay width Γ(φ e → γγ) given by For a CP -even scalar boson φ e = h, H, F φeγγ receive contributions from charged fermions, the charged W gauge boson, and the charged scalar boson H ± : with τ χ = 4m 2 χ /m 2 φe .The F φeγγ χ (x) functions can be obtained from the results for the reducible diagrams presented in Appendix B by setting s = m 2 φe .They are given by for φ e = h, H f (x) is given by On the other hand, when the intermediary scalar boson is the CP -odd one A, we only have the contributions of charged fermions As for the φ i → φ e Z decay width, it is given as follows Note that τ φe = 4m 2 φe /m 2 φi , thus τ φe = 4m 2 φ /m 2 A for A → Zγγ and τ φe = 4m 2 A /m 2 φ for φ → Zγγ (φ = h, H).A similar expression with the corresponding replacements is obeyed by the φ e → W ± H ∓ decays if kinematically allowed.
All the necessary coupling constants g φ f f , g φW W , g φZZ , g φAZ , g φW ± H ∓ , g φH − H + (φ = h, H), along with g A f f and g AW ± H ∓ are shown in Table II of Appendix A. Other coupling constants involved in decays such as H → hh and φ → AA can be found in Ref. [3,4] for instance.To obtain the branching ratio BR(φ e → γγ), we need the main decay widths of both CP -even and CP -odd scalar bosons, which have already been studied in the literature considerably [3].For completeness we present in Appendix D all the necessary formulas, which can also be helpful to obtain the branching ratio for the φ i → Zγγ decay in type-II THDM and make a comparison with that of other decay channels.

IV. NUMERICAL ANALYSIS AND RESULTS
We now turn to the numerical analysis.To begin with, we will analyze the current constraints on the parameter space of type-II THDM.
A. Allowed parameter space of type-II THDM After the Higgs boson discovery, several studies have been devoted to explore the implications on the parameter space of THDMs [22][23][24][25][26]. From the recent analyses of the ATLAS and CMS collaborations [27], it is inferred that the properties of the 125 GeV scalar boson found at the LHC are highly consistent with the SM predictions, thereby imposing strong constraints on the scalar sector of SM extensions.If one of such theories predicts several CP -even physical Higgs bosons, one of them must correspond to the SM one and reproduce its couplings to fermions and gauge bosons.In type-II THDM, the scalar boson h is usually assumed to be the lightest one and so is identified with the SM Higgs boson, which constrains the parameter space of the model to a region very close to the alignment limit sin(β − α) = 1, where the heavy Higgs H does not couple to the gauge bosons and the coupling hZA is absent at tree-level [26,28,29].The couplings of the h Higgs boson to the fermions involve the mixing angles α and β.Therefore, the LHC data can impose strong constraints on both parameters.Other constraints can be obtained from theoretical requirements such as vacuum stability and unitarity of the scalar potential as well as perturbativity of the Higgs couplings.Also, the oblique parameters S, T and U can impose strong constraints on the masses of the new Higgs bosons A and H, requiring that at least one of them is very heavy: a CP -odd scalar with m A ∼ 200 GeV requires a heavy CP -even scalar with m H ≥ 600 GeV and viceversa.As for the charged Higgs boson mass, it can be constrained through experimental measurements on low energy FCNC processes.
All of the above constraints can be complemented with the direct searches of additional Higgs bosons at LEP and the LHC.Below we present the constraints most relevant for our numerical analysis.
• Mixing angles β and β − α: since the h Higgs boson is identified with the SM Higgs boson, the LHC data restrict β − α to lie very close to π/2, namely, | sin(β − α)| > 0.999, with a small interval around t β = 1 where such a constraint is less stringent.Furthermore, in type-II THDM, for β − α π/2, the scalar couplings to the top quark (bottom quark) behaves as 1/t β (t β ), thus FCNCs process are very sensitive to small and large values of t β , which will impose stringent constraints on this parameter.We can thus consider values of t β in the range 1-30.
• Mass of the charged Higgs boson m H ± : while the direct search at LEP imposed the constraint m H ± > 80 GeV [30], the measurement of the B → X s γ branching ratio imposes the very stringent bound m H ± > 570 GeV, for t β ∼ 1.5 [31].
• Mass of the CP -odd scalar m A : the authors of Ref. [25] examine the scenarios where either m A or m H is set to a large value about 600-700 GeV while the other one is bounded via theory constraints and experimental data, with the remaining free parameters set to the values mentioned above.We will follow closely this analysis as it is of interest for the present work.We first examine the case of a light CP -odd scalar and a heavy CP -even scalar with mass m H = 600 GeV.In this scenario the searches for the decays A → τ τ , A → γγ, and A → hZ exclude the region m A < 350 GeV, whereas the LHC data on the Higgs boson require m A to be larger than 220 GeV.On the other hand, the search for the channel b b → A → τ τ allows for m A values in the range 350-700 GeV and impose the upper limit t β < 2 for m A ≤ 500 GeV, whereas t β < 15 for m A ≤ 500 GeV.
• Mass of the heavy CP -even scalar m H : we now examine the scenario with a light CP -even scalar and a heavy CP -odd scalar with m A = 700 GeV.In this case the whole constraints require m H > 300 GeV, whereas the b b → H/A → τ τ channel imposes an upper bound on t β as a function of m H .For instance, for m H = 200 (600) GeV t β < 6 (15).On the other hand, the searches for the H decays into τ τ , W W , ZZ, γγ, and hh require tan β > 2.5 for m H < 380 GeV.For a lighter m A = 600 GeV, the search for the A → HZ channel can exclude the region m H < 270 GeV.
We now turn to study the behavior of the A → Zγγ and φ → Zγγ (φ = h, H) branching ratios as functions of the parameters t β , β − α, m H ± , m A and m H .We stick to the still allowed values for these parameters, whereas for the SM parameters we take the values given in Ref. [32].For our analysis we used the LoopTools package [33,34] for the numerical evaluation of the Passarino-Veltman scalar functions appearing in the decay amplitudes.The dominant decay widths of the CP -odd and CP -even scalar bosons were evaluated by our own Mathematica code that implements the formulas of Appendix D, including the QCD corrections for the decays into light quarks.

B. A → Zγγ branching ratio
We work in a region close to the alignment limit and use sin(β − α) = 0.999.In this scenario, the strength of the hZA vertex is negligible and so the contributions to the A → Zγγ decay only arise from box diagrams and reducible diagrams with H exchange, which receive their main contributions from the top and bottom quarks.The contribution of the loops with W gauge bosons turns out to be negligibly small as it is proportional to cos 2 (β − α), whereas the charged scalar boson also gives a very small contribution for m H ± of the order of a few hundred GeVs.
We can distinguish two scenarios of interest: m A < m H + m Z and m A > m H + m Z .Below we examine the behavior of the A → Zγγ branching ratio in such scenarios.We consider the scenario with m H = 600 GeV and analyze the behavior of BR(A → Zγγ) as a function of m A in the range 350 − 650 GeV.For the mixing angle β we consider two values: t β = 2 and t β = 10, which are allowed for m A < 500 GeV and 500 GeV ≤ m A ≤ 700 GeV, respectively.In the upper plots of Fig. 4 we show the behavior of the A → Zγγ branching ratio as function of m A for the two chosen values of t β .We also show the main decay modes of the CP -odd scalar boson: A → b b, A → W − H + , t t, gg, γγ, and Zγ.The decay A → Zh has a negligible branching ratio in the region close to the alignment limit and is not shown in the plots.We note that the main contribution to BR(A → Zγγ) arises from the reducible diagrams with top quarks, whereas the contribution of the loops with charged scalar bosons is negligible.Since in this scenario the intermediary scalar boson H is far from the resonance, the reducible diagram contribution is very small, though it is larger than the box diagram contribution by almost two orders of magnitude.Therefore, the Z → Aγγ branching ratio is thus very small.For instance, for t β = 2, BR(A → Zγγ) is of the order of 10 −11 for m A = 300 GeV with a small increase as m A increases.When t β increases up to 10, BR(A → Zγγ) decreases about one order of magnitude as the top quark contribution is suppressed by a factor of 1/t β .In this region of the parameter space of the type-II THDM, the A → Zγγ branching ratio is considerably smaller than those of the one-loop induced decays A → γγ and A → Zγ.The branching ratios for the main A decay channels are also shown.

Scenario with mA > mH + mZ
We now turn to analyze the scenario where the CP -even scalar is relatively light, with a mass m H = 270 GeV along with a heavier CP -odd scalar with a mass in the range 600-1000 GeV.We use t β = 5 and t β = 10, which are allowed for m A = 600 GeV and m A = 700 GeV, respectively.In this scenario the intermediate scalar boson H is on resonance and the CP -odd scalar can decay as A → ZH with a large branching ratio.The decay A → Zγγ would then proceed in two stages: after the CP -odd scalar boson decays as A → ZH, the on-shell CP -even scalar boson decays into a photon pair H → γγ, namely, A → HZ → Zγγ.The enhancement of BR(A → Zγγ) becomes evident in the lower plots of Fig. 4, where we show its behavior as a function of m A , along with that of the branching ratios of other decay modes of the CP -odd scalar boson.We observe that BR(A → Zγγ) increases up to four orders of magnitude with respect to the result obtained in the scenario with m A < m H + m Z and can reach values of the order of 10 −6 − 10 −5 when m A is in the 600-800 GeV range.In this mass regime, the main decay is A → HZ, which explains why the A → Zγγ decay has such an enhanced branching ratio.For illustrative purpose we also show the branching ratios for the decays A → Z bb and A → Zgg, which arise from the decay A → HZ followed by the decays H → bb and H → gg.We note that the dominant decay channel is A → Z bb.
The above-described behavior of BR(A → Zγγ) is best illustrated in the contour plot on the m A vs m H plane shown in Fig. 5 for two values of t β .We observe that BR(A → Zγγ) can reach its largest values, of the order of 10 −5 , in the region where m A > m A + m Z , whereas it is negligible when m A < m A + m Z .Since the A → Zγγ decay receives its main contribution from the loops with top quark, it decreases as t β increases.

C. H → Zγγ branching ratio
We now analyze the behavior of the H → Zγγ branching ratio as a function of m H in scenarios analogue to those discussed for the CP -odd scalar boson.For sin(β − α) = 0.999, apart from the box diagram contribution, the only contribution from reducible diagrams is that with an intermediary CP -odd scalar boson A, which receives contributions mainly from the top and bottom quarks.The diagram mediated by the Z gauge boson gives a negligible contribution since the HZZ vertex is proportional to cos(β − α).

Scenario with mH < mA + mZ
We consider a heavy CP -odd scalar with a mass m A = 600 GeV and take m H in the range 300-600 GeV.For t β we use the values 3 and 10.In the upper plots of Fig. 6 we show the branching ratios for the main decay channels of the H scalar boson.We note that the H → Zγγ decay has a very suppressed branching ratio up to five orders of magnitude smaller than the branching ratios of the one-loop induced decays H → γγ and H → Zγ.It increases for smaller t β but it seems still beyond the reach of detection.

Scenario with mH > mA + mZ
In this scenario we consider m A = 350 GeV and take m H in the range 600-1000 GeV.We also use t β = 2 and t β = 10.For the mass of the charged scalar boson we use m H ± = 575 GeV as we do not need to assume that m H ± > m H since the H → W − H + decay channel has a negligible branching ratio proportional to cos(β − α) 2 .In the lower plots of Fig. 6 we show the H → Zγγ branching ratio along with those of the main decay channels.We note that there is a considerable enhancement of BR(H → Zγγ), up to 5 orders of magnitude, now that the H → ZA decay is allowed, thus BR(H → Zγγ) can be as large as 10 −3 for t β = 2. Again we include the decays H → ZA → Z bb and H → ZA → Zgg, with the dominant decay being the H → ZA → Z bb.
In Fig. 7, we also show the contour plot of BR(H → Zγγ) in the m H vs m A plane for two values of t β .Again it is evident that BR(H → Zγγ) can reach its largest values when m H > m A + m Z and it is negligible when Finally, we would like to comment shortly on the potential detection of the A → Zγγ and φ → Zγγ (φ = h, H) decays in the scenario we are considering in type-II THDM.Although there is a considerable enhancement of the corresponding branching ratios, it still seems not enough to put these decays at the reach of experimental detection at the LHC in the near future.In Fig. 8 we show the leading order production cross section for the CP -even and CP -odd scalar bosons via gluon fusion at the LHC at √ s = 14 TeV as a function of the scalar boson mass.It turns out that with an integrated luminosity of 300 fb −1 , to be achieved in LHC run 3, we would have about 1.64 × 10 5 (3.2 × 10 5 ) CP -even (CP -odd) scalar bosons with a mass m φ = 500 GeV produced per year, but these numbers drop by one order of magnitude when m φ = 700 GeV.For BR(H → Zγγ) O(10 −3 ), we would only have about 164 H → Zγγ events prior to imposing the kinematic cuts, which would render this decay hard to detect.The situation might be more promising at a future high-luminosity 100 TeV pp collider, where we could have thousands of H → Zγγ events prior to imposing the kinematic cuts.As discussed below, this event number would increase in type-I THDM by one order of magnitude as the respective branching ratios would have such an enhancement in that model.

D. h → Zγγ branching ratio
We now briefly discuss the lightest CP -even scalar boson decay h → Zγγ.Since h must mimic the properties of the SM Higgs boson, it is expected that the h → Zγγ branching ratio does not deviate considerably from its SM value.For sin(β − α) 1, the only contributions arise from box diagrams and the reducible diagram mediated by the Z gauge boson.As mentioned before, the hZA vertex is considerably suppressed, whereas the hHZ one is forbidden due to CP invariance.Even more, there is no enhancement due to the Z-mediated reducible diagram since m h < 2m Z .For m h = 125 GeV and t β = 10 GeV we obtain BR(h → Zγγ) 10 −9 , which does not deviate significantly from the SM value [16].Therefore the new physics effects provided by the THDM does not give a significant enhancement to this decay and seem very far from the reach of detection.

E. Kinematic distributions
In the scenario in which the intermediary scalar boson is off-shell, the analysis of the behavior of some kinematic distributions could be helpful to disentangle the decay signal from the potential background.The Z gauge boson energy distribution dΓ(φ i → Zγγ)/dE Z and the photon invariant mass distribution dΓ(φ i → Zγγ)/dm m γγ could be useful for this task.To obtain the former, one can plug the relation dx 1 = (2/m φi )dE Z into Eq.( 21) to obtain where the Z gauge boson energy is defined in the interval (m Z , (m 2 φi + m 2 Z )/(2m φi )).On the other hand, the expression for the photon invariant mass distribution dΓ(φ i → Zγγ)/dm m γγ is obtained using the relation where m γγ is defined in the interval (0, m φi − m Z ).
For illustrative purpose we show in Fig. 9 the energy distribution dΓ(A → Zγγ)/dE Z and the photon invariant mass distribution dΓ(A → Zγγ)/dm γγ in the alignment limit for m H = 700 GeV, m H ± = 750 GeV, t β = 15, and a few values of m A .We observe that in the rest frame of the CP -odd scalar, the Z gauge boson energy is peaked at about one half of m A .A similar situation is observed for the invariant mass m γγ .We now briefly analyze these decays in the framework of type-I THDM, where the charged scalar boson mass has a lower bound.Since the main contribution to the decays A → γγ and H → γγ arises from the top quark, the effect of a charged Higgs scalar boson with a mass less than 570 GeV would not have a considerably effect on the decays we are interested in.However, in type-I THDM the couplings of the CP -odd scalar boson to both quark types are now proportional to cot β, and the same is true for the couplings of CP -even scalar bosons in the cos(α − β) → 0 limit [37].Therefore, the decay widths of the scalar bosons into the bb pair would be suppressed for large t β , which can do have an effect on our decays indeed.Consider for instance the decay A → Zγγ in the scenario where m A > m H + m Z .In type-II THDM the main decay channel is A → ZH → Z bb, but for large t β this decay would get suppressed in type-I THDM as the H → bb decay gets suppressed.This can translate into an enhancement of the A → ZH → Zγγ decay width.To analyze this scenario we have performed the explicit calculation of the A → Zγγ and H → Zγγ decay widths in type-I THDM in the same scenarios considered in type-II THDM.The results for the respective branching ratios and those of the main decay channels are shown in Fig. 10, where we observe that the A → Zγγ and H → Zγγ decays can have an enhancement of about one order of magnitude with respect to the values obtained on type-II THDM.

V. CONCLUSIONS
In this work we have calculated the one-loop contributions to the decays of the CP -odd and CP -even scalar bosons A → Zγγ and φ → Zγγ (φ = h, H) in the framework of THDMs.We have presented analytical expressions for both box and reducible diagrams in terms of Passarino-Veltman scalar functions, though the main contributions arise from reducible diagrams.We first discuss the A → Zγγ decay, which has not been discussed previously in the literature to our knowledge.For the numerical analysis we worked within the type-II THDM and considered a region of the parameter space still consistent with experimental data, with sin(β − α) 1, where the lightest CP -even scalar boson h is identified with the SM Higgs boson, the hZA vertex has a negligibly small strenght, and the heavy CP -even scalar does not couple to the weak gauge bosons.It was found that the A → Zγγ branching ratio is only relevant in the scenario where m A > m H + m Z , when the intermediary H boson is on-shell.For m A > 600 GeV and t β close to 1, BR(A → Zγγ) can reach values of the order of 10 −5 − 10 −4 , but it decreases by about one order of magnitude as t β increases up to 10, which stems from the fact that the dominant contribution arises from the loops with the top quark, which couples to the scalar boson with a strength proportional to 1/t β .On the other hand, when m A < m H + m Z , BR(A → Zγγ) is negligibly small, of the order of 10 −10 .As far as the H → Zγγ decay is concerned, it exhibits a similar behavior and its branching ratio is non-negligible only in the scenario where m H > m A + m Z , when the CP -odd scalar is now on-shell.In this region of the parameter space, BR(H → Zγγ) can reach the level 10 -6

Zhh
FIG. 10: Branching ratios for the A → Zγγ and H → Zγγ decays in type-I THDM as a function of the scalar boson masses for m H ± = 570 GeV, sin(β − α) = 0.999, and two values of t β allowed by theory and experimental constraints.The branching ratios for the main decay channels are also shown.
of 10 −4 − 10 −3 for m H > 600 GeV and t β 1, but it decreases for larger t β .We also discussed the h → Zγγ decay, which receives contribution from box diagrams and a reducible diagram mediated by the Z gauge boson.Since the properties of the h scalar boson are nearly identical to the SM Higgs boson, it is found that the h → Zγγ branching ratio does not deviates significantly from the SM prediction and it is of the order of 10 −9 .Our calculation is in agreement with previous evaluations.The new physics effects of THDMs are thus not relevant for this decay.Finally we also estimated these rare decays in the framework of type-I THDM, where we find that the respective branching ratios can be enhanced buy about one order of magnitude with respect to those of type-II THDM.
12: Feynman rules for the couplings of the scalar bosons to fermions in THDMs.The corresponding coupling constants for type-II THDM are shown in Table II.
TABLE II: Constants for the couplings of the scalar bosons to fermions and gauge bosons in type-II THDM as described in Figs. 12 and 13.We have used the short-hand notation sa = sin a and ca = cos a.The g φZZ couplings obey g φZZ = 1 c 2 W g φW W [4].
CP -even scalar bosons to a pair of charged scalar boson φH − H + , which emerge from the Higgs potential (3) once it is diagonalized.We now present the form factors of Eqs. ( 17) and (19) in terms of Passarino-Veltman scalar functions.
(c)  17) and where the kinematical invariant variables s 1 , s 2 and s were defined in Eqs. ( 13)- (15).In addition, we use the following auxiliary variables: for i = 1, 2 and j = A, Z.As for the three-and four-point Passarino-Veltman scalar functions C i and D i , they are defined as (B6) As we can see for Eq. ( B2)-(B4) the box diagrams amplitudes are free of ultraviolet divergences since are free of two-point Passarino-Veltman scalar functions.

b. Reducible diagram contribution
The reducible diagram of Fig. 2 only contribute to the form factor F 1 of Eq. ( 17).The Passarino-Veltman technique allowed us to obtain the following results for the fermion and W gauge boson contributions where the three-point scalar function C(s, m 2 χ ) can be written in terms of elementary functions as follows where f (x) is given in Eq. (30).The box diagram contributions to the form factors of Eq. ( 19) are given as follows and with the two-point Passarino-Veltman scalar functions defined as ∆B(r 2 1 , r It is also evident that ultraviolet divergence cancel out.

b. Reducible diagram contribution
The reducible diagrams related to the processes φ → Zχ * → Zγγ, with χ = A, Z, yield the following contribution to the form factor of Eq. ( 20) Appendix C: Squared average amplitudes From the general form of the invariant amplitudes for the A → Zγγ and φ → Zγγ (φ = h, H) decays presented in Eqs. ( 17) and ( 19), respectively, we can readily obtain the square amplitudes averaged over photon and Z polarizations, which are required for the calculation of the decay width (21).The results can be written as follows.

1.
A → Zγγ decay and From Eq. ( 19) we obtain with Gi (s, s 1 , s 2 ) = G i (s, s 2 , s 1 ) and Appendix D: Decay widths of CP -even and CP -odd scalar bosons For completeness, we present the expressions for the most relevant A → X and φ → X (φ = h, H) decays, with X a final multiparticle state.These formulas have been summarized for instance in [3,35,36].We use the notation introduced in the Feynman rules shown in Figs. 12 and 13.

CP -even scalar boson decays
The tree-level two-body decay width into fermion pairs is with f φ f f = gm f g φ f f /(2m W ), where the g φ f f constants are shown in Table II for type-II THDM.Also, we use the definition τ a = 4m 2 a /m 2 φ and N f c stands for the fermion color number.The widths of the decays into a pair of on-shell gauge bosons V = W Z, when kinematically allowed, are given by with n V = 1 (2) for V = W (Z).Here f φW W = gm W g φW W and f φZZ = gm W g φW W /c 2 W , where again the g φV V constants are shown in Table II for type-II THDM.
For the present work another relevant decay is φ → ZA, whose decay width was already presented in Eq. (32), which can also be useful to compute the φ → W ∓ H ± decay when kinematically allowed.On the other hand, we will assume that other tree-level decays such as φ → AA and φ → H − H + are not kinematically allowed and we refrain from presenting the respective decay widths here.
One-loop decays can also be important for Higgs boson phenomenology: while the decay φ → γγ has a clean signature, the decay φ → gg is important for the cross section of Higgs boson production via gluon fusion.As for the φ → γγ decay width, it is given in Eqs. ( 27)- (29), which can also used for the two-gluon decay width by taking the quark contribution only and making the replacements α 2 → 2α 2 S and N f c Q 2 f → 1.The φ → Zγ decay has also been largely studied in the literature.The decay width can be written as with F φZγ = F φZγ f (τ f , ξ f ) + F φZγ W (τ f , ξ W ) + F φZγ H ± (τ H ± , ξ H ± ).The contributions of charged fermions, the W gauge boson, and the charged scalar are given by where we introduced the definition ξ i = 4m 2 i /m 2 Z .

CP -odd scalar boson decays
The decay of a CP -odd scalar boson A into a pair of fermions of distinct flavor is given by where now we use the definition τ a = 4m 2 a /m 2 A .There are no decays into pairs of electroweak gauge bosons at the tree-level, but the A → φZ (φ = h, H) decay can be kinematically allowed.Its decay width is given in Eq. ( 32) and a similar expression with the corresponding replacements is obeyed by the A → W ± H ∓ decay if kinematically allowed.
As far as one-loop decays are concerned, the two-photon decay proceeds via charged fermion loops and its decay width is presented in Eqs.(27) and (31), whereas the two-gluon decay width can be obtained from these equations by summing over quarks only and making the additional replacements α 2 → 2α 2 S and N f c Q 2 f → 1.The A → Zγ decay also receives contribution from charged fermions only and its decay width is given by Eq. (D3), with φ → A and 3. QCD radiative corrrections for the decays φ → qq For light quarks, the running mass mq at the scale m φ must be used in Eqs (D1) and (D5) to take into account the next-to-leading order QCD corrections.As for higher order QCD corrections, they are important and must be also included.They are summarized in [36] and we include them here for completeness.For light quarks we have where p = 1 (3) for CP -even (CP -odd) scalar boson and the running quark mass mq is defined at the scale m φ .As for ∆ qq , it is the same for both CP -even and CP -odd scalar bosons for m φ m q .In the MS renormalization scheme it is given by ∆ qq = 5.For the top quark, the leading order QCD corrections are give by [36] Γ(φ → tt) = 3g 2 g 2 φ tt m t m φ 32πm 2 W (1 − τ t ) p/2 1 + 4 3 with β t = 1 − τ t , whereas ∆ t φ (β) is given in Ref. [36].However, these corrections are small compared to the case of the b and c quarks.

FIG. 3 :
FIG.3: Feynman diagram that also contributes to the φ → Zγγ (φ = h, H) decay in the THDM, in addition to Feynman diagrams analogue to those of Figs.1 and 2. The diagram obtained by exchanging the photons is not shown.

FIG. 4 :
FIG.4: Branching ratio for the A → Zγγ decay in type-II THDM as a function of mA for m H ± = 570 GeV, sin(β − α) = 0.999, and two values of t β allowed by theory and experimental constraints.In the upper (lower) plots we use mH = 600 (270) GeV.The branching ratios for the main A decay channels are also shown.

FIG. 6 :- 6 7. 4 × 10 - 6 FIG. 7 :
FIG.6: Branching ratio for the H → Zγγ decay in type-II THDM as a function of mH for m H ± = 570 GeV, sin(β − α) = 0.999, and two values of t β allowed by theory and experimental constraints.In the upper (lower) plots we use mA = 600 (350) GeV.The branching ratios for the main H decay channels are also shown.

2 FIG. 8 :
FIG.8: CP -even and CP -odd scalar boson production cross section via gluon fusion as a function of the scalar boson mass at the LHC at √ s = 14 TeV in type-II THDM.We use sin(α − β) 0.999 and t β = 2.The right axis shows the annual event number achieved with an integrated luminosity of 300 fb −1 .

FIG. 13 :
FIG. 13: Feynman rules necessary for our calculation in THDMs.All the four-momenta are incoming.Here φ = h, H in diagram (b) but φ = h, H, A in diagram (f).Also, in diagram (e) g AH + H − = sW and g ZH + H − = c2W /2cW .The remaining coupling constants for type-II THDM are presented in TableII.

2 .
φ → Zγγ (φ = h, H) decay a. Box diagram contribution 67 ᾱsπ + (35.94 − 1.36N f ) ᾱ2 s π 2 + . . .,(D8)where N f is the number of flavors of light quarks and ᾱs is the strong coupling constant defined at m φ scale.As for ∆ φ , it differs for CP -even or CP -odd scalar bosons and it is given at order ᾱ2

TABLE I :
Couplings of quarks and leptons to the Higgs doublets Φi in THDMs with natural flavor conservation.The superscript i stands for the generation index.