#  CP violation in $B_{s}$ decays mediated by two nearly-degenerate Majorana
neutrinos111Phys.Rev.D 105 (2022) 073001,
https://link.aps.org/doi/10.1103/PhysRevD.105.073001

Guangshuai Zhang1 Bo-Qiang Ma1,2,3<EMAIL_ADDRESS>1 School of Physics and
State Key Laboratory of Nuclear Physics and Technology, Peking University,
Beijing 100871, China
2Center for High Energy Physics, Peking University, Beijing 100871, China
3Collaborative Innovation Center of Quantum Matter, Beijing, China

###### Abstract

We study the CP violation in the lepton-number-violating process
$B_{s}^{0}\rightarrow D_{s}^{-}\mu^{+}\mu^{+}\pi^{-}$ and its CP-conjugate
process $\bar{B_{s}^{0}}\rightarrow D_{s}^{+}\mu^{-}\mu^{-}\pi^{+}$ that are
induced by two GeV-scale Majorana neutrinos with nearly-degenerate masses. Our
result shows that the size of the CP violation could become considerable if
the mass difference between the two Majorana neutrinos is around the decay
width of the neutrinos. We perform experimental analysis on the CP-violating
processes at LHCb within its upgrade II. The analysis shows that under current
constraint on the heavy-light neutrino mixing parameter, it is possible that
such CP violation is observed with the LHCb experimental ability. We also give
the upper bound on the heavy-light mixing parameter $|U_{\mu N}|^{2}$ under
the assumption that no positive signal of the process is observed. The result
shows that such modes can give complementary constraint compared with previous
experiments such as NuTeV, BEBC, etc. in the mass region 1 GeV $<m_{N}<$ 3
GeV.

Majorana neutrino; CP violation; $B_{s}$ decay; lepton-number-violating
process

##  I Introduction

The neutrino oscillation experiments have confirmed that at least two of the
three generations of neutrinos have nonzero masses Super-Kamiokande:1998kpq ;
Super-Kamiokande:2010orq . The relation between the flavor eigenstates and
mass eigenstates of three generations of neutrinos is described by the
Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix. Thus, two questions about the
nature of the neutrino arise. The first is the origin of the neutrino mass. As
a well-known approach for the generation of neutrino mass, the seesaw
mechanism offers a natural explanation for the tininess of the neutrino mass,
for which the introduction of one or more generations of heavy Majorana
neutrinos is necessary Minkowski:1977sc ; GellMann:1980vs ; Mohapatra:1979ia ;
Yanagida:1980xy . Studies on the cosmological effect of the heavy neutrinos
have shown that these heavy neutrinos can generate the observed baryon
asymmetry of the Universe through leptogenesis Canetti:2012kh , and such heavy
neutrinos can also serve as natural candidates for dark matter Dodelson:1993je
; Shi:1998km ; Asaka:2005pn ; Canetti:2012vf .

Another question is whether a neutrino is a Majorana particle, i.e., its
antiparticle is identical to itself, or not. If the neutrino is a Dirac
particle, then the lepton number is conserved ($\Delta L=0$). While if it is a
Majorana particle, then the conservation of the lepton number can be violated
by 2 units ($\Delta L=2$). Thus, the most important approach to establish the
Majorana nature of the neutrino is to search for the lepton-number-violating
(LNV) processes. Up to now, the neutrinoless double beta ($0\nu\beta\beta$)
decay is the most explored LNV channel DellOro:2016tmg ; Dolinski:2019nrj ;
Engel:2016xgb . Another way is to look for LNV processes that are induced by
Majorana neutrino in the hardron or $\tau$ lepton decays. In literature, the
LNV processes in the decays of mesons ($K,D,D_{s},B,B_{c}$) Abad:1984gh ;
Littenberg:1991ek ; Littenberg:2000fg ; Ali:2001gsa ; Ivanov:2004ch ;
Dib:2000wm ; Atre:2005eb ; Atre:2009rg ; Helo:2010cw ; Cvetic:2010rw ;
Cvetic:2016fbv ; Cvetic:2017vwl ; Cvetic:2019shl ; Milanes:2016rzr ; Mejia-
Guisao:2017gqp ; Milanes:2018aku ; Asaka:2016rwd ; Chun:2019nwi ;
Quintero:2011yh , baryons ($\Sigma^{-},\Xi^{-},\Lambda_{b}$) Barbero:2002wm ;
Barbero:2007zm ; Barbero:2013fc ; Mejia-Guisao:2017nzx ; Zhang:2021wjj , and
$\tau$ lepton Castro:2012gi ; Gribanov:2001vv ; Cvetic:2002jy ; Atre:2005eb
were extensively investigated. If the mass of the hypothetical Majorana
neutrino lies in the range from hundreds of keV to several GeV, the decay
widths of the corresponding LNV hadron/$\tau$ decays can be resonantly
enhanced due to the on shellness of the intermediate Majorana neutrino. Thus,
it is possible for these LNV decays to be observed by current or future
hardron collider experiments. On the other hand, the nonobservation of these
LNV hadron or $\tau$ decays can set strong constraints on the heavy-light
neutrino mixing parameters in the resonant range.

The neutrino oscillation experiments also showed that the $\theta_{13}$ angle
in the PMNS matrix has a nonzero value T2K:2011ypd ; DayaBay:2012fng ;
RENO:2012mkc ; thus, the CP violation in the neutrino sector is still
possible. Moreover, the introduction of one or more generations of heavy
neutrinos brings in more CP-violating phases in the extended PMNS matrix,
which describes the mixing between the flavor eigenstates and mass eigenstates
of three normal neutrinos as well as the hypothetical sterile ones. Suppose
there only exists one generation of sterile neutrino, the CP-violating phases
cause no observable effect in the sterile-neutrino-induced LNV processes.
However, if there exist two generations of GeV-scale sterile neutrinos that
have nearly-degenerate masses, the CP-violating phases in the extended-PMNS
matrix could cause observable CP violation in the sterile-neutrino-induced LNV
meson decay. The idea was first pointed out in Ref. Cvetic:2013eza . The CP
violation in such Majorana-neutrino-induced LNV processes arises through two
different mechanisms. The first is the interference between the amplitudes
contributed by the two generations of sterile neutrinos. And the second is the
neutrino oscillation during the propagation of the two on shell sterile
neutrinos, which was first investigated in Ref. Cvetic:2015ura  in LNV $B$
meson decays. CP-violating LNV processes induced by similar mechanisms in the
decays of other mesons Cvetic:2014nla ; Dib:2014pga ; Cvetic:2015naa ;
Cvetic:2015ura ; Cvetic:2020lyh ; Zhang:2020hwj  and $\tau$ leptons Zamora-
Saa:2016ito ; Tapia:2019coy  as well as $W$ boson Najafi:2020dkp  were also
studied in detail in the literature. Besides, the existence of two nearly-
degenerate Majorana neutrinos together with the CP violating phases can lead
to a significant difference between the decay widths of the LNV meson decays
and those of the lepton-number-conserving ones Abada:2019bac ; Das:2017hmg ,
which is contrary to the common hypothesis. A remarkable conclusion about the
CP violation in LNV meson decays due to intermediate sterile neutrino
interference is that the relative size of the CP violation is independent of
the neutrino mass while the parameter $\Delta m_{N}/\Gamma_{N}$ remains
unchanged, where $\Delta m_{N}$ is the masses difference and $\Gamma_{N}$ the
decay width of the neutrino. It is still unclear whether the conclusion holds
true for four-body decays. Thus in this paper, we would apply the mechanism to
four-body decays of $B_{s}$ mesons, of which the Feynman diagram is shown in
Fig. 1. To our knowledge, no previous studies have explored the Majorana-
neutrino-induced LNV decays of $B_{s}$ meson except Ref. Mejia-Guisao:2017gqp
. Neither has any experimental research about the process been done. The mass
difference between $B_{s}$ and $D_{s}$ reaches 3.4 GeV, which may extend the
constraining mass region provided by previous channels such as three-body $D$
and $B$ meson decays. On the other hand, the branching fraction of
$B_{s}\rightarrow D_{s}l\nu$ in $B_{s}$ decay reaches $8.1\%$, which means the
branching fraction of the related LNV processes may be considerable since the
latter is proportional to the former in case that $m_{N}$ equals zero.

The motivation for the existence of two heavy neutrinos with nearly-degenerate
masses comes from the $\nu$MSM model Asaka:2005pn ; Gorbunov:2007ak , which
proposes the existence of two generations of Majorana neutrinos with almost
degenerate masses between 100 MeV and a few GeV as well as a light Majorana
neutrino of mass $\sim$ 10 keV. The $\nu$MSM allows one to explain
simultaneously neutrino oscillations, dark matter and baryon asymmetry of the
Universe Gorbunov:2007ak .

In this work we study the CP violation between the Majorana-neutrino-mediated
LNV process $B_{s}^{0}\rightarrow D_{s}^{-}\mu^{+}\mu^{+}\pi^{-}$ and its CP-
conjugate process $\bar{B_{s}^{0}}\rightarrow D_{s}^{+}\mu^{-}\mu^{-}\pi^{+}$.
These processes are induced by two Majorana neutrinos that have nearly
degenerate but not equal masses. We focus on the neutrino mass region between
0.5 and 3.5 GeV where the resonant enhancement appears. We deal with the CP
violation that arises from the interference between the amplitudes contributed
by the two Majorana neutrinos. Moreover, we investigate the possibility of the
detection of CP asymmetries in such decays during the LHCb upgrade II
LHCb:2018roe .

The structure of this paper is as follows. In Sec. II, we introduce the
formalism we used for calculation and give the expression for the size of CP
violation. In Sec. III, we perform an experimental analysis on the processes
under the experimental ability of LHCb during its upgrade II and give the
upper bounds on $|U_{\mu N}|^{2}$ under the assumption that such modes were
not observed in experiments. Sec. IV gives the summary and the conclusion of
our work.

##  II Formalism

###  II.1 Decay widths for the two processes

We define the light neutrino flavor eigenstates as

$\nu_{l}=\sum_{i=1}^{3}U_{l\nu_{i}}\nu_{i}+U_{lN_{1}}N_{1}+U_{lN_{2}}N_{2}.$
(1)

Here, $\nu_{i}~{}(i=1,2,3)$ and $N_{j}~{}(j=1,2)$ represent the mass
eigenstates of a light neutrino and heavy neutrino, respectively.
$U_{lN_{j}}~{}(j=1,2)$ is the heavy-light neutrino mixing elements of the
extended PMNS matrix (between $l$ lepton and the $j$th heavy neutrino). We
parametrize $U_{lN_{i}}$ as

$U_{lN_{i}}=|U_{lN_{i}}|e^{i\phi_{lN_{i}}}.$ (2)

Here, $\phi_{lN_{i}}$ is the CP-odd phase angle. We also require that the
masses of $N_{1}$ and $N_{2}$ satisfy that

$\displaystyle m_{\mu}+m_{\pi}\leq m_{N_{1}}$ $\displaystyle\leq
m_{B_{s}}-m_{D_{s}}-m_{\mu},$ (3) $\displaystyle~{}m_{\mu}+m_{\pi}\leq
m_{N_{2}}$ $\displaystyle\leq m_{B_{s}}-m_{D_{s}}-m_{\mu},$

so that the intermediate neutrinos are nearly on shell and thus, the resonant
enhancement appears. Another assumption we make is that the masses of the two
neutrinos are nearly degenerate (suppose $m_{N_{1}}<m_{N_{2}}$),

$\Delta m_{N}=m_{N_{2}}-m_{N_{1}}\ll m_{N_{i}}~{}(i=1,2).$ (4)

Figure 1: Feynman diagrams for $B_{s}^{0}\rightarrow
D_{s}^{-}\mu^{+}\mu^{+}\pi^{-}$ and $\bar{B}_{s}^{0}\rightarrow
D_{s}^{+}\mu^{-}\mu^{-}\pi^{+}$

The Feynman diagrams for the two processes are shown in Fig. 1.222There exists
a “crossed” diagram where $\mu_{1}$ and $\mu_{2}$ are exchanged. Since when
the lifetime of the intermediate heavy neutrino is long enough, the two
leptons appear at displaced vertices and the corresponding processes of
“direct” channel and “cross” channel can be distinguished from each other
Dib:2014iga . Thus, there is no interference between them and the
corresponding decay width can be added directly. This brings in a factor 2 in
the final result, which cancels out with the factor 1/2! due to the
indistinguishability between the two muons .. We denote the amplitude of
$B_{s}^{0}\rightarrow D_{s}^{-}\mu^{+}\mu^{+}\pi^{-}$ as $\mathcal{M}^{+}$ and
that of its CP-conjugate process as $\mathcal{M}^{-}$. The amplitudes can be
written explicitly from the Feynman diagram

$\displaystyle
i\mathcal{M}^{+}=G_{F}^{2}V_{cb}^{*}V_{ud}^{*}f_{\pi}^{*}\sum_{i=1}^{2}U_{\mu
N_{i}}^{*2}m_{N_{i}}P_{N_{i}}\bar{u}(\mu_{2})\not{p}_{\pi}\gamma_{\mu}(1-\gamma_{5})v(\mu_{1})\langle
D_{s}^{-}|c(0)\gamma^{\mu}b(0)|B_{s}^{0}\rangle,$ (5) $\displaystyle
i\mathcal{M}^{-}=G_{F}^{2}V_{cb}V_{ud}f_{\pi}\sum_{i=1}^{2}U_{\mu
N_{i}}^{2}m_{N_{i}}P_{N_{i}}\bar{u}(\mu_{2})\not{p}_{\pi}\gamma_{\nu}(1+\gamma_{5})v(\mu_{1})\langle
D_{s}^{+}|\bar{c}(0)\gamma^{\nu}\bar{b}(0)|\bar{B}_{s}^{0}\rangle.$

Here, $G_{F}=1.1664\times 10^{-5}\mathrm{GeV}^{-2}$ is the Fermi coupling
constant, and $f_{\pi}$ = 0.1304 GeV is the pion decay constant Zyla:2020zbs .
$V_{cb}$ and $V_{ud}$ are the charm-bottom and upper-down Cabibbo-Kobayashi-
Maskawa (CKM) matrix element. In this work we take that $|V_{ud}|=0.974$ and
$|V_{cb}|=4.012\times 10^{-2}$ Zyla:2020zbs . $m_{N_{i}}$ is the mass of the
$i$th heavy neutrino. $\bar{u}(\mu_{2})$ and $v(\mu_{1})$ are the spinors of
the two muons. The propagator $P_{N_{i}}$ is defined as

$P_{N_{i}}=\frac{1}{(p_{B_{s}}-p_{D_{s}}-p_{\mu_{1}})^{2}-m_{N_{i}}^{2}+i\Gamma_{N_{i}}M_{N_{i}}}.$
(6)

Also, $\langle D_{s}^{-}|c(0)\gamma_{\mu}b(0)|B_{s}^{0}\rangle=\langle
D_{s}^{+}|\bar{c}(0)\gamma_{\mu}\bar{b}(0)|\bar{B}_{s}^{0}\rangle$ is the
$B_{s}-D_{s}$ transition matrix element that can be parametrized as
Monahan:2017uby

$\langle
D_{s}^{-}|c(0)\gamma_{\mu}b(0)|B_{s}^{0}\rangle=f_{0}(t_{1})\frac{m_{B_{s}}^{2}-m_{D_{s}}^{2}}{t_{1}}q^{\mu}+f_{+}(t_{1})\Big{[}p_{B_{s}}^{\mu}+p_{D_{s}}^{\mu}-\frac{m_{B_{s}}^{2}-m_{D_{s}}^{2}}{t_{1}}q^{\mu}\Big{]},$
(7)

where $q^{\mu}=p_{B_{s}}^{\mu}-p_{D_{s}}^{\mu}$ is the transferred momentum
and $t_{1}=q^{2}$. In this work, we use the numerical result of the form
factors $f_{0}(t_{1})$ and $f_{+}(t_{1})$ from Ref. Monahan:2017uby  which is
calculated by lattice QCD. Equation (5) shows that two CP-odd factors appear
in the amplitudes that are nontrivial. The first is the CP-odd phases of the
heavy-light mixing parameter $U_{\mu N_{i}}$ and the second comes from the
weak interaction vertex $\gamma^{\mu}(1-\gamma_{5})$. In case there exists
only one generation of heavy neutrino, such a CP-odd term would not result in
a CP violation in the observable quantity since after being squared, the
difference between the two amplitudes vanishes. While in case there are more
than one generation of heavy neutrino, a CP violation would arise in the decay
width. For simplicity, we write the amplitudes as

$\displaystyle\mathcal{M^{+}}$ $\displaystyle=U_{\mu
N_{1}}^{*2}\bar{\mathcal{M}}^{+}_{1}+U_{\mu
N_{2}}^{*2}\bar{\mathcal{M}}^{+}_{2},$ (8) $\displaystyle\mathcal{M^{-}}$
$\displaystyle=U_{\mu N_{1}}^{2}\bar{\mathcal{M}}^{-}_{1}+U_{\mu
N_{2}}^{2}\bar{\mathcal{M}}^{-}_{2}.$

Here $\bar{\mathcal{M}}^{\pm}_{i}$ represents the canonical amplitude
contributed by the $i$th heavy neutrino

$\displaystyle\bar{\mathcal{M}}^{+}_{i}$
$\displaystyle=G_{F}^{2}V_{cb}^{*}V_{ud}^{*}f_{\pi}^{*}m_{N_{i}}P_{N_{i}}\bar{u}(\mu_{2})\not{p}_{\pi}\gamma_{\mu}(1-\gamma_{5})v(\mu_{1})\langle
D_{s}^{-}|c(0)\gamma^{\mu}b(0)|B_{s}^{0}\rangle,$ (9)
$\displaystyle\bar{\mathcal{M}}^{-}_{i}$
$\displaystyle=G_{F}^{2}V_{cb}V_{ud}f_{\pi}m_{N_{i}}P_{N_{i}}\bar{u}(\mu_{2})\not{p}_{\pi}\gamma_{\mu}(1+\gamma_{5})v(\mu_{1})\langle
D_{s}^{-}|c(0)\gamma^{\mu}b(0)|B_{s}^{0}\rangle.$

It is straightforward to prove that

$\bar{\mathcal{M}}^{+}_{i}\bar{\mathcal{M}}^{+*}_{j}=\bar{\mathcal{M}}^{-}_{i}\bar{\mathcal{M}}^{-*}_{j}.$
(10)

We then define the squared amplitude matrix $M$ as

$M_{ij}=\bar{\mathcal{M}}^{+}_{i}\bar{\mathcal{M}}^{+*}_{j}=\bar{\mathcal{M}}^{-}_{i}\bar{\mathcal{M}}^{-*}_{j}.$
(11)

The explicit form of $M$ is shown in Appendix B. From Equation (49) in
Appendix B we can verify that

$M_{12}=M_{21}^{*}.$ (12)

Then the decay width for the two processes can be written as

$\displaystyle\Gamma^{+}$ $\displaystyle=\int d\Phi_{4}\Big{[}|U_{\mu
N_{1}}|^{4}M_{11}+|U_{\mu N_{2}}|^{4}M_{22}+U_{\mu N_{1}}^{*2}U_{\mu
N_{2}}^{2}M_{12}+U_{\mu N_{1}}^{2}U_{\mu N_{2}}^{*2}M_{21}\Big{]},$ (13)
$\displaystyle\Gamma^{-}$ $\displaystyle=\int d\Phi_{4}\Big{[}|U_{\mu
N_{1}}|^{4}M_{11}+|U_{\mu N_{2}}|^{4}M_{22}+U_{\mu N_{1}}^{2}U_{\mu
N_{2}}^{*2}M_{12}+U_{\mu N_{1}}^{*2}U_{\mu N_{2}}^{2}M_{21}\Big{]},$

where $d\Phi_{4}$ is the four-body phase space

$d\Phi_{4}=\frac{1}{2m_{B_{s}}}\frac{1}{(2\pi)^{8}}\frac{d^{3}\bm{p}_{D_{s}}}{(2E_{D_{s}})}\frac{d^{3}\bm{p}_{1}}{(2E_{1})}\frac{d^{3}\bm{p}_{2}}{(2E_{2})}\frac{d^{3}\bm{p}_{\pi}}{(2E_{\pi})}\delta^{4}(p_{B_{s}}-p_{D_{s}}-p_{\mu_{1}}-p_{\mu_{2}}-p_{\pi}).$
(14)

The explicit form of $d\Phi_{4}$ and the reference frame we use to define the
kinematic variables for numerical calculation are given in Appendix A. Note
that the factors $|P_{N_{i}}|^{2}$ appears in Equation (49), which can be
approximated as

$|P_{N_{i}}|^{2}=|\frac{1}{p_{N}^{2}-m_{N_{i}}^{2}+i\Gamma_{N_{i}}m_{N_{i}}}|^{2}\approx\frac{\pi}{m_{N_{i}}\Gamma_{N_{i}}}\delta(p_{N}^{2}-m_{N_{i}}^{2})~{}(i=1,2),$
(15)

when it is satisfied that $\Gamma_{N_{i}}\ll m_{N_{i}}~{}(i=1,2)$. Thus with
the narrow width approximation, $\int d\Phi_{4}M_{ii}$ can be simplified as

$\int d\Phi_{4}M_{ii}=\bar{\Gamma}(B_{s}\rightarrow
D_{s}\mu_{1}N_{i})\frac{\bar{\Gamma}(N_{i}\rightarrow\mu_{2}\pi)}{\Gamma_{N_{i}}}.$
(16)

Here $\bar{\Gamma}(B_{s}\rightarrow D_{s}\mu_{1}N_{i})$ and
$\bar{\Gamma}(N_{i}\rightarrow\mu_{2}\pi)$ are the canonical decay widths for
the subprocesses $B_{s}\rightarrow D_{s}\mu_{1}N_{i}$ and
$N_{i}\rightarrow\mu_{2}\pi$, respectively,

$\bar{\Gamma}(B_{s}\rightarrow
D_{s}\mu_{1}N_{i})=\frac{\Gamma(B_{s}\rightarrow D_{s}\mu_{1}N_{i})}{|U_{\mu
N_{i}}|^{2}},\quad\bar{\Gamma}(N_{i}\rightarrow\mu_{2}\pi)=\frac{\Gamma(N_{i}\rightarrow\mu_{2}\pi)}{|U_{\mu
N_{i}}|^{2}}.$ (17)

The last part that is yet unknown is the decay width for the heavy neutrinos
$\Gamma_{N_{i}}$. In the literature, $\Gamma_{N_{i}}$ can be calculated by
summing up all possible decaying channels of the Majorana neutrino, which is
explained in Appendix D. From Equation 63 as well as the nearly-degenerate
condition $m_{N_{1}}\approx m_{N_{2}}$, we can simplify the ratio between
$\Gamma_{N_{1}}$ and $\Gamma_{N_{2}}$ as

$\frac{\Gamma_{N_{2}}}{\Gamma_{N_{1}}}\approx\frac{|U_{eN_{2}}|^{2}a_{e}(m_{N_{1}})+|U_{\mu
N_{2}}|^{2}a_{\mu}(m_{N_{1}})+|U_{\tau
N_{2}}|^{2}a_{\tau}(m_{N_{1}})}{|U_{eN_{1}}|^{2}a_{e}(m_{N_{1}})+|U_{\mu
N_{1}}|^{2}a_{\mu}(m_{N_{1}})+|U_{\tau N_{1}}|^{2}a_{\tau}(m_{N_{1}})},$ (18)

where the meaning of the factors $a_{l}(l=e,\mu,\tau)$ is explained in
Appendix D. The equation will be used in the analytical analysis in the next
subsection. However, here we emphasize that we treat the lifetime of the
sterile neutrino $\tau_{N}=\hbar/\Gamma_{N}$ as a free parameter that can be
measured by LHCb experiments in our experimental analysis. In Sec. III we
explain in detail our treatment and justify that the treatment is consistent
with the result in Appendix D.

###  II.2 CP violation

In the following part we define that

$\mathcal{A}_{CP}=\frac{S^{-}}{S^{+}}$ (19)

to measure the size of the CP violation, where $S^{\pm}$ are defined as

$S^{\pm}=\Gamma^{+}\pm\Gamma^{-}.$ (20)

Then, following Eqs (13) and (2), we have

$S^{-}=4|U_{\mu N_{1}}|^{2}|U_{\mu N_{2}}|^{2}\sin(2\phi_{\mu
N_{1}}-2\phi_{\mu N_{2}})\int d\Phi_{4}\mathrm{Im}M_{12}.$ (21)

Here $\mathrm{Im}M_{12}$ represents the imaginary part of $M_{12}$. As shown
by Appendix B, the imaginary part of $M_{12}$ is proportional to the imaginary
part of $P_{N_{1}}P_{N_{2}}^{*}$,

$\mathrm{Im}P_{N_{1}}P_{N_{2}}^{*}=\frac{m_{N_{2}}\Gamma_{N_{2}}(p_{N}^{2}-m_{N_{1}}^{2})-m_{N_{1}}\Gamma_{N_{1}}(p_{N}^{2}-m_{N_{2}}^{2})}{[(p_{N}^{2}-m_{N_{1}}^{2})^{2}+(m_{N_{1}}\Gamma_{N_{1}})^{2}][(p_{N}^{2}-m_{N_{2}}^{2})^{2}+(m_{N_{2}}\Gamma_{N_{2}})^{2}]}.$
(22)

The physics meaning of Equation (21) includes two parts. First, the appearance
of the factor $\sin(2\phi_{\mu N_{1}}-2\phi_{\mu N_{2}})$ implies that the
physics origin of the CP violation is the difference between the CP-odd phases
of mixing parameters between the two heavy neutrinos with the common ones. We
define the CP-odd phase difference $\vartheta_{12}$ as

$\vartheta_{12}=2(\phi_{\mu N_{1}}-\phi_{\mu N_{2}}),$ (23)

which is the key parameter in deciding on the size of the CP violation.
Second, the CP violation of the processes in consideration arises as a result
of the interference of the contributions of the two generations of Majorana
neutrinos, which is shown explicitly by the factor
$\mathrm{Im}P_{N_{1}}P_{N_{2}}^{*}$. In case that $m_{N_{1}}=m_{N_{2}}$,
$\mathrm{Im}P_{N_{1}}P_{N_{2}}^{*}$ vanishes and so does the CP-violating term
$S^{-}$.

The sum $S^{+}$ can be written as

$S^{+}=2\int d\Phi_{4}\Big{[}|U_{\mu N_{1}}|^{4}M_{11}+|U_{\mu
N_{2}}|^{4}M_{22}+2|U_{\mu N_{1}}|^{2}|U_{\mu
N_{2}}|^{2}\mathrm{Re}M_{12}\cos\vartheta_{12}\Big{]},$ (24)

where $\mathrm{Re}M_{12}$ represents the real part of $M_{12}$, of which the
corresponding key factor is the real part of $P_{N_{1}}P_{N_{2}}^{*}$,

$\mathrm{Re}P_{N_{1}}P_{N_{2}}^{*}=\frac{(p_{N}^{2}-m_{N_{1}}^{2})(p_{N}^{2}-m_{N_{2}}^{2})+m_{N_{1}}m_{N_{2}}\Gamma_{N_{1}}\Gamma_{N_{2}}}{[(p_{N}^{2}-m_{N_{1}}^{2})^{2}+(m_{N_{1}}\Gamma_{N_{1}})^{2}][(p_{N}^{2}-m_{N_{2}}^{2})^{2}+(m_{N_{2}}\Gamma_{N_{2}})^{2}]}.$
(25)

Note that Eqs (16) and (17) show that $\Gamma_{N_{1}}\times\int
d\Phi_{4}M_{11}$($\Gamma_{N_{2}}\times\int d\Phi_{4}M_{22}$) is only the
function of the heavy neutrino mass $m_{N_{1}}$($m_{N_{2}}$) and the two
functions are of the same form. Considering the nearly-degenerate situation
$m_{N_{1}}\approx m_{N_{2}}$, we have

$\frac{\int d\Phi_{4}M_{22}}{\int
d\Phi_{4}M_{11}}\approx\frac{\Gamma_{N_{1}}}{\Gamma_{N_{2}}}.$ (26)

Then, the relative size of the CP violation $\mathcal{A}_{CP}$ can be written
in the following simple form

$\mathcal{A}_{CP}=\frac{2\delta_{I}\kappa\sin\vartheta_{12}}{1+\kappa^{2}\frac{\Gamma_{N_{1}}}{\Gamma_{N_{2}}}+2\kappa\delta_{R}\cos\vartheta_{12}}.$
(27)

Here, we define that

$\delta_{R}=\frac{\int d\Phi_{4}\mathrm{Re}M_{12}}{\int
d\Phi_{4}M_{11}},~{}\delta_{I}=\frac{\int d\Phi_{4}\mathrm{Im}M_{12}}{\int
d\Phi_{4}M_{11}},$ (28)

to measure the relative size of the interference term of the two heavy
neutrinos. The parameter $\kappa$ is defined as

$\kappa=\frac{|U_{\mu N_{2}}|^{2}}{|U_{\mu N_{1}}|^{2}}.$ (29)

Figure 2: The parameters $\delta_{I}$ and $\delta_{R}$ as a function of $y$
for different values of $m_{N_{1}}$. Note that since the deviation of
$\delta_{I}/\delta_{R}$ between different choices of $m_{N_{1}}$ is very small
(detailed numerical result shows that the difference is less than 1%), to
distinguish them, we use different styles of plot, i.e., line plot and scatter
point plot, to represent the results for different choices of $m_{N_{1}}$.

In order to make a more detailed analysis of this result, we make another
assumption about the properties of the heavy neutrinos that the two Majorana
neutrinos have approximately the same mixing parameters with the three
generations of normal neutrinos, namely,

$\frac{|U_{lN_{2}}|^{2}}{|U_{lN_{1}}|^{2}}\approx 1~{}(l=e,\mu,\tau).$ (30)

Thus, the right-hand side of Equation (18) simply gives one. The size of the
CP violation $\mathcal{A}_{CP}$ can be further simplified as

$\mathcal{A}_{CP}=\frac{\delta_{I}\sin\vartheta_{12}}{1+\delta_{R}\cos\vartheta_{12}}.$
(31)

The final result Equation (31) is only the function of $m_{N_{1}}$,
$m_{N_{2}}$ and the angle $\vartheta_{12}$, or equivalently, $m_{N_{1}}$,
$\Delta m_{N}$ and $\vartheta_{12}$. A natural way to measure the size of
$\Delta m_{N}$ is to compare it with the decay width of the heavy neutrino
$\Gamma_{N}$; thus, we define that

$y=\frac{\Delta m_{N}}{\Gamma_{N}}.$ (32)

Note that $\delta_{I}$ and $\delta_{R}$ rely both on the values of $m_{N_{1}}$
and $y$. However, numerical result shows that in the mass region under
consideration, namely, 0.5 GeV $<m_{N_{1}}<$ 3.5 GeV, the variation of
$\delta_{I}$ and $\delta_{R}$ due to different choices of $m_{N}$ is very
small (less than $1\%$). In Fig. 2, we present the numerical results of
$\delta_{I}$ and $\delta_{R}$ as a function of $y$ under the condition that
$m_{N_{1}}$ equals 1, 2 and 3 GeV. In conclusion, $\delta_{I}$ and
$\delta_{R}$ can be considered as nearly independent of the choice of
$m_{N_{1}}$ in the mass region 0.5 GeV $<m_{N_{1}}<$ 3.5 GeV. Thus, the size
of the CP violation $\mathcal{A}_{CP}$ is also independent of $m_{N_{1}}$. The
relative size of $\mathcal{A}_{CP}$ reaches its maximum when $y$ is around 1.
The above conclusions are all in accordance with the three-body case studied
in Refs. Cvetic:2013eza ; Cvetic:2015ura . Note that the conclusion does not
hold true when the mass interval is larger than several tens of GeV. Reference
Najafi:2020dkp  studies CP violation of a similar mechanism in $W$ boson
decays and it shows that under the condition that $m_{N_{1}}$ equals 5, 10, 20
and 60 GeV, $\mathcal{A}_{CP}$ shows considerable differences. Another
significant difference is that in our result $\mathcal{A}_{CP}$ reaches
maximum when $y$ is around 1 while in Ref. Najafi:2020dkp  $y$ is between 0.3
and 0.5. See Ref. Najafi:2020dkp  for details. In Fig. 3, we give the
numerical result of $\mathcal{A}_{CP}$ as a function of $y$ under different
choices of $\vartheta_{12}$.

Finally, the decay width for the process $B_{s}^{0}\rightarrow
D_{s}^{-}\mu^{+}\mu^{+}\pi^{-}$ and its CP-conjugate process can be written in
the more compact form,

$\Gamma^{\pm}=2|U_{\mu
N_{1}}|^{4}[1+\delta_{R}\cos\vartheta_{12}\pm\delta_{I}\sin\vartheta_{12}]\int
d\Phi^{4}M_{11},$ (33)

We define the averaged branching ratio $\mathcal{B}r_{avr}$ for experimental
analysis in the next section,

$\mathcal{B}r_{avr}=\frac{1}{2\Gamma_{B_{s}}}(\Gamma^{+}+\Gamma^{-}),$ (34)

here, $\Gamma_{B_{s}}=4.362\times 10^{-13}~{}\mathrm{GeV}$ is the total decay
width of $B_{s}$ meson Zyla:2020zbs .

Figure 3: The size of the CP violation $\mathcal{A}_{C}P$ as a function of $y$
for different choices of $\vartheta_{12}$.

##  III Experimental analysis

In this part, we evaluate the possibility that such CP violation is observed
by LHCb in its upgrade II LHCb:2018roe . In the experiment, the most important
quantity is the absolute size of the averaged branching ratio
$\mathcal{B}r_{avr}$. In this section, we do not distinguish the two sterile
neutrinos due to the degeneracy and use $N$ to represent both of them.

The expected number of $B_{s}^{0}$ mesons produced on LHCb can be estimated as

$N_{B_{s}}=\mathcal{L}\times\sigma_{b\bar{b}}\times f(b\rightarrow
B_{s}^{0}),$ (35)

where $\mathcal{L}\approx 200~{}\mathrm{fb^{-1}}$ is the expected integrated
luminosity of LHCb until 2035 LHCb:2018roe , $\sigma_{b\bar{b}}\simeq
144~{}\mathrm{\mu b}$ the $b\bar{b}$ cross section within the LHCb covered
$\eta$ range ($2<\eta<5$) Aaij:2016avz , and $f(b\rightarrow B_{s}^{0})\simeq
4.4\%$ the hardronization factor of b quark to $B_{s}$ which is estimated from
Ref. Aaij:2019pqz 333Ref. Aaij:2019pqz  shows that the ratio between the
production fraction of $\Lambda_{b}$ hardrons and the sum of the fraction of
$B^{-}$ and $\bar{B}^{0}$ is around 0.259 (averaged between
$4~{}\mathrm{GeV}<p_{T}<25~{}\mathrm{GeV}$ and $2<\eta<5$), while that between
$\bar{B}^{0}_{s}$ and the sum of $B^{-}$ and $\bar{B}^{0}$ is 0.122. Since
$B^{\pm}$, $B^{0}/\bar{B}^{0}$, $\bar{B}^{0}_{s}/B^{0}_{s}$ and
$\Lambda_{b}^{0}/\bar{\Lambda}^{0}_{b}$ make the majority of $b\bar{b}$
products, the fraction $f(b\bar{b}\rightarrow\Lambda_{b})$ is estimated as
$0.122/(1+0.259+0.122)\times 0.5\sim 0.044$.. The final result is
$N_{B_{s}}\approx 1.9\times 10^{12}$.

In order to get the expected sensitivity of LHCb on the processes in
consideration, we also need to know the detection efficiency of the LHCb
detector $\epsilon~{}(B_{s}\rightarrow D_{s}\mu\mu\pi)$, which contains the
contribution from geometrical acceptance, trigger and selection requirements
and particle identification LHCb:2013svv . A systematic evaluation of the
detection efficiency requires a Monte Carlo simulation under a LHCb
configuration as well as considering final state radiation generation and
interaction of the produced particles with the detector and its response
LHCb:2013svv . Here, however, we use an indirect approach to give an
approximation to the detection efficiency, whose accuracy is enough for our
calculation. Reference LHCb:2015wdu  shows that the simulated detection
efficiency of the process $B_{s}^{0}\rightarrow\phi(K^{+}K^{-})\mu^{+}\mu^{-}$
is $1.1\%$. The main difference between
$B_{s}^{0}\rightarrow\phi(K^{+}K^{-})\mu^{+}\mu^{-}$ and $B_{s}^{0}\rightarrow
D_{s}^{-}\pi^{-}\mu^{+}\mu^{+}$ comes from the replacement of $K^{-}$ with
$D_{s}^{-}$ 444The difference between the detection of $K$ and $\pi$ as well
as between $\mu^{+}$ and $\mu^{-}$ is very small and can be overlooked here..
$D_{s}$ is reconstructed by the golden mode $D_{s}\rightarrow KK\pi$, which
requires two additional charged tracks and thus would reduce the detection
efficiency. From Ref. LHCb:2020ist , we know that the ratio between the
detection efficiency of $B_{s}^{0}\rightarrow K^{-}\mu^{+}\nu$ and
$B_{s}^{0}\rightarrow D_{s}^{-}\mu^{+}\nu$ is around 0.733 (averaged over the
full $q^{2}$ range). Thus, we estimate the detection efficiency of
$B_{s}\rightarrow D_{s}\pi\mu\mu$ to be $1.1\%\times 0.733\approx 0.81\%$. We
note that the accuracy of this approximation is enough for magnitude
estimation.

Another factor we need to consider is the efficiency loss due to the flight of
the long-lived particle, i.e., the intermediate sterile neutrinos $N$, in the
detector. The sterile neutrinos are produced nearly on shell and would travel
for certain distance before decaying into its aftermath Dib:2014iga . We
include this effect by adding another factor $\mathcal{P}_{N}$ to the total
detection efficiency, which relies on the lifetime of the sterile neutrino
$\tau_{N}$ as

$\mathcal{P}_{N}=1-\exp{(-L_{D}/L_{N})},$ (36)

where $L_{D}\sim$ 1 m is the length of the detector and

$L_{N}=c\tau_{N}\gamma_{N}\beta_{N}$ (37)

is the decay length of the sterile neutrino Dib:2014iga . Here $c$ represents
the light speed and $\gamma_{N}\beta_{N}$ is the Lorentz time dilation factor
of $N$. To our knowledge, the choice of $\gamma_{N}\beta_{N}$ between 1 and 10
is common in the literature such as Refs. Cvetic:2016fbv ; Dib:2014iga  based
on the realistic condition of collider experiments. In this study, we take
that $\gamma_{N}\beta_{N}$ equals 4 instead of doing detailed calculation of
$\gamma_{N}\beta_{N}$ since the former is enough for magnitude estimation. For
$\tau_{N}=1000$ ps , the factor $\mathcal{P}_{N}$ is about $56\%$, while for
$\tau_{N}\leq 100$ ps $\mathcal{P}_{N}$ is very close to 1 and the effect is
almost negligible.

Following the practice in experimental research of sterile neutrinos such as
Refs. LHCb:2014osd ; LHCb:2016inz , we take $\tau_{N}$ as a free parameter
that can be measured by LHCb experiment to avoid more complication, instead of
calculating $\tau_{N}$ through $\tau_{N}=\hbar/\Gamma_{N}$. We assume that the
sterile neutrino lifetime $\tau_{N}=[100,1000]$ ps, which is within the
acceptance of LHCb. Here, we justify that this assumption is consistent with
Equation (63) under a certain choice of the size of $|U_{lN}|^{2}$. It can be
read from Ref. Deppisch:2015qwa  that around 1 GeV the currently known upper
bound on $|U_{eN}|^{2}$ is between $10^{-7}$ and $10^{-8}$ , and the upper
bound on $|U_{\mu N}|^{2}$ is between $10^{-4}$ and $10^{-5}$, while that on
$|U_{\tau N}|^{2}$ is between $10^{-3}$ and $10^{-2}$. The possible region
that the corresponding lifetime of the sterile neutrino can lie in is drawn in
Fig. 4. The plot shows that within the mass range
$[1~{}\mathrm{GeV},~{}4~{}\mathrm{GeV}]$, the choice that
$\tau_{N}=[100,1000]$ ps is acceptable.

In Table. 1, we present the expected number of events for certain choices of
related parameters, based on the experimental ability discussed above. We use
$\mathcal{N}_{+}$/$\mathcal{N}_{-}$ to represent the event numbers of
$B_{s}^{0}\rightarrow
D_{s}^{-}\mu^{+}\mu^{+}\pi^{-}$/$\bar{B_{s}^{0}}\rightarrow
D_{s}^{+}\mu^{-}\mu^{-}\pi^{+}$. The expected number of events at LHCb upgrade
II is estimated as

$\mathcal{N}_{\pm}=N_{B_{s}}\times\epsilon~{}(B_{s}\rightarrow
D_{s}\mu\mu\pi)\times\mathcal{P}_{N}\times\frac{\Gamma^{\pm}}{\Gamma_{B_{s}}}.$
(38)

For $|U_{\mu N}|^{2}=10^{-4}$, several hundreds of events can be expected, and
there is significant difference between $\mathcal{N}_{+}$ and
$\mathcal{N}_{-}$. For $|U_{\mu N}|^{2}=10^{-5}$, event numbers decrease to a
few, and the observable CP violation is not that significant. However, it
should be noted that our result is based on a relatively conservative
estimation on the experimental ability of LHCb since we use the detection
efficiency of previous LHCb experiments. During the upgrade II, certain
detection ability of LHCb detectors will be improved LHCb:2018roe  and it is
possible that more events can be observed.

|  $\tau_{N}=$ 1000 ps |  $\tau_{N}=$ 100 ps
---|---|---
|  $\mathcal{N}_{+}$ |  $\mathcal{N}_{-}$ |  $\mathcal{N}_{+}$ |  $\mathcal{N}_{-}$
$|U_{\mu N}|^{2}=10^{-4}$ |  665 |  390 |  119 |  70
$|U_{\mu N}|^{2}=10^{-5}$ |  7 |  4 |  1 |  0
Table 1: Expected event numbers for $B_{s}^{0}\rightarrow
D_{s}^{-}\mu^{+}\mu^{+}\pi^{-}$ ($\mathcal{N}_{+}$) and
$\bar{B_{s}^{0}}\rightarrow D_{s}^{+}\mu^{-}\mu^{-}\pi^{+}$
($\mathcal{N}_{+}$) under the assumption that $|U_{\mu N}|^{2}$ equals
$10^{-4}$ and $10^{-5}$. The relative parameters are chosen as follows:
$m_{N}=2~{}\mathrm{GeV}$, $y=1$, $\vartheta_{12}=\pi/4$. Figure 4: The shaded
region represent the possible region of the lifetime of the sterile neutrino
$N$. The black line represents the case when $|U_{eN}|^{2}=10^{-7},|U_{\mu
N}|^{2}=10^{-4},|U_{\tau N}|^{2}=10^{-2}$, and the blue line represents the
case when $|U_{eN}|^{2}=10^{-8},|U_{\mu N}|^{2}=10^{-5},|U_{\tau
N}|^{2}=10^{-3}$. Figure 5: The upper bound on the heavy-light mixing
parameter $|U_{\mu N}|^{2}$ under the assumption that no positive signal about
the processes is observed $B_{s}^{0}\rightarrow
D_{s}^{-}\mu^{+}\mu^{+}\pi^{-}$ on LHCb. The shaded regions represents the
excluded region given by previous experiments including NuTeV Vaitaitis:1999wq
, BEBC CooperSarkar:1985nh , Belle Liventsev:2013zz  and Delphi Abreu:1996pa .

On the other hand, suppose that such modes were not observed on the LHCb
upgrade II, we can constrain the upper bound on the heavy-light mixing
parameter $|U_{\mu N}|^{2}$ by requiring the total number of events to be
lower than some threshold. Reference Feldman:1997qc  shows that the
experimental sensitivity on $|U_{\mu N}|^{2}$ at $95\%$ confidence level
(C.L.) under a background-free environment is obtained for
$\mathcal{N}_{\text{events}}=3.09$. Thus, the upper bound on $|U_{\mu N}|^{2}$
at $95\%$ C.L. is obtained by requiring that

$N_{B_{s}}\times\epsilon~{}(B_{s}\rightarrow
D_{s}\mu\mu\pi)\times\mathcal{P}_{N}\times\mathcal{B}r_{\mathrm{avr}}=3.09.$
(39)

The numerical result of Equation (39) is shown in Fig. 5. For comparison, we
also include currently known upper bounds on $|U_{\mu N}|^{2}$ given by other
experiments including NuTeV Vaitaitis:1999wq , BEBC CooperSarkar:1985nh ,
Belle Liventsev:2013zz , and Delphi Abreu:1996pa  in the plot. The plot shows
that during the mass region $[1~{}\mathrm{GeV},3~{}\mathrm{GeV}]$, the channel
gives a comparable or even stronger constraint on the size of $|U_{\mu
N}|^{2}$. Specifically, for $m_{N}$ equals 2 GeV, the constraint on $|U_{\mu
N}|^{2}$ goes as low as $6\times 10^{-6}$.

##  IV Summary and conclusions

In this paper, we study the lepton-number-violating (LNV) process
$B_{s}^{0}\rightarrow D_{s}^{-}\mu^{+}\mu^{+}\pi^{-}$ and its CP-conjugate
process $\bar{B_{s}^{0}}\rightarrow D_{s}^{+}\mu^{-}\mu^{-}\pi^{+}$ that are
induced by two nearly-degenerate Majorana neutrinos and explore the
possibility for searching for the CP violation in such processes. We point out
that the physics origin of the CP violation is $\vartheta_{12}$, which is
defined as the difference between the CP-odd phases of mixing parameters
between two generations of heavy neutrinos with the normal ones. The CP
violation becomes considerable when the masses of the two generations of heavy
neutrinos are nearly degenerate but have a nonzero difference $\Delta m_{N}$.
The numerical result draws the following conclusion. First, the relative size
of the CP violation $\mathcal{A}_{CP}$ is only the function of the mass
difference $\Delta m_{N}$ and $\vartheta_{12}$ and $\mathcal{A}_{CP}$ is
nearly independent of the absolute mass of the lighter heavy neutrino in the
mass region we considered. Second, $\mathcal{A}_{CP}$ reaches its maximum when
$\Delta m_{N}$ is around the size of the decay width of the intermediate
sterile neutrino and the maximum value depends on the CP-odd phase difference
$\vartheta_{12}$. It should be noted that the above conclusion is drawn under
the assumption that the two Majorana neutrinos have approximately the same
mixing parameters with the three normal neutrinos.

We also analyze the possibility that such CP violation is observed by LHCb
during its upgrade II. It is shown that under a current constraint on the
heavy-light neutrino mixing parameter, namely, $10^{-4}<|U_{\mu
N}|^{2}<10^{-5}$, it is possible that such a CP violation can be observed with
the LHCb experimental ability. We also give the upper bound on the heavy-light
mixing parameter under the assumption that no positive signal of the processes
is observed. The result shows that such modes can give a complementary
constraint on the heavy-light mixing parameter compared with previous
experiments including NuTeV, BEBC, Belle, and Delphi in the mass region 1 GeV
$<m_{N}<$ 3 GeV. Thus, we note that it is worthwhile searching for such modes
on LHCb due to the possibility of both observing new types of CP violation in
$B_{s}$ meson decays and setting complementary constraints on $|U_{\mu
N}|^{2}$.

##  V Acknowledgement

This work is supported by National Natural Science Foundation of China (Grant
No. 12075003).

## APPENDIX

##  Appendix A THE FOUR-BODY PHASE SPACE $d\Phi_{4}$

The four-body phase space $d\Phi_{4}$ is

$d\Phi_{4}=\frac{1}{2m_{B_{s}}}\frac{1}{(2\pi)^{8}}d_{4},$ (40)

where the factor $d_{4}$ is defined as

$d_{4}=\frac{\mathrm{d}^{3}\bm{p}_{D_{s}}}{(2E_{D_{s}})}\frac{d^{3}\bm{p}_{1}}{(2E_{1})}\frac{d^{3}\bm{p}_{2}}{(2E_{2})}\frac{d^{3}\bm{p}_{\pi}}{(2E_{\pi})}\delta^{4}(p_{B_{s}}-p_{D_{s}}-p_{\mu_{1}}-p_{\mu_{2}}-p_{\pi}).$
(41)

Remember the fact that

$\frac{d^{3}\bm{p}_{N}}{(2E_{N})}=d^{4}p_{N}\delta_{+}(p_{N}^{2}-m_{N}^{2}),$
(42)

then it is straightforward to prove that $d_{4}$ can be factorized as

$d_{4}=d_{3}(B_{s}\rightarrow
D_{s}\mu_{1}N)d_{2}(N\rightarrow\mu_{2}\pi)dp_{N}^{2}.$ (43)

Here, $d_{3}(B_{s}\rightarrow D_{s}\mu_{1}N)$ is the corresponding factor in
the three-body phase space for the subprocess $B_{s}\rightarrow D_{s}\mu_{1}N$
and $d_{2}(N\rightarrow\mu_{2}\pi)$ is that for $N\rightarrow\mu_{2}\pi$.
$p_{N}$ is the four-momenta of the intermediate sterile neutrino. In the
following we set $t_{2}=p_{N}^{2}$. In the rest frame of $\mu_{1}-N$,
$d_{3}(B_{s}\rightarrow D_{s}\mu_{1}N)$ can be simplified as

$d_{3}(B_{s}\rightarrow
D_{s}\mu_{1}N)=\frac{1}{64m_{B_{s}}^{2}}\lambda^{\frac{1}{2}}(m_{B_{s}}^{2},m_{D_{s}}^{2},t_{1})\lambda^{\frac{1}{2}}(t_{1},m_{\mu}^{2},m_{N}^{2})\frac{1}{t_{1}}dt_{1}d\Omega^{*}_{\mu_{1}}d\Omega_{D_{s}},$
(44)

where $t_{1}=(p_{N}+p_{\mu_{1}})^{2}=(p_{B_{s}}-p_{D_{s}})^{2}$, and
$d\Omega^{*}_{\mu_{1}}=d\cos\theta_{1}d\varphi_{1}$ is the solid angle of
$\mu_{1}$ in the $\mu_{1}-N$ rest frame, while
$d\Omega_{D_{s}}=d\cos\theta_{2}\varphi_{2}$ is that of $D_{s}$ in the rest
frame of $B_{s}$. The function $\lambda(x,y,z)$ is the kinematic Källen
function, $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2yz-2xz$. In the rest frame of
$N$, $d_{2}(N\rightarrow\mu_{2}\pi)$ is written as

$d_{2}(N\rightarrow\mu_{2}\pi)=\frac{1}{8t_{2}}\lambda^{\frac{1}{2}}(t_{2},m_{\mu}^{2},m_{\pi}^{2})d\Omega_{\mu_{2}},$
(45)

where $d\Omega_{\mu_{2}}=d\cos\theta_{3}d\varphi_{3}$ is the solid angle of
$\mu_{2}$ in the rest frame of $N$. As is shown in the next part, the square
of the amplitude $|\mathcal{M}|^{2}$ is independent of $\varphi_{1}$ and
$\Omega_{D_{s}}$; thus, they can be integrated and give a factor $8\pi^{2}$.
As a result, the four-body phase space $d_{4}$ is

$d_{4}=\frac{\pi^{2}}{64m_{B_{s}}^{2}}\lambda^{\frac{1}{2}}(m_{B_{s}}^{2},m_{D_{s}}^{2},t_{1})\lambda^{\frac{1}{2}}(t_{1},m_{\mu}^{2},m_{N}^{2})\lambda^{\frac{1}{2}}(t_{2},m_{\mu}^{2},m_{\pi}^{2})\frac{1}{t_{1}t_{2}}dt_{1}dt_{2}d\cos\theta_{1}d\cos\theta_{3}d\varphi_{3}.$
(46)

Here, $t_{1}$ and $t_{2}$ satisfy that $(\sqrt{t_{2}}+m_{\mu})^{2}\leq
t_{1}\leq(m_{B_{s}}-m_{D_{s}})^{2}$ and $t_{2}\geq(m_{\mu}+m_{\pi})^{2}$. The
definitions of $\theta_{1},\phi_{1},\theta_{3},\phi_{3}$ and the reference
frames we used are shown in Fig. 6.

Figure 6: The reference frames we use and the definitions of the solid angles.

##  Appendix B EXPLICIT FORM OF THE SQUARED MATRIX $M$

For simplicity, we define in the following that

$K=G_{F}^{2}V_{cb}V_{ud}f_{\pi}$ (47)

and

$T^{\pm}=\bar{u}(\mu_{2})\not{p}_{\pi}\gamma_{\mu}(1\mp\gamma_{5})v(\mu_{1})\langle
D_{s}^{-}|c(0)\gamma_{\mu}b(0)|B_{s}^{0}\rangle.$ (48)

$T^{\pm}$ contains all the spinor structures of the decaying processes. It can
be checked that $|T^{+}|^{2}=|T^{-}|^{2}\equiv|T|^{2}$. Then the elements of
the squared amplitude matrix $M$ are written as

$M_{ij}=|K|^{2}m_{N_{i}}m_{N_{j}}P_{N_{i}}P_{N_{j}}^{*}|T|^{2}$ (49)

Apart from $P_{N_{i}}P_{N_{j}}^{*}$, other factors in $M_{ij}$ are always
real. Thus, the real/imaginary part of $M_{12}$ is proportional to the
real/imaginary part of $P_{N_{1}}P_{N_{2}}^{*}$. The completed form of
$|T|^{2}$ is

$|T|^{2}=C_{00}f_{0}^{2}(t_{1})+C_{01}f_{0}(t_{1})f_{+}(t_{1})+C_{11}f_{+}^{2}(t_{1}).$
(50)

Here, the coefficients are

$\displaystyle C_{00}$
$\displaystyle=\frac{8}{t_{1}^{2}}(m_{B_{s}}^{2}-m_{D_{s}}^{2})^{2}\Big{\\{}m_{\pi}^{2}[t_{1}(p_{\mu_{1}}\cdot
p_{\mu_{2}})-2(q\cdot p_{\mu_{1}})(q\cdot p_{\mu_{2}})]+4(p_{\pi}\cdot
p_{\mu_{2}})(q\cdot p_{\pi})(q\cdot p_{\mu_{1}})-2t_{1}(p_{\pi}\cdot
p_{\mu_{1}})(p_{\pi}\cdot p_{\mu_{2}})\Big{\\}},$ (51) $\displaystyle C_{01}$
$\displaystyle=\frac{16}{t}(m_{B_{s}}^{2}-m_{D_{s}}^{2})\Big{\\{}(q\cdot
p_{\mu_{1}})\big{[}2(p_{\pi}\cdot p_{\mu_{2}})(P\cdot
p_{\pi})-m_{\pi}^{2}(P\cdot p_{\mu_{2}})\big{]}+(P\cdot
p_{\mu_{1}})\big{[}2(p_{\pi}\cdot p_{\mu_{2}})(q\cdot
p_{\pi})-m_{\pi}^{2}(q\cdot p_{\mu_{2}})\big{]}\Big{\\}},$ $\displaystyle
C_{11}$ $\displaystyle=8\Big{\\{}P^{2}[m_{\pi}^{2}(p_{\mu_{1}}\cdot
p_{\mu_{2}})-2(p_{\pi}\cdot p_{\mu_{1}})(p_{\pi}\cdot
p_{\mu_{2}})]-2m_{\pi}^{2}(P\cdot p_{\mu_{1}})(P\cdot
p_{\mu_{2}})+4(p_{\pi}\cdot p_{\mu_{2}})(P\cdot p_{\pi})(P\cdot
p_{\mu_{1}})\Big{\\}},$

where we define that

$P=p_{B_{s}}+p_{D_{s}}-\frac{m_{B_{s}}^{2}-m_{D_{s}}^{2}}{t_{1}}q.$ (52)

Most inner products in Equation (51) can be written as functions of
$t_{1},t_{2},\theta_{1},\theta_{3},\varphi_{3}$ directly, namely,

$\displaystyle p_{D_{s}}\cdot
q=\frac{1}{2}(m_{B_{s}}^{2}-m_{D_{s}}^{2}-t_{1}),~{}p_{D_{s}}\cdot
p_{\mu_{1}}=E_{D_{s}}^{*}E_{\mu_{1}}^{*}-|\bm{p}_{D_{s}}^{*}||\bm{p}_{\mu_{1}}^{*}|\cos\theta_{1},~{}q\cdot
p_{\mu_{1}}=\frac{1}{2}(t_{1}+m_{\mu}^{2}-t_{2}),$ (53) $\displaystyle
p_{\mu_{2}}\cdot
p_{\pi}=\frac{1}{2}(t_{2}-m_{\mu}^{2}-m_{\pi}^{2}),~{}~{}~{}p_{\mu_{1}}\cdot
p_{\pi}=E_{\mu_{1}}E_{\pi}+|\bm{p}_{\mu_{1}}||\bm{p}_{\mu_{2}}|\cos\theta_{3},~{}~{}~{}~{}p_{\mu_{1}}\cdot
p_{\mu_{2}}=E_{\mu_{1}}E_{\mu_{2}}-|\bm{p}_{\mu_{1}}||\bm{p}_{\mu_{2}}|\cos\theta_{3}.$

Here, we list the explicit forms for the energies and three-momenta that
appear in the above terms (the superscript “$*$” represents the value in the
$W^{*}$ rest frame),

$\displaystyle
E_{D_{s}}^{*}=\frac{1}{2\sqrt{t_{1}}}(m_{B_{s}}^{2}-m_{D_{s}}^{2}-t_{1}),~{}E_{\mu_{1}}^{*}=\frac{1}{2\sqrt{t_{1}}}(t_{1}+m_{\mu}^{2}-t_{2}),~{}|\bm{p}_{D_{s}}^{*}|=\frac{1}{2\sqrt{t_{1}}}\lambda^{\frac{1}{2}}(m_{B_{s}}^{2},m_{D_{s}}^{2},t_{1}),$
(54)
$\displaystyle|\bm{p}_{\mu_{1}}^{*}|=\frac{1}{2\sqrt{t_{1}}}\lambda^{\frac{1}{2}}(t_{1},m_{\mu}^{2},t_{2}),\qquad
E_{\mu_{1}}=\frac{1}{2\sqrt{t_{2}}}(t_{1}-m_{\mu}^{2}-t_{2}),~{}~{}E_{\mu_{2}}=\frac{1}{2\sqrt{t_{2}}}(t_{2}-m_{\mu}^{2}+m_{\pi}^{2}),$
$\displaystyle
E_{\pi}=\frac{1}{2\sqrt{t_{2}}}(t_{2}-m_{\mu}^{2}+m_{\pi}^{2}),\qquad|\bm{p}_{\mu_{1}}|=\frac{1}{2\sqrt{t_{2}}}\lambda^{\frac{1}{2}}(t_{1},m_{\mu}^{2},t_{2}),~{}~{}|\bm{p}_{\mu_{2}}|=\frac{1}{2\sqrt{t_{2}}}\lambda^{\frac{1}{2}}(t_{2},m_{\mu}^{2},m_{\pi}^{2}),$
$\displaystyle
E_{N}^{*}=\frac{1}{2\sqrt{t_{1}}}(t_{1}+t_{2}-m_{\mu}^{2}),\qquad~{}|\bm{p}_{N}^{*}|=\frac{1}{2\sqrt{t_{1}}}\lambda^{\frac{1}{2}}(t_{1},t_{2},m_{\mu}^{2}).$

Note that $q=p_{\mu_{1}}+p_{\mu_{2}}+p_{\pi}$, so $q\cdot p_{\mu_{2}}$ and
$q\cdot p_{\pi}$ can also be written directly. In order to get $p_{D_{s}}\cdot
p_{\mu_{2}}$ and $p_{D_{s}}\cdot p_{\pi}$, we need to do a Lorentz
transformation from the $W*$ rest frame to the $N$ rest frame on $p_{D_{s}}$,
since only the four-momenta vector of $D_{s}$ in the $W*$ rest frame can be
written directly. We need to first rotate the $z$ axis of $\Sigma_{W^{*}}$ to
the direction of $\bm{p}_{\mu_{1}}$, then boost the four-momenta vector of
$p_{D_{s}}$ from the $W^{*}$ rest frame to the $N$ rest from. In the $W^{*}$
rest frame where the $z$ axis points to the direction of $\bm{p}_{\mu_{1}}$,
the four-momenta vector of $D_{s}$ is

$p_{D_{s}}^{*}=(E_{D_{s}}^{*},~{}|\bm{p}^{*}_{D_{s}}|\sin\theta_{1},~{}0,~{}|\bm{p}^{*}_{D_{s}}|\cos\theta_{1}).$
(55)

The Lorentz boost from the $W^{*}$ rest frame to the $N$ rest frame is

$B(\bm{v}_{N})=\left(\begin{array}[]{cccc}\gamma_{N}&~{}0&~{}0&~{}\gamma_{N}|\bm{v}_{N}|\\\
0&~{}1&~{}0&~{}0\\\ 0&~{}0&~{}1&~{}0\\\
\gamma_{N}|\bm{v}_{N}|&~{}0&~{}0&~{}\gamma_{N}\end{array}\right),$ (56)

where $\gamma_{N}$ and $\bm{v}_{N}$ are the Lorentz factor and the speed of
$N$ in the $W^{*}$ rest frame, respectively,

$\gamma_{N}=\frac{E_{N}^{*}}{m_{N}},~{}~{}~{}|\bm{v}_{N}|=\frac{|\bm{p}_{N}^{*}|}{E_{N}^{*}}.$
(57)

Here, $E_{N}^{*}$ and $|\bm{p}_{N}^{*}|$ are given in Equation (54). As a
result, $p_{D_{s}}\cdot p_{\mu_{2}}$ and $p_{D_{s}}\cdot p_{\pi}$ are written
as

$\displaystyle p_{D_{s}}\cdot
p_{\mu_{2}}=\gamma_{N}(E_{D_{s}}^{*}+|\bm{v}_{N}||\bm{p}_{D_{s}}^{*}|\cos\theta_{1})E_{\mu_{2}}-|\bm{p}_{D_{s}}^{*}||\bm{p}_{\mu_{2}}|\sin\theta_{1}\sin\theta_{3}\cos\varphi_{3}-\gamma_{N}(|\bm{p}_{D_{s}}^{*}|\cos\theta_{1}+|\bm{v}_{N}|E_{D_{s}}^{*})|\bm{p}_{\mu_{2}}|\cos\theta_{3},$
(58) $\displaystyle p_{D_{s}}\cdot
p_{\pi}~{}=\gamma_{N}(E_{D_{s}}^{*}+|\bm{v}_{N}||\bm{p}_{D_{s}}^{*}|\cos\theta_{1})E_{\pi}+|\bm{p}_{D_{s}}^{*}||\bm{p}_{\mu_{2}}|\sin\theta_{1}\sin\theta_{3}\cos\varphi_{3}+\gamma_{N}(|\bm{p}_{D_{s}}^{*}|\cos\theta_{1}+|\bm{v}_{N}|E_{D_{s}}^{*})|\bm{p}_{\mu_{2}}|\cos\theta_{3}$

##  Appendix C DECAY WIDTH FOR $B_{s}\rightarrow D_{s}\mu_{1}N_{i}$ AND FOR
$N_{i}\rightarrow\mu_{2}\pi$

Figure 7: Feynman diagram for $B_{s}\rightarrow D_{s}\mu_{1}N_{i}$

The Feynman diagram for $B_{s}\rightarrow D_{s}\mu_{1}N_{i}$ is shown in Fig.
7. The amplitude for the process is

$i\mathcal{M}(B_{s}\rightarrow D_{s}\mu_{1}N_{i})=\frac{G_{F}}{\sqrt{2}}U_{\mu
N_{i}}V_{cb}\bar{u}(p_{N_{i}})\gamma_{\mu}(1-\gamma_{5})v(p_{1})\langle
D_{s}^{-}|c(0)\gamma^{\mu}b(0)|B_{s}^{0}\rangle,$ (59)

Integrating the amplitude over three-body phase space, we have the decay width
for $B_{s}\rightarrow D_{s}\mu_{1}N_{i}$,

$\Gamma(B_{s}\rightarrow
D_{s}\mu_{1}N_{i})=\frac{G_{F}^{2}V_{cb}^{2}}{384\pi^{3}}\frac{1}{m_{B_{s}}^{3}}\int_{t_{min}}^{t_{max}}dt\frac{1}{t^{2}}\lambda^{\frac{1}{2}}(m_{B_{s}}^{2},m_{D_{s}}^{2},t)\lambda^{\frac{1}{2}}(m_{\mu}^{2},m_{N_{i}}^{2},t)\Big{[}f_{+}^{2}(t)D_{1}(t)+f_{0}^{2}(t)D_{0}(t)\Big{]},$
(60)

where $t_{min}=(m_{N_{i}}+m_{\mu})^{2}$, $t_{max}=(m_{B_{s}}-m_{D_{s}})^{2}$
and coefficients $D_{1}(t)$ and $D_{0}(t)$ are

$\displaystyle D_{1}(t)$
$\displaystyle=\Big{[}(t-m_{D_{s}}^{2})^{2}-2m_{B_{s}}^{2}(t+m_{D_{s}}^{2})+m_{B_{s}}^{4}\Big{]}\Big{[}2t^{2}-tm_{N_{i}}^{2}+m_{\mu}^{2}(2m_{N_{i}}^{2}-t)-m_{N_{i}}^{2}-m_{\mu}^{4}\Big{]}$
(61) $\displaystyle D_{0}(t)$
$\displaystyle=3(m_{B_{s}}^{2}-m_{D_{s}}^{2})\Big{[}tm_{N_{i}}^{2}+m_{\mu}^{2}(2m_{N_{i}}^{2}+t)-m_{N_{i}}^{4}-m_{\mu}^{4}\Big{]}.$

As for the process $N_{i}\rightarrow\mu_{2}\pi$, the decay width is well known
in the literature,

$\Gamma(N_{i}\rightarrow\mu\pi)=\frac{G_{F}^{2}}{16\pi}|V_{ud}|^{2}|U_{\mu
N_{i}}|^{2}f_{\pi}^{2}m_{N_{i}}\lambda^{1/2}(m_{N_{i}}^{2},m_{\mu}^{2},m_{\pi}^{2})\left[\left(1-\frac{m_{\mu}^{2}}{m_{N_{i}}^{2}}\right)^{2}-\frac{m_{\pi}^{2}}{m_{N_{i}}^{2}}\left(1+\frac{m_{\mu}^{2}}{m_{N_{i}}^{2}}\right)\right].$
(62)

##  Appendix D TOTAL DECAY WIDTH FOR THE INTERMEDIATE MAJORANA NEUTRINO

Although $\Gamma_{N_{i}}$ can be calculated through the channel-by-channel
approach, which sums up the partial decay rates for all the leptonic and
semileptonic decay modes of $N$ Atre:2009rg , for neutrino mass larger than 1
GeV the uncertainties of the hadronic parameters such as the decay constants
of the final hadronic states are large. Since we are interested in the mass
range between 1 and 4 GeV, here we use the inclusive approach introduced in
Ref. Helo:2010cw , which approximates the semileptonic decays of $N$ by its
decays into free quark-antiquark pairs and leptons and the approximation is
better than the channel-by-channel method for neutrino mass more than 1 GeV
Bondarenko:2018ptm . We refer to Ref. Helo:2010cw  for details of the
calculation. Then all the decay channels of $N$ are three-body decays and the
partial decay rates are proportional to $m_{N}^{5}$. Thus, the total decay
width can be written as

$\Gamma_{N_{i}}=\frac{G_{F}^{2}m_{N_{i}}^{5}}{96\pi^{3}}\sum_{l=e,\mu,\tau}|U_{lN}|^{2}a_{l}(m_{N_{i}}),$
(63)

where $a_{l}(m_{N_{i}})$ are dimensionless functions of $m_{N_{i}}$, and the
numerical values of $a_{l}(m_{N_{i}})$ are presented in Fig. 8.

Figure 8: Numerical values of $a_{l}(m_{N})$

## References

  * (1) Y. Fukuda et al. [Super-Kamiokande], Phys. Rev. Lett. 81, 1562-1567 (1998) doi:10.1103/PhysRevLett.81.1562 [arXiv:hep-ex/9807003 [hep-ex]].
  * (2) R. Wendell et al. [Super-Kamiokande], Phys. Rev. D 81, 092004 (2010) doi:10.1103/PhysRevD.81.092004 [arXiv:1002.3471 [hep-ex]].
  * (3) P. Minkowski, Phys. Lett. B 67, 421-428 (1977) doi:10.1016/0370-2693(77)90435-X
  * (4) M. Gell-Mann, P. Ramond and R. Slansky, Conf. Proc. C 790927, 315-321 (1979) [arXiv:1306.4669 [hep-th]].
  * (5) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980) doi:10.1103/PhysRevLett.44.912
  * (6) T. Yanagida, Prog. Theor. Phys. 64, 1103 (1980) doi:10.1143/PTP.64.1103
  * (7) L. Canetti, M. Drewes, T. Frossard and M. Shaposhnikov, Phys. Rev. D 87, 093006 (2013) doi:10.1103/PhysRevD.87.093006 [arXiv:1208.4607 [hep-ph]].
  * (8) S. Dodelson and L. M. Widrow, Phys. Rev. Lett. 72, 17-20 (1994) doi:10.1103/PhysRevLett.72.17 [arXiv:hep-ph/9303287 [hep-ph]].
  * (9) X. D. Shi and G. M. Fuller, Phys. Rev. Lett. 82, 2832-2835 (1999) doi:10.1103/PhysRevLett.82.2832 [arXiv:astro-ph/9810076 [astro-ph]].
  * (10) T. Asaka and M. Shaposhnikov, Phys. Lett. B 620, 17-26 (2005) doi:10.1016/j.physletb.2005.06.020 [arXiv:hep-ph/0505013 [hep-ph]].
  * (11) L. Canetti, M. Drewes and M. Shaposhnikov, Phys. Rev. Lett. 110, no.6, 061801 (2013) doi:10.1103/PhysRevLett.110.061801 [arXiv:1204.3902 [hep-ph]].
  * (12) S. Dell’Oro, S. Marcocci, M. Viel and F. Vissani, Adv. High Energy Phys. 2016, 2162659 (2016) doi:10.1155/2016/2162659 [arXiv:1601.07512 [hep-ph]].
  * (13) M. J. Dolinski, A. W. P. Poon and W. Rodejohann, Ann. Rev. Nucl. Part. Sci. 69, 219-251 (2019) doi:10.1146/annurev-nucl-101918-023407 [arXiv:1902.04097 [nucl-ex]].
  * (14) J. Engel and J. Menéndez, Rept. Prog. Phys. 80, no.4, 046301 (2017) doi:10.1088/1361-6633/aa5bc5 [arXiv:1610.06548 [nucl-th]].
  * (15) J. Abad, J. G. Esteve and A. F. Pacheco, Phys. Rev. D 30, 1488 (1984) doi:10.1103/PhysRevD.30.1488
  * (16) L. S. Littenberg and R. E. Shrock, Phys. Rev. Lett. 68, 443-446 (1992) doi:10.1103/PhysRevLett.68.443
  * (17) L. S. Littenberg and R. Shrock, Phys. Lett. B 491, 285-290 (2000) doi:10.1016/S0370-2693(00)01041-8 [arXiv:hep-ph/0005285 [hep-ph]].
  * (18) A. Ali, A. V. Borisov and N. B. Zamorin, Eur. Phys. J. C 21, 123-132 (2001) doi:10.1007/s100520100702 [arXiv:hep-ph/0104123 [hep-ph]].
  * (19) M. A. Ivanov and S. G. Kovalenko, Phys. Rev. D 71, 053004 (2005) doi:10.1103/PhysRevD.71.053004 [arXiv:hep-ph/0412198 [hep-ph]].
  * (20) C. Dib, V. Gribanov, S. Kovalenko and I. Schmidt, Phys. Lett. B 493, 82-87 (2000) doi:10.1016/S0370-2693(00)01134-5 [arXiv:hep-ph/0006277 [hep-ph]].
  * (21) A. Atre, V. Barger and T. Han, Phys. Rev. D 71, 113014 (2005) doi:10.1103/PhysRevD.71.113014 [arXiv:hep-ph/0502163 [hep-ph]].
  * (22) A. Atre, T. Han, S. Pascoli and B. Zhang, JHEP 05, 030 (2009) doi:10.1088/1126-6708/2009/05/030 [arXiv:0901.3589 [hep-ph]].
  * (23) J. C. Helo, S. Kovalenko and I. Schmidt, Nucl. Phys. B 853, 80-104 (2011) doi:10.1016/j.nuclphysb.2011.07.020 [arXiv:1005.1607 [hep-ph]].
  * (24) G. Cvetic, C. Dib, S. K. Kang and C. S. Kim, Phys. Rev. D 82, 053010 (2010) doi:10.1103/PhysRevD.82.053010 [arXiv:1005.4282 [hep-ph]].
  * (25) G. Cvetic and C. S. Kim, Phys. Rev. D 94, no.5, 053001 (2016) [erratum: Phys. Rev. D 95, no.3, 039901 (2017)] doi:10.1103/PhysRevD.94.053001 [arXiv:1606.04140 [hep-ph]].
  * (26) G. Cvetic and C. S. Kim, Phys. Rev. D 96, no.3, 035025 (2017) [erratum: Phys. Rev. D 102, no.1, 019903 (2020); erratum: Phys. Rev. D 102, no.3, 039902 (2020)] doi:10.1103/PhysRevD.96.035025 [arXiv:1705.09403 [hep-ph]].
  * (27) G. Cvetič and C. S. Kim, Phys. Rev. D 100, no.1, 015014 (2019) doi:10.1103/PhysRevD.100.015014 [arXiv:1904.12858 [hep-ph]].
  * (28) D. Milanes, N. Quintero and C. E. Vera, Phys. Rev. D 93, no.9, 094026 (2016) doi:10.1103/PhysRevD.93.094026 [arXiv:1604.03177 [hep-ph]].
  * (29) J. Mejia-Guisao, D. Milanés, N. Quintero and J. D. Ruiz-Alvarez, Phys. Rev. D 97, no.7, 075018 (2018) doi:10.1103/PhysRevD.97.075018 [arXiv:1708.01516 [hep-ph]].
  * (30) D. Milanés and N. Quintero, Phys. Rev. D 98, no.9, 096004 (2018) doi:10.1103/PhysRevD.98.096004 [arXiv:1808.06017 [hep-ph]].
  * (31) T. Asaka and H. Ishida, Phys. Lett. B 763, 393-396 (2016) doi:10.1016/j.physletb.2016.10.070 [arXiv:1609.06113 [hep-ph]].
  * (32) E. J. Chun, A. Das, S. Mandal, M. Mitra and N. Sinha, Phys. Rev. D 100, no.9, 095022 (2019) doi:10.1103/PhysRevD.100.095022 [arXiv:1908.09562 [hep-ph]].
  * (33) N. Quintero, G. Lopez Castro and D. Delepine, Phys. Rev. D 84, 096011 (2011) [erratum: Phys. Rev. D 86, 079905 (2012)] doi:10.1103/PhysRevD.84.096011 [arXiv:1108.6009 [hep-ph]].
  * (34) C. Barbero, G. Lopez Castro and A. Mariano, Phys. Lett. B 566, 98-107 (2003) doi:10.1016/S0370-2693(03)00773-1 [arXiv:nucl-th/0212083 [nucl-th]].
  * (35) C. Barbero, L. F. Li, G. Lopez Castro and A. Mariano, Phys. Rev. D 76, 116008 (2007) doi:10.1103/PhysRevD.76.116008 [arXiv:0709.2431 [hep-ph]].
  * (36) C. Barbero, L. F. Li, G. López Castro and A. Mariano, Phys. Rev. D 87, no.3, 036010 (2013) doi:10.1103/PhysRevD.87.036010 [arXiv:1301.3448 [hep-ph]].
  * (37) J. Mejia-Guisao, D. Milanes, N. Quintero and J. D. Ruiz-Alvarez, Phys. Rev. D 96, no.1, 015039 (2017) doi:10.1103/PhysRevD.96.015039 [arXiv:1705.10606 [hep-ph]].
  * (38) G. Zhang and B.-Q. Ma, Phys. Rev. D 103, no.3, 033004 (2021) doi:10.1103/PhysRevD.103.033004 [arXiv:2101.05566 [hep-ph]].
  * (39) G. Lopez Castro and N. Quintero, Phys. Rev. D 85, 076006 (2012) [erratum: Phys. Rev. D 86, 079904 (2012)] doi:10.1103/PhysRevD.85.076006 [arXiv:1203.0537 [hep-ph]].
  * (40) V. Gribanov, S. Kovalenko and I. Schmidt, Nucl. Phys. B 607, 355-368 (2001) doi:10.1016/S0550-3213(01)00169-9 [arXiv:hep-ph/0102155 [hep-ph]].
  * (41) G. Cvetic, C. Dib, C. S. Kim and J. D. Kim, Phys. Rev. D 66, 034008 (2002) [erratum: Phys. Rev. D 68, 059901 (2003)] doi:10.1103/PhysRevD.66.034008 [arXiv:hep-ph/0202212 [hep-ph]].
  * (42) K. Abe et al. [T2K], Phys. Rev. Lett. 107, 041801 (2011) doi:10.1103/PhysRevLett.107.041801 [arXiv:1106.2822 [hep-ex]].
  * (43) F. P. An et al. [Daya Bay], Phys. Rev. Lett. 108, 171803 (2012) doi:10.1103/PhysRevLett.108.171803 [arXiv:1203.1669 [hep-ex]].
  * (44) J. K. Ahn et al. [RENO], Phys. Rev. Lett. 108, 191802 (2012) doi:10.1103/PhysRevLett.108.191802 [arXiv:1204.0626 [hep-ex]].
  * (45) G. Cvetič, C. S. Kim and J. Zamora-Saá, J. Phys. G 41, 075004 (2014) doi:10.1088/0954-3899/41/7/075004 [arXiv:1311.7554 [hep-ph]].
  * (46) G. Cvetic, C. S. Kim, R. Kogerler and J. Zamora-Saa, Phys. Rev. D 92, 013015 (2015) doi:10.1103/PhysRevD.92.013015 [arXiv:1505.04749 [hep-ph]].
  * (47) G. Cvetič, C. S. Kim and J. Zamora-Saá, Phys. Rev. D 89, no.9, 093012 (2014) doi:10.1103/PhysRevD.89.093012 [arXiv:1403.2555 [hep-ph]].
  * (48) C. O. Dib, M. Campos and C. S. Kim, JHEP 02, 108 (2015) doi:10.1007/JHEP02(2015)108 [arXiv:1403.8009 [hep-ph]].
  * (49) G. Cvetic, C. Dib, C. S. Kim and J. Zamora-Saa, Symmetry 7, 726-773 (2015) doi:10.3390/sym7020726 [arXiv:1503.01358 [hep-ph]].
  * (50) G. Cvetic, C. S. Kim, S. Mendizabal and J. Zamora-Saa, Eur. Phys. J. C 80, no.11, 1052 (2020) doi:10.1140/epjc/s10052-020-08625-0 [arXiv:2007.04115 [hep-ph]].
  * (51) J. Zhang, T. Wang, G. Li, Y. Jiang and G. L. Wang, Phys. Rev. D 103, no.3, 035015 (2021) doi:10.1103/PhysRevD.103.035015 [arXiv:2010.13286 [hep-ph]].
  * (52) J. Zamora-Saa, JHEP 05, 110 (2017) doi:10.1007/JHEP05(2017)110 [arXiv:1612.07656 [hep-ph]].
  * (53) S. Tapia and J. Zamora-Saá, Nucl. Phys. B 952, 114936 (2020) doi:10.1016/j.nuclphysb.2020.114936 [arXiv:1906.09470 [hep-ph]].
  * (54) F. Najafi, J. Kumar and D. London, JHEP 04, 021 (2021) doi:10.1007/JHEP04(2021)021 [arXiv:2011.03686 [hep-ph]].
  * (55) A. Abada, C. Hati, X. Marcano and A. M. Teixeira, JHEP 09, 017 (2019) doi:10.1007/JHEP09(2019)017 [arXiv:1904.05367 [hep-ph]].
  * (56) A. Das, P. S. B. Dev and R. N. Mohapatra, Phys. Rev. D 97, no.1, 015018 (2018) doi:10.1103/PhysRevD.97.015018 [arXiv:1709.06553 [hep-ph]].
  * (57) D. Gorbunov and M. Shaposhnikov, JHEP 10, 015 (2007) [erratum: JHEP 11, 101 (2013)] doi:10.1088/1126-6708/2007/10/015 [arXiv:0705.1729 [hep-ph]].
  * (58) R. Aaij et al. [LHCb], [arXiv:1808.08865 [hep-ex]].
  * (59) C. Dib and C. S. Kim, Phys. Rev. D 89, no.7, 077301 (2014) doi:10.1103/PhysRevD.89.077301 [arXiv:1403.1985 [hep-ph]].
  * (60) P. A. Zyla et al. [Particle Data Group], PTEP 2020, no.8, 083C01 (2020) doi:10.1093/ptep/ptaa104
  * (61) C. J. Monahan, H. Na, C. M. Bouchard, G. P. Lepage and J. Shigemitsu, Phys. Rev. D 95, no.11, 114506 (2017) doi:10.1103/PhysRevD.95.114506 [arXiv:1703.09728 [hep-lat]].
  * (62) R. Aaij et al. [LHCb], Phys. Rev. Lett. 118, no.5, 052002 (2017) [erratum: Phys. Rev. Lett. 119, no.16, 169901 (2017)] doi:10.1103/PhysRevLett.118.052002 [arXiv:1612.05140 [hep-ex]].
  * (63) R. Aaij et al. [LHCb], Phys. Rev. D 100, no.3, 031102 (2019) doi:10.1103/PhysRevD.100.031102 [arXiv:1902.06794 [hep-ex]].
  * (64) R. Aaij et al. [LHCb], Phys. Rev. D 87, no.11, 112009 (2013) doi:10.1103/PhysRevD.87.112009 [arXiv:1304.6317 [hep-ex]].
  * (65) R. Aaij et al. [LHCb], JHEP 09, 179 (2015) doi:10.1007/JHEP09(2015)179 [arXiv:1506.08777 [hep-ex]].
  * (66) R. Aaij et al. [LHCb], Phys. Rev. Lett. 126, no.8, 081804 (2021) doi:10.1103/PhysRevLett.126.081804 [arXiv:2012.05143 [hep-ex]].
  * (67) R. Aaij et al. [LHCb], Phys. Rev. Lett. 112, no.13, 131802 (2014) doi:10.1103/PhysRevLett.112.131802 [arXiv:1401.5361 [hep-ex]].
  * (68) R. Aaij et al. [LHCb], Eur. Phys. J. C 77, no.4, 224 (2017) doi:10.1140/epjc/s10052-017-4744-6 [arXiv:1612.00945 [hep-ex]].
  * (69) F. F. Deppisch, P. S. Bhupal Dev and A. Pilaftsis, New J. Phys. 17, no.7, 075019 (2015) doi:10.1088/1367-2630/17/7/075019 [arXiv:1502.06541 [hep-ph]].
  * (70) A. Vaitaitis et al. [NuTeV and E815], Phys. Rev. Lett. 83, 4943-4946 (1999) doi:10.1103/PhysRevLett.83.4943 [arXiv:hep-ex/9908011 [hep-ex]].
  * (71) A. M. Cooper-Sarkar et al. [WA66], Phys. Lett. B 160, 207-211 (1985) doi:10.1016/0370-2693(85)91493-5
  * (72) D. Liventsev et al. [Belle], Phys. Rev. D 87, no.7, 071102 (2013) [erratum: Phys. Rev. D 95, no.9, 099903 (2017)] doi:10.1103/PhysRevD.87.071102 [arXiv:1301.1105 [hep-ex]].
  * (73) P. Abreu et al. [DELPHI], Z. Phys. C 74, 57-71 (1997) [erratum: Z. Phys. C 75, 580 (1997)] doi:10.1007/s002880050370
  * (74) G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873-3889 (1998) doi:10.1103/PhysRevD.57.3873 [arXiv:physics/9711021 [physics.data-an]].
  * (75) K. Bondarenko, A. Boyarsky, D. Gorbunov and O. Ruchayskiy, JHEP 11, 032 (2018) doi:10.1007/JHEP11(2018)032 [arXiv:1805.08567 [hep-ph]].


11institutetext: Inner Mongolia Normal University, 81 Zhaowuda Road, Saihan
District, Hohhot, Inner Mongolia, China
11email: {cisha1573,limeng20131015}<EMAIL_ADDRESS>
# LSTM-RPA: A Simple but Effective Long Sequence Prediction Algorithm for
Music Popularity Prediction††thanks: Supported by the National Natural Science
Foundation of China under Grant Nos. 61806103, 61562068; ;Inner Mongolia
discipline inspection and supervision big data laboratory open project fund
No. imdbd2020013;Science and technology planning project of Inner Mongolia
Autonomous Region No. jh20180175;Science and technology research project of
colleges and universities in Inner Mongolia Autonomous Region No. njzy21578,
No. njzy21551; the Youth Innovation and Entrepreneurship Talents of Inner
Mongolia “grassland talents” Project; the Graduate Innovation Fund of Inner
Mongolia Normal University under Grant No. CXJJS18112; ††thanks: This work has
been published on CNKI and the paper has been published in _Journal of
Computer Engineering and Applications_.
https://kns.cnki.net/kcms/detail/11.2127.TP.20210915.1447.010.html

Kun Li 11 Meng Li 22 Yanling Li 33 Min Lin 44

###### Abstract

The big data about music history contains information about time and users’
behavior. Researchers could predict the trend of popular songs accurately by
analyzing this data. The traditional trend prediction models can better
predict the short trend than the long trend. In this paper, we proposed the
improved LSTM Rolling Prediction Algorithm (LSTM-RPA), which combines LSTM
historical input with current prediction results as model input for next time
prediction. Meanwhile, this algorithm converts the long trend prediction task
into multiple short trend prediction tasks. The evaluation results show that
the LSTM-RPA model increased F score by 13.03%, 16.74%, 11.91%, 18.52%,
compared with LSTM, BiLSTM, GRU and RNN. And our method outperforms tradi-
tional sequence models, which are ARIMA and SMA, by 10.67% and 3.43%
improvement in F score.Code:https://github.com/maliaosaide/lstm-rpa

###### Keywords:

Long trend prediction LSTM Music popularity prediction Big data

##  1 Introduction

With the development of network technology, music as one of the streaming
media has also been developed rapidly, such as QQ Music 111https://y.qq.com/ ,
Spotify222https://www.spotify.com/ , Ali music 333https://www.xiami.com/ and
NetEase Cloud Music 444https://music.163.com/. People could listen, share,
download, collect music, and communicate their music ideas online. At the same
time, a large number of users’ behavior data are generated. The data from
different users have implicit relationships, although they seem like
irrelevant. And it also attracted the attention of many researchers and
Internet companies. With the help of music popularity prediction, enterprises
can formulate accurate propaganda strategies. The researchers recommend songs
for users accurately by analyzing the active users’ music play sequence[1],
musical complexity (chroma, rhythm, timbre, and arousal)[2], and Twitter
users’ music-listening behaviors[3]. Meanwhile, some studies analyzed the
music semantics constructs (genre, mood, instrumental, theme) [4], the music
frequency spectrogram image [5], and the music raw audio signal [5] to predict
the popularity of songs. We could predict the trend of music popularity
precisely from the analysis of different aspects of music.

In this paper, we investigated the relationship between the historical
information of audiences’ behavior and the trend of music popularity, and
proposed an algorithm to improve the accuracy of long trend prediction. This
historical information comes from the competition dataset of ” Popular Music
Prediction”555DataSet:
https://tianchi.aliyun.com/competition/entrance/231531/information-hosted by
Tsinghua University and Alicloud Tianchi Big Data Platform. The goal of the
contest is to predict music trends for the next month by analyzing this
historical data. This competition attracted 5,475 teams all over the world,
and this topic is of great research significance. We proposed the LSTM Rolling
Prediction Algorithm (LSTM-RPA) model with different features by analyzing the
time series of the data. This research could play an important role in
business decision making. We believe that using historical behavioral
information of audiences could improve the performance of trend prediction.

##  2 Related Work

From movie popularity prediction[6], stock market trend prediction[7], to
short-trend traffic prediction[8], trend prediction and time series modeling
have a long research history. Research [9] used the Auto Regressive Integrated
Moving Average (ARIMA) model to predict the Dez dam reservoir inflow monthly,
and their experimental results show that the Mean Square Error (MSE) of ARIMA
is smaller than that of the Auto Regressive Moving Average (ARMA). In[10], the
research found that the ARIMA model cannot predict the non-linear time series
accurately. So, they proposed the hybrid model based on Artificial Neural
Network (ANN) and Genetic Programming (GP), which predicted more accurately
than ARIMA on non-linear time series.

Research [11] established the Artificial “Music Market” and found a successful
song was not only related to the quality of the song itself but also to social
preferences, which increased the difficulty of prediction. In[5], the
researchers used the frequency spectrogram images of music raw audio signals
as the CNN model input feature to predict popularity music, and the accuracy
of this model is 61% in the test set. Research [12] established HitMusicNet
based on multimodal end-to-end Deep Learning (DL) architecture. HitMusicNet
used autoencoder to compress high-dimensional features, which consist of
audio, text lyrics, and meta-data to improve the accuracy of the model
prediction. Research [13] divided the time features of music play-count into a
basic trend and incremental trend, then they proposed Time Series based Music
Prediction (TSMP) which is based on category optimal value selection. When
songs’ play-count increased exponentially, the accuracy of the TSMP model
prediction would be reduced. Therefore, they proposed an Extend-Time Series
Music Prediction (E-TSMP) algorithm, which is combining Sub-Sequence Pattern
Matching Method (SSPMM) and Additional Processing (AD), to improve the
prediction accuracy.

In DL, Recurrent Neural Network (RNN) could learn key information very well
from sequences [14]. The Long Short-Term Memory (LSTM) network, the
improvement of RNN, can solve the problem of RNN gradient explosion and
gradient disappearance effectively in long sequences, and can transmit key
information to further cells [15]. Thus, we used LSTM to predict songs’ play-
count. We proposed a Rolling Prediction Algorithm (RPA) to solve the problem
of low accuracy in the long-trend prediction task. In this paper, we build two
kinds of models, that are traditional sequence models (ARIMA and Simple Moving
Average (SMA)[16]) and DL sequence models (LSTM, Bidirectional Long-Short Term
Memory (BiLSTM) [17], Gated Recurrent Unit (GRU) [18] and RNN model).

##  3 Methodology

###  3.1 LSTM model

Assume a time series
$X={x_{1},x_{2},x_{3},\dots,x_{n}},x_{1},x_{2},x_{3},\dots,x_{n}$ represent
each point in time that contains the feature information. (1)$\sim$(6) shows
the LSTM formula:

$\displaystyle f_{t}=\sigma(W_{f}\cdot h_{t-1}+W_{f}\cdot x_{t}+b_{f})$ (1)
$\displaystyle i_{t}=\sigma(W_{i}\cdot h_{t-1}+W_{i}\cdot x_{t}+b_{i})$ (2)
$\displaystyle\tilde{C}=\tanh(W_{c}\cdot h_{t-1}+W_{c}\cdot x_{t}+b_{c})$ (3)
$\displaystyle C_{t}=f_{t}\cdot C_{t-1}+i_{t}\cdot\tilde{C}$ (4)
$\displaystyle o_{t}=\sigma(W_{o}\cdot h_{t-1}+W_{o}\cdot x_{t}+b_{o})$ (5)
$\displaystyle h_{t}=o_{t}\cdot\tanh(C_{t})$ (6)

In LSTM block,$\sigma$ is the sigmoid function, $\tilde{C}$ is the new memory
of the current layer. In time t, $C_{t}$ is the memory cell block, $h_{t}$ is
the hidden state.

###  3.2 Single Feature LSTM Model

In the music popularity prediction task, one of the most important features is
the songs’ historical play-count. For this single feature, we built the Single
Feature LSTM Model (SF-LSTM) as shown in Fig. 1. This model inputs time series
data of different time steps, then these data respective go through two LSTM
layers which consist of 64 and 32 neurons respectively. We use ReLU as an
activation function to reduce information loss. After that, The SF-LSTM uses a
full-connected layer to reduce dimension, and outputs the prediction results.
$p$ is the input time step, $q$ is the output time step, and $H_{init}$ is the
initial hidden state.

Figure 1: Single Feature LSTM Model.

###  3.3 Multiple Features LSTM Model

Besides songs’ historical play-count, there are many other important features,
such as weekly play-count, monthly play-count, download-count, collect-count,
etc. We build the Multiple Features LSTM Model (MF-LSTM) as shown in Fig. 2,
which consists of three LSTM layers. The first two layers are also LSTM layers
that consist of 64 and 32 neurons respectively, ReLU as the activation
function. But the last layer is the LSTM layer that consisted of q neurons
with output step, sigmoid as the activation function. The particularity of
model is that when the input step and the output step are equal, the input and
output dimensions are same.

Figure 2: Multiple Features LSTM Model.

###  3.4 Rolling Prediction Algorithm (RPA)

Most of trend prediction algorithms only predict once[20]. That is, if we
predict the next five days, only change the output layer into five neurons. In
the case of mini-batch sets or lacking features, the model cannot get accuracy
highly in long-trend prediction task. So, we proposed the Rolling Prediction
Algorithm (RPA) which is based on moving average and Shortcut ideas [21] to
solve this problem. This algorithm optimizes the long trend prediction for
sequence models with mini-batch data sets or lacking features, and turn the
long trend prediction task into multiple short trend prediction tasks.

The RPA considers transforming a sequence with time $T$ step into multiple
prediction tasks with p length. When the RPA construct the input of the model
in predicting time $t+1$, it brings in the time t-1 and t information
together, so that the model can consider the information of $t-1$ and $t$ in
predicting time $t+1$. (7)$\sim$(8) shows algorithm formula:

$\displaystyle
X_{t}^{*}=X_{t-1}[x_{1}^{{}^{\prime}},x_{2}^{{}^{\prime}},\ldots,x_{l}^{{}^{\prime}}]+Y_{t}[y_{1},y_{2},\ldots,y_{p-l}]$
(7) $\displaystyle Y_{t+1}^{*}=unit(X_{t}^{*})$ (8)

Set the model input sequence $X$, length $p$; output sequence $Y$, length $q$;
and $l$ is rolling step, so $0\leq l\leq p,0\leq p-l\leq q$. In (7), $X_{t-1}$
is the input sequence for the model to predict $Y$ at time $t$, we reversely
order select $l$ input values:
$x_{1}^{{}^{\prime}},x_{2}^{{}^{\prime}},\ldots,x_{l}^{{}^{\prime}}$. $Y_{t}$
is the model prediction result at the time $t$, we order select $p-l$
prediction values: $y_{1},y_{2},\ldots,y_{p-l}$. Above all, at the prediction
$t+1$ time, the model input sequence $X_{t}^{*}$ consists of $X_{t-1}$ and
$Y_{t}$. In (8), $unit$ is the trained sequence neural network, and the
combination sequence $X_{t}^{*}$ gets the prediction result $Y_{t+1}^{*}$ by
$unit$. We build the LSTM-RPA, the network structure is shown as Fig. 3.

Figure 3: LSTM-RPA model. We use SF-LSTM and MF-LSTM as $unit$ to make up of
the LSTM-RPA model. $x_{t},x_{t-1},\ldots,x_{t-p}$ is input for time $t-1$,
and $y_{t+1},y_{t+2},\ldots,Y_{t+q}$ is prediction result for time $t$. At
time $t+1$, The model prediction input $X_{t}^{*}$ consist of
$X_{t-1}[x_{1}^{{}^{\prime}},x_{2}^{{}^{\prime}},\ldots,x_{l}^{{}^{\prime}}]$
and $Y_{t}[y_{1},y_{2},\ldots,y_{p-l}]$ sequences, then the predicted result
$Y_{t}^{*}[Y_{t+1}^{*},Y_{t+2}^{*},\ldots,Y_{t+n}^{*}]$ is acquired by the
trained LSTM model.

##  4 Experiments

###  4.1 Data preparation and analysis

The dataset includes 183 days, 50 artists, 349946 users and 10842 songs,
collected from March 1 to August 30, 2015. The sensitive information (Song_id,
Artist_id, and User_id) has been replaced by the fixed-length string, and one
song belongs to only one artist. The data format shows in Table 1, 2.

Table 1: User action | Name | Type | Description
---|---|---|---
| Song_id | String | Track name
| Artist_id | String | Artists’ name
| Ds | String | Date of data recording
| Gmt_create | String | Time of user’s playing
| Action_type | String | Type of songs’ action:1. play; 2. download; 3. collect
Table 2: Relations between artists and songs | Name | Type | Description
---|---|---|---
| Song_id | String | Track name
| Artist_id | String | Artists’ name
| Publish_time | String | Publish time
| Gender | String | Gender: 1. Male; 2. Female; 3. Combination
| Language | String | Language:1. Chinese; 2. French; 3. English

Table 1 shows users’ behavior, such as Track name, Artist, Date of data
recording, Type of songs’ action. Table 2 shows the relationship between the
artists and the songs, such as the song belongs to which one artist, the
song’s language, and the initial play-count. After preprocessing the data set,
we found that one artist’s data was lost seriously. So we selected the other
49 artists’ data as the experimental data.

###  4.2 Feature analysis and data slicing

Considering the influence of the different features on model accuracy, we
select daily songs’ play-count as the input feature in a single feature
experiment, and select daily artist songs’ play-count, download-count and
collect-count as the input feature in multiple features experiment.

At the same time, we chose the March to July dataset as the train and dev set
divided by 122:31, and the August data as the test set. The experiments’ goal
is to use the data from the first five months to predict daily songs’ play-
count in August.

(a) A (b) B (c) C (d) D

Figure 4: From March to July, the trend of songs playback by artists A, B, C,
D. The x-axis denotes the date, its range is 1-122; the y-axis denotes the
play-count

###  4.3 Music popularity trend analysis

Fig. 4 shows the daily play-count trend of 4 artists, which selected randomly
from 49 artists from March to July. The x-axis denotes the date, its range is
1-122. The y-axis denotes the play-count. Artist A’s play-count fluctuated a
lot in the 1 to 40 days, but after 40 days the play-count gradually stabilized
around 1000 daily. It is speculated that before 40 days this fluctuation may
be the recommended period after the artist releases a new album, or the
season’s top TV series, or movie songs produced by the artist. Artist B’s
total play-count trended to jitter upwards, with a rapid growth period of
nearly 10 days after 30 days. This was possibly influenced by external
platform recommendations or similar song recommendations. After 80 days, the
play-count trended to fluctuate slowly. Artist C had a steady trend in 1 70
days, but on the 75th day there have huge fluctuations; we guess: (1) Error
data. (2) Special holidays impact. (3) Influenced by sensational entertainment
news. Artist D’s daily play-count is extremely unstable.

Each artist’s historical play-count data contains their unique time
information, and there is a large gap between the different artists’ data. It
is impossible to use one model to predict all artist songs’ play-count trend,
so we build LSTM-RPA for 49 artists independently.

###  4.4 Evaluation

The competition organizer gives the evaluation function: F score, which is
used to measure the difference between the predicted value and the actual
value. a denotes artists $a\in W$, and $W$ is artists set. $Y_{a,d}$ is actual
play-count on $d$ day, and $X_{a,d}$ is model prediction value on $d$ day.
$\sigma_{a}$ is normalized variance $a$:

$\displaystyle\sigma_{a}=\sqrt{\frac{1}{N}\Sigma_{d=1}^{N}(\frac{X_{a,d}-Y_{a,d}}{Y_{a,d}})^{2}}$
(9)

In (9), N is the total number of days predicted. $\sigma_{a}$ indicate the gap
between the $X_{a,d}$ and $Y_{a,d}$. If $\sigma_{a}$ is small,
$(1-\sigma_{a})$ is large which means the model’s prediction very accurate.
$\phi_{a}$ is artist weight.

$\displaystyle\phi_{a}=\sqrt{\Sigma_{d=1}^{N}Y_{a,d}}$ (10)

Finally, the F score is defined.

$\displaystyle F=\Sigma_{a\in W}(1-\sigma_{a})\cdot\phi_{a}$ (11)

In this paper, we use the F score as a measure of prediction accuracy. The
larger the F score, the more accurate the prediction.

###  4.5 Experimental result

The music trend prediction is not only influenced by different features, but
also by different rolling steps and time steps. When the input is single or
multiple features, we researched the effect of different rolling and time
steps on the accuracy of long trend prediction, and it is based on the Keras
framework.

###  4.6 Single feature LSTM-RPA experiment

This experiment established and compared four kinds of neural networks: LSTM,
BiLSTM, GRU, and RNN. The GRU and RNN network structure is the same as the SF-
LSTM model. The BiLSTM model modifies the first two LSTM layers to
bidirectional LSTM layers, which means data flows forward and backward only in
the current bidirectional LSTM layer. The baseline result (black line) is the
prediction score of the corresponding model without using the RPA. In this
experiment, we only consider a positive score and change the negative score to
0.

Figure 5: F scores of different models with different time steps and rolling
steps in a single feature experiment. The x-axis denotes the F score, the
y-axis denotes the time step.

Fig. 5 shows the predicted F score of 49 artists by different rolling
prediction models under different time and rolling steps in the single
feature. When the time step is less than 3, the LSTM-RPA prediction accuracy
is better than the LSTM model (baseline). However, when the time step is
greater than 3, the LSTM-RPA model prediction accuracy is greatly reduced due
to the overfitting. In first four time steps, the BiLSTM-RPA model has an
underfitting (F score is 276.9514) at the time step 2, and other prediction
scores are higher than the baseline score. At the time step 4, different
rolling steps have different effects for the model’s prediction accuracy. When
rolling step increase, the model’s accuracy is also improved. As a variant of
the LSTM-RPA, the GRU-RPA model also has a high prediction accuracy when the
time step is less than 3. However, different from the LSTM-RPA model, the GRU-
RPA model still has prediction accuracy when the time step is greater than 3.
More than the half time-step, the RNN-RPA model’s prediction score is greater
than the baseline score.

###  4.7 Multiple features LSTM-RPA experiment

Because songs’ play-count trend is not only influenced by historical data but
also influenced by songs’ download-count and collect-count. Based on the
single feature LSTM-RPA experiment, we modified four neural networks
structure. GRU and RNN network structures are same as the MF-LSTM. The output
layer of BiLSTM model is the LSTM layer and others are the bidirectional LSTM
layer. The black line is the baseline score.

Figure 6: F scores of different models with different time step and rolling
step in the multi-feature experiment. The x-axis denotes the F score, the
y-axis denotes the time step.

Fig. 6 shows the F score of different models with different time and rolling
steps in the multiple feature experiment. When increasing the features, the
LSTM-RPA model prediction accuracy is better than the LSTM model without RPA
at all times steps. When the time step is greater than 5, the BiLSTM-RPA model
prediction accuracy decreases with the increase of time step. This is because
that the amount of training data increases when the time step increases. At
the same time, the amount of forward and reverse training data increases in
the bidirectional network structure, which makes the model strengthen the
memory of some non-necessary data, and thus leads to the decrease of the
prediction accuracy. In the GRU-RPA model experiment, the RPA could optimize
the baseline which has best score in all model’s baselines. The prediction
accuracy of the RNN-RPA model is greatly influenced by different rolling steps
in the case of the same time step, but some prediction F score is still higher
than the baseline.

###  4.8 Experimental analysis

The experiment explored the effects of single and multiple features on the
RPA, then selected and compared each RPA model’s best F scores and
corresponding baseline under different features. Table 3 shows the best F
scores and baseline of each RPA model in the single feature experiment. The
LSTM-RPA model has the highest F score and the lowest average error among the
four models. And the F score of GRU-RPA and LSTM-RPA model are similar. The
RPA could improve the prediction accuracy of RNN model better than other
models. This experiment shows that the RPA has more than 10% optimization
effect for these sequence models in the single feature experiment.

Table 3: The F score and baseline of each rolling prediction model in the single feature experiment | Model using RPA | Best F score | Baseline score | Optimum ratio | Average error
---|---|---|---|---|---
| LSTM | 4366.81 | 3863.12 | 13.03% | 26.34
| BiLSTM | 4326.58 | 3740.34 | 15.65% | 27.19
| GRU | 4360.07 | 3902.00 | 11.73% | 26.49
| RNN | 4287.89 | 3684.44 | 16.37% | 27.96

Table 4 shows the best F scores and baseline of different rolling prediction
models in the multi-feature experiment. The LSTM-RPA, BiLSTM-RPA, and GRU-RPA
have little difference in the best F score. And compared with their baselines,
the RPA has the highest improvement effect on the LSTM-RPA model. But this
optimization effect less than the RNN-RPA model. In all models, the GRU-RPA
model has the highest F score and the lowest average error. These experimental
results show that the RPA still has more than a 4% optimization effect in the
multi-feature experiment.

Table 4: The F score and baseline of each rolling prediction model in the multi-feature experiment | Model using RPA | Best F score | Baseline score | Optimum ratio | Average error
---|---|---|---|---|---
| LSTM | 4273.16 | 3813.11 | 12.06% | 28.26
| BiLSTM | 4263.43 | 4027.17 | 5.89% | 28.46
| GRU | 4291.26 | 4117.23 | 4.22% | 27.89
| RNN | 4198.39 | 3685.90 |  12.81% | 30.60

Based on single-feature and multi-feature experiments, we selected the rolling
prediction model with the highest predictive accuracy: single-feature LSTM-
RPA, and compared with the baseline of different models and traditional trend
prediction algorithms: ARIMA and SMA. Table 5 shows the best F scores and
average errors of traditional models and sequence models. The experimental
results show that the LSTM-RPA model has a better F score than other models.
The LSTM-RPA compared with ARIMA and SMA, F score increases by 10.67% and
3.43%, average error decreases by 32.64% and 11.23%, average error variance
decreases by 368.730% and 35.29%. And it compared with LSTM, BiLSTM, GRU, and
RNN, F score increases by 13.03%, 16.74%, 11.91% and 18.52%, average error
decreases by 39.02%, 48.55%, 36.02% and 52.88%, average error variance
decreases by 183.69%, 252.83%, 175.23% and 225.79%.

Table 5: The F scores and mean errors of best models | Evaluation | ARIMA | SMA | LSTM-RPA | LSTM | BiLSTM | GRU | RNN
---|---|---|---|---|---|---|---|---
| F score | 3945.53 | 4221.99 | 4366.81 | 3863.12 | 3740.34 | 3902.00 | 3684.44
| Average error | 34.94 | 29.30 | 26.34 | 36.62 | 39.13 | 35.83 | 40.27
| Mean error variance | 45.42 | 13.11 | 9.69 | 27.49 | 34.19 | 26.67 | 31.57

In summary, the model using RPA has higher prediction accuracy and lower
average error than baseline in single-feature and multi-feature experiments.
In the proportion of improvement of algorithm optimization, the RPA
optimization effect on multi-feature experiments is lower than that on single-
feature experiments. At the same time, the LSTM-RPA model has higher
prediction accuracy than the traditional trend prediction models. And the
LSTM-RPA model has a lower average error variance than the LSTM model. All in
all, the RPA has some optimization effect in the long trend prediction task.

##  5 Conclusion

In this paper, we analyzed the data about information of artist and user’s
behavior, selected different features (play-count, download-count and collect-
count) for experiments. The goal of experiments is to predict artists’ play-
count daily for 30 days. We proposed the LSTM Rolling Prediction Algorithm
(LSTM-RPA) based on the moving average and Shortcut idea to improve the
accuracy of the model in long trend prediction task, by analyzing the
relationship between the model’s historical input and current prediction
results. At time t+1, the LSTM-RPA model can consider the information of time
t-1 and t, and turn the long trend prediction task into multiple short trend
prediction tasks. Experiment results show that the prediction accuracy of the
single-feature LSTM-RPA model is better than the traditional model and DL
sequence models. The LSTM-RPA compared with ARIMA and SMA, the F score
increases by 10.67% and 3.43%, the average error decreases by 32.64% and
11.23%. And The LSTM-RPA increased F score by 13.03%, 16.74%, 11.91% and
18.52%, decreased average error by 39.02%, 48.55%, 36.02%, and 52.88%,
compared with LSTM, BiLSTM, GRU, and RNN. So, our method could improve
prediction accuracy of sequence model in long trend prediction task.

In the future, we will introduce the Self-Attention Mechanism to increase
model’s accuracy in long trend prediction task.

## References

  * [1] Cheng Z, Shen J, Zhu L, et al. Exploiting Music Play Sequence for Music Recommendation[C]//IJCAI. 2017, 17: 3654-3660.
  * [2] Lee J, Lee J S. Predicting music popularity patterns based on musical complexity and early stage popularity[C]//Proceedings of the Third Edition Workshop on Speech, Language & Audio in Multimedia. 2015: 3-6.
  * [3] Kim Y, Suh B, Lee K. # nowplaying the future Billboard: mining music listening behaviors of Twitter users for hit song prediction[C]//Proceedings of the first international workshop on Social media retrieval and analysis. 2014: 51-56.
  * [4] Ren J, KAUFFMAN R J. Understanding music track popularity in a social network[C]. AIS, 2017.
  * [5] Pleus B R M. Music Popularity Prediction via Techniques in Deep Supervised Learning[J].
  * [6] Latif M H, Afzal H. Prediction of movies popularity using machine learning techniques[J]. International Journal of Computer Science and Network Security (IJCSNS), 2016, 16(8): 127-131.
  * [7] Li Zhenzhen, Wu Qun. Research on stock forecasting algorithm based on LSTM neural network [J]. Fujian computer, 2019,35 (07): 41-43.
  * [8] Song Yujia, Zhang Jian, Xing Bin. Building short-term traffic prediction model based on long-term and short-term memory network [J]. Highway, 2019,64(07): 224-229.
  * [9] Mohammad Valipour, Mohammad Ebrahim Banihabib, Seyyed Mahmood Reza Behbahani. Parameters Estimate of Autoregressive Moving Average and Autoregressive Integrated Moving Average Models and Compare Their Ability for Inflow Forecasting[J]. Journal of Mathematics and Statistics,2012,8(3).
  * [10] Yi-Shian Lee, Lee-Ing Tong. Forecasting time series using a methodology based on autoregressive integrated moving average and genetic programming[J]. Knowledge-Based Systems,2010,24(1).
  * [11] Salganik, Matthew J., Peter Sheridan Dodds, and Duncan J. Watts. ”Experimental study of inequality and unpredictability in an artificial cultural market.” science 311.5762 (2006): 854-856.
  * [12] Martín-Gutiérrez D, Peñaloza G H, Belmonte-Hernández A, et al. A Multimodal End-to-End Deep Learning Architecture for Music Popularity Prediction[J]. IEEE Access, 2020, 8: 39361-39374.
  * [13] Yu Weisheng, Deng Weisheng, Zhang Yao, Li Shuyu. Prediction Study of Popular Trend of Music Based on Time Series [J]. Computer Engineering and Science, 2018,40 (09): 1703-1709.
  * [14] Mishra S, Rizoiu M A, Xie L. Modeling popularity in asynchronous social media streams with recurrent neural networks[J]. arXiv preprint arXiv:1804.02101, 2018.
  * [15] Youru Li, Zhenfeng Zhu, Deqiang Kong, Hua Han,Yao Zhao. EA-LSTM: Evolutionary attention-based LSTM for time series prediction[J]. Knowledge-Based Systems,2019,181.
  * [16] Lauren S , Harlili S D . Stock trend prediction using simple moving average supported by news classification[C]// International Conference of Advanced Informatics: Concept. IEEE, 2014.
  * [17] Xu Xianfeng, Liu ahui, Chen Yulu, Cai lulu. Short term power prediction of bilstm photovoltaic power generation based on meteorological factors [J]. Computer system application, 2020,29 (07): 205-211.
  * [18] Zhang Jinlei, Luo Yuling, Fu Qiang. Financial time series prediction based on gated cyclic unit neural network [J]. Journal of Guangxi Normal University (NATURAL SCIENCE EDITION), 2019,37 (02): 82-89.
  * [19] Hochreiter S, Schmidhuber J. Long short-term memory [J]. Neural Computation, 1997, 9(8): 1735-1780.
  * [20] Yu H, Li Y, Zhang S, et al. Popularity Prediction for Artists Based on User Songs Dataset[C]//Proceedings of the 2019 5th International Conference on Computing and Artificial Intelligence. 2019: 17-24.
  * [21] Wang Yuandong. Image Classification Method and Application Based on ResNet Model [D]. East China Jiaotong University, 2019.