KIAS-P21043, CTP-SCU/2021031, APCTP Pre2021-015

#  Muon $g-2$, $B\to K^{(*)}\mu^{+}\mu^{-}$ anomalies, and leptophilic dark
matter in $U(1)_{\mu-\tau}$ gauge symmetry

P. Ko<EMAIL_ADDRESS>School of Physics, KIAS, Seoul 02455, Korea  Quantum
Universe Center, KIAS, Seoul 02455, Korea  Takaaki Nomura<EMAIL_ADDRESS>College of Physics, Sichuan University, Chengdu 610065, China  Hiroshi Okada
<EMAIL_ADDRESS>Asia Pacific Center for Theoretical Physics, Pohang
37673, Republic of Korea  Department of Physics, Pohang University of Science
and Technology, Pohang 37673, Republic of Korea

###### Abstract

We propose a new class of $U(1)_{\mu-\tau}$ gauged model that can explain
recent flavor anomalies such as muon $g-2$, $b\to s\mu^{+}\mu^{-}$, as well as
include a scalar dark matter candidate while satisfying all the
phenomenological constraints. For this purpose, we add new vectorlike quarks
and leptons as well as two inert singlet scalars, one leptophilic and the
other leptophobic. In our model $b\to s\mu^{+}\mu^{-}$ anomalies can be
explained by loop induced interactions among the SM and exotic quarks and
$Z^{\prime}$ gauge boson. In our numerical analysis, we show allowed region of
our model to be narrow, and it would be tested soon, for example by searching
for 4 muons and dimuon + missing transverse momentum from
$pp\to\mu^{+}\mu^{-}Z^{\prime}$ and
$pp\to\nu_{\mu,\tau}\bar{\nu}_{\mu,\tau}Z^{\prime}$ followed by
$Z^{\prime}\to\mu^{+}\mu^{-},\nu_{\mu,\tau}\bar{\nu}_{\mu,\tau}$ at the LHC.

##  I Introduction

Even after discovery of the standard model (SM) Higgs boson, there are a
number of reasons for physics beyond the SM (BSM), such as Baryon Asymmetry of
Universe (BAU), neutrino oscillations Esteban:2020cvm , nonbaryonic dark
matter Planck:2015fie , all of which are well established. Also there are some
anomalies in the flavor sector: the ratio between rare $B$ meson decays of
$B^{+}\to K^{+}\mu^{+}\mu^{-}$ and $B^{+}\to K^{+}e\bar{e}$ ($R_{K}$)
Hiller:2003js ; Bobeth:2007dw ; Aaij:2014ora ; Aaij:2019wad , and muon
anomalous magnetic moment (muon $g-2$) Bennett:2006fi ; Abi:2021gix . These
anomalies may call for explanations within some new physics beyond the SM
(BSM), although some anomalies may (partly) originate from poorly understood
SM predictions involving nonperturbative QCD.

Especially, muon $g-2$ is known as longstanding anomaly that was initially
discovered by E821 experiment at Brookhaven National Lab (BNL) two decades ago
Bennett:2006fi , and has recently been confirmed by E989 Run 1 experiment at
Fermilab(FNAL) Abi:2021gix ; combining these measurements we obtain

$\displaystyle\Delta a_{\mu}=(25.1\pm 5.9)\times 10^{-10},$ (1)

whose deviation from SM prediction Aoyama:2012wk ; Aoyama:2019ryr ;
Czarnecki:2002nt ; Gnendiger:2013pva ; Davier:2017zfy ; Keshavarzi:2018mgv ;
Colangelo:2018mtw ; Hoferichter:2019mqg ; Davier:2019can ; Keshavarzi:2019abf
; Kurz:2014wya ; Melnikov:2003xd ; Masjuan:2017tvw ; Colangelo:2017fiz ;
Hoferichter:2018kwz ; Gerardin:2019vio ; Bijnens:2019ghy ; Colangelo:2019uex ;
Blum:2019ugy ; Colangelo:2014qya ; Hagiwara:2011af  is 4.2 $\sigma$.

Also $R_{K}$ anomaly has been recently updated by the LHCb collaboration
Aaij:2021vac :

$\displaystyle R_{K}=0.846^{+0.042+0.013}_{-0.039-0.012}\quad(1.1{\rm
GeV^{2}}<q^{2}<6{\rm GeV^{2}}),$ (2)

where the first (second) uncertainty is statical (systematic) one and $q^{2}$
is the dilepton invariant mass squared. Its deviation from the SM prediction
is about 3.1 $\sigma$. There are also discrepancies in the measurements of the
angular observable $P^{\prime}_{5}$ in the $B$ meson decay ($B\to
K^{*}\mu^{+}\mu^{-}$) DescotesGenon:2012zf ; Aaij:2015oid ; Aaij:2013qta ;
Abdesselam:2016llu ; Wehle:2016yoi , as well as in $R_{K^{*}}=BR(B\to
K^{*}\mu^{+}\mu^{-})/BR(B\to K^{*}e^{+}e^{-})$ Aaij:2017vbb ; we call these as
$b\to s\mu^{+}\mu^{-}$ anomalies.

$U(1)_{\mu-\tau}$ gauge symmetry is anomaly free even without additional
chiral fermions beyond the SM fermions He:1990pn ; He:1991qd . It has been
frequently applied to models that explains some of these flavor anomalies. Due
to feature of flavor dependent symmetry, we often obtain flavor specific
signals that can be more testable through various experiments than the SM.
Moreover, since a new gauge boson does not directly interact with electrons up
to small knetic mixing, we can evade stringent constraints on the gauge
coupling and gauge boson mass from LEP (or ILC) easier than the other gauged
models. 111Notice here that we have to consider a constraint of neutrino
trident production in gauged $U(1)_{\mu-\tau}$ models Altmannshofer:2014pba ,
and we will discuss it in the main text.  Also the muon $g-2$ within the
$U(1)_{\mu-\tau}$ gauge models was studied right after the BNL announcement
Baek:2001kca . And the PAMELA $e^{+}$ excess was studied in the
$U(1)_{\mu-\tau}$-charged dark matter with mass $\sim O(1)$ TeV in Ref.
Baek:2008nz .

In fact, there is a wide variety of recent applications of $U(1)_{\mu-\tau}$
models, which are distinguished by additional particle contents to address
flavor anomalies; Refs Baek:2015mna ; Nomura:2018cle ; Nomura:2018vfz ;
Asai:2018ocx ; Asai:2017ryy ; Lee:2017ekw  for neutrino oscillations to
control lepton flavor structure, Refs Davier:2010nc ; Davier:2017zfy ;
Davier:2019can ; Borah:2021mri ; Qi:2021rhh ; Singirala:2021gok ;
Buras:2021btx ; Zhou:2021vnf ; Borah:2021jzu ; Chen:2021vzk ; Zu:2021odn ;
Huang:2021nkl ; Patra:2016shz ; Altmannshofer:2016oaq  for muon $g-2$, Refs
Altmannshofer:2014cfa ; Crivellin:2015mga ; Crivellin:2015lwa ; Ko:2017yrd ;
Kumar:2020web ; Han:2019diw ; Chao:2021qxq ; Baek:2019qte ; Chen:2017usq ;
Borah:2021khc  for $b\to s\mu^{+}\mu^{-}$ anomalies, and (or) DM. Qualitative
features in $U(1)_{\mu-\tau}$ models for flavor anomalies are the following.
The muon $g-2$ anomaly is explained by light $Z^{\prime}$ around
$\mathcal{O}(10)$ to $\mathcal{O}(100)$ MeV scale or by Yukawa interactions
with exotic fermions controlled by $U(1)_{\mu-\tau}$ symmetry. For explaining
$b\to s\mu^{+}\mu^{-}$ anomalies with $U(1)_{\mu-\tau}$, generally we need to
add new particles since the corresponding gauge boson $Z^{\prime}$ does not
couple to quarks. Mixing among the SM quarks and new vector-like quarks are
considered to generate flavor violating $Z^{\prime}$-quark coupling in refs
Crivellin:2015mga ; Altmannshofer:2014cfa ; Han:2019diw  while flavor
dependent charge assignment in the SM quark sector is considered in ref.
Crivellin:2015lwa .

In ref. Ko:2017yrd , we proposed a gauged $U(1)_{\mu-\tau}$ model in which
interactions among $Z^{\prime}$ and quarks are radiatively generated via one-
loop diagram with vector-like quarks($Q^{\prime}$) and scalar dark matter(DM)
candidate propagating inside loop 222See also refs. Kumar:2020web ;
Chao:2021qxq ; Baek:2019qte ; Chen:2017usq  for similar mechanism and
analysis. Also scenarios in which vector-like fermions are mixed with the SM
fermions are also found e.g. in refs. Altmannshofer:2014cfa ; Navarro:2021sfb
.. In this work, we focused on a loop diagram associated with
$Z^{\prime}Q^{\prime}\bar{Q}^{\prime}$ interaction. But there are
contributions associated with $ZQ^{\prime}\bar{Q}^{\prime}$ and
photon$-Q^{\prime}\bar{Q}^{\prime}$ interactions, since vector-like quarks
interact with $Z$ and photon that should be taken into account to check
consistency of the scenario. Note also that we could not explain muon $g-2$ in
this scenario since $Z^{\prime}$ from $U(1)_{\mu-\tau}$ should be heavier than
$\mathcal{O}(10)$ GeV. Also in Ref. Ko:2017yrd  we introduced only vector-like
quarks and not vector-like leptons, and the model looks non symmetric. Thus it
is interesting to add vector-like leptons in addition to vector-like quarks in
the previous scenario to make the model look more symmetrical, to explain muon
$g-2$ and consider loop diagrams associated with $Z$ and photon interactions
to check consistency for $b\rightarrow s\mu^{+}\mu^{-}$ anomalies. We will
also introduce two inert complex scalars $\chi_{q}$ (leptophobic) and
$\chi_{l}$ (leptophilic) that are charged under $U(1)_{\mu-\tau}$ in such a
way that Yukawa couplings among the vector-like quarks (leptons), SM quarks
(leptons) and the inert scalar $\chi_{q}$ ($\chi_{l}$) are allowed by gauge
symmetries.

In the present work, we extend the previous model Ko:2017yrd  by introducing
vector-like leptons to explains muon $g-2$, $R_{K}$, and DM in a gauged
$U(1)_{\mu-\tau}$ symmetry, where DM is assumed to be a weak interacting
massive particle (WIMP). We then estimate one loop diagrams considering
$Z^{\prime}$, $Z$ and photon vertices with vector-like quarks that affects
$b\to s\mu^{+}\mu^{-}$ process to explain $R_{K}$ anomaly and to check
consistency of the model. We also estimate muon $g-2$ and DM relic density to
search for parameters explaining them. In addition we take all the valid
constraints such as neutral meson mixings and $B_{s(d)}^{0}\to\mu^{+}\mu^{-}$
into consideration and show the allowed region in our numerical analysis.

This letter is organized as follows. In Sec. II, we present our model and show
several constraints from $B\to K^{(*)}\ell^{+}\ell^{-}$, neutral meson
mixings, muon $g-2$, and (leptophilic) dark matter. In Sec.III we carry out
the numerical analysis and demonstrate our allowed region. Finally, Sec. IV is
devoted to the summary of our results and conclusion.

##  II Model setup and Constraints

|  $Q^{\prime}$ |  $\chi_{q}$ |  $L^{\prime}$ |  $E$ |  $\chi_{\ell}$ |  $\varphi$
---|---|---|---|---|---|---
$SU(3)_{C}$ | $\bm{3}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$
$SU(2)_{L}$ | $\bm{2}$ | $\bm{1}$ | $\bm{2}$ | $\bm{1}$ | $\bm{1}$ | $\bm{1}$
$U(1)_{Y}$ | $\frac{1}{6}$ | $0$ | $-\frac{1}{2}$ | $-1$ | $0$ | $0$
$U(1)_{\mu-\tau}$ | $q_{x}$ | $q_{x}$ | $\frac{q_{x}}{2}$ | $\frac{q_{x}}{2}$ | $\frac{q_{x}}{2}-1$ | ${q_{\varphi}}$
Table 1:  Charge assignments of the new fields $Q^{\prime}$,
$L^{\prime}\equiv[N^{\prime},E^{\prime}]^{T}$, $E$, $\chi_{q}$ and
$\chi_{\ell}$ under $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times
U(1)_{\mu-\tau}$. Here $Q^{\prime}$, $L^{\prime}$ and $E$ are vector-like
fermions. $\chi_{q}$ and $\chi_{\ell}$ are complex inert scalar bosons either
of which can be considered as a DM candidate. $\varphi$ plays a role in
breaking the $U(1)_{\mu-\tau}$ symmetry spontaneously by the nonzero VEV;
$\langle\varphi\rangle=v_{\varphi}/\sqrt{2}$.

Figure 1: The diagrams introducing effective coupling for
$Z^{\prime}_{\mu}(Z_{\mu})\bar{b}\gamma^{\mu}P_{L}s+h.c.$ and
$F^{\mu\nu}\bar{b}\sigma_{\mu\nu}P_{L,R}s+h.c.$ interaction.

In this section we define our model. We introduce a vector-like exotic quark
$Q^{\prime}\equiv[U^{\prime},D^{\prime}]^{T}$ and lepton
$L^{\prime}\equiv[N^{\prime},E^{\prime}]^{T}$ under isospin doublet and a
singly charged-lepton $E$ under isospin singlet, two complex scalar bosons
$\chi_{q}$ and $\chi_{\ell}$, and a boson $\varphi$, where
$\chi_{q}$($\chi_{\ell}$) couples only to quark(lepton) sector respectively.
The singlet scalar $\varphi$ plays a role in breaking the $U(1)_{\mu-\tau}$
symmetry spontaneously developing the nonzero vacuum expectation value (VEV);
$\langle\varphi\rangle=v_{\varphi}/\sqrt{2}$. Here we take the
$U(1)_{\mu-\tau}$ charge of $\varphi$ so as to retain the inert features of
$\chi_{q,\ell}$; we forbid terms $\chi_{l,q}^{(*)}(\varphi^{(*)})^{n}$ with
$n$ being integer. This can be achieved by imposing conditions for the charge
of $\varphi$ as follows:

$\displaystyle q_{\varphi}\neq\pm\frac{1}{n}\left(\frac{q_{x}}{2}-1\right),\
\pm\frac{1}{n}q_{x}.$ (3)

In addition we forbid these operators,
$\bar{Q}_{L}Q^{\prime}_{R}\varphi^{(*)}$,
$\bar{L}_{L_{\ell}}L^{\prime}_{L}\varphi^{(*)}$ and
$\bar{\ell}_{R}E_{L}\varphi^{(*)}$ $(\ell=e,\mu,\tau)$ and their Hermitian
conjugates, so that newly added vector-like fermions and the SM ones do not
mix. In fact, $\bar{Q}_{L}Q^{\prime}_{R}\varphi^{(*)}$,
$\bar{L}_{L_{e,\mu}}L^{\prime}_{L}\varphi^{(*)}$ and
$\bar{e}_{R}(\bar{\mu}_{R})E_{L}\varphi^{(*)}$ are already forbidden by the
conditions in Eq. (3). To forbid remaining operators we require additional
conditions,

$q_{\varphi}\neq\pm\left(\frac{q_{x}}{2}+1\right).$ (4)

The SM Higgs is denoted by $H$ and its VEV is given by
$[0,v/_{H}/\sqrt{2}]^{T}$. $Q^{\prime}$ as well as $\chi_{q}$ is requested to
induce sizable $B\to K^{(*)}\ell\bar{\ell}$ anomalies, $E$ and $L^{\prime}$ as
well as $\chi_{\ell}$ are required to get sizable lepton anomalous magnetic
moment. 333With $E$ only, one cannot get sizable lepton anomalous magnetic
moment, since the mass of $E$ has to be very small at about 200 GeV, which
would be ruled out by LHC. $L^{\prime}$ plays a role in evading chiral
suppression and thus sizable values are obtained. Charge assignments of these
new fields are summarized in Table 1. Then, the valid Lagrangian under these
symmetries is given by

$\displaystyle-\mathcal{L}_{\rm VLF+\chi}=$ $\displaystyle
M_{Q^{\prime}}\bar{Q}^{\prime}Q^{\prime}+M_{E}\bar{E}^{\prime}E^{\prime}+{M_{L^{\prime}}\bar{L}^{\prime}_{L}L^{\prime}_{R}}+y\overline{L^{\prime}_{L}}HE_{R}+m_{\chi_{q}}^{2}|\chi_{q}|^{2}+m_{\chi_{\ell}}^{2}|\chi_{\ell}|^{2}$
$\displaystyle+(f_{i}\overline{Q^{\prime}_{R}}Q_{L_{i}}\chi_{q}+g_{\mu}\overline{\mu_{R}}E_{L}\chi_{\ell}^{*}+h_{\mu}\overline{L_{L_{\mu}}}L^{\prime}_{R}\chi_{\ell}^{*}+{\rm
h.c.}),$ (5)

where $i=1-3$ are generation indices, $Q_{L}$’s are the SM quark doublets and
$\mu_{R}$ is the right-handed muon. We have abbreviated kinetic terms for
simplicity.

The scalar potential in our model is given by

$\displaystyle
V=-M_{H}^{2}H^{\dagger}H-M^{2}_{\varphi}\varphi^{*}\varphi+\lambda_{H}(H^{\dagger}H)^{2}+\lambda_{\varphi}(\varphi^{*}\varphi)^{2}+\lambda_{H\varphi}(H^{\dagger}H)(\varphi^{*}\varphi),$
(6)

where we have abbreviated scalar potential associated with $\chi_{q,\ell}$ for
simplicity. We expand scalar fields around its VEVs as

$\displaystyle H=\begin{pmatrix}w^{+}\\\
\frac{1}{\sqrt{2}}(v+\tilde{h}+iG_{Z})\end{pmatrix},\quad\varphi$
$\displaystyle=\frac{1}{\sqrt{2}}(v_{\varphi}+\tilde{\phi}+iG_{Z^{\prime}}),$
(7)

where $\tilde{h}$ and $\tilde{\phi}$ correspond to CP-even physical scalar
boson states while $w^{+}$, $G_{Z}$ and $G_{Z^{\prime}}$ are the Nambu-
Goldstone(NG) bosons which are absorbed by the weak gauge bosons $Z,~{}W$ and
$Z^{\prime}$. The mass matrix for physical scalar bosons is obtained as

$\displaystyle
M^{2}=\begin{pmatrix}\lambda_{H}v^{2}&\lambda_{H\varphi}vv_{\varphi}\\\
\lambda_{H\varphi}vv_{\varphi}&\lambda_{\varphi}v_{\varphi}^{2}\end{pmatrix}.$
(8)

The matrix is diagonalized by acting an orthogonal matrix $U$ such that

$\displaystyle U^{T}M^{2}U=\mathrm{diag}(m_{h}^{2},m_{H}^{2}),$ (9)

where

$\displaystyle U=\begin{pmatrix}\cos\alpha&\sin\alpha\\\
-\sin\alpha&\cos\alpha\end{pmatrix},$ (10)

and the scalar mixing angle $\alpha$ is given by

$\displaystyle\tan
2\alpha=\frac{2\lambda_{H\varphi}vv_{\varphi}}{\lambda_{H}v^{2}-\lambda_{\varphi}v_{\varphi}^{2}}.$
(11)

Then mass eigenstates are given by

$\displaystyle\begin{pmatrix}h\\\
\phi\end{pmatrix}=U^{T}\begin{pmatrix}\tilde{h}\\\
\tilde{\phi}\end{pmatrix}=\begin{pmatrix}\cos\alpha&-\sin\alpha\\\
\sin\alpha&\cos\alpha\end{pmatrix}\begin{pmatrix}\tilde{h}\\\
\tilde{\phi}\end{pmatrix},$ (12)

where $h$ corresponds to the SM-like Higgs boson.

Exotic singly-charged fermions $(E,E^{\prime})$ mix each other after
spontaneous electroweak symmetry breaking, and its mass matrix is given by

$\displaystyle{\cal M}_{E}=\left(\begin{array}[]{cc}M_{E}&m_{E}\\\
m_{E}&M_{L^{\prime}}\\\ \end{array}\right),$ (15)

where $m_{E}\equiv yv/\sqrt{2}$. It is diagonalized by a unitary matrix
$V_{E}$ as ${\rm diag}(M_{1},M_{2})=V_{E}^{\dagger}{\cal M}_{E}V_{E}$, and
$V_{E}$ is given by

$\displaystyle V_{E}=\left(\begin{array}[]{cc}c_{c}&-s_{c}\\\ s_{c}&c_{c}\\\
\end{array}\right),\quad\tan 2\theta_{c}=\frac{2m_{E}}{M_{E}-M_{L^{\prime}}},$
(18)

where $c_{c}(s_{c})$ is shorthand symbol of $\cos\theta_{c}(\sin\theta_{c})$.
We write mass eigenstates as $E_{1}$ and $E_{2}$ whose masses are $M_{1}$ and
$M_{2}$.

##  III Phenomenological Constraints

###  III.1 Wilson coefficients for $B\to K^{(*)}\ell^{+}\ell^{-}$ decay

Here we derive Wilson coefficients associated with $B\to
K^{(*)}\ell^{+}\ell^{-}$ processes derived from diagrams in Fig. 1. Anomalies
observed in $B\to K^{(*)}\ell^{+}\ell^{-}$ decays can be explained by the
shift of the Wilson coefficient $C_{9}$ which corresponds to the operator
$(\bar{s}\gamma_{\mu}P_{L}b)(\mu\gamma^{\mu}\mu)$. In our model, the effective
coupling for $Z^{\prime}_{\mu}\bar{b}\gamma^{\mu}P_{L}s+h.c.$ is induced at
one loop level as shown in Fig. 1 with the Yukawa coupling in Eq. (5). Then
the effective Hamiltonian
$(\bar{s}\gamma_{\mu}P_{L}b)(\bar{\mu}\gamma^{\mu}\mu)$ arises from
$Z^{\prime}$ mediation and the contribution to Wilson coefficient $\Delta
C_{9}^{\mu\mu}(Z^{\prime})$ is obtained as:

$\displaystyle\Delta C_{9}^{\mu\mu}(Z^{\prime})\simeq-\frac{q_{x}g^{\prime
2}}{4(4\pi)^{2}m_{Z^{\prime}}^{2}C_{\rm SM}}f^{*}_{3}f_{2},$
$\displaystyle{C_{\rm SM}\equiv\frac{V_{tb}V^{*}_{ts}G_{\rm F}\alpha_{\rm
em}}{\sqrt{2}\pi},}$ (19)

where $V_{tb}\approx 0.999$, $V_{ts}\approx-0.040$ are the 3-3 and 3-2
elements of CKM matrix respectively, $G_{\rm F}\approx 1.17\times 10^{-5}$ GeV
is the Fermi constant, $\alpha_{\rm em}\approx 1/137$ is the electromagnetic
fine-structure constant, $m_{b}\approx 4.18$ GeV and $m_{s}\approx{0.095}$ GeV
are respectively the bottom and strange quark masses given in the
$\overline{MS}$ scheme at a renormalization scale $\mu=2$ GeV pdg . Notice
here that we have ignored $m_{b}$ and $m_{s}$ in the formula for $C_{9}$ in
Eq. (19), since they are much smaller than $M_{Q^{\prime}}$ and $m_{\chi_{q}}$
in our scenario. The $\Delta C_{9}^{\mu\mu}$ allowed within 2$\sigma$ range is
given by Alguero:2021anc :

$\displaystyle 0.83\leq-\Delta C_{9}^{\mu\mu}\leq 1.52.$ (20)

We also obtain lepton universal Wilson coefficient $\Delta C_{9}(Z)$ from the
diagram where $Z^{\prime}$ is replaced to $Z$. Another lepton universal Wilson
coefficient $\Delta C_{10}(Z)$ associated with the corresponding operator
$(\bar{s}\gamma_{\mu}P_{L}b)(\ell\gamma^{\mu}\gamma_{5}\ell)$ is also arisen
from $Z$ mediation. The formulas of them are given by

$\displaystyle\Delta
C_{9}(Z)\simeq\frac{g_{2}^{2}f^{*}_{3}f_{2}}{4(4\pi)^{2}m_{Z}^{2}c_{w}^{2}C_{\rm
SM}}\left(-\frac{1}{2}+2s^{2}_{w}\right)\left(-\frac{1}{2}+\frac{1}{3}s^{2}_{w}\right),$
(21) $\displaystyle\Delta
C_{10}(Z)\simeq\frac{g_{2}^{2}f^{*}_{3}f_{2}}{8(4\pi)^{2}m_{Z}^{2}c_{w}^{2}C_{\rm
SM}}\left(-\frac{1}{2}+\frac{1}{3}s^{2}_{w}\right).$ (22)

The $\Delta C_{10}$ contributes to $B_{s}^{0}\to\mu^{+}\mu^{-}$ process. Here
we consider measurement of $B_{s}^{0}\to\mu^{+}\mu^{-}$ branching ratio
LHCb:2021vsc ; LHCb:2021awg

$BR(B_{s}^{0}\to\mu^{+}\mu^{-})^{\rm
exp}=(3.09^{+0.46+0.15}_{-0.43-0.11})\times 10^{-9},$ (23)

where the first and second uncertainties are statistical and systematic
respectively. We compare the experimental value with the theoretical one
estimated by Hiller:2014yaa

$BR(B_{s}^{0}\to\mu^{+}\mu^{-})^{\rm th}=|1-0.24\Delta
C_{10}^{\mu\mu}|^{2}BR(B_{s}^{0}\to\mu^{+}\mu^{-})^{\rm SM},$ (24)

where $BR(B_{s}^{0}\to\mu^{+}\mu^{-})^{\rm SM}=(3.65\pm 0.23)\times 10^{-9}$
is the SM prediction Bobeth:2013uxa . $\Delta C_{10}^{\mu\mu}$ is restricted
by the data of the process and $B\to K^{(*)}\ell^{+}\ell^{-}$. For example,
allowed range at 2$\sigma$ is given by Alguero:2021anc  including data of
$B_{s}^{0}\to\mu^{+}\mu^{-}$ decay: $-0.11\leq\Delta C_{10}^{\mu\mu}\leq
0.49$. In fact this range will be modified as we have lepton flavor universal
$\Delta C_{9}(Z)$ and $\Delta C_{10}(Z)$ effects. Here we impose $0<\Delta
C_{10}(Z)<0.79$ to satisfy measurement of $BR(B_{s}^{0}\to\mu^{+}\mu^{-})$; we
assume $\Delta C_{10}(Z)>0$ since it is preferred by data.

Furthermore the Wilson coefficient $C_{7}$ is modified by $\Delta C_{7}$:

$\displaystyle\Delta
C_{7}=\frac{f_{2}f^{*}_{3}}{2}\int[dX]^{3}\frac{xy}{x(x-1)m_{b}^{2}+xm^{2}_{\chi_{q}}+(y+z)M^{2}_{Q^{\prime}}}\left(\frac{4G_{F}}{\sqrt{2}}V_{tb}V^{*}_{ts}\right)^{-1},$
(25)

where $\int[dX]^{3}\equiv\int_{0}^{1}dxdydz\delta(1-x-y-z)$. In our numerical
analysis we apply the constraints Alguero:2021anc :

$\displaystyle-0.01\leq\Delta C_{7}\leq 0.03.$ (26)

Note that Wilson coefficient $\Delta C_{7^{\prime}}$ is also induced but it is
suppressed by $m_{s}/m_{b}$ factor compared to $C_{7}$.

###  III.2 $M-\overline{M}$ mixing

The exotic vector-like quarks and the complex scalar DM $\chi$ induce the
neutral meson ($M$)-antimeson ($\overline{M}$) mixings such as
$K^{0}-\bar{K}^{0}$, $B_{d}-\bar{B}_{d}$, $B_{s}-\bar{B}_{s}$, and
$D^{0}-\bar{D}^{0}$ from the box type one-loop diagrams. The formulae for the
mass splitting are respectively given by Gabbiani:1996hi

$\displaystyle\Delta m_{K}$
$\displaystyle\approx\sum_{a,b=1}^{3}|f_{1}f_{2}|^{2}G^{K}_{box}[m_{\chi_{q}},M_{Q^{\prime}}]\lesssim
3.48\times 10^{-15}\ [{\rm GeV}],$ (27) $\displaystyle\Delta m_{B_{d}}$
$\displaystyle\approx\sum_{a,b=1}^{3}|f_{1}f_{3}|^{2}G^{B_{d}}_{box}[m_{\chi_{q}},M_{Q^{\prime}}]\lesssim
3.36\times 10^{-13}\ [{\rm GeV}],$ (28) $\displaystyle\Delta m_{B_{s}}$
$\displaystyle\approx\sum_{a,b=1}^{3}|f_{2}f_{3}|^{2}G^{B_{s}}_{box}[m_{\chi_{q}},M_{Q^{\prime}}]\lesssim
1.17\times 10^{-11}\ [{\rm GeV}],$ (29) $\displaystyle\Delta m_{D}$
$\displaystyle\approx\sum_{a,b=1}^{3}|f_{1}f_{2}|^{2}G^{D}_{box}[m_{\chi_{q}},M_{Q^{\prime}}]\lesssim
6.25\times 10^{-15}\ [{\rm GeV}],$ (30) $\displaystyle
G^{M}_{box}(m_{1},m_{2})$
$\displaystyle=\frac{m_{M}f_{M}^{2}}{3(4\pi)^{2}}\int[dX]^{3}\frac{x}{xm^{2}_{1}+(y+z)m_{2}^{2}}$
$\displaystyle=\frac{m_{M}f_{M}^{2}}{6(4\pi)^{2}}\left[\frac{m_{1}^{4}-m_{2}^{4}+4m_{1}^{2}m_{2}^{2}\ln\left[\frac{m_{2}}{m_{1}}\right]}{(m_{1}^{2}-m_{2}^{2})^{3}}\right]\hskip
113.81102pt(m_{1}\neq m_{2}),$ (31)

where relevant quarks $(q,q^{\prime})$ are respectively $(d,s)$ for $K$,
$(b,d)$ for ${B_{d}}$, $(b,s)$ for ${B_{s}}$, and $(u,c)$ for $D$, each of the
last inequalities of the above equations represent the upper bound from the
experimental values pdg , and $f_{K}\approx 0.156$ GeV,
$f_{B_{d}(B_{s})}\approx 0.191(0.200)$ GeV, $f_{D}\approx 0.212$ GeV,
$m_{K}\approx 0.498$ GeV, $m_{B_{d}(B_{s})}\approx 5.280(5.367)$ GeV, and
$m_{D}\approx 1.865$ GeV.

###  III.3  Muon anomalous magnetic dipole moment

Recently, E989 Collaboration at Fermilab reported the new result on the muon
$(g-2)$ Abi:2021gix :

$\displaystyle a^{\rm FNAL}_{\mu}=116592040(54)\times 10^{-11},$ (32)

Combining it with the previous BNL result, the new result on the muon $(g-2)$
shows a deviation from the SM prediction at 4.2 $\sigma$ level:

$\displaystyle\Delta a^{\rm new}_{\mu}=(25.1\pm 5.9)\times 10^{-10},$ (33)

which may be a herald for new physics beyond the SM.

In our model, the new dominant contribution arises from the loop diagrams
involving the vector-like leptons and leptophilic DM $\chi_{l}$. The dominant
terms involve a product of $g_{\mu}$ and $h_{\mu}$ 444The contributions
proportional to $|g_{\mu}|^{2}$ or $|h_{\mu}|^{2}$ cannot be dominant because
of the chiral suppression. Thus, we neglect those terms for simplicity., and
its form is given by

$\displaystyle\Delta a^{\rm new}_{\mu}$
$\displaystyle=\frac{m_{\mu}}{(4\pi)^{2}}s_{c}c_{c}g_{\mu}h_{\mu}\left[\frac{F(q_{1},r_{1})}{M_{1}}-\frac{F(q_{2},r_{2})}{M_{2}}\right],$
(34) $\displaystyle F(q,r)$
$\displaystyle\equiv\int_{0}^{1}dx\int_{0}^{1-x}dy\frac{1-2y}{1-x+q^{2}(x^{2}-x)+r^{2}x}\approx\frac{-1+r^{4}-2r^{2}\ln(r^{2})}{2(-1+r^{2})^{3}}+{\cal
O}(q^{2}),$ (35)

where $q_{i}\equiv\frac{m_{\mu}}{M_{i}}(<<1)$ and
$r_{i}\equiv\frac{m_{\chi_{\ell}}}{M_{i}}$. And $M_{i}(i=1,2)$ and
$s_{c}(c_{c})$ are the mass eigenstates for the exotic singly-charged fermions
and their mixings given in Eq. (18), respectively.

###  III.4 Oblique parameters

The vector-like fermions and scalars contribute to oblique parameters through
vacuum polarization diagram for electroweak gauge bosons. Here we consider
$S$\- and $T$-parameters which would constrain our parameter space. They are
given by Peskin:1990zt

$\displaystyle\alpha_{em}S=4e^{2}\left(\sum_{\rm
VLF}\left[\frac{d}{dq^{2}}\Pi_{33}^{\rm VLF}-\frac{d}{dq^{2}}\Pi_{3Q}^{\rm
VLF}\right]_{q^{2}=0}+\left[\frac{d}{dq^{2}}\Pi_{33}^{h,H}-\frac{d}{dq^{2}}\Pi_{3Q}^{h,H}\right]_{q^{2}=0}\right),$
(36)
$\displaystyle\alpha_{em}T=\frac{e^{2}}{s_{W}^{2}c_{W}^{2}m_{Z}^{2}}\left(\sum_{\rm
VLF}\left[\Pi^{\rm VLF}_{\pm}(q^{2})-\Pi^{\rm
VLF}_{33}(q^{2})\right]_{q^{2}=0}+\left[\Pi^{h,H}_{\pm}(q^{2})-\Pi^{h,H}_{33}(q^{2})\right]_{q^{2}=0}\right),$
(37)

where superscript VLF distinguishes contributions from different combination
of our vector-like fermions $\\{E_{1},E_{2},N,U,D\\}$, and $\Pi_{33,3Q,\pm}$
are obtained from vacuum polarization diagrams which are summarized in the
Appendix; here $\Pi^{h,H}(33,3Q,\pm)$ indicates beyond the standard model
contribution. We impose the constraint on new physics contributions to the
$S$\- and $T$-parameters given by pdg

$\displaystyle\Delta S\leq 0.00\pm 0.07$ (38) $\displaystyle\Delta T\leq
0.05\pm 0.06$ (39)

where we fix $U$-parameter to be zero.

###  III.5  Dark matter

In our scenario, $U(1)_{\mu-\tau}$ is a good global symmetry of the model
Lagrangian at renormalizable level, and the lightest $U(1)_{\mu-\tau}$
particle with spin-0 would make a good DM candidate in addition to the
$\nu_{\mu}$ and $\nu_{\tau}$. Then the lighter field of two complex scalars
$\chi_{q}$ (hadro-philic or leptophobic) and $\chi_{\ell}$ (leptophilic) can
be DM candidate.

In fact DM physics in our model is very rich. For generic choices of $q_{x}$
and $q_{\varphi}$, $\chi_{q}$ and $\chi_{l}$ would be separately stable, and
we are with two-component DM models. However this would be no longer true for
special choices of $q_{\varphi}$ and $q_{x}$. For example, the following
operators would be gauge invariant while satisfying Eqs. (3) and (4):

$\displaystyle\chi_{l}\chi_{q}\varphi,~{}{\rm
if}~{}q_{\varphi}=1+\frac{3q_{x}}{2},$ (40)
$\displaystyle\chi_{l}\chi_{q}^{\dagger}\varphi^{2},~{}{\rm
if}~{}q_{\varphi}=\frac{1}{2}\left(1+\frac{q_{x}}{2}\right).$ (41)

Then $\chi_{q}$ and $\chi_{l}$ will mix with each other, and we would end up
with a single component DM model. In fact the mixing can occur without
$\varphi$: $\chi_{l}\chi_{q}$ is allowed if $q_{x}=2/3$. Another interesting
possibility is an operator $\chi_{l}^{\dagger 2}\chi_{q}\varphi^{\dagger}$
that is allowed for $q_{\varphi}=2$. Then after $U(1)_{\mu-\tau}$ symmetry
breaking, the effective $\chi_{l}^{\dagger 2}\chi_{q}$ operator is induced, so
that leptophobic DM $\chi_{q}$ will decay into a pair of leptophilic DM
$\chi_{l}$.

In this paper, we shall ignore all these interesting possibilities, and focus
on the leptophilic $\chi_{\ell}$ as a DM candidate that is simply denoted by
$\chi$ 555See Ref. Ko:2017yrd  for the case where $\chi_{q}$ is DM.. Then, it
annihilates into the SM leptons via $\chi\chi^{*}\to\ell^{+}\ell^{-}$ with
$\psi_{1,2}^{+}$ propagator, and $\chi\chi^{*}\to
Z^{\prime}\to\\{\mu^{+}\mu^{-},\tau^{+}\tau^{-},\nu_{\mu,\tau}\bar{\nu}_{\mu,\tau}\\}$,
666Here, the Yukawa contribution gives an $s$-wave dominant annihilation cross
section, while the $Z^{\prime}$ mediated cross section is the $p$-wave
dominant. indicating our DM is leptophilic. In fact there are annihilation
mode of $\chi\chi^{*}\to Z^{\prime}Z^{\prime}$ if kinematically allowed. In
our analysis we do not consider this mode assuming $m_{\chi}<m_{Z^{\prime}}$
to simplify our analysis. The relic density of DM is given by

$\displaystyle\Omega h^{2}\approx\frac{1.07\times
10^{9}}{\sqrt{g_{*}(x_{f})}M_{Pl}J(x_{f})[{\rm GeV}]},$ (42)

where $g^{*}(x_{f}\approx 25)\approx 100$, $M_{Pl}\approx 1.22\times 10^{19}$,
and $J(x_{f})(\equiv\int_{x_{f}}^{\infty}dx\frac{\langle\sigma v_{\rm
rel}\rangle}{x^{2}})$ is found as

$\displaystyle
J(x_{f})=\int_{x_{f}}^{\infty}dx\left[\frac{\int_{4m_{\chi}^{2}}^{\infty}ds\sqrt{s-4m_{\chi}^{2}}(\sigma
v_{\rm
rel})K_{1}\left(\frac{\sqrt{s}}{m_{\chi}}x\right)}{16m_{\chi}^{5}x[K_{2}(x)]^{2}}\right],$
$\displaystyle(\sigma v_{\rm rel})=(\sigma v_{\rm rel})_{\mu\bar{\mu}}+(\sigma
v_{\rm rel})_{\tau\bar{\tau}}+(\sigma v_{\rm
rel})_{\nu_{\mu}\bar{\nu}_{\mu}+\nu_{\tau}\bar{\nu}_{\tau}},$ (43)
$\displaystyle(\sigma v_{\rm rel})_{\mu\bar{\mu}}\simeq\frac{g^{\prime
4}q_{x}^{2}s(s-4m_{\chi}^{2})}{24\pi(s-m_{Z^{\prime}}^{2})^{2}}+\frac{(s_{c}c_{c}h_{\mu}g_{\mu})^{2}}{16\pi}\int_{0}^{\pi}d\theta\sin\theta\left(\frac{M_{1}}{t-M_{1}^{2}}-\frac{M_{2}}{t-M_{2}^{2}}\right)^{2}s,$
(44) $\displaystyle(\sigma v_{\rm rel})_{\tau\bar{\tau}}=\frac{g^{\prime
4}q_{x}^{2}(s+2m_{\tau}^{2})(s-4m_{\chi}^{2})}{24\pi(s-m_{Z^{\prime}}^{2})^{2}}\sqrt{1-\frac{4m_{\tau}^{2}}{s}},$
(45) $\displaystyle(\sigma v_{\rm
rel})_{\nu_{\mu}\bar{\nu}_{\mu}+\nu_{\tau}\bar{\nu}_{\tau}}=\frac{g^{\prime
4}q_{x}^{2}s(s-4m_{\chi}^{2})}{24\pi(s-m_{Z^{\prime}}^{2})^{2}}.$ (46)

Here $s$ and $t$ are Mandelstam variables, and $K_{1,2}$ are the modified
Bessel functions of the second kind of order 1 and 2, respectively. We write
second term of Eq. (44) in terms of integration by angle $\theta$ since
analytic form after integration is complicated as $t$ depends on $\theta$.
Note that integration by $\theta$ is carried out since corresponding
amplitudes only depend on $s$. In addition we ignore interference term between
$Z^{\prime}$ and $\Psi^{\pm}_{1,2}$ mediated diagrams since it is suppressed
by $m_{\mu}/m_{\chi}$. In our numerical analysis below, we use the current
experimental range for the relic density: $0.11\leq\Omega h^{2}\leq 0.13$
Ade:2013zuv .

Direct detection: We have a spin independent scattering cross section that is
constrained by several experiments such as XENON1T XENON:2015gkh . The
dominant contribution is Higgs portal via interactions in the scalar
potential. However we simply assume the Higgs portal coupling
$c_{\chi_{\ell}\chi_{\ell}h}$ for $\chi_{l}-\chi_{l}-h_{\rm SM}$ to be tiny
enough; $c_{\chi_{\ell}\chi_{\ell}h}\lesssim 10^{-3}$ Arcadi:2019lka , in
order to evade the direct detection via Higgs exchanges.

##  IV Numerical analysis

In this section, we perform the numerical analysis. Relevant free parameters
are scanned to search for parameter points satisfying phenomenological
requirements discussed in the previous section. The ranges of the other input
parameters are set to be as follows:

$\displaystyle g^{\prime}\in[0.001,0.1],\
|q_{x}|\in[0.1,\sqrt{4\pi}/g^{\prime}],\
\\{f_{1,2,3},g_{e,\mu},h_{e,\mu}\\}\in[10^{-2},1],\
m_{Z^{\prime}}\in[m_{\chi_{l}},150]\ [{\rm GeV}],\ $ $\displaystyle
m_{\chi_{\ell}}\in[1,100]\ [{\rm GeV}],\
m_{\chi_{q}}\in[1.2m_{\chi_{\ell}},10^{4}]\ [{\rm GeV}],\
\\{M_{E,L^{\prime}}\\}\in[10^{2},10^{4}]\ [{\rm GeV}],\ $ (47) $\displaystyle
m_{E}\in[10,100]\ [{\rm GeV}],\ M_{Q^{\prime}}\in[10^{3},10^{4}]\ [{\rm
GeV}],$

where we take perturbation limit $|q_{x}g^{\prime}|\leq\sqrt{4\pi}$, and
choose $M_{1}<M_{2}$, and take $m_{\chi_{\ell}}<1.2\times(M_{1},m_{\chi_{q}})$
for simplicity so that we can evade contributions from co-annihilation
processes. Here we require 1000 GeV$\leq M_{Q^{\prime}}$ to avoid constraints
from direct search of vector-like quarks at the LHC. For vector-like lepton,
constraint $m_{VLL}\gtrsim 700$ GeV is given for vector-like lepton decaying
into the SM lepton with $W/Z$ boson Sirunyan:2019ofn . But here we just
consider LEP limit of $m_{VLL}\gtrsim 100$ GeV Achard:2001qw  since our
vector-like lepton decays into lepton and $\chi_{l}$. In addition we impose
the constraint of $550$ GeV $\lesssim m_{Z^{\prime}}/g^{\prime}$ which arises
from the neutrino trident production Altmannshofer:2014pba . Here we take
scalar mixing $\sin\alpha$ to be sufficiently small so that scalar boson
contribution to oblique parameters can avoid constraints since it is
proportional to $\sin^{2}\alpha$. We find that new scalar contribution to $S$
and $T$ parameter is sufficiently small for $\sin\alpha<0.1$. Moreover, taking
$m_{E}<M_{E,L^{\prime}}$, we find that constraint from $S$\- and
$T$-parameters can be easily avoided where $S$ and $T$ are found to be less
than $\sim 0.05$ and $0.01$, respectively.

Firstly, in Fig. 2, we show $C_{9,10}$ values for parameter sets satisfying
the constraints of neutral meson-antimeson mixings $\Delta m_{K},\Delta
m_{B_{d}},\Delta m_{B_{s}},\Delta m_{D}$, correct relic density
$0.11\leq\Omega h^{2}\leq 0.13$ and $0<\Delta C_{10}(Z)<0.79$; here we show
region with $\Delta C_{9}^{\mu\mu}(Z^{\prime})<0$ and $|\Delta C_{9,10}(Z)|$
where region with $\Delta C_{9}^{\mu\mu}(Z^{\prime})>0$ and $\Delta
C_{9,10}>(<)0$ similarly appear. We find that $|\Delta C_{9}(Z)|$ tends to be
small due to suppression factor of $(2s^{2}_{w}-1/2)$. Then the value of
$\Delta C_{9}^{\mu\mu}$ can be $0.83\leq-\Delta C_{9}^{\mu\mu}\leq 1.52$
indicated by region between dashed line which is obtained from global analysis
in Ref. Alguero:2021anc  (see also Ref. Carvunis:2021jga ). Note that non-zero
$\Delta C_{10}$ helps to fit $B_{s}\to\mu\mu$ data. In fact non-zero $\Delta
C^{ee}_{9,10}(Z)$ affect the region from global fit but we still use the
region as our reference. We also find that $C_{7}$ is tiny in our allowed
region as $|C_{7}|<0.001$ and the corresponding plot is omitted.

In Fig. 3, we show our values of $\\{M_{1},a_{\mu}\\}$ for the parameter sets
the same as Fig. 2. We can obtain $\Delta a^{\rm new}_{\mu}=(25.1\pm
5.9)\times 10^{-10}$ at 1(2)$\sigma$ for $M_{1}\lesssim 2000(2700)$ GeV. Here,
the dashed line denotes the allowed region within 2$\sigma$.

Figure 2: Left: scattering plot on $\\{\Delta
C_{9}^{\mu\mu}(Z^{\prime}),\Delta C_{10}(Z)\\}$ plain where each point
satisfies phenomenological constraints except for $\Delta C_{9,10}$. Right the
same plot on $\\{\Delta C_{9}^{\mu\mu}(Z^{\prime}),|\Delta C_{9}(Z)|\\}$
plain. Figure 3: $M_{1}$ versus $\Delta a_{\mu}$ for the parameter sets same
as the previous plots.

In Fig. 4, we show the allowed parameter region for DM mass $m_{\chi_{\ell}}$
and $m_{Z^{\prime}}$. The relic density is explained by $\chi\chi^{*}\to
Z^{\prime}\to\\{\mu^{+}\mu^{-},\tau^{+}\tau^{-},\nu_{\mu,\tau}\bar{\nu}_{\mu,\tau}\\}$
process for the points around $m_{Z^{\prime}}\sim 2m_{\chi}$ due to
enhancement of the annihilation cross section. On the other hand
$\chi\chi^{*}\to\mu^{+}\mu^{-}$ via Yukawa interaction explains the relic
density for the points which do not give resonant enhancement of the
$Z^{\prime}$ exchange process.

Figure 4: $m_{\chi_{\ell}}$ versus $m_{Z^{\prime}}$ in unit of GeV for the
parameter sets same as the previous plots.

Figure 5:  Left: Parameter points on $(m_{Z^{\prime}},g^{\prime})$ plane
satisfying $0.83\leq-\Delta C_{9}^{\mu\mu}\leq 1.52$, $\Delta a^{\rm
new}_{\mu}=(25.1\pm 5.9)\times 10^{-10}$ and other phenomenological
constraints. Shaded region is excluded by the LHC constraint from
$pp\to\mu^{+}\mu^{-}Z^{\prime}(\to\mu^{+}\mu^{-})$ process Sirunyan:2018nnz .
Right: the same points on $(g^{\prime},|q_{x}|)$ plane.

Finally, in left plot of Fig. 5, we show parameter points on
$\\{m_{Z^{\prime}},g_{X}\\}$ plane satisfying $0.83\leq-\Delta
C_{9}^{\mu\mu}\leq 1.52$, $\Delta a^{\rm new}_{\mu}=(25.1\pm 5.9)\times
10^{-10}$ and other phenomenological constraints. In this region the most
important one is the LHC constraint from
$pp\to\mu^{+}\mu^{-}Z^{\prime}(\to\mu^{+}\mu^{-})$ process Sirunyan:2018nnz .
We find that $m_{Z^{\prime}}\lesssim m_{Z}$ region is preferred to explain
experimental anomalies where some parameter region is excluded by the LHC
constraint for $m_{Z^{\prime}}<60$ GeV. The gauge coupling should be small for
small $m_{Z^{\prime}}$ region to avoid the LHC constraint, and large $|q_{X}|$
is required to explain $B$-anomalies in the region. In the right plot in Fig.
5 we show correlation between $g^{\prime}$ and $|q_{x}|$ for the same
parameter sets as the left plot; both negative and positive values of $q_{x}$
appear as same amount and we only show absolute values. We then find
$g^{\prime}|q_{x}|\sim\mathcal{O}(1)$ is required to obtain sizable $\Delta
C_{9}^{\mu\mu}$. Thus the region with $m_{Z^{\prime}}>60$ GeV would be more
natural to explain $b\to s\mu^{+}\mu^{-}$ anomalies since required $|q_{x}|$
is not too large. The remaining allowed region would be further tested in
future at the LHC from the search for
$pp\to\mu^{+}\mu^{-}Z^{\prime}(\to\mu^{+}\mu^{-},\nu_{\mu,\tau}\bar{\nu}_{\mu,\tau})$
and $pp\to\nu_{\mu,\tau}\bar{\nu}_{\mu,\tau}Z^{\prime}(\to\mu\bar{\mu})$:
namely 4 muons or dimuon + missing transverse momentum.

##  V Summary and Conclusions

In this paper, we have proposed a $U(1)_{\mu-\tau}$ gauged model to explain
recent anomalies; the muon anomalous magnetic moment and $b\to
s\mu^{+}\mu^{-}$ at the same time, satisfying several constraints. Also, we
have considered a complex scalar boson DM candidate ($\chi_{l}$) that mainly
couples to leptons. In our model $b\to s\mu^{+}\mu^{-}$ anomalies can be
explained by loop induced interactions among quarks and $Z^{\prime}$. We also
estimated one-loop diagrams associated with vector-like quark and $Z(\gamma)$
vertices that contribute to Wilson coefficients $C_{9(10)}$ and $C_{7}$ to
check consistency of our model. We then find that $C_{9(10)}$ from $Z$
mediation and $C_{7}$ induced by vector-like quark loop can be small, and we
can explain $b\to s\mu^{+}\mu^{-}$ anomalies consistently. In our numerical
analysis, we have found allowed regions to explain anomalies while satisfying
phenomenological constraints. This model would be test soon in future
measurement such as $pp\to\mu^{+}\mu^{-}Z^{\prime}(\to\mu^{+}\mu^{-})$ signal
search at the LHC. We also have the processes
$pp\to\mu^{+}\mu^{-}Z^{\prime}(\to\nu_{\mu,\tau}\bar{\nu}_{\mu,\tau})$ and
$pp\to\nu_{\mu,\tau}\bar{\nu}_{\mu,\tau}Z^{\prime}(\to\mu^{+}\mu^{-})$ giving
signal of $\mu^{+}\mu^{-}$ with missing transverse momentum that is also
searched for at the LHC. In addition our model can be tested by searching for
vector-like fermions at the LHC where vector-like quarks and leptons can be
pair produced providing ”dijet + missing transverse momentum” and ”dilepton +
missing transverse momentum” respectively.

## Acknowledgments

The work of P.K. is supported in part by KIAS Individual Grants under Grant
No.PG021403, and by National Research Foundation of Korea (NRF) Research Grant
NRF- 2019R1A2C3005009. The work of H.O. is supported by an appointment to the
JRG Program at the APCTP through the Science and Technology Promotion Fund and
Lottery Fund of the Korean Government, and also by the Korean Local
Governments - Gyeongsangbuk-do Province and Pohang City H. O. is sincerely
grateful for the KIAS members.

##  Appendix A Vector-like fermion contributions to electroweak vacuum
polarization diagrams

We calculate electroweak vacuum polarization diagrams with vector-like
fermions and scalar bosons, and summarize analytic form of their contributions
to estimate oblique parameters. We can write the vacuum polarizations for $Z$
and $W$ contributing to oblique parameters such that

$\displaystyle\Pi_{Z}^{\mu\nu}=g^{\mu\nu}\frac{e^{2}}{c_{W}^{2}s_{W}^{2}}(\Pi_{33}(q^{2})-2s_{W}^{2}\Pi_{3Q}(q^{2})-s_{W}^{4}\Pi_{QQ}(q^{2})),$
(48)
$\displaystyle\Pi^{\mu\nu}_{W}=g^{\mu\nu}\frac{e^{2}}{s_{W}^{2}}\Pi_{\pm}(q^{2}),$
(49)

where $q$ is four momentum carried by gauge bosons.

Non-zero contributions to $\Pi_{33}(q^{2})$ are summarized as follows;

$\displaystyle\Pi^{E_{1}E_{1}}_{33}(q^{2})=-\frac{s_{c}^{4}}{16\pi^{2}}F(q^{2},M_{1}^{2},M_{1}^{2}),$
(50)
$\displaystyle\Pi^{E_{2}E_{2}}_{33}(q^{2})=-\frac{c_{c}^{4}}{16\pi^{2}}F(q^{2},M_{2}^{2},M_{2}^{2}),$
(51)
$\displaystyle\Pi^{E_{1}E_{2}}_{33}(q^{2})=-\frac{s_{c}^{2}c_{c}^{2}}{8\pi^{2}}F(q^{2},M_{1}^{2},M_{2}^{2}),$
(52)
$\displaystyle\Pi^{N^{\prime}N^{\prime}}_{33}(q^{2})=-\frac{1}{16\pi^{2}}F(q^{2},M_{N^{\prime}}^{2},M_{N^{\prime}}^{2}),$
(53)
$\displaystyle\Pi^{U^{\prime}U^{\prime}}_{33}(q^{2})=-\frac{1}{16\pi^{2}}F(q^{2},M_{U^{\prime}}^{2},M_{U^{\prime}}^{2}),$
(54)
$\displaystyle\Pi^{D^{\prime}D^{\prime}}_{33}(q^{2})=-\frac{1}{16\pi^{2}}F(q^{2},M_{D^{\prime}}^{2},M_{D^{\prime}}^{2}),$
(55)

where the superscripts in the left sides indicate particles inside vacuum
polarization diagrams, and $F(q^{2},m^{2},m^{\prime 2})$ is the loop function
given by

$\displaystyle F(q^{2},m^{2},m^{\prime 2})=$
$\displaystyle\int_{0}^{1}dxdy\delta(1-x-y)\left(\frac{1}{\epsilon_{\rm\overline{MS}}}-\ln\left(\frac{\Delta}{\mu^{2}}\right)\right)$
$\displaystyle\times(2x(1-x)q^{2}-xm^{2}-ym^{\prime 2}+mm^{\prime}),$ (56)
$\displaystyle\Delta=$ $\displaystyle-x(1-x)q^{2}+xm^{2}+ym^{\prime
2},\quad\frac{1}{\epsilon_{\rm\overline{MS}}}\equiv\frac{2}{\epsilon}-\gamma-\ln(4\pi),$
(57)

where $\mu$ is auxiliary parameter having mass dimension. Dependence on $\mu$
is canceled when we calculate $S,T$-parameter from Eqs. (28) and (29) in Sec.
III D.

Similarly we obtain non-zero contributions to $\Pi_{3Q,QQ,\pm}$:

$\displaystyle\Pi^{E_{1}E_{1}}_{3Q}(q^{2})=-\frac{s_{c}^{2}}{8\pi^{2}}F(q^{2},M_{1}^{2},M_{1}^{2}),$
(58)
$\displaystyle\Pi^{E_{2}E_{2}}_{3Q}(q^{2})=-\frac{c_{c}^{2}}{8\pi^{2}}F(q^{2},M_{2}^{2},M_{2}^{2}),$
(59)
$\displaystyle\Pi^{U^{\prime}U^{\prime}}_{3Q}(q^{2})=-\frac{1}{12\pi^{2}}F(q^{2},M_{U^{\prime}}^{2},M_{U^{\prime}}^{2}),$
(60)
$\displaystyle\Pi^{D^{\prime}D^{\prime}}_{3Q}(q^{2})=-\frac{1}{24\pi^{2}}F(q^{2},M_{D^{\prime}}^{2},M_{D^{\prime}}^{2}),$
(61)
$\displaystyle\Pi^{E_{1}E_{1}}_{QQ}(q^{2})=-\frac{1}{4\pi^{2}}F(q^{2},M_{1}^{2},M_{1}^{2}),$
(62)
$\displaystyle\Pi^{E_{2}E_{2}}_{QQ}(q^{2})=-\frac{1}{4\pi^{2}}F(q^{2},M_{2}^{2},M_{2}^{2}),$
(63)
$\displaystyle\Pi^{U^{\prime}U^{\prime}}_{QQ}(q^{2})=-\frac{1}{9\pi^{2}}F(q^{2},M_{U^{\prime}}^{2},M_{U^{\prime}}^{2}),$
(64)
$\displaystyle\Pi^{D^{\prime}D^{\prime}}_{QQ}(q^{2})=-\frac{1}{36\pi^{2}}F(q^{2},M_{D^{\prime}}^{2},M_{D^{\prime}}^{2}),$
(65)
$\displaystyle\Pi^{E_{1}N^{\prime}}_{\pm}(q^{2})=-\frac{s_{c}^{2}}{8\pi^{2}}F(q^{2},M_{1}^{2},M_{N^{\prime}}^{2}),$
(66)
$\displaystyle\Pi^{E_{2}N^{\prime}}_{\pm}(q^{2})=-\frac{c_{c}^{2}}{8\pi^{2}}F(q^{2},M_{2}^{2},M_{N^{\prime}}^{2}),$
(67)
$\displaystyle\Pi^{U^{\prime}D^{\prime}}_{\pm}(q^{2})=-\frac{1}{8\pi^{2}}F(q^{2},M_{U^{\prime}}^{2},M_{D^{\prime}}^{2}).$
(68)

Finally non-zero contributions from new scalar particles beyond the SM are
given by

$\displaystyle\Pi_{33}^{h,H}=-\frac{\sin^{2}\alpha}{32\pi^{2}}\int_{0}^{1}dx\left[2m_{Z}^{2}\ln\left(\frac{\Delta_{h}}{\Delta_{H}}\right)+\Delta_{h}\left(\frac{1}{\epsilon_{\rm\overline{MS}}}-\ln\left(\frac{\Delta_{h}}{\mu^{2}}\right)\right)-\Delta_{H}\left(\frac{1}{\epsilon_{\rm\overline{MS}}}-\left(\frac{\Delta_{H}}{\mu^{2}}\right)\right)\right],$
$\displaystyle\Pi_{\pm}^{h,H}=-\frac{\sin^{2}\alpha}{32\pi^{2}}\int_{0}^{1}dx\left[2m_{W}^{2}\ln\left(\frac{\Delta_{h}}{\Delta^{\prime}_{H}}\right)+\Delta^{\prime}_{h}\left(\frac{1}{\epsilon_{\rm\overline{MS}}}-\left(\frac{\Delta^{\prime}_{h}}{\mu^{2}}\right)\right)-\Delta^{\prime}_{H}\left(\frac{1}{\epsilon_{\rm\overline{MS}}}-\left(\frac{\Delta^{\prime}_{H}}{\mu^{2}}\right)\right)\right],$
$\displaystyle\Delta_{h(H)}=x(1-x)q^{2}+xm^{2}_{h(H)}+(1-x)m_{Z}^{2},$
$\displaystyle\Delta^{\prime}_{h(H)}=x(1-x)q^{2}+xm^{2}_{h(H)}+(1-x)m_{W}^{2}.$
(69)

This completes the new contributions to the vacuum polarization tensors of
electroweak gauge bosons in the $U(1)_{\mu-\tau}$ model considered in this
paper. Then the $S$\- and $T$-parameters can be obtained from Eqs. (28) and
(29) in Sec. III D; divergent part proportional to $\epsilon_{\rm SM}$ will be
canceled when we calculate oblique parameters by those equations.

## References

  * (1) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, JHEP 09, 178 (2020) doi:10.1007/JHEP09(2020)178 [arXiv:2007.14792 [hep-ph]].
  * (2) P. A. R. Ade et al. [Planck], Astron. Astrophys. 594, A13 (2016) doi:10.1051/0004-6361/201525830 [arXiv:1502.01589 [astro-ph.CO]].
  * (3) G. Hiller and F. Kruger, Phys. Rev. D 69, 074020 (2004) [hep-ph/0310219].
  * (4) C. Bobeth, G. Hiller and G. Piranishvili, JHEP 0712, 040 (2007) [arXiv:0709.4174 [hep-ph]].
  * (5) R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 113, 151601 (2014) [arXiv:1406.6482 [hep-ex]].
  * (6) R. Aaij et al. [LHCb], Phys. Rev. Lett. 122 (2019) no.19, 191801 [arXiv:1903.09252 [hep-ex]].
  * (7) G. W. Bennett et al. [Muon g-2], Phys. Rev. D 73, 072003 (2006) [arXiv:hep-ex/0602035 [hep-ex]].
  * (8) B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126, 141801 (2021) [arXiv:2104.03281 [hep-ex]].
  * (9) K. Hagiwara, R. Liao, A. D. Martin, D. Nomura and T. Teubner, J. Phys. G 38, 085003 (2011) [arXiv:1105.3149 [hep-ph]].
  * (10) T. Aoyama, M. Hayakawa, T. Kinoshita and M. Nio, Phys. Rev. Lett. 109, 111808 (2012) doi:10.1103/PhysRevLett.109.111808 [arXiv:1205.5370 [hep-ph]].
  * (11) T. Aoyama, T. Kinoshita and M. Nio, Atoms 7, no.1, 28 (2019) doi:10.3390/atoms7010028
  * (12) A. Czarnecki, W. J. Marciano and A. Vainshtein, Phys. Rev. D 67, 073006 (2003) [erratum: Phys. Rev. D 73, 119901 (2006)] doi:10.1103/PhysRevD.67.073006 [arXiv:hep-ph/0212229 [hep-ph]].
  * (13) C. Gnendiger, D. Stöckinger and H. Stöckinger-Kim, Phys. Rev. D 88, 053005 (2013) doi:10.1103/PhysRevD.88.053005 [arXiv:1306.5546 [hep-ph]].
  * (14) A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 97, no.11, 114025 (2018) doi:10.1103/PhysRevD.97.114025 [arXiv:1802.02995 [hep-ph]].
  * (15) G. Colangelo, M. Hoferichter and P. Stoffer, JHEP 02, 006 (2019) doi:10.1007/JHEP02(2019)006 [arXiv:1810.00007 [hep-ph]].
  * (16) M. Hoferichter, B. L. Hoid and B. Kubis, JHEP 08, 137 (2019) doi:10.1007/JHEP08(2019)137 [arXiv:1907.01556 [hep-ph]].
  * (17) A. Keshavarzi, D. Nomura and T. Teubner, Phys. Rev. D 101, no.1, 014029 (2020) doi:10.1103/PhysRevD.101.014029 [arXiv:1911.00367 [hep-ph]].
  * (18) A. Kurz, T. Liu, P. Marquard and M. Steinhauser, Phys. Lett. B 734, 144-147 (2014) doi:10.1016/j.physletb.2014.05.043 [arXiv:1403.6400 [hep-ph]].
  * (19) K. Melnikov and A. Vainshtein, Phys. Rev. D 70, 113006 (2004) doi:10.1103/PhysRevD.70.113006 [arXiv:hep-ph/0312226 [hep-ph]].
  * (20) P. Masjuan and P. Sanchez-Puertas, Phys. Rev. D 95, no.5, 054026 (2017) doi:10.1103/PhysRevD.95.054026 [arXiv:1701.05829 [hep-ph]].
  * (21) G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, JHEP 04, 161 (2017) doi:10.1007/JHEP04(2017)161 [arXiv:1702.07347 [hep-ph]].
  * (22) M. Hoferichter, B. L. Hoid, B. Kubis, S. Leupold and S. P. Schneider, JHEP 10, 141 (2018) doi:10.1007/JHEP10(2018)141 [arXiv:1808.04823 [hep-ph]].
  * (23) A. Gérardin, H. B. Meyer and A. Nyffeler, Phys. Rev. D 100, no.3, 034520 (2019) doi:10.1103/PhysRevD.100.034520 [arXiv:1903.09471 [hep-lat]].
  * (24) J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez, Phys. Lett. B 798, 134994 (2019) doi:10.1016/j.physletb.2019.134994 [arXiv:1908.03331 [hep-ph]].
  * (25) G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, JHEP 03, 101 (2020) doi:10.1007/JHEP03(2020)101 [arXiv:1910.13432 [hep-ph]].
  * (26) T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung and C. Lehner, Phys. Rev. Lett. 124, no.13, 132002 (2020) doi:10.1103/PhysRevLett.124.132002 [arXiv:1911.08123 [hep-lat]].
  * (27) G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera and P. Stoffer, Phys. Lett. B 735, 90-91 (2014) doi:10.1016/j.physletb.2014.06.012 [arXiv:1403.7512 [hep-ph]].
  * (28) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 77, no.12, 827 (2017) doi:10.1140/epjc/s10052-017-5161-6 [arXiv:1706.09436 [hep-ph]].
  * (29) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 80, no.3, 241 (2020) [erratum: Eur. Phys. J. C 80, no.5, 410 (2020)] doi:10.1140/epjc/s10052-020-7792-2 [arXiv:1908.00921 [hep-ph]].
  * (30) R. Aaij et al. [LHCb], [arXiv:2103.11769 [hep-ex]].
  * (31) S. Descotes-Genon, J. Matias, M. Ramon and J. Virto, JHEP 1301, 048 (2013) [arXiv:1207.2753 [hep-ph]].
  * (32) R. Aaij et al. [LHCb Collaboration], JHEP 1602, 104 (2016) [arXiv:1512.04442 [hep-ex]].
  * (33) R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 111, 191801 (2013) [arXiv:1308.1707 [hep-ex]].
  * (34) A. Abdesselam et al. [Belle Collaboration], arXiv:1604.04042 [hep-ex].
  * (35) S. Wehle et al. [Belle Collaboration], arXiv:1612.05014 [hep-ex].
  * (36) R. Aaij et al. [LHCb], JHEP 08 (2017), 055 [arXiv:1705.05802 [hep-ex]].
  * (37) X. G. He, G. C. Joshi, H. Lew and R. R. Volkas, Phys. Rev. D 43, 22-24 (1991) doi:10.1103/PhysRevD.43.R22
  * (38) X. G. He, G. C. Joshi, H. Lew and R. R. Volkas, Phys. Rev. D 44, 2118-2132 (1991) doi:10.1103/PhysRevD.44.2118
  * (39) W. Altmannshofer, S. Gori, M. Pospelov and I. Yavin, Phys. Rev. Lett. 113, 091801 (2014) doi:10.1103/PhysRevLett.113.091801 [arXiv:1406.2332 [hep-ph]].
  * (40) S. Baek, N. G. Deshpande, X. G. He and P. Ko, Phys. Rev. D 64, 055006 (2001) doi:10.1103/PhysRevD.64.055006 [arXiv:hep-ph/0104141 [hep-ph]].
  * (41) S. Baek and P. Ko, JCAP 10, 011 (2009) doi:10.1088/1475-7516/2009/10/011 [arXiv:0811.1646 [hep-ph]].
  * (42) S. Baek, H. Okada and K. Yagyu, JHEP 04 (2015), 049 [arXiv:1501.01530 [hep-ph]].
  * (43) T. Nomura and H. Okada, Phys. Lett. B 783 (2018), 381-386 [arXiv:1805.03942 [hep-ph]].
  * (44) T. Nomura and H. Okada, Phys. Rev. D 97 (2018) no.9, 095023 [arXiv:1803.04795 [hep-ph]].
  * (45) K. Asai, K. Hamaguchi, N. Nagata, S. Y. Tseng and K. Tsumura, Phys. Rev. D 99 (2019) no.5, 055029 [arXiv:1811.07571 [hep-ph]].
  * (46) K. Asai, K. Hamaguchi and N. Nagata, Eur. Phys. J. C 77 (2017) no.11, 763 [arXiv:1705.00419 [hep-ph]].
  * (47) S. Lee, T. Nomura and H. Okada, Nucl. Phys. B 931 (2018), 179-191 [arXiv:1702.03733 [hep-ph]].
  * (48) M. Davier, A. Hoecker, B. Malaescu and Z. Zhang, Eur. Phys. J. C 71, 1515 (2011) [erratum: Eur. Phys. J. C 72, 1874 (2012)] doi:10.1140/epjc/s10052-012-1874-8 [arXiv:1010.4180 [hep-ph]].
  * (49) D. Borah, A. Dasgupta and D. Mahanta, [arXiv:2106.14410 [hep-ph]].
  * (50) X. Qi, A. Yang, W. Liu and H. Sun, [arXiv:2106.14134 [hep-ph]].
  * (51) S. Singirala, S. Sahoo and R. Mohanta, [arXiv:2106.03735 [hep-ph]].
  * (52) A. J. Buras, A. Crivellin, F. Kirk, C. A. Manzari and M. Montull, JHEP 06 (2021), 068 [arXiv:2104.07680 [hep-ph]].
  * (53) S. Zhou, [arXiv:2104.06858 [hep-ph]].
  * (54) D. Borah, M. Dutta, S. Mahapatra and N. Sahu, [arXiv:2104.05656 [hep-ph]].
  * (55) J. Chen, Q. Wen, F. Xu and M. Zhang, [arXiv:2104.03699 [hep-ph]].
  * (56) L. Zu, X. Pan, L. Feng, Q. Yuan and Y. Z. Fan, [arXiv:2104.03340 [hep-ph]].
  * (57) G. y. Huang, F. S. Queiroz and W. Rodejohann, Phys. Rev. D 103 (2021) no.9, 095005 [arXiv:2101.04956 [hep-ph]].
  * (58) S. Patra, S. Rao, N. Sahoo and N. Sahu, Nucl. Phys. B 917, 317-336 (2017) doi:10.1016/j.nuclphysb.2017.02.010 [arXiv:1607.04046 [hep-ph]].
  * (59) W. Altmannshofer, M. Carena and A. Crivellin, Phys. Rev. D 94, no.9, 095026 (2016) doi:10.1103/PhysRevD.94.095026 [arXiv:1604.08221 [hep-ph]].
  * (60) W. Altmannshofer, S. Gori, M. Pospelov and I. Yavin, Phys. Rev. D 89, 095033 (2014) doi:10.1103/PhysRevD.89.095033 [arXiv:1403.1269 [hep-ph]].
  * (61) A. Crivellin, G. D’Ambrosio and J. Heeck, Phys. Rev. Lett. 114, 151801 (2015) doi:10.1103/PhysRevLett.114.151801 [arXiv:1501.00993 [hep-ph]].
  * (62) A. Crivellin, G. D’Ambrosio and J. Heeck, Phys. Rev. D 91, no.7, 075006 (2015) doi:10.1103/PhysRevD.91.075006 [arXiv:1503.03477 [hep-ph]].
  * (63) P. Ko, T. Nomura and H. Okada, Phys. Rev. D 95 (2017) no.11, 111701 [arXiv:1702.02699 [hep-ph]].
  * (64) C. H. Chen and T. Nomura, Phys. Lett. B 777, 420-427 (2018) doi:10.1016/j.physletb.2017.12.062 [arXiv:1707.03249 [hep-ph]].
  * (65) S. Baek, JHEP 05, 104 (2019) doi:10.1007/JHEP05(2019)104 [arXiv:1901.04761 [hep-ph]].
  * (66) Z. L. Han, R. Ding, S. J. Lin and B. Zhu, Eur. Phys. J. C 79, no.12, 1007 (2019) doi:10.1140/epjc/s10052-019-7526-5 [arXiv:1908.07192 [hep-ph]].
  * (67) N. Kumar, T. Nomura and H. Okada, [arXiv:2002.12218 [hep-ph]].
  * (68) W. Chao, H. Wang, L. Wang and Y. Zhang, Chin. Phys. C 45, no.8, 083105 (2021) doi:10.1088/1674-1137/ac06ba [arXiv:2102.07518 [hep-ph]].
  * (69) D. Borah, M. Dutta, S. Mahapatra and N. Sahu, [arXiv:2109.02699 [hep-ph]].
  * (70) R. Aaij et al. [LHCb], [arXiv:2108.09283 [hep-ex]].
  * (71) R. Aaij et al. [LHCb], [arXiv:2108.09284 [hep-ex]].
  * (72) G. Hiller and M. Schmaltz, Phys. Rev. D 90, 054014 (2014) doi:10.1103/PhysRevD.90.054014 [arXiv:1408.1627 [hep-ph]].
  * (73) C. Bobeth, M. Gorbahn, T. Hermann, M. Misiak, E. Stamou and M. Steinhauser, Phys. Rev. Lett. 112, 101801 (2014) doi:10.1103/PhysRevLett.112.101801 [arXiv:1311.0903 [hep-ph]].
  * (74) M. F. Navarro and S. F. King, [arXiv:2109.08729 [hep-ph]].
  * (75) P. A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * (76) M. E. Peskin and T. Takeuchi, Phys. Rev. Lett. 65, 964-967 (1990)
  * (77) M. Algueró, B. Capdevila, S. Descotes-Genon, J. Matias and M. Novoa-Brunet, [arXiv:2104.08921 [hep-ph]].
  * (78) A. Carvunis, F. Dettori, S. Gangal, D. Guadagnoli and C. Normand, [arXiv:2102.13390 [hep-ph]]. Copy to ClipboardDownload
  * (79) F. Gabbiani, E. Gabrielli, A. Masiero and L. Silvestrini, Nucl. Phys. B 477, 321 (1996) [hep-ph/9604387].
  * (80) P. A. R. Ade et al. [Planck Collaboration], Astron. Astrophys. 571, A16 (2014) [arXiv:1303.5076 [astro-ph.CO]].
  * (81) E. Aprile et al. [XENON], JCAP 04, 027 (2016) doi:10.1088/1475-7516/2016/04/027 [arXiv:1512.07501 [physics.ins-det]].
  * (82) G. Arcadi, A. Djouadi and M. Raidal, Phys. Rept. 842, 1-180 (2020) doi:10.1016/j.physrep.2019.11.003 [arXiv:1903.03616 [hep-ph]].
  * (83) P. Achard et al. [L3 Collaboration], Phys. Lett. B 517, 75 (2001) [hep-ex/0107015].
  * (84) A. M. Sirunyan et al. [CMS Collaboration], Phys. Rev. D 100, no. 5, 052003 (2019) [arXiv:1905.10853 [hep-ex]].
  * (85) A. M. Sirunyan et al. [CMS Collaboration], Phys. Lett. B 792, 345 (2019) [arXiv:1808.03684 [hep-ex]].


# Predicting Generalization in Deep Learning via Metric Learning – PGDL Shared
task

Sebastian Mežnar<EMAIL_ADDRESS>
Jožef Stefan Institute
Jamova 39, 1000 Ljubljana, Slovenia Blaž Škrlj<EMAIL_ADDRESS>
Jožef Stefan Institute
Jamova 39, 1000 Ljubljana, Slovenia
Jožef Stefan International Postgraduate School
Jamova 39,1000 Ljubljana, Slovenia

###### Abstract

The competition “Predicting Generalization in Deep Learning (PGDL)” aims to
provide a platform for rigorous study of generalization of deep learning
models and offer insight into the progress of understanding and explaining
these models. This report presents the solution that was submitted by the user
_smeznar_ which achieved the eight place in the competition. In the proposed
approach, we create simple metrics and find their best combination with
automatic testing on the provided dataset, exploring how combinations of
various properties of the input neural network architectures can be used for
the prediction of their generalization.

##  1 Introduction

Generalization is a difficult problem in deep learning since neural models are
a complex, not directly interpretable branch of models. Due to a large amount
of parameters in neural network models, training objective has commonly
multiple global minima that all minimize the error on the training set well,
but can fail to generalize well on new data. Because of this property,
choosing the right minima is crucial to achieve great performance across all
data Neyshabur et al. (2017).

Recently, an extensive study of complexity measures for generalization was
conducted, where authors trained over 10,000 convolutional networks and tested
more than 40 complexity measures Jiang et al. (2020b). During testing, they
found that PAC-Bayesian bounds are a promising direction for solving
generalization. Another relevant observation is that the evaluation of
generalization complexity measures can be misleading when the number of
studied models is small. When considering how instances differ from one
another, the notion of _distance_ between them becomes of crucial relevance.
This work also builds on some of the key ideas of _metric learning_ Li and
Tian (2018), i.e. the process of automatically determining the appropriate
distance score for a given problem.

The goal of the competition is to create a complexity measure that scores how
well a model generalizes. The competition scores complexity measures based on
the correlation between the ranking obtained from the complexity measure and
the ground-truth ranking obtained empirically by ranking models based on their
generalization gap. The generalization gap is defined as

$g=\bigg{|}\frac{1}{|X|}\sum_{x_{i}\in
X}\mathbbm{1}[f(x_{i})=y_{i}]-\frac{1}{|\overline{X}|}\sum_{\overline{x}_{i}\in\overline{X}}\mathbbm{1}[f(\overline{x_{i}})=\overline{y_{i}}]\bigg{|},$
(1)

where $x_{i}\in X$ is the training data, $\overline{x}_{i}\in\overline{X}$ is
the validation data, $f$ is the classification model and $y$ and
$\overline{y}$ are the ground-truth values of training and validation data
Jiang et al. (2020a). The metric used for calculating the correlation between
two rankings is the Kendall’s Rank-Correlation Coefficient Kendall (1938).
This coefficient is calculated on pairs representing the same rank using
formula:

$\tau=\frac{\text{concordant pairs}-\text{discordant pairs}}{{n\choose 2}}.$
(2)

Two pairs of points ($x_{i},y_{i}$) and ($x_{j},y_{j}$), where $i<j$ are said
to be concordant if $x_{i}<x_{j}$ and $y_{i}<y_{j}$ holds or $x_{i}>x_{j}$ and
$y_{i}>y_{j}$ holds, otherwise two pairs are discordant.

In our work we investigated whether, because of the complexity of deep
learning models, there does not exist a complexity measure that exploits only
some part of a model structure and gives a very accurate generalization score.
Instead we focus on finding a collection of different metrics, that, when
combined, give an accurate approximation of the generalization. For this we
create simple metrics that represent different aspects of the model (and the
data set) structure and combine them into the final approximation. Since this
combination might not be trivial, we explore the combination space
_automatically_.

##  2 Methodology

We created seven simple metrics that give a score between 0 and 1. These
metrics are the following:

  * •

Non-uniform weights: The percentage of weights that changed for more than the
standard deviation of the initial weights. This can be defined as

$c=\frac{1}{|W|}\sum_{w\in
W}\mathbbm{1}[|w_{i}^{t}-w_{i}|>\beta];\beta=\sigma(\\{|w_{i}^{t}-w_{i}|;w_{i}\in
W\\})\textbf{,}$

where $W$ is a list of weights, $w_{i}^{t}$ is the trained weight $w_{i}$ and
$\sigma(A)$ is the standard deviation of $A$.

  * •

Normalized difference: We take the absolute difference between the end weights
and the initial weights. We then divide these values by the maximum and take
the mean value. This can be defined as

$c=\frac{1}{|W|\cdot\tau}\cdot\sum_{w_{i}\in
W}|w_{i}-w_{i}^{t}|;\tau=\max\limits_{w_{i}\in W}(|w_{i}-w_{i}^{t}|),$

where $W$ is a list of weights and $w_{i}^{t}$ is the trained weight $w_{i}$.

  * •

Distance to a simple classifier: The difference between the accuracy of the
deep learning model and the logistic regression model on the training set.
This can be defined as

$c=\frac{1}{|X|}\bigg{|}\sum_{x_{i}\in
X}\mathbbm{1}[\textrm{DNN}(x_{i})=y_{i}]-\sum_{x_{i}\in
X}\mathbbm{1}[\textrm{LR}(x_{i})=y_{i}]\bigg{|}$

where $X$ is the training set, $y_{i}$ is the ground truth of $x_{i}$,
$\textrm{DNN}(x_{i})$ is the prediction of our deep learning model and
$\textrm{LR}(x_{i})$ is the prediction of our logistic regression model.

  * •

Margin: Percentage of predictions where the difference between the two classes
with the highest score is greater than $0.2$. This can be defined as

$c=\frac{1}{|X|}\cdot\sum_{x_{i}\in
X}\mathbbm{1}[(\textrm{DNN}_{1}(x)-\textrm{DNN}_{2}(x))\geq 0.2],$

where $X$ is the training data and
$\textrm{DNN}_{1}(x_{i}),\textrm{DNN}_{2}(x_{i})$ the prediction scores for
the best and second best class for instance $x_{i}$.

  * •

Learning rate mean: We train the model with different learning rates (0.0001,
0.001, 0.01) three times for five epochs and store accuracies we get after
each epoch. We take the variance of accuracies for each different value of the
learning rate and divide it by 0.002 (we assume variance is less than this).
The final score is the mean value of these variances. This is defined as

$c=\frac{1}{3\cdot 2\cdot
10^{-3}}\cdot\sum_{i=1}^{3}\sigma^{2}(\textrm{lr}_{i}),$

where $\textrm{lr}_{i}$ are the accuracy scores for the $i$-th learning rate
tested and $\sigma^{2}(A)$ is the variance of $A$.

  * •

Learning rate distance: We train the model with different learning rates
(0.0001, 0.001, 0.01) three times for five epochs and store accuracies we get
after each epoch. We take the variance of accuracies for each different value
of the learning rate and divide it by 0.002 (we assume variance is less than
this). The final score is the difference between the variances and the mean
value of variances normalized. This is defined as

$c=\frac{2}{3\cdot 2\cdot
10^{-3}}\cdot\sum_{i=1}^{3}\big{|}\sigma^{2}(\textrm{lr}_{i})-\gamma\big{|};\gamma=\frac{1}{3}\sum_{j=1}^{3}(\sigma^{2}(\textrm{lr}_{j})),$

where $\textrm{lr}_{i}$ are the accuracy scores for the $i$-th learning rate
tested and $\sigma^{2}(A)$ is the variance of $A$.

  * •

Categorical cross entropy: The categorical cross entropy of predictions on the
training set. This is defined as

$c=-\frac{1}{|X|}\cdot\sum_{i=1}^{|X|}\sum_{j=1}^{C}y_{i,j}\cdot\log~{}\hat{y}_{i,j},$

where $X$ is the training set, $C$ is the number of classes, $y_{i,j}$ is the
ground-truth of the prediction $\hat{y}_{i,j}$.

We created the the prediction by summing the transformed scores of this
metrics; $h(x,f)=\sum_{i=i}^{7}t(x_{i},f_{i})$, where $f\in\\{-1,0,1\\}^{7}$
is a vector of metric inclusion and $t$ is the transformation function shown
in the Equation 3.

$t(x_{i},f_{i})=\begin{cases*}x_{i}&if $f_{i}=1$\\\ 1-x_{i}&if $f_{i}=-1$\\\
0&if $f_{i}=0$\end{cases*}$ (3)

This gives us $3^{7}=2187$ different combinations of simple metrics. We chose
the best combination by automatically testing them on the provided dataset and
choosing the one with the highest score. In our tests, the highest score was
achieved by the combination where metrics "Learning rate distance" and "Non-
uniform weights" were given the metric inclusion value $1$, "Learning rate
mean" feature inclusion value $-1$ and all others feature inclusion value $0$.

By observing the scores of the combinations, we saw that those achieving the
best results were very similar. This gives us further confidence that the
chosen combination of metrics gives some insight into the model’s structure
and how it generalizes.

##  3 Conclusion

We introduced our methodology where simple metrics are combined to generate a
complexity score. Because this combination might not be trivial, different
combinations are tested automatically and the best one is chosen. Because of
this, our approach might be a good starting point for creating a good method
for predicting generalization.

During testing we observed that similar combination usually score similarly.
For each metric we can average scores where metric inclusion value is the
same. By comparing these values we can observe how metrics rank and which
inclusion value is the best for a given metric.

The presented methodology has a few advantages. One of these is the modality.
New metrics can easily be added or removed to improve the performance and
better predict generalization. The automatic approach makes it possible to
explore non-trivial combinations of simple metrics and rank them based on the
overall contribution.

The main disadvantage of this approach is that we need to incorporate already
existing metrics into it. The way we implement how scores from different
metrics are merged also limits the possible metrics that can be added, since
they need to give a score between 0 and 1, otherwise one metric would give
more to the final value than the others. This can probably be avoided by using
hyperparameter tuning approaches such Bayesian Optimization. Another
disadvantage of our approach is that adding new metrics without discarding
others might greatly impact time needed to calculate the best combination.

In further work, we want to add new metrics and discard some that did not
perform well. We also want to improve automatic search for the combination of
metrics using more advanced methods for parameter tuning such as Bayesian
Optimization or genetic algorithms.

## References

  * Jiang et al. (2020a) Yiding Jiang, Pierre Foret, Yak, and Scott et al. Yak.  Predicting generalization in deep learning competition documentation, 2020a.  URL https://sites.google.com/view/pgdl2020/.
  * Jiang et al. (2020b) Yiding Jiang, Behnam Neyshabur, Dilip Krishnan, Hossein Mobahi, and Samy Bengio.  Fantastic generalization measures and where to find them.  _international conference on learning representations_ , 2020b.
  * Kendall (1938) M. G. Kendall.  A new measure of rank correlation.  _Biometrika_ , 30(1-2):81–93, 06 1938.  ISSN 0006-3444.  doi: 10.1093/biomet/30.1-2.81.  URL https://doi.org/10.1093/biomet/30.1-2.81.
  * Li and Tian (2018) Dewei Li and Yingjie Tian.  Survey and experimental study on metric learning methods.  _Neural networks_ , 105:447–462, 2018.
  * Neyshabur et al. (2017) Behnam Neyshabur, Srinadh Bhojanapalli, David Mcallester, and Nati Srebro.  Exploring generalization in deep learning.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 30, pages 5947–5956. Curran Associates, Inc., 2017.  URL https://proceedings.neurips.cc/paper/2017/file/10ce03a1ed01077e3e289f3e53c72813-Paper.pdf.


# A convenient implementation of the overlap between arbitrary Hartree-Fock-
Bogoliubov vacua for projection

Zao-Chun Gao<EMAIL_ADDRESS>Qing-Li Hu  Y. S. Chen  China Institute of
Atomic Energy, P.O. Box 275 (10), Beijing 102413, P.R. China

###### Abstract

Overlap between Hartree-Fock-Bogoliubov(HFB) vacua is very important in the
beyond mean-field calculations. However, in the HFB transformation, the $U,V$
matrices are sometimes singular due to the exact emptiness ($v_{i}=0$) or full
occupation ($u_{i}=0$) of some single-particle orbits. This singularity may
cause some problem in evaluating the overlap between HFB vacua through
Pfaffian. We found that this problem can be well avoided by setting those zero
occupation numbers $u_{i},v_{i}$ to some tiny values denoted by
$\varepsilon(>0)$, which numerically satisfies $1+\varepsilon^{2}=1$ (e.g.,
$\varepsilon=10^{-8}$ when using the double precision data type). This
treatment does not change the HFB vacuum state because
$u_{i}^{2},v_{i}^{2}=\varepsilon^{2}$ are numerically zero relative to 1.
Therefore, for arbitrary HFB transformation, we say that the $U,V$ matrices
can always be nonsingular. From this standpoint, we present a new convenient
Pfaffian formula for the overlap between arbitrary HFB vacua, which is
especially suitable for symmetry restoration. Testing calculations have been
performed for this new formula. It turns out that our method is reliable and
accurate in evaluating the overlap between arbitrary HFB vacua.

###### keywords:

Hartree-Fock-Bogoliubov method, beyond mean-field method, Pfaffian

††journal: Physics Letter B

##  1 Introduction

The Hartree-Fock-Bogoliubov (HFB) approximation has been a great success in
understanding interacting many-body quantum systems in all fields of physics.
However, the beyond mean-field effects (e.g., the nuclear vibration and
rotation) are missing in the HFB calculations. Methods that go beyond mean-
field, such as the Generator Coordinate Method(GCM) and the projection method,
are expected to take those missing effects into consideration and present
better description of the many-body quantum system. In the beyond mean-field
calculations, operator matrix elements and overlaps between multi-
quasiparticle HFB states are basic blocks. These matrix elements and overlaps
can be evaluated using the generalized Wick’s theorem (GWT)[1, 2], or
equivalently using Pfaffian [3, 4, 5, 6, 7], or using the compact formula in
Ref.[8]. However, in the efficient calculations (e.g., see [5]), all of the
matrix elements and overlaps require the value of the overlap between HFB
vacua.

Thus, the reliable and accurate evaluation of the overlap between HFB vacua is
very important for the stability and the efficiency of the beyond mean-field
calculations. Especially in cases near to the Egido pole [9], the overlap
between HFB vacua is very tiny, and a small error could lead to a large
uncertainty of the matrix elements. In the past, numerical calculations of the
overlap were performed with the Onishi formula [10]. Unfortunately, the Onishi
formula leaves the sign of the overlap undefined due to the square root of a
determinant. Several efforts have been made to overcome this sign problem [11,
12, 13, 14, 15, 16]. In 2009, Robledo proposed a different overlap formula
with the Pfaffian rather than the determinant [17]. This formula completely
solves the sign problem but requires the inversion of the matrix $U$ in the
Bogoliubov transformation. To avoid the singularity of $U$, the formula for
the limit when several orbits are fully occupied is given in Ref. [18].
Simultaneously, the limit when some orbits are exact empty was also considered
to reduce the computational cost. Meanwhile, various Pfaffian formulae for the
overlap between HFB vacua have been proposed by several authors [3, 6, 7]. In
Ref. [7], the overlap formula does not require the inversion of $U$, but the
empty orbits in the Fock space should be omitted.

In practical calculation, one should first identify the singularity of the
matrices $U$ and $V$ in the Bogoliubov transformation. This can be easily
tested with the Bloch-Messiah theorem (see details in Ref. [19]). The matrices
$U$ and $V$ can be decomposed as $U=D\bar{U}C$ and $V=D^{*}\bar{V}C$. Here,
$D$ and $C$ are unitary matrices. $\bar{U}$ and $\bar{V}$ refer to the BCS-
transformation and are constructed from the occupation numbers $u_{i},v_{i}$
with $0\leq u_{i},v_{i}\leq 1$ and $u_{i}^{2}+v_{i}^{2}=1$ (see Eq. (7.9) and
Eq. (7.12) in Ref.[19]). The limits of fully occupied ($u_{i}=0,v_{i}=1$) and
fully empty ($u_{i}=1,v_{i}=0$) levels have been carefully treated in Refs.
[6, 7, 18] to avoid the collapse of the overlap computation.

However, we note that in most realistic cases the $v_{i}$’s can be extremely
close to 0 or 1 but not exact 0 or 1. Strictly speaking, these levels with
such extreme emptiness or occupation should be considered but may lead to
exotic values ( extremely huge or extremely tiny) of the Pfaffian in the
proposed formulae. What is worse, the Pfaffian values are easily out of the
scope of the double precision data type and cause the computation collapsed.

Careful treatment must be made to avoid such data overflow. In this paper, we
implement an accurate and reliable calculation for the overlap between
arbitrary HFB vacua in a unified way. For the cases of ($u_{i}=1,v_{i}=0$) and
($u_{i}=0,v_{i}=1$), we treat them as the cases of
($u_{i}=1,v_{i}=\varepsilon$) and ($u_{i}=\varepsilon,v_{i}=1$), respectively.
The tiny quantity $\varepsilon>0$ is chosen such that $\varepsilon^{2}$ should
be numerical zero relative to 1 in the practical calculation. In other words,
$\varepsilon$ should numerically satisfy $1+\varepsilon^{2}=1$. Under this
condition, $\varepsilon$ may be chosen as large as possible so that the
calculated Pfaffian values are not necessarily too huge or too tiny. For
instance, one can choose $\varepsilon=10^{-8}$ when using double precision.
Because $v_{i}^{2}(u_{i}^{2})=\varepsilon^{2}$ is actually zero relative to
$u_{i}^{2}(v_{i}^{2})=1$ in practical calculations, this treatment does not
change the HFB vacuum at all. Therefore, without losing the generality, we
assume that all levels in the Fock space are partly occupied, but some of
their $u_{i},v_{i}$ values are allowed to be extremely close to 0 or 1.
Ideally, $U,V$ are nonsingular in our assumption, and we can derive a new
formula for the overlap between the HFB vacua based on the work of Bertsch and
Robledo [7]. This formula is especially convenient for the symmetry
restoration. Numerical calculations have been carried out for heavy nuclear
system to test the precision of the new formula by comparing with the Onishi
formula.

In section 2, the formalism of the new overlap formula is given. Section 3
provides an example of numerical calculation. A summary is given in section 4.

##  2 The overlap between the HFB vacua

We denote $\hat{c}^{\dagger}_{i}$ and $\hat{c}_{i}$ as the creation and
annihilation operators defined in an $M$-dimensional Fock-space. The Hartree-
Fock-Bogoliubov(HFB) transformation is

$\left(\begin{array}[]{c}\hat{\beta}\\\
\hat{\beta}^{\dagger}\end{array}\right)=\left(\begin{array}[]{cc}U^{\dagger}&V^{\dagger}\\\
V^{T}&U^{T}\end{array}\right)\left(\begin{array}[]{c}\hat{c}\\\
\hat{c}^{\dagger}\end{array}\right).$ (1)

Here, we assume $U$ and $V$ are nonsingular matrices, and their shapes are
$M\times M$. The HFB vacuum (unnormalized) can be written as

$\displaystyle|\phi\rangle=\hat{\beta}_{1}\hat{\beta}_{2}...\hat{\beta}_{M}|-\rangle,$
(2)

where $|-\rangle$ is the true vacuum. By definition, one has

$\displaystyle\hat{\beta}_{i}|\phi\rangle=0\quad\mathrm{for}\quad 1\leq i\leq
M.$ (3)

The second HFB vacuum $|\phi^{\prime}\rangle$ is defined in the same way, but
the prime, ‘′’, is attached to the corresponding symbols to show difference.

The overlap between $|\phi\rangle$ and$|\phi^{\prime}\rangle$ is given by

$\displaystyle\langle\phi|\phi^{\prime}\rangle$ $\displaystyle=$
$\displaystyle\langle-|\hat{\beta}^{\dagger}_{M}\hat{\beta}^{\dagger}_{M-1}...\hat{\beta}^{\dagger}_{1}\hat{\beta}^{\prime}_{1}\hat{\beta}^{\prime}_{2}...\hat{\beta}^{\prime}_{M}|-\rangle$
(4) $\displaystyle=$ $\displaystyle
s_{M}\langle-|\hat{\beta}^{\dagger}_{1}\hat{\beta}^{\dagger}_{2}...\hat{\beta}^{\dagger}_{M}\hat{\beta}^{\prime}_{1}\hat{\beta}^{\prime}_{2}...\hat{\beta}^{\prime}_{M}|-\rangle,$

where, $s_{M}=(-1)^{[M(M-1)/2]}$. If $M$ is even, $s_{M}=(-1)^{M/2}$.
Following the technique of Bertsch and Robledo [7], one can obtain

$\displaystyle\langle\phi|\phi^{\prime}\rangle$ $\displaystyle=$
$\displaystyle
s_{M}\mathrm{pf}\left(\begin{array}[]{cc}V^{T}U&V^{T}V^{\prime*}\\\
-V^{\prime\dagger}V&U^{\prime\dagger}V^{\prime*}\end{array}\right).$ (7)

The shape of the matrix in Eq. (7) is $2M\times 2M$, and no empty levels are
omitted. For the norm overlap $\langle\phi|\phi\rangle$, it is real and
positive. From Eq.(7) and the Bloch-Messiah theorem, one can get

$\displaystyle\langle\phi|\phi\rangle=s_{M}\mathrm{pf}\left(\begin{array}[]{cc}V^{T}U&V^{T}V^{*}\\\
-V^{\dagger}V&U^{\dagger}V^{*}\end{array}\right)=\prod_{i=1}^{M/2}v_{i}^{2}.$
(10)

Denoting $\prod_{i=1}^{M/2}v_{i}$ by $\mathfrak{N}$, the normalized quasi-
particle vacuum, $|\psi\rangle$, can be written as

$\displaystyle|\psi\rangle=\frac{|\phi\rangle}{\mathfrak{N}}.$ (11)

Then, one finds that

$\displaystyle\langle\psi|\psi^{\prime}\rangle=\frac{s_{M}}{\mathfrak{NN^{\prime}}}\mathrm{pf}\left(\begin{array}[]{cc}V^{T}U&V^{T}V^{\prime*}\\\
-V^{\prime\dagger}V&U^{\prime\dagger}V^{\prime*}\end{array}\right).$ (14)

In the symmetry restoration, the general rotational operator, involving the
spin and particle number projection, may be written as

$\displaystyle\hat{\mathbb{R}}(\Xi)=\hat{R}(\Omega)e^{-i\hat{N}\phi_{n}}e^{-i\hat{Z}\phi_{p}},$
(15)

where $\hat{R}(\Omega)$ is the rotation operator, and $\Omega$ refers to the
three Euler angles $\alpha,\beta,\gamma$. $e^{-i\hat{N}\phi_{n}}$ and
$e^{-i\hat{Z}\phi_{p}}$ are ‘gauge’ rotational operators induced by the
neutron and proton number projection. $\hat{N}$ and $\hat{Z}$ are neutron and
proton number operators, respectively. $\phi_{n}$ and $\phi_{p}$ are ”gauge”
angles for neutron and proton, respectively. $\Xi$ refers to
$(\Omega,\phi_{n},\phi_{p})$. The matrix element
$\langle\psi|\hat{\mathbb{R}}(\Xi)|\psi^{\prime}\rangle$ needs to be
calculated. Let’s define the general rotation transformation for symmetry
restoration,

$\displaystyle\hat{\mathbb{R}}(\Xi)\left(\begin{array}[]{cc}\hat{c}\\\
\hat{c}^{\dagger}\end{array}\right)\hat{\mathbb{R}}^{\dagger}(\Xi)=\left(\begin{array}[]{cc}\mathbb{\mathbb{D}}^{\dagger}(\Xi)&0\\\
0&\mathbb{\mathbb{D}}^{T}(\Xi)\end{array}\right)\left(\begin{array}[]{cc}\hat{c}\\\
\hat{c}^{\dagger}\end{array}\right),$ (22)

where $\mathbb{D}_{ij}(\Xi)=\langle i|\hat{\mathbb{R}}(\Xi)|j\rangle$, and
$|i(j)\rangle=\hat{c}^{\dagger}_{i(j)}|-\rangle$. The $\mathbb{D}(\Xi)$ matrix
has the dimension $M\times M$. One can get

$\displaystyle\hat{\mathbb{R}}(\Xi)\left(\begin{array}[]{cc}\hat{\beta}^{\prime}\\\
\hat{\beta}^{\prime\dagger}\end{array}\right)\hat{\mathbb{R}}^{\dagger}(\Xi)=\mathcal{D}(\Xi)\left(\begin{array}[]{cc}\hat{c}\\\
\hat{c}^{\dagger}\end{array}\right),$ (27)

where

$\displaystyle\mathcal{D}(\Xi)=\left(\begin{array}[]{cc}[\mathbb{D}(\Xi)U^{\prime}]^{\dagger}&[\mathbb{D}^{*}(\Xi)V^{\prime}]^{\dagger}\\\
{[\mathbb{D}^{*}(\Xi)V^{\prime}]^{T}}&{[\mathbb{D}(\Xi)U^{\prime}]^{T}}\end{array}\right).$
(30)

By comparing Eq.(27) with Eq.(1), one can obtain the rotated overlap by
replacing $U^{\prime}$ and $V^{\prime}$ in Eq.(14) with
$\mathbb{D}(\Xi)U^{\prime}$ and $\mathbb{D}^{*}(\Xi)V^{\prime}$, respectively.
Thus

$\displaystyle\mathbb{N}_{\mathrm{pf}}(\Xi)=\langle\psi|\hat{\mathbb{R}}(\Xi)|\psi^{\prime}\rangle=\frac{s_{M}}{\mathfrak{NN^{\prime}}}\mathrm{pf}[\mathcal{M}(\Xi)],$
(31)

where

$\mathcal{M}(\Xi)=\left(\begin{array}[]{cc}V^{T}U&V^{T}\mathbb{D}(\Xi)V^{\prime*}\\\
-V^{\prime\dagger}\mathbb{D}^{T}(\Xi)V&U^{\prime\dagger}V^{\prime*}\end{array}\right).$
(32)

This formula is essentially the same as the one proposed by Bertsch and
Robledo [7], but we will transform it into a new form. Supposing that there is
a $\Xi_{0}$ satisfying $\mathbb{N}_{\mathrm{pf}}(\Xi_{0})\neq 0$, we have

$\displaystyle\frac{\mathbb{N}_{\mathrm{pf}}(\Xi)}{\mathbb{N}_{\mathrm{pf}}(\Xi_{0})}$
$\displaystyle=$
$\displaystyle\frac{\mathrm{pf}[\mathcal{M}(\Xi)]}{\mathrm{pf}[\mathcal{M}(\Xi_{0})]}=\frac{\mathrm{pf}\left[P\mathcal{M}(\Xi)P^{T}\right]}{\mathrm{pf}\left[P\mathcal{M}(\Xi_{0})P^{T}\right]}$
(33) $\displaystyle=$
$\displaystyle\frac{\mathrm{pf}[\mathcal{W}(\Xi)]}{\mathrm{pf}[\mathcal{W}(\Xi_{0})]},$

where

$\displaystyle\mathcal{W}(\Xi)=\left(\begin{array}[]{cc}[U^{\prime}V^{\prime-1}]^{\dagger}&-\mathbb{D}^{T}(\Xi)\\\
\mathbb{D}(\Xi)&UV^{-1}\end{array}\right),$ (36)

and $P$ is

$\displaystyle P=\left(\begin{array}[]{cc}0&(V^{\prime\dagger})^{-1}\\\
(V^{T})^{-1}&0\end{array}\right).$ (39)

Therefore, one can get

$\displaystyle{\mathbb{N}_{\mathrm{pf}}(\Xi)}=\mathcal{C}{\mathrm{pf}[\mathcal{W}(\Xi)]},$
(40)

where, the coefficient $\mathcal{C}$ is actually independent of $\Xi_{0}$, and
can be written as

$\displaystyle\mathcal{C}=\frac{\mathbb{N}_{\mathrm{pf}}(\Xi_{0})}{\mathrm{pf}[\mathcal{W}(\Xi_{0})]}=\frac{s_{M}}{\mathfrak{NN^{\prime}}\mathrm{det}P}={s_{M}\Delta\mathfrak{NN^{\prime}}}.$
(41)

Here, $\Delta$ is a phase determined by

$\displaystyle\Delta=\mathrm{det}D^{*}\mathrm{det}D^{\prime}\mathrm{det}C\mathrm{det}C^{\prime*}.$
(42)

In Eq.(40), we have used the Bloch-Messiah theorem and the following equation

$\displaystyle\mathrm{det}P=\mathrm{det}[(V^{\prime\dagger})^{-1}]\mathrm{det}[(V^{T})^{-1}].$
(43)

Eq.(40) looks more convenient to be implemented and may save some computing
time in contrast to Eq.(31), where extra evaluation of
$V^{T}\mathbb{D}(\Xi)V^{\prime*}$ is required for each mesh point in the
integral of projection.

For comparison, let us present a brief introduction of the overlap of the
Onishi formula [10]. The unitary transformation of the quasi-particles under
rotation $\hat{\mathbb{R}}(\Xi)$ can be written as

$\displaystyle\hat{\mathbb{R}}(\Xi)\left(\begin{array}[]{cc}\hat{\beta}^{\prime}\\\
\hat{\beta}^{\prime\dagger}\end{array}\right)\hat{\mathbb{R}}^{\dagger}(\Xi)=\left(\begin{array}[]{cc}\mathbb{X}(\Xi)&\mathbb{Y}(\Xi)\\\
\mathbb{Y}^{*}(\Xi)&\mathbb{X}^{*}(\Xi)\end{array}\right)\left(\begin{array}[]{cc}\hat{\beta}\\\
\hat{\beta}^{\dagger}\end{array}\right),$ (50)

where

$\displaystyle\mathbb{X}(\Xi)$ $\displaystyle=$ $\displaystyle
U^{\prime\dagger}\mathbb{D}^{\dagger}(\Xi)U+V^{\prime\dagger}\mathbb{D}^{T}(\Xi)V,$
$\displaystyle\mathbb{Y}(\Xi)$ $\displaystyle=$ $\displaystyle
U^{\prime\dagger}\mathbb{D}^{\dagger}(\Xi)V^{*}+V^{\prime\dagger}\mathbb{D}^{T}(\Xi)U^{*}.$
(51)

The Onishi formula is then expressed as (see Ref.[20]),

$\displaystyle\mathbb{N}_{\mathrm{Onishi}}(\Xi)=\langle\psi|\hat{\mathbb{R}}(\Xi)|\psi^{\prime}\rangle$
(52) $\displaystyle=$
$\displaystyle(\pm)\sqrt{\mathrm{det}[\mathbb{X}(\Xi)]}e^{-i(M_{n}\phi_{n}+M_{p}\phi_{p})/2},$

where $M_{n}$ and $M_{p}$ are the numbers of neutron and proton orbits in the
Fock space, respectively. The value of $\mathrm{det}[\mathbb{X}(\Xi)]$ is a
complex number, and the sign of the square root is left undefined. Extra
efforts must be made to determine the sign before the application of the
Onishi formula. For instance, in the Projected Shell Model [21] without
particle number projection, the overlap between the BCS vacua is real and
positive, thus there is no sign ambiguity and the Onishi formula works.

##  3 Numerical test of the overlap formulae

Although the sign problem is solved in Eq.(31) and Eq.(40), one can imagine
that $\mathfrak{N}$, $\mathfrak{N}^{\prime}$ are extremely tiny numbers by
definition. Thus $\mathrm{pf}[\mathcal{M}(\Xi)]$ is also very tiny, but
$\mathrm{pf}[\mathcal{W}(\Xi)]$ should be huge. Numerical accuracy of Eqs.
(31) and (40) needs to be carefully tested. It is believed that the Onishi
formula is accurate except for its undetermined sign. So, it is helpful to
compare the numerical values of the overlaps using Eq.(31), Eq.(40) and
Eq.(52).

To demonstrate the accuracy and the reliability of the Eqs. (31) and (40),
numerical calculations are performed for the typical example of the deformed
heavy nucleus 226Th. For projection, we should take
$|\psi\rangle=|\psi^{\prime}\rangle$, and then $\Delta=1$.

The $U,V$ matrices are obtained from the Nilsson+BCS method. The single
particle levels are generated from the Nilsson Hamiltonian with the standard
parameters [22]. The single-particle model space contains 5 neutron major
shells with $N=$4-8 and 5 proton major shells with $N=$3-7, i.e., the Fock
space has 145 neutron levels ($M_{n}=290$) and 110 proton levels
($M_{p}=220$). The numbers of the active neutrons and protons are 96 and 70,
respectively. The quadrupole deformation is taken to be $\epsilon_{2}=0.2$.
Here, we only consider the axial symmetry for simplicity.

In the no pairing case, the BCS vacuum becomes a pure slater determinant,
which is a challenge for Eq.(40) because all $v_{i}$’s above the Fermi surface
are zero. Consequently, $\mathfrak{N}=0$ and $\mathcal{W}(\Xi)$ is meaningless
due to the singularity of $V$. Here, we use the double precision data type and
set $v_{i}=\varepsilon=10^{-8}$ for those $v_{i}=0$ orbits to avoid the
collapse of calculation. Therefore we have

$\displaystyle\langle\phi_{n}|\phi_{n}\rangle=(10^{-16})^{\frac{290-96}{2}}=10^{-1552},$
$\displaystyle\langle\phi_{p}|\phi_{p}\rangle=(10^{-16})^{\frac{220-70}{2}}=10^{-1200},$

where, $|\phi_{n}\rangle$ and $|\phi_{p}\rangle$ are BCS vacua for neutrons
and protons, respectively, and
$|\phi\rangle=|\phi_{n}\rangle|\phi_{p}\rangle$. The tiny numbers $10^{-1552}$
and $10^{-1200}$ are too far out of the scope of the double precision data
($\sim 10^{\pm 307}$). To avoid the data overflow, we multiply the tiny
variable by $10^{200}$ several times until the scaled absolute value falls
into the interval $[10^{-200},10^{200}]$. In other words, we use a number $y$
and an integer number $k$ to express a tiny number $x$ through
$x=y\times(10^{-200})^{k}$. If $x$ is a huge number, then $k$ is negative.

Figure 1: (Color online) Overlaps of the ground state neutron slater
determinant for 226Th as functions of $\phi_{n}$ with Euler angles
$\alpha=\gamma=0^{\circ},\beta=10^{\circ}$, calculated with present formula
[Eq.(40)] and the Onishi formula [Eq.(52) with ‘+’ sign].
Re$[\mathbb{N}(\phi_{n})]$ and Im$[\mathbb{N}(\phi_{n})]$ are the real and
imaginary parts of the overlap.

However, for the Onishi formula of Eq.(52), we do not need to change $v_{i}=0$
to $v_{i}=\varepsilon$. The overlaps for the neutron part, calculated with
Eq.(40) and Eq.(52), are compared in Fig.1. The curves of Eq.(40) are
continuous, but the sign uncertainty of Eq.(52) causes the discontinuity.
However, if one copies the sign of Eq.(40) to Eq.(52), one can compare
numerical difference between Eq.(40) and Eq.(52) using the following quantity,
$R$,

$\displaystyle
R=\left|\frac{\mathbb{N}_{\mathrm{Onishi}}(\phi_{n})}{\mathbb{N}_{\mathrm{pf}}(\phi_{n})}-1\right|.$
(53)

In all calculations, we found that $R<10^{-12}$ with double precision. This
confirms that a small change of $v_{i}$ from zero to $\varepsilon$ almost does
not affect numerical accuracy. However, it is crucial to keep Eq.(40) valid.
Yet notice that $\mathbb{N}_{\mathrm{pf}}(\Xi)$ in Eq.(40) is obtained from a
product of tiny and giant numbers. The same calculations have also been done
with Eq.(31), and we also get $R<10^{-12}$. Thus we have presented an
alternative way of using Eq.(31), where we set $v_{i}=\varepsilon$ for those
empty orbits rather than omitting them[7].

Once the overlap is available, it is straightforward to perform the symmetry
restoration. The deformed BCS vacuum of 226Th has been projected onto good
particle number and spin. Therefore, one can test how precise the numerical
calculations with Eq.(40) satisfy

$\displaystyle\sum_{N,Z,I}\langle\psi|\hat{P}^{N}\hat{P}^{Z}\hat{P}^{I}_{00}|\psi\rangle=1,$
(54)

where $\hat{P}^{N}$, $\hat{P}^{Z}$, and $\hat{P}^{I}_{MK}$ are neutron-number,
proton-number, and spin projection operators, respectively. For the above
vacuum state without pairing (i.e. the ground state slater determinant), the
particle numbers of both neutrons and protons are good. Indeed, our particle
number projection (using 16 mesh points in the integral) shows that
$\langle\psi_{n}|\hat{P}^{N}|\psi_{n}\rangle=1$ ($N=96$), or 0 ($N\neq 96$)
with numerical errors less than $10^{-13}$. Calculations for the protons also
have the same accuracy. This again shows the reliability of Eq. (40). Angular
momentum projection is also performed on the same state in addition to the
particle number projection. The amplitude of
$\langle\psi|\hat{P}^{N}\hat{P}^{Z}\hat{P}^{I}_{00}|\psi\rangle$ with
($N=96,Z=70$) is plotted as a function of spin $I$ in Fig.2. In the integral
of the spin projection, 100 mesh points are taken, and the range of spin is
$0\leq I\leq 70$, and we indeed reproduced Eq.(54) with numerical error around
$10^{-12}$.

Figure 2:  The amplitude of projection,
$\langle\psi|\hat{P}^{N}\hat{P}^{Z}\hat{P}^{I}_{00}|\psi\rangle$, as a
function of spin $I$ at $N=96$ and $Z=70$ using Eq.(40). $|\psi\rangle$ is the
axially deformed BCS vacuum but without pairing.

We also have tested Eq.(40) in the projection of the triaxially deformed
vacuum with normal pairing, which seems more convenient to use Eq.(40). With
the present method, similar accuracy has also been achieved.

##  4 Summary

Following the strategy of Bertsch and Robledo [7], we have proposed a new
formula of the overlap between HFB vacua by using the Pfaffian identity and
assuming that the inverse of the $V$ matrix exists. This formula is especially
convenient and efficient in the symmetry restoration, and has the same high
accuracy as the Onishi formula as well as the correct sign. The reliability of
the present formula has been tested by carrying out the calculations of the
overlap and the quantum number projection for the heavy nucleus 226Th. In the
testing calculations, one has to be faced with two numerical problems: (1) The
extreme (huge or tiny) quantities are certainly encountered, and we have
properly treated this situation to avoid data overflow (see the text). (2) For
those empty orbits with $v_{i}=0$, which make Eq.(40) invalid, one can change
$v_{i}$ to a small quantity $\varepsilon(>0)$ to avoid the singularity of $V$
matrix. It turns out that such treatments work very well. Testing calculations
have confirmed that the present formula is even applicable to the pure slater
determinant without losing the numerical accuracy. Thus it is promising that
Eq.(40) may be applicable in evaluating the overlap between arbitrary HFB
vacua.

Acknowledgements Z. G. thanks Prof. Y. Sun and Dr. F. Q. Chen for the
stimulating and fruitful discussions. The authors acknowledge support from the
National Natural Science Foundation of China under Contract Nos. 11175258,
11021504 and 11275068.

## References

  * [1] R. Balian and E. Brezin, Nuovo Cimento B 64 (1969) 37.
  * [2] K. Hara, S. Iwasaki, Nucl. Phys. A 332 (1979) 61.
  * [3] M. Oi, T. Mizusaki, Phys. Lett. B 707 (2012) 305.
  * [4] T. Mizusaki, M. Oi, Phys. Lett. B 715 (2012) 219.
  * [5] T. Mizusaki, M. Oi, Fang-Qi Chen, Yang Sun, Phys. Lett. B 725 (2013) 175.
  * [6] B. Avez, M. Bender, Phys. Rev. C. 85 (2012) 034325.
  * [7] G.F. Bertsch and L.M. Robledo, Phys. Rev. Lett. 108 (2012) 042505.
  * [8] S. Perez-Martin and L.M. Robledo, Phys. Rev. C 76 (2007) 064314.
  * [9] M. Anguiano, J.L. Egido, L.M. Robledo, Nucl. Phys. 696 (2001) 467.
  * [10] N. Onishi and S. Yoshida, Nucl. Phys. 80 (1966) 367.
  * [11] K. Neergård and E. Wüst, Nucl. Phys. A 402 (1983) 311.
  * [12] Q. Haider and D. Gogny, J. Phys. G: Nucl. Part. Phys. 18 (1992) 993.
  * [13] F. Dönau, Phys. Rev. C 58 (1998) 872.
  * [14] M. Bender and Paul-Henri Heenen, Phys. Rev. C 78 (2008) 024309.
  * [15] K. Hara, A. Hayashi and P. Ring, Nucl. Phys. A 385 (1982) 14.
  * [16] M. Oi and N. Tajima, Phys. Lett. B 606 (2005) 43.
  * [17] L.M. Robledo, Phys. Rev. C 79 (2009) 021302(R).
  * [18] L.M. Robledo, Phys. Rev. C 84 (2011) 014307.
  * [19] P. Ring and P. Schuck, The Nuclear Many-Body Problem, Springer-Verlag, 1980.
  * [20] K.W. Schmid, Prog. Part. Nucl. Phys. 52 (2004) 565.
  * [21] K. Hara and Y. Sun, Int. J. Mod. Phys. E 4 (1995) 637.
  * [22] T. Bengtsson and I. Ragnarsson, Nucl. Phys. A 436 (1985) 14.