Zhu Luo ^1,
, Cheng-feng Cai ^1,
, Zhao-feng Kang ^2,
, Zhao-huan Yu ^1,3,,
, Hong-hao Zhang ^1,,
,

Corresponding author: Zhao-huan Yu, yuzhaoh5@mail.sysu.edu.cn

Corresponding author: Hong-hao Zhang, zhh98@mail.sysu.edu.cn

1.School of Physics, Sun Yat-Sen University, Guangzhou 510275, China
2.School of Physics, Huazhong University of Science and Technology, Wuhan 430074, China
3.ARC Centre of Excellence for Particle Physics at the Terascale, School of Physics, The University of Melbourne, Victoria 3010, Australia
Received Date:2018-11-08
Available Online:2019-02-01
Abstract:We confront the perturbativity problem in the real scalar quintuplet minimal dark matter model. In the original model, the quintuplet quartic self-coupling inevitably hits a Landau pole at a scale ~10¹⁴ GeV, far below the Planck scale. In order to push up this Landau pole scale, we extend the model with a fermionic quintuplet and three fermionic singlets which couple to the scalar quintuplet via Yukawa interactions. Involving such Yukawa interactions at a scale ~10¹⁰ GeV can not only keep all couplings perturbative up to the Planck scale, but can also explain the smallness of neutrino masses via the type-I seesaw mechanism. Furthermore, we identify the parameter regions favored by the condition that perturbativity and vacuum stability are both maintained up to the Planck scale.

HTML

--> --> -->

One of the biggest mysteries of nature, dark matter (DM) has drawn much attention from astrophysicists, cosmologists, and particle physicists. Among various guesses at the identity of the DM particle, the most extensively studied class of DM candidates is weakly interacting massive particles (WIMPs), because they can naturally explain the observed DM relic abundance via the thermal production mechanism in the early Universe [1-4]. WIMP models can be easily constructed by introducing a dark sector with electroweak SU(2)_L multiplets. Introducing one nontrivial SU(2)_L multiplet leads to the so-called minimal dark matter (MDM) models [5-18], which only involve the minimal content of new fields. Introducing more than one SU(2)_L multiplet results in a richer phenomenology, but the models are much more complicated [19-42].
The philosophy of the MDM models is to extend the standard model (SM) in a minimal way to involve dark matter [5]. For this purpose, a fermionic or scalar SU(2)_L × U(1)_Y multiplet in a representation $(n, Y)$ is introduced. The potential DM candidate would be the electrically neutral component that should be the lightest new state. If the dimension of the SU(2)_L representation $n$ is large enough to forbid dangerous decay operators, this neutral state would be able to play the role of the DM particle without imposing an artificial $Z_2$ symmetry. In other words, the stability of DM can be protected by an accidental $Z_2$ symmetry due to a large dimension of the SU(2)_L representation and the Lorentz invariance. As proposed in the original paper [5], the minimal choices are a quintuplet ($n=5$) for the fermionic case [6, 8, 14-16] and a septuplet ($n=7$) for the scalar case [9, 11, 12, 16, 17]. Adding new electroweak multiplets will push up the SU(2)_L gauge coupling $g_2$ when it runs to high scales. Nevertheless, the MDM model with a quintuplet fermion or a septuplet scalar can keep $g_2$ perturbative up to ~10¹⁹ GeV [5] (~10¹⁶-10¹⁷ GeV [43]) based on one-loop (two-loop) $\beta$ functions.
Scalar MDM models are quite different from fermionic ones, since scalars will bring in more coupling terms. Such complexity has caused the neglect of a dangerous decay operator for the septuplet scalar model in the original consideration: the dimension-5 operator $\Phi^3 H^{\dagger } H$ will induce DM decay at loop level [43]. Therefore, the accidental $Z_2$ symmetry in the septuplet scalar model is not totally strict. Nonetheless, one can always introduce an artificial $Z_2$ symmetry to make the model work again, but considering $n=7$ would not be special any more. In this case, discussing a triplet ($n=3$) or quintuplet ($n=5$) real^①) scalar multiplet with $Y=0$ would be more economic. The real scalar triplet model has been studied thoroughly in Refs. [44-46], while the real scalar quintuplet model is less explored and will be the main topic of this paper.
In such scalar MDM models, scalar coupling terms may lead to another problem. Solutions to the renormalization group equations (RGEs) show that the scalar self-interaction couplings will go to infinity, i.e., a Landau pole (LP) will show up, at an energy scale far below the Planck scale [17, 46, 47]. With two independent septuplet self-interaction terms, the real scalar septuplet model hits a Landau pole at a scale around 10⁸ GeV if the DM particle mass is fixed to satisfy the observed relic abundance. In our previous work [17], we attempted to push up the LP scale via introducing Yukawa couplings between the scalar septuplet and extra fermionic multiplets. A bonus of these extra fermions is to explain the smallness of neutrino masses through the type-III seesaw mechanism [48]. We found that such a setup can push up the LP scale to ~10¹⁴ GeV at best.
On the other hand, a real scalar quintuplet lives in a smaller representation and has only one independent self-interaction term. Consequently, the quintuplet couplings should evolve slower and reach a Landau pole at a higher scale. If extra fermionic multiplets are introduced, we may even push the LP scale above the Planck scale. Besides, such fermions could be used to explain the tiny neutrino masses via the type-I seesaw mechanism [49-52]. In this work, we will explore these possibilities in the real scalar quintuplet MDM model. For completeness, we will also discuss the constraint from the observed DM relic abundance, the bounds from direct and indirect detection experiments, and the stability of the electroweak vacuum.
The paper is organized as follows. In Section 2, we introduce the quintuplet MDM model, and discuss its phenomenological constraints and the LP scale. In Section 3, we study an extension with extra fermions for pushing up the LP scale, and discuss the constraints from perturbativity and vacuum stability. Conclusions and discussions are given in Section 4. Appendix A gives the $\beta$ functions and initial values of SM couplings, while Appendix B gives the detailed calculation of the Sommerfeld enhancement effect.

2

-->

2.1.Model details

In the real scalar quintuplet MDM model, the dark sector only involves a real scalar quintuplet $\Phi$ with $Y=0$, which can be expressed as

$ \Phi=\frac{1}{\sqrt{2}}(\Delta^{(2)}, \ \Delta^{(1)}, \ \Delta^{(0)}, \ \Delta^{(-1)}, \ \Delta^{(-2)})^{\rm{T}}. $

(1)

The self-conjugate condition implies $(\Delta^{(Q)})^\ast=\Delta^{(-Q)}$. The real scalar $\Delta^{(0)}$ is a viable DM candidate. The gauge covariant derivative of $\Phi$ is

$ D_\mu \Phi=\partial_\mu \Phi-ig_2W^a_\mu\tau^a \Phi, $

(2)

where $\tau^a$ are generators for the SU(2)_L representation 5:

$\begin{split} &{\tau ^1} = \left( {\begin{array}{*{20}{c}} {} & { - 1} & {} & {} & {}\\ { - 1} & {} & { - \sqrt 6 /2} & {} & {}\\ {} & { - \sqrt 6 /2} & {} & {\sqrt 6 /2} & {}\\ {} & {} & {\sqrt 6 /2} & {} & 1\\ {} & {} & {} & 1 & {} \end{array}} \right), \\ &{\tau ^2} = \left( {\begin{array}{*{20}{c}} {} & i & {} & {} & {}\\ { - i} & {} & {\sqrt 6 i/2} & {} & {}\\ {} & { - \sqrt 6 i/2} & {} & { - \sqrt 6 i/2} & {}\\ {} & {} & {\sqrt 6 i/2} & {} & { - i}\\ {} & {} & {} & i & {} \end{array}} \right), \end{split}$

(3)

$ \tau^3 = {\rm{diag}}(2, 1, 0, -1, -2). $

(4)

Thus, the covariant kinetic term for $\Phi$ can be expanded as

$ \begin{aligned} \cal{L}_{\rm{kin}}&=(D_\mu \Phi)^{\dagger } D^\mu \Phi \\ &=\frac{1}{2}(\partial_\mu\Delta^{(0)})^2+\sum_{Q=1}^2\left(\partial_\mu\Delta^{(Q)}\right)\left(\partial^\mu\Delta^{(-Q)}\right) \\ &\quad+\sum_{Q=1}^2\left(QeA^\mu+Qg_2c_{\rm{W}} Z^\mu\right)\Delta^{(-Q)}i\overleftrightarrow{\partial_\mu}\Delta^{(Q)}\\ &\quad -g_2 \left[W^{+, \mu}\left(\sqrt{2}\Delta^{(-2)}i\overleftrightarrow{\partial_\mu}\Delta^{(1)} +\sqrt{3}\Delta^{(-1)}i\overleftrightarrow{\partial_\mu}\Delta^{(0)}\right)+{\rm{h.c.}}\right] \\ &\quad \!+\!\left(e^2A_\mu A^\mu\!+\!g_2^2c_{\rm{W}}^2 Z_\mu Z^\mu\!+\!2eg_2c_{\rm{W}} A_\mu Z^\mu\right) \sum_{Q=1}^2Q^2\Delta^{(Q)}\Delta^{(-Q)} \\ &\quad +g_2^2 W^+_\mu W^{-, \mu}\bigg[3(\Delta^{(0)})^2+5\Delta^{(1)}\Delta^{(-1)} +2\Delta^{(2)}\Delta^{(-2)}\bigg] \\ &\quad -g_2^2\bigg\{ W^+_\mu\left(s_{\rm{W}} A^\mu + c_{\rm{W}} Z^\mu\right) \left(\sqrt{3}\Delta^{(0)}\Delta^{(-1)}+3\sqrt{2}\Delta^{(1)}\Delta^{(-2)}\right) \\ &\quad + W^+_\mu W^{+, \mu} \bigg[\frac{3}{2}(\Delta^{(-1)})^2 -\sqrt{6}\Delta^{(0)}\Delta^{(-2)}\bigg] +{\rm{h.c.}}\bigg\}, \end{aligned} $

(5)

where $s_{\rm{W}}\equiv\sin\theta_{\rm{W}}$, $c_{\rm{W}}\equiv\cos\theta_{\rm{W}}$, and $\overleftrightarrow{\partial_\mu}$ is defined as $F\overleftrightarrow{\partial_\mu} G=F\partial_\mu G-G\partial_\mu F$.
In order to protect the stability of $\Delta^{(0)}$, we require that $\Phi$ is odd under a $Z_2$ symmetry, while all SM fields are even. The scalar potential is constructed by $\Phi$ and the SM Higgs doublet $H$. Since the operator $\Phi^{\dagger }\tau^a\Phi$ vanishes due to the self-conjugation condition, the general form of the potential respecting the $Z_2$ symmetry is given by only five independent terms:

$ \begin{split} V &= {\mu ^2}{H^\dagger }H + {m^2}{\Phi ^\dagger }\Phi + \lambda {({H^\dagger }H)^2} + {\lambda _2}{({\Phi ^\dagger }\Phi )^2} \\ &\quad+{\lambda _3}({H^\dagger }H)({\Phi ^\dagger }\Phi ). \end{split} $

(6)

Therefore, this model just brings in two couplings, $\lambda_2$ and $\lambda_3$, and one mass parameter $m$ as new free parameters. We assume that the vacuum expectation value (VEV) of the Higgs field is nonzero, while the VEV of $\Phi$ remains zero. Then the minimization of the potential implies two conditions, $\mu^2<0$ and $m^2-{\lambda_3\mu^2}/{(2\lambda)}\geqslant 0$. As in the SM, the VEV of the Higgs doublet is $\langle H\rangle=(0, {v}/{\sqrt{2}})^{\rm{T}}$ with $ v\equiv$ $\sqrt{{-\mu^2}/{\lambda}}$=246.22 GeV.
After the Higgs field acquires a VEV, the $\lambda_3(H^{\dagger } H)$ $ (\Phi^{\dagger }\Phi)$ term contributes equally to the masses of all the $\Phi$ components. Therefore, at the tree level all components are degenerate with a shifted mass $m_0$, given by

$ m_{0}^2=m^2+\frac{\lambda_3}{2}v^2. $

(7)

Electroweak one-loop corrections break this degeneracy, making $\Delta^{(1)}$ and $\Delta^{(2)}$ slightly heavier than $\Delta^{(0)}$. When $m_0 \gg m_Z$, the mass difference between $\Delta^{(Q)}$ and $\Delta^{(0)}$ is [5]

$ m_Q - m_0 = Q^2 \Delta m, $

(8)

where $\Delta m=\alpha_2 m_W\sin^2({\theta_{\rm{W}}}/{2})\simeq $167 MeV, with $\alpha_2\equiv g_2^2/{(4\pi)}$.
Vacuum stability (VS) sets a stringent constraint on the model. The philosophy is that the potential should remain bounded from below as the couplings evolve to high energies. The VS conditions can be obtained by means of the copositive criteria [53]:

$ \lambda\geqslant 0, \quad \lambda_2\geqslant 0, \quad \lambda_3+2\sqrt{\lambda\lambda_2}\geqslant 0. $

(9)

2 -->

2.2.Experimental constraints

The observation of the DM relic abundance sets a constraint on the $\Delta^{(0)}$ mass $m_0$. Assuming DM is thermally produced in the early Universe, its relic abundance can be expressed as [54]

$ \begin{split} {\Omega _{{\rm{DM}}}}{h^2} \simeq & \frac{{1.07 \times {{10}^9}\;{\rm{Ge}}{{\rm{V}}^{{\rm{ - 1}}}}}}{{J({x_{F}})\sqrt {{g_ * }} {M_{{\rm{Pl}}}}}}\;{\rm{with}}\;\; \\ J({x_{F}}) =& \int_{{x_{F}}}^\infty {\frac{{\langle {\sigma _{{\rm{eff}}}}v\rangle }}{{{x^2}}}} {\rm d}x, \end{split} $

(10)

where $M_{\rm{Pl}}$ is the Planck mass, $x_{F}$ is the freeze-out parameter, $g_\ast$ is the total number of effectively relativistic degrees of freedom, and $\left<\sigma_{{\rm{eff}}}v\right>$ the effective thermally averaged annihilation cross section accounting for the coannihilation effect.
For $m_0\gg m_h$, annihilation and coannihilation into gauge and Higgs bosons in the s-wave are dominant, leading to the following result [9]:

$ \langle\sigma_{{\rm{eff}}}v\rangle\simeq \frac{66\pi \alpha_2^2}{5m_0^2}+\frac{\lambda_3^2}{80\pi m_0^2}. $

(11)

We take $x_{F}\simeq 25$ and $\sqrt{g_\ast}\simeq 10.33$ for $T\sim \cal{O}$ (TeV), and calculate the prediction of the relic abundance in the quintuplet MDM model.
A more proper treatment is to consider the Sommerfeld enhancement (SE) effect for DM freeze-out, following the strategy in Refs. [6, 15]. Annihilation and coannihilation channels in the dark sector are categorized by the electric charges of the two-body states. Enhancement factors can be computed in various categories via numerically solving the Schr?dinger equations for the two-body states. The inclusion of the SE effect would increase the effective annihilation cross section, and hence reduce the relic abundance for fixed model parameters. Details of the calculation are summarized in Appendix B.
Assuming thermally produced $\Delta^{(0)}$ particles in the early Universe fulfill the observed relic abundance, $\Omega_{\rm{DM}}h^2=$$0.1193\pm 0.0014 $ [55], the $\Delta^{(0)}$ mass $m_0$ can be constrained within a narrow range, depending on the coupling $\lambda_3$. As shown in Fig. 1, a slightly larger $m_0$ is favored for a larger $|\lambda_3|$. If the SE effect is not taken into account, the favored $m_0$ is ~5 TeV. After including the SE effect, the favored value^②) is $m_0$~7 TeV.
Other constraints come from DM direct and indirect detection experiments. Direct detection uses the response to DM-nucleon scattering. The only tree-level diagram of the spin-independent $\Delta^{(0)}$-nucleon scattering process is mediated by the Higgs boson, but its cross section is suppressed by $m_0^2$ and thus is insignificant unless $\lambda_3$ is very large. On the other hand, according to the discussions in Ref. [58], the gauge loops lead to a DM-nucleon cross section of ~2×10^?46 cm² for a fermionic quintuplet. As the gauge interactions of the scalar quintuplet are similar to those of the fermionic one, we may expect that the cross section for the scalar case would also be around this value^③). Current direct detection experiments, such as PandaX-II [59] and XENON1T [60], have set an upper bound of ~5×10^?45 cm² on the DM-nucleon cross section for a DM particle mass of 7 TeV. Thus, the scalar quintuplet model can evade current direct searches, but should be well tested in near future experiments.
For indirect detection of DM annihilation in space, the dominant process is $\Delta^{(0)}\Delta^{(0)}\rightarrow W^+W^-$. The SE effect on such a process is even more significant, since the velocities of Galactic DM particles are much lower than those at the freeze-out epoch. Details of the related calculation are also given in Appendix B. Using the 95% CL exclusion limit on annihilation cross section obtained by the combined MAGIC and Fermi-LAT $\gamma$-ray observations of dwarf satellite galaxies, assuming DM totally annihilating into $W^+W^-$ [56], we find that indirect detection experiments have given very stringent constraints, as shown by the red regions in Fig. 1. Nonetheless, three windows around $m_0 \simeq $7 TeV, 8-10 TeV, and 13-14 TeV survive. As a result, the parameter regions suggested by the observation of relic abundance are still available.
2 -->

2.3.$\beta$ functions and the Landau pole scale

Renormalization group equation (RGE) evolution of couplings are determined by $\beta$ functions. In Appendix A, we list the $\beta$ functions in the SM, as well as initial values of SM couplings. When the renormalization scale $\mu$ goes above a scale $\Lambda_{\rm{s}} \sim m_0$, the effect of the real scalar quintuplet should be involved. Hereafter we will adopt $\Lambda_{\rm{s}} $=10 TeV, as suggested by the observed relic abundance. At one-loop level, the real scalar quintuplet MDM model introduces new terms in the $\beta$ functions of the SU(2)_L gauge coupling $g_2$ and the Higgs quartic coupling $\lambda$, while the $\beta$ functions of the other SM couplings do not change. Here we list the $\beta$ functions that are relevant in the following analysis: [47]

$ \begin{split} {\beta _{{g_1}}} &= \beta _{{g_1}}^{{\rm{SM}}},\;\;\;{\kern 1pt} {\beta _{{g_2}}} = \beta _{{g_2}}^{{\rm{SM}}} + \frac{1}{{16{\pi ^2}}}\frac{5}{3}g_2^3,\;\;\;{\kern 1pt}\\ {\beta _{{g_3}}} &= \beta _{{g_3}}^{{\rm{SM}}},\;\;\;{\kern 1pt} {\beta _{{y_t}}} = \beta _{{y_t}}^{{\rm{SM}}}, \end{split} $

(12)

$ \begin{split} {\beta _\lambda } &= \beta _\lambda ^{{\rm{SM}}} + \frac{1}{{16{\pi ^2}}}\frac{5}{2}\lambda _3^2,\;\\ {\beta _{{\lambda _2}}} &= \frac{1}{{16{\pi ^2}}}\left[26\lambda _2^2 + 108g_2^4 - 72g_2^2{\lambda _2} + 2\lambda _3^2\right], \end{split} $

(13)

$ \begin{split} {\beta _{{\lambda _3}}}= & \frac{1}{{16{\pi ^2}}}\Bigg[ {12\lambda {\lambda _3} + 14{\lambda _2}{\lambda _3} + 4\lambda _3^2 + 18g_2^4} \Bigg.\\ &\Bigg. { - {\lambda _3}\left( {\frac{{81}}{2}g_2^2 + \frac{9}{{10}}g_1^2 - 6y_t^2} \right)} \Bigg]. \end{split} $

(14)

Note that $g_1$ is related to the U(1)_Y gauge coupling $g'$ by $g_1\equiv \sqrt{5/3}g'$. The solution to $\alpha_2 = g_2^2/(4\pi)$ for $\mu > \Lambda_{\rm{s}}$ is just

$ \alpha_2^{-1}(\mu) = \alpha_2^{-1}(m_Z) -\frac{b_2^{\rm{SM}}}{2\pi}\ln\frac{\Lambda_{\rm{s}}}{m_Z} -\frac{b_2^{\rm{s}}}{2\pi}\ln\frac{\mu}{\Lambda_{\rm{s}}}, $

(15)

where $b_2^{\rm{SM}}=-19/6$ and $b_2^{\rm{s}}=b_2^{\rm{SM}}+5/3=-3/2$ are coefficients corresponding to $\beta_{g_2}^{\rm{SM}}$ in Eq. (A1) and $\beta_{g_2}$ in Eq. (12), respectively.
Below we analyze the Landau pole scale. The large coefficients of the $\lambda_2^2$ and $g_2^4$ terms in the beta function of $\lambda_2$ drive $\lambda_2$ to divergence at high energy scales. At one-loop level, $g_1$ and $y_t$ do not enter the beta function of $\lambda_2$, while $\lambda_2$ and $\lambda_3$ do not contribute to the beta functions of $g_1$ and $y_t$. Therefore, $g_1$ and $y_t$ remain small at high energy scales. For an approximate analysis, we can safely neglect $g_1$ and $y_t$ in Eqs. (13) and (14). Thus, the RGEs for $\lambda$, $\lambda_2$, and $\lambda_3$ become

$ \frac{{\rm d}\lambda}{{\rm d}t'} = 24\lambda^2+\frac{9}{8}g_2^4-9g_2^2\lambda+\frac{5}{2}\lambda_3^2, $

(16)

$ \frac{{\rm d}\lambda_2}{{\rm d}t'} = 26\lambda_2^2+108g_2^4-72g_2^2\lambda_2+2\lambda_3^2, $

(17)

$ \frac{{\rm d}\lambda_3}{{\rm d}t'} = 4\lambda_3^2+18g_2^4-\frac{81}{2}g_2^2\lambda_3+12\lambda\lambda_3+14\lambda_2\lambda_3, $

(18)

where $t'\equiv (4\pi)^{-2} \ln(\mu/\Lambda_{\rm{s}})$. For convenience, we define three functions $f_i$ by $\lambda=f_1g_2^2$, $\lambda_2=f_2g_2^2$, and $\lambda_3=f_3g_2^2$, and obtain the equations for them:

$ \frac{{\rm d}f_1}{{\rm d}G} = 24f_1^2+\frac{9}{8}-(9+2b_2^{\rm{s}})f_1+\frac{5}{2}f_3^2, $

(19)

$ \frac{{\rm d}f_2}{{\rm d}G} = 26f_2^2+108-(72+2b_2^{\rm{s}})f_2+2f_3^2, $

(20)

$ \frac{{\rm d}f_3}{{\rm d}G} = 4f_3^2+18-\left(\frac{81}{2}+2b_2^{\rm{s}}\right)f_3+12f_1f_3+14f_2f_3. $

(21)

where $G(t')=(b_2^{\rm{s}})^{-1}\ln[g_2(t')/g_2(0)]$.
In order to eliminate the linear terms, we further define functions $\hat{f}_1$ and $\hat{f}_2$ by the shifts $f_1=\hat{f}_1+1/8$ and $f_2=\hat{f}_2+69/52$. Then we have

$ \frac{{\rm d}\hat{f}_1}{{\rm d}G}=24\hat{f}_1^2+\frac{3}{4}+\frac{5}{2}f_3^2, \quad \frac{{\rm d}\hat{f}_2}{{\rm d}G}=26\hat{f}_2^2+\frac{6471}{104}+2f_3^2. $

(22)

Noting that $\hat{f}_2$ runs much faster than $\hat{f}_1$ and $f_3$, we can simply neglect $f_3$ in the second equation and find its solution as

$ \hat{f}_2(t')=\hat{d}\tan\left[\hat{c}_2\hat{d}G(t')+\tan^{-1}\left(\frac{\hat{f}_2(0)}{\hat{d}}\right)\right], $

(23)

where $\hat{c}_2=26$ and $\hat{d}=\sqrt{6471/(104c_2)}\approx1.55$. The Landau pole is reached when $\hat{c}_2\hat{d}G(t')+\tan^{-1}[\hat{f}_2(0)/\hat{d}]=$$\pi/2 $. Thus, the corresponding scale is

$ \begin{split} \Lambda _{{\rm{LP}}}^{({f_2})} &= {\Lambda _{\rm{s}}}\exp \left[ { - \frac{{2\pi }}{{b_2^{\rm{s}}{\alpha _2}({\Lambda _{\rm{s}}})}}} \right.\\ &\quad\times\left. {\left( {\exp \left( { - \frac{{b_2^{\rm{s}}\pi }}{{{{\hat c}_2}\hat d}}\left[ {1 - \frac{2}{\pi }{{\tan }^{ - 1}}\left( {\frac{{{{\hat f}_2}(0)}}{{\hat d}}} \right)} \right]} \right) - 1} \right)} \right]. \end{split} $

(24)

Setting $\hat{f}_2(0)=-69/52$, which corresponds to $\lambda_2=$0 at $\mu = \Lambda_{\rm{s}}$, we find that the maximal LP scale for $\lambda_2$ is $\Lambda_{\rm{LP}}^{(\lambda_2)}$=5.6×10¹⁴ GeV, which is far below the Planck scale. Such a Landau pole implies that other new physics may exist between the quintuplet mass scale and the Planck scale, rendering all the couplings finite.

In this section, we will attempt to push up the LP scale obtained above. A lesson we can learn from the standard model is that the top Yukawa coupling gives a negative contribution to the self-coupling of the Higgs boson. As the Landau pole is induced by the self-coupling of the quintuplet, it is straightforward to introduce extra fermions with a Yukawa coupling to the quintuplet to shift the Landau pole. Such a motivation leads to the 5-5-1 model studied below.
2

-->

3.1.Yukawa interactions

There are three minimal ways to construct Yukawa interactions with the quintuplet scalar: introducing fermions in $({\bf{1}}, 0)\oplus({\bf{5}}, 0)$, $({\bf{3}}, 0)\oplus({\bf{3}}, 0)$, and $({\bf{4}}, \pm1/2)\oplus\; ({\bf{2}}, \mp1/2)$The first and second options have potential for explaining the tiny neutrino masses via the type-I and type-III seesaw mechanisms, respectively. In order to keep $\Delta^{(0)}$ stable, one of the two fermions participating the Yukawa interaction should be odd under the $Z_2$ symmetry. Particularly, in the second option, the two (3, 0) fermions should be different: one is $Z_2$-odd, and the other one is $Z_2$-even. Consequently, in order to give correct neutrino oscillation properties, we have to introduce at least one more triplet which is $Z_2$-even. On the other hand, in the third option, the $({\bf{4}}, \pm 1/2)$ representation should correspond to a $Z_2$-odd Dirac fermion for avoiding anomalies [61-64]. The $({\bf{2}}, \mp 1/2)$ fermions can just be the SM lepton doublets which are $Z_2$-even. However, this case is less interesting for us, as it cannot explain the neutrino masses.
For these reasons, in this work we only concentrate on the $({\bf{1}}, 0)\oplus({\bf{5}}, 0)$ case. The resulting model is dubbed the “5-5-1” model. Minimally, we introduce a left-handed self-conjugated fermionic quintuplet $\Psi_{\rm{L}}$ and several right-handed fermionic singlets $N_{a, {\rm{R}}}$. It is convenient to use the tensor notation for writing down the interaction terms. The tensor notation can be translated to the familiar vector notation using the following dictionaries of $\Phi$ and $\Psi$:

$\begin{split} &\Phi = \frac{1}{{\sqrt 2 }}\left( {\begin{array}{*{20}{c}} {{\Delta ^{(2)}}}\\ {{\Delta ^{(1)}}}\\ {{\Delta ^{(0)}}}\\ {{\Delta ^{( - 1)}}}\\ {{\Delta ^{( - 2)}}} \end{array}} \right) = \frac{1}{{\sqrt 2 }}\left( {\begin{array}{*{20}{c}} {{\Phi _{1111}}}\\ {2{\Phi _{1112}}}\\ {\sqrt 6 {\Phi _{1122}}}\\ { - 2{\Phi _{1222}}}\\ {{\Phi _{2222}}} \end{array}} \right),\\ &{\Psi _{\rm{L}}} = \left( {\begin{array}{*{20}{c}} {{\Psi _{ + 2, {\rm{L}}}}}\\ {{\Psi _{ + 1, {\rm{L}}}}}\\ {{\Psi _{0, {\rm{L}}}}}\\ {{\Psi _{ - 1, {\rm{L}}}}}\\ {{\Psi _{ - 2, {\rm{L}}}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{\Psi _{1111, {\rm{L}}}}}\\ {2{\Psi _{1112, {\rm{L}}}}}\\ {\sqrt 6 {\Psi _{1122, {\rm{L}}}}}\\ { - 2{\Phi _{1222, {\rm{L}}}}}\\ {{\Phi _{2222, {\rm{L}}}}} \end{array}} \right). \end{split}$

(25)

At the renormalizable level, the Yukawa interactions can be expressed as

$ {\cal{L}}_{\rm{yuk}}= -(y_\nu)_{ab} \overline{\ell_{a, {\rm{L}}}^i}H_iN_{b, {\rm{R}}} -\sqrt{6}y\Phi_{ijkl}\overline{\Psi_{\rm{L}}^{ijkl}}N_{3, {\rm{R}}} + {\rm{h.c.}}, $

(26)

where $\ell_{a, {\rm{L}}}$ denotes the SM lepton doublets. $i, j, k, l=1, 2$ are totally symmetric SU(2)_L indices. $a$ and $b$ are family indices and at least two singlets are required for generating the realistic neutrino mixing. Equation (26) respects the $Z_2$ symmetry, with $\Psi_{\rm{L}}$ being $Z_2$-odd and $N_{a, {\rm{R}}}$ being $Z_2$-even. Consequently, the new terms do not endanger stability of the scalar quintuplet, which has mass far below the fermionic one.
The above Yukawa interactions involve many new parameters. To illustrate the point, we adopt some working simplifications. In the first, we assume that there are three generations of $N_R$, and the third generation has the largest Yukawa coupling, which is close to the $\tau$ Yukawa coupling, i.e., $(y_\nu)_{33}\simeq y_\tau$. Other elements in the $y_\nu$ matrix are much smaller but should be consistent with neutrino oscillation data. Next, we assume all the new fermions, $\Psi_{\rm{L}}$ and $N_{a, {\rm{R}}}$ have the same mass $M_{51}$, which is a characteristic scale of the 5-5-1 model. Thus, the neutrino masses given by the seesaw mechanism are $\sim m_\tau^2/M_{51}$. According to the current cosmological constraint [65], $m_\nu\sim m_\tau^2/M_{51}\lesssim$0.2 eV, which implies that $M_{51}$~10¹⁰ GeV. Hereafter, $M_{51}$=10¹⁰ GeV will be set as a benchmark scale of the 5-5-1 model. In the concrete numerical analysis, we will comment on the situation of deviation from this scale setup. Finally, in the 5-5-1 Yukawa interaction term, i.e., the second term of $\cal{L}_{\rm{yuk}}$, we only consider the coupling to $N_{3, {\rm{R}}}$, neglecting the other two couplings. Thus, we just need to deal with one 5-5-1 Yukawa coupling, $y$. This 5-5-1 Yukawa term can be expanded as

$ \begin{split} &{L_{{\rm{yuk}}}} \supset - \sqrt 3 y\left[{\Delta ^{(0)}}\overline {{\Psi _0}} {N_3}\right . \\ &+\left.\left({\Delta ^{( - 1)}}\overline {{\Psi _{ - 1}}} {N_3} + {\Delta ^{( - 2)}}\overline {{\Psi _{ - 2}}} {N_3} + {\rm{h}}.{\rm{c}}.\right)\right], \end{split} $

(27)

where $\Psi_{-Q}=\Psi_{-Q, {\rm{L}}}+(\Psi_{+Q, {\rm{L}}})^{\rm{c}}$ and $N_3=N_{3, {\rm{R}}}+(N_{3, {\rm{R}}})^{\rm{c}}$. The gauge couplings of the fermionic quintuplet $\Psi$ are given by

$ \begin{split} \mathcal{L}_\Psi=&g_2\big(\sqrt{3}W^+_\mu\overline{\Psi_0}\gamma^\mu\Psi_{-1}+\sqrt{2}W^+_\mu\overline{\Psi_{-1}}\gamma^\mu\Psi_{-2}+{\rm{h.c.}}\big)\\ &-\big(eA_\mu + g_2 c_{\rm{W}} Z_\mu\big) \big(2\overline{\Psi_{-2}}\gamma^\mu\Psi_{-2}+\overline{\Psi_{-1}}\gamma^\mu\Psi_{-1}\big). \end{split} $

(28)

2 -->

3.2.$\beta$ functions and analytic results

The contributions of the quintuplet and singlet fermions at $\mu > M_{51}$ further modify the one-loop $\beta$ functions of the SU(2)_L gauge coupling and scalar couplings:

$ \delta \beta_{g_2} = \frac{1}{16\pi^2}\frac{20}{3}g_2^3, \quad \delta \beta_{\lambda_2}=\frac{1}{16\pi^2}\left(-72y^4+16y^2\lambda_2\right), $

(29)

$ \delta \beta_{\lambda_3} = \frac{1}{16\pi^2}8y^2\lambda_3, \quad \beta_y =\frac{1}{16\pi^2}y\left(19y^2-18g_2^2\right). $

(30)

In these expressions, we have neglected the effect of the neutrino Yukawa couplings $(y_\nu)_{ab}$. This will be justified later.
Note that the $\beta$ function of $g_2$ becomes positive for $\mu>M_{51}$. Thus, one may worry about up to which scale $g_2$ can remain perturbative. By requiring $\alpha_2=4\pi$, we find the non-perturbative scale of $g_2$ in the 5-5-1 model as

$ \Lambda_{g_2}^{\rm{NP}}=M_{51}\left(\frac{\Lambda_{\rm{s}}}{M_{51}}\right)^{b_2^{\rm{s}}/b_2^{\rm{tot}}}\exp\left[\frac{2\pi}{b_2^{\rm{tot}}}\left(\alpha_2^{-1}(\Lambda_{\rm{s}})-\frac{1}{4\pi}\right)\right], $

(31)

where $b_2^{\rm{tot}}=b_2^{\rm{s}}+20/3=31/6$. Setting $\Lambda_{\rm{s}}$=10 TeV, we find that almost any $M_{51}>\Lambda_{\rm{s}}$ will give $\Lambda_{g_2}^{\rm{NP}}>M_{\rm{Pl}}$. Thus, $g_2$ would still be perturbative at the Planck scale.
By solving the RGEs, we obtain the exact values of $g_2(\mu)$ and $y(\mu)$ at $\mu>M_{51}$ as

$ g_2^2(\mu) = g_2^2(M_{51}) \left[1-\frac{1}{8\pi^2}b_2^{\rm{tot}}g_2^2(M_{51})\ln\frac{\mu}{M_{51}}\right]^{-1}, $

(32)

$ y^2(\mu) \!=\! \left(18+b_2^{\rm{tot}}\right)g_2^2(\mu) \left[19\!+\!F_0\left(\frac{g_2^2(\mu)}{g_2^2(M_{51})}\right)^{(18+b_2^{\rm{tot}})/b_2^{\rm{tot}}}\right]^{-1}, $

(33)

where

$ F_0 \equiv \left(18+b_2^{\rm{tot}}\right)\frac{g_2^2(M_{51})}{y^2(M_{51})}-19. $

(34)

$F_0 = 0$ leads to a critical value of $y(M_{51})$,

$ y_{\rm{c}}(M_{51}) \equiv \sqrt{\frac{18+b_2^{\rm{tot}}}{19}} g_2(M_{51}). $

(35)

If $y(M_{51})>y_c(M_{51})$, i.e., $F_0<0$, the 5-5-1 Yukawa coupling $y$ will reach a Landau pole at

$ \begin{split} \Lambda _{{\rm{LP}}}^{(y)} =& {M_{51}}\exp \left[ {\frac{{8{\pi ^2}}}{{b_2^{{\rm{tot}}}g_2^2({M_{51}})}}} \right.\\ &\times\left. {\left( {1 - {{\left[ {1 - \frac{{18 + b_2^{{\rm{tot}}}}}{{19}}\frac{{g_2^2({M_{51}})}}{{{y^2}({M_{51}})}}} \right]}^{b_2^{{\rm{tot}}}/(18 + b_2^{{\rm{tot}}})}}} \right)} \right]. \end{split} $

(36)

Then the condition for $\Lambda_{\rm{LP}}^{(y)}>M_{\rm{Pl}}$ is

$ \begin{split} &{y^2}({M_{51}})<\frac{{18 + b_2^{{\rm{tot}}}}}{{19}}g_2^2({M_{51}})\\ &\quad\times{\left[ {1 - {{\left( {1 - \frac{{b_2^{{\rm{tot}}}g_2^2({M_{51}})}}{{8{\pi ^2}}}\ln \frac{{{M_{{\rm{Pl}}}}}}{{{M_{51}}}}} \right)}^{(18 + b_2^{{\rm{tot}}})/b_2^{{\rm{tot}}}}}} \right]^{ - 1}}. \end{split} $

(37)

Note that the perturbative condition $y\leqslant4\pi$ will give a smaller upper bound on $y(M_{51})$.
As $y$ grows at high scales, the effect of $(y_\nu)_{ab}$ may become important. The $\beta$ function of $y_\nu$ with all Yukawa coupling included is

$ \beta_{y_\nu}=\frac{y_\nu}{16\pi^2}\left[\frac{3}{2}y_\nu^{\dagger } y_\nu+T+\frac{15}{2}y^2-\frac{3}{2}y_e^{\dagger } y_e-\frac{9}{4}g_2^2-\frac{9}{20}g_1^2\right], $

(38)

where $T=3{\rm{tr}}(y_u^{\dagger } y_u)+3{\rm{tr}}(y_d^{\dagger } y_d)+{\rm{tr}}(y_e^{\dagger } y_e)+{\rm{tr}}(y_\nu^{\dagger } y_\nu)$, and $y_u$, $y_d$, $y_e$, $y_\nu$ are the Yukawa coupling matrices for quarks and leptons. For simplicity we consider that only one real element of $y_\nu$ on the diagonal dominates, and denote it $\hat{y}_\nu$. At high scales, $y$ and $g_2$ are large, and the above equation can be approximated as

$ \beta_{\hat{y}_\nu}\simeq\frac{\hat{y}_\nu}{16\pi^2}\left[\frac{5}{2}\hat{y}_\nu^2+\frac{15}{2}y^2-\frac{9}{4}g_2^2\right]. $

(39)

As $\hat{y}_\nu$ contributes to the self-energy of $N_{3, {\rm{R}}}$, it will modify $\beta_y$ by

$ \delta\beta_{y}=\frac{y}{16\pi^2}\hat{y}_\nu^2. $

(40)

Therefore, the growing $\hat{y}_\nu$ will boost the running of $y$, and hence the LP scale of $y$ becomes lower. Nonetheless, $\hat{y}_\nu$ is important only when $y$ is very large (not far from its LP), so the LP scale would not change too much. Thus, it is still reasonable to neglect the effect of $(y_\nu)_{ab}$.
On the other hand, if $y(M_{51})< y_{\rm{c}}(M_{51})$, i.e., $F_0>0$, as $\mu$ goes up $y(\mu)$ will increase at the beginning, and then turn its direction at some scale, and then exponentially drop down to zero. If the decreasing behavior happens at a scale lower than the Planck scale, the effect of the 5-5-1 Yukawa coupling is not significant. In this case, the model is quite similar to the original quintuplet MDM model, where $\lambda_2$ blows up before the Planck scale. If the decrease happens at some scale higher than the Planck scale, then $y$ might be large enough to keep the $\lambda_2$ finite.
2 -->

3.3.Numerical calculation

The above analysis is based on analytic calculations. Below we present the results obtained by solving the RGEs numerically .
Firstly, we investigate the impact of $y(M_{51})$ on the running of couplings. Figure 2 shows the evolution of the couplings with different values of $y(M_{51})$ for fixing $\Lambda_{\rm{s}}$=10 TeV, $M_{51}$=10¹⁰ GeV, $\lambda_2(\Lambda_{\rm{s}})=0.015$, and $\lambda_3(\Lambda_{\rm{s}})=0.2$. We can see that the evolution behavior dramatically depends on the delicate input value of $y(M_{51})$.
Figure2. (color online) Evolution of the couplings ${{\lambda_2}}$, ${{\lambda_3}}$, and ${{y}}$ in the 5-5-1 model with different values of ${{y(M_{51})}}$. We have fixed ${{\Lambda_{\rm{s}}}}$=10 TeV, ${{M_{51}}}$=10¹⁰ GeV, ${{\lambda_2(\Lambda_{\rm{s}})=0.015}}$, and ${{\lambda_3(\Lambda_{\rm{s}})=0.2}}$. The red, blue, and green lines correspond to the evolution of ${{\lambda_2}}$, ${{\lambda_3}}$, and ${{y}}$, respectively.

1) If $y(M_{51})=0.662$, it will be unable to slow down the growing of $\lambda_2$, which reaches a Landau pole at a scale lower than the Planck scale.
2) If $y(M_{51})=0.666$, the Landau pole scale of $\lambda_2$ will be push up to near the Planck scale.
3) If $y(M_{51})=0.667$, all couplings will remain perturbative up to the Planck scale, and the VS conditions will be satisfied at the same time.
4) If $y(M_{51})=0.672$, although all couplings will remain perturbative up to the Planck scale, $\lambda_2$ will become negative, leading to an unstable vacuum before the Planck scale.
The above fine-tuning of $y(M_{51})$ is expected. In the $\beta$ function of $\lambda_2$ one has to arrange a delicate cancellation between $108g_2^4$ and $72y^4$ at $M_{51}$ so that $\beta_{\lambda_2}$ is well under control in a wide energy region. Otherwise, the effect of Yukawa damping on $\lambda_2$ is insufficient for a slightly smaller $y(M_{51})$ and too much for a slightly larger $y(M_{51})$ that renders a negative $\beta_{\lambda_2}$ too early.
Secondly, we study the parameter regions where the perturbativity and VS conditions are satisfied. We choose $y(M_{51})=0.662, 0.666, 0.670$ as three typical inputs, and perform scans in the $\lambda_2-\lambda_3$ plane with $\Lambda_{\rm{s}}$=10 TeV and $M_{51}$=10¹⁰ GeV. The results are shown in Fig. 3. In the blue regions, all parameters can remain perturbative up to the Planck scale, while in the orange region, the vacuum remains stable up to the Planck scale. The overlap regions simultaneously satisfy the perturbativity and VS conditions.
Figure3. (color online) Regions favored by the perturbativity and VS conditions in the ${{\lambda_2(\Lambda_{\rm{s}})-\lambda_3(\Lambda_{\rm{s}})}}$ plane for the 5-5-1 model with different values of ${{y(M_{51})}}$. We have fixed ${{\Lambda_{\rm{s}}}}$=10 TeV and ${{M_{51}}}$=10¹⁰ GeV. The blue and orange regions correspond to the parameter regions satisfying the perturbativity and VS conditions, respectively. The overlap regions are favored by both conditions.

As we vary the initial value of the Yukawa coupling $y(M_{51})$, the overlap region varies. We find that the perturbativity and VS conditions constrain $\lambda_3(\Lambda_{\rm{s}})$ within a range of $ 0.14 < \lambda_3(\Lambda_{\rm{s}}) < 0.4$, and $ \lambda_2(\Lambda_{\rm{s}})$ within a range of $0 < \lambda_2(\Lambda_{\rm{s}}) < 0.053. $ This favored range of $\lambda_3(\Lambda_{\rm{s}}) $ is also indicated as a green strip in Fig. 1, for comparing with other phenomenological constraints.
Figure1. (color online) Constraints on the real scalar quintuplet MDM model in the $\lambda_3-m_0$ plane. The purple (blue) strip corresponds to the 1$\sigma$ range of the DM relic abundance measured by the Planck experiment [55] without (with) the Sommerfeld enhancement effect. The red regions are excluded at 95% CL by the combined result of the MAGIC and Fermi-LAT indirect detection experiments [56]. The green band indicates the region satisfying the vacuum stability and perturbativity conditions in the 5-5-1 model described in Section 3.

The VS condition $\lambda_3+2\sqrt{\lambda\lambda_2}\geqslant0$ seems to allow a small negative $\lambda_3(\Lambda_{\rm{s}})$. However, the results in Fig. 3 exclude the whole regions with negative $\lambda_3(\Lambda_{\rm{s}})$. The reason is that the vacuum stability also requires the Higgs quartic coupling $\lambda$ remaining positive when running to higher scales. In the SM, $\lambda$ will turn negative at a scale ~10 ⁹ GeV. The existence of the quintuplet scalar could change this behavior, because the $\lambda_3$ coupling has a positive contribution to the $\beta$ function of $\lambda$, as shown in Eq. (13). If $|\lambda_3|$ is not large enough to turn around the trend of $\lambda$, however, the vacuum will be unstable at some scale. We set $\lambda_3(\Lambda_{\rm{s}})=-0.02$, $\lambda_2(\Lambda_{\rm{s}})=0.015$, and $y(M_{51})=0.666$, and demonstrate the evolution of $\lambda(\mu)$ in Fig. 4. We can see that $\lambda$ goes to negative values at scales ~10⁹-10¹⁴ GeV. In this case, the vacuum stability cannot be ensured along the whole way to the Planck scale.
Figure4. (color online) Evolution of ${{\lambda}}$ in the 5-5-1 model with ${{\Lambda_{\rm{s}}}}$=10 TeV, ${{M_{51}}}$=10¹⁰ GeV, ${{\lambda_2(\Lambda_{\rm{s}})=0.015}}$, ${{\lambda_3(\Lambda_{\rm{s}})=-0.02}}$, and ${{y(M_{51})=}}$0.666.

To end up this section, we would like to make a comment on the possible impacts from a different scale set for $M_{51}$ other than the benchmark scale 10¹⁰ GeV. If $M_{51}$ lies below this scale, the Yukawa term can slow down the running of $\lambda_2$ because of the lower scale, thus requiring a smaller $y(M_{51})$. But $M_{51}$ cannot be many orders of magnitude below, because $\beta_{\lambda_2}$ is about to change the sign and render $\lambda_2$ negative even far below the Planck scale. $M_{51}\gg 10^{10}$ GeV is also disfavored, because $\lambda_2(M_{51})$ would be too large to be stopped by a perturbative $y$. One can employ a parallel analysis for other choices of $M_{51}$, and the configuration of the viable parameter space will not change significantly as long as $M_{51}$ is not very far from the benchmark value.

Perturbativity puts a strong constraint on scalar MDM models, especially when the multiplet lives in a large SU(2)_L representation. The scalar self-couplings usually reach a Landau pole at an energy scale far below the Planck scale, in spite of the initial values. There are two reasons leading to such a disaster. One is that the quadratic self-interaction Lagrangians result in terms with large coefficients in the $\beta$ functions of the self-couplings. Once the self-couplings obtains a modest value, such terms will drive the self-couplings to grow exponentially and soon violate perturbativity. Another reason is that there are also significant terms contributed by other couplings in the $\beta$ functions of the self-couplings, e.g., a $g_2^4$ term with a large coefficient. These terms ensure that even when their initial values are very tiny, the self-couplings will quickly obtain modest values after a short journey of running.
In a previous work [17], we found that the perturbativity problem in the real septuplet scalar MDM model is quite stubborn: even after some fermions with Yukawa interactions are introduced to slow down the running of scalar self-couplings, the model is still unable to remain perturbative up to the Planck scale. Nonetheless, this may be achievable if the MDM scalar lives in a smaller representation. Therefore, we have studied the real scalar quintuplet MDM model in this work. The observed relic abundance implies that the scale for introducing such a quintuplet should be near 10 TeV. Our calculation suggests that the quintuplet self-coupling $\lambda_2$ will hit a Landau pole at a scale of 5.6×10¹⁴ GeV, which is consistent with Ref. [47].
In order to push up this LP scale, we have extended the model with Yukawa couplings of the scalar quintuplet to a fermionic quintuplet and three fermionic singlets, resulting in the so-called 5-5-1 model. The new singlets can also play the role of right-handed neutrinos, explaining the smallness of neutrino masses by the type-I seesaw mechanism. We have found that if such Yukawa couplings are involved after a scale of $M_{51}$~10¹⁰ GeV, all couplings can remain perturbative up to the Planck scale. The reason is that the Yukawa couplings contribute a large negative term to the $\beta$ function of $\lambda_2$ and hence slow down the growing of $\lambda_2$ at high scales.
We have also investigated the parameter regions favored by the perturbativity and vacuum stability conditions up to the Planck scale. It has been found that these conditions constrain the Higgs-quintuplet coupling $\lambda_3(\Lambda_{\rm{s}})$ at the quintuplet scale $\Lambda_{\rm{s}}$ within a range of $0.14 < \lambda_3(\Lambda_{\rm{s}}) < $0.4, and the quintuplet self-coupling $\lambda_2(\Lambda_{\rm{s}})$ within a range of $0 < \lambda_2(\Lambda_{\rm{s}}) < 0.53$.

At one-loop level, the $\beta $ functions of SM couplings are give by

$ \tag{A1} \begin{split} \beta _{{g_1}}^{{\rm{SM}}} &= \frac{1}{{16{\pi ^2}}}\frac{{41}}{{10}}g_1^3,\;\;\beta _{{g_2}}^{{\rm{SM}}} = \frac{1}{{16{\pi ^2}}}\left( { - \frac{{19}}{6}} \right)g_2^3,\;\;\\ \beta _{{g_3}}^{{\rm{SM}}}& = \frac{1}{{16{\pi ^2}}}( - 7)g_3^3, \end{split} $

$ \tag{A2} \begin{aligned} \beta_{y_t}^{\rm{SM}}=\frac{1}{16\pi^2}y_t\left(\frac{9}{2}y_t^2-\frac{9}{4}g_2^2-\frac{17}{20}g_1^2-8g_3^2\right), \end{aligned} $

$ \tag{A3} \begin{split} \beta _\lambda ^{{\rm{SM}}} =& \frac{1}{{16{\pi ^2}}}\left\{ {24{\lambda ^2} - 6y_t^4 + \frac{3}{8}\left[ {2g_2^4 + {{\left( {g_2^2 + \frac{3}{5}g_1^2} \right)}^2}} \right]} \right.\\& \left. { + \lambda \left( { - 9g_2^2 - \frac{9}{5}g_1^2 + 12y_t^2} \right)} \right\}. \end{split} $

As most of the Yukawa couplings are negligible, only the top Yukawa coupling $y_t$ is considered in the above expressions.
In the RGE calculation, we use the following $\overline{{\rm{MS}}}$ values at $m_Z$ as initial values for gauge couplings [66]:

$ \tag{A4} \alpha_{\rm{s}}(m_Z) = \frac{1}{4\pi}\left[g_3(m_Z)\right]^2=0.1184, $

$ \tag{A5} \alpha(m_Z) = \frac{1}{4\pi}\left[g_2(m_Z)s_{\rm{W}}(m_Z)\right]^2=\frac{1}{127.926}, $

$ \tag{A6} s_{\rm{W}}^2 = \sin^2\theta_{\rm{W}}(m_Z)=0.2312. $

The measured values of $y_t$ and $\lambda$ are obtained from the pole masses of the top quark and the Higgs boson, $m_t$ and $m_h$, respectively. Therefore, we need to derive their $\overline{{\rm{MS}}}$ values $y_t(m_Z)$ and $\lambda(m_Z)$ at $m_Z$ by the matching conditions [67]

$ \tag{A7} y_t(\mu_0)=\frac{\sqrt{2}m_t}{v}\left[1+\delta_t(\mu_0)\right], \;\; \lambda(\mu_0)=\frac{m_h^2}{2v^2}\left[1+\delta_h(\mu_0)\right], $

setting $\mu_0 = m_Z$. The related functions are:

$ \tag{A8} \delta_t(\mu_0) = \left(-\frac{4\alpha_s}{4\pi}-\frac{4}{3}\frac{\alpha}{4\pi}+\frac{9}{4}\frac{m_t^2}{16\pi^2v^2}\right)\ln\frac{\mu_0^2}{m_t^2}+c_t, $

$ \tag{A9} \delta_h(\mu_0) = \frac{2v^2}{m_h^2}\frac{1}{32\pi^2v^4}\left[h_0(\mu_0)+m_h^2h_1(\mu_0)+m_h^4h_2(\mu_0)\right], $

$ \tag{A10} h_0(\mu_0) = -24m_t^4\ln\frac{\mu_0^2}{m_t^2}+6m_Z^4\ln\frac{\mu_0^2}{m_Z^2}+12m_W^4\ln\frac{\mu_0^2}{m_W^2}+c_0, $

$ \tag{A11} h_1(\mu_0) = 12m_t^2\ln\frac{\mu_0^2}{m_t^2}-6m_Z^2\ln\frac{\mu_0^2}{m_Z^2}-12m_W^2\ln\frac{\mu_0^2}{m_W^2}+c_1, $

$ \tag{A12} h_2(\mu_0) = \frac{9}{2}\ln\frac{\mu_0^2}{m_h^2}+\frac{1}{2}\ln\frac{\mu_0^2}{m_Z^2}+\ln\frac{\mu_0^2}{m_W^2}+c_2. $

Here the constants $c_0$, $c_1$, and $c_2$ are independent of $\mu_0$. Their contributions to $\delta_h$ are less than 0.02 and can be neglected. The constant $c_t$ lies in a range of $-0.052\leqslant c_t\leqslant$?0.042. We take$c_t=-0.052$ in the calculation, but choosing another value within the range would not essentially change our results.

In order to account for the SE effect on the DM relic abundance and the annihilation cross section related in indirect detection, we consider the two-body Schr?dinger equations for (co-)annihilation pairs of dark sector particles, following Refs. [6, 15]. The two-body states can be categorized by their charges, denoted as $Q$. For the scalar quintuplet, there are three types of initial states ($Q=0, 1, 2$), and each of them is endowed with an annihilation rate matrix and a potential matrix. The potential matrices account for the non-relativistic effect of exchanging electroweak gauge bosons between the two incoming particles, with the following forms:

$ \tag{B1} \begin{split} &{V_{Q = 0}} = \left( {\begin{array}{*{20}{c}} {8\Delta m - 4A}&{ - 2B}&0\\ { - 2B}&{2\Delta m - A}&{ - 3\sqrt 2 B}\\ 0&{ - 3\sqrt 2 B}&0 \end{array}} \right), \\ &{V_{Q = 1}} = \left( {\begin{array}{*{20}{c}} {5\Delta m - 2A}&{ - \sqrt 6 B}\\ { - \sqrt 6 B}&{2\Delta m - 3B} \end{array}} \right), \\ &{V_{Q = 2}} = \left( {\begin{array}{*{20}{c}} {4\Delta m}&{ - 2\sqrt 3 B}\\ { - 2\sqrt 3 B}&{2\Delta m + A} \end{array}} \right). \end{split} $

Here $A=\alpha/r+\alpha_2c_{W}^2e^{-m_Zr}/r$ and $B=\alpha_2e^{-m_Wr}/r$ corresponds to neutral and charged gauge bosons, respectively.
The Schr?dinger equations for the two-body wave functions $\psi_i^{(j)}$ are given by

$ \tag{B2} -\frac{1}{m_0}\frac{\partial^2\psi_i^{(j)}}{\partial r^2}+\sum\limits_{k} V_{ik}\psi_{k}^{(j)}=K\psi_i^{(j)}, $

where $K=m_0\beta^2$ with $\beta$ a characteristic velocity of DM particles. We adopt $\beta=10^{-3}$ and 0.28 for DM particles in the Galaxy and at the freeze-out epoch, respectively. The boundary conditions are

$ \tag{B3} \psi_{i}^{(j)}(0)=\delta_{i}^j, \quad \frac{\partial\psi_i^{(j)}}{\partial r}(\infty)=i\sqrt{m_0(K-V_{ii}(\infty))}\, \psi_i^{(j)}(\infty). $

After solving these equations, annihilation and coannihilation cross sections with the SE effect are given by

$ \tag{B4} \sigma_i v=(A\Gamma A^{\dagger })_{ii} $

where $A_{ij}=\psi_i^{(j)}(\infty)$, and $\Gamma$ is the tree-level annihilation rate matrix.
The annihilation rate matrices account for annihilation and coannihilation of dark sector particles into SM particles at tree level. Dominant final states are pairs of gauge bosons and of Higgs bosons. Using the optical theorem, the annihilation rate matrices are given by

$ \tag{B5} \begin{split} &\begin{aligned} {{\Gamma _{Q = 0}}} =&{\displaystyle\frac{{6\pi \alpha _2^2}}{{25m_0^2}}}\left( {\begin{array}{*{20}{c}} {12}&6&{2\sqrt 2 }\\ 6&9&{5\sqrt 2 }\\ {2\sqrt 2 }&{5\sqrt 2 }&6 \end{array}} \right)\\ &+ {\displaystyle\frac{{\lambda _3^2}}{{200\pi m_0^2}}}\left( {\begin{array}{*{20}{c}} 1&1&{\displaystyle\frac{1}{{\sqrt 2 }}}\\ 1&1&{\displaystyle\frac{1}{{\sqrt 2 }}}\\ {\displaystyle\frac{1}{{\sqrt 2 }}}&{\displaystyle\frac{1}{{\sqrt 2 }}}&{\displaystyle\frac{1}{2}} \end{array}} \right), \end{aligned}\\ &{{\Gamma _{Q = 1}}} = {\displaystyle\frac{{6\pi \alpha _2^2}}{{25m_0^2}}}\left( {\begin{array}{*{20}{c}} 6&{\sqrt 6 }\\ {\sqrt 6 }&1 \end{array}} \right), \\ &{{\Gamma _{Q = 2}}} = {\displaystyle\frac{{6\pi \alpha _2^2}}{{25m_0^2}}}\left( {\begin{array}{*{20}{c}} 4&{ - \sqrt {12} }\\ { - \sqrt {12} }&3 \end{array}} \right). \end{split} $

These matrices are utilized in the calculation of relic abundance. One may check that the sum of all the diagonal elements can reproduce the tree-level effective annihilation cross section (11). Note that when computing annihilation cross sections, the states with $Q=1, 2$ contribute twice since there are contributions from the (++, ?) pair as well as from the (??, +) pair.
In order to estimate the constraint from the MAGIC and Fermi-LAT experiments, whose result was obtained by assuming DM totally annihilating into $W^+W^-$ [56], we also utilize the annihilation rate matrix for $\Delta^{(0)}\Delta^{(0)}\rightarrow W^+W^-$,

$ \tag{B6} \begin{split} {\Gamma _{WW}^{Q = 0}}=& {\displaystyle\frac{{2\pi \alpha _2^2}}{{m_0^2}}}\left( {\begin{array}{*{20}{c}} 4&{10}&{6\sqrt 2 }\\ {10}&{25}&{15\sqrt 2 }\\ {6\sqrt 2 }&{15\sqrt 2 }&{18} \end{array}} \right) \\ & +{\displaystyle\frac{{\lambda _3^2}}{{16\pi m_0^2}}}\left( {\begin{array}{*{20}{c}} 1&1&{\displaystyle\frac{1}{{\sqrt 2 }}}\\ 1&1&{\displaystyle\frac{1}{{\sqrt 2 }}}\\ {\displaystyle\frac{1}{{\sqrt 2 }}}&{\displaystyle\frac{1}{{\sqrt 2 }}}&{\displaystyle\frac{1}{2}} \end{array}} \right). \end{split} $