# Supplemental Material:
Universal Bound on Energy Cost of Bit Reset in Finite Time

Yi-Zheng Zhen  Shenzhen Institute for Quantum Science and Engineering and
Department of Physics,Southern University of Science and Technology, Shenzhen
518055, China  Hefei National Laboratory for Physical Sciences at Microscale
and Department of Modern Physics, University of Science and Technology of
China, Hefei, 230026, China  Dario Egloff  Institute of Theoretical Physics,
Technische Universität Dresden, D-01062 Dresden, Germany  Max Planck Institute
for the Physics of Complex Systems, Nöthnitzer Strasse 38, 01187 Dresden,
Germany  Kavan Modi  School of Physics and Astronomy, Monash University,
Clayton, Victoria 3800, Australia  Shenzhen Institute for Quantum Science and
Engineering and Department of Physics,Southern University of Science and
Technology, Shenzhen 518055, China  Oscar Dahlsten<EMAIL_ADDRESS>Shenzhen Institute for Quantum Science and Engineering and Department of
Physics,Southern University of Science and Technology, Shenzhen 518055, China

##  I Proof of the work penalty lower bound

###  I.1 Step I: Geometric formulation of work penalty

For the relative entropy

$D\left[p\|\gamma\right]=\sum_{i}p_{i}\ln\frac{p_{i}}{\gamma_{i}},$ (1)

here we show that the two parts in its time derivative

$\dot{D}\left[p\|\gamma\right]=\dot{p}\partial_{p}D\left[p\|\gamma\right]+\dot{\gamma}\partial_{\gamma}D\left[p\|\gamma\right]$
(2)

correspond to entropy production rate and work penalty rate, respectively.

For the first term,

$\displaystyle\dot{p}\partial_{p}D\left[p\|\gamma\right]$
$\displaystyle=\sum_{i}\dot{p_{i}}\partial_{p_{i}}D\left[p\|\gamma\right]$ (3)
$\displaystyle=\sum_{i}\dot{p_{i}}\left(\ln\frac{p_{i}}{\gamma_{i}}+1\right)$
(4) $\displaystyle=\sum_{i}\dot{p_{i}}\left(\ln p_{i}+\beta E_{i}+\ln
Z\right)$ (5) $\displaystyle=\sum_{i}\dot{p_{i}}\ln
p_{i}+\beta\sum_{i}\dot{p_{i}}E_{i}$ (6)
$\displaystyle=-k_{B}^{-1}\dot{S}+\beta\dot{Q}$ (7)
$\displaystyle=-k_{B}^{-1}\dot{\Sigma},$ (8)

where we have used the thermal state expression $\gamma_{i}=\exp(-\beta
E_{i})/Z$, the definition of entropy production $\Sigma=\Delta S-Q/T$, and the
fact that $\sum_{i}\dot{p}_{i}=0$.

For the second term, recall that $W_{\rm pn}=W-W_{\rm qs}$ with
$W(\tau)=\int_{0}^{\tau}dt\sum_{i}p_{i}\dot{E_{i}}$ and $W_{\rm
qs}(\tau)=-k_{B}T\ln Z(\tau)/Z(0)$. We have

$\displaystyle\dot{W}_{\rm pn}$
$\displaystyle=\sum_{i}p_{i}\dot{E_{i}}+k_{B}T\frac{d}{dt}\ln Z$ (9)
$\displaystyle=k_{B}T\sum_{i}p_{i}\frac{d}{dt}\left(\beta E_{i}+\ln Z\right)$
(10) $\displaystyle=-k_{B}T\sum_{i}p_{i}\frac{d}{dt}\ln\frac{e^{-\beta
E_{i}}}{Z}$ (11)
$\displaystyle=-k_{B}T\sum_{i}p_{i}\frac{d}{dt}\ln\gamma_{i},$ (12)

such that the second term in Eq.(2) is

$\dot{\gamma}\partial_{\gamma}D\left[p\|\gamma\right]=-\sum_{i}p_{i}\frac{d}{dt}\ln\gamma_{i}=\beta\dot{W}_{\rm
pn}.$ (13)

Substituting Eqs.(8) and (13) into Eq.(2), we obtain

$\dot{D}\left[p\|\gamma\right]=-k_{B}^{-1}\dot{\Sigma}+\beta\dot{W}_{\rm pn}.$
(14)

The integration over time $\tau$ yields

$W_{\rm pn}\left(\tau\right)=k_{B}T\Delta
D\left[p\|\gamma\right]+T\Sigma\left(\tau\right),$ (15)

with $\Delta D[p\|\gamma]=D[p(\tau)\|\gamma(\tau)]-D[p(0)\|\gamma(0)]$.

##### Generalization to continuous case

Making the physically justified assumptions that both $p$ and $\gamma$ are in
Schwartz’ space (that is, the probabilities are reasonably well localised) and
that $\gamma$ has full support (the energy does not diverge), one directly
sees that simply replacing the sum by the corresponding integral will give the
same result for the continuous case. Of course there are well-motivated, but
strictly speaking unphysical, models that do not follow the above assumptions
and more care must be taken in those cases. For instance, if the energy does
diverge in some volume in space at some time, then $p$ needs to go to zero
there before that happens, since otherwise the model introduces infinite work
on the system. See Ref. [1] for details.

###  I.2 Step 2: Coarse-graining relative entropy, entropy production and
dynamics

We apply the coarse-graining to a generic system to create a bit. Precisely,
if the system is observed in one of states in
$\Omega_{0}=\\{i_{1},i_{2},\dots,i_{n}\\}$, we say the bit is in logical state
0; otherwise, if the system is in $\Omega_{1}=\\{j_{1},j_{2},\dots,j_{m}\\}$,
we say the bit is in logical state 1. For the bit to be logical, $\Omega_{0}$
and $\Omega_{1}$ must be disjoint and complementary. Then, the distribution of
bit states is denoted by $P^{\rm bit}=[P_{0}^{\rm bit},P_{1}^{\rm bit}]$ with

$P_{a}=\sum_{i\in\Omega_{a}}p_{i},\quad\forall a\in\left\\{0,1\right\\}.$ (16)

The master equation of coarse-grained states can be written as [2]

$\frac{d}{dt}P_{a}^{\rm bit}\left(t\right)=\Gamma^{\rm
bit}_{a\overline{a}}\left(t\right)P_{\overline{a}}^{\rm
bit}\left(t\right)-\Gamma^{\rm bit}_{\overline{a}a}\left(t\right)P_{a}^{\rm
bit}\left(t\right),$ (17)

where $\overline{a}$ is the negation of bit value $a\in\\{0,1\\}$, and
$\Gamma^{\rm
bit}_{a\overline{a}}=\sum_{i\in{\Omega}_{a}}\sum_{j\in{\Omega}_{\overline{a}}}\Gamma_{ij}p_{j}/P_{\overline{a}}$
can be shown to be the transition rate between coarse-grained states.

####  I.2.1 Relative entropy

We directly apply the log-sum inequality

$\sum_{i}u_{i}\ln\frac{u_{i}}{v_{i}}\geqslant\left(\sum_{i}u_{i}\right)\ln\frac{\sum_{i}u_{i}}{\sum_{i}v_{i}},\quad\forall
u_{i},v_{i}\geqslant 0,$ (18)

where the equality holds when $u_{i}/v_{i}$ is a constant for each $i$. The
relative entropy can be bounded as

$\displaystyle D[p\|\gamma]$
$\displaystyle=\sum_{i}p_{i}\ln\frac{p_{i}}{\gamma_{i}}=\sum_{a}\sum_{i\in\Omega_{a}}p_{i}\ln\frac{p_{i}}{\gamma_{i}}$
(19)
$\displaystyle\geqslant\sum_{a}\left(\sum_{i\in\Omega_{a}}p_{i}\right)\ln\frac{\sum_{i\in\Omega_{a}}p_{i}}{\sum_{i\in\Omega_{a}}\gamma_{i}}$
(20) $\displaystyle=\sum_{a}P_{a}^{\rm bit}\ln\frac{P_{a}^{\rm
bit}}{\gamma_{a}^{\rm bit}}$ (21) $\displaystyle=D\left[P^{\rm
bit}\|\gamma^{\rm bit}\right].$ (22)

For convenience, we use

$D_{\epsilon}(\tau)=D\left[P^{\rm bit}(\tau)\|\gamma^{\rm bit}(\tau)\right]$
(23)

to denote the relative entropy between the final state and final thermal
state, where $\epsilon=P_{1}^{\rm bit}(\tau)$.

####  I.2.2 Entropy production

From the master equation

$\dot{p}_{i}=\sum_{j\left(\neq i\right)}\Gamma_{ij}p_{j}-\Gamma_{ji}p_{i},$
(24)

where $\Gamma_{ij}$ are transition rates satisfying $\Gamma_{ij}\geqslant 0$
for $i\neq j$ and $\sum_{i}\Gamma_{ij}=0$ for all $j$, the entropy production
has a rate

$\displaystyle\dot{\Sigma}$ $\displaystyle=\dot{S}-\frac{\dot{Q}}{T}$ (25)
$\displaystyle=-k_{B}\sum_{i}\dot{p}_{i}\left(\ln p_{i}+\beta E_{i}\right)$
(26) $\displaystyle=k_{B}\sum_{i\neq
j}\left(\Gamma_{ij}p_{j}-\Gamma_{ji}p_{i}\right)\ln\frac{\gamma_{i}}{p_{i}}$
(27) $\displaystyle=\frac{k_{B}}{2}\sum_{i\neq
j}\left(\Gamma_{ij}p_{j}-\Gamma_{ji}p_{i}\right)\ln\frac{\gamma_{i}p_{j}}{p_{i}\gamma_{j}}$
(28) $\displaystyle=\frac{k_{B}}{2}\sum_{i\neq
j}\left(\Gamma_{ij}p_{j}-\Gamma_{ji}p_{i}\right)\ln\frac{\Gamma_{ij}p_{j}}{\Gamma_{ji}p_{i}}.$
(29)

Here, we have used the thermal state expression $\gamma_{i}=e^{-\beta
E_{i}}/Z$ in the third line and the detailed balance condition
$\Gamma_{ij}\gamma_{j}=\Gamma_{ji}\gamma_{i}$ in the last line.

Applying the time derivative on Eq.(16) gives

$\displaystyle\dot{P}_{a}^{\rm bit}$
$\displaystyle=\sum_{i\in\Omega_{a}}\dot{p}_{i}$ (30)
$\displaystyle=\sum_{i\in\Omega_{a}}\sum_{j\left(\neq
i\right)}\left(\Gamma_{ij}p_{j}-\Gamma_{ji}p_{i}\right)$ (31)
$\displaystyle=\sum_{i\in\Omega_{a}}\sum_{j\left(\neq
i\right)}\left(\Gamma_{ij}\frac{p_{j}}{P_{\overline{a}}^{\rm
bit}}P_{\overline{a}}^{\rm bit}-\Gamma_{ji}\frac{p_{i}}{P_{a}^{\rm
bit}}P_{a}^{\rm bit}\right)$ (32) $\displaystyle=\Gamma^{\rm
bit}_{a\overline{a}}P_{\overline{a}}^{\rm bit}-\Gamma^{\rm
bit}_{\overline{a}a}P_{a}^{\rm bit}.$ (33)

Here,

$\Gamma^{\rm
bit}_{a\overline{a}}=\sum_{i\in\Omega_{a}}\sum_{j\in\Omega_{\overline{a}}}\Gamma_{ij}\frac{p_{j}}{P_{\overline{a}}^{\rm
bit}}$ (34)

are transition rates between coarse-grained states, i.e. bit states 0 and 1.

With the coarse-grained transition rates, we can define the coarse-grained
entropy production as a similar form of Eq.(29), i.e.

$\dot{\Sigma}^{\rm bit}=\frac{k_{B}}{2}\sum_{a}\left(\Gamma^{\rm
bit}_{a\overline{a}}P^{\rm bit}_{\overline{a}}-\Gamma^{\rm
bit}_{\overline{a}a}P^{\rm bit}_{a}\right)\ln\frac{\Gamma^{\rm
bit}_{a\overline{a}}P^{\rm bit}_{\overline{a}}}{\Gamma^{\rm
bit}_{\overline{a}a}P^{\rm bit}_{a}}.$ (35)

Now we show that $\Sigma\geqslant\Sigma^{\rm bit}$ by proving
$\dot{\Sigma}\geqslant\dot{\Sigma}^{\rm bit}\geqslant 0$.

In Eq.(29), we can always drop the terms in the sum with $i,j$ in the same
$\Omega_{a}$. This is due to the fact that $\Gamma_{ij}p_{j}\geqslant 0$
$\forall i\neq j$ and therefore
$(\Gamma_{ij}p_{j}-\Gamma_{ji}p_{i})\ln(\Gamma_{ij}p_{j}/\Gamma_{ji}p_{i})\geqslant
0$. Then, using the log sum inequality again, we have

$\displaystyle\dot{\Sigma}$
$\displaystyle\geqslant\frac{k_{B}}{2}\sum_{a}\sum_{i\in\Omega_{a}}\sum_{j\in\Omega_{\overline{a}}}\left(\Gamma_{ij}p_{j}-\Gamma_{ji}p_{i}\right)\ln\frac{\Gamma_{ij}p_{j}}{\Gamma_{ji}p_{i}}$
(36)
$\displaystyle\geqslant\frac{k_{B}}{2}\sum_{a}\left(\sum_{i\in\Omega_{a},j\in\Omega_{\overline{a}}}\Gamma_{ij}p_{j}\ln\frac{\sum_{i\in\Omega_{a},j\in\Omega_{\overline{a}}}\Gamma_{ij}p_{j}}{\sum_{i\in\Omega_{a},j\in\Omega_{\overline{a}}}\Gamma_{ji}p_{i}}\right.$
(37)
$\displaystyle\qquad\qquad+\left.\sum_{i\in\Omega_{a},j\in\Omega_{\overline{a}}}\Gamma_{ji}p_{i}\ln\frac{\sum_{i\in\Omega_{a},j\in\Omega_{\overline{a}}}\Gamma_{ji}p_{i}}{\sum_{i\in\Omega_{a},j\in\Omega_{\overline{a}}}\Gamma_{ij}p_{j}}\right)$
(38) $\displaystyle=\frac{k_{B}}{2}\sum_{a}\left(\Gamma^{\rm
bit}_{a\overline{a}}P^{\rm bit}_{\overline{a}}-\Gamma^{\rm
bit}_{\overline{a}a}P^{\rm bit}_{a}\right)\ln\frac{\Gamma^{\rm
bit}_{a\overline{a}}P^{\rm bit}_{\overline{a}}}{\Gamma^{\rm
bit}_{\overline{a}a}P^{\rm bit}_{a}}$ (39) $\displaystyle=\dot{\Sigma}^{\rm
bit}.$ (40)

Note that $\dot{\Sigma}\geqslant\dot{\Sigma}^{\rm bit}\geqslant 0$, as each
term in the sum is non-negative due to the relation $(u-v)\log(u/v)\geqslant
0$ for $u,v>0$. This also recovers the second law as it shows that the entropy
production and its rate is non-negative. In the bit reset case, since
$\Sigma(0)=\Sigma^{\rm bit}(0)=0$ at the beginning, we finally obtain

$\Sigma(\tau)\geqslant\Sigma^{\rm bit}(\tau),$ (41)

i.e., the coarse-grained entropy production is lower than the fine-grained
one.

####  I.2.3 Tightness of the bound

We have show that the work penalty can be lower bounded by

$W_{\rm pn}\left(\tau\right)\geqslant k_{B}TD\left[P^{\rm
bit}\left(\tau\right)\|\gamma^{\rm bit}\left(\tau\right)\right]+T\Sigma^{\rm
bit}\left(\tau\right).$ (42)

As we used log-sum inequality in the proof, the minimal work penalty is
obtained by letting $p_{i}/\gamma_{i}$ be a constant for all $i$ in the same
coarse-grained states, i.e.

$\frac{p_{i}}{\gamma_{i}}=c_{a},\quad\forall i\in\Omega_{a}.$ (43)

As we will show in the next subsection, this equation leads to the assumption
of local equilibrium, based on which Eq. (42) becomes an inequality.

###  I.3 Effective two-level system in local equilibrium

We show that by assuming the local equilibrium, the coarse-grained master
equation respects detailed balance, such that the coarse-grained bit is
reduced to an effective two-level system. The local equilibrium assumption can
be written as

$\frac{p_{i}}{\gamma_{i}}=\frac{p_{j}}{\gamma_{j}}$ (44)

whenever $i$ and $j$ are in the same composite bit. Calling $\Omega_{0}$ the
first bit and $\Omega_{1}$ the second, this is the case if either
$i,j\in\Omega_{0}$ or $i,j\in\Omega_{1}$. For this special case of local
equilibrium, we have a very simple formula for the thermal state $\gamma^{\rm
bit}$ of the coarse grained bits, namely

$\gamma^{\rm bit}_{a}:=\sum_{i\in\Omega_{a}}\gamma_{i}.$ (45)

For the thermal state, we need to have that it is stationary under the
(undriven) master equation, such that

$\frac{\partial}{\partial t}\gamma^{\rm bit}_{0}=0.$ (46)

As the coarse grained system only has two levels, this is also the only
condition (unless the dynamics is trivial). Using the very useful property
that

$\displaystyle\frac{P^{\rm bit}_{a}}{\gamma^{\rm
bit}_{a}}=\frac{\sum_{i\in\Omega_{a}}p_{i}}{\sum_{i\in\Omega_{a}}\gamma_{i}}=\frac{\sum_{i\in\Omega_{a}}\frac{p_{k}}{\gamma_{k}}\gamma_{i}}{\sum_{i\in\Omega_{a}}\gamma_{i}}=\frac{p_{k}}{\gamma_{k}},$
(47)

for any $k\in\Omega_{a}$, we get indeed

$\displaystyle\frac{\partial}{\partial t}\gamma^{\rm bit}_{0}$
$\displaystyle=\Gamma^{\rm bit}_{00}\gamma_{0}^{\rm bit}+\Gamma^{\rm
bit}_{01}\gamma_{1}^{\rm bit}$
$\displaystyle=\sum_{i\in\Omega_{0}}\left(\sum_{j\in\Omega_{0}}\Gamma_{ij}\frac{p_{j}}{\sum_{k\in\Omega_{0}}p_{k}}\gamma^{\rm
bit}_{0}\right.$
$\displaystyle\qquad\qquad\left.+\sum_{j\in\Omega_{1}}\Gamma_{ij}\frac{p_{j}}{\sum_{k\in\Omega_{1}}p_{k}}\gamma^{\rm
bit}_{1}\right)$
$\displaystyle=\sum_{i\in\Omega_{0}}\left(\sum_{j\in\Omega_{0}}\Gamma_{ij}\gamma_{j}+\sum_{j\in\Omega_{1}}\Gamma_{ij}\gamma_{j}\right)$
$\displaystyle=\sum_{i\in\Omega_{0}}\frac{\partial}{\partial t}\gamma_{i}=0,$
(48)

as the $\gamma_{i}$ are the stationary components of the fine grained thermal
state. This also means that the $\gamma^{\rm bit}_{a}$ satisfy the equilibrium
condition (i.e. the detailed balance)

$\displaystyle\gamma^{\rm bit}_{0}=\frac{\Gamma^{\rm bit}_{01}}{\Gamma^{\rm
bit}_{01}-\Gamma^{\rm bit}_{00}}.$ (49)

With this in mind we can directly see that the coarse grained entropy
production is the same as the fine grained one. We start with the observation
that

$\displaystyle\ln\frac{\Gamma_{10}^{\rm bit}P_{0}}{\Gamma_{01}^{\rm
bit}P_{1}}=\ln\frac{P_{0}(\Gamma_{10}-\Gamma_{00})}{\Gamma_{01}^{\rm
bit}}+\ln\frac{\Gamma_{10}^{\rm bit}}{P_{1}(\Gamma_{10}-\Gamma_{00})}$ (50)
$\displaystyle=\ln\frac{P_{0}}{\gamma^{\rm
bit}_{0}}-\ln\frac{P_{1}}{\gamma^{\rm bit}_{1}}$ (51)

and, from Eq.(33), we have that

$\displaystyle\dot{P}_{a}$ $\displaystyle=\Gamma^{\rm
bit}_{a\overline{a}}P_{\overline{a}}-\Gamma^{\rm
bit}_{\overline{a}a}P_{a}=-\dot{P}_{\overline{a}}.$ (52)

It follows that

$\displaystyle\frac{d}{dt}\Sigma^{\rm bit}$ $\displaystyle=\frac{k_{{\rm
B}}}{2}\sum_{a\in\left\\{0,1\right\\}}\left(\Gamma_{a\overline{a}}^{{\rm
bit}}P_{\overline{a}}-\Gamma_{\overline{a}a}^{{\rm
bit}}P_{a}\right)\ln\frac{\Gamma_{a\overline{a}}^{\rm
bit}P_{\overline{a}}}{\Gamma_{\overline{a}a}^{\rm bit}P_{a}}$ (53)
$\displaystyle=k_{{\rm
B}}\sum_{a\in\left\\{0,1\right\\}}\dot{P}_{a}\ln\frac{P_{a}}{\gamma^{\rm
bit}_{a}}$ (54) $\displaystyle=k_{{\rm
B}}\sum_{a\in\left\\{0,1\right\\}}\sum_{i\in\Omega_{a}}\dot{p}_{i}\ln\frac{p_{i}}{\gamma_{i}}$
(55) $\displaystyle=k_{{\rm B}}\sum_{i}\dot{p}_{i}\ln\frac{p_{i}}{\gamma_{i}}$
(56) $\displaystyle=\frac{d}{dt}\Sigma,$ (57)

where the first line is the definition, in the second line we used Eq.(50) and
(52), in the third line we used linearity of the derivative and Eq.(47).

As for the relative entropy, let

$\frac{p_{i}}{\gamma_{i}}=c_{a},\quad\forall i\in\Omega_{a}.$ (58)

We can compare

$\displaystyle D\left[p\|\gamma\right]$
$\displaystyle=\sum_{a}\sum_{i\in\Omega_{a}}p_{i}\ln\frac{p_{i}}{\gamma_{i}}=\sum_{a}\sum_{i\in\Omega_{a}}p_{i}\ln
c_{a}$ (59) $\displaystyle=\sum_{a}\sum_{i\in\Omega_{a}}\gamma_{i}c_{a}\ln
c_{a}=\sum_{a}\gamma_{a}^{{\rm bit}}c_{a}\ln c_{a}$ (60)

with

$\displaystyle D\left[P^{{\rm bit}}\|\gamma^{{\rm bit}}\right]$
$\displaystyle=\sum_{a}P_{a}^{{\rm bit}}\ln\frac{P_{a}^{{\rm
bit}}}{\gamma_{a}^{{\rm bit}}}=\sum_{a}\gamma_{a}^{{\rm bit}}c_{a}\ln c_{a}$
(61) $\displaystyle=D\left[p\|\gamma\right].$ (62)

Therefore, for the system satisfying local equilibrium, we can prove that

$W_{\rm pn}(\tau)=D_{\epsilon}(\tau)+\Sigma^{\rm bit}(\tau)$ (63)

such that the penalty inequality is saturated.

###  I.4 Partial swap model

When the system is coarse-grained to a bit system, the dynamics of the bit can
always be described by the partial swap model. To see this, we introduce a
coarse-grained stationary state $P^{\rm st}$ such that

$\Gamma^{\rm bit}_{01}P^{\rm st}_{1}=\Gamma^{\rm bit}_{10}P^{\rm st}_{0},$
(64)

and let

$\mu=\Gamma^{\rm bit}_{01}+\Gamma^{\rm bit}_{10}.$ (65)

Then, we have $\Gamma^{\rm bit}_{01}=\mu P_{0}^{\rm st}$ and $\Gamma^{\rm
bit}_{10}=\mu P_{1}^{\rm st}$. Substitute it into Eq.(33) and we have

$\displaystyle\dot{P}^{\rm bit}_{0}$ $\displaystyle=\mu P_{0}^{\rm
st}P_{1}^{\rm bit}-\mu P_{1}^{\rm st}P_{0}^{\rm bit}=\mu\left(P_{0}^{\rm
st}-P_{0}^{\rm bit}\right)$ (66) $\displaystyle\dot{P}{\rm bit}$
$\displaystyle=\mu P_{1}^{\rm st}P_{0}^{\rm bit}-\mu P_{0}^{\rm st}P_{1}^{\rm
bit}=\mu\left(P_{1}^{\rm st}-P_{1}^{\rm bit}\right).$ (67)

Finally, Eq.(33) can be equivalently written as

$\dot{P}^{\rm bit}(t)=\mu(t)\left(P^{\rm st}(t)-P^{\rm bit}(t)\right).$ (68)

To see the meaning of partial swap, consider an infinitesimal time $dt$. A
state $P^{\rm bit}(t)$ evolves to $P^{\rm bit}(t+dt)=(1-\mu(t)dt)P^{\rm
bit}(t)+\mu(t)dtP^{\rm st}(t)$. That is, $P^{\rm bit}(t+dt)$ is obtained by
swapping the state $P^{\rm bit}(t)$ with the stationary state $P^{\rm st}$ at
a probability $\mu(t)dt$. Therefore, $\mu(t)$ is called the partial swap rate.

###  I.5 Step 3: Entropy production time-scaling via speed

Here, we prove the entropy production lower bound via the speed limit on the
partial swap model.

Take Eq.(68) into Eq.(35), we have the coarse-grained entropy production rate

$\displaystyle\dot{\Sigma}^{\rm bit}$ $\displaystyle=\mu
k_{B}\sum_{a\in\left(0,1\right)}\left(P_{a}^{\rm bit}-P_{a}^{\rm
st}\right)\ln\frac{P_{a}^{\rm bit}}{P_{a}^{\rm st}}$ (69)
$\displaystyle\geqslant 2\mu
k_{B}\sum_{a\in\left(0,1\right)}\frac{\left(P_{a}^{\rm bit}-P_{a}^{\rm
st}\right)^{2}}{P^{\rm bit}_{a}+P_{a}^{\rm st}},$ (70)

where we have used $(u-v)\ln(u/v)\geqslant 2(u-v)^{2}/(u+v)$ for any
$u,v\geqslant 0$. To obtain the speed limit, consider the norm-1 distance

$\displaystyle L\left[P^{\rm bit}\left(\tau\right),P^{\rm
bit}\left(0\right)\right]$ (71) $\displaystyle=$
$\displaystyle\sum_{a}\left|P^{\rm bit}_{a}\left(\tau\right)-P^{\rm
bit}_{a}\left(0\right)\right|=\sum_{a}\left|\int_{0}^{\tau}dt\dot{P}^{\rm
bit}_{a}\right|$ (72) $\displaystyle\leqslant$
$\displaystyle\int_{0}^{\tau}dt\sum_{a}\left|\dot{P}^{\rm
bit}_{a}\right|=\int_{0}^{\tau}dt\sum_{a}\left|\mu\left(P_{a}^{{\rm
st}}-P_{a}^{\rm bit}\right)\right|$ (73) $\displaystyle=$
$\displaystyle\int_{0}^{\tau}dt\sum_{a}\sqrt{\mu^{2}\frac{\left(P_{a}^{{\rm
st}}-P_{a}^{\rm bit}\right)^{2}}{P_{a}^{{\rm st}}+P^{\rm
bit}_{a}}}\sqrt{P_{a}^{{\rm st}}+P_{a}^{\rm bit}}$ (74)
$\displaystyle\leqslant$
$\displaystyle\int_{0}^{\tau}dt\sqrt{\mu^{2}\sum_{a}\frac{\left(P_{a}^{{\rm
st}}-P_{a}^{\rm bit}\right)^{2}}{P_{a}^{{\rm st}}+P^{\rm
bit}_{a}}}\sqrt{\sum_{a}\left(P_{a}^{{\rm st}}+P^{\rm bit}_{a}\right)}$ (75)
$\displaystyle=$ $\displaystyle\int_{0}^{\tau}dt\sqrt{\frac{\mu}{k_{{\rm
B}}}\dot{\Sigma}^{\rm bit}}$ (76) $\displaystyle\leqslant$
$\displaystyle\sqrt{\frac{\langle\mu\rangle_{\tau}\tau}{k_{{\rm
B}}}\Sigma^{\rm bit}(\tau)}.$ (77)

Here, we have used the Cauchy-Schwartz inequality twice, and
$\langle\mu\rangle_{\tau}=\tau^{-1}\int_{0}^{\tau}dt\mu(t)$ is the time
average of partial swap rate. The speed limit is therefore expressed as

$\tau\geqslant\frac{k_{B}L\left[P^{\rm bit}\left(\tau\right),P^{\rm
bit}\left(0\right)\right]^{2}}{\langle\mu\rangle_{\tau}\Sigma^{\rm bit}}.$
(78)

Combined with $\Sigma(\tau)\geqslant\Sigma^{\rm bit}(\tau)$, we have

$\Sigma(\tau)\geqslant\frac{k_{B}L\left[P^{\rm bit}\left(\tau\right),P^{\rm
bit}\left(0\right)\right]^{2}}{\langle\mu\rangle_{\tau}\tau}.$ (79)

###  I.6 Determining the partial swap rate

We now show how to determine or bound the coarse-grained) partial swap rate
$\mu(t)=\Gamma_{01}^{{\rm bit}}(t)+\Gamma_{10}^{{\rm bit}}(t)$ given the fine-
grained master equation. We firstly show how to determine it exactly for two
cases of master equations and then give an upper bound for the general case.

####  I.6.1 Case of partial swap

We consider firstly the case where the fine-grained system satisfies the
partial swap model, i.e.

$\dot{p}_{i}(t)=\mu(t)\left(\gamma_{i}(t)-p_{i}(t)\right),\qquad\text{for
}i=1,2,\dots,N.$ (80)

It can be written as the master equation

$\dot{p}_{i}(t)=\sum_{j}\Gamma_{ij}(t)p_{j}(t),\qquad\Gamma_{ij}(t)=\mu(t)\left(\gamma_{i}(t)-\delta_{ij}\right).$
(81)

Then, for the coarse-grained states

$P_{a}^{{\rm bit}}(t)=\sum_{i\in\Omega_{a}}p_{i}(t),$ (82)

the coarse-grained transition rates can be derived as

$\displaystyle\Gamma_{a\overline{a}}^{{\rm bit}}(t)$
$\displaystyle=\sum_{i\in\Omega_{a}}\sum_{j\in\Omega_{\overline{a}}}\frac{\Gamma_{ij}(t)p_{j}(t)}{P_{\overline{a}}^{{\rm
bit}}(t)}$ (83)
$\displaystyle=\mu(t)\sum_{i\in\Omega_{a}}\sum_{j\in\Omega_{\overline{a}}}\frac{\left(\gamma_{i}(t)-\delta_{ij}\right)p_{j}(t)}{P_{\overline{a}}^{{\rm
bit}}(t)}$ (84)
$\displaystyle=\mu(t)\sum_{i\in\Omega_{a}}\gamma_{i}\sum_{j\in\Omega_{\overline{a}}}\frac{p_{j}(t)}{P_{\overline{a}}^{{\rm
bit}}(t)}$ (85) $\displaystyle=\mu(t)\gamma_{a}^{{\rm bit}}(t),$ (86)
$\displaystyle\Gamma_{aa}^{{\rm bit}}(t)$
$\displaystyle=-\Gamma_{\overline{a}a}^{{\rm bit}}(t)$ (87)
$\displaystyle=\mu(t)\left(\gamma_{a}^{{\rm bit}}(t)-1\right),$ (88)

for $a\in{0,1}$. Therefore, the coarse-grained master equation can be
transferred to the partial swap model with the same partial swap rate $\mu$,
i.e.

$\dot{P}_{a}^{{\rm bit}}(t)=\mu(t)\left(\gamma_{a}^{{\rm bit}}(t)-P_{a}^{{\rm
bit}}(t)\right).$ (89)

The time averaged partial swap rate is
$\langle\mu\rangle_{\tau}=\tau^{-1}\int_{0}^{\tau}\mu(t)dt$.

####  I.6.2 Case of local equilibrium

In the second case, we consider the local equilibrium where the system state
in a coarse-grained set is proportional to the thermal state, i.e.
$p_{i}(t)/\gamma_{i}(t)=p_{j}(t)/\gamma_{j}(t)$ for $i,j\in\Omega_{a}$. Then,
the coarse-grained transition rates can be expressed as

$\displaystyle\Gamma_{a\overline{a}}^{{\rm bit}}(t)$
$\displaystyle=\sum_{i\in\Omega_{a}}\sum_{j\in\Omega_{\overline{a}}}\frac{\Gamma_{ij}(t)p_{j}(t)}{P^{\rm
bit}_{\overline{a}}(t)}$ (90)
$\displaystyle=\sum_{i\in\Omega_{a}}\sum_{j\in\Omega_{\overline{a}}}\frac{\Gamma_{ij}(t)\gamma_{j}(t)}{\gamma^{\rm
bit}_{\overline{a}}(t)}.$ (91)

The corresponding partial swap rate is then as

$\displaystyle\mu(t)$ $\displaystyle=\Gamma_{01}^{{\rm
bit}}(t)+\Gamma_{10}^{{\rm bit}}(t)$ (92)
$\displaystyle=\sum_{i\in\Omega_{0}}\sum_{j\in\Omega_{1}}\frac{\Gamma_{ij}(t)\gamma_{j}(t)}{\gamma^{\rm
bit}_{1}(t)}+\frac{\Gamma_{ij}(t)\gamma_{j}(t)}{\gamma^{\rm bit}_{0}(t)}.$
(93)

Therefore, $\mu(t)$ can be obtained if knowing the fine-grained thermal state
$\gamma_{i}(t)$ and the transition rates $\Gamma_{ij}(t)$. In this case, one
does not need to solve the master equation to obtain the system states
$p_{i}(t)$.

In fact, one can further simplify the derivation since the coarse-grained
transition rates satisfy the detailed balance. As we have shown in Sec. I C.,

$\Gamma_{01}^{{\rm bit}}\gamma_{1}^{\rm bit}(t)=\Gamma_{10}^{{\rm
bit}}\gamma_{0}^{\rm bit}(t).$ (94)

We define a time-dependent parameter $\omega(t)$ describing the thermalization
speed in terms of the coarse-grained transition rates by the expression

$\omega(t)=\Gamma_{a\overline{a}}^{\rm bit}(t)/\gamma_{a}^{\rm
bit}(t),\qquad\text{for }a=0,1.$ (95)

This gives the partial swap rate

$\mu(t)=\Gamma_{01}^{{\rm bit}}(t)+\Gamma_{10}^{{\rm bit}}(t)=\omega(t).$ (96)

Therefore, it is sufficient to obtain $\mu(t)$ by solely looking at
$\omega(t)=\Gamma_{01}^{\rm bit}(t)/\gamma_{0}^{\rm bit}(t)$, and the time
averaged partial swap rate is
$\langle\mu\rangle_{\tau}=\tau^{-1}\int_{0}^{\tau}\omega(t)dt$.

####  I.6.3 General discrete case

For the third case, we consider the general discrete system and upper bound
$\mu$ by the underlying dynamics, i.e.

$\displaystyle\mu(t)$ $\displaystyle=\Gamma_{01}^{{\rm
bit}}(t)+\Gamma_{10}^{{\rm bit}}(t)$ (97)
$\displaystyle=\sum_{i\in\Omega_{0}}\sum_{j\in\Omega_{1}}\frac{\Gamma_{ij}(t)p_{j}(t)}{P_{1}^{{\rm
bit}}(t)}+\frac{\Gamma_{ji}(t)p_{i}(t)}{P_{0}^{{\rm bit}}(t)}$ (98)
$\displaystyle\leqslant\sum_{i\in\Omega_{0}}\max_{j\in\Omega_{1}}\Gamma_{ij}(t)+\sum_{j\in\Omega_{1}}\max_{i\in\Omega_{0}}\Gamma_{ji}(t).$
(99)

In this manner, the time-averaged partial swap rate can be also upper bounded
in terms of the fine-grained master equation transition rates.

Finally we remark that for the continuous system, the transition rate
$\Gamma_{x,x^{\prime}}$, where $x,x^{\prime}$ are the continuous variable, may
diverge so one needs to consider a region around $x$. We will develop this
analysis in follow-up work.

##  II Time-scaling of relative entropy

As the partial swap model in Eq.(68) is a first-order partial differential
equation, we have a solution $P^{{\rm bit}}\left(t\right)=\left[P_{0}^{{\rm
bit}}\left(t\right),P_{1}^{{\rm bit}}\left(t\right)\right]$ with

$P_{1}^{{\rm bit}}\left(t\right)={\rm
e}^{-\left\langle\mu\right\rangle_{t}t}\left[P_{1}^{{\rm
bit}}\left(0\right)+\int_{0}^{t}ds\,{\rm
e}^{\left\langle\mu\right\rangle_{s}s}\mu\left(s\right)P_{1}^{{\rm
st}}\left(s\right)\right].$ (100)

During the bit reset, $P_{1}^{{\rm st}}\left(s\right)\geqslant P_{1}^{{\rm
st}}\left(t\right)$ for any two times $s\leqslant t$. Then, $P_{1}^{{\rm
bit}}\left(t\right)$ can be lower bounded by

$\displaystyle P_{1}^{{\rm{\rm bit}}}\left(t\right)$
$\displaystyle\geqslant{\rm
e}^{-\left\langle\mu\right\rangle_{t}t}P_{1}^{{\rm{\rm{\rm
bit}}}}\left(0\right)+{\rm e}^{-\left\langle\mu\right\rangle_{t}t}P_{1}^{{\rm
st}}\left(t\right)\int_{0}^{t}ds\,{\rm
e}^{\left\langle\mu\right\rangle_{s}s}\mu\left(s\right)$ (101)
$\displaystyle={\rm e}^{-\left\langle\mu\right\rangle_{t}t}P_{1}^{{\rm{\rm{\rm
bit}}}}\left(0\right)+\left(1-{\rm
e}^{-\left\langle\mu\right\rangle_{t}t}\right)P_{1}^{{\rm st}}\left(t\right),$
(102)

where $\int_{0}^{t}ds\,{\rm
e}^{\left\langle\mu\right\rangle_{s}s}\mu\left(s\right)={\rm
e}^{\left\langle\mu\right\rangle_{t}t}-1$ has been used. To bound
$D\left[P^{{\rm bit}}\|\gamma^{{\rm bit}}\right]$, we need to know how
$D\left[P^{{\rm bit}}\|\gamma^{{\rm bit}}\right]$ varies with $P_{1}^{{\rm
bit}}$. We prove the following relation.

Proposition._For three binary distributions $p=\left[p_{0},p_{1}\right]$,
$x=\left[x_{0},x_{1}\right]$ and $y=\left[y_{0},y_{1}\right]$, if
$p_{0}\geqslant p_{1}$ and $p_{1}\geqslant\alpha
x_{1}+\left(1-\alpha\right)y_{1}$ holds for some $\alpha\in\left[0,1\right]$,
then, for another binary distribution $q=\left[q_{0},q_{1}\right]$ with
$q_{1}\leqslant p_{1},$the following inequality holds,_

$\displaystyle D\left[p\|q\right]\geqslant$ $\displaystyle\alpha
D\left[x\|q\right]+\left(1-\alpha\right)D\left[y\|q\right]$
$\displaystyle+\left(\alpha^{2}-\alpha\right)\left(D\left[x\|y\right]+D\left[y\|x\right]\right).$
(103)

Proof. Note that $D\left[p\|q\right]$ is monotonicly increasing with $p_{1}$.
This can be checked by deriving $\partial_{p_{1}}D\left[p\|q\right]$ and
considering its boundaries $q_{1}\leqslant p_{1}\leqslant 1/2$. Then, let
$\beta=1-\alpha$ and we have

$\displaystyle D\left[p\|q\right]\geqslant$ $\displaystyle D\left[\alpha
x+\beta y\|q\right]$ $\displaystyle=$ $\displaystyle\sum_{k=0,1}\left(\alpha
x_{k}+\beta y_{k}\right)\log\frac{\alpha x_{k}+\beta y_{k}}{q_{k}}$
$\displaystyle=$ $\displaystyle\alpha\sum_{k}x_{k}\log\frac{\alpha x_{k}+\beta
y_{k}}{x_{k}}\frac{x_{k}}{q_{k}}$
$\displaystyle+\beta\sum_{k}y_{k}\log\frac{\alpha x_{k}+\beta
y_{k}}{y_{k}}\frac{y_{k}}{q_{k}}$ $\displaystyle=$ $\displaystyle-\alpha
D\left[x\|\alpha x+\beta y\right]+\alpha D\left[x\|q\right]$
$\displaystyle-\beta D\left[y\|\alpha x+\beta y\right]+\beta
D\left[y\|q\right]$ $\displaystyle\geqslant$ $\displaystyle\ \alpha
D\left[x\|q\right]+\beta D\left[y\|q\right]$
$\displaystyle-\alpha\beta\left(D\left[x\|y\right]+D\left[y\|x\right]\right).$

In the last inequality, we have used the convexity of relative entropy.
_[QED.]_

Now, substitute Eq.(102) into Eq.(103), we have

$\displaystyle D_{\epsilon}\left(\tau\right)=$ $\displaystyle D\left[P^{{\rm
bit}}\left(\tau\right)\|\gamma^{{\rm bit}}\left(\tau\right)\right]$
$\displaystyle\geqslant$ $\displaystyle qD\left[\gamma^{{\rm{\rm{\rm
bit}}}}\left(0\right)\|\gamma^{{\rm bit}}\left(\tau\right)\right]$
$\displaystyle+\left(1-q\right)D\left[P^{{\rm
st}}\left(\tau\right)\|\gamma^{{\rm bit}}\left(\tau\right)\right]$
$\displaystyle-q\left(1-q\right)J\left[\gamma^{\rm bit}\left(0\right),P^{{\rm
st}}\left(\tau\right)\right],$ (104)

where $q={\rm e}^{-\langle\mu\rangle_{\tau}\tau}$ and $J[r,s]=D[r\|s]+D[r\|s]$
is the symmetric relative entropy. We have the following cases:

  * •

Case 1: for small enough time $\tau$ such that $1-{\rm
e}^{-\left\langle\mu\right\rangle_{\tau}\tau}$ almost vanishes, up to the
first order of $\left(1-e^{-\left\langle\mu\right\rangle_{\tau}\tau}\right)$,
we have

$\displaystyle D_{\epsilon}\left(\tau\right)\geqslant$ $\displaystyle
e^{-\left\langle\mu\right\rangle_{\tau}\tau}D\left[\gamma^{{\rm{\rm{\rm
bit}}}}\left(0\right)\|\gamma^{{\rm bit}}\left(\tau\right)\right]$
$\displaystyle-\left(1-e^{-\left\langle\mu\right\rangle_{\tau}\tau}\right)D\left[\gamma^{{\rm{\rm{\rm
bit}}}}\left(0\right)\|P^{{\rm st}}\left(\tau\right)\right]$
$\displaystyle+{\cal
O}\left(\left(1-e^{-\left\langle\mu\right\rangle_{\tau}\tau}\right)^{2}\right).$
(105)

As in this case $1\gtrsim e^{-\left\langle\mu\right\rangle_{\tau}\tau}\gg
1-{\rm e}^{-\left\langle\mu\right\rangle_{\tau}\tau}$,
$D_{\epsilon}\left(\tau\right)$ is mainly in an exponential scaling.

  * •

Case 2: if $P^{{\rm st}}\left(\tau\right)=\gamma^{{\rm
bit}}\left(\tau\right)$, we have

$\displaystyle D_{\epsilon}\left(\tau\right)\geqslant$ $\displaystyle
q^{2}J\left[\gamma^{{\rm{\rm{\rm bit}}}}\left(0\right),\gamma^{{\rm
bit}}\left(\tau\right)\right]$ $\displaystyle-qD\left[\gamma^{{\rm
bit}}\left(\tau\right)\|\gamma^{{\rm{\rm{\rm bit}}}}\left(0\right)\right]$
(106) $\displaystyle\geqslant$ $\displaystyle
q\left(2q-1\right)D\left[\gamma^{{\rm
bit}}\left(0\right)\|\gamma^{{\rm{\rm{\rm bit}}}}\left(\tau\right)\right]$
(107)

where we have used $D[r\|s]\geqslant D[s\|r]$ for two binary distribution with
$1/2\geqslant r_{1}\geqslant s_{1}$.

To show this, one can consider a function
$g\left(r_{1}\right)=D\left[r\|s\right]-D\left[s\|r\right]$ and check its
derivatives on $r_{1}$. It yields $\partial_{r_{1}}g\geqslant 0$ and
$\partial_{r_{1}}^{2}g\geqslant 0$ for $s_{1}\leqslant r_{1}\leqslant 1/2$.
Thus, $g\left(r_{1}\right)\geqslant 0$ leads to $D\left[r\|s\right]\geqslant
D\left[s\|r\right]$.

##  III Constant-shifting protocol on a two-level system

As described in the main text, the constant-shifting protocol has $N$ steps.
Initially, $E_{0}(0)=E_{1}(0)=0$ and $P_{0}(0)=P_{1}(0)=1/2$. Then,
$E_{0}(t)=0$ is fixed during a finite time $\tau$, while at $t=(k-1)\tau/N$
for each $k=1,2,\dots,N$, the energy $E_{1}$ is instantly lifted with an
increment $E_{\max}/N$. The system undergoes a thermalization in the following
$\Delta\tau=\tau/N$ time. At the $k$-th step, the energy $E_{1}$ is denoted by
$E^{k}_{1}$, the thermal state is denoted by $\gamma^{k}$, and the system
state is denoted by $P^{k}$.

We have shown that the relative entropy change is equivalent to the sum of
work penalty and the negative entropy production. Thus, for each energy lift
step, entropy production is zero and relative entropy change is exactly work
penalty, while for each thermalization step, work penalty is zero and relative
entropy change is exactly the negative entropy production. To summarize, the
entire entropy production is

$\Sigma=\sum_{k=1}^{N}\Sigma^{k}=\sum_{k=1}^{N}D\left[P^{k-1}\|\gamma^{k}\right]-D\left[P^{k}\|\gamma^{k}\right],$
(108)

while the entire work penalty is

$\displaystyle W_{{\rm pn}}$ $\displaystyle=\sum_{k=1}^{N}W_{{\rm
pn}}^{k}=\sum_{k=1}^{N}D\left[P^{k-1}\|\gamma^{k}\right]-D\left[P^{k-1}\|\gamma^{k-1}\right]$
(109) $\displaystyle=D\left[P^{N}\|\gamma^{N}\right]+\Sigma.$ (110)

###  III.1 Relative entropy production bound with a fixed thermal state

We first look at the partial swap model $\dot{p}=\mu\left(\gamma-p\right)$
with a time-independent $\gamma$. This corresponds to each thermalsation step.
The $p(t)$ can be solved as

$p\left(t\right)=e^{-\mu t}p\left(0\right)+\left(1-e^{-\mu t}\right)\gamma.$
(111)

Using the convexity of relative entropy, we have

$D[p(t)\|\gamma]\leqslant e^{-\mu t}D[p(0)\|\gamma]$ (112)

For the later use, we also derive a lower bound on $D[p(t)\|\gamma]$. We check
the time derivative of $D\left[p\|\gamma\right]$

$\displaystyle\frac{d}{dt}D\left[p\|\gamma\right]$
$\displaystyle=\frac{d}{dt}\sum_{i}p_{i}\left(\ln p_{i}-\ln\gamma_{i}\right)$
$\displaystyle=\sum_{i}\frac{d}{dt}p_{i}\left(\ln p_{i}-\ln\gamma_{i}\right)$
$\displaystyle=\mu\sum_{i}\left(\mu_{i}-p_{i}\right)\ln\frac{p_{i}}{\gamma_{i}}$
$\displaystyle=-\mu\left(D\left[p\|\gamma\right]+D\left[\gamma\|p\right]\right)$
$\displaystyle\geqslant-2\mu D\left[p\|\gamma\right],$ (113)

where we have used $D\left[p\|\gamma\right]\geqslant D\left[\gamma\|p\right]$.
Because relative entropy is always non-negative, the integration yields

$D[p(t)\|\gamma]\geqslant e^{-2\mu t}D[p(0)\|\gamma].$ (114)

###  III.2 Upper bound of relative entropy

Let

$q=e^{-\mu\tau/N}.$ (115)

The $D(\tau)$ in the constant-shifting protocol can be upper bounded by

$\displaystyle D(\tau)$ $\displaystyle=D\left[P^{N}\|\gamma^{N}\right]$
$\displaystyle\leqslant qD\left[P^{N-1}\|\gamma^{N}\right]$
$\displaystyle\leqslant
q^{2}D\left[P^{N-2}\|\gamma^{N}\right]+q\left(1-q\right)D\left[\gamma^{N-1}\|\gamma^{N}\right]$
$\displaystyle\leqslant\dots$ $\displaystyle\leqslant
q^{N}D\left[P^{0}\|\gamma^{N}\right]+\sum_{k=1}^{N-1}q^{k}\left(1-q\right)D\left[\gamma^{N-k}\|\gamma^{N}\right]$
$\displaystyle\leqslant
q^{N}D\left[P^{0}\|\gamma^{N}\right]+\sum_{k=1}^{N-1}q^{k}\left(1-q\right)D\left[\gamma^{0}\|\gamma^{N}\right]$
$\displaystyle=\left[q^{N}+\sum_{k=1}^{N-1}q^{k}\left(1-q\right)\right]D\left[\gamma^{0}\|\gamma^{N}\right]$
$\displaystyle=qD\left[\gamma^{0}\|\gamma^{N}\right],$ (116)

where we have used $D\left[P^{k-1}\|\gamma^{k}\right]\geqslant
D\left[P^{k-1}\|\gamma^{k-1}\right]$,
$D\left[\gamma^{k}\|\gamma^{N}\right]\leqslant
D\left[\gamma^{0}\|\gamma^{N}\right]$, and $P^{0}=\gamma^{0}$. Therefore,
$D(\tau)$ has an exponential scaling

$D\left[P^{N}\|\gamma^{N}\right]\leqslant
e^{-\mu\frac{\tau}{N}}D\left[\gamma^{0}\|\gamma^{N}\right],$ (117)

when $N$ is fixed.

###  III.3 Exponential scaling of the relative entropy for fixed error

The Eq.(117) shows that $D(\tau)$ of the constant-shifting protocol is
exponentially decaying up to a factor $D[\gamma(0)\|\gamma(\tau)]$. This
factor is a constant for the case of fixed energy. Here, we also show that for
the case of fixed error, although $\gamma(\tau)$ is changing for different
$\tau$, still $D(\tau)$ decays exponentially.

Consider a protocol where the reset error is fixed. This means that for
different times the final energy of the state is dependent on the total time
$\tau$ of the protocol. Let $E_{\max}(\tau_{1})$ be the final energy for the
protocol with length $\tau$ and error $\epsilon$. We have shown that for the
same upper energy at longer total times the relative entropy term decays
exponentially. As the protocol is longer, additionally, the final state will
be nearer to the thermal state (that keeps fixed, as it only depends on the
final energy). The final error will thus be smaller than $\epsilon$ for that
protocol. This means that, if instead of the energy, we keep the final error
fixed, we can drive the energy of the level to be set to zero slower, which
will, at any time of the protocol with fixed error, drive the thermal state
slower to the reset state, in comparison to the case of the protocol with
fixed final energy. It then directly follows that for the protocol with fixed
error $\epsilon_{1}$ the final thermal state for any length of the protocol
$\tau_{2}>\tau_{1}$ will be farther away of the reset state, than for the
protocol with fixed maximal energy $E_{\max}$, while the final state will
obviously be the same (as the error is fixed). This means that the final state
will be nearer (by any distance measure) to the final thermal state.
Concluding, we find that the decay of the relative entropy term needs to be at
least as fast in the case of protocols with fixed error as in the case of
protocols with fixed final energy, where it is exponential; which yields the
claim.

###  III.4 Exponential dacay of work penalty

We check the total work cost

$W=\sum_{k=0}^{N-1}P_{1}^{k}{\cal E}.$ (118)

For each time interval of thermalization $\Delta\tau=\tau/N$, the partial swap
model leads to

$P_{1}^{k}=e^{-\mu\Delta\tau}P_{1}^{k-1}+\left(1-e^{-\mu\Delta\tau}\right)\gamma^{k}_{1}.$
(119)

The total work cost is then

$\displaystyle W$ $\displaystyle=\sum_{k=1}^{N}P_{1}^{k-1}{\cal E}$ (120)
$\displaystyle=\sum_{k=1}^{N}\gamma_{1}^{k-1}{\cal
E}+\sum_{k=1}^{N}\left(P_{1}^{k-1}-\gamma_{1}^{k-1}\right){\cal E}$ (121)
$\displaystyle=\sum_{k=1}^{N}\gamma_{1}^{k-1}{\cal
E}+e^{-\mu\Delta\tau}\sum_{k=2}^{N}\left(P_{1}^{k-2}-\gamma_{1}^{k-1}\right){\cal
E}.$ (122)

To calculate the exact value of $W$ is a difficult work. Instead, we calculate
the lower and upper bounds of $W$. The first term
$\sum_{k=1}^{N}\gamma_{1}^{k-1}{\cal E}$ can be bounded by (see Fig.1

$\displaystyle\sum_{k=1}^{N}\gamma_{1}^{k-1}{\cal E}$
$\displaystyle\geqslant\int_{0}^{E_{\max}}\gamma_{1}\left(E\right)dE$ (123)
$\displaystyle\sum_{k=1}^{N}\gamma_{1}^{k-1}{\cal E}$
$\displaystyle\leqslant\int_{0}^{E_{\max}}\gamma_{1}\left(E-{\cal
E}\right)dE.$ (124)

Substituting the above two equations into the work cost in Eq.(122), together
with $\gamma_{1}^{k}\leqslant P_{1}^{k}\leqslant 1/2$, we have bounds as

$\displaystyle W\geqslant$
$\displaystyle\int_{0}^{E_{\max}}\gamma_{1}\left(E\right)dE+e^{-\mu\Delta\tau}\sum_{k=2}^{N}\left(\gamma_{1}^{k-2}-\gamma_{1}^{k-1}\right){\cal
E},$ (125) $\displaystyle\geqslant$ $\displaystyle W_{\rm
qs}+e^{-\mu\Delta\tau}\left(P_{0}^{{\rm eq}}-P_{N-1}^{{\rm eq}}\right){\cal
E}$ (126) $\displaystyle=$ $\displaystyle W_{\rm
qs}+e^{-\mu\Delta\tau}\left(\frac{1}{2}-\frac{1}{1+e^{\beta\left(E_{\max}-{\cal
E}\right)}}\right){\cal E},$ (127)

$\displaystyle W\leqslant$
$\displaystyle\int_{0}^{E_{\max}}\gamma_{1}\left(E-{\cal
E}\right)dE+e^{-\mu\Delta\tau}\sum_{k=2}^{N}\left(\frac{1}{2}-\gamma_{1}^{k-1}\right){\cal
E}=\int_{-{\cal E}}^{E_{\max}-{\cal
E}}\gamma_{1}\left(E\right)dE+e^{-\mu\Delta\tau}\left[\frac{N-1}{2}-\sum_{k=2}^{N}\gamma_{1}^{k-1}\right]{\cal
E}$ (128) $\displaystyle\leqslant$ $\displaystyle W_{\rm
qs}+\left(\int_{-{\cal E}}^{0}-\int_{E_{\max}-{\cal
E}}^{E_{\max}}\right)\gamma_{1}\left(E\right)dE+e^{-\mu\Delta\tau}\left[\frac{E_{\max}-{\cal
E}}{2}-\int_{{\cal E}}^{E_{\max}}\gamma_{1}\left(E\right)dE\right]$ (129)
$\displaystyle=$ $\displaystyle W_{\rm
qs}+T\ln\left[\frac{1}{2}+\frac{1}{2}\frac{e^{\beta{\cal E}}+e^{-\beta
E_{\max}}}{1+e^{-\beta\left(E_{\max}-{\cal
E}\right)}}\right]+e^{-\mu\Delta\tau}\left[\frac{E_{\max}-{\cal E}}{2}-W_{\rm
qs}+\int_{0}^{{\cal E}}\gamma_{1}dE\right]$ (130) $\displaystyle=$
$\displaystyle W_{\rm qs}+e^{-\mu\Delta\tau}\left[\frac{E_{\max}}{2}-W_{\rm
qs}\right]+T\ln\left[\frac{1}{2}+\frac{1}{2}\frac{e^{\beta{\cal E}}+e^{-\beta
E_{\max}}}{1+e^{-\beta\left(E_{\max}-{\cal
E}\right)}}\right]+e^{-\mu\Delta\tau}\left(T\ln\frac{2}{1+e^{-\beta{\cal
E}}}-\frac{{\cal E}}{2}\right).$ (131)

The last term is always less than $0$. For the second term, let
$x=e^{\beta{\cal E}}\geqslant 1$ and $y=e^{-\beta E_{\max}}=x^{-N}$, using
$\ln x\leqslant x-1$,

$\displaystyle\ln\left[\frac{1}{2}+\frac{1}{2}\frac{e^{\beta{\cal
E}}+e^{-\beta E_{\max}}}{1+e^{-\beta\left(E_{\max}-{\cal E}\right)}}\right]$
$\displaystyle\leqslant$
$\displaystyle\frac{1}{2}\frac{x+y}{1+xy}-\frac{1}{2}=\frac{\left(x-1\right)\left(x^{N}-1\right)}{2\left(x^{N}+x\right)}$
(132) $\displaystyle\leqslant\frac{x-1}{2}.$ (133)

Finally, the work penalty can be characterized by $W_{\rm pn}^{L}\leqslant
W_{\rm pn}\leqslant W_{\rm pn}^{U}$ with

$\displaystyle W_{\rm pn}^{L}$
$\displaystyle=e^{-\mu\tau/N}\left[\frac{1}{2}-\frac{1}{1+e^{\beta\left(E_{\max}-{\cal
E}\right)}}\right]{\cal E}$ (134) $\displaystyle W_{\rm pn}^{U}$
$\displaystyle=T\frac{e^{\beta{\cal
E}}-1}{2}+e^{-\mu\tau/N}\left[\frac{E_{\max}}{2}-W_{\rm qs}\right].$ (135)

Therefore, the work penalty decays exponentially to a non-zero constant.

Figure 1: Restrict sum by integral.

##  IV The partial swap rate of overdamped Fokker-Plank equation

The dynamics of a Brownian particle in a double-well potential can be
described by the Fokker-Planck equation

$\partial_{t}p\left(x,t\right)=D\partial_{x}\left[\beta\partial_{x}V\left(x,t\right)+\partial_{x}\right]p\left(x,t\right).$
(136)

where $p(x,t)$ is the distribution of the Brownian particle position at time
$t$, $D$ is the diffusion constant, and $V(x,t)$ is the double-well potential.
Compared with the discrete case, the position $x$ casts as the continuous
states label and the potential $V(x,t)$ the internal energies. Therefore, by
writing $\sum_{i}p_{i}f_{i}$ as $\int_{-\infty}^{+\infty}dxp(x)f(x)$, where
$f$ can be internal energies (potentials) or stochastic entropies, we
generalize the definitions in the discrete case to the continuous. For
instance, for $V(x,t)$ at time $t$, the corresponding thermal distribution is

$\gamma\left(x,t\right)=\frac{e^{-\beta V\left(x,t\right)}}{Z},\quad
Z\left(t\right)=\int_{-\infty}^{+\infty}dxe^{-\beta V(x,t)}.$ (137)

Then, it can be verified that the work penalty relation

$\displaystyle W_{\rm pn}\left(\tau\right)$
$\displaystyle=W\left(\tau\right)-W_{\rm pn}$ (138) $\displaystyle=T\Delta
D\left[p\|\gamma\right]+T\Sigma\left(\tau\right),$ (139)

naturally holds.

To perform the coarse-graining analysis, it is necessary to write the Fokker-
Planck equation in the form of continuous master equation

$\partial_{t}p\left(x,t\right)=\int_{\mathbb{R}}dx^{\prime}\Gamma\left(x|x^{\prime};t\right)p\left(x^{\prime},t\right).$
(140)

Here, $\Gamma(x|x^{\prime};t)$ is the transition rate from $x^{\prime}$ to
$x$. It satisfies $\Gamma(x|x^{\prime};t)\geqslant 0$ $\forall x\neq
x^{\prime},t$ and $\int dx\Gamma(x|x^{\prime};t)=0$ $\forall t$. To write the
differential form of Eq.(136) to the integral form of Eq.(140), we consider a
small time interval

$\delta t=t-t^{\prime}.$ (141)

Because the evolution is Markovian, the distribution $p(x,t)$ can be viewed as
a one-dimensional random walk from $p(x,t^{\prime})$,

$p\left(x,t\right)=\int_{\mathbb{R}}dx^{\prime}P\left[x,t|x^{\prime},t^{\prime}\right]p\left(x^{\prime},t^{\prime}\right),$
(142)

where $P[x,t|x^{\prime},t^{\prime}]$ is the jump probability from position
$x^{\prime}$ at time $t^{\prime}$ to position $x$ at time $t$. Then, we can
approximate $\partial_{t}p(x,t)$ via the left derivative

$\frac{p\left(x,t\right)-p\left(x,t^{\prime}\right)}{\delta t}=\frac{1}{\delta
t}\int_{\mathbb{R}}dx^{\prime}\left(P\left[x,t|x^{\prime},t^{\prime}\right]-\delta_{xx^{\prime}}\right)p\left(x^{\prime},t^{\prime}\right).$
(143)

Compared with Eq.(140), we can approximate $\Gamma(x|x^{\prime};t)$ by

$\Gamma\left(x|x^{\prime};t\right)=\lim_{\delta t\rightarrow
0}\frac{P\left[x,t|x^{\prime},t^{\prime}\right]-\delta_{xx^{\prime}}}{\delta
t}.$ (144)

Here we assume that the above limitation is well-defined. It this is not the
case, we could always use $P[x,t|x^{\prime},t^{\prime}]$ of small time $\delta
t$ to evaluate $\Gamma(x|x^{\prime};t)$.

From the Fokker-Planck equation in Eq.(136), $P[x,t|x^{\prime},t^{\prime}]$
can be solved as [3]

$P\left[x,t|x^{\prime},t^{\prime}\right]=\frac{1}{\sqrt{4\pi D\delta
t}}e^{-\frac{\left(\delta x+\beta
D\partial_{x^{\prime}}V\left(x^{\prime},t^{\prime}\right)\delta
t\right)^{2}}{4D\delta t}}.$ (145)

Upon Eq.(144), we can verify that the detailed balanced condition holds, i.e.

$\Gamma\left(x|x^{\prime};t\right)\gamma\left(x^{\prime},t\right)=\Gamma\left(x^{\prime}|x;t\right)\gamma\left(x,t\right).$
(146)

This allows us to apply the coarse-grained entropy production in Eq.(35) as a
lower bound of the actual entropy production, as shown in Eq.(41). Here, the
coarse-grained states are

$\displaystyle P_{0}\left(t\right)$
$\displaystyle=\int_{-\infty}^{0}dxp\left(x,t\right),$ (147) $\displaystyle
P_{1}\left(t\right)$ $\displaystyle=\int^{+\infty}_{0}dxp\left(x,t\right).$
(148)

That is, when using the Brownian particle in double-well potential as a bit,
we choose position in the left well as bit value 0 while position in the right
well as bit value 1.

Analogy to the discrete case, the coarse-grained bit system satisfies mater
equation with transition rate

$\displaystyle\Gamma_{01}^{\rm bit}$
$\displaystyle=\int_{-\infty}^{0}du\int_{0}^{+\infty}dv\Gamma\left(u|v;t\right)\frac{p\left(v,t\right)}{P_{1}\left(t\right)},$
(149) $\displaystyle\Gamma_{10}^{\rm bit}$
$\displaystyle=\int_{0}^{+\infty}dv\int_{-\infty}^{0}du\Gamma\left(v|u;t\right)\frac{p\left(u,t\right)}{P_{0}\left(t\right)}.$
(150)

The dynamics can also be described by the partial swap model as in Eq.(68).
Particularly, the partial swap rate $\mu=\Gamma_{01}^{\rm
bit}+\Gamma_{10}^{\rm bit}$ can be well estimated.

##  V Implied bound on information throughput

As a two level system optimizes the work penalty, here we consider the bit is
constituted by truly two-level systems for the information throughput. By
substituting

$\displaystyle\ln\frac{Z(0)}{Z(\tau)}+D(\tau)$ (151) $\displaystyle=$
$\displaystyle\ln 2-\ln(1+e^{-\beta
E_{\max}})+\epsilon\ln\frac{\epsilon}{1/(1+e^{\beta E_{\max})}}$ (152)
$\displaystyle+(1-\epsilon)\ln\frac{1-\epsilon}{1/(1+e^{-\beta E_{\max}})}$
(153) $\displaystyle=$ $\displaystyle\ln
2-H_{b}(\epsilon)+\epsilon\ln\frac{1+e^{\beta E_{\max}}}{1+e^{-\beta
E_{\max}}}$ (154) $\displaystyle=$ $\displaystyle\ln
2-H_{b}(\epsilon)+\epsilon\beta E_{\max},$ (155)

where $H_{b}(\epsilon)=\epsilon\ln\epsilon+(1-\epsilon)\ln(1-\epsilon)$ is the
binary entropy, into the binary switch energy

$\displaystyle E_{\rm bit}$ $\displaystyle=W_{\rm qs}+W_{\rm pn}$ (156)
$\displaystyle=k_{B}T\ln\frac{Z(0)}{Z(\tau_{\rm SW})}+k_{B}TD(\tau_{\rm
SW})+T\Sigma(\tau_{\rm SW}),$ (157)

we obtain

$\displaystyle P$ $\displaystyle=\frac{E_{\rm bit}}{B}$ (158)
$\displaystyle\geqslant\frac{k_{B}Tn}{\tau_{\rm SW}}\left(\ln
2-H_{b}(\epsilon)+\epsilon\beta E_{\max}+\frac{(1-2\epsilon)^{2}}{\mu\tau_{\rm
SW}}\right).$ (159)

## References

  * Gardiner [2004] C. W. Gardiner,  _Handbook of Stochastic Methods: for Physics, Chemistry and the Natural Sciences_ , 3rd ed. (Springer, New York, 2004).
  * Esposito [2012] M. Esposito, Stochastic Thermodynamics under Coarse Graining, Phys. Rev. E 85, 041125 (2012).
  * Seifert [2012] U. Seifert, Stochastic Thermodynamics, Fluctuation Theorems and Molecular Machines, Rep. Prog. Phys. 75, 126001 (2012).


# Sterile Neutrino Portal Dark Matter with $Z_{3}$ Symmetry

An Liu1 Zhi-Long Han1<EMAIL_ADDRESS>Yi Jin1,2 Honglei Li1 1School of
Physics and Technology, University of Jinan, Jinan, Shandong 250022, China
2Guangxi Key Laboratory of Nuclear Physics and Nuclear Technology, Guangxi
Normal University, Guilin, Guangxi 541004, China

###### Abstract

In this paper, we consider the sterile neutrino portal dark matter with
$Z_{3}$ symmetry. This model further extends the canonical type-I seesaw with
a fermion singlet $\chi$ and a scalar singlet $\phi$. Under the $Z_{3}$
symmetry, the dark sector transforms as $\chi\to e^{i2\pi/3}\chi,\phi\to
e^{i2\pi/3}\phi$, while the standard model particles and the sterile neutrino
$N$ transform trivially. Besides the interactions as $y_{N}\phi\bar{\chi}N$
and $\lambda_{H\phi}(H^{\dagger}H)(\phi^{\dagger}\phi)$ allowed in the $Z_{2}$
symmetry, the $Z_{3}$ symmetry also introduces two new terms, i.e.,
$y_{\chi}\phi\overline{\chi^{c}}\chi$ and $\mu\phi^{3}/2$. These new
interactions induce additional semi-annihilation processes as $\chi\chi\to
N\chi$ and $\phi\phi\to h\phi$ for WIMP dark matter. We then perform a
comprehensive analysis of the phenomenology of this $Z_{3}$ symmetric model.
Viable parameter space is explored under the constraints from dark matter
relic density, Higgs invisible decay, indirect and direct detection for both
fermion and scalar dark matter. We find that the semi-annihilation channels
$\chi\chi\to N\chi$ and $\phi\phi\to N\chi$ can lead to quite different
phenomena from the $Z_{2}$ symmetric model, which provides a viable pathway to
distinguish these two kinds of model.

##  I Introduction

The identity of particle dark matter (DM) and the explanation for the tiny
mass of neutrinos remain outstanding questions in particle physics, garnering
attention as crucial topics in current research. Their common origin presents
the possibility of future exploration into new physics beyond the standard
model. The Weakly Interacting Massive Particle (WIMP) is the most promising
dark matter candidate Arcadi:2017kky . However, this scenario, such as the
extensively studied Higgs portal McDonald:1993ex ; Burgess:2000yq ;
Cline:2013gha  and $Z^{\prime}$ portal model Alves:2013tqa ; Alves:2015pea ;
Rodejohann:2015lca ; Klasen:2016qux , usually suffer stringent constraints
from direct detection. Therefore, a new interaction portal for the WIMP dark
sector should be considered.

Sterile neutrinos $N$ are introduced to generate the tiny neutrino mass via
the canonical Type-I seesaw mechanism Minkowski:1977sc ; Mohapatra:1979ia .
For a proper mixing angle with active neutrinos, the keV scale sterile
neutrino can serve as a decaying dark matter Dodelson:1993je ;
Boyarsky:2018tvu . Then the radiative decay $N\to\nu\gamma$ leads to an
observable signature at X-ray telescopes Roach:2019ctw , which is able to
explain the tentative 3.5 keV line signal Boyarsky:2014ska . However, the
parameter space for sterile neutrino dark matter now is tightly constrained
Roach:2022lgo . If the sterile neutrinos $N$ are charged under the dark group,
the lightest sterile neutrino becomes a stable dark matter. In this scenario,
the tree-level Type-I seesaw is also forbidden by the dark group, then light
neutrino mass could be generated via the radiative mechanism Ma:2006km ;
Ma:2007gq ; Ding:2016wbd ; Liu:2022byu .

On the other hand, the electroweak scale sterile neutrino $N$ is an ideal
messenger between the dark sector and the standard model Ballett:2019pyw ;
Blennow:2019fhy ; Falkowski:2009yz ; Falkowski:2011xh ; GonzalezMacias:2015rxl
; Gonzalez-Macias:2016vxy ; Liu:2020mxj ; Coito:2022kif ; Liu:2022rst . This
is facilitated through the new Yukawa coupling $y_{N}\phi\bar{\chi}N$, which
enables the secluded channel $\phi\phi/\chi\chi\to NN$, providing an
additional annihilation pathway for the WIMP dark matter. In particular, this
scenario features a relatively small nucleon scattering cross section, and
permits the indirect detection of observable gamma-ray signals Pospelov:2007mp
; Folgado:2018qlv ; Tang:2015coo ; Campos:2017odj ; Batell:2017rol ;
Batell:2017cmf . For an electroweak scale dark matter annihilating via the
sterile neutrino portal, it is still allowed by direct detection and is
hopefully probed by indirect detection in the near future. Meanwhile, the
sterile neutrino portal dark matter produced through the freeze-in mechanism
is also extensively studied Falkowski:2017uya ; Becker:2018rve ;
Chianese:2018dsz ; Bandyopadhyay:2020qpn ; Cheng:2020gut ; Chang:2021ose ;
Barman:2022scg ; Liu:2022cct .

The interactions between the dark sector and the standard model particles are
typically governed by the dark group, such as the well-studied $Z_{2}$
Escudero:2016ksa  or $U(1)_{B-L}$ symmetry Escudero:2016tzx . Although the
simplest $Z_{2}$ symmetry has demonstrated success in dark matter with a
simplified phenomenology, more sophisticated dark groups, e.g. $Z_{N}(N\geq
3)$ Batell:2010bp , $A_{4}$ Hirsch:2010ru  and $SU(2)$ symmetry Diaz-
Cruz:2010czr , are also options. For instance, an alternative explanation for
the observed relic abundance of dark matter through semi-annihilation may be
achieved by introducing the next simplest $Z_{3}$ symmetry, leading to a lower
bound on the direct detection cross section Belanger:2012vp ; Belanger:2012zr
.

In this paper, we consider the sterile neutrino portal dark matter with
$Z_{3}$ symmetry Bandyopadhyay:2022tsf ; Ghosh:2023ocl . This model includes a
fermion singlet $\chi$ and a scalar singlet $\phi$, both of which transform
non-trivially under the exact $Z_{3}$ symmetry as $\chi\to
e^{i2\pi/3}\chi,\phi\to e^{i2\pi/3}\phi$. The standard model particles and the
sterile neutrinos are not charged under the imposed $Z_{3}$ symmetry. Compared
with the $Z_{2}$ symmetry, the $Z_{3}$ symmetry allows two new interaction
terms, i.e., $y_{\chi}\phi\overline{\chi^{c}}\chi$ and $\mu\phi^{3}/2$, which
would lead to new annihilation channels of dark matter. The semi-annihilation
of fermion dark matter via the process $\chi\chi\to N\chi$ in the framework of
effective field theory has been considered in Ref. Bandyopadhyay:2022tsf .
Focusing on the self-interaction of dark scalar $\phi$, Ref. Ghosh:2023ocl
studies the non-thermal production of dark matter $\chi$ by the late time
decay $\phi\to\chi\nu$. In this paper, we perform a comprehensive analysis of
WIMP dark matter for both scalar and fermion scenarios. The non-thermal
production of dark matter will be considered in a separate paper Liu:2023hzl .

This paper is structured into several sections. In Sec. II, we introduce the
sterile neutrino portal dark matter model with $Z_{3}$ symmetry. In Sec. III,
we illustrate the evolution of the dark matter abundance under certain
scenarios. Then we perform a random scan to obtain the viable parameter space
for correct relic abundance. In Sec. IV, we calculate the branching ratio of
Higgs invisible decay. In Sec. V, we explore the indirect detection
constraints of dark matter. In Sec. VI, we compute the direct detection cross
section of dark matter. Finally, in Sec. VII, we provide concluding remarks
for our study.

##  II The Model

This model further extends the Type-I seesaw with a dark sector under the
$Z_{3}$ symmetry. Sterile neutrinos $N$ are introduced to generate tiny
neutrino mass. In this paper, we consider the electroweak scale $N$ in order
to accommodate WIMP dark matter. The dark sector consists of a scalar singlet
$\phi$ and a fermion singlet $\chi$. Under the dark $Z_{3}$ symmetry, the dark
sector transforms as $\chi\to e^{i2\pi/3}\chi,\phi\to e^{i2\pi/3}\phi$, while
the standard model particles and the sterile neutrinos transform trivially.
The lightest particle in the dark sector serves as dark matter. In this paper,
both fermion and scalar dark matter will be considered.

The Yukawa interaction takes the form of

$-\mathcal{L}_{Y}=\left(y_{\nu}\overline{L}\tilde{H}N+y_{N}\phi\bar{\chi}N+h.c.\right)+y_{\chi}\phi\overline{\chi^{c}}\chi,$
(1)

where $L$ is the left-handed lepton doublet and $H$ is the Higgs doublet with
$\tilde{H}=i\sigma_{2}H^{*}$. Light neutrino mass is generated by the Type-I
seesaw as

$m_{\nu}=-\frac{v^{2}}{2}y_{\nu}m_{N}^{-1}y_{\nu}^{T},$ (2)

where $v=246$ GeV is the vacuum expectation value of the Higgs field. For
electroweak scale sterile neutrinos, the mixing angle with light neutrino
$\theta$ is at the order of $\sqrt{m_{\nu}/m_{N}}\sim 10^{-6}$, which is far
below current collider limits Abdullahi:2022jlv .

The scalar potential under the exact $Z_{3}$ symmetry is

$V=-\mu_{H}^{2}H^{\dagger}H+\mu_{\phi}^{2}\phi^{\dagger}\phi+\lambda_{H}(H^{\dagger}H)^{2}+\lambda_{\phi}(\phi^{\dagger}\phi)^{2}+\lambda_{H\phi}(H^{\dagger}H)(\phi^{\dagger}\phi)+\left(\frac{\mu}{2}\phi^{3}+h.c.\right),$
(3)

where all the parameters are taken to be real. After the electroweak symmetry
breaking, the physical mass of dark scalar $\phi$ is
$m_{\phi}^{2}=\mu^{2}_{\phi}+\lambda_{H\phi}v^{2}/2$. The scalar potential in
Equation (3) must have a finite minimum to prevent unbounded energy, which
requires Belanger:2012zr

$\lambda_{H}>0,\quad\lambda_{\phi}>0,\quad\lambda_{H\phi}+2\sqrt{\lambda_{H}\lambda_{\phi}}>0.$
(4)

Meanwhile, the stability of electroweak vacuum sets an upper bound on the
cubic coupling $\mu$ as

$\mu\leq 2\sqrt{\lambda_{\phi}}m_{\phi},$ (5)

in the limit of small $\lambda_{H\phi}$. In order to maintain the validity of
perturbation theory, $\left|\lambda_{\phi}\right|\leqslant\pi$ and
$\left|\lambda_{H\phi}\right|\leqslant 4\pi$ should be further satisfied. In
the following studies, we assume $\mu<3m_{\phi}$, which is also allowed by the
unitarity constraints Hektor:2019ote .

##  III Relic density

Figure 1: The dominant annihilation channels of fermion dark matter. Panel (a)
is for the secluded channel $\chi\chi\to NN$, while panels (b)-(d) are for the
semi-annihilation channel $\chi\chi\to N\chi$.  Figure 2: The dominant
annihilation channels of scalar dark matter. Panel (a) and (b) denote the
annihilation channels to SM final states. Panel (c) is for the secluded
channel $\phi\phi\to NN$. Panels (d)-(f) are for the fermion semi-annihilation
channel $\phi\phi\to N\chi$. Panels (g)-(i) are for the scalar semi-
annihilation channel $\phi\phi\to h\phi$.

In this section, we first discuss the annihilation channels of dark matter. As
shown in Figure 1 , there are two dominant annihilation channels for the
fermion dark matter $\chi$. One is the secluded channel $\chi\chi\to NN$,
which also exists in the $Z_{2}$ symmetry model. The other one is the semi-
annihilation channel $\chi\chi\to N\chi$ DEramo:2010keq , which is induced by
the new Yukawa coupling $y_{\chi}\phi\overline{\chi^{c}}\chi$ under the
$Z_{3}$ symmetry. In principle, there is also the scalar semi-annihilation
channel $\chi\chi\to\phi h$, however the annihilation cross section is
$p$-wave suppressed for the Majorana-like Yukawa coupling
$y_{\chi}\phi\overline{\chi^{c}}\chi$ Cai:2015zza ; Guo:2021rre . For the
scalar dark matter $\phi$, there are four kinds of annihilation channels as
depicted in Figure 2 . Apart from the extensively studied Higgs portal
channels $\phi\phi\to\text{SM}$, the secluded channel $\phi\phi\to NN$ is also
allowed. Meanwhile, the cubic term $\mu\phi^{3}$ and the Yukawa coupling
$y_{\chi}\phi\overline{\chi^{c}}\chi$ induce two additional semi-annihilation
channels $\phi\phi\to N\chi$ and $\phi\phi\to h\phi$ if kinematically allowed.
Comparing with the simplest $Z_{3}$ scalar singlet dark matter Hektor:2019ote
, the fermion channel $\phi\phi\to N\chi$ is unique in this model. Therefore,
the sterile neutrino portal semi-annihilation channels $\chi\chi\to N\chi$ and
$\phi\phi\to N\chi$ provide a viable pathway to distinguish from the other
models. It is notable that when masses of the dark sector are nearly
degenerate, the co-annihilation channels as $\phi\chi\to h\chi/\phi N$ are
also possible. For simplicity, we do not consider such co-annihilation
channels in this paper.

Figure 3: The evolution of fermion dark matter abundance in different major
annihilation channels. The orange horizontal lines correspond to the Planck
observed abundance for $m_{\mathrm{DM}}=500$ GeV.

As the WIMP dark matter candidate, it is in thermal equilibrium at the very
beginning, and then decouples from the thermal bath when the temperature is
low enough. Defining the variable $z=m_{\text{DM}}/T$, the evolution of
fermion dark matter abundance $Y_{\chi}$ is determined by the Boltzmann
equation

$\displaystyle\frac{\mathrm{d}Y_{\chi}}{\mathrm{d}z}=-\frac{\lambda}{z^{2}}\langle\sigma
v\rangle_{\chi\chi\to
NN}\left(Y_{\chi}^{2}-(Y_{\chi}^{\mathrm{eq}})^{2}\right)-\frac{\lambda}{2z^{2}}\langle\sigma
v\rangle_{\chi\chi\to
N\chi}\left(Y_{\chi}^{2}-Y_{\chi}^{\mathrm{eq}}Y_{\chi}\right),$ (6)

where $\lambda$ is defined as $\lambda\equiv\sqrt{\pi
g_{*}/{45}}m_{\mathrm{DM}}M_{\mathrm{Pl}}$. Here, $g_{*}$ is the effective
number of degrees of freedom of the relativistic species and
$M_{\mathrm{Pl}}=1.2\times 10^{19}$ GeV is the Planck mass. The sterile
neutrino $N$ is assumed in thermal equilibrium Li:2022bpp . Similarly, the
evolution of scalar dark matter is calculated as

$\displaystyle\frac{\mathrm{d}Y_{\phi}}{\mathrm{d}z}$ $\displaystyle=$
$\displaystyle-\frac{\lambda}{z^{2}}\langle\sigma
v\rangle_{\phi\phi\to\mathrm{SM}}\left(Y_{\phi}^{2}-(Y_{\phi}^{\mathrm{eq}})^{2}\right)-\frac{\lambda}{2z^{2}}\langle\sigma
v\rangle_{\phi\phi\to
h\phi}\left(Y_{\phi}^{2}-Y_{\phi}^{\mathrm{eq}}Y_{\phi}\right)$
$\displaystyle-\frac{\lambda}{z^{2}}\langle\sigma v\rangle_{\phi\phi\to
NN}\left(Y_{\phi}^{2}-(Y_{\phi}^{\mathrm{eq}})^{2}\right)-\frac{\lambda}{2z^{2}}\langle\sigma
v\rangle_{\phi\phi\to
N\chi}\left(Y_{\phi}^{2}-\frac{(Y_{\phi}^{\mathrm{eq}})^{2}}{Y_{\chi}^{\mathrm{eq}}}Y_{\chi}\right).$

The thermally averaged cross sections $\langle\sigma v\rangle$ are calculated
with MicrOMEGAs Belanger:2018ccd .

Figure 4: Same as Figure 3, but for scalar dark matter.

Figure 3 and Figure 4 present the evolution of the dark matter relic abundance
via various annihilation channels during the early universe. According to the
general definition, the lighter of fermion $\chi$ and scalar $\phi$ is the
dark matter candidate. For illustration, we set the dark matter mass
$m_{\mathrm{DM}}=500$ GeV, the other heavier particle mass
$m_{\mathrm{Heavier}}=800$ GeV, and the sterile neutrino mass $m_{N}=180$ GeV.
The secluded channel $\chi\chi\to NN$ only involves the Yukawa coupling
$y_{N}$, whose impacts on the abundance is shown in panel (a) of Figure 3. The
contribution of semi-annihilation channel $\chi\chi\to N\chi$ is turned off
simply by setting $y_{\chi}=0$. To obtain correct relic abundance,
$y_{\chi}\sim\mathcal{O}(0.1)$ is required when $\chi\chi\to NN$ is the only
annihilation channel. For the process $\chi\chi\to N\chi$, both the Yukawa
coupling $y_{N}$ and $y_{\chi}$ contribute. We then fix $y_{N}=0.2$, and show
the impact of $y_{\chi}$ in panel (b). The observed relic abundance is
reproduced with $y_{\chi}\sim\mathcal{O}(1)$. Since $\chi\chi\to NN$ is also
kinematically allowed, it is clear that when $y_{\chi}\ll y_{N}$, the relic
abundance is actually determined by the secluded channel.

For the scalar dark matter $\phi$, we first show the impact of
$\lambda_{H\phi}$ on the canonical Higgs portal annihilation channels in panel
(a) of Figure 4. Contributions of other kinds of annihilation channels are
forbidden by fixing $y_{\chi}=y_{N}=0$ and $\mu=0$ GeV. These Higgs portal
channels are efficient to obtain the desired abundance with
$\lambda_{H\phi}\gtrsim\mathcal{O}(0.01)$. The contribution of scalar semi-
annihilation $\phi\phi\to h\phi$ is depicted in panel (b) while setting
$\lambda_{H\phi}=0.05$. A relatively large cubic coupling $\mu\gtrsim 100$ GeV
is required to make this channel the dominant one. In panels (c) to (f), we
consider the secluded channel $\phi\phi\to NN$ and fermion semi-annihilation
channel $\phi\phi\to N\chi$. Similar to the fermion dark matter scenario,
correct abundance is achieved with $y_{N}\sim\mathcal{O}(0.1)$ or
$y_{\chi}\sim\mathcal{O}(1)$ when $\phi\phi\to NN$ or $\phi\phi\to N\chi$ is
the dominant annihilation channel respectively. Different from the fermion
dark matter, the $s$-channel of semi-annihilation $\phi\phi\to N\chi$ is
induced by the cubic term $\mu\phi^{3}$ but not the Yukawa coupling
$y_{\chi}\phi\overline{\chi^{c}}\chi$. The contributions of the $s$\- and
$t/u$-channels are then separately shown in panels (d) and (e). Since the
contribution of $\phi\phi\to N\chi$ is suppressed by the final states phase
space in the benchmark points, a larger contribution is possible with lighter
$\chi,N$, which is illustrated in panel (f).

Generating the appropriate cosmological relic density, as determined with high
accuracy by the Planck experiment $\Omega_{\rm DM}h^{2}=0.120\pm 0.001$
Planck:2018vyg , is a crucial prerequisite for a viable dark matter candidate.
With a seesaw related mixing angle $\theta\sim\sqrt{m_{\nu}/m_{N}}$, the
lifetime of sterile neutrino $\tau_{N}$ would be longer than
$\mathcal{O}(0.1)$ s for $m_{N}<1$ GeV, which is excluded by Big Bang
Nucleosynthesis Ruchayskiy:2012si . So we assume $m_{N}>1$ GeV in this paper,
and then perform a random scan to explore the following dark sector parameter
space:

$y_{\chi,N}\in[10^{-4},1],\lambda_{H\phi}\in[10^{-6},1],\mu/m_{\phi}\in[0,3],m_{\chi,\phi}\in[1,10^{3}]~{}\text{GeV}.$
(8)

Figure 5: Distributions of samples with correct relic density for fermion dark
matter. The purple and blue points denote that the dominant annihilation
channel is $\chi\chi\to NN$ and $\chi\chi\to N\chi$ respectively.

Samples with correct relic density in the $3\sigma$ range of the Planck value
are kept for later study. Survived samples are then classified by the dominant
annihilation channels. The results are shown in Figure 5 and Figure 6 for the
case of fermion and scalar dark matter respectively. When $\chi\chi\to NN$ is
the dominant channel, it is clear that a lower bound on $y_{N}$ exists, which
is approximately $y_{N}\gtrsim(m_{\chi}/10^{4}~{}\text{GeV})^{1/2}$. Including
the contribution of semi-annihilation channel $\chi\chi\to N\chi$ would allow
$y_{N}$ to be about two orders of magnitude smaller. Similarly, the lower
limit $y_{\chi}\gtrsim m_{\chi}/10^{4}~{}\text{GeV}$ should be satisfied when
$\chi\chi\to N\chi$ is the dominant channel. It is also clear that for large
enough $y_{\chi}$, i.e., $y_{\chi}\gtrsim 0.3$, the semi-annihilation channel
will always have the largest contribution. These two channels are well
separated in the $y_{N}-y_{\chi}$ plane, namely $\chi\chi\to N\chi$ is the
leading one when $y_{\chi}\gtrsim 0.6y_{N}$. As shown in Figure 1 (b), there
is a $s$-channel for the $\chi\chi\to N\chi$ annihilation. Therefore, the
contribution of this channel is enhanced when $m_{\phi}\simeq 2m_{\chi}$,
which leads to the deep cusp in panel (d) of Figure 5. Apart from the
resonance region, the larger the ratio $m_{\phi}/m_{\chi}$ is, the bigger the
factor $\sqrt{y_{\chi}y_{N}}$ is required.

Figure 6:  Distributions of samples with correct relic density for scalar dark
matter. The purple, blue, yellow and green points denote that the dominant
annihilation channel is $\phi\phi\to NN$, $\phi\phi\to N\chi$, $\phi\phi\to$
SM, and $\phi\phi\to h\phi$ respectively.

For the scalar dark matter, $\lambda_{H\phi}\lesssim 0.1$ is enough to realize
the correct relic density when 10 GeV $\lesssim m_{\phi}\lesssim 10^{3}$ GeV.
The sharp dip around $m_{\phi}\sim m_{h}/2$ corresponds to the on-shell
production of $h$ in the $s$-channel, where $\lambda_{H\phi}$ can be as small
as $10^{-4}$. The SM channel dominant samples via the Higgs portal distribute
mainly on the upper edge of the allowed region of $\lambda_{H\phi}$. For the
scalar semi-annihilation dominant samples, $\lambda_{H\phi}\gtrsim 10^{-3}$ is
required. The lower limit on $\lambda_{H\phi}$ for the $h\phi$ dominant
channel grows as $m_{\phi}$ increases. So the contribution of the $\phi\phi\to
h\phi$ channel is only important in the range of $[10^{2},10^{3}]$ GeV. Since
the $h\phi$ channel also involves the cubic coupling $\mu$, the $h\phi$
dominant samples have the largest value of $\mu\lambda_{H\phi}$ in the allowed
region. For the secluded channel $\phi\phi\to NN$ and the semi-annihilation
channel $\phi\phi\to N\chi$, lower limits on $y_{N}$ are required, which is
similar to the fermion dark matter. However, due to additional contributions
from the $\phi\phi\to\text{SM}$ and $\phi\phi\to h\phi$ channels, $y_{N}$ can
be as small as $10^{-4}$ when $m_{\phi}$ is above $\mathcal{O}(10)$ GeV, which
is different from the fermion dark matter. In panels (d) and (e) of Figure 6,
we show the corresponding parameters that are involved in the $\phi\phi\to
N\chi$ channel. Approximate lower limits on the factors $\sqrt{y_{\chi}y_{N}}$
and $\mu y_{N}$ exist for $\phi\phi\to N\chi$ dominant samples. To make sure
this semi-annihilation channel is kinematically allowed,
$m_{\chi}/m_{\phi}\lesssim(2m_{\phi}-m_{N})/m_{\phi}\lesssim 2$ should be
satisfied, which is depicted in panel (f). There is also a fake upper limit
$m_{\chi}/m_{\phi}<10$ for the $\phi\phi\to h\phi$ dominant samples, because
it happens to be $m_{\chi}^{\text{max}}\simeq 10^{3}$ GeV and
$m_{\phi}^{\text{min}}\simeq 10^{2}$ GeV.

##  IV Higgs Invisible Decay

In principle, the light sterile neutrino could induce additional Higgs decay
mode $h\to\nu N$ when $m_{N}\lesssim m_{h}$. However, with the seesaw
predicted Yukawa coupling
$y_{\nu}\sim\sqrt{2m_{\nu}m_{N}}/v\sim\mathcal{O}(10^{-6})$, the corresponding
decay width is heavily suppressed. Meanwhile, for a light enough dark matter,
it can contribute to the Higgs invisible decay. The corresponding branching
ratio has been constrained by the ATLAS experiment with ATLAS:2020kdi :

$\displaystyle\mathrm{Br}_{\mathrm{inv}}$
$\displaystyle=\frac{\Gamma_{\mathrm{inv}}}{\Gamma_{\mathrm{inv}}+\Gamma_{\mathrm{SM}}}<0.11,$
(9)

where $\Gamma_{\mathrm{SM}}\simeq 4$ MeV is the standard Higgs width. The
theoretical Higgs invisible decay widths into the dark matter are given by
Escudero:2016ksa :

$\displaystyle\Gamma(h\rightarrow\phi\phi)$ $\displaystyle=$
$\displaystyle\frac{\lambda_{H\phi}^{2}v^{2}}{8\pi
m_{h}}\sqrt{1-\frac{4m_{\phi}^{2}}{m_{h}^{2}}},$ (10)
$\displaystyle\Gamma(h\rightarrow\bar{\chi}\chi)$ $\displaystyle=$
$\displaystyle\frac{m_{h}(\lambda_{H\chi}^{eff})^{2}}{8\pi}\left(1-\frac{4m_{\chi}^{2}}{m_{h}^{2}}\right)^{3/2},$
(11)

where the one-loop effective $h\bar{\chi}\chi$ coupling is

$\displaystyle\lambda_{H\chi}^{eff}$ $\displaystyle=$
$\displaystyle\lambda_{H\phi}\frac{y_{N}^{2}}{16\pi^{2}}\frac{m_{N}}{(m_{\phi}^{2}-m_{N}^{2})^{2}}\left(m_{\phi}^{2}-m_{N}^{2}+m_{N}^{2}\log\frac{m_{N}^{2}}{m_{\phi}^{2}}\right).$
(12)

Figure 7: Branching ratios of invisible Higgs decay induced by fermion (left
panels) and scalar dark matter (right panels). In panels (a)-(d), the purple,
slate blue, and red points are excluded by ATLAS search of Higgs invisible
decay ATLAS:2020kdi , direct detection of LZ LZ:2022ufs , and indirect
detection Campos:2017odj , respectively. The blue points satisfy all current
constraints but are within the reach of the future LZ experiment LZ:2015kxe ,
meanwhile the green points are unconstrained. In panels (e) and (f), we also
classify the samples according to the dominant annihilation channels. Labels
of the samples are the same as in Figure 5 and Figure 6.

Figure 7 shows the theoretical branching ratios of Higgs invisible decay
induced by dark matter and various constraints. For fermion dark matter, the
branching ratio of Higgs invisible decay is below $10^{-10}$, which is because
of the loop suppression of effective Higgs coupling $\lambda_{H\chi}^{eff}$.
The predicted branching ratio is proportional to the coupling
$\lambda_{H\phi}$. For maximum Br${}_{\text{inv}}\simeq 10^{-10}$, the
dominant annihilation channel is found to be $NN$. Meanwhile, for extremely
tiny Br${}_{\text{inv}}\lesssim 10^{-22}$, the $N\chi$ channel is the dominant
one. Because the future HL-LHC could only probe Br${}_{\text{inv}}\gtrsim
0.033$ ATLAS:2018jlh , this negligible tiny branching ratio induced by fermion
dark matter is thus far beyond the reach of direct collider measurement. As
will show later, most samples in the light dark matter region below about 60
GeV are excluded by indirect search.

The scenario for scalar dark matter is quite different, where
Br${}_{\text{inv}}$ could be the dominant decay channel of Higgs when
$\lambda_{H\phi}>0.1$. Under current constraints, the region with
$m_{\phi}\lesssim 53$ GeV and Br${}_{\text{inv}}>0.11$ is excluded by the
ATLAS experiment, which corresponds to $\lambda_{H\phi}\gtrsim 10^{-2}$
disallowed in this region. According to panel (f), most excluded samples in
this region are dominantly annihilation via the $\phi\phi\to$ SM channel. For
$m_{\phi}\gtrsim 10$ GeV, the direct detection experiment as LZ LZ:2022ufs
sets a more stringent constrain than the Higgs invisible decay, where
Br${}_{\text{inv}}\gtrsim 10^{-3}$ and $\lambda_{H\phi}\gtrsim 4\times
10^{-4}$ might be excluded. Such a small branching ratio is also beyond the
reach of future HL-LHC ATLAS:2018jlh . As for scalar dark matter lighter than
10 GeV, the Higgs invisible decay leads to more strict constrain than direct
detection. The annihilation channels of samples in the region which escapes
the limits from Higgs invisible decay and direct detection are $\phi\phi\to
NN$ and $\phi\phi\to N\chi$. Although this region is also tightly constrained
by indirect detection, there are still some samples that satisfy all current
constraints. We have checked that most of these light-mass survived samples
annihilate via $\phi\phi\to N\chi$ with the special requirement
$2m_{\phi}\lesssim m_{N}+m_{\chi}$. So if HL-LHC discovers a relatively large
Br${}_{\text{inv}}$, the dark matter candidate should be a scalar with mass
around a few GeV, meanwhile the sterile neutrino is also at the GeV scale.

##  V Indirect Detection

Figure 8: Exclusion limits of indirect detection experiments. In panels (a)
and (b), the yellow and orange curves represent the bounds of the antiproton-
to-proton flux ratio from AMS-02 AMS:2016oqu  and gamma-rays in the Milky Way
dSphs from FermiLAT Fermi-LAT:2015att  with typical thermal annihilation cross
section $\langle\sigma v\rangle=2.2\times
10^{-26}~{}\text{cm}^{3}\text{s}^{-1}$. In panels (c) to (f), the blue and
purple curves illustrate the limits on annihilation cross section from Fermi-
LAT and H.E.S.S. Campos:2017odj . The yellow line shows the sensitivity of
future CTA experiment for the $W^{+}W^{-}$ annihilation modes CTA:2015yxo .
Other labels are the same as in Figure 7.

The indirect detection experiments aim to search for various types of
particles produced in dark matter annihilation. The differential flux arising
from the annihilation of dark matter is calculated as

$\frac{d\Phi}{dE}=\frac{1}{4\pi}\frac{\langle\sigma
v\rangle}{2m_{\mathrm{DM}}^{2}}\frac{dN}{dE}\cdot\int_{\Delta\Omega}d\Omega\int\rho_{\mathrm{DM}}^{2}(s)ds,$
(13)

where $\rho_{\mathrm{DM}}$ is the dark matter density of the observed object.
The energy spectrum $dN/dE$ describes the distribution of observed particles
from dark annihilation. According to previous studies, current indirect limits
could constrain dark matter mass below about 50 GeV Batell:2017rol . In this
region, the dominant annihilation channels for both fermion and scalar dark
matter are $NN$ and $N\chi$ final states, after imposing the constraint from
Higgs invisible decay. The resulting spectrum $dN/dE$ depends on both the
masses $m_{\text{DM}},m_{N}$ and decay modes of sterile neutrino $N$. With
light $m_{N}<m_{W}$, the three-body decay width via off-shell $W$ and $Z$ can
be estimated as

$\Gamma_{N}\approx\frac{G_{f}^{2}}{192\pi^{3}}|\theta_{\alpha}|^{2}m_{N}^{5},$
(14)

where $\theta_{\alpha}(\alpha=e,\mu,\tau)$ describes the mixing angle between
sterile and active neutrinos for different flavors. For a heavier sterile
neutrino, the two-body decays $N\to W^{\pm}\ell^{\mp},Z\nu,h\nu$ are the
dominant channels. The continuum spectra of muon flavor $N$ is similar to the
electron flavor scenario, while the tau flavor $N$ produces a slightly
stronger gamma-ray spectrum Campos:2017odj . In this paper, we consider an
electron flavor $N$ for a conservative study.

Using the observed spectrum, the exclusion limits on dark matter annihilation
can be derived by performing a likelihood analysis, although the large
astrophysical uncertainties would affect these limits. In Figure 8, we show
the indirect detection constraints from the antiproton observations of AMS-02
and the gamma-rays observations of Fermi-LAT Batell:2017rol . For light
$m_{N}\sim\mathcal{O}$(GeV), the Fermi-LAT experiment could exclude
$m_{\text{DM}}\lesssim 60$ GeV. Meanwhile for heavier $m_{N}$, the AMS-02
result would exclude $m_{\text{DM}}\lesssim 80$ GeV. The combination of these
two limits excludes most samples in the region below $m_{\text{DM}}\lesssim
50\mathrm{GeV}$. It is notable that these two limits in Figure 8 (a) and (b)
are obtained with fixed thermal annihilation cross section $\langle\sigma
v\rangle=2.2\times 10^{-26}~{}\text{cm}^{3}\text{s}^{-1}$. In panel (c)-(f),
we show the theoretically predicted annihilation cross section for indirect
detection, where the cross sections of the semi-annihilation process as
$\chi\chi/\phi\phi\to N\chi$, $\phi\phi\to h\phi$ are multiplied by a factor
of $1/2$. Since the annihilation cross section today can be much smaller than
the thermal target $\langle\sigma v\rangle=2.2\times
10^{-26}~{}\text{cm}^{3}\text{s}^{-1}$, we require that the samples are
excluded by indirect detection when the corresponding cross sections are also
above the Fermi-LAT limit on $\langle\sigma v\rangle$.

In the low mass region below 50 GeV, the two annihilation channels of fermion
dark matter lead to quite different results. For the $\chi\chi\to NN$ dominant
samples, the corresponding annihilation cross sections are at the typical
value $\langle\sigma v\rangle=2.2\times
10^{-26}~{}\text{cm}^{3}\text{s}^{-1}$, so these samples are excluded by
indirect detection. However for the $\chi\chi\to N\chi$ dominant samples, due
to the existence of $s$-channel contribution via the dark scalar $\phi$, the
annihilation cross section could be much smaller than $2.2\times
10^{-26}~{}\text{cm}^{3}\text{s}^{-1}$ when $2m_{\chi}\simeq m_{\phi}$. In
this special scenario, the $\chi\chi\to N\chi$ dominant samples also satisfy
the indirect constraints. For the scalar dark matter, although there are three
annihilation modes in these low mass region, the $\phi\phi\to\text{SM}$
channel is tightly constrained by Higgs invisible decay and direct detection.

Similar to the fermion dark matter, most of the $\phi\phi\to NN$ dominant
samples could be excluded by indirect detection in the low mass region, while
some of the $\phi\phi\to N\chi$ dominant samples are still allowed. Although
there is no on-shell $s$-channel contribution of the $\phi\phi\to N\chi$
channel, we find that the allowed samples satisfy $2m_{\phi}\lesssim
m_{\chi}+m_{N}$ as shown in Figure 9. So these samples fall into the forbidden
region, where the non-relativistic velocity of dark matter today can not
overcome the mass splitting $m_{\chi}+m_{N}-2m_{\phi}$ DAgnolo:2015ujb . The
scalar dark matter annihilates into SM particles also has on-shell $s$-channel
contributions when $2m_{\phi}\simeq m_{h}$. Meanwhile, the scalar semi-
annihilation $\phi\phi\to h\phi$ meets the forbidden condition when
$m_{\phi}\lesssim m_{h}$. The resulting annihilation cross sections of these
two kinds are much smaller than the typical value $\langle\sigma
v\rangle=2.2\times 10^{-26}~{}\text{cm}^{3}\text{s}^{-1}$, so these samples
are hard to probe by indirect detection.

Figure 9: Annihilation cross section for indirect detection. The horizontal
dash-dotted line is the minimum value of the Fermi-LAT limit in Figure 8

For dark matter around the TeV scale, the most stringent constraint is from
H.E.S.S. observation HESS:2011zpk . Currently, no samples could be excluded by
H.E.S.S.. If there is no positive signal at future direct detection
experiments, the scalar semi-annihilation $\phi\phi\to h\phi$ would be
excluded. So in this heavy mass region, we only consider the sterile neutrino
portal annihilation channels $\phi\phi/\chi\chi\to NN$. Because the photon
spectrum from $NN$ final state is similar to $W^{+}W^{-}$ final state, we also
show the future limits from the Cherenkov Telescope Array (CTA) CTA:2015yxo ,
which will cover most samples above 200 GeV.

##  VI Direct Detection

In this model, the dark matter scatters off the atomic nucleus elastically via
the $t$-channel exchange of Higgs boson $h$. For scalar dark matter, this
scattering happens at tree-level, while for fermion dark matter, it is induced
at one-loop level Escudero:2016ksa . The dark matter direct detection
experiments measure the nuclear recoil energy and set constraints on the dark
matter-nucleon scattering cross section. Until now, no concrete signal is
observed by direct detection experiments, such as PandaX-4T PandaX-4T:2021bab
, XENONnT XENON:2023sxq , and LZ LZ:2022ufs . In this paper, we consider the
most stringent limit from LZ and Darkside-50 DarkSide-50:2022qzh  at present
and the future projected limit from LZ LZ:2015kxe . For light dark matter
below about 10 GeV, the Darkside-50 experiment sets the most stringent limit,
which excludes $\sigma_{\rm SI}\gtrsim 10^{-43}~{}\text{cm}^{2}$. For heavier
dark matter, the LZ limit is the most tight one, where the minimum is at
$m_{\text{DM}}=30$ GeV with $\sigma_{\rm SI}=5.9\times
10^{-48}~{}\text{cm}^{2}$.

Within the context of the Higgs-portal effective scenarios Arcadi:2021mag ,
the spin-independent cross section for dark matter collision with nucleon can
be expressed as

$\displaystyle\sigma_{\rm SI}^{\phi n}=\frac{\lambda_{H\phi}^{2}}{\pi
m_{h}^{4}}\frac{m_{n}^{4}f_{n}^{2}}{(m_{\phi}+m_{n})^{2}}\;,$ (15)
$\displaystyle\sigma_{\rm SI}^{\chi n}=\frac{(\lambda_{H\chi}^{eff})^{2}}{\pi
m_{h}^{4}}\frac{m_{n}^{4}m_{\chi}^{2}f_{n}^{2}}{(m_{\chi}+m_{n})^{2}}\;.$ (16)

The nucleon mass is denoted as $m_{n}$, and the parameter $f_{n}\simeq 0.3$
are used to parameterize the Higgs-nucleon interactions Cline:2013gha . The
effective coupling $\lambda_{H\chi}^{eff}$ is calculated in Equation (12).

Figure 10:  The predicted spin-independent cross section and various exclusion
limits. The purple and red lines are the current exclusion limits from
Darkside-50 and LZ respectively. The yellow line is the projected LZ limit.
Other labels are the same as in Figure 7.

The scanned results are shown in Figure 10 for both fermion and scalar dark
matter. Because the scattering cross section for fermion dark matter is
suppressed by the one-loop factor, the predicted values are usually far below
current experimental limits. Within the parameter space scanned in Equation
(8), even the future projection of the LZ experiment could not have a positive
signature when taking into account the limits from indirect detection. We also
find that samples with relatively large scattering cross section $\sigma_{\rm
SI}\gtrsim 10^{-50}~{}\text{cm}^{2}$ are dominant by the semi-annihilation
channel $\chi\chi\to N\chi$. In principle, increasing the maximum values of
related coupling $\lambda_{H\phi}$ and $y_{N}$ to the perturbation limits
could lead to the predicted scattering cross section above current limits
Escudero:2016ksa . These detectable samples are expected to annihilate via
$\chi\chi\to N\chi$ channel in this model.

As for the scalar dark matter, the coupling $\lambda_{H\phi}$ induces quite a
large scattering cross section with correct relic density. The direct
detection experiments, such as LZ, could now exclude samples with
$\lambda_{H\phi}\gtrsim 4\times 10^{-4}$. As already discussed in Section IV,
the Higgs invisible decay has excluded $\lambda_{H\phi}\gtrsim 10^{-2}$ for
$m_{\phi}<m_{h}/2$. It is clear in panel (d) of Figure 10 that the exclusion
limit from Higgs invisible decay is more stringent than the direct limit from
Darkside-50, therefore we do not consider the Darkside-50 limit in this study.
For $m_{\phi}\gtrsim 10$ GeV, the direct limit from LZ is over two orders of
magnitudes tighter than the Higgs invisible limit. Since the indirect
detection has already excluded most samples with $m_{\phi}\lesssim 50$ GeV,
the future projection of LZ is also hard to probe the allowed samples in this
region. Together with panel (f) of Figure 10, the only viable region of
$\phi\phi\to\text{SM}$ dominant samples is the narrow resonance region at
$m_{\phi}\lesssim m_{h}/2$. The predicted cross section of SM dominant samples
can be as small as $\sigma_{\rm SI}\sim 10^{-49}~{}\text{cm}^{2}$, which is
also beyond the reach of future LZ.

For heavier scalar dark matter above 100 GeV, only the direct detection
experiments set the corresponding limit now. The $\phi\phi\to\text{SM}$
dominant samples lead to the largest predicted cross section $\sigma_{\rm
SI}\simeq 10^{-45}~{}\text{cm}^{2}$, which is disfavored by the current LZ
limit when $m_{\phi}<1$ TeV. It is still possible for $\phi\phi\to\text{SM}$
dominant samples with $m_{\phi}$ above 1 TeV Athron:2018ipf , but is out of
the parameter space we scanned in Equation (8). The new contribution of the
semi-annihilation $\phi\phi\to h\phi$ channel would induce a smaller
scattering cross section. Since this semi-annihilation channel also involves
the cubic coupling, the stability and unitarity bounds $\mu<3m_{\phi}$ then
lead to a lower bound for the predicted cross section. For instance, the
minimum predicted cross section is about $3\times 10^{-47}~{}\text{cm}^{2}$
with $m_{\phi}\sim 130$ GeV. Once $\phi\phi\to h\phi$ is kinematically
allowed, this lower bound on $\sigma_{\rm SI}$ increases as $m_{\phi}$ is
larger. Under the current LZ limit, some of the $\phi\phi\to h\phi$ dominant
samples are still allowed. In the future, the projected LZ limit could probe
all the $\phi\phi\to h\phi$ dominant samples. So if there is still no positive
signature at future LZ, the allowed samples will be dominant by $\phi\phi\to
NN$ and $\phi\phi\to N\chi$ channels. We then expect observable signatures at
indirect detection experiments from these two channels for most of the allowed
parameter space.

##  VII Conclusion

Besides generating tiny neutrino mass via the type-I seesaw mechanism, the
electroweak scale sterile neutrino $N$ can also mediate the interaction
between the dark sector and the standard model. Beyond the simplest $Z_{2}$
symmetry, we extend the sterile neutrino portal dark matter with $Z_{3}$
symmetry in this paper. We introduce a scalar singlet $\phi$ and a fermion
singlet $\chi$ to the dark sector. Under the dark $Z_{3}$ symmetry, the dark
sector transforms as $\chi\to e^{i2\pi/3}\chi,\phi\to e^{i2\pi/3}\phi$, while
the standard model particles and the sterile neutrinos transform trivially. In
this paper, we consider WIMP dark matter for both scalar and fermion
scenarios. The $Z_{3}$ symmetry introduces two new interactions, i.e.,
$y_{\chi}\phi\overline{\chi^{c}}\chi$ and $\mu\phi^{3}/2$, which lead to semi-
annihilation channels as $\chi\chi\to N\chi$ and $\phi\phi\to h\phi$.

For the fermion dark matter $\chi$, the annihilation channels are secluded
$\chi\chi\to NN$ and the semi-annihilation $\chi\chi\to N\chi$. Because the
effective $h\bar{\chi}\chi$ coupling is induced at the one-loop level, the
contribution of fermion dark matter to Higgs invisible decay is negligible
tiny. The resulting dark matter-nucleon scattering cross section is also
beyond future experimental reach. Currently, the indirect detection could
exclude most of the samples with $m_{\chi}\lesssim 50$ GeV. In the future, the
CTA experiment is expected to probe the high mass region. However, due to the
$s$-channel contribution of dark scalar $\phi$ to the semi-annihilation
$\chi\chi\to N\chi$, the corresponding annihilation cross section is much
smaller than the usual thermal value $\langle\sigma v\rangle=2\times
10^{-26}~{}\text{cm}^{3}\text{s}^{-1}$ when $2m_{\chi}\simeq m_{\phi}$. In
this special scenario, even the indirect detection can not have a positive
signature.

For the scalar dark matter $\phi$, there are four kinds of annihilation
channels, i.e., the Higgs portal $\phi\phi\to\text{SM}$, the secluded channel
$\phi\phi\to NN$, and the semi-annihilations $\phi\phi\to N\chi,h\phi$. The
direct Higgs portal interaction
$\lambda_{H\phi}(H^{\dagger}H)(\phi^{\dagger}\phi)$ generates observable
signatures from Higgs invisible decay, indirect and direct detection. Under
these constraints, the Higgs portal $\phi\phi\to\text{SM}$ dominant samples
are only viable at the resonance region $m_{\phi}\lesssim m_{h}/2$ for dark
scalar below 1 TeV. The semi-annihilation $\phi\phi\to h\phi$ might be the
dominant one in the range of $[10^{2},10^{3}]$ GeV. This channel predicts a
lower bound on dark matter-nucleon scattering cross section, which can be
fully detected by the future LZ experiment. Meanwhile, the secluded
$\phi\phi\to NN$ and semi-annihilation $\phi\phi\to N\chi$ channels can easily
satisfy current bounds and are promising at indirect detection experiments.
Although the light mass region is tightly constrained, we find that when the
forbidden relation $2m_{\phi}\lesssim m_{\chi}+m_{N}$ is satisfied, the semi-
annihilation channel $\phi\phi\to N\chi$ also has a suppressed annihilation
cross section for indirect detection.

Compared with the $Z_{2}$ symmetric model, the new interactions introduced in
the $Z_{3}$ symmetric model enlarge the viable parameter space. For instance,
light dark matter below about 50 GeV is entirely excluded by indirect
detection in the $Z_{2}$ symmetric model. However, in the $Z_{3}$ symmetric
model, the semi-annihilation $\chi\chi\to N\chi$ and $\phi\phi\to N\chi$ are
still allowed under certain circumstances. Therefore, once light dark matter
is discovered, the $Z_{3}$ symmetric model would be preferred.

## Acknowledgments

This work is supported by the National Natural Science Foundation of China
under Grant No. 11805081 and 11635009, Natural Science Foundation of Shandong
Province under Grant No. ZR2019QA021 and ZR2022MA056, the Open Project of
Guangxi Key Laboratory of Nuclear Physics and Nuclear Technology under Grant
No. NLK2021-07.

## References

  * (1) 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]].
  * (2) C. P. Burgess, M. Pospelov and T. ter Veldhuis, Nucl. Phys. B 619, 709-728 (2001) [arXiv:hep-ph/0011335 [hep-ph]].
  * (3) J. McDonald, Phys. Rev. D 50, 3637-3649 (1994) [arXiv:hep-ph/0702143 [hep-ph]].
  * (4) J. M. Cline, K. Kainulainen, P. Scott and C. Weniger, Phys. Rev. D 88, 055025 (2013) [erratum: Phys. Rev. D 92, no.3, 039906 (2015)] [arXiv:1306.4710 [hep-ph]].
  * (5) A. Alves, S. Profumo and F. S. Queiroz, JHEP 04, 063 (2014) [arXiv:1312.5281 [hep-ph]].
  * (6) A. Alves, A. Berlin, S. Profumo and F. S. Queiroz, Phys. Rev. D 92, no.8, 083004 (2015) [arXiv:1501.03490 [hep-ph]].
  * (7) W. Rodejohann and C. E. Yaguna, JCAP 12, 032 (2015) [arXiv:1509.04036 [hep-ph]].
  * (8) M. Klasen, F. Lyonnet and F. S. Queiroz, Eur. Phys. J. C 77, no.5, 348 (2017) [arXiv:1607.06468 [hep-ph]].
  * (9) P. Minkowski, Phys. Lett. B 67, 421-428 (1977)
  * (10) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980)
  * (11) S. Dodelson and L. M. Widrow, Phys. Rev. Lett. 72 (1994), 17-20 [arXiv:hep-ph/9303287 [hep-ph]].
  * (12) A. Boyarsky, M. Drewes, T. Lasserre, S. Mertens and O. Ruchayskiy, Prog. Part. Nucl. Phys. 104, 1-45 (2019) [arXiv:1807.07938 [hep-ph]].
  * (13) B. M. Roach, K. C. Y. Ng, K. Perez, J. F. Beacom, S. Horiuchi, R. Krivonos and D. R. Wik, Phys. Rev. D 101, no.10, 103011 (2020) [arXiv:1908.09037 [astro-ph.HE]].
  * (14) A. Boyarsky, J. Franse, D. Iakubovskyi and O. Ruchayskiy, Phys. Rev. Lett. 115, 161301 (2015) [arXiv:1408.2503 [astro-ph.CO]].
  * (15) B. M. Roach, S. Rossland, K. C. Y. Ng, K. Perez, J. F. Beacom, B. W. Grefenstette, S. Horiuchi, R. Krivonos and D. R. Wik, Phys. Rev. D 107, no.2, 023009 (2023) [arXiv:2207.04572 [astro-ph.HE]].
  * (16) E. Ma, Phys. Rev. D 73 (2006), 077301 [arXiv:hep-ph/0601225 [hep-ph]].
  * (17) E. Ma, Phys. Lett. B 662 (2008), 49-52 [arXiv:0708.3371 [hep-ph]].
  * (18) R. Ding, Z. L. Han, Y. Liao and W. P. Xie, JHEP 05, 030 (2016) [arXiv:1601.06355 [hep-ph]].
  * (19) J. Liu, Z. L. Han, Y. Jin and H. Li, JHEP 12, 057 (2022) [arXiv:2207.07382 [hep-ph]].
  * (20) A. Falkowski, J. Juknevich and J. Shelton, [arXiv:0908.1790 [hep-ph]].
  * (21) A. Falkowski, J. T. Ruderman and T. Volansky, JHEP 05 (2011), 106 [arXiv:1101.4936 [hep-ph]].
  * (22) V. Gonzalez Macias and J. Wudka, JHEP 07, 161 (2015) [arXiv:1506.03825 [hep-ph]].
  * (23) V. González-Macías, J. I. Illana and J. Wudka, JHEP 05, 171 (2016) [arXiv:1601.05051 [hep-ph]].
  * (24) M. Blennow, E. Fernandez-Martinez, A. Olivares-Del Campo, S. Pascoli, S. Rosauro-Alcaraz and A. V. Titov, Eur. Phys. J. C 79 (2019) no.7, 555 [arXiv:1903.00006 [hep-ph]].
  * (25) P. Ballett, M. Hostert and S. Pascoli, Phys. Rev. D 101 (2020) no.11, 115025 [arXiv:1903.07589 [hep-ph]].
  * (26) A. Liu, Z. L. Han, Y. Jin and F. X. Yang, Phys. Rev. D 101, no.9, 095005 (2020) [arXiv:2001.04085 [hep-ph]].
  * (27) L. Coito, C. Faubel, J. Herrero-García, A. Santamaria and A. Titov, JHEP 08, 085 (2022) [arXiv:2203.01946 [hep-ph]].
  * (28) A. Liu, F. L. Shao, Z. L. Han, Y. Jin and H. Li, Eur. Phys. J. C 83, no.5, 423 (2023) [arXiv:2205.11846 [hep-ph]].
  * (29) M. Pospelov, A. Ritz and M. B. Voloshin, Phys. Lett. B 662 (2008), 53-61 [arXiv:0711.4866 [hep-ph]].
  * (30) Y. L. Tang and S. h. Zhu, JHEP 03 (2016), 043 [arXiv:1512.02899 [hep-ph]].
  * (31) M. D. Campos, F. S. Queiroz, C. E. Yaguna and C. Weniger, JCAP 07 (2017), 016 [arXiv:1702.06145 [hep-ph]].
  * (32) B. Batell, T. Han and B. Shams Es Haghi, Phys. Rev. D 97 (2018) no.9, 095020 [arXiv:1704.08708 [hep-ph]].
  * (33) B. Batell, T. Han, D. McKeen and B. Shams Es Haghi, Phys. Rev. D 97 (2018) no.7, 075016 [arXiv:1709.07001 [hep-ph]].
  * (34) M. G. Folgado, G. A. Gómez-Vargas, N. Rius and R. Ruiz De Austri, JCAP 08 (2018), 002 [arXiv:1803.08934 [hep-ph]].
  * (35) A. Falkowski, E. Kuflik, N. Levi and T. Volansky, Phys. Rev. D 99, no.1, 015022 (2019) [arXiv:1712.07652 [hep-ph]].
  * (36) M. Becker, Eur. Phys. J. C 79, no.7, 611 (2019) [arXiv:1806.08579 [hep-ph]].
  * (37) M. Chianese and S. F. King, JCAP 09, 027 (2018) [arXiv:1806.10606 [hep-ph]].
  * (38) P. Bandyopadhyay, E. J. Chun and R. Mandal, JCAP 08, 019 (2020) [arXiv:2005.13933 [hep-ph]].
  * (39) Y. Cheng and W. Liao, Phys. Lett. B 815, 136118 (2021) [arXiv:2012.01875 [hep-ph]].
  * (40) Z. F. Chang, Z. X. Chen, J. S. Xu and Z. L. Han, JCAP 06, 006 (2021) [arXiv:2104.02364 [hep-ph]].
  * (41) B. Barman, P. S. Bhupal Dev and A. Ghoshal, [arXiv:2210.07739 [hep-ph]].
  * (42) A. Liu, F. L. Shao, Z. L. Han, Y. Jin and H. Li, [arXiv:2212.10043 [hep-ph]].
  * (43) M. Escudero, N. Rius and V. Sanz, Eur. Phys. J. C 77, no.6, 397 (2017) [arXiv:1607.02373 [hep-ph]].
  * (44) M. Escudero, N. Rius and V. Sanz, JHEP 02, 045 (2017) [arXiv:1606.01258 [hep-ph]].
  * (45) B. Batell, Phys. Rev. D 83, 035006 (2011) [arXiv:1007.0045 [hep-ph]].
  * (46) M. Hirsch, S. Morisi, E. Peinado and J. W. F. Valle, Phys. Rev. D 82, 116003 (2010) [arXiv:1007.0871 [hep-ph]].
  * (47) J. L. Diaz-Cruz and E. Ma, Phys. Lett. B 695, 264-267 (2011) [arXiv:1007.2631 [hep-ph]].
  * (48) G. Belanger, K. Kannike, A. Pukhov and M. Raidal, JCAP 04, 010 (2012) [arXiv:1202.2962 [hep-ph]].
  * (49) G. Belanger, K. Kannike, A. Pukhov and M. Raidal, JCAP 01 (2013), 022 [arXiv:1211.1014 [hep-ph]].
  * (50) P. Bandyopadhyay, D. Choudhury and D. Sachdeva, Phys. Rev. D 107, no.1, 015020 (2023) [arXiv:2206.05811 [hep-ph]].
  * (51) D. K. Ghosh, P. Ghosh and S. Jeesun, [arXiv:2301.13754 [hep-ph]].
  * (52) A. Liu, F. L. Shao, Z. L. Han, Y. Jin and H. Li, in preparation
  * (53) A. M. Abdullahi, P. B. Alzas, B. Batell, J. Beacham, A. Boyarsky, S. Carbajal, A. Chatterjee, J. I. Crespo-Anadon, F. F. Deppisch and A. De Roeck, et al. J. Phys. G 50, no.2, 020501 (2023) [arXiv:2203.08039 [hep-ph]].
  * (54) A. Hektor, A. Hryczuk and K. Kannike, JHEP 03, 204 (2019) [arXiv:1901.08074 [hep-ph]].
  * (55) F. D’Eramo and J. Thaler, JHEP 06 (2010), 109 [arXiv:1003.5912 [hep-ph]].
  * (56) Y. Cai and A. P. Spray, JHEP 01, 087 (2016) [arXiv:1509.08481 [hep-ph]].
  * (57) J. Guo, Z. Kang and P. Zhang, Phys. Lett. B 826, 136885 (2022) [arXiv:2108.12964 [hep-ph]].
  * (58) S. P. Li and X. J. Xu, [arXiv:2212.09109 [hep-ph]].
  * (59) G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov and B. Zaldivar, Comput. Phys. Commun. 231, 173-186 (2018) [arXiv:1801.03509 [hep-ph]].
  * (60) N. Aghanim et al. [Planck], Astron. Astrophys. 641 (2020), A6 [erratum: Astron. Astrophys. 652 (2021), C4] [arXiv:1807.06209 [astro-ph.CO]].
  * (61) O. Ruchayskiy and A. Ivashko, JCAP 10, 014 (2012) [arXiv:1202.2841 [hep-ph]].
  * (62) [ATLAS], ATLAS-CONF-2020-052.
  * (63) [ATLAS], ATL-PHYS-PUB-2018-054.
  * (64) J. Aalbers et al. [LZ], [arXiv:2207.03764 [hep-ex]].
  * (65) D. S. Akerib et al. [LZ], [arXiv:1509.02910 [physics.ins-det]].
  * (66) M. Aguilar et al. [AMS], Phys. Rev. Lett. 117, no.9, 091103 (2016)
  * (67) M. Ackermann et al. [Fermi-LAT], Phys. Rev. Lett. 115, no.23, 231301 (2015) [arXiv:1503.02641 [astro-ph.HE]].
  * (68) J. Carr et al. [CTA], PoS ICRC2015 (2016), 1203 [arXiv:1508.06128 [astro-ph.HE]].
  * (69) R. T. D’Agnolo and J. T. Ruderman, Phys. Rev. Lett. 115, no.6, 061301 (2015) [arXiv:1505.07107 [hep-ph]].
  * (70) A. Abramowski et al. [H.E.S.S.], Phys. Rev. Lett. 106, 161301 (2011) [arXiv:1103.3266 [astro-ph.HE]].
  * (71) Y. Meng et al. [PandaX-4T], Phys. Rev. Lett. 127, no.26, 261802 (2021) [arXiv:2107.13438 [hep-ex]].
  * (72) E. Aprile et al. [XENON], [arXiv:2303.14729 [hep-ex]].
  * (73) P. Agnes et al. [DarkSide-50], Phys. Rev. D 107, no.6, 063001 (2023) [arXiv:2207.11966 [hep-ex]].
  * (74) G. Arcadi, A. Djouadi and M. Kado, Eur. Phys. J. C 81, no.7, 653 (2021) [arXiv:2101.02507 [hep-ph]].
  * (75) P. Athron, J. M. Cornell, F. Kahlhoefer, J. Mckay, P. Scott and S. Wild, Eur. Phys. J. C 78, no.10, 830 (2018) [arXiv:1806.11281 [hep-ph]].