CERN-TH-2022-022

††thanks: Corresponding Author††thanks: Corresponding Author

# A Two-Component Dark Matter Model and
its Associated Gravitational Waves

Francesco Costa<EMAIL_ADDRESS>Institute
for Theoretical Physics, Georg-August University Göttingen, Friedrich-Hund-
Platz 1, Göttingen, D-37077 Germany  Sarif Khan<EMAIL_ADDRESS>Institute for Theoretical Physics, Georg-August University Göttingen,
Friedrich-Hund-Platz 1, Göttingen, D-37077 Germany  Jinsu Kim
<EMAIL_ADDRESS>Theoretical Physics Department, CERN, 1211 Geneva 23,
Switzerland

###### Abstract

We consider an extension of the Standard Model that accounts for the muon
$g-2$ tension and neutrino masses and study in detail dark matter
phenomenology. The model under consideration includes a WIMP and a FIMP scalar
dark matter candidates and thus gives rise to two-component dark matter
scenarios. We discuss different regimes and mechanisms of production,
including the novel freeze-in semi-production, and show that the WIMP and FIMP
together compose the observed relic density today. The presence of the extra
scalar fields allows phase transitions of the first order. We examine the
evolution of the vacuum state and discuss stochastic gravitational wave
signals associated with the first-order phase transition. We show that the
gravitational wave signals may be probed by future gravitational wave
experiments which may serve as a complementary detection signal.

##  I Introduction

The Standard Model (SM) of particle physics proved to be very precise in
describing the nature of the physical world. However, some of its problems
were highlighted and studied in the past decades, including the neutrino
masses, the existence of dark matter (DM), and the muon $g-2$ tension. In the
SM, neutrinos are massless. However, the evidence of neutrino oscillations
indicates otherwise Super-Kamiokande:1998kpq ; Gonzalez-Garcia:2002bkq . The
mass splitting from neutrino oscillation experiments is constrained to be
$|\Delta m_{21}^{2}|=7.42_{-0.20}^{+0.21}\times 10^{-5}\,{\rm eV}^{2}$ between
the first and the second mass eigenstates, while it is $\Delta
m_{32}^{2}=2.517_{-0.028}^{+0.026}\times 10^{-3}\,{\rm eV}^{2}$ for the second
and the third Esteban:2020cvm . Moreover, from cosmological data, we also have
a bound on the sum of their masses $\sum_{i}m_{\nu_{i}}<0.3\,{\rm eV}$
Goobar:2006xz .

The recent data coming from Fermilab Muong-2:2021ojo  increased the tension
between the SM theoretical prediction for the muon anomalous magnetic moment,
the $g-2$ factor, and the experimental data. At the moment there is a
$4.2\sigma$ discrepancy,

$\displaystyle\Delta a_{\mu}=a_{\mu}^{\rm exp}-a_{\mu}^{\rm SM}=(2.51\pm
0.59)\times 10^{-9}\,,$ (1)

suggesting the presence of new physics at a scale of hundreds of GeV.

Finally, the SM fails to accommodate one or more particles that may play the
role of the DM. Ever since the proposal by Zwicky for a dark, collision-less,
and matter-like component of the energy budget of the universe Zwicky:1933gu ;
Bertone:2016nfn , evidences from different sources for a cold, particle-like
DM have cumulated Ostriker:1973uit ; Planck:2018vyg ; Corbelli:1999af . The
most promising and studied solution to this problem is the Weekly Interacting
Massive Particle (WIMP) Gunn:1978gr ; Hut:1977zn ; Lee:1977ua ; Bertone:2004pz
. The WIMP DM is, however, strongly constrained by experimental data
XENON:2018voc ; CMS:2016lcl ; MAGIC:2016xys ; Arcadi:2017kky ; PandaX-
II:2016vec ; LUX:2016ggv . Thus, more attention has been drawn to alternative
DM production mechanisms. For example, the freeze-in mechanism has gained
increasing interest McDonald:2001vt ; Choi:2005vq ; Kusenko:2006rh ;
Hall:2009bx ; Cheung:2011nn ; Elahi:2014fsa ; Arcadi:2015ffa ; Bernal:2017kxu
; Benakli:2017whb ; Bernal:2018qlk ; Bernal:2019mhf ; Covi:2020pch ;
Khan:2020pso ; Garcia:2020hyo ; Bernal:2020qyu . The abundance of a WIMP DM is
produced through the freeze-out mechanism, which is a thermal process, and it
is generally inversely proportional to the thermal cross section. Instead, the
DM produced via the freeze-in mechanism, called Feebly Interacting Massive
Particle (FIMP), is out of equilibrium with respect to the thermal bath of the
SM particles. A small coupling between the visible sector and the DM is
predicted, making this candidate more difficult to detect and to constrain
with direct detection experiments.111 In Refs. Kim:2017mtc ; Kim:2018xsp ;
Goudelis:2018xqi , it was pointed out that such a small coupling can naturally
arise in a clockwork framework Choi:2015fiu ; Kaplan:2015fuy ; Giudice:2016yja
.  Nonetheless, both freeze-out and freeze-in production mechanisms are
physically viable and not mutually exclusive. It is thus worth exploring the
possibility of multi-component DM scenarios, where both the WIMP and FIMP DM
contribute to the current relic density $\Omega_{\rm DM}h^{2}=0.120\pm 0.001$
observed by the Planck experiment Planck:2018vyg . Recent studies on multi-
component DM scenarios include Refs. Zurek:2008qg ; Profumo:2009tb ;
Feldman:2010wy ; Bian:2013wna ; Biswas:2013nn ; Bhattacharya:2013hva ;
Bian:2014cja ; Belanger:2014vza ; Esch:2014jpa ; Arcadi:2016kmk ;
Bhattacharya:2016ysw ; DuttaBanik:2016jzv ; Bhattacharya:2017fid ;
Ahmed:2017dbb ; Bernal:2018aon ; Bhattacharya:2018cgx ; Elahi:2019jeo ;
Borah:2019aeq ; Bhattacharya:2019fgs ; Yaguna:2019cvp ; Abdallah:2019svm ;
Belanger:2020hyh ; Choi:2021yps ; DiazSaez:2021pfw ; DiazSaez:2021pmg ;
Saez:2021qta ; Belanger:2021lwd ; Ho:2021ojb ; Ho:2022erb ;
Bhattacharya:2022wtr ; Das:2022oyx .

In this paper, we consider an extension of the SM and explain the
aforementioned three problems of the SM in a single unified framework. A novel
set-up is proposed where we introduce three massive right-handed (RH)
neutrinos $N_{R}$ that, through the standard type-I seesaw mechanism
Minkowski:1977sc ; Gell-Mann:1979vob , provide a mass to the SM neutrinos. The
SM is then also extended with two SM-singlet scalar fields $\phi_{1}$ and
$\phi_{2}$ that play the role of the WIMP-like DM and FIMP-like DM,
respectively. Finally, we introduce an extra $U(1)_{L_{\mu}-L_{\tau}}$ gauge
symmetry with a related gauge boson $Z_{\mu\tau}$ which receives its mass from
a second Higgs field $\phi_{H}$. The presence of a new massive gauge boson,
with the vacuum expectation value (VEV) of $\phi_{H}$ around hundreds of GeV
can solve the $g-2$ tension Abdallah:2011ew ; Khalil:2015wua ; Lindner:2016bgg
; Chun:2016hzs ; Calibbi:2018rzv ; Arnan:2019uhr ; Calibbi:2020emz ;
Athron:2021iuf .

Appropriately assigning $U(1)_{L_{\mu}-L_{\tau}}$ charges for the DM
particles, we also address a novel production mechanism, namely the freeze-in
mechanism by semi-production processes Bringmann:2021tjr ; Hryczuk:2021qtz
like $\phi_{1}\phi_{2}\leftrightarrow\phi_{2}\phi_{2}$, that is the inverse of
the semi-annihilation process DEramo:2010keq . This mechanism produces an
exponentially increasing DM yield, and it typically requires a larger coupling
than the standard freeze-in scenarios do, since the DM abundance is also
suppressed by the small initial abundance which is generically required for
the freeze-in production mechanism.

The evolution of the vacuum state of the scalar potential becomes non-trivial
due to the three extra scalar fields. First-order phase transitions (FOPTs)
may thus arise, producing stochastic gravitational wave (GW) signals
Kamionkowski:1993fg  detectable by future GW experiments such as LISA
Baker:2019nia  which is a space-based detector comprising of three spacecraft,
utilising laser interferometry, DECIGO Seto:2001qf ; Kawamura:2006up ;
Sato:2017dkf ; Isoyama:2018rjb ; Kawamura:2020pcg  which is a proposed GW
antenna in space designed to observe GWs in the 0.1 – 10 Hz frequency range,
consisting of four clusters of LISA-like three spacecraft, and BBO
Corbin:2005ny ; Crowder:2005nr ; Harry:2006fi  which is a proposed follow-up
of the LISA experiment, aiming to form a triangular shape consisting of four
LISA-like detectors, similar to the DECIGO. For recent studies on this
subject, see, e.g., Refs. Grojean:2006bp ; Huber:2008hg ; Espinosa:2008kw ;
Caprini:2015zlo ; Artymowski:2016tme ; Baldes:2017rcu ; Beniwal:2018hyi ;
Hashino:2018zsi ; Caprini:2018mtu ; Bian:2018mkl ; Bian:2018bxr ; Bian:2019szo
; Bian:2019kmg ; Caprini:2019egz ; Di:2020ivg ; Zhou:2021cfu ;
Mohamadnejad:2021tke ; Bian:2021dmp . This possibility gives a complementary
detection signal to the standard (in-)direct detection and collider searches
that potentially can probe our model and unveil the nature of the DM. We
present a region of the model parameter space that produces detectable GW
signals from a FOPT, relieves the muon $g-2$ tension, gives masses to the SM
neutrinos, and explains the correct DM abundance by a two-component WIMP-FIMP
relic density.

The rest of the paper is organised as follows. We set up our model in Section
II, introducing the particle content, the gauge groups, and the mass spectrum
of the theory. We also present the standard type-I seesaw mechanism adopted to
explain the neutrino masses, and we give a brief explanation of the muon $g-2$
tension. In Section III, we discuss possible DM scenarios. We divide the
parameter space into three regimes and study both one-component and two-
component scenarios. In Section IV, the FOPT and its associated GWs are
studied. We showcase four benchmark points that explain the muon $g-2$,
neutrino masses, and correct DM relic density. The benchmark points predict GW
signals within the detectability of future GW experiments, in particular
Ultimate-DECIGO, which is an ultimate, idealised version of the DECIGO, whose
sensitivity is only limited by quantum noises. We conclude in Section V.

##  II Model

We consider the following Lagrangian:

$\displaystyle\mathcal{L}=\mathcal{L}_{\rm
SM}+\mathcal{L}_{\phi_{H}}+\mathcal{L}_{N}+\mathcal{L}_{\rm
DM}+\mathcal{L}_{\rm
int}-\frac{1}{4}F_{\mu\tau}^{\alpha\beta}{F_{\mu\tau}}_{\alpha\beta}\,,$ (2)

which obeys the symmetry of the complete gauge group $SU(3)_{c}\times
SU(2)_{L}\times U(1)_{Y}\times U(1)_{L_{\mu}-L_{\tau}}$, where
$\mathcal{L}_{\rm SM}$ is the SM Lagrangian including the SM Higgs field
$\phi_{h}$, $\mathcal{L}_{\phi_{H}}$ is the Lagrangian for the
$U(1)_{L_{\mu}-L_{\tau}}$ Higgs field $\phi_{H}$,

$\displaystyle\mathcal{L}_{\phi_{H}}=(D_{\mu}\phi_{H})^{\dagger}(D^{\mu}\phi_{H})+\mu_{H}^{2}|\phi_{H}|^{2}-\lambda_{H}|\phi_{H}|^{4}\,,$
(3)

and $\mathcal{L}_{N}$ is the Lagrangian for the RH neutrinos containing their
kinetic terms, mass terms, and Yukawa terms with the SM lepton doublets,

$\displaystyle\mathcal{L}_{N}$
$\displaystyle=\sum_{i=e,\mu,\tau}\frac{i}{2}\bar{N_{i}}\gamma^{\mu}D^{N}_{\mu}N_{i}-\frac{1}{2}M_{ee}\bar{N_{e}^{c}}N_{e}-\frac{1}{2}M_{\mu\tau}(\bar{N_{\mu}^{c}}N_{\tau}+\bar{N_{\tau}^{c}}N_{\mu})$
$\displaystyle\quad-
h_{e\mu}(\bar{N_{e}^{c}}N_{\mu}+\bar{N_{\mu}^{c}}N_{e})\phi_{H}^{\dagger}-h_{e\tau}(\bar{N_{e}^{c}}N_{\tau}+\bar{N_{\tau}^{c}}N_{e})\phi_{H}-\sum_{i=e,\mu,\tau}y_{i}\bar{L_{i}}\tilde{\phi}_{h}N_{i}+{\rm
h.c.}\,,$ (4)

where $\tilde{\phi}_{h}=i\sigma_{2}\phi^{*}_{h}$, and $M_{ee}$ and
$M_{\mu\tau}$ are constants whose mass-dimension is one, while $h_{e\mu}$,
$h_{e\tau}$, and $y_{i}$ are dimensionless coupling constants. In Eq. (2),
$\mathcal{L}_{\rm DM}$ is the DM Lagrangian that is given by

$\displaystyle\mathcal{L}_{\rm
DM}=\sum_{i=1,2}(D^{\mu}\phi_{i})^{\dagger}(D_{\mu}\phi_{i})-\sum_{i=1,2}\mu_{i}^{2}\phi_{i}^{\dagger}\phi_{i}-\sum_{i=1,2}\lambda_{i}(\phi_{i}^{\dagger}\phi_{i})^{2}-\lambda_{12}(\phi^{\dagger}_{1}\phi_{1})(\phi^{\dagger}_{2}\phi_{2})-\mu(\phi^{\dagger}_{1}\phi^{3}_{2}+{\rm
h.c.})\,.$ (5)

Furthermore, $\mathcal{L}_{\rm int}$ in Eq. (2) contains all the interactions
between the SM Higgs field $\phi_{h}$, the second $U(1)_{L_{\mu}-L_{\tau}}$
Higgs field $\phi_{H}$, and the DM fields $\phi_{1,2}$,

$\displaystyle\mathcal{L}_{\rm
int}=-\lambda_{hH}(\phi_{h}^{\dagger}\phi_{h})(\phi_{H}^{\dagger}\phi_{H})-\sum_{i=1,2\;j=h,H}\lambda_{ij}(\phi_{i}^{\dagger}\phi_{i})(\phi_{j}^{\dagger}\phi_{j})\,.$
(6)

The covariant derivatives in Eqs. (2) – (5) can generically be written as
$D_{\nu}X=(\partial_{\nu}+ig_{\mu\tau}Q_{\mu\tau}(X){Z_{\mu\tau}}_{\nu})X$,
where $X$ is a SM-singlet field whose $U(1)_{L_{\mu}-L_{\tau}}$ charge is
$Q_{\mu\tau}(X)$ (see Table 2), and $g_{\mu\tau}$ is the
$U(1)_{L_{\mu}-L_{\tau}}$ gauge coupling. Finally, the kinetic term for the
extra gauge boson $Z_{\mu\tau}$ is given by the last term in Eq. (2) with its
field strength tensor
$F_{\mu\tau}^{\alpha\beta}=\partial^{\alpha}Z_{\mu\tau}^{\beta}-\partial^{\beta}Z_{\mu\tau}^{\alpha}$.

In general, the Lagrangian (2) may include the gauge kinetic mixing term
Holdom:1985ag ,

$\displaystyle\mathcal{L}\supset\frac{\zeta}{2}F_{\mu\tau}^{\alpha\beta}F_{\alpha\beta}\,,$
(7)

between the $U(1)_{L_{\mu}-L_{\tau}}$ gauge boson and the SM $U(1)_{Y}$ gauge
boson whose field-strength tensor is denoted by $F_{\alpha\beta}$. In the
presence of the gauge kinetic mixing term, one may work with the physical
gauge boson states instead of the original gauge boson states by diagonalising
the mass matrix of the gauge bosons Babu:1997st . Furthermore, as we shall see
shortly, the DM phenomenology as well as the FOPT-associated GWs are
qualitatively indifferent to the gauge kinetic mixing term. Therefore, since
the kinetic mixing term does not play an important role in our discussion, we
assume, for simplicity, that $\zeta\ll 1$ in this work.222 In Refs.
Altmannshofer:2019zhy ; Biswas:2021dan , it was shown that small values of the
kinetic mixing parameter $\zeta$ are favoured from the muon $g-2$ aspect when
taking into account the experimental constraint of Borexino Harnik:2012ni ;
Borexino:2017rsf . See also, e.g., Ref. Bauer:2018onh  for a comprehensive
study on experimental constraints on the kinetic mixing parameter $\zeta$.

| Gauge
---
Group
${\rm SU(2)}_{\rm L}$
${\rm U(1)}_{\rm Y}$
| Baryon Fields
---
$Q_{L}^{i}=(u_{L}^{i},d_{L}^{i})^{T}$ | $u_{R}^{i}$ | $d_{R}^{i}$
$2$ | $1$ | $1$
$1/6$ | $2/3$ | $-1/3$
| Lepton Fields
---
$L_{L}^{i}=(\nu_{L}^{i},e_{L}^{i})^{T}$ | $e_{R}^{i}$ | $N_{R}^{i}$
$2$ | $1$ | $1$
$-1/2$ | $-1$ | $0$
| Scalar Fields
---
$\phi_{h}$ | $\phi_{H}$ | $\phi_{1}$ | $\phi_{2}$
$2$ | $1$ | $1$ | $1$
$1/2$ | $0$ | $0$ | $0$
Table 1: Particle contents and their corresponding charges under the SM gauge group. | Gauge
---
Group
$U(1)_{L_{\mu}-L_{\tau}}$
| Baryon Fields
---
$(Q^{i}_{L},u^{i}_{R},d^{i}_{R})$
$0$
| Lepton Fields
---
$(L_{L}^{e},e_{R},N_{R}^{e})$ | $(L_{L}^{\mu},\mu_{R},N_{R}^{\mu})$ | $(L_{L}^{\tau},\tau_{R},N_{R}^{\tau})$
$0$ | $1$ | $-1$
| Scalar Fields
---
$\phi_{h}$ | $\phi_{H}$ | $\phi_{1}$ | $\phi_{2}$
$0$ | $1$ | $3n_{\mu\tau}$ | $n_{\mu\tau}$
Table 2: Particle contents and their corresponding charges under
$U(1)_{L_{\mu}-L_{\tau}}$.

The presence of the interaction term between $\phi_{h}$ and $\phi_{H}$ in Eq.
(6) introduces a mass mixing. In unitary gauge, the Higgs fields $\phi_{h}$
and $\phi_{H}$ after the spontaneous breaking of $SU(2)_{L}\times
U(1)_{Y}\times U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry may be expressed as

$\displaystyle\phi_{h}=\begin{pmatrix}0\\\
\frac{v+H}{\sqrt{2}}\end{pmatrix}\qquad\text{and}\qquad\phi_{H}=\begin{pmatrix}\frac{v_{\mu\tau}+H_{\mu\tau}}{\sqrt{2}}\end{pmatrix}\,,$
(8)

where $v$ and $v_{\mu\tau}$ are the VEVs of the Higgs fields $\phi_{h}$ and
$\phi_{H}$, respectively. The scalar mass matrix is then given by

$\displaystyle\mathcal{M}^{2}_{\rm
scalar}=\left(\begin{array}[]{cc}2\lambda_{h}v^{2}&\lambda_{hH}v_{\mu\tau}v\\\
\lambda_{hH}v_{\mu\tau}v&2\lambda_{H}v_{\mu\tau}^{2}\end{array}\right)\,.$
(11)

In the presence of the Higgs-portal coupling $\lambda_{hH}$, the physical
states are obtained after diagonalising the matrix $\mathcal{M}^{2}_{\rm
scalar}$. The mass eigenstates $h_{1}$ and $h_{2}$ can be written as

$\displaystyle h_{1}=H\cos\theta+H_{\mu\tau}\sin\theta\,,\quad
h_{2}=-H\sin\theta+H_{\mu\tau}\cos\theta\,.$ (12)

The mixing angle $\theta$ and the mass eigenvalues $M^{2}_{h_{1}}$ and
$M^{2}_{h_{2}}$ are given by

$\displaystyle\tan 2\theta$
$\displaystyle=\frac{\lambda_{hH}v_{\mu\tau}v}{\lambda_{h}v^{2}-\lambda_{H}v_{\mu\tau}^{2}}\,,$
(13) $\displaystyle M^{2}_{h_{1}}$
$\displaystyle=\lambda_{h}v^{2}+\lambda_{H}v_{\mu\tau}^{2}-\sqrt{(\lambda_{h}v^{2}-\lambda_{H}v_{\mu\tau}^{2})^{2}+(\lambda_{hH}vv_{\mu\tau})^{2}}\,,$
(14) $\displaystyle M^{2}_{h_{2}}$
$\displaystyle=\lambda_{h}v^{2}+\lambda_{H}v_{\mu\tau}^{2}+\sqrt{(\lambda_{h}v^{2}-\lambda_{H}v_{\mu\tau}^{2})^{2}+(\lambda_{hH}vv_{\mu\tau})^{2}}\,.$
(15)

We identify the lighter scalar field $h_{1}$ with the observed SM Higgs field.

For the masses of the WIMP and FIMP, we obtain, with $\mu^{2}_{1,2}>0$, as

$\displaystyle
M^{2}_{1}=\mu^{2}_{1}+\lambda_{1h}\frac{v^{2}}{2}+\lambda_{1H}\frac{v^{2}_{\mu\tau}}{2}\,,\quad
M^{2}_{2}=\mu^{2}_{2}+\lambda_{2h}\frac{v^{2}}{2}+\lambda_{2H}\frac{v^{2}_{\mu\tau}}{2}\,.$
(16)

We summarise the particle contents of our model and their corresponding
charges in Table 1 and Table 2. In the remaining part of this section, we
present the standard type-I seesaw mechanism that we adopt to explain the
neutrino masses, and we briefly explain how the muon $g-2$ tension can be
relieved in our model. For a detailed explanation, readers may refer to e.g.
Refs. Biswas:2016yan ; Biswas:2016yjr .

###  II.1 Neutrino masses

Once the SM and $U(1)_{L_{\mu}-L_{\tau}}$ Higgs fields develop VEVs, the RH
neutrino mass matrix can be expressed as

$\displaystyle
M_{R}=\begin{pmatrix}M_{ee}&\frac{h_{e\mu}v_{\mu\tau}}{\sqrt{2}}&\frac{h_{e\tau}v_{\mu\tau}}{\sqrt{2}}\\\
\frac{h_{e\mu}v_{\mu\tau}}{\sqrt{2}}&0&M_{\mu\tau}e^{i\eta}\\\
\frac{h_{e\tau}v_{\mu\tau}}{\sqrt{2}}&M_{\mu\tau}e^{i\eta}&0\end{pmatrix}\,,$
(17)

where $\eta$ is the only fermionic phase factor that cannot be absorbed by
field redefinitions, and we see from the Yukawa terms in Eq. (4) that the
Dirac mass matrix can be written as

$\displaystyle M_{D}=\begin{pmatrix}\frac{y_{e}v}{\sqrt{2}}&0&0\\\
0&\frac{y_{\mu}v}{\sqrt{2}}&0\\\
0&0&\frac{y_{\tau}v}{\sqrt{2}}\end{pmatrix}\,.$ (18)

Therefore, the complete neutrino mass matrix is a $6\times 6$ matrix in the
basis $(\nu_{l},N_{l})$,

$\displaystyle M_{\nu}=\begin{pmatrix}0&M_{D}\\\
M^{T}_{D}&M_{R}\end{pmatrix}\,,$ (19)

After diagonalisation, we can obtain the mass matrix for the mass eigenstates.
Then, we can write the light neutrino mass and heavy mass matrix as follows:

$\displaystyle m^{\rm light}_{\nu}=-M^{T}_{D}M^{-1}_{R}M_{D}\,,\quad M^{\rm
heavy}_{R}=M_{R}\,.$ (20)

With the RH neutrino mass matrix elements in GeV range and the Dirac mass
matrix in keV range, one may easily obtain the neutrino mass in the correct
experimental range Esteban:2020cvm ; see also Refs. Biswas:2016yan ;
Biswas:2016yjr  for details.

The RH neutrino mass matrix squared, $(M_{R}^{\rm heavy})^{2}$, can be
diagonalised analytically when $M_{ee}=M_{\mu\tau}$ and $h_{e\tau}=h_{e\mu}$,
and we obtain the eigenvalues as

$\displaystyle(M_{ee}-h_{e\tau}v_{\mu\tau})^{2}\,,\quad
M_{ee}^{2}\,,\quad(M_{ee}+h_{e\tau}v_{\mu\tau})^{2}\,.$ (21)

In the following, we assume that this is the case.

###  II.2 Muon $g-2$

Figure 1: Constraints on $g_{\mu\tau}$ and $M_{Z_{\mu\tau}}$. The cyan region
relieves the muon $g-2$ tension. The magenta region is excluded by the
$Z\rightarrow 4\mu$ searches from LHC CMS:2012bw ; ATLAS:2014jlg ; CMS:2018yxg
, while the grey region is excluded by the neutrino trident experiments CHARM-
II and CCFR CHARM-II:1990dvf ; CCFR:1991lpl ; Altmannshofer:2014pba . See also
Ref. Chun:2018ibr  for constraints coming from the lepton universality test
and LEP searches. The colour of the scan points represents $M_{h_{2}}$. The
star ($*$) corresponds to our benchmark point 1 (see Table 3 in Section IV).

The presence of additional gauge boson $Z_{\mu\tau}$ can alleviate the
$(g-2)_{\mu}$ anomaly through the one-loop contribution, resulting in
Gninenko:2001hx ; Baek:2001kca

$\displaystyle\Delta
a_{\mu}=\frac{g^{2}_{\mu\tau}}{8\pi^{2}}\int^{1}_{0}\frac{2x(1-x)^{2}}{(1-x)^{2}+rx}dx\,,$
(22)

where $r=M^{2}_{Z_{\mu\tau}}/m^{2}_{\mu}$ . Figure 1 shows the region that
addresses the discrepancy between the experimental and theoretical values of
muon $g-2$, together with constraints from the neutrino trident experiments
such as CHARM-II CHARM-II:1990dvf  and CCFR CCFR:1991lpl ;
Altmannshofer:2014pba  and the LHC $Z\rightarrow 4\mu$ searches CMS:2012bw ;
ATLAS:2014jlg ; CMS:2018yxg . We observe that $M_{Z_{\mu\tau}}\lesssim 0.1$
GeV region with $4\times 10^{-4}\lesssim g_{\mu\tau}\lesssim 1\times 10^{-3}$
successfully explains the muon $g-2$ tension. We also present our scan points
whose colour represents the value of $M_{h_{2}}$. One may clearly see from
Fig. 1 that $M_{h_{2}}\gtrsim 1.1$ TeV is disfavoured from the muon $g-2$
point of view as long as the quartic couplings are in the perturbative regime.
Therefore, throughout the paper, we consider $M_{h_{2}}\lesssim 1.1$ TeV. The
star ($*$) in Fig. 1 depicts our benchmark point 1 (see Table 3). Strong GW
signals can be emitted from a region that explains the muon $g-2$ tension. In
Section IV, we discuss possible GW signals in detail.

##  III Two-Component Dark Matter

In this section, we examine the possibility of having DM component(s) in the
present model. As the model contains two scalar DM candidates, we may have a
single-component or two-component DM scenario depending on the mass range of
the WIMP and FIMP DM particles. We first discuss the production of DM when the
$\mu$ term in Eq. (5) is dominant and the other quartic terms associated with
the DM are also significant. We also look at the scenario when the $\mu$
parameter is less significant and quartic terms are the ones which take part
in the DM productions. In regime I and regime II, we study the effect of the
$\mu$ term on the production of DM. In these cases, depending on the mass
range of the WIMP DM, we obtain both the single-component and two-component DM
scenarios. In regime III, we study DM productions when the $\mu$ term is small
and the quartic terms are relevant. In this regime, we have a two-component DM
scenario where one component is WIMP-type DM and another component is FIMP-
type DM. In studying DM phenomenology, we have implemented our model in
FeynRules Alloul:2013bka  and generated the CalcHEP files Belyaev:2012qa . We
have then used micrOMEGAs Belanger:2018ccd  to solve the coupled Boltzmann
equations relevant for our study. Some useful analytical expressions are
derived and summarised in Appendix A. We discuss the different regimes in
detail below.

###  III.1 Regime I ($M_{h_{2}}<2M_{2}$ and $M_{1}>3M_{2}$)

As the WIMP DM mass is larger than three times the mass of the FIMP DM, a
three-body decay channel from the WIMP DM to the FIMP DM is open in this
regime. Since the WIMP DM decays into the FIMP DM, this regime gives us a
single-component DM scenario, unless the lifetime is larger than the age of
the Universe. It would require extremely small couplings to make the lifetime
larger than the age of the Universe, and we do not consider such a scenario.
Additionally, we assume that the SM and BSM Higgs masses are such that the
decay production of the FIMP DM is kinematically forbidden. Nevertheless, a
freeze-in contribution through annihilation processes,
$AB\rightarrow\phi^{\dagger}_{2}\phi_{2}$, where $A$ and $B$ are the SM
particles, will be there.

The Boltzmann equations associated with the WIMP and FIMP DM productions are
given by

$\displaystyle\frac{dY_{1}}{dx}$
$\displaystyle=-\frac{2\pi^{2}}{45}\frac{M_{\rm
Pl}M_{1}\sqrt{g_{*}(x)}}{1.66x^{2}}\langle\sigma v\rangle_{\rm
th}\left(Y^{2}_{1}-Y^{{\rm eq}2}_{1}\right)-\frac{3M_{\rm
Pl}x\sqrt{g_{*}(x)}}{1.66M^{2}_{1}g_{s}(x)}\langle\Gamma\rangle\left(Y_{1}-Y^{3}_{2}\right)\,,$
$\displaystyle\frac{dY_{2}}{dx}$
$\displaystyle=\frac{4\pi^{2}}{45}\frac{M_{\rm
Pl}M_{1}\sqrt{g_{*}}}{1.66x^{2}}\sum_{i,j\in{\rm SM},\phi_{1}}\langle\sigma
v\rangle_{ij}\left(Y^{{\rm eq}}_{i}Y^{{\rm
eq}}_{j}-Y^{2}_{2}\right)+\frac{3M_{{\rm
Pl}}x\sqrt{g_{*}(x)}}{1.66M^{2}_{1}g_{s}(x)}\langle\Gamma\rangle\left(Y_{1}-Y^{3}_{2}\right)\,,$
(23)

where $M_{\rm Pl}=1.22\times 10^{19}$ GeV is the Planck mass, and $g_{*}(x)$
and $g_{s}(x)$ are the effective and entropic degrees of freedom of the
Universe. Here, $Y_{1,2}\equiv n_{\phi_{1,2}}/S$ are the yields, with
$n_{\phi_{1,2}}$ being the number densities and $S$ the entropy density. The
first equation corresponds to the evolution of the WIMP DM and the second
equation represents the production of the FIMP DM. In the right hand side of
the first equation, the first term is the annihilation of the WIMP DM to the
SM particles. Here, $\langle\sigma v\rangle_{\rm th}$ is the thermal average
of cross section times velocity of DM annihilating to the SM particles. The
second term implies the three-body decay of the WIMP DM to the FIMP DM, where
$\langle\Gamma\rangle$ is the thermal average of the decay rate $\Gamma$,
defined as $\langle\Gamma\rangle=\Gamma K_{1}(x)/K_{2}(x)$ with $K_{1,2}$
being the modified Bessel functions of the second kind. The analytical
expression for the three-body decay is provided in Appendix A; see Eq. (63).
Similarly, the first term in the right hand side of the second equation
represents the annihilation contribution to the FIMP DM and the second term is
the decay contribution of the WIMP DM to the FIMP DM. Here, $\langle\sigma
v\rangle_{ij}$ is the thermal average of annihilations of $i,j$ particles to
FIMP DM. Due to the allowed decay term of the WIMP to the FIMP, we see that
the WIMP DM eventually decays to the FIMP DM before big bang nucleosynthesis
(BBN). We still do not have any contribution to visible energy even when the
WIMP DM decays after BBN, and thus, our model remains safe from the
constraints which come from light elements abundances Kawasaki:2017bqm .

Figure 2: Evolution of DM relic density produced by different mechanisms:
freeze-in production from three-body decay (red dashed line), thermal freeze-
out of WIMP (green double-dot-dashed line), freeze-in production from
annihilation (blue dot-dashed line), and the total contribution (black solid
line). The model parameters are chosen as follows: $M_{1}=1650$ GeV,
$M_{2}=500$ GeV, $M_{h_{2}}=500$ GeV, $M_{Z_{\mu\tau}}=0.1$ GeV,
$g_{\mu\tau}=9\times 10^{-4}$, $\lambda_{2i}=6\times 10^{-12}$ ($i=h,H$),
$\lambda_{12}=6\times 10^{-12}$, $\lambda_{1i}=0.5$ ($i=h,H$), and
$\mu=7.5\times 10^{-11}$. The magenta dot-dashed line corresponds to the
correct value of DM relic density.

In Fig. 2, we show the evolution of the DM relic density in terms of
$x=M_{1}/T$. We examine the contributions of different production mechanisms.
The blue dot-dashed line in Fig. 2 corresponds to the production of the FIMP
DM through annihilation processes which saturate at $x\simeq 1$. The
dominating processes in the annihilation contribution are the four-point
contact terms which are $AB\rightarrow\phi^{\dagger}_{2}\phi_{2}$
($A,B=h_{1,2},\phi_{1}$) and are not propagator-suppressed. The green double-
dot-dashed line represents the evolution of the WIMP DM which freezes out at
$x\simeq 20$ and starts to decay into the FIMP DM at $x\simeq 10^{5}$. Since
the WIMP DM decays into the FIMP DM, this regime corresponds to a single-
component DM scenario. The red dashed line is the freeze-in production from
the three-body decay of the WIMP DM i.e.
$\phi_{1}\rightarrow\phi_{2}\phi_{2}\phi_{2}$ that happens at $x\simeq 1$
which means that the WIMP DM is in thermal equilibrium with the cosmic soup.
Moreover, there is also another contribution that is superWIMP (SW)
contribution appears at $x\simeq 10^{5}$ Feng:2003xh . It comes from the
three-body decay of the WIMP DM. Finally, the black solid line corresponds to
the total DM relic density which comes after summing all the contributions.
The magenta dot-dashed line corresponds to the correct value of the DM relic
density. We see that, for the choice of model parameters outlined in the
caption of Fig. 2, our model correctly produces the exact amount of DM relic
density.

Figure 3: Evolutions of DM relic density for three different values of the
WIMP DM mass $M_{1}$ (left) and for three different values of the FIMP DM mass
$M_{2}$ (right). For the rest of the model parameters, see Fig. 2.

Changes in the DM relic density with respect to the masses of the WIMP DM and
FIMP DM are shown respectively in the left panel and the right panel of Fig.
3. The freeze-in production of the FIMP DM due to the annihilation is
insensitive to the WIMP DM mass. This is consistent with the observation that
there is no direct effect of the WIMP DM mass on the annihilation production
of the FIMP DM apart from the annihilation process
$\phi^{\dagger}_{1}\phi_{1}\rightarrow\phi^{\dagger}_{2}\phi_{2}$ which has
negligible dependence on the WIMP DM mass. The freeze-in production from the
decay of the WIMP DM when it is in thermal equilibrium, i.e., $x\lesssim 20$,
is in general inversely proportional to the WIMP DM mass in this regime.
However, we observe the opposite behaviour, i.e., we get more production as
the WIMP DM mass increases. This is due to the fact that, for a low value of
the WIMP DM mass, $M_{1}=1650$ GeV, we have a phase-space suppression in the
decay. Thus, we get less amount of the FIMP DM, and when we increase the WIMP
DM mass, the effect of phase space gets reduced, and we obtain more DM from
decay. The WIMP DM freezes out at $x\simeq 20$ and starts to decay into the
FIMP DM at $x\simeq 10^{5}$. We see that the WIMP DM starts to decay into the
FIMP DM earlier as the WIMP DM mass increases as the double-dot-dashed lines
in Fig. 3 indicate. This is because the decay width is linearly proportional
to the WIMP DM mass. On top of that, there is also the phase-space suppression
which further reduces the decay width and delays the WIMP decay. When the
decay of the WIMP DM happens, we see a rise in the production of the FIMP DM
at $x\simeq 10^{5}$ which is similar to the superWIMP production mechanism.

In the right panel of Fig. 3, we show the dependence on the FIMP DM mass. For
the FIMP production due to the annihilation which is represented by the dot-
dashed lines, we see a slight variation in the relic density. This is because
$h_{2}h_{2}\rightarrow\phi^{\dagger}_{2}\phi_{2}$ is the dominant process, and
we have taken the BSM Higgs mass to be $M_{h_{2}}=500$ GeV which is comparable
to the FIMP DM mass. Therefore, suppression due to the phase-space factor and
increment due to mass compensate each other. The FIMP DM production due to the
three-body decay is shown by the dashed lines. We observe one order of
magnitude difference in the DM production when we vary the FIMP DM mass from
500 GeV to 400 GeV. This happens purely because the effect of phase space is
small, and the same effect continues when we decrease the FIMP DM mass
further. Finally, let us discuss the production of the WIMP DM which decays
into the FIMP DM at $x\simeq 10^{5}$. Again, we see that, as the FIMP DM mass
decreases from 500 GeV to 400 GeV, WIMP decay width increases due to lower
phase-space suppression which indicates that the WIMP decays earlier. This is
visible by the double-dot-dashed lines. When the WIMP DM decay happens, we
have further production of the FIMP DM similar to the superWIMP production. In
both the left and right panels of Fig. 3, solid lines represent the total DM
relic density after summing all the production contributions.

Figure 4: Dependence of the DM relic density on the $\mu$ (left) and
$\lambda_{12}$ (right) parameters. For the rest of the model parameters, see
Fig. 2.

The left and right panels of Fig. 4, shows the dependence of DM relic density
for three different values of the $\mu$ and $\lambda_{12}$ parameters,
respectively. In the left panel, from the freeze-in contribution through
annihilation, we see that there is no change in the relic density coming from
the annihilation contribution which is represented by the dot-dashed line.
This is because the annihilation process associated with the $\mu$ term is
proportional to $Y^{{\rm eq}}_{1}Y_{2}$ ($Y_{2}=0$ at initial value of $x$),
while other annihilation terms are proportional to $Y^{{\rm eq}}_{A}Y^{{\rm
eq}}_{B}$ ($A,B$ are the annihilating particles). Let us turn to the
production of the FIMP DM from the three-body decay of the WIMP DM. The
production of the FIMP DM before $x\simeq 10$ occurs in the domain when the
WIMP DM is still in thermal equilibrium, and from its decay, the FIMP DM is
produced. We clearly see that the production has a quadratic dependence on the
$\mu$ parameter which is perfectly consistent with the analytical expression
given in the Appendix A; see Eq. (63). We note that the freeze-out temperature
of the WIMP DM does not depend on the $\mu$ parameter while the WIMP DM decay
does. Decay of the WIMP happens earlier (later) for a higher (lower) value of
$\mu$. This is consistent with the analytical expression; see Appendix A for
details. Depending on the decay occurrence, the superWIMP contribution to the
FIMP DM happens earlier or later and has an equal contribution for all three
values, as the freeze-out contributions do not depend on the $\mu$ parameter,
and this contribution is equal to $\Omega^{2}_{\rm DM}h^{2}=\Omega^{1}_{\rm
DM}h^{2}(M_{2}/M_{1})$.

From the right panel of Fig. 4, one may clearly see that the three-body decay
and the superWIMP production do not depend on $\lambda_{12}$. However, we see
changes in the production coming from the annihilation process,
$\phi^{\dagger}_{1}\phi_{1}\rightarrow\phi^{\dagger}_{2}\phi_{2}$. In both the
left and right panels of Fig. 4, the solid lines correspond to the total
contribution in DM relic density.

Figure 5: Dependence of the DM relic density on the quartic coupling of the
WIMP with the Higgses $\lambda_{1i}$ ($i=h,H$) (left) and the quartic coupling
of the FIMP with Higgses $\lambda_{2i}$ ($i=h,H$) (right). For the rest of the
model parameters, see Fig. 2.

Finally, the dependence of the DM relic density on the quartic couplings of
the DM with the Higgses, $\lambda_{ij}$ ($i=1,2$ and $j=h,H$), is shown in
Fig. 5. From the left panel, we see that the production of the FIMP DM from
the freeze-in by the three-body decay and annihilation does not change; see
the red dashed and dot-dashed lines. We can also see that there is a
significant change in the final value when we consider the freeze-out
production of the WIMP DM. This can be explained in a very simple way. The
WIMP relic density is determined from the inverse of thermal average of cross
section times velocity, namely $\Omega^{1}_{\rm DM}h^{2}\propto
1/\lambda^{2}_{1i}$ ($i=h,H$). Therefore, larger values of $\lambda_{1j}$
($j=h,H$) imply that we face the situation when most of the particles
annihilate away, and thus, we have less abundance for WIMP. The relative
strength of the WIMP DM relic density due to different values of quartic
couplings is given by $\Omega^{1}_{\rm DM}h^{2}|_{A}/\Omega^{1}_{\rm
DM}h^{2}|_{B}=\lambda^{2}_{1i}|_{B}/\lambda^{2}_{1i}|_{A}$ ($i=h,H$). This is
consistent with our numerical results as one may see from the blue, red, and
green lines for the freeze-out production of the WIMP DM. Since there is no
variation of the FIMP DM production from the three-body decay and
annihilation, most of the changes in the FIMP DM production comes after
$x\simeq 10^{5}$ when the WIMP DM decays into the FIMP DM. Solid lines
correspond to the total contribution after taking into account all the
production mechanisms.

On the other hand, from the right panel of Fig. 5, we see that the FIMP DM
production from the three-body decay of the WIMP when the WIMP is in thermal
equilibrium does not change as the quartic couplings vary. We can also see
that the freeze-out production of the WIMP DM shown by the double-dot-dashed
line is not affected by the change of the quartic couplings. However, the FIMP
production from annihilation changes. This is indeed consistent with the
observation that the freeze-in contribution by annihilation processes is
proportional to the quartic couplings. In this case, the production is
directly proportional to the quartic coupling which is visible by the dot-
dashed lines; see Eqs. (70) and (71) in Appendix A.

###  III.2 Regime II ($M_{h_{2}}<2M_{2}$ and $M_{1}<3M_{2}$)

When the WIMP DM mass is less than the three times the FIMP DM mass, the
three-body decay channel of the WIMP DM to the FIMP DM is kinematically
forbidden. Therefore, in this case, we have a two-component DM scenario, one
WIMP-type DM and one FIMP-type DM.

The WIMP DM freezes out at $x\simeq 20$, and we obtain relic density in the
experimentally allowed range put by Planck Planck:2018vyg  near the Higgs
resonance region. For the FIMP DM, we examine the effect of the $\mu$ term on
the FIMP DM production and choose the FIMP DM mass in such a way that the
decay channel $h_{2}\rightarrow\phi^{\dagger}_{2}\phi_{2}$ is kinematically
forbidden. Nevertheless, we have an annihilation contribution in the
production of the FIMP DM through the freeze-in mechanism. Due to the presence
of the $\mu$ term, there exists $\phi_{1}\phi_{2}\rightarrow\phi_{2}\phi_{2}$
annihilation process, and it will exhibits an exponential growth. In the
production of the FIMP, at $x\simeq 0.01$, we have a tiny amount of FIMP DM
produced from the annihilation processes of the SM particles, and at $x\simeq
1$, the exponential enhancement will take place which will be discussed in
detail below.

The governing Boltzmann equations in this regime are given by

$\displaystyle\frac{dY_{1}}{dx}$
$\displaystyle=-\frac{2\pi^{2}}{45}\frac{M_{{\rm
Pl}}M_{1}}{1.66x^{2}}\sqrt{g_{*}}\langle\sigma v\rangle_{\rm
th}\left(Y^{2}_{1}-Y^{{\rm eq}2}_{1}\right)\,,$
$\displaystyle\frac{dY_{2}}{dx}$
$\displaystyle=-\frac{2\pi^{2}}{45}\frac{M_{{\rm
Pl}}M_{1}}{1.66x^{2}}\sqrt{g_{*}}\langle\sigma v\rangle_{\rm
exp}\left(Y^{2}_{2}-Y^{{\rm
eq}}_{1}Y_{2}\right)+\frac{4\pi^{2}}{45}\frac{M_{{\rm
Pl}}M_{1}}{1.66x^{2}}\sqrt{g_{*}}\sum_{i,j\in{\rm SM},\phi_{1}}\langle\sigma
v\rangle_{ij}\left(Y^{{\rm eq}}_{i}Y^{{\rm eq}}_{j}-Y^{2}_{2}\right)\,.$ (24)

The first (second) equation represents the evolution of the WIMP (FIMP) DM.
Here, $\langle\sigma v\rangle_{\rm exp}$ is the thermal average associated
with the exponential growth computed using the prescription described in
Appendix A.1.1. We now discuss the effect of model parameters on the
production of WIMP and FIMP DM by different mechanisms.

Figure 6: Evolution of the DM relic density in regime II. In this regime, we
have a two-component DM scenario with one WIMP DM and one FIMP DM. The
evolution of the WIMP DM is shown in the green double-dot-dashed line, while
the freeze-in productions of FIMP DM through annihilation and exponential
growth are shown in the blue dot-dashed line and red dashed line,
respectively. The black solid line corresponds to the sum of the WIMP and FIMP
DM relic densities, and the magenta dot-dashed line indicates the observed
value of DM relic density. The model parameters are chosen as $M_{1}=500$ GeV,
$M_{2}=700$ GeV, $M_{h_{2}}=1000$ GeV, $\lambda_{2i}=6\times 10^{-12}$
($i=h,H$), $\lambda_{12}=6\times 10^{-12}$, $\lambda_{1h}=0.05$,
$\lambda_{1H}=0.125$, and $\mu=10^{-6}$.

In Fig. 6, the evolution of the WIMP and FIMP DM relic densities is shown. The
model parameters are chosen in such a way that the WIMP and FIMP DM relic
densities contribute equally and generate a total DM relic density in the
correct ballpark value as referred by the Planck collaboration Planck:2018vyg
. The red dashed line in Fig. 6 corresponds to the exponential growth of FIMP
DM due to the presence of the process
$\phi_{1}\phi_{2}\rightarrow\phi_{2}\phi_{2}$. This kind of process can be
solved analytically. The co-moving number density can be expressed as

$\displaystyle Y_{2}=Y^{\rm ini}_{2}e^{\int_{x_{\rm
ini}}^{\infty}\frac{2\pi^{2}}{45}\frac{M_{{\rm
Pl}}M_{1}}{1.66x^{2}}\sqrt{g_{*}}\langle\sigma v\rangle Y^{{\rm
eq}}_{1}dx}\,,$ (25)

where $Y^{\rm ini}_{2}=Y_{2}(x=x_{\rm ini})$. Thus, we see an exponential
enhancement of the FIMP DM. The blue dot-dashed line represents the FIMP DM
production through the annihilation processes
$AB\rightarrow\phi^{\dagger}_{2}\phi_{2}$ where $A$ and $B$ are the particles
in thermal equilibrium. The green double-dot-dashed line indicates the WIMP DM
production through the freeze-out mechanism which happens at $x\simeq 20$. The
total sum of the WIMP and FIMP DM contributions is depicted by the black solid
line which matches with the correct value of DM relic density $\Omega_{\rm
DM}h^{2}=0.12$.

Figure 7: Dependence of the DM relic density on the WIMP DM mass $M_{1}$
(left) and FIMP DM mass $M_{2}$ (right). For the rest of the model parameters,
see Fig. 6.

Figure 7 shows the dependence of DM relic density on the WIMP and FIMP DM
masses. In the left panel, we see that the WIMP DM mass has no observable
effect on the freeze-in production of FIMP DM through annihilation. However,
the WIMP DM mass affects the exponential growth of FIMP DM as the thermal
average of cross section times velocity is inversely proportional to mass of
the initial state particle, which is the WIMP DM in the present case. In the
case of the WIMP DM production, when we increase or decrease the WIMP DM mass
around $M_{1}=500$ GeV, we get more abundance for the WIMP DM. This is
understood from the resonance behaviour of the Higgs-mediated diagram. In Fig.
7, we have considered $M_{h_{2}}=1000$ GeV which is the resonance region for
$M_{1}=500$ GeV DM. Thus, we get a large annihilation cross section which
results in the reduction in WIMP abundance. If the WIMP DM mass deviates from
500 GeV, we get a smaller value of annihilation cross section and higher WIMP
DM abundance. The solid lines correspond to the total contribution in DM relic
density both from the WIMP and FIMP contributions.

In the right panel of Fig. 7, we may observe the effect of FIMP DM mass on the
production of WIMP and FIMP DM by different mechanisms. The FIMP DM mass has
little impact on the FIMP DM production through annihilation, while it has an
observable effect on the exponentially enhanced production of FIMP DM. For
$x\lesssim 1$, we see that the change in the FIMP DM relic density is
proportional to the FIMP DM mass. On the other hand, for $x\gtrsim 1$, we get
a similar kind of enhancement as discussed in the previous paragraph, due to
the dependence of the thermal average of the cross section on mass. As the
chosen masses are in a large range, we see no big difference in the produced
relic densities like before. There is also no effect of the FIMP DM mass on
the production of WIMP DM through the freeze-out process. The solid lines are
the total sum of WIMP and FIMP DM relic densities, and they all match the
correct value of the DM relic density given by Planck.

Figure 8: Dependence of the DM relic density on the $\mu$ parameter (left) and
the quartic coupling $\lambda_{12}$ between the WIMP DM and FIMP DM (right).
For the rest of the model parameters, see Fig. 6.

The dependences of the production of WIMP and FIMP DM on the $\mu$ parameter
and the quartic coupling between the WIMP DM and FIMP DM $\lambda_{12}$ are
respectively shown in the left and right panels of Fig. 8. The $\mu$ parameter
only affects the process $\phi_{1}\phi_{2}\rightarrow\phi_{2}\phi_{2}$, and
thus, other DM productions do not change. Looking at the freeze-in production
of FIMP DM from annihilation and WIMP DM production through the freeze-out
process, we easily see that these production mechanisms do not vary when $\mu$
changes. However, we can see a strong dependence of the exponential
enhancement on the $\mu$ parameter. If we take $\mu\lesssim 10^{-7}$, then the
exponential enhancement is absent, while for $\mu\gtrsim 5\times 10^{-6}$,
there exists a tremendous exponential enhancement in the production which
overproduces the DM. Thus, higher values of $\mu$ are disfavoured. The solid
lines are again the total sum of WIMP and FIMP contributions, and the
variation in their values are solely due to the effect of the exponential
enhancement.

On the other hand, since the quartic coupling $\lambda_{12}$ is in the feeble
regime, it does not contribute to the WIMP DM production which is clearly
visible by the double-dot-dashed line which is same for all the three values
of $\lambda_{12}$. The dashed line, which accounts for the FIMP DM production
through the exponential enhancement, is also unchanged for different values of
$\lambda_{12}$. However, the freeze-in contribution through annihilation
depends on the $\lambda_{12}$ parameter as the dot-dashed lines indicate. The
amount of DM production through annihilation depends quadratically on the
$\lambda_{12}$ parameter.

Figure 9:  Dependence of the DM relic density on the quartic coupling between
the WIMP DM and the SM Higgs (left) and the quartic coupling between the WIMP
DM and the BSM Higgs (right). For the rest of the model parameters, see Fig.
6.

The left panel and the right panel of Fig. 9 show the dependence of the DM
relic density on the quartic coupling between the WIMP DM and the SM Higgs
$\lambda_{1h}$ and the quartic coupling between the WIMP DM and the BSM Higgs
$\lambda_{1H}$, respectively. The coupling $\lambda_{1h}$ connects the WIMP DM
to the visible sector through the SM Higgs. Since this quartic coupling does
not affect the exponential growth of the FIMP DM and has a negligible effect
on the FIMP DM through annihilation, there is no change in the FIMP DM
production for different values of $\lambda_{1h}$. We can see, however,
changes in the WIMP DM production, although the difference is small. The small
dependence on $\lambda_{1h}$ is due to the fact that the WIMP DM mass is
chosen in such a way that it lies in the BSM Higgs resonance regime. Moreover,
we have kept $\lambda_{1h}$ below $0.1$. Otherwise, the WIMP DM will be ruled
out by the direct detection experiments.

The right panel of Fig. 9 indicates that the quartic coupling $\lambda_{1H}$
has no effect on the FIMP DM production as well. However, we see that a change
in $\lambda_{1H}$ results in an order of magnitude variation in the WIMP DM
relic density. This is because our parameters are chosen such that the WIMP
mass is in the resonance region for the second, BSM Higgs, $M_{1}\simeq
M_{h_{2}}/2$. Therefore, a change in $\lambda_{1H}$ which measures the
coupling strength for $\phi^{\dagger}_{1}\phi_{1}h_{2}^{2}$ has a direct
impact on the thermal DM relic density.

Figure 10:  Dependence of the DM relic density on the quartic coupling between
the FIMP DM and the SM and BSM Higgses (left) and the BSM Higgs mass (right).
For the rest of the model parameters, see Fig. 6.

The left panel of Fig. 10 shows the dependence of the DM relic density on the
quartic coupling between the FIMP DM and the SM and BSM Higgses $\lambda_{2i}$
($i=h,H$). Since $\lambda_{2h}$ and $\lambda_{2H}$ are associated with the
FIMP DM, it does not affect the WIMP DM production as one may easily see from
the figure. The quartic coupling $\lambda_{2i}$ ($i=h,H$) also does not alter
the FIMP DM production by the exponential enhancement. On the other hand, we
see that the FIMP production by annihilation gets affected due to the
variation of $\lambda_{2h}$ and $\lambda_{2H}$. This is because the associated
annihilation processes $AB\rightarrow\phi^{\dagger}_{2}\phi_{2}$, where $A$
and $B$ belong to the SM and BSM particles, directly depend on the strength of
the $\lambda_{2h}$ and $\lambda_{2H}$ couplings. The changes in the solid
line, which is the total sum of both the FIMP and WIMP contributions, are due
to the variation in FIMP DM relic density coming from the annihilation part.

The right panel of Fig. 10 shows the dependence of the DM relic density on the
BSM Higgs mass $M_{h_{2}}$. Since in this regime, the decay process
$h_{2}\rightarrow\phi^{\dagger}_{2}\phi_{2}$ is not allowed, we do not see any
observable effect on the FIMP DM production. However, we see an effect on the
production of WIMP DM. The reason is exactly the same as the one we discussed
earlier for the left panel of Fig. 7. Here as well, since $M_{1}=500$ GeV, if
$M_{h_{2}}$ deviates from $M_{h_{2}}=1000$ GeV, we are basically going away
from the resonance region. This means that DM freezes out earlier due to the
reduction in the thermal cross section, and we get higher WIMP DM relic
density. The changes in the solid line are purely due to variation in the WIMP
contribution to the DM relic density.

###  III.3 Regime III ($M_{h_{2}}>2M_{2}$ with $\mu$ negligible)

In this regime, one should take into account the FIMP DM production from the
decay of the Higgses as well. Throughout the discussion, we assume that $\mu$
is negligible and focus on two-component DM scenarios. Since $\mu$ is
negligible, we may neglect the exponential enhancement in the FIMP DM
production. We note that this scenario is different from the individual study
of WIMP Rodejohann:2015lca ; Biswas:2016ewm ; Biswas:2016yan  and FIMP
Biswas:2016yjr  as the FIMP DM can also be produced from the annihilation of
the WIMP DM through the process
$\phi^{\dagger}_{1}\phi_{1}\rightarrow\phi^{\dagger}_{2}\phi_{2}$. This
annihilation contribution can be increased or decreased with the strength of
the $\lambda_{12}$ parameter as discussed in the right panel of Fig. 8 and can
produce the FIMP DM with the correct DM relic density. Therefore, our study on
the two-component DM scenario in the regime III is new and interesting. We
provide analytical expressions for the decay and $2\rightarrow 2$ contact
annihilation processes in Appendix A; see Eqs. (68)–(71).

The Boltzmann equations associated with the WIMP and FIMP DM are given by

$\displaystyle\frac{dY_{1}}{dx}$
$\displaystyle=-\frac{2\pi^{2}}{45}\frac{M_{{\rm
Pl}}M_{1}}{1.66x^{2}}\sqrt{g_{*}}\langle\sigma v\rangle_{\rm
th}\left(Y^{2}_{1}-Y^{{\rm eq}2}_{1}\right)\,,$ (26)
$\displaystyle\frac{dY_{2}}{dx}$ $\displaystyle=\frac{2M_{{\rm
Pl}}}{1.66M^{2}_{1}}\frac{x\sqrt{g_{*}(x)}}{g_{s}(x)}\sum_{i=1,2}\langle\Gamma_{h_{i}\rightarrow\phi^{\dagger}_{1}\phi_{1}}\rangle\left(Y^{{\rm
eq}}_{h_{i}}-Y^{2}_{2}\right)+\frac{4\pi^{2}}{45}\frac{M_{{\rm
Pl}}M_{1}\sqrt{g_{*}(x)}}{1.66x^{2}}\sum_{i,j\in{\rm
SM},\phi_{1}}\langle\sigma v\rangle_{ij}\left(Y^{{\rm eq}}_{i}Y^{{\rm
eq}}_{j}-Y^{2}_{2}\right)\,.$

In the following, we solve the above Boltzmann equations and discuss the
correlation between the model parameters by performing scans with the
following range:

$\displaystyle 10^{-3}\leq\theta\leq 10^{-1}\,,\qquad
10^{-3}\leq\lambda_{1h},\lambda_{1H}\leq 10^{-1}\,,\qquad
10^{-12}\leq\lambda_{2h},\lambda_{2H},\lambda_{12}\leq 10^{-10}\,,$
$\displaystyle 10^{-10}\leq n_{\mu\tau}\leq 10^{-8}\,,\qquad 10^{-4}\leq
g_{\mu\tau}\leq 10^{-2}\,,$ (27) $\displaystyle 10^{-3}\leq
M_{Z_{\mu\tau}}[{\rm GeV}]\leq 1\,,\qquad 200\leq M_{h_{2}}[{\rm GeV}]\leq
1100\,,\qquad 1\leq M_{1,2}[{\rm GeV}]\leq 1000\,.$

When performing the scans, we demand the total DM relic density to be in the
range $0.01\leq\Omega_{\rm DM}h^{2}\leq 0.12$. We stress that, when the sum of
the WIMP and FIMP DM relic densities is smaller than $\Omega_{\rm
DM}h^{2}=0.12$, the rest of the amount can easily be obtained by suitably
adjusting the $\mu$ parameter which we neglect at the moment.

Figure 11: Allowed parameter space in the $M_{2}$ – $\lambda_{2h}$ plane
(left) and in the $M_{2}$ – $\lambda_{2H}$ plane (right) with the total DM
relic density in the range $0.01\leq\Omega_{\rm DM}h^{2}\leq 0.12$. The
parameter ranges used in the scan are given in Eq. (27).

The allowed parameter regions in the $M_{2}$ – $\lambda_{2h}$ and $M_{2}$ –
$\lambda_{2H}$ planes after imposing $0.01\leq\Omega_{\rm DM}h^{2}\leq 0.12$
are shown in Fig. 11. From the left panel, we see a sharp correlation between
the FIMP DM mass $M_{2}$ and the coupling $\lambda_{2h}$ for $M_{2}\lesssim
62$ GeV. The sharp correlation may be understood as follows. For
$M_{2}\lesssim 62$ GeV, the DM can be produced from the SM Higgs decay. The
decay mode $h_{1}\rightarrow\phi^{\dagger}_{2}\phi_{2}$ is proportional to
$\lambda_{1h}$ and the phase-space factor $\sqrt{1-4M^{2}_{2}/M_{h_{1}}^{2}}$.
Moreover, the DM relic density is proportional to the DM mass as well. Thus,
as the DM mass increases, the $\lambda_{2h}$ coupling needs to be decreased in
order to obtain the correct DM relic density. For DM mass in the range $50$ –
$60$ GeV, we have the phase-space suppression. Hence, in order to get DM in
this range, we need a larger value of $\lambda_{2h}$. For $M_{2}\gtrsim 62.5$
GeV, we do not have the decay channel of the SM Higgs into the FIMP DM, and
annihilation processes take over. We thus do not have a sharp correlation in
the large $M_{2}$ region.

From the right panel of Fig. 11, we also observe a similar kind of behaviour
between $M_{2}$ and $\lambda_{2H}$. However, for $M_{2}\gtrsim 62.5$ GeV, the
$\lambda_{2H}$ parameter cannot be arbitrarily large as the $h_{2}$ decay mode
is present. Moreover, we see that $\lambda_{2H}$ may become as large as
$10^{-10}$ whereas $\lambda_{2h}$ can go only up to $10^{-11}$. The reason for
this is that we varied $M_{h_{2}}$ up to 1.1 TeV and that the FIMP DM relic
density through decay is proportional to $\lambda^{2}_{2H}/M_{h_{2}}$.
Therefore, $\lambda_{2H}$ may become larger as $M_{h_{2}}$ takes a larger
value, which is impossible for the SM Higgs case.

Figure 12:  Allowed parameter space in the $\lambda_{2h}$ – $\lambda_{2H}$
plane (left) and in the $M_{h_{2}}$ – $\Delta a_{\mu}$ plane (right) with the
total DM relic density in the range $0.01\leq\Omega_{\rm DM}h^{2}\leq 0.12$.
The muon $g-2$ tension may be relieved in the magenta-coloured region. The
parameter ranges used in the scan are given in Eq. (27).

The left panel of Fig. 12 shows the allowed parameter region in the
$\lambda_{2h}$ – $\lambda_{2H}$ plane. Since we have considered both the WIMP
and FIMP DM in the DM relic density bound, it is hard to bound the quartic
couplings from below as there will always be a contribution from the WIMP DM.
However, we may obtain an upper bound on $\lambda_{2h}$ and $\lambda_{2H}$
above which the DM is overproduced. We find the upper limits as
$\lambda_{2h}\lesssim 10^{-11}$ and $\lambda_{2H}\lesssim 10^{-10}$ for the
choice of model parameters used in the scan (27).

The right panel of Fig. 12 shows the allowed parameter region in the
$M_{h_{2}}$ – $\Delta a_{\mu}$ plane. The magenta-coloured region corresponds
to the correct experimental range of muon $g-2$. The parameter space above the
allowed magenta-coloured band is ruled out, and the points below the band
demand additional positive contributions in $(g-2)_{\mu}$ to match the
experimental range. One may see a correlation between $\Delta a_{\mu}$ and
$M_{h_{2}}$. For a larger value of $M_{h_{2}}$, we get a lower value of
$\Delta a_{\mu}$. This is due to the fact that a higher value of $M_{h_{2}}$
indicates a higher VEV of the BSM Higgs, $v_{\mu\tau}$. Since
$v_{\mu\tau}=M_{Z_{\mu\tau}}/g_{\mu\tau}$, a higher VEV implies a lower value
of $g_{\mu\tau}$ which reduces the $\Delta a_{\mu}$ contribution.

Figure 13:  Left: Results of the scans in the $M_{1}$ – $\langle\sigma
v\rangle_{b\bar{b}}$ plane. The points above the red or blue lines are ruled
out by indirect detection experiments MAGIC:2016xys ; Reinert:2017aga . Right:
Results of the scans in the $M_{1}$ – $\sigma_{\rm SI}$ plane. The points
above the red line is ruled out from the direct detection experiment Xenon-1T
XENON:2018voc . The blue line indicates the sensitivity of the future direct
detection experiment PandaX PandaX:2018wtu . The parameter ranges used in the
scan are given in Eq. (27).

In Fig. 13, we present the indirect and direct detection bounds on the mass of
the WIMP DM. In the left panel, the AMS 02 indirect detection bound coming
from the WIMP DM annihilation to $b\bar{b}$ MAGIC:2016xys ; Reinert:2017aga
is indicated with a red line. We see a sharp rise around $M_{1}\simeq 62$ GeV
which corresponds to the SM Higgs resonance region. A part of the region is
ruled out by the indirect detection bound and the rest of the region is to be
probed shortly by different ongoing indirect detection experiments
MAGIC:2016xys . In the right panel, the spin-independent direct detection
XENON-1T bound XENON:2018voc  on the WIMP DM is shown. One may see from Fig.
13 that a part of the $M_{1}\leq 500$ GeV region is already ruled out from the
direct detection experiments. The rest of the region will be explored in the
future by different proposed experiments like Darwin DARWIN:2016hyl  and
PandaX PandaX:2018wtu .

##  IV Gravitational Waves from Phase Transition

The presence of the extra scalar fields in addition to the SM Higgs field not
only makes the phenomenology of DM much richer, but it also makes the
evolution dynamics of the vacuum state non-trivial and may lead to a FOPT in
the early universe as opposed to the SM case whose phase transition is of the
cross-over type Kajantie:1996mn . See, e.g., Ref. Hindmarsh:2020hop  for a
recent review on the FOPT. As a consequence, stochastic GW signals may be
emitted. The produced stochastic GW signals have a potential to be detected by
future GW experiments such as LISA Baker:2019nia , DECIGO Seto:2001qf , and
BBO Corbin:2005ny , and this possibility gives a complementary detection
signal to the standard (in-)direct detection and collider searches. In this
section, we examine stochastic GW signals from a FOPT in our model and compare
them with the sensitivity curves of future GW experiments. Furthermore, we
present benchmark points that explain the muon $g-2$ tension, neutrino masses,
and correct DM relic density, while producing strong GW signals that are
within the detectability of Ultimate-DECIGO.

We closely follow Ref. Caprini:2015zlo  to estimate the stochastic GW signal
from a FOPT. The three main sources of the GWs produced by a FOPT include the
collision of bubble walls, the sound wave in the plasma, and the magneto-
hydrodynamic turbulence in the plasma, and thus,

$\displaystyle\Omega_{\rm GW}h^{2}\simeq\Omega_{\rm col}h^{2}+\Omega_{\rm
sw}h^{2}+\Omega_{\rm turb}h^{2}\,,$ (28)

where

$\displaystyle\Omega_{\rm col}h^{2}=1.67\times
10^{-5}\left(\frac{H_{*}}{\beta}\right)^{2}\left(\frac{\kappa_{\phi}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{\frac{1}{3}}\left(\frac{0.11v_{w}^{3}}{0.42+v_{w}^{2}}\right)\left(\frac{3.8\left(f/f_{\rm
col}\right)^{2.8}}{1+2.8\left(f/f_{\rm col}\right)^{3.8}}\right)\,,$ (29)
$\displaystyle\Omega_{\rm sw}h^{2}=2.65\times
10^{-6}\left(\frac{H_{*}}{\beta}\right)\left(\frac{\kappa_{v}\alpha}{1+\alpha}\right)^{2}\left(\frac{100}{g_{*}}\right)^{\frac{1}{3}}v_{w}\left(f/f_{\rm
sw}\right)^{3}\left(\frac{7}{4+3\left(f/f_{\rm
sw}\right)^{2}}\right)^{\frac{7}{2}}\,,$ (30)

and

$\displaystyle\Omega_{\rm turb}h^{2}=3.35\times
10^{-4}\left(\frac{H_{*}}{\beta}\right)\left(\frac{\kappa_{\rm
turb}\alpha}{1+\alpha}\right)^{\frac{3}{2}}\left(\frac{100}{g_{*}}\right)^{\frac{1}{3}}\left(\frac{v_{w}\left(f/f_{\rm
turb}\right)^{3}}{\left[1+\left(f/f_{\rm
turb}\right)\right]^{\frac{11}{3}}\left(1+8\pi f/h_{*}\right)}\right)\,,$ (31)

with

$\displaystyle h_{*}=1.65\times 10^{-5}\,{\rm Hz}\left(\frac{T_{*}}{100{\rm
GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}\,.$ (32)

The expressions for $f_{\rm col}$, $f_{\rm sw}$, and $f_{\rm turb}$ are given
as follows:

$\displaystyle f_{\rm col}=1.65\times 10^{-5}\,{\rm
Hz}\left(\frac{0.62}{1.8-0.1v_{w}+v_{w}^{2}}\right)\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100{\rm
GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}\,,$ (33)
$\displaystyle f_{\rm sw}=1.9\times 10^{-5}\,{\rm
Hz}\left(\frac{1}{v_{w}}\right)\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100{\rm
GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}\,,$ (34)

and

$\displaystyle f_{\rm turb}=2.7\times 10^{-5}\,{\rm
Hz}\left(\frac{1}{v_{w}}\right)\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100{\rm
GeV}}\right)\left(\frac{g_{*}}{100}\right)^{\frac{1}{6}}\,.$ (35)

Here, $g_{*}$ is the number of effective degrees of freedom at $T=T_{*}$. For
the bubble wall velocity $v_{w}$, we use Steinhardt:1981ct

$\displaystyle
v_{w}=\frac{\sqrt{1/3}+\sqrt{\alpha^{2}+2\alpha/3}}{1+\alpha}\,,$ (36)

and we adopt Kamionkowski:1993fg

$\displaystyle\kappa=\frac{0.715\alpha+(4/27)\sqrt{3\alpha/2}}{1+0.715\alpha}\,,\qquad\kappa_{v}=\frac{\alpha}{0.73+0.083\sqrt{\alpha}+\alpha}\,,\qquad\kappa_{\rm
turb}=0.1\kappa_{v}\,,$ (37)

In estimating the sound-wave contribution to the GW signal, we have ignored
the possible suppression factor associated with the lifetime of the sound-wave
source.333We thank the anonymous referee for pointing out this. The
suppression factor may be estimated as Ellis:2018mja ; Ellis:2019oqb ;
Ellis:2020awk ; Guo:2020grp

$\displaystyle\mathcal{S}_{\rm sw}\simeq 1.81\times{\rm
min}\left\\{1,\frac{2(8\pi)^{1/3}}{\sqrt{3}}v_{w}\left(\frac{H_{*}}{\beta}\right)\sqrt{\frac{1+\alpha}{\kappa_{v}\alpha}}\right\\}\,,$
(38)

which corresponds to $\mathcal{O}\left(10^{-2}-10^{-1}\right)$ for the
benchmark points (BPs) presented in Table 3. It is also important to note that
such a suppression may be followed by a possible enhancement in the turbulence
contribution to the GW signal Ellis:2018mja . The precise determination
requires dedicated and sophisticated numerical simulations which go beyond the
scope of the present work.

From Eqs. (29) – (35), one may see that the key parameters that control the GW
signal are $\alpha$, $\beta/H_{*}$, and $T_{n}$, where

$\displaystyle\alpha=\frac{\rho_{\rm vac}}{\rho_{\rm
rad}^{*}}\,,\qquad\frac{\beta}{H_{*}}=T_{*}\frac{dS_{\rm
E}}{dT}\bigg{|}_{T_{*}}\,,$ (39)

with $S_{\rm E}$ being the Euclidean action of a bubble and
$\rho_{\text{vac}}$ the energy density released during the FOPT. We note that
$\rho_{\rm rad}^{*}=g_{*}\pi^{2}T_{*}^{4}/30$. Throughout the section, we take
$T_{*}$ to be the nucleation temperature $T_{n}$, i.e., $T_{*}=T_{n}$.

To understand the dynamics of the FOPT, we use the one-loop effective
potential,444 For a gauge dependence issue, readers may refer to Refs.
Nielsen:1975fs ; Fukuda:1975di ; Patel:2011th ; Chiang:2017zbz ; Croon:2020cgk
.

$\displaystyle V_{\rm eff}^{\rm 1-loop}=V^{\rm tree}+V^{\rm 1-loop}_{{\rm
eff},T=0}+V^{\rm 1-loop}_{{\rm eff},T\neq 0}\,.$ (40)

Here, $V^{\rm tree}$ is the tree-level potential and $V^{\rm 1-loop}_{{\rm
eff},T=0}$ is the zero-temperature one-loop Coleman-Weinberg contribution
Coleman:1973jx  which, in the $\overline{\rm MS}$ scheme, is given by

$\displaystyle V^{\rm 1-loop}_{{\rm
eff},T=0}=\pm\sum_{i}n_{i}\frac{M_{i}^{4}}{64\pi^{2}}\left[\ln\frac{M_{i}^{2}}{\Lambda^{2}}-c_{i}\right]\,,$
(41)

where $\Lambda$ is the renormalisation scale which we take to be
$\Lambda^{2}=(v^{2}+v_{\mu\tau}^{2})/2$, $n_{i}$ is the number of degrees of
freedom of the particle with field-dependent mass $M_{i}$, the constants
$c_{i}$ are $1/2$ ($3/2$) for transverse gauge bosons (all other particles),
and $+$ ($-$) is for bosons (fermions). The last correction, $V^{\rm
1-loop}_{{\rm eff},T\neq 0}$, is the finite-temperature one-loop correction
given by Dolan:1973qd

$\displaystyle V^{\rm 1-loop}_{{\rm eff},T\neq
0}=\sum_{i}\frac{T^{4}}{2\pi^{2}}n_{i}I_{\pm}\left(\frac{m_{i}^{2}}{T^{2}}\right)\,,$
(42)

with

$\displaystyle I_{\pm}(x)=\pm\int_{0}^{\infty}dy\,y^{2}\ln\left(1\mp
e^{-\sqrt{y^{2}+x}}\right)\,,$ (43)

where $+$ ($-$) is for fermions (bosons). To take into account the re-summed
ring diagrams, we replace the field-dependent masses as

$\displaystyle
M_{i}^{2}\rightarrow\widetilde{M}_{i}^{2}=M_{i}^{2}+\Pi_{i}(T)\,,$ (44)

where $\Pi_{i}(T)$ are the thermal masses Carrington:1991hz . For the scalars
in our model these corrections are

$\displaystyle\Pi_{h}$
$\displaystyle=\frac{T^{2}}{48}\left(3g_{1}^{2}+9g_{2}^{2}+12y_{t}^{2}+12\lambda_{h}+2\lambda_{hH}+4\lambda_{1h}+4\lambda_{2h}\right)\,,$
(45) $\displaystyle\Pi_{H}$
$\displaystyle=\frac{T^{2}}{24}\left(6g_{\mu\tau}^{2}+6\lambda_{H}+\lambda_{hH}+2\lambda_{1H}+2\lambda_{2H}+4h_{e\tau}^{2}\right)\,,$
(46) $\displaystyle\Pi_{\varphi_{1}}$
$\displaystyle=\Pi_{\eta_{1}}=\frac{T^{2}}{24}\left(8\lambda_{1}+2\lambda_{12}+\lambda_{1h}+\lambda_{1H}\right)\,,$
(47) $\displaystyle\Pi_{\varphi_{2}}$
$\displaystyle=\Pi_{\eta_{2}}=\frac{T^{2}}{24}\left(8\lambda_{2}+2\lambda_{12}+\lambda_{2h}+\lambda_{2H}\right)\,,$
(48)

where $\varphi_{1,2}$ ($\eta_{1,2}$) are the real (imaginary) components of
the DM candidate $\phi_{1,2}$, and for the gauge boson, for which only the
longitudinal mode receive corrections,

$\displaystyle\Pi_{W_{L}^{1,2,3}}=\frac{11}{6}g_{2}^{2}T^{2}\,,\quad\Pi_{B_{L}}=\frac{11}{6}g_{1}^{2}T^{2}\,,\quad\Pi_{Z_{\mu\tau
L}}=\frac{1}{3}g_{\mu\tau}^{2}T^{2}\,.$ (49)

Fermions do not receive any corrections.

In order to estimate the $\alpha$ and $\beta/H_{*}$ parameters in our model,
we performed a numerical analysis by using a modified version of
CosmoTransitions Wainwright:2011kj  together with the mass spectra given
above. We restrict our focus on the case where only the SM and the BSM Higgses
develop VEVs, taking zero VEVs for the DM candidates $\phi_{1}$ and $\phi_{2}$
throughout the temperature evolution of the system. We work with the following
input parameters:

$\displaystyle v_{\mu\tau}\,,\quad M_{h_{2}}\,,\quad\theta\,,\quad
M_{ee}\,,\quad M_{1}\,,\quad M_{2}\,,\quad\mu\,,\quad g_{\mu\tau}\,,\quad
h_{e\tau}\,,$
$\displaystyle\lambda_{1}\,,\quad\lambda_{2}\,,\quad\lambda_{1h}\,,\quad\lambda_{1H}\,,\quad\lambda_{12}\,,\quad\lambda_{2h}\,,\quad\lambda_{2H}\,,$
(50)

with the assumptions $M_{ee}=M_{\mu\tau}$ and $h_{e\tau}=h_{e\mu}$ which allow
us to analytically diagonalise the RH neutrino mass matrix as we discussed in
Sec. II. In the following, we take $h_{e\tau}=\sqrt{2}M_{ee}/v_{\mu\tau}$. We
impose the vacuum stability conditions,

$\displaystyle\lambda_{h}>0\,,\quad\lambda_{H}>0\,,\quad
4\lambda_{h}\lambda_{H}-\lambda_{hH}^{2}\geq 0\,,$ (51)

as well as the perturbativity and unitarity bounds,

$\displaystyle\left|\lambda_{h}\right|<4\pi,\quad\left|\lambda_{H}\right|<4\pi,\quad\left|\lambda_{hH}\right|<8\pi,\quad
3\lambda_{h}+2\lambda_{H}+\sqrt{\left(3\lambda_{h}-2\lambda_{H}\right)^{2}+2\lambda_{hH}^{2}}<8\pi\,.$
(52)

We focus on the following range of the parameters:

$\displaystyle 20\leq v_{\mu\tau}\,[{\rm GeV}]\leq 250\,,\quad 130\leq
M_{h_{2}}\,[{\rm GeV}]\leq 1000\,,\quad 0\leq\theta\leq 0.4\,,$ $\displaystyle
1\leq M_{ee}\,[{\rm GeV}]\leq 100\,,\quad 0\leq\lambda_{1h}\leq 0.5\,,\quad
0\leq\lambda_{1H}\leq 0.5\,,\quad 10^{-12}\leq\lambda_{2h}\leq 10^{-10}\,,$
(53) $\displaystyle 100\leq M_{1}\,[{\rm GeV}]\leq 300\,,\quad 50\leq
M_{2}\,[{\rm GeV}]\leq 250\,.$ (54)

while fixing the other parameters as follows:

$\displaystyle\lambda_{1}=\lambda_{2}=0.1\,,\quad\lambda_{2H}=\lambda_{12}=6\times
10^{-12}\,,\quad g_{\mu\tau}=6\times 10^{-4}\,,\quad\mu=10^{-7}\,.$ (55)

We note that the upper bound of the mixing angle, $\theta=0.4$, is chosen by
considering the LHC constraints on the $hVV$ couplings coming from the
measurements of the Higgs decay into gauge bosons ATLAS:2016neq . We observe
FOPTs and its associated GW signals for a wide range of the parameter values,
including the mixing angle. To show that GWs can be accompanied with both
small and large values of the mixing angle $\theta$, we present two BPs with a
large value of $\theta$ and two BPs with a small value of $\theta$.

In Fig. 14, we show the associated GW signals together with the sensitivity
curves of future GW experiments. We select four BPs and present the results in
Table 3. From Fig. 14, we see that all of our four BPs, and many other
signals, are well within the reach of detectability of Ultimate-DECIGO, while
their signal strengths are below the sensitivity curves of BBO and DECIGO.
Taking into account the suppression factor (38), we see that some of the GW
signals for the chosen BPs fall below the sensitivity curve of the Ultimate-
DECIGO, while some stay marginally within the sensitivity curve. However, we
stress that the conclusion that the GW signals associated with the FOPT within
our model, which simultaneously accounts for the muon $g-2$ tension, neutrino
masses, and two-component DM scenarios, are within the reach of the Ultimate-
DECIGO sensitivity curve remains intact.

One may see from Table 3 that, in the parameter space that can solve the muon
$g-2$ tension, generate the neutrinos masses, and produce the correct DM relic
density, GWs are also expected whose signals are strong enough to be seen by
Ultimate-DECIGO. The presented four BPs clearly showcase that the muon $g-2$,
neutrino masses, and two-component DM scenarios are accounted for in our model
which, at the same time, predicts stochastic GWs associated with the FOPT in a
single unified framework.

Figure 14: GW spectrum from the FOPT together with the sensitivity curves of future GW experiments. The predicted GW signals span in the frequency range $0.01\,{\rm Hz}\lesssim f\lesssim 100\,{\rm Hz}$ with the magnitude as large as $\Omega_{\rm GW}h^{2}\simeq 10^{-18}$. While the strengths of the signals are below the sensitivity curves of BBO and DECIGO, they are well within the reach of detectability of Ultimate-DECIGO. The data for the sensitivity curves of LISA, BBO, and DECIGO are obtained from Ref. Schmitz:2020syl , and for the data for the sensitivity curve of Ultimate-DECIGO, we used Ref. Ringwald:2020vei . Earlier work on the sensitivity curves includes e.g. Refs. Larson:1999we ; Kudoh:2005as ; Yagi:2011wg ; Thrane:2013oya ; Moore:2014lga ; Kuroyanagi:2014qza ; Saikawa:2018rcs ; Robson:2018ifk . The red dotted, green dashed, blue solid, and cyan dot-dahsed lines correspond to our four benchmark points 1, 2, 3, and 4, respectively, that are summarised in Table 3.  BP | $v_{\mu\tau}$ | $M_{h_{2}}$ | $M_{ee}$ | $M_{1}$ | $M_{2}$ | $\theta$ | $\lambda_{2h}$ | $\lambda_{1h}$ | $\lambda_{1H}$ | $\alpha$ | $\beta/H$ | $T_{n}$ | $\frac{\Omega_{\rm DM}^{1}}{\Omega_{\rm DM}}$ | $\frac{\Omega_{\rm DM}^{2}}{\Omega_{\rm DM}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 80.14 | 408.06 | 98.45 | 250.31 | 169.75 | 0.388 | $6.0\times 10^{-12}$ | 0.1 | 0.1 | 0.0036 | 4994.4 | 235.0 | 0.64 | 0.36
2 | 81.69 | 415.34 | 99.83 | 244.57 | 99.39 | 0.387 | $4.5\times 10^{-12}$ | 0.1 | 0.1 | 0.0037 | 15293.2 | 238.8 | 0.58 | 0.42
3 | 81.66 | 398.97 | 98.19 | 210.0 | 209.9 | 0.002 | $4.5\times 10^{-12}$ | 0.1 | 0.1 | 0.0068 | 4884.0 | 178.7 | 0.79 | 0.21
4 | 83.28 | 366.09 | 67.07 | 249.51 | 119.54 | 0.034 | $1.7\times 10^{-11}$ | 0.289 | 0.228 | 0.0056 | 47146.8 | 189.2 | 0.15 | 0.85
Table 3: Four BPs. Values of the mass-dimensionful parameters are given in
units of GeV. We present the model parameters, the GW-related parameters,
$\alpha$, $\beta/H$, and $T_{n}$, and the DM relic densities $\Omega_{\rm
DM}^{1,2}$. For all of these four BPs, the muon $g-2$ tension and the neutrino
masses can be accounted for. The GW signals corresponding to the four BPs are
highlighted in Fig. 14. The other parameters are chosen as
$\lambda_{2H}=6\times 10^{-12}$, $\mu=10^{-7}$, $\lambda_{12}=6\times
10^{-12}$, $g_{\mu\tau}=6\times 10^{-4}$, and $n_{\mu\tau}=10^{-8}$. The
neutrino masses can be obtained for the choice of parameters as studied in
detail in Ref. Biswas:2016yan , and $(g-2)_{\mu}$ will also be obtained in the
correct ballpark value as shown in Fig. 1. In determining the WIMP and FIMP DM
individual contribution, we have considered the total DM relic density of
$\Omega_{\rm DM}h^{2}=0.12$.

##  V Conclusion

In this paper, we studied an extension of the Standard Model that accounts for
the dark matter, the muon $g-2$ tension, and the neutrino masses, in a single
unified framework. We introduced three massive right-handed neutrinos which,
through the type-I seesaw mechanism, provide a mass to the Standard Model
neutrinos. We then extended the Standard Model by introducing two scalar
fields that play the role of the dark matter. Finally, an extra
$U(1)_{L_{\mu}-L_{\tau}}$ gauge symmetry is imposed, where the associated
gauge boson $Z_{\mu\tau}$ alleviates the muon $g-2$ tension.

As the model we considered contains two Standard Model-singlet scalar dark
matter candidates, $\phi_{1}$ and $\phi_{2}$, we examined the possibility of a
single-component as well as two-component dark matter scenarios. Focusing on
three different regimes, we showed how a single-component or two-component
dark matter scenario can emerge by numerically solving the coupled Boltzmann
equations. We found that, when the $\mu\phi_{1}^{\dagger}\phi_{2}^{3}$ term is
not negligible, both the single- and two-component scenarios may be obtained,
depending on the mass range. When the $\mu$ parameter is small or absent, we
showed that a two-component dark matter scenario naturally arises without
dependence on the mass range of the WIMP and FIMP DM. In the case of a two-
component scenario, one component becomes the WIMP-type dark matter and the
other component is the FIMP-type dark matter. We performed a numerical scan
and presented viable parameter spaces which are compatible with the current
experimental bounds such as the direct and indirect detections, relieving the
muon $g-2$ tension at the same time.

The presence of the extra scalar fields not only makes the dark matter
phenomenology richer. It also affects the evolution dynamics of the vacuum
state. As opposed to the Standard Model case whose phase transition is of the
cross-over type, a first-order phase transition may be realised in our model.
Consequently, stochastic GW signals may be emitted. We investigated the
parameter space where the first-order phase transition occurs and scrutinised
the associated stochastic gravitational wave signals. Performing a numerical
scan, we showed that the predicted gravitational waves are strong enough to be
probed by future gravitational wave experiments such as Ultimate-DECIGO.

We explicitly demonstrated that our model is capable of accommodating the
three problems of the Standard Model, namely the dark matter, neutrino masses,
and the muon $g-2$ tension, by presenting four benchmark points. The chosen
four benchmark points give rise to the first-order phase transition, and
consequently, we observe the associated gravitational wave signals. All of the
four benchmark points are within the reach of detectability of Ultimate-
DECIGO. Furthermore, the chosen benchmark points realise two-component dark
matter scenarios. We expect that the gravitational wave feature of our model
may serve as a complementary detection signal to the standard (in-)direct
detection and collider searches.

## Acknowledgements

The work of F.C. is supported by the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreement No
860881-HIDDeN. This work used the Scientific Compute Cluster at GWDG, the
joint data center of Max Planck Society for the Advancement of Science (MPG)
and University of Göttingen.

##  Appendix A Analytical Expressions for the Freeze-in Dark Matter

###  A.1 FIMP-WIMP interactions

We summarise different channels for the FIMP production, considering the
interaction between the FIMP DM and the WIMP DM.

####  A.1.1 Exponential yield

In the regime where $M_{1}<3M_{2}$, the decay of the WIMP DM is kinetically
forbidden. For the scattering process
$\phi_{2}\phi_{2}\leftrightarrow\phi_{1}\phi_{2}$, the Boltzmann equation is
given by

$\displaystyle\dot{n}_{\phi_{2}}+3Hn_{\phi_{2}}=\langle\sigma v\rangle
n_{\phi_{2}}n^{{\rm eq}}_{\phi_{1}}\,,$ (56)

to a good approximation.555The viability of the use of number densities in the
Boltzmann equations is questioned and checked in, for example, Ref. Du:2021jcj
by considering the backreaction effects and solving the Boltzmann equations at
the level of the phase-space distribution. Since $\phi_{2}$ is in a FIMP
regime, i.e., out of equilibrium, its number density is considerably low, and
we can thus neglect the quadratic term in $n_{\phi_{2}}$. Notice that if the
decay were allowed, the decay channel would become the dominant process as it
may produce the total relic density of DM with a coupling orders of magnitude
smaller than the scattering process as it is shown in Sec. III. In terms of
the yield $Y_{2}=n_{\phi_{2}}/S$, where $S$ is the entropy density, the
Boltzmann equation can be re-written as

$\displaystyle\frac{dY_{2}}{dT}=-\frac{1}{HT}\langle\sigma v\rangle n^{{\rm
eq}}_{\phi_{1}}Y_{2}\,.$ (57)

The solution has an exponential behaviour,

$\displaystyle Y_{2}=Y_{0}\exp\left[\int dT\,\frac{1}{HT}\langle\sigma
v\rangle n^{{\rm eq}}_{\phi_{1}}\right]\,.$ (58)

The thermal average of cross section times velocity, $\langle\sigma v\rangle$,
can be obtained by

$\displaystyle\langle\sigma v\rangle_{A\,B\rightarrow
C\,D}=\frac{1}{8M_{A}^{2}M_{B}^{2}K_{2}\left(M_{A}/T\right)K_{2}\left(M_{B}/T\right)}\int_{\left(M_{A}+M_{B}\right)^{2}}^{\infty}ds\,\frac{\sigma_{AB\rightarrow
CD}}{\sqrt{s}}p_{AB}K_{1}\left(\frac{\sqrt{s}}{T}\right)\,,$ (59)

where $M_{A,B}$ are the masses of $A$ and $B$, $s$ is the centre-of-mass
energy, $p_{AB}=[s-(M_{A}-M_{B})^{2}][s-(M_{A}+M_{B})^{2}]$, and $K_{1,2}$ are
the modified Bessel functions of the second kind. In our model, the cross
section is given by

$\displaystyle\sigma=\frac{|\mathcal{M}|^{2}}{32\pi
s}\sqrt{\frac{s(s-4M_{2}^{2})}{(s-(M_{1}+M_{2})^{2})(s-(M_{1}-M_{2})^{2})}}\,,$
(60)

where the matrix element is given by $|\mathcal{M}|^{2}=36\mu^{2}$.

####  A.1.2 Three-body decay

In the opposite regime, $M_{1}>3M_{2}$, the decay channel of the WIMP DM to
the FIMP DM is open. The differential decay rate is given by

$\displaystyle
d\Gamma=\frac{M_{1}\pi}{32(2\pi)^{5}}|\mathcal{M}|^{2}\,dx_{3}\,dx_{1}\,d\cos\theta\,d\phi,$
(61)

where

$\displaystyle 0\leq\phi\leq 2\pi\,,\quad-1\leq\cos\theta\leq 1\,,\quad
2\sqrt{c}\leq x_{3}\leq 1+c-a-b-2\sqrt{a}\sqrt{b}\,,\quad x_{1}^{-}\leq
x_{1}\leq x_{1}^{+}\,,$ $\displaystyle
x_{1}^{\pm}=1+a-b+c-x_{3}-\frac{1}{2}\left(2c-x_{3}\right)\left(1+\frac{a-b}{1+c-x_{3}}\right)\pm\frac{1}{2}y_{3}\sqrt{1-2\frac{a+b}{1+c-x_{3}}+\frac{(a-b)^{2}}{\left(1+c-x_{3}\right)^{2}}}\,,$
(62)

with $a=b=c=M^{2}_{2}/M^{2}_{1}$ and $y_{3}=\sqrt{x_{3}-4c}$. Therefore, we
obtain the decay rate as follows:

$\displaystyle\Gamma=\frac{M_{1}}{256\pi^{3}}|\mathcal{M}|^{2}\int^{1-3a}_{2\sqrt{a}}dx_{3}\sqrt{\frac{(x_{3}^{2}-4a)(1-3a-x_{3})}{1+a-x_{3}}}\,,$
(63)

with $|\mathcal{M}|^{2}=36\mu^{2}$.

###  A.2 FIMP-SM/BSM Higgs interactions

We consider now different FIMP production channels through the interactions
with the SM and BSM Higgses. Rotating into the mass eigenstates and
considering the electroweak broken phase, we have the following interaction:

$\displaystyle\mathcal{L}$ $\displaystyle\supset
2v|\phi_{2}|^{2}h_{2}(\lambda_{2H}\cos{\theta}-\lambda_{2h}\sin{\theta})+2v_{\mu\tau}|\phi_{2}|^{2}h_{1}(\lambda_{2H}\sin{\theta}+\lambda_{2h}\cos{\theta})$
$\displaystyle\quad+|\phi_{2}|^{2}h_{1}h_{2}(\lambda_{2H}\cos{\theta}\sin{\theta}-\lambda_{2h}\cos{\theta}\sin{\theta})+\frac{1}{2}|\phi_{2}|^{2}h_{2}^{2}(\lambda_{2H}\cos^{2}{\theta}-\lambda_{2h}\sin^{2}{\theta})$
$\displaystyle\quad+\frac{1}{2}|\phi_{2}|^{2}h_{1}^{2}(\lambda_{2H}\cos^{2}{\theta}-\lambda_{2h}\sin^{2}{\theta})$
(64)
$\displaystyle=2\lambda_{a}v|\phi_{2}|^{2}h_{2}+2\lambda_{b}v_{\mu\tau}|\phi_{2}|^{2}h_{1}+\lambda_{c}|\phi_{2}|^{2}h_{1}h_{2}+\frac{1}{2}\lambda_{d}|\phi_{2}|^{2}h_{2}^{2}+\frac{1}{2}\lambda_{e}|\phi_{2}|^{2}h_{1}^{2}\,,$
(65)

where we have introduced new coupling constants as

$\displaystyle\lambda_{a}$
$\displaystyle\equiv\lambda_{2H}\cos{\theta}-\lambda_{2h}\sin{\theta}\,,$
$\displaystyle\lambda_{b}$
$\displaystyle\equiv\lambda_{2H}\sin{\theta}+\lambda_{2h}\cos{\theta}\,,$
$\displaystyle\lambda_{c}$
$\displaystyle\equiv\lambda_{2H}\cos{\theta}\sin{\theta}-\lambda_{2h}\cos{\theta}\sin{\theta}\,,$
(66) $\displaystyle\lambda_{d}$
$\displaystyle\equiv\lambda_{2H}\cos^{2}{\theta}-\lambda_{2h}\sin^{2}{\theta}\,,$
$\displaystyle\lambda_{e}$
$\displaystyle\equiv\lambda_{2H}\cos^{2}{\theta}-\lambda_{2h}\sin^{2}{\theta}\,.$

####  A.2.1 Decay contribution

In the parameter space where $M_{h_{1},h_{2}}>2M_{2}$, the decays of the SM
and the BSM Higgs fields are allowed, with the decay rates

$\displaystyle\Gamma_{\rm SM}=\frac{v_{\mu\tau}^{2}}{4\pi
M_{h_{1}}}\lambda_{b}\sqrt{M_{h_{1}}^{2}-4M_{2}^{2}}\,,\qquad\Gamma_{\rm
BSM}=\frac{v^{2}}{4\pi
M_{h_{2}}}\lambda_{a}\sqrt{M_{h_{2}}^{2}-4M_{2}^{2}}\,.$ (67)

Solving the Boltzmann equation, we obtain

$\displaystyle Y_{2}$ $\displaystyle=\int
dT\,\frac{M_{h_{1}}^{2}}{2\pi^{2}HS}K_{1}\left(\frac{M_{h_{1}}}{T}\right)\Gamma_{\rm
SM}+\int
dT\,\frac{M_{h_{2}}^{2}}{2\pi^{2}HS}K_{1}\left(\frac{M_{h_{2}}}{T}\right)\Gamma_{\rm
BSM}$ $\displaystyle=\frac{135M_{\rm Pl}}{1.66\cdot
8\pi^{3}g_{*}^{S}\sqrt{g_{*}}}\left(\frac{\Gamma_{\rm
SM}}{M_{h_{1}}^{2}}+\frac{\Gamma_{\rm BSM}}{M_{h_{2}}^{2}}\right)\,.$ (68)

We note that, when the decay channel is kinematically open, it dominates the
production.

####  A.2.2 Scattering contribution

When the decay is inactive, the dominant contribution to the FIMP production
comes from the scattering process. The yield in this case is given by

$\displaystyle\frac{dY_{2}}{dT}=-\frac{1}{32\pi^{4}HS}\int_{4M_{h_{12}}^{2}}^{\infty}ds\,\sigma\left(s-4M_{h_{12}}^{2}\right)K_{1}\left(\frac{\sqrt{s}}{T}\right)\,,$
(69)

where $M_{h_{12}}$ is $M_{h_{1}}$ ($M_{h_{2}}$) if the relevant interaction is
governed by the coupling $\lambda_{e}$ ($\lambda_{d}$). In the $M_{h_{1}}\ll
M_{2}$ limit, we get

$\displaystyle Y_{2}=\frac{135\lambda_{e}^{2}M_{\rm Pl}}{1.66\cdot
4096\pi^{4}g_{*}^{S}\sqrt{g_{*}}M_{2}^{2}}\,,$ (70)

while in the $M_{h_{2}}\ll M_{2}$ limit, we find

$\displaystyle Y_{2}=\frac{135\lambda_{d}^{2}M_{\rm Pl}}{1.66\cdot
4096\pi^{4}g_{*}^{S}\sqrt{g_{*}}M_{2}^{2}}\,.$ (71)

As $h_{1}$ represents the SM Higgs in our consideration, the first case,
namely $M_{h_{1}}\ll M_{2}$, always holds to be the case.

## References

  * (1) Y. Fukuda et al. [Super-Kamiokande], Phys. Rev. Lett. 81, 1562-1567 (1998) [arXiv:hep-ex/9807003 [hep-ex]].
  * (2) M. C. Gonzalez-Garcia and Y. Nir, Rev. Mod. Phys. 75, 345-402 (2003) [arXiv:hep-ph/0202058 [hep-ph]].
  * (3) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, JHEP 09, 178 (2020) [arXiv:2007.14792 [hep-ph]].
  * (4) A. Goobar, S. Hannestad, E. Mortsell and H. Tu, JCAP 06, 019 (2006) [arXiv:astro-ph/0602155 [astro-ph]].
  * (5) B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126, no.14, 141801 (2021) [arXiv:2104.03281 [hep-ex]].
  * (6) F. Zwicky, Helv. Phys. Acta 6, 110-127 (1933).
  * (7) G. Bertone and D. Hooper, Rev. Mod. Phys. 90, no.4, 045002 (2018) [arXiv:1605.04909 [astro-ph.CO]].
  * (8) J. P. Ostriker and P. J. E. Peebles, Astrophys. J. 186, 467-480 (1973).
  * (9) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * (10) E. Corbelli and P. Salucci, Mon. Not. Roy. Astron. Soc. 311, 441-447 (2000) [arXiv:astro-ph/9909252 [astro-ph]].
  * (11) J. E. Gunn, B. W. Lee, I. Lerche, D. N. Schramm and G. Steigman, Astrophys. J. 223, 1015-1031 (1978).
  * (12) P. Hut, Phys. Lett. B 69, 85 (1977).
  * (13) B. W. Lee and S. Weinberg, Phys. Rev. Lett. 39, 165-168 (1977).
  * (14) G. Bertone, D. Hooper and J. Silk, Phys. Rept. 405, 279-390 (2005) [arXiv:hep-ph/0404175 [hep-ph]].
  * (15) E. Aprile et al. [XENON], Phys. Rev. Lett. 121, no.11, 111302 (2018) [arXiv:1805.12562 [astro-ph.CO]].
  * (16) V. Khachatryan et al. [CMS], JHEP 10, 129 (2016) [arXiv:1606.03577 [hep-ex]].
  * (17) M. L. Ahnen et al. [MAGIC and Fermi-LAT], JCAP 02, 039 (2016) [arXiv:1601.06590 [astro-ph.HE]].
  * (18) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo and F. S. Queiroz, Eur. Phys. J. C 78, no.3, 203 (2018) [arXiv:1703.07364 [hep-ph]].
  * (19) A. Tan et al. [PandaX-II], Phys. Rev. Lett. 117, no.12, 121303 (2016) [arXiv:1607.07400 [hep-ex]].
  * (20) D. S. Akerib et al. [LUX], Phys. Rev. Lett. 118, no.2, 021303 (2017) [arXiv:1608.07648 [astro-ph.CO]].
  * (21) J. McDonald, Phys. Rev. Lett. 88, 091304 (2002) [arXiv:hep-ph/0106249 [hep-ph]].
  * (22) K. Y. Choi and L. Roszkowski, AIP Conf. Proc. 805, no.1, 30-36 (2005) [arXiv:hep-ph/0511003 [hep-ph]].
  * (23) A. Kusenko, Phys. Rev. Lett. 97, 241301 (2006) [arXiv:hep-ph/0609081 [hep-ph]].
  * (24) L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, JHEP 03, 080 (2010) [arXiv:0911.1120 [hep-ph]].
  * (25) C. Cheung, G. Elor and L. Hall, Phys. Rev. D 84, 115021 (2011) [arXiv:1103.4394 [hep-ph]].
  * (26) F. Elahi, C. Kolda and J. Unwin, JHEP 03, 048 (2015) [arXiv:1410.6157 [hep-ph]].
  * (27) G. Arcadi, L. Covi and M. Nardecchia, Phys. Rev. D 92, no.11, 115006 (2015) [arXiv:1507.05584 [hep-ph]].
  * (28) N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen, Int. J. Mod. Phys. A 32, no.27, 1730023 (2017) [arXiv:1706.07442 [hep-ph]].
  * (29) K. Benakli, Y. Chen, E. Dudas and Y. Mambrini, Phys. Rev. D 95, no.9, 095002 (2017) [arXiv:1701.06574 [hep-ph]].
  * (30) N. Bernal, M. Dutra, Y. Mambrini, K. Olive, M. Peloso and M. Pierre, Phys. Rev. D 97, no.11, 115020 (2018) [arXiv:1803.01866 [hep-ph]].
  * (31) N. Bernal, F. Elahi, C. Maldonado and J. Unwin, JCAP 11, 026 (2019) [arXiv:1909.07992 [hep-ph]].
  * (32) L. Covi, A. Ghosh, T. Mondal and B. Mukhopadhyaya, [arXiv:2008.12550 [hep-ph]].
  * (33) S. Khan, Eur. Phys. J. C 81, no.7, 598 (2021) [arXiv:2007.13008 [hep-ph]].
  * (34) M. A. G. Garcia, Y. Mambrini, K. A. Olive and S. Verner, Phys. Rev. D 102, no.8, 083533 (2020) [arXiv:2006.03325 [hep-ph]].
  * (35) N. Bernal, J. Rubio and H. Veermäe, JCAP 10, 021 (2020) [arXiv:2006.02442 [hep-ph]].
  * (36) J. Kim and J. McDonald, Phys. Rev. D 98, no.2, 023533 (2018) [arXiv:1709.04105 [hep-ph]].
  * (37) J. Kim and J. McDonald, Phys. Rev. D 98, no.12, 123503 (2018) [arXiv:1804.02661 [hep-ph]].
  * (38) A. Goudelis, K. A. Mohan and D. Sengupta, JHEP 10, 014 (2018) [arXiv:1807.06642 [hep-ph]].
  * (39) K. Choi and S. H. Im, JHEP 01, 149 (2016) [arXiv:1511.00132 [hep-ph]].
  * (40) D. E. Kaplan and R. Rattazzi, Phys. Rev. D 93, no.8, 085007 (2016) [arXiv:1511.01827 [hep-ph]].
  * (41) G. F. Giudice and M. McCullough, JHEP 02, 036 (2017) [arXiv:1610.07962 [hep-ph]].
  * (42) K. M. Zurek, Phys. Rev. D 79, 115002 (2009) [arXiv:0811.4429 [hep-ph]].
  * (43) S. Profumo, K. Sigurdson and L. Ubaldi, JCAP 12, 016 (2009) [arXiv:0907.4374 [hep-ph]].
  * (44) D. Feldman, Z. Liu, P. Nath and G. Peim, Phys. Rev. D 81, 095017 (2010) [arXiv:1004.0649 [hep-ph]].
  * (45) L. Bian, R. Ding and B. Zhu, Phys. Lett. B 728, 105-113 (2014) [arXiv:1308.3851 [hep-ph]].
  * (46) A. Biswas, D. Majumdar, A. Sil and P. Bhattacharjee, JCAP 12, 049 (2013) [arXiv:1301.3668 [hep-ph]].
  * (47) S. Bhattacharya, A. Drozd, B. Grzadkowski and J. Wudka, JHEP 10, 158 (2013) [arXiv:1309.2986 [hep-ph]].
  * (48) L. Bian, T. Li, J. Shu and X. C. Wang, JHEP 03, 126 (2015) [arXiv:1412.5443 [hep-ph]].
  * (49) G. Bélanger, F. Boudjema, A. Pukhov and A. Semenov, Comput. Phys. Commun. 192, 322-329 (2015) [arXiv:1407.6129 [hep-ph]].
  * (50) S. Esch, M. Klasen and C. E. Yaguna, JHEP 09, 108 (2014) [arXiv:1406.0617 [hep-ph]].
  * (51) S. Bhattacharya, P. Poulose and P. Ghosh, JCAP 04, 043 (2017) [arXiv:1607.08461 [hep-ph]].
  * (52) G. Arcadi, C. Gross, O. Lebedev, Y. Mambrini, S. Pokorski and T. Toma, JHEP 12, 081 (2016) [arXiv:1611.00365 [hep-ph]].
  * (53) A. Dutta Banik, M. Pandey, D. Majumdar and A. Biswas, Eur. Phys. J. C 77, no.10, 657 (2017) [arXiv:1612.08621 [hep-ph]].
  * (54) S. Bhattacharya, P. Ghosh, T. N. Maity and T. S. Ray, JHEP 10, 088 (2017) [arXiv:1706.04699 [hep-ph]].
  * (55) A. Ahmed, M. Duch, B. Grzadkowski and M. Iglicki, Eur. Phys. J. C 78, no.11, 905 (2018) [arXiv:1710.01853 [hep-ph]].
  * (56) N. Bernal, D. Restrepo, C. Yaguna and Ó. Zapata, Phys. Rev. D 99, no.1, 015038 (2019) [arXiv:1808.03352 [hep-ph]].
  * (57) S. Bhattacharya, P. Ghosh and N. Sahu, JHEP 02, 059 (2019) [arXiv:1809.07474 [hep-ph]].
  * (58) F. Elahi and S. Khatibi, Phys. Rev. D 100, no.1, 015019 (2019) [arXiv:1902.04384 [hep-ph]].
  * (59) D. Borah, R. Roshan and A. Sil, Phys. Rev. D 100, no.5, 055027 (2019) [arXiv:1904.04837 [hep-ph]].
  * (60) S. Bhattacharya, P. Ghosh, A. K. Saha and A. Sil, JHEP 03, 090 (2020) [arXiv:1905.12583 [hep-ph]].
  * (61) C. E. Yaguna and Ó. Zapata, JHEP 03, 109 (2020) [arXiv:1911.05515 [hep-ph]].
  * (62) W. Abdallah, S. Choubey and S. Khan, JHEP 06, 095 (2019) [arXiv:1904.10015 [hep-ph]].
  * (63) G. Bélanger, A. Pukhov, C. E. Yaguna and Ó. Zapata, JHEP 09, 030 (2020) [arXiv:2006.14922 [hep-ph]].
  * (64) S. M. Choi, J. Kim, P. Ko and J. Li, JHEP 09, 028 (2021) [arXiv:2103.05956 [hep-ph]].
  * (65) B. Díaz Sáez, K. Möhling and D. Stöckinger, JCAP 10, 027 (2021) [arXiv:2103.17064 [hep-ph]].
  * (66) B. Díaz Sáez, P. Escalona, S. Norero and A. R. Zerwekh, JHEP 10, 233 (2021) [arXiv:2105.04255 [hep-ph]].
  * (67) B. D. Sáez and K. Ghorbani, Phys. Lett. B 823, 136750 (2021) [arXiv:2107.08945 [hep-ph]].
  * (68) G. Belanger, A. Mjallal and A. Pukhov, Phys. Rev. D 105, no.3, 035018 (2022) [arXiv:2108.08061 [hep-ph]].
  * (69) S. Y. Ho, P. Ko and C. T. Lu, [arXiv:2107.04375 [hep-ph]].
  * (70) S. Y. Ho, P. Ko and C. T. Lu, [arXiv:2201.06856 [hep-ph]].
  * (71) S. Bhattacharya, P. Ghosh, J. Lahiri and B. Mukhopadhyaya, [arXiv:2202.12097 [hep-ph]].
  * (72) A. Das, S. Gola, S. Mandal and N. Sinha, [arXiv:2202.01443 [hep-ph]].
  * (73) P. Minkowski, Phys. Lett. B 67, 421-428 (1977).
  * (74) M. Gell-Mann, P. Ramond and R. Slansky, Conf. Proc. C 790927, 315-321 (1979) [arXiv:1306.4669 [hep-th]].
  * (75) W. Abdallah, A. Awad, S. Khalil and H. Okada, Eur. Phys. J. C 72, 2108 (2012) [arXiv:1105.1047 [hep-ph]].
  * (76) S. Khalil and C. S. Un, Phys. Lett. B 763, 164-168 (2016) [arXiv:1509.05391 [hep-ph]].
  * (77) M. Lindner, M. Platscher and F. S. Queiroz, Phys. Rept. 731, 1-82 (2018) [arXiv:1610.06587 [hep-ph]].
  * (78) E. J. Chun and J. Kim, JHEP 07, 110 (2016) [arXiv:1605.06298 [hep-ph]].
  * (79) L. Calibbi, R. Ziegler and J. Zupan, JHEP 07, 046 (2018) [arXiv:1804.00009 [hep-ph]].
  * (80) P. Arnan, A. Crivellin, M. Fedele and F. Mescia, JHEP 06, 118 (2019) [arXiv:1904.05890 [hep-ph]].
  * (81) L. Calibbi, M. L. López-Ibáñez, A. Melis and O. Vives, JHEP 06, 087 (2020) [arXiv:2003.06633 [hep-ph]].
  * (82) P. Athron, C. Balázs, D. H. J. Jacob, W. Kotlarski, D. Stöckinger and H. Stöckinger-Kim, JHEP 09, 080 (2021) [arXiv:2104.03691 [hep-ph]].
  * (83) T. Bringmann, P. F. Depta, M. Hufnagel, J. T. Ruderman and K. Schmidt-Hoberg, Phys. Rev. Lett. 127, no.19, 19 (2021) [arXiv:2103.16572 [hep-ph]].
  * (84) A. Hryczuk and M. Laletin, JHEP 06, 026 (2021) [arXiv:2104.05684 [hep-ph]].
  * (85) F. D’Eramo and J. Thaler, JHEP 06, 109 (2010) [arXiv:1003.5912 [hep-ph]].
  * (86) M. Kamionkowski, A. Kosowsky and M. S. Turner, Phys. Rev. D 49, 2837-2851 (1994) [arXiv:astro-ph/9310044 [astro-ph]].
  * (87) J. Baker, J. Bellovary, P. L. Bender, E. Berti, R. Caldwell, J. Camp, J. W. Conklin, N. Cornish, C. Cutler and R. DeRosa, et al. [arXiv:1907.06482 [astro-ph.IM]].
  * (88) N. Seto, S. Kawamura and T. Nakamura, Phys. Rev. Lett. 87, 221103 (2001) [arXiv:astro-ph/0108011 [astro-ph]].
  * (89) S. Kawamura, T. Nakamura, M. Ando, N. Seto, K. Tsubono, K. Numata, R. Takahashi, S. Nagano, T. Ishikawa and M. Musha, et al. Class. Quant. Grav. 23, S125-S132 (2006).
  * (90) S. Sato, S. Kawamura, M. Ando, T. Nakamura, K. Tsubono, A. Araya, I. Funaki, K. Ioka, N. Kanda and S. Moriwaki, et al. J. Phys. Conf. Ser. 840, no.1, 012010 (2017).
  * (91) S. Isoyama, H. Nakano and T. Nakamura, PTEP 2018, no.7, 073E01 (2018) [arXiv:1802.06977 [gr-qc]].
  * (92) S. Kawamura, M. Ando, N. Seto, S. Sato, M. Musha, I. Kawano, J. Yokoyama, T. Tanaka, K. Ioka and T. Akutsu, et al. PTEP 2021, no.5, 05A105 (2021) [arXiv:2006.13545 [gr-qc]].
  * (93) V. Corbin and N. J. Cornish, Class. Quant. Grav. 23, 2435-2446 (2006) [arXiv:gr-qc/0512039 [gr-qc]].
  * (94) J. Crowder and N. J. Cornish, Phys. Rev. D 72, 083005 (2005) [arXiv:gr-qc/0506015 [gr-qc]].
  * (95) G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner and E. S. Phinney, Class. Quant. Grav. 23, 4887-4894 (2006) [erratum: Class. Quant. Grav. 23, 7361 (2006)].
  * (96) C. Grojean and G. Servant, Phys. Rev. D 75, 043507 (2007) [arXiv:hep-ph/0607107 [hep-ph]].
  * (97) S. J. Huber and T. Konstandin, JCAP 09, 022 (2008) [arXiv:0806.1828 [hep-ph]].
  * (98) J. R. Espinosa, T. Konstandin, J. M. No and M. Quiros, Phys. Rev. D 78, 123528 (2008) [arXiv:0809.3215 [hep-ph]].
  * (99) C. Caprini, M. Hindmarsh, S. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No, A. Petiteau, P. Schwaller and G. Servant, et al. JCAP 04, 001 (2016) [arXiv:1512.06239 [astro-ph.CO]].
  * (100) M. Artymowski, M. Lewicki and J. D. Wells, JHEP 03, 066 (2017) [arXiv:1609.07143 [hep-ph]].
  * (101) I. Baldes, JCAP 05, 028 (2017) [arXiv:1702.02117 [hep-ph]].
  * (102) A. Beniwal, M. Lewicki, M. White and A. G. Williams, JHEP 02, 183 (2019) [arXiv:1810.02380 [hep-ph]].
  * (103) K. Hashino, M. Kakizaki, S. Kanemura, P. Ko and T. Matsui, JHEP 06, 088 (2018) [arXiv:1802.02947 [hep-ph]].
  * (104) C. Caprini and D. G. Figueroa, Class. Quant. Grav. 35, no.16, 163001 (2018) [arXiv:1801.04268 [astro-ph.CO]].
  * (105) L. Bian and Y. L. Tang, JHEP 12, 006 (2018) [arXiv:1810.03172 [hep-ph]].
  * (106) L. Bian and X. Liu, Phys. Rev. D 99, no.5, 055003 (2019) [arXiv:1811.03279 [hep-ph]].
  * (107) L. Bian, W. Cheng, H. K. Guo and Y. Zhang, Chin. Phys. C 45, no.11, 113104 (2021) [arXiv:1907.13589 [hep-ph]].
  * (108) L. Bian, Y. Wu and K. P. Xie, JHEP 12, 028 (2019) [arXiv:1909.02014 [hep-ph]].
  * (109) C. Caprini, M. Chala, G. C. Dorsch, M. Hindmarsh, S. J. Huber, T. Konstandin, J. Kozaczuk, G. Nardini, J. M. No and K. Rummukainen, et al. JCAP 03, 024 (2020) [arXiv:1910.13125 [astro-ph.CO]].
  * (110) Y. Di, J. Wang, R. Zhou, L. Bian, R. G. Cai and J. Liu, Phys. Rev. Lett. 126, no.25, 251102 (2021) [arXiv:2012.15625 [astro-ph.CO]].
  * (111) R. Zhou, L. Bian and J. Shu, [arXiv:2104.03519 [hep-ph]].
  * (112) A. Mohamadnejad, JHEP 03, 188 (2022) [arXiv:2111.04342 [hep-ph]].
  * (113) L. Bian, Y. L. Tang and R. Zhou, [arXiv:2111.10608 [hep-ph]].
  * (114) B. Holdom, Phys. Lett. B 166, 196-198 (1986).
  * (115) K. S. Babu, C. F. Kolda and J. March-Russell, Phys. Rev. D 57, 6788-6792 (1998) [arXiv:hep-ph/9710441 [hep-ph]].
  * (116) W. Altmannshofer, S. Gori, J. Martín-Albo, A. Sousa and M. Wallbank, Phys. Rev. D 100, no.11, 115029 (2019) [arXiv:1902.06765 [hep-ph]].
  * (117) A. Biswas and S. Khan, [arXiv:2112.08393 [hep-ph]].
  * (118) R. Harnik, J. Kopp and P. A. N. Machado, JCAP 07, 026 (2012) [arXiv:1202.6073 [hep-ph]].
  * (119) M. Agostini et al. [Borexino], Phys. Rev. D 100, no.8, 082004 (2019) [arXiv:1707.09279 [hep-ex]].
  * (120) M. Bauer, P. Foldenauer and J. Jaeckel, JHEP 07, 094 (2018) [arXiv:1803.05466 [hep-ph]].
  * (121) A. Biswas, S. Choubey and S. Khan, JHEP 09, 147 (2016) [arXiv:1608.04194 [hep-ph]].
  * (122) A. Biswas, S. Choubey and S. Khan, JHEP 02, 123 (2017) [arXiv:1612.03067 [hep-ph]].
  * (123) S. Chatrchyan et al. [CMS], JHEP 12, 034 (2012) [arXiv:1210.3844 [hep-ex]].
  * (124) G. Aad et al. [ATLAS], Phys. Rev. Lett. 112, no.23, 231806 (2014) [arXiv:1403.5657 [hep-ex]].
  * (125) A. M. Sirunyan et al. [CMS], Phys. Lett. B 792, 345-368 (2019) [arXiv:1808.03684 [hep-ex]].
  * (126) D. Geiregat et al. [CHARM-II], Phys. Lett. B 245, 271-275 (1990).
  * (127) S. R. Mishra et al. [CCFR], Phys. Rev. Lett. 66, 3117-3120 (1991).
  * (128) W. Altmannshofer, S. Gori, M. Pospelov and I. Yavin, Phys. Rev. Lett. 113, 091801 (2014) [arXiv:1406.2332 [hep-ph]].
  * (129) E. J. Chun, A. Das, J. Kim and J. Kim, JHEP 02, 093 (2019) [erratum: JHEP 07, 024 (2019)] [arXiv:1811.04320 [hep-ph]].
  * (130) S. N. Gninenko and N. V. Krasnikov, Phys. Lett. B 513, 119 (2001) [arXiv:hep-ph/0102222 [hep-ph]].
  * (131) S. Baek, N. G. Deshpande, X. G. He and P. Ko, Phys. Rev. D 64, 055006 (2001) [arXiv:hep-ph/0104141 [hep-ph]].
  * (132) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, Comput. Phys. Commun. 185, 2250-2300 (2014) [arXiv:1310.1921 [hep-ph]].
  * (133) A. Belyaev, N. D. Christensen and A. Pukhov, Comput. Phys. Commun. 184, 1729-1769 (2013) [arXiv:1207.6082 [hep-ph]].
  * (134) G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov and B. Zaldivar, Comput. Phys. Commun. 231, 173-186 (2018) [arXiv:1801.03509 [hep-ph]].
  * (135) M. Kawasaki, K. Kohri, T. Moroi and Y. Takaesu, Phys. Rev. D 97, no.2, 023502 (2018) [arXiv:1709.01211 [hep-ph]].
  * (136) J. L. Feng, A. Rajaraman and F. Takayama, Phys. Rev. Lett. 91, 011302 (2003) [arXiv:hep-ph/0302215 [hep-ph]].
  * (137) W. Rodejohann and C. E. Yaguna, JCAP 12, 032 (2015) [arXiv:1509.04036 [hep-ph]].
  * (138) A. Biswas, S. Choubey and S. Khan, JHEP 08, 114 (2016) [arXiv:1604.06566 [hep-ph]].
  * (139) A. Reinert and M. W. Winkler, JCAP 01, 055 (2018) [arXiv:1712.00002 [astro-ph.HE]].
  * (140) J. Aalbers et al. [DARWIN], JCAP 11, 017 (2016) [arXiv:1606.07001 [astro-ph.IM]].
  * (141) H. Zhang et al. [PandaX], Sci. China Phys. Mech. Astron. 62, no.3, 31011 (2019) [arXiv:1806.02229 [physics.ins-det]].
  * (142) K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, Phys. Rev. Lett. 77, 2887-2890 (1996) [arXiv:hep-ph/9605288 [hep-ph]].
  * (143) M. B. Hindmarsh, M. Lüben, J. Lumma and M. Pauly, SciPost Phys. Lect. Notes 24, 1 (2021) [arXiv:2008.09136 [astro-ph.CO]].
  * (144) P. J. Steinhardt, Phys. Rev. D 25, 2074 (1982).
  * (145) J. Ellis, M. Lewicki and J. M. No, JCAP 04, 003 (2019) [arXiv:1809.08242 [hep-ph]].
  * (146) J. Ellis, M. Lewicki, J. M. No and V. Vaskonen, JCAP 06, 024 (2019) [arXiv:1903.09642 [hep-ph]].
  * (147) J. Ellis, M. Lewicki and J. M. No, JCAP 07, 050 (2020) [arXiv:2003.07360 [hep-ph]].
  * (148) H. K. Guo, K. Sinha, D. Vagie and G. White, JCAP 01, 001 (2021) [arXiv:2007.08537 [hep-ph]].
  * (149) N. K. Nielsen, Nucl. Phys. B 101, 173-188 (1975).
  * (150) R. Fukuda and T. Kugo, Phys. Rev. D 13, 3469 (1976).
  * (151) H. H. Patel and M. J. Ramsey-Musolf, JHEP 07, 029 (2011) [arXiv:1101.4665 [hep-ph]].
  * (152) C. W. Chiang and E. Senaha, Phys. Lett. B 774, 489-493 (2017) [arXiv:1707.06765 [hep-ph]].
  * (153) D. Croon, O. Gould, P. Schicho, T. V. I. Tenkanen and G. White, JHEP 04, 055 (2021) [arXiv:2009.10080 [hep-ph]].
  * (154) S. R. Coleman and E. J. Weinberg, Phys. Rev. D 7, 1888-1910 (1973).
  * (155) L. Dolan and R. Jackiw, Phys. Rev. D 9, 3320-3341 (1974).
  * (156) M. E. Carrington, Phys. Rev. D 45, 2933-2944 (1992).
  * (157) C. L. Wainwright, Comput. Phys. Commun. 183, 2006-2013 (2012) [arXiv:1109.4189 [hep-ph]].
  * (158) G. Aad et al. [ATLAS and CMS], JHEP 08, 045 (2016) [arXiv:1606.02266 [hep-ex]].
  * (159) K. Schmitz, JHEP 01, 097 (2021) [arXiv:2002.04615 [hep-ph]].
  * (160) A. Ringwald, K. Saikawa and C. Tamarit, JCAP 02, 046 (2021) [arXiv:2009.02050 [hep-ph]].
  * (161) S. L. Larson, W. A. Hiscock and R. W. Hellings, Phys. Rev. D 62, 062001 (2000) [arXiv:gr-qc/9909080 [gr-qc]].
  * (162) H. Kudoh, A. Taruya, T. Hiramatsu and Y. Himemoto, Phys. Rev. D 73, 064006 (2006) [arXiv:gr-qc/0511145 [gr-qc]].
  * (163) K. Yagi and N. Seto, Phys. Rev. D 83, 044011 (2011) [erratum: Phys. Rev. D 95, no.10, 109901 (2017)] [arXiv:1101.3940 [astro-ph.CO]].
  * (164) E. Thrane and J. D. Romano, Phys. Rev. D 88, no.12, 124032 (2013) [arXiv:1310.5300 [astro-ph.IM]].
  * (165) C. J. Moore, R. H. Cole and C. P. L. Berry, Class. Quant. Grav. 32, no.1, 015014 (2015) [arXiv:1408.0740 [gr-qc]].
  * (166) S. Kuroyanagi, K. Nakayama and J. Yokoyama, PTEP 2015, no.1, 013E02 (2015) [arXiv:1410.6618 [astro-ph.CO]].
  * (167) K. Saikawa and S. Shirai, JCAP 05, 035 (2018) [arXiv:1803.01038 [hep-ph]].
  * (168) T. Robson, N. J. Cornish and C. Liu, Class. Quant. Grav. 36, no.10, 105011 (2019) [arXiv:1803.01944 [astro-ph.HE]].
  * (169) Y. Du, F. Huang, H. L. Li, Y. Z. Li and J. H. Yu, [arXiv:2111.01267 [hep-ph]].