Decay $\phi\to Z \gamma\gamma$ ($\phi=h, H,A$) in the minimal supersymmetric standard model

The decays of the CP-even ($h,\,H$) and the CP-odd ($A$) Higgs bosons $h,\,H,\,A\to Z\gamma\gamma$ are calculated in the context of the minimal supersymmetric standard model (MSSM), where they are induced at the one-loop level via box and reducible Feynman diagrams. For the numerical evaluation of the branching ratios we employ the decoupling limit and consider values for the MSSM parameters still consistent with the current experimental data. We consider a particular scenario for the radiative corrections to the Higgs boson masses, dubbed hMSSM, which leads to only two free parameters, namely $m_A$ and $\tan\beta$. We found that the branching ratio of the $CP$-odd Higgs boson decay $A\to Z\gamma\gamma$ can reach values up to $10^{-8}$ for $m_A=750$ GeV, which is three order of magnitude smaller than that of the $A\to\gamma\gamma$ channel. On the other hand, the branching ratio of the $H\to Z\gamma\gamma$ decay is of the order of $10^{-7}$ for large values of $m_H$ and small values of $\tan\beta$, which is comparable with the branching ratio of the $H\to Z\gamma$ channel. As far as the branching ratio of the SM-like Higgs boson decay $h\to Z\gamma\gamma$ is concerned, it is considerably suppressed, of the order of $10^{-11}$.


I. INTRODUCTION
The discovery of a Higgs boson with a 125 GeV mass in the Run 1 of the Large Hadron Collider (LHC) by the ATLAS and CMS Collaborations [1,2] is a milestone in particle physics and the standard model (SM), which only predicts one of such particles. Nonetheless, the search for additional Higgs bosons predicted by theories with extended scalar sectors is still underway. A plethora of beyond the SM theories (BSMTs) of this class have been proposed to improve the SM but also to explore effects absent in this theory. High energy experiments have looked for signals of new physics (NP), either through the direct production of new particles or the search for evidences of anomalous contributions to the couplings of the SM particles. However, up to now no evidences along these lines have been found, thereby shifting the bounds on the associated scales into the TeV scale. In spite of this scenario, the search for new Higgs bosons may provide signals of NP and give a hint to the path to extend the SM. Hopefully any additional scalar boson would show up at the LHC in a near future.
Among the most popular BSMTs, the minimal supersymmetric standard model (MSSM) [3,4] stands out and has been the focus of considerable attention in the literature during the last decades. Supersymmetry (SUSY) associate bosons with fermions, predicting the existence of fermion partners for all SM bosons as well as scalar partners for all SM fermions. Unlike the SM scalar sector, which requires only one Higgs doublet to trigger the electroweak symmetry breaking, two Higgs doublets are necessary in the framework of the MSSM. After the electroweak symmetry is broken, five physical Higgs bosons emerge: two CP -even Higgs bosons H and h, with m H > m h by convention, a CP -odd Higgs boson A, often called pseudoscalar, and a pair of charged Higgs bosons H ± . The properties of one of the neutral CP -even scalar bosons must match those of the 125 GeV scalar boson detected at the LHC. Of particular importance for our work are the superpartners of the Higgs bosons and gauge bosons (Higgsinos and gauginos), which mix to give rise to neutralino and chargino states χ ± . Extensive studies have been carried out to constraint the allowed region of the parameter space of the MSSM via low energy observables along with the data from direct searches at the LHC [5,6]. Furthermore, the LHC data on the 125 GeV Higgs boson obtained at √ s = 7 and 8 GeV are in agreement with the SM predictions for the Higgs boson couplings to the heavier SM particles, thereby placing tight constraints on the parameter space of BSMTs. Along this line, it has been found that LHC data gives rise to a strong tension with the naturalness within the MSSM since the mass of the stop sector is required to be rather heavy, typically above 1 TeV. One-loop induced Higgs boson decays into vector bosons such as φ → γγ and φ → Zγ (φ = h, H, A) have been long studied in the literature within the context of several BSMTs. These decays are interesting since they are induced by loops carrying heavy charged particles and are thus sensitive to the charge and spin of new charged particles. Within the MSSM, a potential enhancement of the φ → γγ decay width may arise from loops carrying charged Higgs bosons, sfermions, and charginos. This possibility was considered in [7] to explain the apparent excess of the diphoton decay rate of the 125 GeV Higgs boson reported by the ATLAS and CMS Collaborations [1,2]. In the scenario with heavy sfermions, the chargino contribution is the only viable option that can induce a significant enhancement of the h → γγ decay rate, which can only reach the 10% level for typical values of the SUSY parameters, but it can be considerably larger for light charginos and low tan β [7]. In particular, there is a narrow area of the MSSM parameter space where the h → γγ rate can be enhanced up to 25%, whereas a less substantial increase of around 10 − 20% can be reached in a much wider area [8]. Similar studies have been carried out for the φZγ (φ = h, H, A) coupling [9]. Although no excess in the h → γγ decay rate was observed in the data collected at the LHC Run 2, it is possible that charginos can give an enhancement to the diphoton decay rate of heavy scalar bosons over the contributions of non-SUSY particles.
If any additional Higgs scalar boson is discovered, it would be interesting to test its exotic decay modes as they may shed light on the underlying theory. In this work we study the one-loop induced three-body decays of the CP -odd and CP -even Higgs bosons φ → Zγγ (φ = h, H, A) in the context of the MSSM, where chargino loops can enhance significantly the decay rate as compared to the fermion and squark contributions. These decays receive contributions from both box and reducible Feynman diagrams. Experimentally this decay channel has the advantage of a relatively low background since the final state contains two photons plus a pair of energetic back-to-back leptons. In addition, as occurs with the diphoton coupling, the AZγγ and HZγγ ones are also sensitive to the charge and spin of the particles running into the loops, such as the charginos of the MSSM. Within the SM framework the decays h → Zγγ [10] and the analogue one h → Zgg [11] have been studied. More recently, the A → Zγγ and H → Zγγ decays were studied in the context of THDMs [12]. Also, the φZγγ (φ = h, H, A) coupling was analyzed via photon fusion γγ → Zφ in the context of the SM and the MSSM [13]. In our study, we consider the existing limits on the parameter space of the MSSM. We work in a scenario for the radiative corrections to the Higgs boson masses, dubbed hMSSM (habemus MSSM? [14]), which leads to only two free parameters, namely m A and tan β. Since the current limit on chargino masses is mχ± > 103.5 GeV [15], it is worth assessing if there is a significant enhancement to the φ → Zγγ (φ = H, A) decay rate from chargino loops.
The organization of this paper is as follows. In section II a brief outline of the MSSM is presented, with focus on the Higgs sector as well as the chargino and squark sectors. Section III is devoted to the calculation of the decay φ → Zγγ (φ = H, A) through the Passarino-Veltman reduction scheme. In Section IV a brief discussion on the allowed parameter space of the MSSM is presented along with the numerical analysis of the φ → Zγγ (φ = H, A) decay width. Finally in Section V we present the conclusions and outlook. The Feynman rules involved in our calculations as well as some lengthy formulas are presented in Appendices A and B.

II. OVERVIEW OF THE MINIMAL SUPERSYMMETRIC STANDARD MODEL
SUSY predicts new scalar particles called squarks and sleptons, which are superpartners of the quarks and leptons. Also, gauge bosons have fermionic partners, called gauginos. In the MSSM [16], which is the minimal realization of SUSY, two complex Higgs doublets with opposite hypercharges are required to give masses to the up-and downtype fermions and ensuring anomaly cancellation [17]. The superpartners of the physical Higgs bosons are known as Higgsinos. Below we present an outline of the MSSM Higgs sector, focusing only on those aspects relevant for our calculation.

A. Higgs sector
The Higgs sector of the MSSM is essentially the same as that of the Type II-THDM [18], where two complex isodoublets are introduced The H 1 doublet couples to the down-type quarks and the charged leptons while the H 2 Higgs doublet only couples to the up-type quarks. After the electroweak symmetry breaking, three of the original eight degrees of freedom of H 1 and H 2 become the longitudinal polarizations of the W ± and Z gauge bosons, which thus acquire masses, whereas the remaining degrees of freedom correspond to the physical Higgs spectrum, which is determined as follows. The neutral components of the two Higgs fields develop vacuum expectations values (VEVs) where ∆M 2 ij stands for the radiative corrections that depend on the SUSY scale, the stop/bottom trilinear couplings A t/b , and the Higgsino mass µ. After the diagonalization of M 2 S and the rotation (4), one obtains the masses of the physical neutral Higgs bosons up to radiative corrections where Note that m h ≤ m Z if radiative corrections are not taken into account. In a simplified version of the model, the so called hMSMS [14], it is assumed that ∆M 22 , which includes the dominant stop corrections, is the only relevant correction term, i.e. ∆M 22 ∆M 11 , ∆M 12 . Thus if the lightest CP -even scalar boson h is assumed to be the one found at the LHC, ∆M 22 can be traded for the already known mass m h as follows [14] Therefore m H can be written as and the mixing angle α obeys There are thus two free parameters that are usually chosen as t β ≡ tan β and m A . Our analysis below will be focused on the hMSSM scenario.

B. Chargino Sector
The superpartners of the charged Higgs boson and the W gauge boson mix to give rise to the chargino mass eigenstatesχ ± 1,2 . The mass eigenvalues mχ± 1,2 , the mixing angles, and the phases are determined by the elements of the chargino mass matrix, which in the (W − ,H − ) basis is given by thereby being determined by the fundamental SUSY parameters: the gaugino mass M 2 , the Higgs mass µ and the mixing angle t β . In several models, the lightest chargino stateχ + 1 is expected to be one of the lightest SUSY particle and might play an important role in the experimental detection of SUSY.
The Lagrangian for the mass terms of the charginos can be written in the gaugino-Higgsino eigenbasis as where Ψ L and Ψ R are two-components Weyl spinors The chargino mass matrix X can be diagonalized by rotating the wino and Higgsino fields via 2×2 unitary matrices, namely, χ L = V Ψ L and χ R = U Ψ R , with χ L,R being the chargino mass eigenstates. This leads to where If CP invariance is assumed, U and V are orthogonal matrices, whereas in CP -violating theories, the parameters M 2 and µ can be complex. By a reparametrization of the fields, M 2 can be taken as real and positive [19], so that the only non-trivial invariant phase is that of µ, namely, µ = |µ|e iθ The unitary matrices U y V must be chosen so that the elements of the diagonal matrix Mχ + are real and nonnegative. The corresponding eigenvalues are given by where We can observe from Eqs. (16) and (17) that the chargino masses mχ+ i depend mainly on µ, M 2 and t β , which are constrained by the LEP bound on the mass of the lightest chargino: mχ± > 103.5 GeV.

C. Sfermion sector
In addition to t β and µ, the sfermion sector can be described by three additional parameters for each sfermion kind: the left-and right-handed soft SUSY-breaking scalar masses mf L and mf R along with the trilinear coupling A f . The sfermion mass matrix can be written as with In order to obtain the mass eigenstatesf 1 andf 2 , the sfermion matrices are diagonalized by a rotation by an angle θ f via a 2 × 2 unitary matrix The mixing angle and sfermion masses are then given by It is usually assumed that the masses of the first and second generation sfermions are zero [20]. For the calculation of the φ → Zγγ (φ = H, A) decay we find it convenient to work in the unitary gauge. The necessary Feynman rules are shown in Appendix A. Below we present the analytical expression for the invariant amplitude and the corresponding decay width.
We turn to the most relevant aspects of our calculation. In the unitary gauge, the decay φ → Zγγ (φ = h, H, A) is induced at the one-loop level by charged particles through the box and reducible Feynman diagrams shown in Figs. 1 and 2. Apart from the loops with SM charged fermions and the W gauge boson, in the MSSM there are contributions from the charged scalar boson, charginos, and sfermions. Due to CP invariance, the W gauge boson, the charged scalar boson, and the sfermions only contribute through reducible diagrams.
We first define the kinematics conditions and introduce a set of invariant variables meant to simplify our analytical results. The nomenclature for the four-momenta of the external particles is as follows with p 2 = m 2 φ , k 2 1 = k 2 2 = 0 and k 2 3 = m 2 Z due to the mass-shell condition. We now introduce the following Lorentz invariant kinematic variables which obey s + s 1 + s 2 = m 2 φ + m 2 Z by four-momentum conservation. All the scalar products between the external four-momenta can be expressed in terms of s, s 1 and s 2 , together with the scaled variable µ Z = m 2 Z /m 2 φ , so that the invariant amplitudes of the φ → Zγγ (φ = h, H, A) decay can be written in terms of these variables. Furthermore, using the transversality condition for the external gauge bosons, i.e., k 1 · µ (k 1 ) = k 2 · ν (k 2 ) = k 3 · α (k 3 ) = 0, all the terms proportional to k µ 1 , k ν 2 and k α 3 can be dropped from the invariant amplitude, which simplifies considerably the calculation.
After writing out the corresponding amplitude for each Feynman diagram, the one-loop integrals were reduced to scalar integrals via the Passarino-Veltman reduction scheme [21], which was carried out with the aid of the Mathematica package FeynCalc [22]. We first present the results for the A → Zγγ decay.
A. A → Zγγ decay

Box diagram contribution
This decay receive contributions of box diagrams and reducible diagrams. The box diagrams are shown in Fig.  1: the particles running into the loops are SM charged fermions and charginos. By CP invariance, the pseudoscalar boson cannot couple to a pair of W gauge bosons or charged scalar bosons. Also, it cannot couple to a pair of identical sfermions by Bose symmetry. Therefore, the only new contribution from the MSSM to box diagrams arises from charginos. Although there can be two charginos of distinct flavor into the box diagrams (there are non-diagonal couplings Zχ + 1 χ + 2 and Aχ + 1 χ + 2 in the MSSM) we only consider the contribution of a single chargino as non-diagonal chargino couplings are much more suppressed than the diagonal ones. In the case of the SM fermions, the dominant contributions arise from the heaviest quarks: for small t β , the top quark box dominates, whereas for large t β the bottom quark box becomes relevant. This is due to the presence of the factors 1/t β and t β in the corresponding Yukawa couplings. Once the Passarino-Veltman reduction was applied to reduce the amplitude of each box diagram into a combination of Passarino-Veltman scalar functions, we verified that the total amplitude is ultraviolet finite and obeys U (1) em gauge invariance as well as Bose symmetry. Although each box diagram has ultraviolet divergences, they are cancelled out when summing over all diagrams. The full amplitude can be arranged in the following manifestly gauge-invariant form with the Lorentz structures given as follows The form factors F Box i depend on the variables s, s 1 , s 2 , and µ Z . The explicit expressions are somewhat cumbersome and are presented in Appendix B.

Reducible diagram contribution
The decay A → Zγγ is also induced by the reducible Feynman diagrams shown in Fig. 2, which are mediated by a CP -even Higgs boson φ = h, H decaying into a photon pair via triangle and bubble diagrams. The decay thus proceeds as follows A → Zφ * → Zγγ. The AZZ coupling is absent at tree level and so are reducible diagrams mediated by a Z gauge boson. Reducible diagrams also yield an ultraviolet-finite and gauge invariant amplitude by themselves. They contribute to the general amplitude of Eq. (27) through the form factor F 1 via the loops with SM charged fermions, the W gauge boson, charginos, sfermions, and the charged scalar boson H ± . Note that two identical sfermions do couple to a CP -even Higgs boson and so they can contribute to the triangle diagrams. The triangle diagram contribution to the F 1 form factor can be written as where φ = h, H and the F φ−X

1
(X = f, W, H ± ,χ + i ,q i ) functions are presented in Appendix B in terms of Passarino-Veltman scalar functions.
The total contribution to the F i form factors is thus given by These diagrams are similar to those shown in Fig. 1 and receive contributions from SM charged fermions and charginos. Although there can be box diagrams with the W gauge bosons, the charged scalar boson, and sfermions, these contributions exactly cancel out: it turns out that by CP invariance the φZγγ vertex function is composed by Levi-Civita tensors, which can only arise from the amplitudes of fermions loops involving a trace of a Dirac chain with a γ 5 matrix.
The most general Lorentz structure for the φ → Zγγ (φ = h, H) decay can be written in the following U (1) em gauge-invariant manifest form where we use the shorthand notation αkpq = αβλρ k β p λ q ρ , etc. To arrive to this expression, we used Schouten's identity. The form factors G i , which depend on the kinematic variables s, s 1 , s 2 , and µ Z , are presented in Appendix B in terms of Passarino-Veltman scalar functions.

Reducible diagram contribution
There are two kinds of reducible Feynman diagrams. Those of the first kind are similar to the ones shown in Fig.  2, but with the virtual CP -even Higgs boson replaced by the CP -odd Higgs boson A, which means that there are only loops with SM charged fermions and charginos. For the reasons explained above, there are no triangle loops with virtual W gauge bosons, charged scalar bosons, or sfermions. There is also another kind of reducible Feynman diagrams mediated by the Z gauge boson as shown in Fig. 3 (the well known triangle anomaly Z * γγ), which receive contributions from charged fermions only: only fermion loops can give rise to the Levi-Civita tensor appearing in the corresponding vertex function via the trace of a Dirac chain with a γ 5 matrix. Z-mediated reducible diagrams can also get chargino contributions, but they are proportional to gχ i A , which exactly cancels in the hMSMM. Reducible diagrams only contribute to the form factor G 3 , which reads where stand for the contributions of the reducible Feynman diagrams mediated by the CP -odd Higgs boson and the Z gauge boson, repectively. The corresponding expressions are presented in Appendix B. The square of the invariant amplitudes of Eqs. (27) and (29), averaged over polarizations of the photons and the Z gauge boson, can be obtained after some lengthy algebra. The results are shown in Appendix C and can be plugged into the following formula to obtain the decay width for φ = h, H, A. In the rest frame of the decaying Higgs boson, the scaled variables s, s 1 and s 2 , defined in Eq. (25), are related to the energy of the Z boson (one of the photons) E Z (E γ ) according to where Below we analyze the behavior of the φ → Zγγ (φ = h, H, A) decay width in the parameter space of the MSSM still allowed by current constraints.

IV. NUMERICAL ANALYSIS AND RESULTS
The MSSM parameters that are relevant to our study are mostly associated to the electroweakino and Higgs sectors, namely, M 2 , µ, m A , t β , as well as the sfermion and chargino masses. We first present an overview of the allowed values of these parameters, according to experimental and theoretical constraints, focusing on the so called hMSSM scenario. Afterwards we evaluate the φ → Zγγ (φ = h, H, A) decay in the allowed region of the parameter space.
A. Allowed parameter space of the hMSSM Until now, most of the observed effects at high energy experiments have been found compatible with the SM: the experimental data agree with the theoretical predictions to a high accuracy. The status of the SM was cemented by the Higgs boson discovery at the LHC [23]. Subsequent measurements [24] leave small room for NP effects and so can have strong implications on the parameter space of BSMTs such as the MSSM. In particular, the parameters associated with the scalar sector of the MSSM have become tightly constrained by the measurement of the couplings of the 125 GeV Higgs boson, whose properties must be reproduced by the lightest CP -even Higgs boson h of the MSSM. This favors regions very close to the alignment limit, sin(β − α) = 1, in which the h couplings to SM particles are SM-like. In such a limit, the mass eigenbasis of the CP -even sector coincides with the one in which the gauge bosons receive their masses from only one of the two Higgs doublets. Also, the parameters α and β are related by tan α = −1/ tan β, and the couplings of the heavy CP -even Higgs boson to the weak gauge bosons vanish, although it can still couple to the fermions. Furthermore, the non-observation of any additional Higgs bosons at the LHC implies that their masses are likely to lie above the elecroweak scale i.e. mφ m Z . For our calculation we consider a region close to the decoupling limit of the MSSM [20] in which sin(β − α) = 1 and consider that the new Higgs bosons are heavy.

SUSY particle masses
The discovery of the Higgs boson at the LHC also has had severe implications on the sfermion sector as the stop mass is required to be of the order of 1-10 TeV to accommodate a 125 GeV scalar boson. However, the tension between the experimental and theoretical values of the muon anomalous magnetic moment can be alleviated by the presence of SUSY particles with masses at the GeV scale, which are also required by the little hierarchy problem. This suggests that SUSY particles such as electroweakinos may have masses at the GeV scale. The search for such particles is one of the main goals of the LHC after the Higgs discovery. Along this line, the non-observation of squarks and gluinos at the early stage of the LHC, hints that the masses of such particles lie in the TeV scale. Also, direct electroweakino searches have imposed strong bounds on their masses [25][26][27] via the chargino-neutralino pair production channel. The most up-to-date analysis constrains the lightest neutralino mass to the 100 − 150 GeV interval, whereas the lower bound on the mass of a degenerate wino-like χ 0 2 and χ ± 1 is about 300 GeV. In view of these constraints, a heavy Higgs boson with a mass below 500 − 600 GeV is not kinematically allowed to decay into an electroweakino final state.
2. CP -odd Higgs scalar boson mass mA and ratio of VEVs t β Direct searches of heavy scalar resonances at the LHC are useful to contrain m A and t β . Recently, the combined LHC data at √ s = 7 and 8 TeV allowed the ATLAS and CMS Collaborations searching for neutral scalar bosons and constraining the parameter space of the MSSM via the following decay channels H → γγ [28], H → W W [29], H → ZZ [30] H → hh → bbbb, H → hh → bbγγ [31], A → Zh [32] and H/A → τ + τ − [33], with the latter channel imposing the most stringent bounds. The most up-to-date constraints arise from the analysis of the LHC data for Higgs boson production in association with bottom quarksbbH/A and via gluon fusion at √ s = 13 TeV [34,35].
Such data were used by the ATLAS Collaboration to perform a search for neutral Higgs bosons in the 200-1200 GeV mass range decaying into the τ + τ − channel [34]. It was found that for m A = 200 GeV (1500 GeV), the region with t β > 1 (t β > 42) is excluded at 95% C.L. A similar study by the CMS Collaboration [35], yields slightly less stringent constraints. Other decay channels such as H → ZZ and A → hZ, are useful to exclude the region with t β 4 and m A 350 GeV [36], but the region with m A 400 GeV and t β 10 is still allowed.
As for indirect constraints on the MSSM parameters, the most stringent bounds can be obtained from the study of the rare B-meson decays B s → X s γ and B s → µ + µ − , whose branching ratios are constrained to the following 95% C.L. intervals [37] 2.82 × 10 −4 < BR(B s → X s γ) < 4.04 × 10 −4 , which allows one to find the allowed region in the m A − t β plane [38]. It turns out that the region with m A < 350 GeV and t β ≥ 25 is excluded by the B s → µ + µ − decay, which is sensitive to the region with large t β and relatively small m A . On the other hand, the region where m A < 350 GeV and t β ≤ 8 is not favored by the B s → X s γ decay, whereas the region with large m A (m A 400 GeV) and small or moderate values of t β is still unconstrained by the experimental data. In our analysis we consider values of t β and m A in the following intervals 1 < t β < 15 and 400 GeV < m A < 750 GeV.

SUSY parameters
The input parameters M 2 and µ are relatively close and hence the charginos are mixed states of bino and Higgsino states. We fix these values according to the current bounds. In particular, we take M 2 = 500 GeV and µ = 400 GeV in our analysis. As far as the CP -violating phase angle Φ µ we take cos(Φ µ ) = 1. With these input parameters, the mass for the lightest chargino is mχ+ 1 = 377.9 for t β = 5. Finally, the U (1) Y gaugino mass parameter, M 1 , is fixed via the GUT relation B. Behavior of the decay φ → Zγγ (φ = h, H, Z) in the hMSSM We first present a general discussion of the most interesting features of the decay φ → Zγγ (φ = h, H, Z). As already mentioned, the decay A → Zγγ receives contributions from loops with SM charged fermions, the W gauge boson, the charged scalar boson, charginos, and squarks. On the other hand, only charged fermion and charginos contribute to the φ → Zγγ (φ = h, H) decay. It turns out that the charged scalar contribution to the A → Zγγ decay is considerably smaller than the fermion contribution and so is the W gauge boson contribution, which is negligible in regions close to the alignment limit as it is proportional to c β−α . Also, squark contributions to both A → Zγγ and φ → Zγγ (φ = h, H) decays are suppressed, which means that the only relevant contributions to both decays are those of SM charged fermions and charginos. Since fermion contributions are proportional to the charged particle mass, the main contributions of SM fermions arise from the top and bottom quarks. For low values of t β the most important contribution comes from the top quark, whereas for large t β the bottom quark contribution becomes relevant. As for the chargino contribution, because of the presence of non-diagonal couplings of both the scalar boson φ and the Z gauge boson to charginos, the box diagrams for the decay φ → Zγγ (φ = h, H) can include two distinct charginos into the loops. However, non-diagonal chargino couplings are very suppressed and we can safely neglect this type of contributions, which in fact are absent in the reducible diagrams. Another point worth mentioning is that the calculation of chargino box diagrams seems to be too cumbersome as the coupling Aχ i +χ j is both scalar and pseudoscalar, as shown in Table II, however, the scalar (pseudoscalar) part vanishes in the SUSY limit. Therefore the corresponding Lorentz structures for the chargino contributions are identical to those induced by SM charged fermions.
In THDMs the decay A → Zγγ can have a considerable enhancement [12] in the region of the parameter space where m A > m Z + m H , when the decaying CP -odd scalar boson is heavy enough to decay as A → ZH, with the CP -even scalar boson decaying subsequently into a photon pair H → γγ. A similar argument goes for the enhancement of the decay H → Zγγ in the scenario with m H > m Z + m A . However, neither scenario is possible in the MSSM as the masses m H and m A are not independent (they are almost degenerate indeed) so the decays A → ZH and H → AZ are not kinematically allowed. Although the decay A → Zh is not kinematically forbidden for very heavy m A , it is very suppressed as its invariant amplitude is proportional to c β−α , which vanishes in the alignment limit. We thus do not expect any enhancement of the decay φ → Zγγ (φ = H, A), let alone the decay h → Zγγ. The latter has a negligibly small branching fraction with no considerable devation from the SM prediction, of the order of 10 −11 , and so we refrain from analyzing it more detailed. We thus focus our study on the decays A → Zγγ and H → Zγγ.
Below we present the numerical analysis of the φ → Zγγ (φ = H, A) decay. The Passarino-Veltman scalar function were evaluated via the LoopTools routines [39,40]. For the sake of completeness we also calculate the main tree-level and one-loop level induced two body-decays into SM particles, such as φ →bb,tt, γγ, gg, and Zγ. Notice that decays such as A → Zh, H → ZZ and H → W W have negligible rates in regions of the parameter space close to the alignment limit. For comparison purposes we also calculate the three-body decays φ → Ztt and Zbb. As far as the decays into SUSY particles are concerned, the only kinematically allowed channel is the one into a neutralino pair φ → χ 0 1 χ 0 1 , which can give a relevant contribution to the total width and was already studied in [41]. Other decays of the heavy neutral scalar bosons into SUSY particles such as squarks or charginos are not kinematically allowed in the region of the parameter space we are considering.

A → Zγγ branching ratio
In the region close to the decoupling limit, the AhZ coupling is highly suppressed, thus the reducible diagram contribution to the A → Zγγ decay arises only from H exchange, which in fact dominates over the box diagrams contribution by almost two orders of magnitude. The behavior of the branching ratio BR(A → Zγγ) as a function of m A in the mass range 350-750 GeV is shown in the top plots of Fig. 4 for two values of t β , whereas the bottom plots show its behavior as a function of t β for two values of m A . We can observe that for t β = 2 (top-left plot), BR(A → Zγγ) is of the order of 10 −8 for m A below 350 GeV, but drops by more than one order of magnitude when thett channel becomes open, though it increases again up to around 10 −8 as m A increases up to 750 GeV. When t β increases up to 15 (top-right plot), there is no significant change in the BR(A → Zγγ) behavior, except around m A = 350 GeV when it is considerably smaller than in the t β = 2 scenario. In the bottom plots of Fig. 4 it is evident that there is very little dependence of BR(A → Zγγ) on t β for m A > 350 GeV. We conclude that the rate for the threebody decay A → Zγγ is well below those for the decays A → Zγ and A → γγ: BR(A → Zγγ)/BR(A → γγ) ∼ 10 −4 and BR(A → Zγγ)/BR(A → Zγ) ∼ 10 −3 when m A = 750 GeV and t β is small. These values decrease by at least one order of magnitude for m A around 400 GeV.

H → Zγγ branching ratio
We now perform an analysis of the behavior of the H → Zγγ branching ratio. As in the case of the A → Zγγ decay, box diagrams give a smaller contribution to the H → Zγγ decay than the ones arising from the reducible diagrams, which in this case can be mediated by a virtual CP -odd scalar boson and the Z gauge boson. However, Z-mediated reducible diagrams give a negligible contribution since its amplitude is proportional to c β−α . Thus the main contribution to this decay arises from the reducible diagram with H exchange via the triangle loops with SM This fact becomes evident in the bottom plots of Fig. 5, where we observe that BR(H → Zγγ) decreases by almost one order of magnitude as t β increases from 2 to 15. As far as the dominant decays are concerned, thett (bb) channel is the dominant one for small (large) t β . The only decay into SUSY particles, the neutralino channel, can also be relevant once it becomes open. In fact for m A = 400 GeV, the rate of theχ 0 1χ 0 1 channel surpasses those of both thett andbb channels in a narrow region around t β = 6, which becomes wider, about 5 t β 11, for m A = 750 GeV.

V. CONCLUSIONS
In this work we have presented an explicit calculation of the decays of the CP -odd and CP -even scalar bosons A → Zγγ and H → Zγγ in the framework of the MSSM, focusing on the so-called hMSSM scenario. These decays have been computed previously in the context of THDMs, where the main contributions arise from loops with the top and bottom quarks. In the MSSM there are new contributions from squarks and charginos. The latter is comparable to those of the top and bottom quarks for certain values of the parameters of the electroweakino sector still consistent with the experimental bounds. As for the contribution from squarks, it is negligibly small since the squark masses lie in the TeV scale. Although both decays receive contributions from box and reducible diagrams, the dominant contribution arises from the latter. We have presented numerical results for the corresponding branching ratios in the context of the decoupling limit and considered values for the free parameters m A and t β still consistent with constraints from experimental data. For the A → Zγγ decay, we found that the branching ratio is suppressed by three orders of magnitude compared with the diphoton channel. In general BR(A → Zγγ) increases as m A increases, though it is not very sensitive to t β . For m A = 750 GeV, BR(A → Zγγ) is of the order of 10 −8 although it can drop up to 10 −10 for small values of m A . As far as the H → Zγγ decay is concerned, its branching ratio exhibits a similar behavior to that of the A → Zγγ decay, although it decreases as t β increases, which stems from the fact that the coupling of the CP -odd Higgs boson to charginos is more suppressed than that of the CP -even Higgs boson. It was found that BR(H → Zγγ) is of the order of 10 −7 for m A = 700 GeV and t β = 2. For these values of the free parameters, BR(H → Zγγ) is close to BR(H → Zγ). As far as the MSSM contribution to the h → Zγγ decay, it is negligible and the corresponding branching ratio is of the order of 10 −11 . Although the studied decays have very suppressed branching ratios, the AZγγ and HZγγ couplings can be of interest in the study of heavy Higgs boson production at a linear collider working in the γγ mode.
FIG. 6: Feynman rules involving only SM particles. All the 4-momenta are incoming and the Lorentz structures are given by Γ αµν (k1, k2, k3) = (k1 − k2) ν g αµ + (k2 − k3) α g µν + (k3 − k1) µ g αν and Σ αβµν = 2g αβ g µν − g αµ g βν − g αν g βµ . The coupling constants read where the indices i, j stand for the chargino generation, and the coupling constants are given by In Fig. 7 we present the Feynman rules for the couplings of the neutral Higgs bosons to fermions and charginos, whereas the ones for the couplings of scalar bosons and squarks to gauge bosons, which are obtained by expanding the covariant derivative in terms of the physical fields, are shown in Fig. 8. Note that the AV V (V = W, Z) and HhZ couplings are absent due to CP conservation.
Finally, as for the trilinear couplings involving scalar bosons and sfermions, which arise from the scalar potential, the corresponding Feynman rules are  Tables I and II.   TABLE II: Constants for the couplings of the neutral Higgs bosons to charginos. We have used the short-hand notation sa = sin a and ca = cos a.
where φ = h, H and the coupling constants g φH − H + are shown in Table I. In the notation of the third generation sfermions, the constants g φqiqj read with the matrices C φqq summarizing the couplings of the Higgs bosons to the squark eigenstates. In particular, the Htt and hqq coupling constants are whereas the couplings of the CP -odd scalar boson to sfermions are  Table I.
Appendix B: A → Zγγ and φ → Zγγ decay amplitudes We now present the form factors F i (i = 1 . . . 7) of Eq. (27) in terms of Passarino-Veltman scalar function. The contribution of SM charged fermions is identical to that arising in type-II THDMs, which was already presented in Ref. [12]. We thus content ourselves with presenting the contributions of charginos and sfermions. The latter only contribute to reducible diagrams.
We first make the following definitions We also use the already defined kinematic variables s, s 1 , and s 2 [Eqs. (32)- (34)] and introduce the auxiliary variables ∆ ij and X defined as

A → Zγγ decay form factors
The box diagrams give the following contributions to the form factors and It is evident that box diagram amplitudes are free of ultraviolet divergences as they are free of two-point Passarino-Veltman scalar functions.
As for the reducible diagrams, they only contribute the form factor F 1 , which is defined in Eq. (28) and receive the following partial contributions from fermions f , the W gauge boson, the charged scalar boson H ± , charginosχ + i and squarksq i : with φ = h, H. The three-point scalar function C 0 (s, m 2 ) can be written in terms of elementary functions as follows with f (x) being given by x < 1. (B10)