Probing the Dark Matter of Three-loop Radiative Neutrino Mass Generation Model with the Cherenkov Telescope Array

We investigate the prospect of detecting the Dark Matter (DM) candidate in the three-loop radiative neutrino mass generation model extended with large electroweak multiplets of the Standard Model (SM) gauge group, at the future imaging atmospheric Cherenkov telescope known as the Cherenkov Telescope Array (CTA). We find that the addition of such large electroweak multiplets leads to a sizable Sommerfeld enhanced annihilation of the DM with $O(\text{TeV})$ mass, into the SM gauge bosons which results in continuum and line-like spectra of very high energy (VHE) gamma-rays, and therefore becomes observable for the CTA. We determine the viable models by setting the upper limit on the $SU(2)_{L}$ isospin of the multiplets from the partial-wave unitarity constraints and the appearance of low-scale Landau pole in the gauge coupling. Afterwards, by considering the continuum VHE gamma-rays produced from the DM annihilation at the Galactic center, we probe the parameter space of the model using the sensitivity reach of the CTA.


Introduction
In recent years, the Imaging Atmospheric Cherenkov Telescopes (IACT) have not only opened new avenues for the ground-based very high energy gamma-ray astronomy [1][2][3] but also offered a testing ground for the Dark Matter (DM) of the universe [4,5]. The parameter space of the DM which has the electroweak quantum numbers and mass in the GeV-TeV range has been probed by currently operating major IACTs: High Energy Stereoscopic System (H.E.S.S.) [6][7][8][9][10], Major Atmospheric Gamma Imaging Cherenkov Telescope (MAGIC) [11,12] and Very Energetic Radiation Imaging Telescope Array (VERITAS) [13]. The Cherenkov Telescope Array (CTA) which is an ongoing international development project for the next-generation IACT, will have the capability to observe gamma-ray of energy from 20 GeV to at least 300 TeV over a large area and wide range of view (up to 10 • ) with more than 100 telescopes located in northern and southern hemispheres. It will allow CTA to achieve a sensitivity about a factor 10 better than current instruments such as H.E.S.S., MAGIC or VERITAS [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28]. The DM charged under the Standard Model (SM) gauge group and with the mass in the multi-TeV range can produce highly energetic diffuse or monochromatic gamma-rays (based on the annihilation channels) which will be within the observational reach of the CTA. Therefore in this paper, we have explored the detection possibility of the TeV DM candidate of the Three-loop Radiative Neutrino Mass Model, known as the Krauss-Nasri-Trodden model, at the CTA.
The Krauss-Nasri-Trodden (KNT) model [29] is one of the early models of radiative neutrino mass generation (for a comprehensive review, please see [30]) which ties the origin and smallness of the neutrino mass with the DM of the universe. The additional Beyond Standard Model (BSM) fields of the model are two single charged singlet scalars, S + 1 , S + 2 and three SM singlet RH neutrinos, N R i , i = 1, 2, 3 with masses lie in the GeV-TeV range. There is a Z 2 symmetry, {S + 2 , N R i } → {−S + 2 , −N R i } which not only omits the tree-level Dirac mass term, L L N R i H, where L L and H are the SM left-handed lepton doublet and Higgs doublet respectively, but also stabilizes the lightest singlet RH neutrino N R 1 to play the role of DM. Besides, the KNT model can be generalized [31] by replacing S + 2 with Φ having integer isospin and hypercharge, Y = 1 and N R i with F i that has integer isospin and Y = 0 under SM gauge group which leaves the three-loop neutrino mass topology invariant. In the generalized KNT model, the lightest neutral fermion component, F 0 1 is the viable DM candidate. Such replacements in KNT model with large electroweak multiplets have been studied for triplet [32], 5-plet [33] and 7-plet [34] cases. In [35], it has been shown that for most of the DM mass range from 1 to 50 TeV, the Sommerfeld enhanced annihilation cross-section to SM gauge bosons in 5-plet and 7-plet are already within reach of H.E.S.S. and the future CTA because more charged component fields of the fermion multiplet contribute to the enhancement and thus lead to a larger annihilation cross-section. Therefore, to probe viable DM candidate in the class of generalized KNT model, we are left with only the minimal model with singlet fermion and next-to-minimal triplet model where DM is the neutral component of a fermion triplet.
The DM of the minimal KNT model, being the singlet under the SM gauge group, has different DM parameter space [36] than the triplet or 5-plet and 7-plet models which have dominant gauge interactions controlling the viable DM parameter space. For this reason, in our CTA sensitivity study, we have addressed the triplet DM because it shares the same region of parameter space as the 5-plet and 7-plet, and left the singlet case for a separate analysis.
The article is organized as follows. In section 2, we present the model and calculate the possible upper bound on the SU (2) L isospin of the electroweak multiplets in the generalized KNT model. In section 3, we describe the Sommerfeld enhanced DM annihilation cross-section, and the gamma-ray flux originated from the DM annihilation in the triplet KNT model. Section 4 delineates the working principle of the CTA, its instrument response functions and possible sources of background, and lastly the observation region we consider in our study. In section 5, we outline our methods to determine the expected gamma-ray counts from the DM signal and backgrounds, and use the Likelihood analysis to set the upper limit on the DM annihilation cross-section in the KNT model. Finally, we conclude in section 6. In appendix A, we present the detailed calculations relevant for our partial-wave unitarity constraints.

The Model
The BSM fields in the generalized KNT model, charged under SM gauge group, SU (3) c × SU (2) L × U (1) Y are Complex scalars: S + 1 ∼ (0, 0, 1), Φ ∼ (0, j φ , 1), and Real fermions: where j φ and j F are integer isospin of SU (2) L . For example the triplet model would contain, Φ ∼ (0, 1, 1) and F 1,2,3 ∼ (0, 1, 0). The generalized KNT Lagrangian that has the additional Z 2 symmetry, is given by, where, c denotes the charge conjugation and dot sign, in shorthand, refers to appropriate SU (2) contractions. Also L α and e Rα are the LH lepton doublet and RH charged leptons respectively and Greek alphabet α stands for the generation index. Moreover, [F ] αβ = f αβ and [G] iα = g iα are 3 × 3 complex antisymmetric and general complex matrices, respectively. H is the SM Higgs doublet. Finally, V (H, Φ, S 1 ) denotes the scalar potential. The mass splittings in fermionic component fields are zero at tree-level and only receive O(100 MeV) splittings due to the radiative correction after the electroweak symmetry breaking [37][38][39]. On the other hand, the mass splittings among component fields of the scalar multiplet Φ is controlled by λ Hφ2 (Φ † .H).(H † .Φ) term of the scalar potential after electroweak symmetry breaking and the splittings, ∆m allowed by the electroweak precision observables only lead to ∆m 2 ij /M 2 0 ∼ 10 −3 if the invariant mass of the scalar multiplet is, M 0 = 10 TeV, and for M 0 > 10 TeV, the ratio becomes even smaller as shown in [40]. Therefore, the component fields of both scalar and fermion multiplets are almost degenerate at the TeV scale. Consequently, such near-degeneracy of fermion component fields does not diminish the Sommerfeld enhancement of DM annihilation.
As we have seen that the large isospin multiplets of the SU (2) L are allowed in the generalized KNT model, one could ask the upper bound on the electroweak quantum numbers for such large multiplets. In the following sections 2.1 and 2.2, we have used the arguments from the Partial-wave unitarity and appearance of the low-scale Landau pole to set an upper limit on the isospin in the generalized KNT model.

Partial-Wave unitarity constraints
To set the partial-wave unitarity constraints on the electroweak quantum numbers of large multiplets in the KNT model, we have looked into the tree-level scattering of the component fields of the fermion multiplets into the SM gauge bosons, F F → V V , i.e. into W W, Z Z, γ γ, γ Z. Similar analysis with SM gauge bosons as the final states was done in [41]. Moreover, the upper bound on the isospin for the scalar multiplet using same final states, was addressed by [42].
Let us consider the 2 → 2 scattering amplitude in momentum space with initial and final 2-particle states |i and |f , respectively, as where, s is the center-of-mass energy. Here T captures the interaction part of the S-matrix, S = 1+iT . From Eq. (3) by using the Jacob-Wick expansion [43][44][45], the corresponding partial-wave amplitude of total angular momentum J is, in terms of the helicities of the initial (λ i 1 , λ i 2 ) and final (λ f 1 , λ f 2 ) states. In addition, β(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2yz − 2zx is the Kallen function. Moreover, a factor of 1/ √ 2 should be multiplied at the right hand side of Eq. (4) for any identical pairs of particles either in the initial or final states. The unitarity condition on the S-matrix, S † S = 1, implies that where the inequality arises because the sum over n is restricted to 2-particle states only. If we consider elastic scattering for which, f = i, we have from Eq.(5), where the sum over k = i are taken over all of the possible 2-particle inelastic channels. By writing, |a J | 2 = (Rea J ) 2 + (Ima J ) 2 , we have, Now that a J ii lies on the unitary circle, the left-hand side of Eq. (7) is bounded by 1/4, we can have in absence of any inelastic channel, Therefore, the possible maximal bound one can set on the real part of an inelastic channel k is, In practice, one can consider the full transition matrix at tree-level that connects all possible 2-particle states and impose the bound Eq. (8) or (9) on the largest eigenvalue of the matrix, |Re a J f i |. In our analysis, we have considered the initial and final states with total electric charge to be zero, F For such channels, the bound Eq. (9) implies that the isospin of the fermionic multiplet in the KNT model (see Appendix A for detailed calculation) has to be, j F ≤ 8 (one generation) and j F ≤ 7 (three generations) .
One could also include the elastic channels, F i F i → F i F i in the coupled-channel matrix. But the helicity amplitude, with helicity transition 0 → 0, in the t-channel via photon exchange has a singularity at θ = 0. Such presence of the Coloumb singularity in the elastic channel renders the partial-wave unitarity analysis somewhat ineffective. Therefore, to put more concrete bound on the electroweak quantum numbers, one has to look into the appearance of the low scale Landau pole in the SM gauge couplings in the presence of large electroweak multiplets.

Appearence of the Landau pole in gauge coupling
The beta-functions of the SM gauge couplings are modified when the larger electroweak fermionic or scalar multiplets are added to its matter content. If the isospin and the hypercharge of the ith fermionic and j-th scalar multiplets are (j F i , Y F i ) and (j S j , Y S j ), respectively, the one-loop betafunctions for the SU (2) L and U (1) Y gauge couplings g and g Y , respectively become, Here, − 19 6 and 41 6 in the brackets are the SM contribution to g and g Y beta-functions, respectively. κ F i = 1, 1 2 if the i-th fermionic multiplet contains Dirac and Weyl fermions, and η S j = 1, 1 2 if the j-th scalar multiplet contains complex and real scalars, respectively. In addition, the Dynkin index, T (R) of the fermionic or scalar multiplet, R with isospin, j is given by, T (R) = j(j + 1)(2j + 1)/3 and its corresponding dimension as D(R) = 2j + 1. As the last two terms of Eq.(11) give a positive contribution to the beta function of g, for large enough isospin, they can overcome the negative SM contribution and make it positively large. As a consequence, such large positive beta-function will lead to the appearance of the Landau pole for the SU (2) L gauge coupling at a scale much lower than the Planck scale. As we have not observed non-perturbative weak interaction at energy scale few orders above the electroweak scale, the appearance of such low scale Landau pole put a severe constraint on the size of the electroweak multiplet.
For the KNT model, the additional three fermionic multiplets with electroweak quantum number, (j, 0) and the scalar multiplet with (j, 1) will contribute to Eq. (11), and depending on the mass scale where the multiplets are integrated in during the renormalization group running; we tabulate the Landau pole for g in Table 1.   From Table 1, we can see that in the KNT model, if the electroweak multiplets are integrated in at M X = M Z , the Landau pole of the coupling g appears at O(TeV) for the 5-plet (j = 2) and at 440 GeV for the 7-plet (j = 3). If M X is set to higher energies, the Landau poles for 5-plet and 7-plet are also being pushed to higher energies but not as high as the Planck scale as we can see for the triplet case (j = 1). Since we have not seen the signature of new physics in the electroweak scale at energies accessible by the LHC, the appearance of the low-scale Landau pole of the SM gauge coupling for electroweak multiplets with j ≥ 2 has made the singlet (j = 0) and the triplet as the viable BSM addition for the KNT model. One could improve the analysis of the size of the electroweak multiplet and the Landau pole by using the two-loop beta functions of the gauge couplings as done in [46,47]. However, we refrain from using the two-loop beta functions because firstly, we do not know the UV completion of the KNT model and secondly, the KNT Yukawa coupling g iα enters in the two-loop beta functions of the gauge couplings, and eventually the two-loop Landau pole will also depend on these couplings. For this reason, the size of the electroweak multiplets will not be connected to the appearance of the low scale Landau pole in a straightforward way as in the one-loop case. However, even if we consider the UV completion of the model with large electroweak multiplets in a Grand Unified Theory (GUT) setup, it would require either enormous representations of the minimal GUT like SO(10) or a much bigger GUT group where the model can be embedded in its fundamental or adjoint representations, both of which are not theoretically appealing.

Sommerfeld Enhanced DM Annihilation
In the triplet KNT model, the lightest neutral component of the fermion triplet, F 0 1 is the viable DM candidate. The Sommerfeld enhancement of the DM annihilation to SM gauge bosons takes places when the DM is non-relativistic, v DM c and the SM gauge bosons follow m W,Z m DM . In this limit, the exchange of W and Z bosons between triplet component fields will lead to Yukawa potentials and γ exchange will lead to a Coulomb potential, respectively which in turn significantly modifies the wavefunction of the incoming DM states and enhances the annihilation cross-sections. The calculation of Sommerfeld enhanced DM annihilation cross-section is a well-studied subject, so here we follow the prescriptions given in [48]. Nevertheless, we briefly review them to set up our notation.
As the Sommerfeld enhancement is considered for 2 → 2 processes, we first define 2-particle states which consist of incoming component fields of fermion triplet. For the DM (co)annihilation channels where the final states consist of W ± , Z and γ bosons, the only relevant 2-particle states are CP-even states with total electric charges Q = 0, ±1, ±2. In the case of DM annihilation in the galaxy halo at present times, only 2-particle states with Q = 0 are applicable.
In general, if the three fermionic triplets of the KNT model are almost mass degenerate, M F 1 ∼ M F 2 ∼ M F 3 , the 2-particle states will contain component fields not only from one multiplet but also from different ones. But, as the hierarchical mass spectrum of F i is consistent with the best-fit experimental values of the neutrino mixing angles and mass-squared differences [49], it allows us to have M F 1 M F 2,3 and M F 1 to be in O(TeV). It enables us to define 2-particle state vector made out of component fields of F 1 as follows, The modification of the wavefunction which generates the Sommerfeld enhancement, is determined by solving the radial Schrodinger equation with effective potential, where r is the magnitude of the relative distance between two component fields in their center-ofmass frame, the kinetic energy of the incoming DM states, i.e. |ii = F 0 The wavefunction Ψ jj ,ii gives the transition amplitude from |ii states to |jj states in the presence of effective potential, V . The double indices ii , jj and kk run over the states of the 2-particle state vector defined in Eq. (13).
We primarily focus on the S-wave annihilation so we set l = 0 and have Here, is the momentum associated with the 2-particle state, |jj and d jj = m j + m j − 2M F 1 denotes the mass differences between DM and other states of the multiplet. Q kk is the electric charge associated with state |kk . Also, α W = α and n W = 1 for W boson exchange and α Z = α/ cos 2 θ W and n Z = 1/ cos θ W for Z boson exchange. Finally, f jj ,kk is the group theoretical factor associated with SU (2). Now by using dimensionless variables defined as we re-write the coupled radial Schrodinger equations as where the dimensionless momentum,k 2 At large x, Ψ jj ,ii behaves as Ψ jj ,ii ∼ T jj ,ii e ik jj x where T jj ,ii is the transition amplitude provided the effective potential is dominated by Yukawa potential. Now if the annihilation matrix for final state f is given by Γ (f ) jj ,ii , the annihilation cross section is, where c = 2 for |F 0 1 F 0 1 state as it consists of identical fields. The Sommerfeld enhancement factor is then given by S F 0 where σ 0 is the tree-level annihilation cross-section.

Dark Matter Constraints and Parameter Space
The relic density of DM in the universe measured by Planck Collaboration as Ω DM h 2 = 0.120 ± 0.001 (68% C. L.) [50] sets an important constraint for KNT triplet DM. In the case of F 0 1 , the standard thermal freeze-out scenario involves the SM gauge interactions and the KNT yukawa interactions given by g iα terms in Eq. (2).
The DM (co)annihilation into the SM final states, controlled by the SM gauge interactions, take place via S-wave and P -wave channels. Moreover, the S-wave DM (co)annihilation cross-sections receive non-negligible Sommerfeld enhancement. In contrast, the DM annihilation into charged leptons that involves the KNT Yukawa coupling g iα , has S-wave and P -wave channels. However, both of them are suppressed by the light charged lepton mass and the velocity of the DM, respectively. Because of the large multiplicity and unsuppressed Sommerfeld enhanced S-wave contribution, the DM (co)annihilation controlled by gauge interaction are more dominant compared to that into charged (Right) The Spin-independent cross section of triplet DM, F 0 1 and nucleon interaction (orange horizontal line). The shaded region is excluded by XENON1T (2018) data [51] leptons involving the KNT Yukawa coupling. Therefore, essentially the gauge interaction determines the relic density of the triplet KNT DM.
Another constraint comes from the DM direct detection. The spin-independent cross-section for the Majarana DM contained in the electroweak multiplet of integer isospin j does not depend on the DM mass, and is given by [37], where, M Nucl is the mass of the target nucleus, f parametrizes nucleon matrix element as n| q m q qq|n = f m n nn and from lattice result, f = 0.347131 [52]. On the other hand, the spin-dependent cross-section is suppressed by the mass of the DM which is of the O(TeV). Therefore, we can see that the constraint from the DM relic density sets the mass of the DM in the triplet KNT model as 2.55 TeV (without SE) and 2.94 TeV (with SE), respectively, if the standard thermal freeze-out mechanism is considered. But, as pointed out in [40], the nonthermal decay of the scalar, φ + , can produce the DM F 1 within the KNT model, although some fine-tuning in the masssplittings and couplings are needed to get the correct relic density for a wider range of DM masses. However, the fine-tuning can be avoided if the KNT model is coupled to an extended dark sector which would assist the non-thermal production of the DM. Hence, we would be able to probe the DM masses up to 100 TeV in our CTA sensitivity study.

Gamma-ray flux from DM annihilation
The gamma-rays from the DM annihilation consist of the prompt gamma-rays and the secondary gamma-rays. The prompt gamma-rays are produced directly from the DM annihilation, or the decays of the SM final states originated in the annihilation. On the other hand, the secondary gamma-rays come from the inverse Compton scattering (ICS) of e ± produced in DM annihilation, again directly or from decays of the SM final states, with the ambient photon background, i.e. mainly coming from Cosmic Microwave Background photons, dust rescattered light and starlight. As our analysis focus on the diffuse prompt gamma-ray flux coming from the Sommerfeld enhanced DM annihilation into gauge bosons (W + W − ) with energy in the TeV range, we neglect the effect of secondary gamma-ray emission from DM annihilation for simplification.
The differential prompt gamma-ray flux from the DM annihilation for a given DM mass, m F 1 , in the generalized KNT model is, where σv is the corresponding velocity averaged annihilation cross-section for m F 1 , dN f /dE γ is the gamma-ray spectra per annihilation for the annihilation channel, F 0 1 F 0 1 → f and B f is the corresponding branching ratio. The integration over the solid angle, ∆Ω and line-of-sight (l.o.s), s of the squared DM mass density, ρ DM is called the astrophysical J-factor. Here, r is the distance between the Galactic Center (GC) and the point in space characterized by galactic coordinates, (b, l), and it is given as r(s, θ) = (r 2 + s 2 − 2r s cos θ) 1/2 , where r is the distance between the Earth and Galactic Center and θ is the angle between the line-of-sight and the axis connecting the Earth and the GC as shown in Fig. 3. In addition, θ and the galactic coordinates (b, l) are related as cos θ = cos b cos l.

Gamma-ray Spectrum from DM Annihilation
As mentioned already, the dominant DM annihilation channels in the triplet KNT model are those with the SM gauge bosons as final states. Hence, two types of gamma-ray spectra arise from SM gauge boson final states, i) Continuum spectra coming from W + W − and ZZ, and ii) line-like spectra , respectively. Besides, the DM annihilation to γ + X can also give line-like spectra at the end-point of photon energy. The sensitivity studies for IACTs involving the continuum and line-like photon spectra require different strategies. In our study, we focus on the continuum spectra given by the W + W − final states as a representative channel. We have used PPPC4DMID [53] to calculate the gamma-ray spectra coming from the W + W − final states which include the electroweak (EW) corrections [54] that become dominant at O(TeV) energies. From Fig. 4, we can see the decrease in the gamma-ray spectra per annihilation with the mass of the DM.
In passing, let us make some remarks on the EW corrections implemented in the PPPC4DMID. At energies much higher than the electroweak scale, the dominant EW corrections for a typical 2 → 2 process arises from the Sudakov double logarithms which are of the form, δ DL ∼ α 4π L log 2L s M 2 W at L-th loop order [55][56][57]. Here, α is SU(2) coupling and for the non-relativistic DM, the center of mass energy is s ∼ 4M 2 DM . As we can see from

Astrophysical J-factor
The astrophysical J-factor plays an essential role in determining the gamma-ray flux from DM annihilation, and therefore, one needs precise information of DM density profile in the Milky Way, especially around the GC which is the focus of this work. The N-body cosmological simulations are often used to parametrize the DM density profile of the Milky Way. The two most frequently used DM density profiles, that assume spherical symmetry, are the Navarro-Frenk-White (NFW) profile and the Einasto profile, where ρ s and r s are the characteristic density and scale radius, respectively. The NFW profile behaves like r −1 at the GC whereas the Einasto profile does not and is smaller in magnitute. The α is the shape parameter that determines the steepness of Einasto profile at the neighbourhood of the GC.
The numerical values of these parameters are for NFW profile, (ρ s , r s ) = (0.184 GeVcm −3 , 24.42 kpc) and for Einasto profile, (ρ s , r s , α) = (0.033 GeVcm −3 , 28.44 kpc, 0.17), respectively. The value of J-factor is subject to considerable uncertainties because there is no adequate observational data and precise simulation of our galaxy with all its baryonic and DM content. Moreover, the DM density profile depends on the nature of the DM itself (for a review, please see [67]). Therefore, the J-factors used in different studies vary significantly [68][69][70]. Hence, we have used the NFW profile, because of its cuspy nature at GC, as the representative in our CTA sensitivity study to detect the DM of the KNT model.

TeV DM at the CTA
Gamma-rays being electrically neutral, do not deviate from their paths due to the galactic magnetic field when propagating through the galaxy. Moreover, gamma-rays with energies ranging from 100 GeV to 100 TeV, have minimal absorbance in the interstellar medium. So, when these very high energy (VHE) gamma-rays enter the Earth's atmosphere, they undergo collisions with atmospheric molecules and generate an air shower of secondary particles known as Extensive Air Shower. These shower particles descend to Earth with almost the speed of light, and due to the Cherenkov effect, these ultra-relativistic charged particles create a faint blue light in the air which typically lasts for a few nanoseconds. Most of this Cherenkov light is emitted at altitudes ranging between 5 to 15 km, and it propagates down to the ground level as a quasi-planar, thin disk of Cherenkov photons orthogonal to the shower axis. It may cover about 50000 m 2 of an area on the ground. By placing arrays of IACTs within the projected Cherenkov light pool, it is possible to detect the air shower provided that the mirror area of the telescope is large enough to catch enough photons. However, it is challenging to reconstruct the exact geometry of the air shower in space with the observation from a single telescope. Hence, multiple telescopes are deployed to take the image of a separate shower from different points which leads a stereoscopic reconstruction of the shower geometry. The images captured by IACTs after removing the backgrounds contributions shows the track of the air shower, which points back to the celestial origin of the incident gamma-rays, and eventually makes it possible to determine the location of its source in the sky along with its spectral and spatial properties. As we have seen in section 3 that when the DM of the triplet KNT model has mass in O(TeV) range, its Sommerfeld enhanced annihilation into SM gauge bosons can produce either diffuse or sharp VHE gamma-rays. Therefore, the IACTs are ideally suited to detect such TeV DM in the galaxy.
There are three major currently operational IACTs: H.E.S.S., MAGIC and VERITAS which have performed well within their capabilities and have discovered more than hundreds of VHE gamma-ray sources. However, the decisive point of the scientific performance of CTA will be its ability to survey the sky over broad energy ranges from 20 GeV to 300 TeV with better angular resolution, energy reconstruction and sensitivity, which is one order of magnitude better than existing IACTs. Besides, to provide the full sky coverage, CTA will be installed at two different sites: one in the northern hemisphere at La Palma (Canary Islands, Spain) and another in the southern hemisphere at Paranal (Chile). The southern observatory (known to be as CTA South) will study the GC, and its northern counterpart (known to be as CTA North) will survey extra-galactic objects. As the search for the DM at the GC will be one of the key focus of the CTA [26], we have chosen to study its detection possibility of the TeV DM of the triplet KNT model and probe its the parameter space.

CTA Instrument Response Function
The instrument response function (IRF) can be considered as the area times probability that a photon with a given set of input parameters is detected as an event with a set of observables (reconstructed parameters). For CTA, IRF is the product of effective area (A eff ), point spread function (PSF), energy resolution and dispersion, and background rate as a function of energy. The key features of these quantities are as follows.
The effective area is the area within which the CTA can observe air showers. It depends on the energy of the primary gamma-ray and the offset angle φ (the angle between the array's normal direction and the actual source position). Moreover, the point spread function (PSF) of CTA IRFs represents the spatial probability distribution of reconstructed event positions for a point source. Besides, in short, the energy dispersion of the CTA gives the resolution between the actual and reconstructed gamma-ray energies of events recorded by it. Finally, the background rate is the residual cosmic-ray background rate per solid angle (here, square degree) as a function of reconstructed gamma-ray energy in the CTA. As the CTA is not constructed yet, the CTA consortium has carried out several Monte Carlo simulations [71][72][73] to compute them. We specifically make use of South z20 50h in CTA-Performance-prod3b-v2-FITS.tar.gz IRF (southern site with zenith angle 20 • and exposure time 50 hours) throughout our entire study. Let us point out a few factors for choosing this IRF. Due to its privileged location, the CTA southern site is most favorable for surveying the γ sources in GC with energy ranges from 20 GeV to 300 TeV. Since our study is focused on the DM signal from the GC, we considered the IRFs designed for the southern site. Since it is assumed that the total operation time of CTA per year will be 1000 hours, we, therefore, expect that a longer exposure period is necessary to probe a possible a DM source region. From Fig. 6, we can see the variation of the effective area, that will enter into our computation of gamma-ray counts due to DM or background in the CTA, with respect to gamma-ray energy and off-set angle.

CTA Backgrounds
In this section, we summarize the dominant background sources at the CTA considered for our study of detecting the DM signal from the KNT model.

Cosmic Rays:
The Cosmic Rays (CRs) are highly energetic protons and nuclei. Apart from the VHE gamma-rays, CR also trigger air showers, which are 10 3 times larger than that of γ rays, and will act as the dominant background for the CTA. The CR, during their flight from their origin, undergo deflection due to galactic magnetic fields; only particles with sufficiently high energies can reach the Earth's atmosphere. While entering the Earth's atmosphere, CR protons undergo inelastic collisions with atmospheric air molecules, causing mixed hadronic and electromagnetic air showers 6 . However, because of the large transverse momentum transfer due to hadronic interactions, the shower components are broad and exhibit irregular patterns compared to those of electromagnetic showers initiated by the primary gamma-ray.
The shape of the captured image induced by the Cherenkov light from the air shower initiated by VHE gamma-ray can be well approximated by an ellipse. But the image gathered from the CR induced air shower has a distorted elliptical shape which can allow one to discriminate events generated due to gamma-rays and cosmic rays. Nevertheless, a small fraction of images remain indistinguishable from the signal gamma-rays, which constitute the irreducible background for the CTA.
2. Galactic Diffuse Emission: Other than CR, Galactic Diffuse Emission (GDE) acts as a potential background for CTA, which originates from the interaction of CR with interstellar molecules and atomic gas filling the Galactic plane. Observation of Fermi-LAT found that up to energy scales of 100 GeV, the GDE is dominated by π 0 decay, inverse Compton scattering and Bremsstrahlung radiation [76]. The contribution due to GDE becomes significant when the region of interest becomes close to the GC, which is relevant for our case since we are studying the possible DM sources near the GC.
Modelling TeV range GDE for existing ground-based IACTs is a challenging task due to the uncertainties associated with the estimated background in the telescopes' field of view. It is also impossible to find a signal-free region in the sky by apriori consideration [74]. Besides, the data available from current experiments is limited to have a proper modelling of GDE at the 100 TeV energy scale. As a consequence, we have taken a simplied approch following [18,19,25] and left more realistic GDE modelling for a future work. We incorporate the GDE in our analysis using the P7V6 model of Fermi team 7 which fits extremely well for 50 MeV and 500 GeV energy range, and use a power-law extrapolation of P7V6 data to 100 TeV. The H.E.S.S. observation of GDE from Galactic ridge for Galactic coordinates |b| < 0.3 • and |l| < 0.8 • [75,77] is not taken into account as it falls within the region that we have excluded from our study.

Region of Interest
Our observation method is based on the Multi-RoI morphological analysis [19,25]. We assumed that all telescopes are pointing toward the GC with Galactic longitude and latitude b with coordinates (0, 0). The Region of Interest (RoI) is divided into five concentric circles, each with width 1 • such that the outermost circle has radius 5 • , as seen in Fig. The expected number of photon counts for signals (DM) and backgrounds (CR, GDE) in each circular RoI has been computed simultaneously. It is known 6 There is a certain minimum threshold energy required in order to trigger the Cherenkov light shower by VHE-γ and CR. For electrons, muons, pions and protons the threshold energies are 21 MeV, 3.4 GeV, 5.6 GeV and 38 GeV, respectively. 7 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html that the central part of the GC is populated with several astrophysical sources including Sagittarius A* (SgrA*), emitting gamma-rays. In order to exclude photon counts coming from this region we disregard the central part of the GC by introducing a rectangular patch with 0.3 • < b < −0.3 • and −5 • < < 5 • within the RoI. The corresponding value of the J-factor, considering the NFW profile, is shown in Fig. 7.

Expected Gamma-ray Counts at the CTA
In this section, we obtain the expected counts for the dark matter, cosmic rays and galactic diffuse emission using the methods described below. The expected differential count for each RoI i and energy bin j for the source X is calculated using the relation, where the differential gamma-ray flux for the component X, dφ X γ dEγ dΩ is weighted by the effective area A eff (E γ ,p), point spead function, P SF (p,p ) and energy dispersion, E disp (E γ , E γ ,p). Here, (E γ ,p) and (E γ ,p ) denote the energy and direction of the actual and reconstructed gamma-rays, respectively. The P SF (p,p ) which is a probability distribution function of the angular separation between the actual and reconstructed gamma-rays, is mostly relevant for a point source in the sky but the sources we are considering, i.e. DM, CR or GDE, are extended so, the PSF can be well approximated by a deltafunction. Moreover, the energy dispersion, E disp is important if we are looking for a specific energy of the gamma-ray, which applys to the line-like spectrum from DM annihilation. But, in our case, the resulting gamma-ray spectrum from DM annihilation has continuum spectrum so again weighting the differential gamma-ray flux with gaussian-like energy dispersion function is not important. Finally the photon count is obtained using where, T obs is the total observation time, and the integration region ∆E j denotes the j-th energy bin. In our analysis, we divide our photon counts into twenty logarithmically spaced energy bins and five RoIs as described in section 4.4. Also, the observation time is set to T obs = 100 hr. For the DM annihilating into W + W − final states, the associated photon count is computed using Eq.(21), (24) and (25) for the NFW profile .
Again, we compute the expected photon count for the GDE using Eq. (24) and (25) from our simplified gamma-ray flux for the GDE which is a power-law extrapolation up to 100 TeV from the Fermi P7V6 model, The expected gamma-ray count associated with the CR is evaluated using ctools version 1.6.2 [78] which is a software toolkit designed for data analysis of the CTA and other IACTs. The ctools consists of a set of binary executable C++ and python tools for performing the necessary data analysis where each of these tools needs a set a parameter like telescope pointing directions, radius of FoV, pixelation, time period of observation, calibration datebase etc. Equipped with IRF-South Z20 50h, we use four components of the ctools in our analysis: ctobssim, ctbin, ctcubemask and csresspec following the sequence, Although csresspec inspects the spectral fit residual in ctools, we use it to extract the count per energy bin. Besides, ctcubemask of the ctools-1.6.2 does not have an option to carry out rectangular masking of the region: 0.3 • < b < −0.3 • and −5 • < < 5 • . To do this task, we exclude the region from longitude −5 • to 5 • using overlapping circles of radius 0.3 • . The output from ctools of the expected gamma-ray counts due to the cosmic-rays from first four RoIs are shown in Fig. 8. In addition, as seen in Fig. 9, for both 1st (innermost) and 5th (outermost) RoI, the Sommerfeld enhancement in the DM annihilation allows its count rate to be comparable with the CR count. Moreover, we can see that in the case of CR and GDE, at low energy, the expected count rates are smaller in the 5th RoI than that in the 1st RoI. On the contrary, when the gamma-ray energy becomes higher, the count rate in the 5th RoI becomes larger compared to the 1st RoI. In Fig. 6 one can see that for low energy, the effective area of the CTA is small and increases with the gamma-energy. Besides, the incident gamma-ray at a large angle with respect to the CTA telescope axis needs to have comparatively higher energy to initiate Cherenkov light that is detectable at the CTA. For this reason, we see a low count rate at low energy in the 5th RoI. On the other hand, with increasing angle towards the outer RoIs, the combination Ω i × A eff becomes large at higher energies, and therefore shows higher gamma-ray counts at high energies for the CR and GDE. In the case of DM, the expected count rate is proportional to the J-factor, which slightly decreases towards outer RoI, as seen in Fig. 7. For that reason, we see that the expected count rate for DM is large to some extent in the 1st RoI compared to the 5th one.  The impact of the Sommerfeld enhancement on the DM annihilation is seen again in Fig. 10 where we consider the total expected count rate in each RoI for the DM, CR and GDE. In this case, the variation of the expected count rate with respect to gamma-ray energy at each RoI is averaged out.

Likelihood Analysis
In this section, we combine the expected gamma-ray counts coming from DM, CR and GDE backgrounds at the CTA, and determine its prospect to detect the DM in the triplet KNT model.
We divide our photon counts into twenty logarithmically spaced energy bins and five RoIs in a procedure similar to [18,25] using the binned Poisson likelihood analysis. The counts are labeled with µ X ij where X indicates the source -DM, CR, or GDE, and the lower indices indicate the count in the i-th energy bin and j-th RoI. The number of expected counts is then given by where, we have introduced the rescaling parameters {R CR , R GDE }. We may now write the likelihood function as the product of the independent Poisson distributions associated with each energy bin and  RoI: where, n ij represents the number of observed counts. However, we do not simulate individual Poisson realizations of the observed counts n ij , but instead obtain an Asimov data set [79]. In this procedure, the entire observed count is considered to be the sum of CR and GDE backgrounds only, and it can be found by setting µ DM ij → 0 and R CR/GDE i → 1 in Eq. (27). Furthermore, we account for the systematic uncertainties in the signals by introducing a new set of Gaussian distributed nuisance parameters α ij , with mean 1 and standard deviation σ α . This approach results in a modification of Eq.(28) so that the new likelihood function is containing the nuisance parameters we denote collectively by the set θ = {α, R CR , R GDE }.
We now use maximum likelihood estimates of the nuisance parameters separately in each of the above three cases to evaluate the test statistic where the subtext MLE indicates maximum likelihood estimates of the underlying parameters and we suppress the implicit dependence on the DM mass, m F 1 and cross-section, σv in µ for notational clarity (Recall that µ DM ij ∼ σv ). λ is asymptotically equivalent to a χ 2 distribution with one degree of freedom (since µ in Eq.(30) is still a free parameter). We bound the rescaling factors in Eq. (27) so that 0.5 ≤ R CR i ≤ 1.5 and 0.2 ≤ R GDE i ≤ 5 (in accordance with [18]) and we require α ≥ 0 in order to keep probabilities non-negative in Eq. (29). While we formally maximize the likelihood function, in practice we minimize the negative log-likelihood function: where we neglect constants and terms involving n ij which ultimately cancel in evaluating λ and hence do not affect our analysis. Furthermore, we reduce the set of nuisance parameters by first analytically maximizing log L with respect to α ij : which, upon imposing the above mentioned constraint α ij ≥ 0, gives: which we substitute back into Eq. (29) and Eq. (31). Note that α ij → 1 in the limit σ α → 0. With the explicit dependence on α eliminated from the likelihood function in Eq.(29), we may regard the nuisance parameters as the reduced set θ = {R CR , R GDE }. At this stage we wish to find an upper limit of the 95% confidence interval for σv . Towards this end, we numerically calculate λ: first we find the denominator in Eq.(30), thereby obtaining µ MLE ; next, we modify the numerator by gradually increasing σv until λ reaches 2.71 (Recall that χ 2 1 ∼ Z 2 with Z ∼ N (0, 1), and hence for finding the upper limit of the 95% confidence interval, we require Z 2 ≈ 1.645 2 ≈ 2.71). We subsequently use this value of σv for the 95% upper limit in the exclusion line for one particular fixed mass m F 1 . We then repeat this method for masses ranging from 1 TeV to 100 TeV, thereby obtaining an exclusion line for σv . Finally, we repeat this entire process in the three distinct scenarios mentioned above, i.e. with the following backgrounds: CR, CR and GDE, as well as CR, GDE and 1% systematic uncertainties. All of these results are illustrated in Fig. 11.

Result and Discussion
In Fig. 11, we can see that the CTA can probe the DM of the triplet KNT model with the mass up to 25.7 TeV. The DM mass range 1.5 TeV ≤ M F 1 ≤ 4.25 TeV, where the direct detection constraint sets the lower limit, can be excluded the CTA observation of continuum gamma-rays coming DM annihilation at the GC. However, the CTA sensitivity derived here is subjected to several issues discussed in the following.
In our study, one prominent uncertainty which influences the expected gamma-ray count from the DM annihilation at the CTA and may degrade its DM detection sensitivity, is the astrophysical J-factor in Eq. 21. As the astrophysical J-factor depends on the DM density profile of the Milky Way, the precise information of its density profile near the GC is essential to reduce the uncertainty that translates into the J-factor used for calculating gamma-ray flux. We have used the NFW profile which is cuspy at the GC to calculate the corresponding J-factors for our annular RoIs centering the GC. However, one could use joint profile likelihood analysis, as done, for example, in [64], to set upper limits on both the DM annihilation cross-section and the J-factor values associated with the considered RoIs. Also, the proper treatment of the background in RoIs contributes to the CTA sensitivity for DM detection. In our case, we consider the gamma-rays from the DM annihilation, and the backgrounds CR and GDE simultaneously from the same the region of the sky. One could have considered two regions of the sky, one where the count from DM signal is expected to be high, known as the ON region, and adjacent DM signal-free but background dominated region taken as OFF. The statistical significance of the DM signal is determined in this scenario by contrasting counts from ON and OFF region. However, in such ON-OFF method, it is not always possible to define the signal-free region beforehand for extended sources like the DM, which could lead to systematic uncertainties in the counts and eventually incorrect substantial significance. The method we use is less prone to systematic uncertainties as all the counts are taken for the same region.
Moreover, our extrapolation of the GDE up to 100 TeV from the Fermi P7V6 model given by Eq. (26) is simplified as the GDE can have non-trivial spatial and spectral properties at such high energies. For this reason, we have derived CTA sensitivity with and without the GDE. If we consider only the CR as the background, we get the most stringent limit (red dashed line in Fig. 11) on the DM annihilation cross-section into W + W − final states in the triplet KNT model at the GC from our analysis. The inclusion of the GDE reduces the sensitivity as seen from the green dashed line in Fig.  11. Finally, as expected, the systematics of 1% further degrades the CTA sensitivity as illustrated by the blue dashed line in Fig. 11.

Conclusion and Outlook
In this work, we have studied the sensitivity of the upcoming Cherenkov Telescope Array to detect the DM candidate, with mass at the TeV range, of the KNT model annihilating into W + W − at the Galactic Center. In the following, we summarize our key results in order.
As shown in section 2, the KNT model can be generalized with large fermionic multiplets, F i of three generations and scalar multiplet, Φ, all of them having the same integer isospin, j of SU (2) L , and hypercharges, Y F i = 0 and Y Φ = 1, respectively without changing the topology of three-loop neutrino mass generation. Here, we set the upper limit on the value of the SU (2) L isospin for the KNT model using the bounds from the partial-wave unitarity on the coupled channels, F (Q) i F (−Q) → W + W − , ZZ, γγ, γ Z. The bound on the SU (2) isospin of fermion multiplet is j F ≤ 8 for one generation and j F ≤ 7 for three generations, respectively. The inclusion of the elastic channels in the coupled-channel analysis makes the partial-wave unitarity bound somewhat ineffective because of the Coulomb singularity arising from photon exchange in F Therefore, a more refined bound on the SU (2) isospin comes from the appearance of the low-scale Landau pole in the gauge coupling in the presence of large electroweak multiplets. From Table 1, we can see that the j = 3 case is almost ruled out and j = 2 is already in tension with the non-observation of non-perturbative SU (2) coupling and the new physics appearing at the LHC. It leaves the minimal KNT with singlets and the triplet KNT with j = 1 as the only viable models within the class of such three-loop radiative neutrino mass generation models.
However, the region of parameter space associated with the dark matter of the singlet and triplet KNT models are different. Because the DM of the singlet KNT does not have any SM gauge interaction, and its thermal freeze-out and annihilation at the present day are controlled by the KNT Yukawa couplings which directly connects them to the neutrino mass generation and the related low-energy constraints. For this reason, the viability of the singlet KNT DM with O(TeV) mass is left for a future study. On the other hand, the DM of the triplet KNT model, being charged under the SM gauge group, can have large Sommerfeld enhancement, and its parameter space is determined mostly by the gauge interactions when the DM mass is in the TeV range.
Besides, the SE annihilation of such heavy DM can produce detectable VHE gamma-rays at the IACTs. We show that the DM of the triplet KNT model annihilating into W + W − at the GC can be probed up to 25.7 TeV by the future CTA experiment considering NFW profile and 100 hours of observation. However, there is some room for improvement of our sensitivity analysis. The triplet KNT model also has SE annihilation channels into γγ and γ Z which would give line-like signature in the gamma-ray spectrum around, E γ ∼ M F 1 . Although the statistical analyses related to the continuum and line-like gamma-ray spectra have some differences, a proper combination of both analyses can lead to a better CTA sensitivity. Also, the consideration of the morphological analyses, adopted in [18,26,27] would result in an improved sensitivity of the CTA to detect the triplet KNT DM.
In addition, one can consider the dwarf spheroidal galaxies (dSphs), that is expected to contain up to O(10 3 ) times more mass in DM than in visible matter and have low astrophysical gamma-ray background than the GC. Such characteristics of dSphs make them an ideal environment to search for the DM at the CTA [26]. Furthermore, the electrons produced from DM annihilation (mostly as secondary) can generate synchrotron radiation in the presence of the magnetic field, which could be observed as a diffuse radio emission in the present and future radio observatories like Square Kilometer Array (SKA) [82], and the sensitivity to detect such radio signal from DM in the dSphs is promising [80,81]. Therefore, a combined sensitivity study of the CTA and SKA to detect the DM in the dSphs will enable us to probe further the parameter space of the generalized KNT model [83] and other variants [84,85] of the three-loop radiative neutrino mass generation models. strument response functions provided by the CTA Consortium and Observatory, see http://www.ctaobservatory.org/science/cta-performance/ (version prod3b-v2) for more details. Also, this research made use of ctools, a community-developed analysis package for Imaging Air Cherenkov Telescope data. ctools is based on GammaLib, a community-developed toolbox for the scientific analysis of astronomical gamma-ray data.
where, the first two entries correspond to fermion's helicity, ± 1 2 and last two entries denote gauge boson's helicity, ±1. Also, t 3 is the value associated with the component field, F (Q) for the SU (2)'s T 3 generator. Here, as the hypercharge of the fermion multiplet is zero, the electric charge of its component field follows, Q = t 3 . Moreover, V ± (t 3 ) = 1