# Can the muon anomalous magnetic moment and the $B$ anomalies be
simultaneously explained in a minimally extended $Z^{\prime}$ model?

H. Lee<EMAIL_ADDRESS>National Centre for Nuclear Research,
Pasteura 7, 02-093 Warsaw, Poland  A. E. Cárcamo Hernández
<EMAIL_ADDRESS>Departamento de Física, Universidad Técnica Federico
Santa María,
Casilla 110-V, Valparaíso, Chile  Centro Científico-Tecnológico de Valparaíso,
Casilla 110-V, Valparaíso, Chile  Millennium Institute for Subatomic Physics
at High-Energy Frontier (SAPHIR), Fernández Concha 700, Santiago, Chile

(August 27, 2024)

###### Abstract

We analyze the implications of an extended 2HDM theory in quark masses and
mixings, electroweak precision tests, charged lepton flavor violating (CLFV)
decays, neutrino trident production, B meson, $g-2$ muon and W anomalies. In
the extended 2HDM theory considered in this work, the SM gauge symmetry is
enlarged by the local $U(1)^{\prime}$ symmetry, whereas its scalar and fermion
sectors are augmented by the inclusion of a gauge singlet scalar and
vectorlike fermions, respectively. Within the framework of this model, we
determine the allowed ranges of the $Z^{\prime}$ gauge boson mass consistent
with the constraints arising from neutrino trident production, $R_{K^{*}}$
anomaly, $B_{s}$ meson oscillation. We conclude the $R_{K^{*}}$ and muon $g-2$
anomalies cannot be simultaneously explained by the same new physics under all
the constraints. Besides that, we found that the new physics sources, the
$Z^{\prime}$ and non SM scalars cannot be connected via the scalar potential.
Furthermore, we have found that the model under consideration can successfully
comply with the constraints arising from oblique $T,S,U$ parameters as well as
with the $W$ mass and $g-2$ muon anomalies. In particular we found that the
muon $g-2$ anomaly, whose dominant contribution arises from the one loop level
scalar exchange, can be successfully accommodated within the $2\sigma$
experimentally allowed range.

##  I INTRODUCTION

The Standard Model (SM) has explained many observables with a great
consistency, so being established as a cornerstone to understand our nature.
In spite of its great successes, there are a few critical limitations in the
SM and one of them is muon and electron anomalies. Until the new result for
muon and electron $g-2$ were released in 2020, explaining the two anomalies
simultaneously had been regarded as a potential new physics signal due to
their relative different overall sign and magnitude. The old experimental
constraints for the muon (the Brookhaven E821 experiment at BNL
Muong-2:2006rrc ) and electron $g-2$ (the Berkeley 2018 experiment
Parker:2018vye ) at $1\sigma$ read:

$\begin{split}\Delta
a_{\mu}&=a_{\mu}^{\mathop{\mathrm{e}xp}}-a_{\mu}^{\mathop{\mathrm{S}M}}=\left(26.1\pm
8.0\right)\times 10^{-10},\\\ \Delta
a_{e}&=a_{e}^{\mathop{\mathrm{e}xp}}-a_{e}^{\mathop{\mathrm{S}M}}=\left(-0.88\pm
0.36\right)\times 10^{-12}.\end{split}$ (1)

The old muon $g-2$ experimental constraint reports $3.3\sigma$, whereas the
old electron $g-2$ reports $-2.4\sigma$ deviation from the SM prediction. One
of our works Hernandez:2021tii  tried to explain both anomalies in a unified
way, based on the old experimental constraints, and showed that the SM $W$
exchange at one-loop level is not possible, however the CP-even and -odd
scalar exchange at one-loop level allows to successfully explain both
anomalies within the $1\sigma$ experimentally allowed range. This situation
has changed significantly when the new experimental constraints for the muon
anomalous magnetic moment was released in 2020. The renewed muon $g-2$
constraint is given by FNAL Muong-2:2021ojo

$\begin{split}\Delta
a_{\mu}&=a_{\mu}^{\mathop{\mathrm{e}xp}}-a_{\mu}^{\mathop{\mathrm{S}M}}=\left(25.1\pm
5.9\right)\times 10^{-10},\end{split}$ (2)

which reports $4.2\sigma$ deviation from the SM prediction Aoyama:2020ynm ;
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 .
Around the same time, the renewed electron $g-2$ constraint is given by the
LKB 2020 experiment Morel:2020dww

$\begin{split}\Delta
a_{e}&=a_{e}^{\mathop{\mathrm{e}xp}}-a_{e}^{\mathop{\mathrm{S}M}}=\left(0.48\pm
0.30\right)\times 10^{-12},\end{split}$ (3)

which reports $1.6\sigma$ SM deviation. Comparing the old constraints for the
muon and electron $g-2$ with the renewed results, one can know that the muon
$g-2$ gets more important whereas the electron $g-2$ gets less important due
to the reduced room for physics beyond the SM. Especially, the muon $g-2$ with
$4.2\sigma$ SM deviation can be studied in the context of $Z^{\prime}$ models
Altmannshofer:2016brv ; CarcamoHernandez:2019ydc ; Belanger:2015nma ;
CarcamoHernandez:2019xkb ; Allanach:2015gkd ; Raby:2017igl ; Kawamura:2019rth
; Kawamura:2019hxp ; Kawamura:2021ygg  or of leptoquarks Cheung:2001ip ;
ColuccioLeskow:2016dox ; Crivellin:2020tsz  or of an augmented scalar sector
with/without vectorlike fermions Arnan:2019uhr ; Crivellin:2018qmi ;
Crivellin:2021rbq ; Hernandez:2021tii . Instead of considering the electron
$g-2$ for new physics, the $B$ anomalies start getting more attention and
interest by the particle physics community due to its sizeable SM deviation
corresponding to $3-4\sigma$ at most. The $B$ anomalies we consider in this
work are the $R_{K}$ anomaly which has been studied by BABAR BaBar:2012mrf ,
Belle BELLE:2019xld , LHC LHCb:2014vgu ; LHCb:2019hip , and LHCb LHCb:2021trn
and $R_{K^{*}}$ by LHCb LHCb:2017avl . Specially, the $R_{K}$ anomaly which
has been studied by LHCb LHCb:2021trn  reports a recent measurement in the
dilepton mass-squared range
$1.1\mathop{\mathrm{G}eV}^{2}<q^{2}<6.0\mathop{\mathrm{G}eV}^{2}$

$R_{K}^{\left[1.1,6.0\right]}=\frac{\mathop{\mathrm{B}R}\left(B\rightarrow
K\mu^{+}\mu^{-}\right)}{\mathop{\mathrm{B}R}\left(B\rightarrow
Ke^{+}e^{-}\right)}=0.846_{-0.041}^{+0.044},$ (4)

which corresponds for $2.4-2.5\sigma$ SM deviation and appears to indicate the
breaking of SM lepton universality. The breaking of SM lepton universality can
also be seen by the $R_{K^{*}}$ anomaly, which has been explored by LHCb
LHCb:2017avl , with increased SM deviation of $3.1\sigma$

$R_{K^{*}}^{\left[1.1,6.0\right]}=\frac{\mathop{\mathrm{B}R}\left(B\rightarrow
K^{*}\mu^{+}\mu^{-}\right)}{\mathop{\mathrm{B}R}\left(B\rightarrow
K^{*}e^{+}e^{-}\right)}=0.690_{-0.120}^{+0.160}.$ (5)

The $B$ anomaly can be accessed at tree level via the virtual exchange of a
neutral and massive $Z^{\prime}$ gauge boson Crivellin:2015mga ;
Crivellin:2015lwa ; Chiang:2017hlj ; King:2017anf ; King:2018fcg ;
Falkowski:2018dsl , or via the leptoquark exchange Becirevic:2017jtw ;
deMedeirosVarzielas:2018bcy ; DeMedeirosVarzielas:2019nob ; King:2021jeo  for
the purpose of enhancing its theoretical sensitivity so that it can be seen by
close future experiments. The sizeable anomalies, muon $g-2$ and $B$
anomalies, are an important clue in the search for physics beyond the SM and
simultaneously explaining both of them in a well-motivated BSM theory can be
an intuitive study for new physics. These attempts to consider both anomalies
in an unified way can be performed via $Z^{\prime}$ gauge boson
Navarro:2021sfb  or leptoquark.

We start from this question: is it possible to explain the muon $g-2$ and $B$
anomalies in an extended $Z^{\prime}$ model while fullfilling all the well-
known SM constraints? Answering to this question is not easy as there can
exist very diverse $Z^{\prime}$ models and one of them, like the one of Ref.
Navarro:2021sfb  which has been a good motivation for this work and a
guideline, tells that it is possible to explain both anomalies in a simplified
fermionic $Z^{\prime}$ model with a fewer number of parameters. However, the
simplified $Z^{\prime}$ model assumes no direct mixing between the SM
fermions, so not considering CLFV constraints such as
$\tau\rightarrow\mu\gamma$ and $\tau\rightarrow 3\mu$, which appears in our
work as one of the most significant constraints, therefore drawing an opposite
conclusion that the muon $g-2$ and $R_{K^{*}}$ anomalies cannot be explained
by the same new physics as there is no overlapped $Z^{\prime}$ mass range in
the model under consideration.

In this work, what we want to achieve is to confirm whether the muon anomalous
magnetic moment $g-2$ and $B$ anomalies can be explained in a unified way,
taking into account many SM constraints as possible, in an extended
$Z^{\prime}$ model. The extended $Z^{\prime}$ model features both a fermionic
$Z^{\prime}$ model CarcamoHernandez:2019ydc ; Belanger:2015nma ;
CarcamoHernandez:2019xkb ; Allanach:2015gkd ; Raby:2017igl ; Kawamura:2019rth
; Kawamura:2019hxp ; Kawamura:2021ygg ; Navarro:2021sfb  and an extend 2HDM
model Hernandez:2021tii ; CarcamoHernandez:2021yev . To be more specific, in
the extended 2HDM theory considered in this work, the SM gauge symmetry is
enhanced by the inclusion of the $U(1)^{\prime}$ local symmetry. Its scalar
sector is the same as the scalar sector of the extended 2HDM theory of
Hernandez:2021tii , thus allowing to simultaneously accommodate both muon and
electron $g-2$ anomalies via the non-SM scalar exchange as in
Hernandez:2021tii . However we consider a $U(1)^{\prime}$ local symmetry
instead of the global $U(1)^{\prime}$ of the model of Hernandez:2021tii . It
is worth mentioning that in the model considered in Hernandez:2021tii , the
CKM mixing matrix can be explained to a good approximation
CarcamoHernandez:2021yev , predicting the mass range of vectorlike charged
leptons as well as vectorlike quarks by charged lepton flavor violation (CLFV)
decays and flavor changing neutral currents (FCNCs). Given that in the
extended 2HDM theory considered in this work, a local $U(1)^{\prime}$ symmetry
is considered, its gauge sector includes a neutral massive fermionic
$Z^{\prime}$ gauge boson. Therefore, the BSM model under consideration has two
new physics sources, which are $Z^{\prime}$ gauge boson and the non-SM scalars
(CP-even and -odd scalars), and then we try to answer whether the muon anomaly
and $B$ anomalies can be simultaneously explained in the BSM model.

For the investigation, the SM constraints that will be analyzed in this work
are:

  * •

The SM lepton sector : Muon $g-2$, Neutrino trident process, CLFV constraint
$\tau\rightarrow\mu\gamma$ and $\tau\rightarrow 3\mu$

  * •

The SM quark sector : $R_{K(K^{*})}$ anomaly, $B_{s}$ meson oscillation,
Collider constraint, The CKM mixing matrix

  * •

Electroweak precision data corresponding to the oblique $T$, $S$ and $U$
parameters as well as the SM $W$ mass anomaly and the tree level vacuum
stability of the scalar potential.

The SM fermion sector in this BSM model is extended by a complete vectorlike
family and one vectorlike family contains two fermionic seesaw mediators that
will trigger an seesaw mechanism to produce the masses of the second and third
generation of the SM fermions Hernandez:2021tii ; CarcamoHernandez:2021yev .
And that is why the SM constraints shown above consist of the second and third
generation of the SM fermions. Adding one more vectorlike family can give rise
to two extra fermionic mediators, allowing the whole generation of the SM
fermions to acquire masses, however this approach makes our analysis much more
complicated since we want to diagonalize the mass matrices in a very concise,
accurate, and complete way without any assumptions to increase predictability
of the model. For this reason, we restrict our attention to the second and
third generation of SM fermions. On top of that, it is worth mentioning that
this BSM model is a simple and economical BSM model as it can explain well the
SM mixing parameters as well as the SM fermion masses, thus allowing a
comprehensive and complete analysis of the anomalies. First of all, we try to
determine the highest order of the $Z^{\prime}$ coupling constant in the well-
motivated BSM model using the charged lepton sector observables, which are the
muon $g-2$, the branching ratio of the CLFV $\tau\rightarrow\mu\gamma$ and
$\tau\rightarrow 3\mu$ decays, and then to find the highest order can reach up
to unity. Based on the determined order of $g_{X}$, $M_{Z^{\prime}}$ ranges
are derived from the neutrino trident production, $R_{K^{*}}$ anomaly, and
$B_{s}$ meson oscillation and we find there is no overlapped region which can
satisfy all constraints, thus considering only one obtained from $R_{K^{*}}$
anomaly and regarding one from $B_{s}$ meson oscillation to be explained by
some other new physics sources at higher energies. At this stage, we separated
$Z^{\prime}$ mass range derived from each fit of $R_{K^{*}}$ anomaly by two
cases where the first is one that the experimental CMS and upper-limit of
$\tau\rightarrow\mu\gamma$ are considered together, and the second is one
considering only the experimental upper-limit of $\tau\rightarrow\mu\gamma$
decay constraint and we call the second “theoretically interesting
$Z^{\prime}$ mass range. The reason why we distinguish two cases is it is
evident that the $R_{K^{*}}$ anomaly and muon $g-2$ can not be explained
simultaneously by the same new physics under all the SM constraints as there
is no overlapped $Z^{\prime}$ mass range (as reflected in the first case),
whereas the $Z^{\prime}$ for the $R_{K^{*}}$ anomaly and non-SM scalars for
the muon $g-2$ can be connected via the scalar potential under consideration
in this work while fitting the oblique parameters as well as the $W$ mass
anomaly in the “theoretically interesting $Z^{\prime}$ mass range”. Finally,
we explore the electroweak precision observables as well as muon $g-2$ anomaly
mediated by one loop level scalar exchange by fitting all the relevant scalar-
and gauge-mediated observables and confirm that they can be explained at their
$2\sigma$ constraints at most, thus concluding that the muon $g-2$ and
$R_{K^{*}}$ anomalies can not be explained by the same new physics, however
their new physics sources $Z^{\prime}$ and non-SM scalars (CP-even and -odd
scalars) are closely interconnected via the scalar potential as well as the
electroweak precision observables in the “theoretically interesting
$Z^{\prime}$ mass range”.

The rest of this paper is organized as follows. In section II, we introduce
the particle spectrum and the relevant Lagrangians. In section III, we
construct the mass matrices for the charged lepton and quark sectors and we
perform their diagonalization. In section IV, we determine the $Z^{\prime}$
couplings in the interaction basis and we derive the $Z^{\prime}$ interactions
in the mass basis. In section V, the lepton sector phenomenology is analyzed
by the massive neutral $Z^{\prime}$ gauge boson, determining the upper bound
of the $Z^{\prime}$ coupling constant $g_{X}$ considering the $Z^{\prime}$
mass as a free parameter. In section VI, the quark sector phenomenology is
explored with the obtained $Z^{\prime}$ coupling constant and the $Z^{\prime}$
mass ranges are derived from the quark sector observables like the $R_{K^{*}}$
anomaly and the $B_{s}$ meson oscillation parameters as well as from the
neutrino trident production. In section VIII, we analyze the electroweak
precision data and fit all the relevant scalar- and gauge-mediated observables
including muon $g-2$ and $W$ mass anomalies to their experimental bounds and
then finally we conclude that the muon $g-2$ and $R_{K^{*}}$ anomaly can not
be explained simultaneously by the same new physics, however they can be
successfully accommodated within their $2\sigma$ experimentally allowed range
in the “theoretically interesting $Z^{\prime}$ mass range”. Finally, we state
our conclusions in section IX. We relegate all the derived theoretical
predictions for the vacuum polarization amplitudes $\Pi$s contributing to the
oblique parameters $T,S,U$ to Appendices A to E.

##  II AN EXTENDED $Z^{\prime}$ MODEL WITH A FOURTH vectorlike FAMILY

The SM, which has survived against many experiments, is an important theory
which has successfully described the electromagnetic, weak and strong
interactions. However the SM is not a complete answer for our nature due to
the fact that it does not provide a mechanism for explaining the tiny neutrino
masses, the strong SM fermion mass hierarchy, the quark and lepton mixings as
well as the flavor anomalies and the $g-2$ anomalies. In order to address
these issues, we consider an extension of the SM motivated to provide a
mechanism that explain the SM fermion mass hierarchy as well as the other
issues previously mentioned. The well-measured SM mass parameters shows a very
strong hierarchy between the SM neutrino mass and the top quark mass.
Especially, the tiny SM neutrino masses are regarded to be explained by the
see-saw mechanism mediated by heavy right-handed neutrinos, rather than by the
Yukawa interaction. As in the SM neutrinos, it can be thought that the SM
charged fermions masses can also be generated by a seesaw mechanism mediated
by heavy charged vectorlike fermions. Therefore, all the SM fermion masses in
the BSM theory considered in this work can be explained by the same dynamical
mechanism. This approach has several advantages: 1) The effective Yukawa
couplings are proportional to a product of two other dimensionless couplings
at each vertex, so a moderate hierarchy in those couplings can yield a
quadratically larger hierarchy in the effective couplings, which together with
the mass of the heavy vectorlike fermions allow to successfully explain the
strong SM fermion mass hierarchy. 2) Given that SM charged fermion masses will
scale with the inverse of heavy vectorlike fermion masses, one can obtain a
range for the mass scale of heavy vectorlike fermions. 3) The vectorlike
charged leptons which mediate the seesaw mechanism that yields the SM charged
lepton masses are also crucial for radiatively generating the muon and
electron anomalous magnetic moments, whose magnitudes are not explained by
Standard Model. For the reasons mentioned, this economical BSM theory is quite
interesting and well-motivated. The particle content with their assigments
under the $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times U(1)^{\prime}$ gauge
symmetry of the BSM theory considered in this work are shown in Table 1:

Field | $Q_{iL}$ | $u_{iR}$ | $d_{iR}$ | $L_{iL}$ | $e_{iR}$ | $Q_{kL}$ | $u_{kR}$ | $d_{kR}$ | $L_{kL}$ | $e_{kR}$ | $\nu_{kR}$ | $\widetilde{Q}_{kR}$ | $\widetilde{u}_{kL}$ | $\widetilde{d}_{kL}$ | $\widetilde{L}_{kR}$ | $\widetilde{e}_{kL}$ | $\widetilde{\nu}_{kR}$ | $\phi$ | $H_{u}$ | $H_{d}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$SU(3)_{C}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$
$SU(2)_{L}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{2}$
$U(1)_{Y}$ | $\frac{1}{6}$ | $\frac{2}{3}$ | $-\frac{1}{3}$ | $-\frac{1}{2}$ | $1$ | $\frac{1}{6}$ | $\frac{2}{3}$ | $-\frac{1}{3}$ | $-\frac{1}{2}$ | $-1$ | $0$ | $\frac{1}{6}$ | $\frac{2}{3}$ | $-\frac{1}{3}$ | $-\frac{1}{2}$ | $-1$ | $0$ | $0$ | $\frac{1}{2}$ | $-\frac{1}{2}$
$U(1)^{\prime}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $1$ | $-1$ | $-1$ | $1$ | $-1$ | $-1$ | $1$ | $-1$ | $-1$ | $1$ | $-1$ | $-1$ | $1$ | $-1$ | $-1$

Table 1: An extended $Z^{\prime}$ model with one complete vectorlike family,
two SM-like Higgses, and one singlet flavon. This BSM theory is exactly the
same as the ones in Refs. Hernandez:2021tii ; CarcamoHernandez:2021yev ,
however in this work the SM gauge symmetry is extended by the local
$U(1)^{\prime}$ symmetry instead of global, therefore this BSM theory features
the massive neutral $Z^{\prime}$ gauge boson. This BSM theory is a combination
of an extended 2HDM model and a $Z^{\prime}$ model. After spontaneous symmetry
breaking (SSB) takes place, this BSM theory gives rise to the effective SM
Yukawa interactions at the electroweak energy scale.

With the particle content shown in table 1, the renormalizable Yukawa
interactions of this model at high energy scale yield the following effective
SM Yukawa interactions at the electroweak energy scale:

$\mathcal{L}_{\mathop{\mathrm{e}ff}}^{\mathop{\mathrm{Y}ukawa}}=y_{ik}^{\psi}(M_{\psi^{\prime}}^{-1})_{kl}{x_{lj}^{\psi^{\prime}}\left\langle\phi\right\rangle}\overline{\psi}_{iL}\widetilde{H}\psi_{jR}+{x_{ik}^{\psi}\left\langle\phi\right\rangle}(M_{\psi}^{-1})_{kl}y_{lj}^{\psi}\overline{\psi}_{iL}\widetilde{H}\psi_{jR}+\mathop{\mathrm{h}.c.}$
(6)

where the indices $i,j=1,2,3$ and $k,l=4$ and $\psi,\psi^{\prime}=Q,u,d,L,e$
and $M$ means heavy vectorlike mass. The corresponding Feynman diagrams are
shown in Figure 1.

Figure 1: Diagrams in this model which lead to the effective Yukawa
interactions, where $\psi,\psi^{\prime}=Q,u,d,L,e$(neutrinos will be treated
separately) $i,j=1,2,3$, $k,l=4$, $M_{lk}$ is vectorlike mass and
$\widetilde{H}=i\sigma_{2}H^{*},H=H_{u,d}$

Notice that the left and right diagrams of Figure 1 are mediated by vectorlike
singlet and vectorlike doublet fermions, respectively. Therefore, a complete
vectorlike family can give rise to two fermionic seesaw operators and this
feature leads to massless first generation of SM fermions, whereas the second
and third families of SM fermions do acquire their masses. Adding one more
vectorlike family will imply the inclusion of two extra fermionic seesaw
operators thus allowing all SM fermions to acquire their masses. However we
only consider one complete vectorlike family for the purpose of keeping this
BSM theory economical and the phenomenological analysis as simple as possible.

###  II.1 Effective Yukawa interactions for the SM fermions

The renormalizable Lagrangian for the quark sector is given by:

$\begin{split}\mathcal{L}_{q}^{\mathop{\mathrm{Y}ukawa+Mass}}&=y_{ik}^{u}\overline{Q}_{iL}\widetilde{H}_{u}u_{kR}+x_{ki}^{u}\phi\overline{\widetilde{u}}_{kL}u_{iR}+x_{ik}^{Q}\phi\overline{Q}_{iL}\widetilde{Q}_{kR}+y_{ki}^{u}\overline{Q}_{kL}\widetilde{H}_{u}u_{iR}\\\
&+y_{ik}^{d}\overline{Q}_{iL}\widetilde{H}_{d}d_{kR}+x_{ki}^{d}\phi\overline{\widetilde{d}}_{kL}d_{iR}+y_{ki}^{d}\overline{Q}_{kL}\widetilde{H}_{d}d_{iR}\\\
&+M_{kl}^{u}\overline{\widetilde{u}}_{lL}u_{kR}+M_{kl}^{d}\overline{\widetilde{d}}_{lL}d_{kR}+M_{kl}^{Q}\overline{Q}_{kL}\widetilde{Q}_{lR}+\mathop{\mathrm{h}.c.}\end{split}$
(7)

where $i,j=1,2,3$ and $k,l=4,5$ and $\widetilde{H}=i\sigma^{2}H^{*}$. After
the singlet flavon $\phi$ develops its vacuum expectation value (vev) and the
vectorlike fermions are integrated out, the renormalizable quark Yukawa terms
give rise to the effective SM quark Yukawa interaction. The corresponding
effective SM Yukawa interactions are given in Figure 2:

Figure 2: Diagrams in this model which lead to the effective Yukawa
interactions for the up quark sector(two above diagrams) and the down quark
sector(two below diagrams) in mass insertion formalism, where $i,j=1,2,3$ and
$k,l=4$ and $M_{lk}$ is vectorlike mass.

As in the quark sector, the renormalizable Yukawa terms for the charged lepton
sector can be written as follows:

$\begin{split}\mathcal{L}_{e}^{\mathop{\mathrm{Y}ukawa+Mass}}&=y_{ik}^{e}\overline{L}_{iL}\widetilde{H}_{d}e_{kR}+x_{ki}^{e}\phi\overline{\widetilde{e}}_{kL}e_{iR}+x_{ik}^{L}\phi\overline{L}_{iL}\widetilde{L}_{kR}+y_{ki}^{e}\overline{L}_{kL}\widetilde{H}_{d}e_{iR}\\\
&+M_{kl}^{e}\overline{\widetilde{e}}_{lL}e_{kR}+M_{kl}^{L}\overline{L}_{kL}\widetilde{L}_{lR}+\mathop{\mathrm{h}.c.},\end{split}$
(8)

which give rise to the effective SM charged lepton Yukawa interactions in
Figure 3:

Figure 3: Diagrams in this model which lead to the effective Yukawa
interactions for the charged lepton sector in mass insertion formalism, where
$i,j=1,2,3$ and $k,l=4$ and $M_{lk}$ is vectorlike mass

As for the SM neutrino sector, it requires a different approach compared to
the SM charged fermions since we assume that the SM left-handed neutrinos are
Majorana particles and they have their masses via expansion through the
vectorlike neutrinos. The renormalizable neutrino Yukawa terms are given by:

$\begin{split}\mathcal{L}_{\nu}^{\mathop{\mathrm{Y}ukawa+Mass}}=y_{ik}^{\nu}\overline{L}_{iL}\widetilde{H}_{u}\nu_{kR}+x_{ik}^{L}\overline{L}_{iL}H_{d}\overline{\widetilde{\nu}}_{kR}+M_{kl}^{M}\overline{\widetilde{\nu}}_{lR}\nu_{kR}+\mathop{\mathrm{h}.c.}\end{split}$
(9)

After the heavy vectorlike neutrinos are integrated out, the renormalizable
Lagrangian induces the Weinberg-like operator which generates the “type 1b
seesaw mechanism” in Figure 4 to differentiate it from the original Weinberg
operator corresponding to the type 1a seesaw mechanism which does not work in
the model under consideration since the SM-like Higgses are charged under the
local $U(1)^{\prime}$ symmetry.

Figure 4: Type Ib seesaw mechanism Hernandez:2021tii ; Hernandez-
Garcia:2019uof  which leads to the effective Yukawa interactions for the
Majorana neutrinos in mass insertion formalism, where $i,j=1,2,3$ and
$k,l=4,5$ and $M_{lk}$ is vectorlike mass.

There are three important features we can read off from the vectorlike mass in
Figure 4. The first one is the nature of the vectorlike mass consisting of two
different vectorlike neutrinos, which is differentiated by the Majorana mass
consisting of the same neutrinos. The second one is the vectorlike mass in
Figure 4 which also violates the lepton number conservation as the Majorana
mass does and this feature can be easily confirmed by checking the lepton
number in the renormalizable neutrino Yukawa interaction of Equation 9. The
third one is that the tiny active neutrino masses are explained from the type
1b seesaw mechanism where each interaction vertex involves two SM Higgs vevs
$v_{u},v_{d}$
($v_{u}^{2}+v_{d}^{2}=\left(246\mathop{\mathrm{G}eV}\right)^{2}$) and allowing
two different vertex values, therefore it is possible to significantly lower
the order of predicted right-handed neutrinos from
$10^{14}\mathop{\mathrm{G}eV}$ up to $\mathop{\mathrm{T}eV}$ scale.

We call the mass insertion formalism introduced from Figure 2 to 4 “low energy
scale seesaw mechanism” and this mechanism can explain both the SM fermion
masses and their mixings as well. On top of that, we also showed that the SM
$Z$ gague boson can cause the FCNCs in both quark and lepton sectors with one
complete vectorlike family and how the unitarity of the first row of the CKM
matrix can be relaxed in an analytical way in one of our works
CarcamoHernandez:2021yev . Based on what we have discussed, we construct the
mass matrices in both quark and charged lepton sectors and then diagonalize
them in the next section.

##  III EFFECTIVE YUKAWA MATRICES USING A MIXING FORMALISM

The mass matrices for the SM fermions and their diagonalization are addressed
in this section and they are well described in section III of Ref.
CarcamoHernandez:2021yev , so we focus on important features of the mass
matrices instead of explaining every details about them. The general mass
matrix for fermions in the flavor basis is given by:

$M^{\psi}=\left(\begin{array}[]{c|ccccc}&\psi_{1R}&\psi_{2R}&\psi_{3R}&\psi_{4R}&\widetilde{\psi}_{4R}\\\\[2.15277pt]
\hline\cr\overline{\psi}_{1L}&0&0&0&y_{14}^{\psi}\langle\widetilde{H}^{0}\rangle&x_{14}^{\psi}\langle\phi\rangle\\\\[4.30554pt]
\overline{\psi}_{2L}&0&0&0&y_{24}^{\psi}\langle\widetilde{H}^{0}\rangle&x_{24}^{\psi}\langle\phi\rangle\\\\[4.30554pt]
\overline{\psi}_{3L}&0&0&0&y_{34}^{\psi}\langle\widetilde{H}^{0}\rangle&x_{34}^{\psi}\langle\phi\rangle\\\\[4.30554pt]
\overline{\psi}_{4L}&y_{41}^{\psi}\langle\widetilde{H}^{0}\rangle&y_{42}^{\psi}\langle\widetilde{H}^{0}\rangle&y_{43}^{\psi}\langle\widetilde{H}^{0}\rangle&0&M_{44}^{\psi}\\\\[4.30554pt]
\overline{\widetilde{\psi}}_{4L}&x_{41}^{\psi^{\prime}}\langle\phi\rangle&x_{42}^{\psi^{\prime}}\langle\phi\rangle&x_{43}^{\psi^{\prime}}\langle\phi\rangle&M_{44}^{\psi^{\prime}}&0\\\
\end{array}\right),$ (10)

where $\psi,\psi^{\prime}=Q,u,d,L,e$ and the zeros appearing in the upper-left
$3\times 3$ block means the SM fermions remain massless without extension of
the SM, which is implemented by assigning $U(1)^{\prime}$ charges to the SM-
like Higgses $H_{u,d}$, and lastly the mass matrix involves three different
mass scales $\langle H^{0}\rangle,\langle\phi\rangle$ and $M$ so it can
naturally accommodate the strong SM hierarchy in this BSM model. In order to
dynamically implement the SM hierarchy, it requires to rotate the mass matrix
of Equation 10 maximally and this feature will be discussed across the next
subsections.

###  III.1 Diagonalizing the charged lepton mass matrix

After all scalars develop its vevs ($v_{d}=\langle H_{d}^{0}\rangle$ and
$v_{\phi}=\langle\phi\rangle$) and then the mass matrix for the charged lepton
is fully rotated, it can be written as follows:

$M^{e}=\left(\begin{array}[]{c|ccccc}&e_{1R}&e_{2R}&e_{3R}&e_{4R}&\widetilde{L}_{4R}\\\\[2.15277pt]
\hline\cr\overline{L}_{1L}&0&0&0&0&0\\\\[4.30554pt]
\overline{L}_{2L}&0&0&0&y_{24}^{e}v_{d}&0\\\\[4.30554pt]
\overline{L}_{3L}&0&0&0&y_{34}^{e}v_{d}&x_{34}^{L}v_{\phi}\\\\[4.30554pt]
\overline{L}_{4L}&0&0&y_{43}^{e}v_{d}&0&M_{44}^{L}\\\\[4.30554pt]
\overline{\widetilde{e}}_{4L}&0&x_{42}^{e}v_{\phi}&x_{43}^{e}v_{\phi}&M_{44}^{e}&0\\\
\end{array}\right),$ (11)

where it is worth mentioning that the fully rotated mass matrix can provide a
dynamical explanation for the SM hierarchy and that with one complete
vectorlike family can give rise to two fermionic seesaw mediators, which means
that two of the SM generations of fermions, i.e., the second and the third one
can acquire masses. Adding one extra vectorlike family can allow for the whole
generations of SM fermions to be massive, as done in one of our works
Hernandez:2021tii , however this approach requires some assumptions to
simplify the resulting mass matrices since they include many fields. For the
purpose of avoiding any assumptions while diagonalizing the mass matrices for
the SM fermions and also keeping those as economical as possible, we constrain
our attention to the massive SM second and third generations with one complete
vectorlike family. It is convenient to rewrite the Yukawa terms by simple mass
parameters and then switch the column of $e_{4R}$ and $\widetilde{L}_{4R}$ in
order for the heavy vectorlike masses to locate at the diagonal position as
follows:

$M^{e}=\left(\begin{array}[]{c|ccccc}&e_{1R}&e_{2R}&e_{3R}&e_{4R}&\widetilde{L}_{4R}\\\\[2.15277pt]
\hline\cr\overline{L}_{1L}&0&0&0&0&0\\\\[4.30554pt]
\overline{L}_{2L}&0&0&0&m_{24}&0\\\\[4.30554pt]
\overline{L}_{3L}&0&0&0&m_{34}&m_{35}\\\\[4.30554pt]
\overline{L}_{4L}&0&0&m_{43}&0&M_{45}^{L}\\\\[4.30554pt]
\overline{\widetilde{e}}_{4L}&0&m_{52}&m_{53}&M_{54}^{e}&0\\\
\end{array}\right)=\left(\begin{array}[]{c|ccccc}&e_{1R}&e_{2R}&e_{3R}&\widetilde{L}_{4R}&e_{4R}\\\\[2.15277pt]
\hline\cr\overline{L}_{1L}&0&0&0&0&0\\\\[4.30554pt]
\overline{L}_{2L}&0&0&0&0&m_{24}\\\\[4.30554pt]
\overline{L}_{3L}&0&0&0&m_{35}&m_{34}\\\\[4.30554pt]
\overline{L}_{4L}&0&0&m_{43}&M_{45}^{L}&0\\\\[4.30554pt]
\overline{\widetilde{e}}_{4L}&0&m_{52}&m_{53}&0&M_{54}^{e}\\\
\end{array}\right),$ (12)

where the left-handed fields follow the order of $12345$ whereas the right-
handed fields follow that of $12354$ (the tilde particles have the index $5$).
The mass matrix of Equation 12 is diagonalized by two methods which are the
numerical singular value decomposition (SVD) and the step-by-step analytical
diagonalization and it is confirmed that subtracting from a result from one
method to that from the other method is quite ignorable
CarcamoHernandez:2021yev . In this work, we mainly make use of the numerical
singular value decomposition method to focus further on the relevant
phenomenology. The rearranged mass matrix for the charged lepton sector is
diagonalized by:

$M_{\mathop{\mathrm{d}iag}}^{e}=\mathop{\mathrm{d}iag}\left(0,m_{\mu},m_{\tau},M_{E_{4}},M_{\widetilde{E}_{4}}\right)=V^{L}M^{e}(V^{e})^{\dagger},$
(13)

where $V^{L}(V^{e})$ is the left-handed(right-handed) mixing matrix for the
charged lepton fields and they connect the flavor basis and the physical basis
as follows:

$\begin{split}\begin{pmatrix}e_{L}\\\\[2.15277pt] \mu_{L}\\\\[2.15277pt]
\tau_{L}\\\\[2.15277pt] E_{4L}\\\\[2.15277pt]
\widetilde{E}_{4L}\end{pmatrix}=V^{L}\begin{pmatrix}e_{1L}\\\\[2.15277pt]
e_{2L}\\\\[2.15277pt] e_{3L}\\\\[2.15277pt] e_{4L}\\\\[2.15277pt]
\widetilde{e}_{4L}\end{pmatrix},\qquad\begin{pmatrix}e_{R}\\\\[2.15277pt]
\mu_{R}\\\\[2.15277pt] \tau_{R}\\\\[2.15277pt]
\widetilde{E}_{4R}\\\\[2.15277pt]
E_{4R}\end{pmatrix}=V^{e}\begin{pmatrix}e_{1R}\\\\[2.15277pt]
e_{2R}\\\\[2.15277pt] e_{3R}\\\\[2.15277pt] \widetilde{e}_{4R}\\\\[2.15277pt]
e_{4R}\end{pmatrix},\end{split}$ (14)

###  III.2 Diagonalizing the up-type quark mass matrix

The mass matrix for the up-type quark sector can be approached as in the
charged lepton sector. The rearranged form is given by:

$M^{u}=\left(\begin{array}[]{c|ccccc}&u_{1R}&u_{2R}&u_{3R}&u_{4R}&\widetilde{Q}_{4R}\\\\[2.15277pt]
\hline\cr\overline{Q}_{1L}&0&0&0&0&0\\\\[4.30554pt]
\overline{Q}_{2L}&0&0&0&y_{24}^{u}v_{u}&0\\\\[4.30554pt]
\overline{Q}_{3L}&0&0&0&y_{34}^{u}v_{u}&x_{34}^{Q}v_{\phi}\\\\[4.30554pt]
\overline{Q}_{4L}&0&0&y_{43}^{u}v_{u}&0&M_{44}^{Q}\\\\[4.30554pt]
\overline{\widetilde{u}}_{4L}&0&x_{42}^{u}v_{\phi}&x_{43}^{u}v_{\phi}&M_{44}^{u}&0\\\
\end{array}\right)=\left(\begin{array}[]{c|ccccc}&u_{1R}&u_{2R}&u_{3R}&\widetilde{Q}_{4R}&u_{4R}\\\\[2.15277pt]
\hline\cr\overline{Q}_{1L}&0&0&0&0&0\\\\[4.30554pt]
\overline{Q}_{2L}&0&0&0&0&m_{24}^{u}\\\\[4.30554pt]
\overline{Q}_{3L}&0&0&0&m_{35}^{u}&m_{34}^{u}\\\\[4.30554pt]
\overline{Q}_{4L}&0&0&m_{43}^{u}&M_{44}^{Q}&0\\\\[4.30554pt]
\overline{\widetilde{u}}_{4L}&0&m_{52}^{u}&m_{53}^{u}&0&M_{44}^{u}\\\
\end{array}\right),$ (15)

Here, the whole form is exactly consistent with the one for the charged lepton
sector except for a few substitutions such as $y^{e}\rightarrow
y^{u},v_{d}\rightarrow v_{u},x^{L}\rightarrow x^{Q}$, and $x^{e}\rightarrow
x^{u}$. Therefore, we can simply reuse the derived result for the charged
lepton sector. The up-type quark mass matrix is diagonalized by:

$M_{\mathop{\mathrm{d}iag}}^{u}=\mathop{\mathrm{d}iag}\left(0,m_{c},m_{t},M_{U_{4}},M_{\widetilde{U}_{4}}\right)=V_{L}^{u}M^{u}(V_{R}^{u})^{\dagger},$
(16)

where the mixing matrices $V_{L,R}^{u}$ are defined as follows:

$\begin{split}\begin{pmatrix}u_{L}\\\\[2.15277pt] c_{L}\\\\[2.15277pt]
t_{L}\\\\[2.15277pt] U_{4L}\\\\[2.15277pt]
\widetilde{U}_{4L}\end{pmatrix}=V_{L}^{u}\begin{pmatrix}u_{1L}\\\\[2.15277pt]
u_{2L}\\\\[2.15277pt] u_{3L}\\\\[2.15277pt] u_{4L}\\\\[2.15277pt]
\widetilde{u}_{4L}\end{pmatrix},\qquad\begin{pmatrix}u_{R}\\\\[2.15277pt]
c_{R}\\\\[2.15277pt] t_{R}\\\\[2.15277pt] \widetilde{U}_{4R}\\\\[2.15277pt]
U_{4R}\end{pmatrix}=V_{R}^{u}\begin{pmatrix}u_{1R}\\\\[2.15277pt]
u_{2R}\\\\[2.15277pt] u_{3R}\\\\[2.15277pt] \widetilde{u}_{4R}\\\\[2.15277pt]
u_{4R}\end{pmatrix}.\end{split}$ (17)

###  III.3 Diagonalizing the down-type quark mass matrix

The down-type quark mass matrix requires some attention as its form is
somewhat different, compared to that for the charged lepton or up-type quark
sector. For comparison, it is convenient to look at the up- and down-type
quark mass matrices together

$M^{u}=\left(\begin{array}[]{c|ccccc}&u_{1R}&u_{2R}&u_{3R}&u_{4R}&\widetilde{Q}_{4R}\\\\[2.15277pt]
\hline\cr\overline{Q}_{1L}&0&0&0&0&0\\\\[4.30554pt]
\overline{Q}_{2L}&0&0&0&y_{24}^{u}v_{u}&0\\\\[4.30554pt]
\overline{Q}_{3L}&0&0&0&y_{34}^{u}v_{u}&x_{34}^{Q}v_{\phi}\\\\[4.30554pt]
\overline{Q}_{4L}&0&0&y_{43}^{u}v_{u}&0&M_{44}^{Q}\\\\[4.30554pt]
\overline{\widetilde{u}}_{4L}&0&x_{42}^{u}v_{\phi}&x_{43}^{u}v_{\phi}&M_{44}^{u}&0\\\
\end{array}\right),\quad
M^{d}=\left(\begin{array}[]{c|ccccc}&d_{1R}&d_{2R}&d_{3R}&d_{4R}&\widetilde{Q}_{4R}\\\\[2.15277pt]
\hline\cr\overline{Q}_{1L}&0&0&0&y_{14}^{d}v_{d}&0\\\\[4.30554pt]
\overline{Q}_{2L}&0&0&0&y_{24}^{d}v_{d}&0\\\\[4.30554pt]
\overline{Q}_{3L}&0&0&0&y_{34}^{d}v_{d}&x_{34}^{Q}v_{\phi}\\\\[4.30554pt]
\overline{Q}_{4L}&0&0&y_{43}^{d}v_{d}&0&M_{44}^{Q}\\\\[4.30554pt]
\overline{\widetilde{d}}_{4L}&0&x_{42}^{d}v_{\phi}&x_{43}^{d}v_{\phi}&M_{44}^{d}&0\\\
\end{array}\right).$ (18)

The first important feature is the fifth column of the up-type mass matrix is
exactly consistent with that of the down-type mass matrix King:2018fcg . What
this implements is the left-handed vectorlike quark doublets contribute to the
right-handed vectorlike doublets equally, so the vectorlike quark doublets get
to have very degenerate masses, and this feature plays a crucial role in
constraining violation of the custodial symmetry as we will see in the oblique
parameter section. The second feature is the presence of the Yukawa term
$y_{14}^{d}v_{d}$. For the up-quark sector, we can rotate $Q_{1L}$ and
$Q_{2L}$ fields further to vanish the up-type Yukawa term $y_{14}^{u}v_{u}$,
however this rotation simply remixes the down-type Yukawa terms
$y_{14}^{d}v_{d}$ and $y_{24}^{d}v_{d}$ and therefore both the down-type
Yukawa terms survive. This rotation between $Q_{1L}$ and $Q_{2L}$ fields can
be done first in the down-type mass matrix, remaining the unvanished up-type
Yukawa term, however this different approach does not change our analyses
since we found that the analyses which will be explored are basis-independent
as confirmed in the SM $Z$ physics in our previous work
CarcamoHernandez:2021yev . Therefore, we consider it is phenomenologically
reasonable to consider the unvanished Yukawa term in the down-type mass matrix
in order to explain the CKM mixing matrix (or equally, the CKM mixing matrix
is mainly based on the down-type quark mixings and this is one of our findings
in the previous work CarcamoHernandez:2021yev ). The last feature to discuss
is the down-type left-handed quark fields can access to all the SM mixings
even thought the first generation of the down-type quarks ($d$ quark) remains
massless in the model under consideration. It can be easily confirmed by
looking at a partially diagonalized mass matrix for the up- and down-type
quark sector. After integrating out the heavy degrees of freedom in both the
mass matrices, the mass matrices have the forms as follows:

$M^{u\prime}=\left(\begin{array}[]{c|ccccc}&u_{1R}&u_{2R}^{\prime}&u_{3R}^{\prime}&\widetilde{Q}_{4R}^{\prime}&u_{4R}^{\prime}\\\\[2.15277pt]
\hline\cr\overline{Q}_{1L}&0&0&0&0&0\\\\[4.30554pt]
\overline{Q}_{2L}^{\prime}&0&m_{22}^{u\prime}&m_{23}^{u\prime}&0&0\\\\[4.30554pt]
\overline{Q}_{3L}^{\prime}&0&m_{32}^{u\prime}&m_{33}^{u\prime}&0&0\\\\[4.30554pt]
\overline{Q}_{4L}^{\prime}&0&0&0&M_{44}^{Q}&0\\\\[4.30554pt]
\overline{\widetilde{u}}_{4L}^{\prime}&0&0&0&0&M_{44}^{u}\\\
\end{array}\right),\quad
M^{d\prime}=\left(\begin{array}[]{c|ccccc}&d_{1R}&d_{2R}^{\prime}&d_{3R}^{\prime}&\widetilde{Q}_{4R}^{\prime}&d_{4R}^{\prime}\\\\[2.15277pt]
\hline\cr\overline{Q}_{1L}^{\prime}&0&m_{12}^{d\prime}&m_{13}^{d\prime}&0&0\\\\[4.30554pt]
\overline{Q}_{2L}^{\prime}&0&m_{22}^{d\prime}&m_{23}^{d\prime}&0&0\\\\[4.30554pt]
\overline{Q}_{3L}^{\prime}&0&m_{32}^{d\prime}&m_{33}^{d\prime}&0&0\\\\[4.30554pt]
\overline{Q}_{4L}^{\prime}&0&0&0&M_{44}^{Q}&0\\\\[4.30554pt]
\overline{\widetilde{d}}_{4L}^{\prime}&0&0&0&0&M_{44}^{d}\\\
\end{array}\right),$ (19)

where the primed field means the field is rotated. The partially diagonalized
down-type mass matrix has two more Yukawa terms $m_{12,13}^{d\prime}$ compared
to the up-type mass matrix and this term allows left-handed mixing with the
first generation of the down-type quarks even though the first generation
remains massless (note that right-handed mixing in the down-type mass matrix
is only 23 mixing as in the other mass matrices). Then the down-type mass
matrix is diagonalized by:

$M_{\mathop{\mathrm{d}iag}}^{d}=\mathop{\mathrm{d}iag}\left(0,m_{s},m_{b},M_{D_{4}},M_{\widetilde{D}_{4}}\right)=V_{L}^{d}M^{d}(V_{R}^{d})^{\dagger},$
(20)

where the mixing matrices $V_{L,R}^{d}$ are defined by:

$\begin{split}\begin{pmatrix}d_{L}\\\ s_{L}\\\ b_{L}\\\ D_{4L}\\\
\widetilde{D}_{4L}\end{pmatrix}=V_{L}^{d}\begin{pmatrix}d_{1L}\\\ d_{2L}\\\
d_{3L}\\\ d_{4L}\\\
\widetilde{d}_{4L}\end{pmatrix},\qquad\begin{pmatrix}d_{R}\\\ s_{R}\\\
b_{R}\\\ D_{4R}\\\
\widetilde{D}_{4R}\end{pmatrix}=V_{R}^{d}\begin{pmatrix}d_{1R}\\\ d_{2R}\\\
d_{3R}\\\ \widetilde{d}_{4R}\\\ d_{4R}\end{pmatrix}.\end{split}$ (21)

##  IV THE $Z^{\prime}$ GAUGE BOSON INTERACTIONS WITH THE VECTORLIKE FAMILY

As the model under consideration features the extended gauge symmetry by the
$U(1)^{\prime}$ local symmetry, this BSM model involves a neutral massive
$Z^{\prime}$ gauge boson. Considering a new physics source $Z^{\prime}$ for
the anomalies is interesting, since it can make a direct connection to the SM
observables and anomalies.

###  IV.1 FCNC mediated by the $Z^{\prime}$ gauge boson with the fourth
vectorlike fermions

The $Z^{\prime}$ interactions with the fourth vectorlike fermions in the
flavor basis are given by:

$\mathcal{L}_{\mathop{\mathrm{f}lavor}}^{Z^{\prime}}=g_{X}Z_{\mu}^{\prime}\left(\overline{Q}_{iL}D_{Q}\gamma^{\mu}Q_{iL}+\overline{u}_{iR}D_{u}\gamma^{\mu}u_{iR}+\overline{d}_{iR}D_{d}\gamma^{\mu}d_{iR}+\overline{L}_{iL}D_{L}\gamma^{\mu}L_{iL}+\overline{e}_{iR}D_{e}\gamma^{\mu}e_{iR}+\overline{\nu}_{iR}D_{\nu}\gamma^{\mu}\nu_{iR}\right)$
(22)

where $i=1,2,3,4,5$ and the $U(1)^{\prime}$ charge matrices, $D$s, are defined
in the model under consideration as follows:

$\begin{split}D_{Q}=\mathop{\mathrm{d}iag}\left(0,0,0,1,-1\right),\quad
D_{u}&=\mathop{\mathrm{d}iag}\left(0,0,0,1,-1\right),\quad\quad
D_{d}=\mathop{\mathrm{d}iag}\left(0,0,0,1,-1\right),\\\
D_{L}=\mathop{\mathrm{d}iag}\left(0,0,0,1,-1\right),\quad
D_{e}&=\mathop{\mathrm{d}iag}\left(0,0,0,1,-1\right),\quad\quad
D_{\nu}=\mathop{\mathrm{d}iag}\left(0,0,0,-1,-1\right),\end{split}$ (23)

where it is worth mentioning that the left-handed fields $Q_{iL},L_{iL}$
follow the order of $12345$ and the right-handed fields
$u_{iR},d_{iR},e_{iR},\nu_{iR}$ follow the order of $12354$ as defined in the
mass matrix diagonalization of section III and the orders determine each $D$
of Equation 23. This BSM model follows the fermionic $Z^{\prime}$ model
CarcamoHernandez:2019ydc , which means that the SM fermions do not interact
with the $Z^{\prime}$ gauge boson in the flavor basis however it starts to
have interactions with the SM fermions via mixings in the physical basis. The
$Z^{\prime}$ interactions in the physical basis are given by:

$\mathcal{L}_{\mathop{\mathrm{p}hysical}}^{Z^{\prime}}=g_{X}Z_{\mu}^{\prime}\left(\overline{Q}_{L}D_{Q}^{\prime}\gamma^{\mu}Q_{L}+\overline{u}_{R}D_{uR}^{\prime}\gamma^{\mu}u_{R}+\overline{d}_{R}D_{dR}^{\prime}\gamma^{\mu}d_{R}+\overline{L}_{L}D_{L}^{\prime}\gamma^{\mu}L_{L}+\overline{e}_{R}D_{eR}^{\prime}\gamma^{\mu}e_{R}+\overline{\nu}_{R}D_{\nu}^{\prime}\gamma^{\mu}\nu_{R}\right),$
(24)

where the primed $D$s are determined as follows ($D_{Q}^{\prime}$ is separated
to $D_{uL}^{\prime}$ and $D_{dL}^{\prime}$ after SSB):

$\begin{split}D_{uL}^{\prime}&=V_{L}^{u}D_{Q}V_{L}^{u\dagger},\quad
D_{uR}^{\prime}=V_{R}^{u}D_{u}V_{R}^{u\dagger},\\\
D_{dL}^{\prime}&=V_{L}^{d}D_{Q}V_{L}^{d\dagger},\quad
D_{dR}^{\prime}=V_{R}^{d}D_{d}V_{R}^{d\dagger},\\\
D_{eL}^{\prime}&=V_{L}^{e}D_{L}V_{L}^{e\dagger},\quad
D_{eR}^{\prime}=V_{R}^{e}D_{e}V_{R}^{e\dagger},\\\
D_{\nu}^{\prime}&=V^{\nu}D_{\nu}V^{\nu\dagger},\\\ \end{split}$ (25)

The derived $D^{\prime}$s of Equation LABEL:eqn:Dprimes induce the flavor
changing neutral currents (FCNCs) so this feature allows to access tree level
contributions to the $\tau\rightarrow\mu\gamma$ decay or $R_{K}$ anomaly, etc.
and those will be discussed in detail in the next sections.

###  IV.2 FCNC mediated by the $Z^{\prime}$ gauge boson in the neutrino sector
with the fourth vectorlike neutrinos

The neutrino sector in the model under consideration requires another approach
compared to the charged fermion sectors. We build the mass matrices for the
charged fermion sectors explicitly since our analyses on the observables which
will be explored are sensitive to mixing parameters coming from their mixing
matrices. As for the neutrino sector, what is required is the $Z^{\prime}$
coupling to $\nu_{\mu}$ pair for the neutrino trident constraint. Therefore,
we do not need to have a full mass matrix for the neutrino sector instead we
consider only the coupling constant of $Z^{\prime}\nu_{\mu}\nu_{\mu}$
interaction. As covered in section III, the neutrino sector in this BSM model
is explained by the type $1$b seesaw mechanism and the analysis was done in
one of our works Hernandez:2021tii . The $D_{\nu}^{\prime}$ is given by
($\alpha=\begin{pmatrix}-1&0\\\ 0&-1\end{pmatrix}$):

$\begin{split}D_{\nu}^{\prime}&=VD_{\nu}V^{T}\\\ &=V\begin{pmatrix}0_{3\times
3}&0_{3\times 2}\\\ 0_{2\times 3}&\alpha\end{pmatrix}V^{T}\\\
&=\begin{pmatrix}U_{\mathop{\mathrm{P}MNS}}&0\\\
0&I\end{pmatrix}\begin{pmatrix}I-\frac{\Theta\Theta^{\dagger}}{2}&\Theta\\\
-\Theta^{\dagger}&I-\frac{\Theta\Theta^{\dagger}}{2}\end{pmatrix}\begin{pmatrix}0_{3\times
3}&0_{3\times 2}\\\ 0_{2\times
3}&\alpha\end{pmatrix}\begin{pmatrix}I-\frac{\Theta\Theta^{\dagger}}{2}&\Theta\\\
-\Theta^{\dagger}&I-\frac{\Theta\Theta^{\dagger}}{2}\end{pmatrix}^{T}\begin{pmatrix}U_{\mathop{\mathrm{P}MNS}}&0\\\
0&I\end{pmatrix}^{T}\\\ &=\begin{pmatrix}U_{\mathop{\mathrm{P}MNS}}&0\\\
0&I\end{pmatrix}\begin{pmatrix}I-\eta_{ij}&\Theta_{ik}\\\
-\Theta_{ki}^{\dagger}&I-\eta_{kl}\end{pmatrix}\begin{pmatrix}0_{3\times
3}&0_{3\times 2}\\\ 0_{2\times
3}&\alpha\end{pmatrix}\begin{pmatrix}I-\eta_{ij}&\Theta_{ik}\\\
-\Theta_{ki}^{\dagger}&I-\eta_{kl}\end{pmatrix}^{T}\begin{pmatrix}U_{\mathop{\mathrm{P}MNS}}&0\\\
0&I\end{pmatrix}^{T}\\\ &=\begin{pmatrix}U_{\mathop{\mathrm{P}MNS}}&0\\\
0&I\end{pmatrix}\begin{pmatrix}\Theta_{ik}\alpha\Theta_{ki}^{T}&\Theta_{ik}\alpha\left(I-\eta_{lk}^{T}\right)\\\
\left(I-\eta_{kl}\right)\alpha\Theta_{ki}^{T}&\left(I-\eta_{kl}\right)\alpha\left(I-\eta_{lk}^{T}\right)\end{pmatrix}\begin{pmatrix}U_{\mathop{\mathrm{P}MNS}}&0\\\
0&I\end{pmatrix}^{T},\end{split}$ (26)

where the indices $i,j=1,2,3$ and $k,l=4,5$ and $U_{\mathop{\mathrm{P}MNS}}$
is the well-known Pontecorvo-Maki-Nakagawa-Sakata matrix and is parametrized
as follows ($c_{ij}=\cos\theta_{ij},s_{ij}=\sin\theta_{ij}$ and
$\delta_{\mathop{\mathrm{C}P}}$ : the Dirac phase, $\alpha^{\prime},\alpha$ :
the Majorana phases) Hernandez-Garcia:2019uof ; Chau:1984fp :

$U_{\mathop{\mathrm{P}MNS}}=\begin{pmatrix}1&0&0\\\ 0&c_{23}&s_{23}\\\
0&-s_{23}&c_{23}\end{pmatrix}\begin{pmatrix}c_{13}&0&s_{13}e^{-i\delta_{\mathop{\mathrm{C}P}}}\\\
0&1&0\\\
-s_{13}e^{i\delta_{\mathop{\mathrm{C}P}}}&0&c_{13}\end{pmatrix}\begin{pmatrix}c_{12}&s_{12}&0\\\
-s_{12}&c_{12}&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}e^{-i\alpha^{\prime}/2}&0&0\\\
0&e^{-i\alpha/2}_{12}&0\\\ 0&0&1\end{pmatrix},$ (27)

and $\Theta$ and the deviation of unitarity $\eta$ coming from the dimension
six operators are defined as follows Blennow:2011vn ; Fernandez-
Martinez:2007iaa ; Broncano:2002rw :

$\Theta\simeq
m_{D}^{\dagger}M_{N}^{-1},\quad\eta=\frac{\Theta\Theta^{\dagger}}{2}=\frac{1}{2}\frac{m_{D}^{\dagger}m_{D}}{M_{N}^{2}}.$
(28)

We do not need to care about all the elements of $D_{\nu}^{\prime}$ instead
focus on the upper-left $3\times 3$ block. Going back to the derivation, it is
given by:

$\begin{split}D_{\nu}^{\prime}&=U_{\mathop{\mathrm{P}MNS}}\Theta_{ik}\alpha\Theta_{ki}^{T}U_{\mathop{\mathrm{P}MNS}}^{T},\quad\text{for
the upper-left $3\times 3$ block}\\\
&=U_{\mathop{\mathrm{P}MNS}}\begin{pmatrix}\Theta_{14}&\Theta_{15}\\\
\Theta_{24}&\Theta_{25}\\\ \Theta_{34}&\Theta_{35}\\\
\end{pmatrix}\begin{pmatrix}-1&0\\\
0&-1\end{pmatrix}\begin{pmatrix}\Theta_{41}^{T}&\Theta_{42}^{T}&\Theta_{43}^{T}\\\
\Theta_{51}^{T}&\Theta_{52}^{T}&\Theta_{53}^{T}\end{pmatrix}U_{\mathop{\mathrm{P}MNS}}^{T}\\\
&=-U_{\mathop{\mathrm{P}MNS}}\begin{pmatrix}\Theta_{14}\Theta_{41}^{T}+\Theta_{15}\Theta_{51}^{T}&\Theta_{14}\Theta_{42}^{T}+\Theta_{15}\Theta_{52}^{T}&\Theta_{14}\Theta_{43}^{T}+\Theta_{15}\Theta_{53}^{T}\\\
\Theta_{24}\Theta_{41}^{T}+\Theta_{25}\Theta_{51}^{T}&\Theta_{24}\Theta_{42}^{T}+\Theta_{25}\Theta_{52}^{T}&\Theta_{24}\Theta_{43}^{T}+\Theta_{25}\Theta_{53}^{T}\\\
\Theta_{34}\Theta_{41}^{T}+\Theta_{35}\Theta_{51}^{T}&\Theta_{34}\Theta_{42}^{T}+\Theta_{35}\Theta_{52}^{T}&\Theta_{34}\Theta_{43}^{T}+\Theta_{35}\Theta_{53}^{T}\\\
\end{pmatrix}U_{\mathop{\mathrm{P}MNS}}^{T}\\\ &\simeq-
U_{\mathop{\mathrm{P}MNS}}\begin{pmatrix}\Theta_{14}\Theta_{41}^{T}&\Theta_{14}\Theta_{42}^{T}&\Theta_{14}\Theta_{43}^{T}\\\
\Theta_{24}\Theta_{41}^{T}&\Theta_{24}\Theta_{42}^{T}&\Theta_{24}\Theta_{43}^{T}\\\
\Theta_{34}\Theta_{41}^{T}&\Theta_{34}\Theta_{42}^{T}&\Theta_{34}\Theta_{43}^{T}\\\
\end{pmatrix}U_{\mathop{\mathrm{P}MNS}}^{T}\\\
&=-2U_{\mathop{\mathrm{P}MNS}}\begin{pmatrix}\eta_{11}&\eta_{12}&\eta_{13}\\\
\eta_{21}&\eta_{22}&\eta_{23}\\\ \eta_{31}&\eta_{32}&\eta_{33}\\\
\end{pmatrix}U_{\mathop{\mathrm{P}MNS}}^{T},\end{split}$ (29)

where at the fourth equality the latter terms including index $5$ can be
safely ignored due to the relative smallness of $v_{d}$ compared to $v_{u}$
and the suppression factor $\epsilon$ Hernandez:2021tii . The experimental
upper limits for the deviation of unitarity $\eta_{ij}$ at $1\sigma$ are given
by Fernandez-Martinez:2016lgt :

| NH | IH
---|---|---
$\eta_{11}$ | $4.2\times 10^{-4}$ | $4.8\times 10^{-4}$
$\eta_{22}$ | $2.9\times 10^{-7}$ | $2.4\times 10^{-7}$
$\eta_{33}$ | $9.2\times 10^{-4}$ | $6.8\times 10^{-4}$
$\eta_{12}$ | $8.4\times 10^{-6}$ | $8.4\times 10^{-6}$
$\eta_{13}$ | $6.5\times 10^{-4}$ | $6.5\times 10^{-4}$
$\eta_{23}$ | $2.5\times 10^{-5}$ | $1.3\times 10^{-5}$

Table 2: Experimental upper limits for the deviation of unitarity $\eta$. The
off-diagonal elements of $\eta_{ij}$ where $i\neq j$ are calculated by the
Schwarz inequality $\eta_{ij}\lesssim\sqrt{\eta_{ii}\eta_{jj}}$

Taking the experimental neutrino mixing angles from NuFIT 5.0 Esteban:2020cvm
, they read off (NH):

$\begin{split}\sin^{2}\theta_{12}&=0.304,\quad\sin^{2}\theta_{23}=0.570,\quad\sin^{2}\theta_{13}=0.02221\\\\[4.30554pt]
U_{\mathop{\mathrm{P}MNS}}&=U_{23}U_{13}U_{12}\\\\[4.30554pt]
U_{\mathop{\mathrm{P}MNS}}&=\left(\begin{array}[]{ccc}0.82495&-0.545205&-0.14903\\\
0.267684&0.609102&-0.746552\\\ 0.497798&0.575975&0.648421\\\
\end{array}\right)\end{split}$ (30)

From Equation 30, we can calculate the $Z^{\prime}$ coupling constant to muon
neutrino pair when $g_{X}$ is assumed to be order of unity:

$g_{\nu_{\mu}\nu_{\mu}}\simeq-5.263\times 10^{-4}.$ (31)

##  V LEPTON SECTOR PHENOMENOLOGY

In this section, we analyze the SM observables in the charged lepton sector
with the $Z^{\prime}$ coupling constants defined in section IV. There are many
$Z^{\prime}$ models in the particle physics community, however most of them do
not lead to a definite prediction as the $Z^{\prime}$ coupling constants and
its mass are not determined experimentally. In spite of this drawback, the
$Z^{\prime}$ has gotten a lot of attention since it can be a new possible FCNC
source and thus might be able to solve the well-known anomalies such as
$R_{K(K^{*})}$. We studied the FCNC interactions mediated by the SM $Z$ gauge
boson with one complete vectorlike family CarcamoHernandez:2021yev  in the
same model frame except that the $U(1)^{\prime}$ local symmetry is replaced by
the global $U(1)^{\prime}$ symmetry and what we found there is the SM $Z$ FCNC
interactions generally give suppressed sensitivities for the diverse CLFV
decays ($\tau\rightarrow\mu\gamma,\tau\rightarrow 3\mu$, and
$Z\rightarrow\mu\tau$), compared to their experimental bounds. In this time,
we assume that the CLFV decays are mainly mediated by the $Z^{\prime}$ gauge
boson and the first task to do is to determine the highest order of the
$Z^{\prime}$ coupling constant $g_{X}$ and the coupling constant can be
constrained by investigating both the muon $g-2$ and the CLFV decays.

###  V.1 Muon anomalous magnetic moment $g-2$

The first observable we consider in the charged lepton sector is the muon
anomalous magnetic moment $g-2$, regarded as one of the most important new
physics signal with $4.2\sigma$ deviation Muong-2:2021ojo . The muon $g-2$ in
the model under consideration gets the leading order contributions from the
one loop diagrams in Figure 5:

Figure 5: Feynman diagrams contributing to the muon anomalous magnetic moment
$g-2$ at one loop level. The cross notation in each diagram means the helicity
flip process.

The corresponding analytic expression for the muon $g-2$ reads off
Raby:2017igl ; CarcamoHernandez:2019ydc :

$\begin{split}\Delta
a_{\mu}^{Z^{\prime}}=-\frac{m_{\mu}^{2}}{8\pi^{2}M_{Z^{\prime}}^{2}}\bigg{[}&\big{(}|g_{\mu\mu}^{L}|^{2}+|g_{\mu\mu}^{R}|^{2}\big{)}F(x_{\mu})+\big{(}|g_{\mu\tau}^{L}|^{2}+|g_{\mu\tau}^{R}|^{2}\big{)}F(x_{\tau})\\\
&+\big{(}|g_{\mu E_{4}}^{L}|^{2}+|g_{\mu
E_{4}}^{R}|^{2}\big{)}F(x_{E_{4}})+\big{(}|g_{\mu\widetilde{E}_{4}}^{L}|^{2}+|g_{\mu\widetilde{E}_{4}}^{R}|^{2}\big{)}F(x_{\widetilde{E}_{4}})\\\
&+\operatorname{Re}\big{(}g_{\mu\mu}^{L}g_{\mu\mu}^{R*}\big{)}\frac{m_{\mu}}{m_{\mu}}G(x_{\mu})+\operatorname{Re}\big{(}g_{\mu\tau}^{L}g_{\mu\tau}^{R*}\big{)}\frac{m_{\tau}}{m_{\mu}}G(x_{\tau})\\\
&+\operatorname{Re}\big{(}g_{\mu
E_{4}}^{L}g_{\mu\widetilde{E}_{4}}^{R*}\big{)}\frac{M_{E_{4}}}{m_{\mu}}G(x_{E_{4}})+\operatorname{Re}\big{(}g_{\mu\widetilde{E}_{4}}^{L}g_{\mu
E_{4}}^{R*}\big{)}\frac{M_{\widetilde{E}_{4}}}{m_{\mu}}G(x_{\widetilde{E}_{4}})\bigg{]},\end{split}$
(32)

where the coupling constant $g_{ij}^{{L}(R)}$ means the left-handed (right-
handed) $Z^{\prime}$ coupling constant to the particle $i$ and $j$, and $F$
and $G$ are the loop functions defined as follows
($m_{\mathop{\mathrm{l}oop}}$ means mass of the particle running in the loop):

$\begin{split}F(x)&=\frac{5x^{4}-14x^{3}+39x^{2}-38x-18x^{2}\mathop{\mathrm{l}n}x+8}{12(1-x)^{4}}\\\\[4.30554pt]
G(x)&=\frac{x^{3}+3x-6x\mathop{\mathrm{l}n}x-4}{2(1-x)^{3}},\qquad
x=\frac{m_{\mathop{\mathrm{l}oop}}^{2}}{M_{Z}^{\prime 2}}.\end{split}$ (33)

The most dominant contributions in the analytic expression of Equation 32 come
from the terms involving the enhancement factors $M_{E_{4}}/m_{\mu}$ or
$M_{\widetilde{E}_{4}}/m_{\mu}$. However, these terms can not keep increasing
as the vectorlike masses get heavier since the relevant coupling constants get
weaker at the same time, so these terms yield a balanced contribution to the
muon $g-2$ prediction and this approach is also applied to the CLFV
$\tau\rightarrow\mu\gamma$ decay as we discuss in the next subsection. The
most recent experimental bound for the muon $g-2$ at $1\sigma$ is given by
FNAL Muong-2:2021ojo :

$\begin{split}\Delta
a_{\mu}&=a_{\mu}^{\mathop{\mathrm{e}xp}}-a_{\mu}^{\mathop{\mathrm{S}M}}=\left(25.1\pm
5.9\right)\times 10^{-10},\end{split}$ (34)

###  V.2 CLFV $\tau\rightarrow\mu\gamma$ decay

As this BSM model features the massive second and third generation of the SM,
the most stringent constraint in the charged lepton sector comes from the CLFV
$\tau\rightarrow\mu\gamma$ decay. The leading order contributions to the CLFV
$\tau\rightarrow\mu\gamma$ decay mainly arise at one loop level as given in
Figure 6:

Figure 6: Feynman diagrams contributing to the CLFV $\tau\rightarrow\mu\gamma$
decay at one loop level. The cross notation in each diagram means the helicity
flip process.

The corresponding analytic expression for the branching ratio of the
$\tau\rightarrow\mu\gamma$ decay reads CarcamoHernandez:2019ydc ; Raby:2017igl
; Lindner:2016bgg ; Lavoura:2003xp ; Chiang:2011cv :

$\begin{split}\mathop{\mathrm{B}R}\left(\tau\rightarrow\mu\gamma\right)&=\frac{\alpha_{\mathop{\mathrm{e}m}}}{1024\pi^{4}}\frac{m_{\tau}^{5}}{M_{Z}^{\prime
4}\Gamma_{\tau}}\Big{(}\lvert
g_{\tau\mu}^{L}g_{\mu\mu}^{L}F(x_{\mu})+g_{\tau\tau}^{L}g_{\tau\mu}^{L}F(x_{\tau})+g_{\tau
E_{4}}^{L}g_{E_{4}\mu}^{L}F(x_{E_{4}})+g_{\tau\widetilde{E}_{4}}^{L}g_{\widetilde{E}_{4}\mu}^{L}F(x_{\widetilde{E}_{4}})\\\
&+\frac{m_{\mu}}{m_{\tau}}g_{\tau\mu}^{L}g_{\mu\mu}^{L}F(x_{\mu})+\frac{m_{\mu}}{m_{\tau}}g_{\tau\tau}^{L}g_{\tau\mu}^{L}F(x_{\tau})+\frac{m_{\mu}}{m_{\tau}}g_{\tau
E_{4}}^{L}g_{E_{4}\mu}^{L}F(x_{E_{4}})+\frac{m_{\mu}}{m_{\tau}}g_{\tau\widetilde{E}_{4}}^{L}g_{\widetilde{E}_{4}\mu}^{L}F(x_{\widetilde{E}_{4}})\\\
&+\frac{m_{\mu}}{m_{\tau}}g_{\tau\mu}^{L}g_{\mu\mu}^{R}G(x_{\mu})+\frac{m_{\tau}}{m_{\tau}}g_{\tau\tau}^{L}g_{\tau\mu}^{R}G(x_{\mu})+\frac{M_{E_{4}}}{m_{\tau}}g_{\tau
E_{4}}^{L}g_{\widetilde{E}_{4}\mu}^{R}G(x_{E_{4}})+\frac{M_{\widetilde{E}_{4}}}{m_{\tau}}g_{\tau\widetilde{E}_{4}}^{L}g_{E_{4}\mu}^{R}G(x_{\widetilde{E}_{4}})\rvert^{2}\\\
&+\lvert
g_{\tau\mu}^{R}g_{\mu\mu}^{R}F(x_{\mu})+g_{\tau\tau}^{R}g_{\tau\mu}^{R}F(x_{\tau})+g_{\tau
E_{4}}^{R}g_{E_{4}\mu}^{R}F(x_{E_{4}})+g_{\tau\widetilde{E}_{4}}^{R}g_{\widetilde{E}_{4}\mu}^{R}F(x_{\widetilde{E}_{4}})\\\
&+\frac{m_{\mu}}{m_{\tau}}g_{\tau\mu}^{R}g_{\mu\mu}^{R}F(x_{\mu})+\frac{m_{\mu}}{m_{\tau}}g_{\tau\tau}^{R}g_{\tau\mu}^{R}F(x_{\tau})+\frac{m_{\mu}}{m_{\tau}}g_{\tau
E_{4}}^{R}g_{E_{4}\mu}^{R}F(x_{E_{4}})+\frac{m_{\mu}}{m_{\tau}}g_{\tau\widetilde{E}_{4}}^{R}g_{\widetilde{E}_{4}\mu}^{R}F(x_{\widetilde{E}_{4}})\\\
&+\frac{m_{\mu}}{m_{\tau}}g_{\tau\mu}^{R}g_{\mu\mu}^{L}G(x_{\mu})+\frac{m_{\tau}}{m_{\tau}}g_{\tau\tau}^{R}g_{\tau\mu}^{L}G(x_{\mu})+\frac{M_{\widetilde{E}_{4}}}{m_{\tau}}g_{\tau
E_{4}}^{R}g_{\widetilde{E}_{4}\mu}^{L}G(x_{\widetilde{E}_{4}})+\frac{M_{E_{4}}}{m_{\tau}}g_{\tau\widetilde{E}_{4}}^{R}g_{E_{4}\mu}^{L}G(x_{E_{4}})\rvert^{2}\Big{)},\end{split}$
(35)

where $\alpha_{\mathop{\mathrm{e}m}}$ is the fine structure constant,
$\Gamma_{\tau}$ is the total decay width of the tau lepton
($\Gamma_{\tau}=5\Gamma\left(\tau_{L}^{-}\rightarrow\nu_{\tau}e_{L}^{-}\overline{\nu}_{e}\right))=2.0\times
10^{-12}\mathop{\mathrm{G}eV}$). The most recent experimental upper-limit for
the CLFV $\tau\rightarrow\mu\gamma$ decay is given by Crivellin:2020ebi ;
MEG:2016leq ; BaBar:2009hkt :

$\mathop{\mathrm{B}R}\left(\tau\rightarrow\mu\gamma\right)_{\mathop{\mathrm{e}xp}}=4.4\times
10^{-8}.$ (36)

###  V.3 CLFV $\tau\rightarrow\mu\mu\mu$ decay

The next observable we consider is the CLFV $\tau\rightarrow\mu\mu\mu$ decay.
The CLFV $\tau\rightarrow 3\mu$ decay mediated by the neutral massive
$Z^{\prime}$ gauge boson can have the leading order contributions at tree
level as given in Figure 7:

Figure 7: Feynman diagrams contributing to the CLFV $\tau\rightarrow 3\mu$
decay at tree level.

We derived the analytic expression for the branching ratio of the CLFV
$\tau\rightarrow 3\mu$ decay mediated by the SM $Z$ gauge boson in one of our
works CarcamoHernandez:2021yev , so we can simply make use of the derived
result after replacing the SM $Z$ gauge boson in the prediction by the neutral
massive $Z^{\prime}$ gauge boson as follows:

$\mathop{\mathrm{B}R}\left(\tau\rightarrow
3\mu\right)=\frac{1}{64\pi^{3}m_{\tau}\Gamma_{\tau}}\int_{m_{\mu}^{2}/m_{\tau}}^{\frac{1}{2}m_{\tau}}\int_{\frac{1}{2}m_{\tau}+\frac{m_{\mu}^{2}}{m_{\tau}}-E_{2}}^{\frac{1}{2}m_{\tau}}\left(\frac{1}{2}\sum_{\mathop{\mathrm{s}pin}}\lvert\mathcal{M}\left(g_{1},g_{2},E_{2},E_{4}\right)\rvert^{2}\right)dE_{4}dE_{2},$
(37)

where the amplitude squared is defined by Hernandez:2021tii
($g_{1}=g_{\mu\tau}^{L2}g_{\mu\mu}^{L2}+g_{\mu\tau}^{R2}g_{\mu\mu}^{R2},g_{2}=g_{\mu\tau}^{L2}g_{\mu\mu}^{R2}+g_{\mu\tau}^{R2}g_{\mu\mu}^{L2}$)

$\begin{split}\frac{1}{2}\sum_{\mathop{\mathrm{s}pin}}\lvert\mathcal{M}\left(g_{1},g_{2},E_{2},E_{4}\right)\rvert^{2}&=\frac{4}{M_{Z^{\prime}}^{4}}\Big{[}g_{1}\left(m_{\tau}^{2}-m_{\tau}(E_{2}+E_{4})\right)\left(-m_{\tau}^{2}-m_{\mu}^{2}+2m_{\tau}(E_{2}+E_{4})\right)\\\
&+g_{2}\left(m_{\tau}E_{4}\right)\left(m_{\tau}^{2}-m_{\mu}^{2}-2m_{\tau}E_{4}\right)\Big{]}\end{split}$
(38)

and $E_{2}$ and $E_{4}$ are the energies defined in the rest frame as follows:

$\begin{split}p_{1}&=(m_{\tau},\vec{0})\\\
p_{2}&=\left(E_{2},\vec{p}_{2}\right)\\\
p_{3}&=\left(m_{\tau}-E_{2}-E_{4},-\vec{p}_{2}-\vec{p}_{4}\right)\\\
p_{4}&=\left(E_{4},\vec{p}_{4}\right)\\\ \end{split}$ (39)

and the energy $E_{4}$ runs from
$\frac{1}{2}m_{\tau}+\frac{m_{\mu}^{2}}{m_{\tau}}-E_{2}$ to
$\frac{1}{2}m_{\tau}$ and the other energy $E_{2}$ runs from
$\frac{m_{\mu}^{2}}{m_{\tau}}$ to $\frac{1}{2}m_{\tau}$ as seen from the
integration range of Equation 37. The experimental upper-limit of the
branching ratio of the CLFV $\tau\rightarrow 3\mu$ decay is given by
Hayasaka:2010np :

$\mathop{\mathrm{B}R}\left(\tau\rightarrow
3\mu\right)_{\mathop{\mathrm{e}xp}}=2.1\times 10^{-8}$ (40)

###  V.4 Neutrino trident production

One more constraint for the $Z^{\prime}$ gauge boson in the lepton sector
arise from the neutrino trident production. The relevant diagram for the
neutrino trident production reads in Figure 8:

Figure 8: Feynman diagram for the neutrino trident production.

The constraint for the neutrino trident production is given by the effective
four fermion interaction by exchanging the neutral massive $Z^{\prime}$ gauge
boson as follows CarcamoHernandez:2019ydc ; CHARM-II:1990dvf ; CCFR:1991lpl ;
Altmannshofer:2014pba :

$\mathcal{L}_{\nu\mathop{\mathrm{T}ri}}^{Z^{\prime},\mathop{\mathrm{e}ff}}=-\frac{g_{\nu_{\mu}\nu_{\mu}}^{L}g_{\mu\mu}^{L}}{M_{Z^{\prime}}^{2}}\left(\overline{\mu}_{L}\gamma_{\mu}\mu_{L}\right)\left(\overline{\nu}_{\mu
L}\gamma^{\mu}\nu_{\mu
L}\right)-\frac{g_{\nu_{\mu}\nu_{\mu}}^{L}g_{\mu\mu}^{R}}{M_{Z^{\prime}}^{2}}\left(\overline{\mu}_{R}\gamma_{\mu}\mu_{R}\right)\left(\overline{\nu}_{\mu
L}\gamma^{\mu}\nu_{\mu L}\right)$ (41)

The experimental bound for the neutrino trident production at $2\sigma$ reads
Falkowski:2017pss :

$\begin{split}-\frac{1}{\left(390\mathop{\mathrm{G}eV}\right)^{2}}&<\frac{g_{\nu_{\mu}\nu_{\mu}}^{L}g_{\mu\mu}^{L}+g_{\nu_{\mu}\nu_{\mu}}^{L}g_{\mu\mu}^{R}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(370\mathop{\mathrm{G}eV}\right)^{2}}\\\
-\frac{1}{\left(390\mathop{\mathrm{G}eV}\right)^{2}}&<\frac{g_{\nu_{\mu}\nu_{\mu}}^{L}\left(g_{\mu\mu}^{L}+g_{\mu\mu}^{R}\right)}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(370\mathop{\mathrm{G}eV}\right)^{2}}\end{split}$
(42)

It is possible to simplify Equation 42 further by using the calculated
coupling constant $g_{\nu_{\mu}\nu_{\mu}}$ of Equation 31:

$-\frac{1}{\left(390\mathop{\mathrm{G}eV}\right)^{2}}<\frac{-5.263\times
10^{-4}\left(g_{\mu\mu}^{L}+g_{\mu\mu}^{R}\right)}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(370\mathop{\mathrm{G}eV}\right)^{2}},$
(43)

where it is worth mentioning that the coupling constant
$g_{\nu_{\mu}\nu_{\mu}}$ is given when we assume the $Z^{\prime}$ coupling
constant $g_{X}$ is the order of unity. As we will see soon in the numerical
analysis for the charged lepton sector, the coupling constant $g_{X}$ can
reach its highest order up to unity (or equally $g_{X}=1$). Therefore, we can
determine the $Z^{\prime}$ mass range from the neutrino trident production
constraint once a numerical range of summing over the left-handed and right-
handed muon coupling constant is given by the numerical analysis and this will
be discussed in the next subsection.

###  V.5 Numerical analysis in the charged lepton sector

In this subsection, we carry out numerical scan in the charged lepton sector.
What we want to achieve via this numerical scan is to determine the highest
order of the $Z^{\prime}$ coupling constant $g_{X}$ and this process is quite
important since the determined order of the coupling constant $g_{X}$ affects
the other analyses which will be carried out in the quark sector as well as in
the electroweak precision observables. The possible range of the $Z^{\prime}$
mass can be constrained in the quark sector via $R_{K^{*}}$ anomaly as well as
$B_{s}$ meson oscillation and we leave constraining the $Z^{\prime}$ mass
until we carry out the quark sector investigation and regard the $Z^{\prime}$
mass as a free parameter in this subsection.

####  V.5.1 Free parameter setup

The free parameter setup given in Table 3 follows the exactly same form
defined in one of our works CarcamoHernandez:2021yev  except for the
$Z^{\prime}$ mass:

Mass parameter | Scanned Region($\mathop{\mathrm{G}eV}$)
---|---
$y_{24}^{e}v_{d}=m_{24}$ | $\pm[1,10]$
$y_{34}^{e}v_{d}=m_{34}$ | $\pm[1,10]$
$y_{43}^{e}v_{d}=m_{43}$ | $\pm[1,10]$
$x_{34}^{L}v_{\phi}=m_{35}$ | $\pm[50,200]$
$x_{42}^{e}v_{\phi}=m_{52}$ | $\pm[50,200]$
$x_{43}^{e}v_{\phi}=m_{53}$ | $\pm[50,200]$
$M_{45}^{L}$ | $\pm[150,2000]$
$M_{54}^{e}$ | $\pm[150,2000]$
$M_{Z^{\prime}}$ | $[75,3000]$

Table 3: Parameter setup to determine the highest order of the $Z^{\prime}$
coupling constant $g_{X}$ by the muon $g-2$ and the CLFV
$\tau\rightarrow\mu\tau$ and $\tau\rightarrow 3\mu$ decays.

The most relevant features about the parameter setup of Table 3 are summarized
as follows:

  1. 1.

As the model under consideration features an extended 2HDM, the up-type Higgs
vev $v_{u}$, very close to $246\mathop{\mathrm{G}eV}$, and the down-type Higgs
vev $v_{d}$, assumed to run from $1$ to $10\mathop{\mathrm{G}eV}$, satisfy the
relation $v_{u}^{2}+v_{d}^{2}=\left(246\mathop{\mathrm{G}eV}\right)^{2}$. The
down-type Higgs $H_{d}$ interacts with the down-type quarks as well as the
charged leptons and thus the mass parameters involving the down-type Higgs vev
$v_{d}$ are determined under the assumption.

  2. 2.

The singlet flavon vev $v_{\phi}$ and vectorlike masses $M_{45}^{L}$ and
$M_{54}^{e}$ and lastly the $Z^{\prime}$ mass $M_{Z^{\prime}}$ are free
parameters and mass parameters including the flavon vev $v_{\phi}$ are
determined to be consistent with the observed charged lepton hierarchy. Here a
main role of the free parameter $M_{Z^{\prime}}$ is to see whether it can be
constrained by the muon $g-2$ and CLFV $\tau\rightarrow\mu\gamma$ and
$\tau\rightarrow 3\mu$ decays. As for the lower bound of the $Z^{\prime}$ mass
in Table 3, it is set to avoid the low mass constraint of $Z^{\prime}$ mass
$5\lesssim M_{Z^{\prime}}\lesssim 70\mathop{\mathrm{G}eV}$ resulting from the
decay process $pp\rightarrow Z\rightarrow 4\mu$ Navarro:2021sfb ;
Falkowski:2018dsl ; Altmannshofer:2014pba ; Altmannshofer:2014cfa ;
Altmannshofer:2016jzy .

  3. 3.

What we constrain in this numerical scan is the muon and tau mass plus the
$23$ mixing angle. The theoretically predicted muon and tau mass should be put
between $\left[1\pm 0.1\right]\times m_{\mu,\tau}^{\mathop{\mathrm{e}xp}}$. On
top of that, we also constrain the $23$ mixing angle and determine the upper-
limit for the mixing angle to be $0.2$ since the sizeable off-diagonal
components of the PMNS mixing matrix mainly arise from the neutrino sector.

####  V.5.2 Numerical scan result for the charged lepton sector

The numerical scan result for the muon $g-2$ and for the CLFV
$\tau\rightarrow\mu\gamma,\tau\rightarrow 3\mu$ decays when the $Z^{\prime}$
coupling constant $g_{X}=1$ is given in Figure 9:

Figure 9: Numerical scan result for $Z^{\prime}$ mass, free parameter, and for
the muon $g-2$ as well as the CLFV $\tau\rightarrow\mu\gamma,\tau\rightarrow
3\mu$ decays when the $Z^{\prime}$ coupling constant $g_{X}$ is fixed by
unity. $M_{E_{4}}$ appearing in $x$ axis is mass of the vectorlike doublet
charged lepton $E_{4}$ and $M_{\widetilde{E}_{4}}$ is mass of the vectorlike
singlet charged lepton $\widetilde{E}_{4}$. The darker blue region in each
plot indicates that the vectorlike doublet or singlet charged lepton mass
should be excluded up to $200\mathop{\mathrm{G}eV}$ Xu:2018pnq ;
Guedes:2021oqx . The brighter blue region in each plot means that the
vectorlike doublet mass $M_{E_{4}}$ should be excluded up to
$790\mathop{\mathrm{G}eV}$ by the CMS experimental constraint
Bhattiprolu:2019vdu ; CMS:2019hsm .

The first feature we discuss in Figure 9 is the $Z^{\prime}$ mass. Here, the
scanned $Z^{\prime}$ mass range appears to be the same as the initial setup
$\left[75,3000\mathop{\mathrm{G}eV}\right]$, which means that the charged
lepton observables (muon $g-2$, the CLFV
$\tau\rightarrow\mu\gamma,\tau\rightarrow 3\mu$ decays) can not constrain the
free parameter $M_{Z^{\prime}}$ at all. However, constraining the $Z^{\prime}$
mass can be done by constraints arising in the quark sector. The second is the
scanned muon $g-2$. Interestingly, the scanned muon $g-2$ in this BSM reports
only negative contributions with order of mostly $10^{-11}$, which implies
that the $Z^{\prime}$ physics is definitely not an answer for the muon anomaly
since it only deteriorates our muon prediction and this feature is one of big
differences, compared to Navarro:2021sfb , as we can not have positive
contributions to the muon $g-2$ regardless of how large the $Z^{\prime}$
coupling constant to right-handed muon pair is. We will try to explain the
muon $g-2$ via scalar exchange and this will be discussed in the next section.
Next we discuss our predictions for the CLFV $\tau\rightarrow\mu\gamma$ decay.
Considering its logarithmic experimental bound,
$\mathop{\mathrm{L}og}_{10}\mathop{\mathrm{B}R}\left(\tau\rightarrow\mu\gamma\right)\approx-7.4$
Crivellin:2020ebi ; MEG:2016leq ; BaBar:2009hkt , our predictions for the
$\tau\rightarrow\mu\gamma$ decay are not significantly constrained by its
experimental constraint, implying that the $Z^{\prime}$ coupling constant
$g_{X}$ can be order of unity at most. The last CLFV observable we discuss is
the $\tau\rightarrow 3\mu$ decay and its situation is similar to the CLFV
$\tau\rightarrow\mu\gamma$ decay in the sense that it can not constrain
$Z^{\prime}$ gauge boson at all. Excluding the benchmark points, which exceed
the experimental upper-limit of the $\tau\rightarrow\mu\gamma$ decay, the
result is given in Figure 10:

Figure 10: Numerical scan result for the charged lepton observables after
excluding the benchmark points which exceed the upper-limit of the
$\tau\rightarrow\mu\gamma$ decay. After excluding the benchmark points which
exceed the experimental upper-limit of $\tau\rightarrow\mu\gamma$ decay, most
of other benchmark points survive and this implies that the highest order of
coupling constant $g_{X}=1$ is not ruled out by the charged lepton
observables.

From the survived benchmark points of Figure 10, we can derive a numerical
range of the $Z^{\prime}$ coupling constant to the left- or right-handed muon
pair $g_{\mu\mu}^{L,R}$ as follows:

Coupling constant | Scanned range
---|---
$g_{\mu\mu}^{L}$ | $4.888\times 10^{-6}<g_{\mu\mu}^{L}<2.113\times 10^{-3}$
$g_{\mu\mu}^{R}$ | $-1.195\times 10^{-2}<g_{\mu\mu}^{R}<-3.756\times 10^{-4}$
$g_{\mu\mu}^{L}+g_{\mu\mu}^{R}$ | $-1.175\times 10^{-2}<g_{\mu\mu}^{L}+g_{\mu\mu}^{R}<1.179\times 10^{-3}$

Table 4: Scanned range of the $Z^{\prime}$ coupling constant to left-handed
(or right-handed) muon pair when the $Z^{\prime}$ coupling constant $g_{X}$ is
set to be order of unity, required to derive mass range of $Z^{\prime}$ gauge
boson from the neutrino trident constraint ($g_{\mu\mu}^{L}+g_{\mu\mu}^{R}$)
as well as the $R_{K^{*}}$ anomaly ($g_{\mu\mu}^{L}$ or $g_{\mu\mu}^{R}$).

where it should be clarified that the upper-limit of the left-handed muon pair
coupling constant $g_{L}^{\mu\mu}$ is given by ruling out the benchmark points
which exceed the CMS and the $\tau\rightarrow\mu\gamma$ experimental
constraints. Without taking the CMS constraint into account, the left-handed
coupling constant $g_{L}^{\mu\mu}$ can reach up to a very similar magnitude of
the lower-limit of the right-handed coupling constant,
$g_{L}^{\mu\mu}=1.339\times 10^{-2}$, and we can connect the physical
$Z^{\prime}$ gauge boson and the CP-odd scalar arising in this BSM model via
the scalar potential in this case as we will see in the electroweak numerical
study. However, the upper-limit of the left-handed coupling constant must be
suppressed from $1.339\times 10^{-2}$ to $2.113\times 10^{-3}$ by the CMS
experimental constraint and what this implements is connecting the
$Z^{\prime}$ mass to the CP-odd scalar becomes much more challenging as
$Z^{\prime}$ mass range derived from the $R_{K^{*}}$ anomaly gets lower
further.

####  V.5.3 $Z^{\prime}$ mass range derived from the neutrino trident
constraint

As we have derived a numerical range of summing over the $Z^{\prime}$ coupling
constants to muon pair in Table 4, it can evaluate $Z^{\prime}$ mass range
derived from the neutrino trident constraint. Reminding the neutrino trident
constraint of Equation 43, it reads:

$-\frac{1}{\left(390\mathop{\mathrm{G}eV}\right)^{2}}<\frac{-5.263\times
10^{-4}\left(g_{\mu\mu}^{L}+g_{\mu\mu}^{R}\right)}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(370\mathop{\mathrm{G}eV}\right)^{2}}.$
(44)

Using the derived result of $g_{\mu\mu}^{L}+g_{\mu\mu}^{R}$, it is possible to
determine $Z^{\prime}$ mass range as given in Table 5:

Coupling constant |  $M_{Z^{\prime}}$ range
---|---
$g_{\mu\mu}^{L}+g_{\mu\mu}^{R}=-1.175\times 10^{-2}$ | $M_{Z^{\prime}}>0.920\mathop{\mathrm{G}eV}$
$g_{\mu\mu}^{L}+g_{\mu\mu}^{R}=1.179\times 10^{-3}$ | $M_{Z^{\prime}}>0.307\mathop{\mathrm{G}eV}$

Table 5: $M_{Z^{\prime}}$ mass range derived from the neutrino trident
production constraint with the $Z^{\prime}$ coupling constants to muon pair
summed over.

Considering both results for the derived $Z^{\prime}$ mass range in Table 5,
it concludes

$M_{Z^{\prime}}>0.920\mathop{\mathrm{G}eV}.$ (45)

##  VI QUARK SECTOR PHENOMENOLOGY

Now that we have defined the preferred order of the $Z^{\prime}$ coupling
constant $g_{X}$, the next task is to constrain mass range of the $Z^{\prime}$
and it can be done by considering constraints arising in the quark sector,
which are the $R_{K^{*}}$ anomaly, $B_{s}$ meson oscillation, collider, and
CKM mixing matrix.

###  VI.1 $R_{K^{*}}$ anomaly

The first observable we consider is the $R_{K^{*}}$ anomaly. The $R_{K^{*}}$
anomaly is also regarded as a potential new physics signal in the quark sector
due to its sizeable SM deviation corresponding to $3.1\sigma$. The $R_{K^{*}}$
anomaly and the $B_{s}\rightarrow\overline{\mu}\mu$ decay, connected to the
$R_{K^{*}}$ anomaly by the crossing symmetry, can be explained by the
effective four fermion operator consisting of left- and right-handed muon pair
and left-handed $b$ and $s$ quarks after integrating out the degree of freedom
of $Z^{\prime}$ gauge boson as follows Navarro:2021sfb :

$\mathcal{L}_{R_{K^{*}}}^{\mathop{\mathrm{e}ff}}=G_{bs\mu}^{L}\left(\overline{s}_{L}\gamma^{\mu}b_{L}\right)\left(\overline{\mu}_{L}\gamma_{\mu}\mu_{L}\right)+G_{bs\mu}^{R}\left(\overline{s}_{L}\gamma^{\mu}b_{L}\right)\left(\overline{\mu}_{R}\gamma_{\mu}\mu_{R}\right)+\mathop{\mathrm{H}.c.},$
(46)

where the coupling constants $G_{bs\mu}^{L,R}$ are known as the Wilson
coefficients and are defined by Navarro:2021sfb :

$\begin{split}G_{bs\mu}^{L}&=-\frac{g_{bs}^{L}g_{\mu\mu}^{L}}{M_{Z^{\prime}}^{2}},\\\
G_{bs\mu}^{R}&=-\frac{g_{bs}^{L}g_{\mu\mu}^{R}}{M_{Z^{\prime}}^{2}}.\end{split}$
(47)

The Feynman diagrams contributing to the $R_{K^{*}}$ anomaly and the
$B_{s}\rightarrow\overline{\mu}\mu$ decay read in Figure 11:

Figure 11: Feynman diagrams contributing to the $R_{K^{*}}$ anomaly (left) and
the $B_{s}\rightarrow\overline{\mu}\mu$ decay (right).

The effective interactions in Equation 46 can be rewritten in terms of the
vector and axial-vector effective operators Navarro:2021sfb ; Geng:2021nhg :

$\begin{split}\mathcal{H}_{\mathop{\mathrm{e}ff}}=N\left[\delta
C_{9}\left(\overline{s}_{L}\gamma^{\mu}b_{L}\right)\left(\overline{\mu}\gamma_{\mu}\mu\right)+\delta
C_{10}\left(\overline{s}_{L}\gamma^{\mu}b_{L}\right)\left(\overline{\mu}\gamma_{\mu}\gamma^{5}\mu\right)\right]+\mathop{\mathrm{H}.c.},\end{split}$
(48)

where $N$ is the normalization factor Navarro:2021sfb :

$N=-\frac{G_{F}}{\sqrt{2}}V_{tb}V_{ts}^{*}\frac{e^{2}}{16\pi^{2}}.$ (49)

The coupling constants $\delta C_{9,10}$ can be rearranged in terms of the
Wilson coefficients $G_{bs\mu}^{L,R}$:

$\begin{split}G_{bs\mu}^{L}&=N\left(\delta C_{9}-\delta C_{10}\right)\\\
G_{bs\mu}^{R}&=N\left(\delta C_{9}+\delta C_{10}\right)\end{split}$ (50)

As the BSM model under consideration can gain access to the coupling constants
$g_{bs}^{L}$ and $g_{\mu\mu}^{L,R}$ in person, it requires a numerical range
of the Wilson coefficients $G_{bs\mu}^{L,R}$ and they can be determined by
using the numerical results of $\delta C_{9,10}$ which has been fitted to
explain the $R_{K^{*}}$ anomaly at $1\sigma$. The coefficients $\delta
C_{9,10}$ has been fitted by two cases: one of which is the “theoretically
clean fit”, meaning that this fitted result can explain the $R_{K^{*}}$
anomaly and the $B_{s}\rightarrow\overline{\mu}\mu$ decay in a unified way
without theoretical uncertainties, and the other is the “global fit” which
includes more substantial theoretical uncertainties when compared to the
$B_{s}\rightarrow\overline{\mu}\mu$ data as well as the ratios of lepton
universality violation and the two fits are displayed in Table 6 and 7:

| Best fit | at $1\sigma$
---|---|---
$(\delta C_{9},\delta C_{10})$ | $(-0.11,0.59)$ | $\left(-0.41<\delta C_{9}<0.17\right),\left(0.38<\delta C_{10}<0.81\right)$
$(G_{bs\mu}^{L},G_{bs\mu}^{R})$ | $\left(\frac{1}{\left(42.50\mathop{\mathrm{T}eV}\right)^{2}},-\frac{1}{\left(51.3\mathop{\mathrm{T}eV}\right)^{2}}\right)$ | $\left(\frac{1}{\left(44.44\mathop{\mathrm{T}eV}\right)^{2}}<G_{bs\mu}^{L}<\frac{1}{\left(40\mathop{\mathrm{T}eV}\right)^{2}}\right),\left(-\frac{1}{\left(35.9\mathop{\mathrm{T}eV}\right)^{2}}<G_{bs\mu}^{R}<\frac{1}{\left(205\mathop{\mathrm{T}eV}\right)^{2}}\right)$

Table 6: “Theoretically clean fit” for the $R_{K^{*}}$ anomaly and the
$B_{s}\rightarrow\overline{\mu}\mu$ decay Navarro:2021sfb ; Geng:2021nhg .

| Best fit | at $1\sigma$
---|---|---
$(\delta C_{9},\delta C_{10})$ | $(-0.56,0.30)$ | $\left(-0.79<\delta C_{9}<-0.31\right),\left(0.15<\delta C_{10}<0.49\right)$
$(G_{bs\mu}^{L},G_{bs\mu}^{R})$ | $\left(\frac{1}{\left(38.34\mathop{\mathrm{T}eV}\right)^{2}},-\frac{1}{\left(69.73\mathop{\mathrm{T}eV}\right)^{2}}\right)$ | $\left(\frac{1}{\left(39.75\mathop{\mathrm{T}eV}\right)^{2}}<G_{bs\mu}^{L}<\frac{1}{\left(36.67\mathop{\mathrm{T}eV}\right)^{2}}\right),\left(-\frac{1}{\left(83.8\mathop{\mathrm{T}eV}\right)^{2}}<G_{bs\mu}^{R}<\frac{1}{\left(44.44\mathop{\mathrm{T}eV}\right)^{2}}\right)$

Table 7: “Global fit” for the $R_{K^{*}}$ anomaly and the
$B_{s}\rightarrow\overline{\mu}\mu$ decay Navarro:2021sfb ; Geng:2021nhg .

What we can confirm via the two fits is the different behavior between
$G_{bs\mu}^{L}$ and $G_{bs\mu}^{R}$: the left-handed Wilson coefficients
$G_{bs\mu}^{L}$ at $1\sigma$ in both fits report the order of
$\left(40\mathop{\mathrm{T}eV}\right)^{-2}$ on average even though the
“theoretically clean fit” is slightly tighter than the “global fit”, whereas
the right-handed Wilson coefficient $G_{bs\mu}^{R}$ in the “theoretically
clean fit” prefers negative contributions more while that in the “global fit”
prefers positive contributions further. The inconsistent behavior of the
right-handed Wilson coefficient $G_{bs\mu}^{R}$ in both fits might be hint at
that the $R_{K^{*}}$ anomaly is explained by the only left-handed Wilson
coefficient $G_{bs\mu}^{L}$ and interactions consisting of the only left-
handed fermionic fields King:2017anf ; King:2018fcg ; Falkowski:2018dsl . We
confirm this feature by applying the scanned range of $Z^{\prime}$ coupling
constant to left- or right-handed muon pair as well as left-handed $b$ and $s$
quarks to the Wilson coefficients $G_{bs\mu}^{L,R}$ and then determine a
possible mass range of $Z^{\prime}$. Considering the constraints coming from
the “theoretically clean fit”, they read

$\begin{split}\frac{1}{\left(44.44\mathop{\mathrm{T}eV}\right)^{2}}&<G_{bs\mu}^{L}=-\frac{g_{bs}^{L}g_{\mu\mu}^{L}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(40\mathop{\mathrm{T}eV}\right)^{2}},\\\
-\frac{1}{\left(35.9\mathop{\mathrm{T}eV}\right)^{2}}&<G_{bs\mu}^{R}=-\frac{g_{bs}^{L}g_{\mu\mu}^{R}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(205\mathop{\mathrm{T}eV}\right)^{2}},\\\
\end{split}$ (51)

and the constraints arising from the “global fit” are given by:

$\begin{split}\frac{1}{\left(39.75\mathop{\mathrm{T}eV}\right)^{2}}&<G_{bs\mu}^{L}=-\frac{g_{bs}^{L}g_{\mu\mu}^{L}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(36.67\mathop{\mathrm{T}eV}\right)^{2}},\\\
-\frac{1}{\left(83.8\mathop{\mathrm{T}eV}\right)^{2}}&<G_{bs\mu}^{R}=-\frac{g_{bs}^{L}g_{\mu\mu}^{R}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(44.44\mathop{\mathrm{T}eV}\right)^{2}}.\\\
\end{split}$ (52)

###  VI.2 $B_{s}-\overline{B}_{s}$ meson mixing oscillation

Next we consider the $B_{s}-\overline{B}_{s}$ meson mixing oscillation and
this observable gives rise to the most strict constraint in the quark sector.
As in the $R_{K^{*}}$ anomaly, the effective four fermion interactions after
integrating out the $Z^{\prime}$ gauge boson take the form Navarro:2021sfb :

$\begin{split}\mathcal{L}_{B_{s}-\overline{B}_{s}}^{\mathop{\mathrm{e}ff}}=-\frac{G_{bs}}{2}\left(\overline{s}\gamma_{\mu}b\right)^{2}+\mathop{\mathrm{H}.c.},\end{split}$
(53)

where the coupling constant $G_{bs}$ reads:

$G_{bs}=\frac{g_{bs}^{L2}+g_{bs}^{L}g_{bs}^{R*}+g_{bs}^{R2}}{M_{Z^{\prime}}^{2}},$
(54)

and it requires to clarify that we count all polarized bases, meaning that not
just left-handed $b$ and $s$ quarks but also right-handed $b$ and $s$ quarks
are considered together unlike Navarro:2021sfb . The effective
$B_{s}-\overline{B}_{s}$ meson mixing operators in Equation 53 yield the
following diagrams at tree level in Figure 12:

Figure 12: Feynman diagrams contributing to the $B_{s}-\overline{B}_{s}$ meson
mixing oscillation in the polarized basis.

The experimental mass difference for the $B_{s}-\overline{B}_{s}$ meson mixing
oscillation is given by two theoretical constraints, which are achieved by
either lattice simulations DiLuzio:2019jyq  or lattice simulations plus sum
rule results DiLuzio:2019jyq

$\begin{split}\Delta
M_{s}^{\mathop{\mathrm{F}LAG}^{\prime}19}&=\left(1.13_{-0.07}^{+0.07}\right)\Delta
M_{s}^{\mathop{\mathrm{e}xp}},\\\ \Delta
M_{s}^{\mathop{\mathrm{A}verage}^{\prime}19}&=\left(1.04_{-0.09}^{+0.04}\right)\Delta
M_{s}^{\mathop{\mathrm{e}xp}},\end{split}$ (55)

The above one of Equation 55 is analyzed by the lattice simulations and gives
the strong bound Navarro:2021sfb

$G_{bs}<\frac{1}{\left(330\mathop{\mathrm{T}eV}\right)^{2}},$ (56)

and the below one is analyzed by both the lattice simulations and sum rule
results and gives the less constraining bound Navarro:2021sfb

$G_{bs}<\frac{1}{\left(220\mathop{\mathrm{T}eV}\right)^{2}},$ (57)

where it is worth mentioning that the less constraining theoretical constraint
is more consistent with the experimental bound. Instead of taking the
theoretical constraints into account, it is also possible to evaluate the
experimental bound in person. The experimental and SM bound for the $B_{s}$
meson oscillation reads Dowdall:2019bea :

$\begin{split}\Delta
M_{s}^{\mathop{\mathrm{e}xp}}&=17.757(21)\mathop{\mathrm{p}s}^{-1}\\\ \Delta
M_{s}^{\mathop{\mathrm{S}M}}&=17.59\begin{pmatrix}+0.33\\\
-1.22\end{pmatrix}\left(0.78\right)\mathop{\mathrm{p}s}^{-1}.\end{split}$ (58)

Then, the new physics effect for the $B_{s}$ meson oscillation can be
evaluated by subtracting the SM bound from the experimental bound as follows:

$\begin{gathered}\Delta M_{s}^{\mathop{\mathrm{N}P}}=\Delta
M_{s}^{\mathop{\mathrm{e}xp}}-\Delta
M_{s}^{\mathop{\mathrm{S}M}}=0.167\begin{pmatrix}+1.22\\\
-0.33\end{pmatrix}\left(-0.78\right)\mathop{\mathrm{p}s}^{-1},\\\
-1.07237\times 10^{-13}\mathop{\mathrm{G}eV}<\Delta
M_{s}^{\mathop{\mathrm{N}P}}<9.125\times
10^{-13}\mathop{\mathrm{G}eV}\end{gathered}$ (59)

We will constrain $Z^{\prime}$ mass range from all the three constraints once
scanned numerical ranges of the coupling constants $g_{bs}^{L,R}$ are obtained
in the quark sector simulation. The final task in this subsection is to
rewrite the new physics contribution in terms of the physical input parameters
Dowdall:2019bea ; Branco:2021vhs ; CarcamoHernandez:2020pnh :

$f_{B_{s}}=0.225\mathop{\mathrm{G}eV},\quad
B_{B_{s}}=1.26,\quad\eta_{B_{s}}=0.55,\quad
m_{B_{s}}=5.3663\mathop{\mathrm{G}eV}.$ (60)

Then, the new physics contribution to the $B_{s}$ meson oscillation reads in
terms of the physical input parameters:

$\begin{split}\Delta
M_{s}^{\mathop{\mathrm{N}P}}&=G_{bs}f_{B_{s}}^{2}B_{B_{s}}\eta_{B_{s}}m_{B_{s}}\\\
&=\frac{g_{bs}^{L2}+g_{bs}^{L}g_{bs}^{R*}+g_{bs}^{R2}}{M_{Z^{\prime}}^{2}}\times\left(0.225\mathop{\mathrm{G}eV}\right)^{2}\times
1.26\times 0.55\times\left(5.3663\mathop{\mathrm{G}eV}\right)\\\
&\simeq\frac{g_{bs}^{L2}+g_{bs}^{L}g_{bs}^{R*}+g_{bs}^{R2}}{M_{Z^{\prime}}^{2}}\times
0.188\mathop{\mathrm{G}eV}^{3}\end{split}$ (61)

###  VI.3 Collider experimental constraints

As the model under consideration features the massive second and third
generation of the SM fermions, the electron collider constraints arising from
the LEP experiments can not be applied to this BSM model. The most relevant
collider constrains arise from LHC experiments. Regarding light $Z^{\prime}$
masses, it can be constrained by the decay process $pp\rightarrow Z\rightarrow
4\mu$ and give rise to the constraint:

$5\mathop{\mathrm{G}eV}<M_{Z^{\prime}}<70\mathop{\mathrm{G}eV},$ (62)

where it is worth mentioning that we already reflected this condition when we
carried out the numerical scan for the charged lepton sector, thus providing
the condition on $Z^{\prime}$ mass which is
$M_{Z^{\prime}}>75\mathop{\mathrm{G}eV}$. The $Z^{\prime}$ mass can be further
constrained by LHC dimuon resonance measurements Abdullah:2017oqj ;
Alonso:2017uky , which are $pp\rightarrow
Z^{\prime}\rightarrow\overline{\mu}\mu$. Following the argument given in
Navarro:2021sfb , the main portal for generating $Z^{\prime}$ is
$\overline{b}b$ channel and this feature can be applied to our BSM model.
Using the narrow width approximation, the cross section for the decay process
$pp\rightarrow Z^{\prime}\rightarrow\overline{\mu}\mu$ is given for any
$g_{bb}$ Navarro:2021sfb :

$\begin{split}\sigma\left(pp\rightarrow
Z^{\prime}\rightarrow\overline{\mu}\mu\right)&\approx\sigma\left(pp\rightarrow
Z^{\prime}\right)\mathop{\mathrm{B}R}\left(Z^{\prime}\rightarrow\overline{\mu}\mu\right),\\\
\mathop{\mathrm{B}R}\left(Z^{\prime}\rightarrow\overline{\mu}\mu\right)&=\frac{\Gamma_{Z^{\prime}\rightarrow\overline{\mu}\mu}}{\Gamma_{Z^{\prime}\rightarrow\overline{\mu}\mu}+\Gamma_{Z^{\prime}\rightarrow\overline{\nu}_{\mu}\nu_{\mu}}+\Gamma_{Z^{\prime}\rightarrow\overline{b}b}+\Gamma_{Z^{\prime}\rightarrow\overline{s}b}+\Gamma_{Z^{\prime}\rightarrow\overline{t}t}},\end{split}$
(63)

where Navarro:2021sfb

$\begin{split}\Gamma\left(Z^{\prime}\rightarrow\overline{\mu}\mu\right)&=\frac{M_{Z^{\prime}}}{24\pi}\left[(g_{\mu\mu}^{L})^{2}+(g_{\mu\mu}^{R})^{2}\right],\\\
\Gamma\left(Z^{\prime}\rightarrow\overline{\nu}_{\mu}\nu_{\mu}\right)&=\frac{M_{Z^{\prime}}}{24\pi}\left[(g_{\nu_{\mu}\nu_{\mu}}^{L})^{2}\right],\\\
\Gamma\left(Z^{\prime}\rightarrow\overline{b}b\right)&=\frac{M_{Z^{\prime}}}{8\pi}\left[(g_{bb}^{L})^{2}+(g_{bb}^{R})^{2}\right],\\\
\Gamma\left(Z^{\prime}\rightarrow\overline{s}b\right)&=\frac{M_{Z^{\prime}}}{8\pi}\left[(g_{sb}^{L})^{2}+(g_{sb}^{R})^{2}\right],\\\
\Gamma\left(Z^{\prime}\rightarrow\overline{t}t\right)&=\frac{M_{Z^{\prime}}}{8\pi}\left[(g_{bb}^{L})^{2}+(g_{bb}^{R})^{2}\right]\left(1-\frac{m_{t}^{2}}{M_{Z^{\prime}}^{2}}\right)\sqrt{1-\frac{4m_{t}^{2}}{M_{Z^{\prime}}^{2}}},\end{split}$
(64)

where the $Z^{\prime}$ coupling constant $g_{bs}$ can be approximated to
$g_{bb}V_{ts}$ and the CKM mixing element $V_{tb}\approx 0.04$ play a role of
suppressing the $Z^{\prime}\overline{s}b$ interactions, therefore the main
decay channel is $Z^{\prime}\overline{b}b$. Given that the dimuon resonance
data from ATLAS ATLAS:2017fih  gives the experimental bound for $Z^{\prime}$
mass, the whole experimental bound reads Navarro:2021sfb :

$150\mathop{\mathrm{G}eV}<M_{Z^{\prime}}<5\mathop{\mathrm{T}eV}.$ (65)

###  VI.4 CKM mixing matrix

The CKM mixing matrix plays a quite crucial role in constraining the quark
sector in the model under consideration. As the investigation for the CKM
mixing matrix was done in one of our works CarcamoHernandez:2021yev  in
detail, we follow the same logic covered there. The renormalizable Lagrangian
for the CKM mixing matrix with the extended fermion spectrum by the vectorlike
family reads CarcamoHernandez:2021yev :

$\begin{split}\mathcal{L}_{\mathop{\mathrm{B}SM}}^{W}&=\frac{g}{\sqrt{2}}\begin{pmatrix}\overline{u}_{1L}&\overline{u}_{2L}&\overline{u}_{3L}&\overline{u}_{4L}&\overline{\widetilde{u}}_{4L}\end{pmatrix}\gamma_{\mu}\begin{pmatrix}1&0&0&0&0\\\
0&1&0&0&0\\\ 0&0&1&0&0\\\ 0&0&0&1&0\\\
0&0&0&0&0\end{pmatrix}\begin{pmatrix}d_{1L}\\\ d_{2L}\\\ d_{3L}\\\ d_{4L}\\\
\widetilde{d}_{4L}\end{pmatrix}W^{\mu+},\\\
&=\frac{g}{\sqrt{2}}\begin{pmatrix}\overline{u}_{1L}&\overline{u}_{2L}&\overline{u}_{3L}&\overline{u}_{4L}&\overline{\widetilde{u}}_{4L}\end{pmatrix}V_{L}^{u\dagger}V_{L}^{u}\gamma_{\mu}\begin{pmatrix}1&0&0&0&0\\\
0&1&0&0&0\\\ 0&0&1&0&0\\\ 0&0&0&1&0\\\
0&0&0&0&0\end{pmatrix}V_{L}^{d\dagger}V_{L}^{d}\begin{pmatrix}d_{1L}\\\
d_{2L}\\\ d_{3L}\\\ d_{4L}\\\ \widetilde{d}_{4L}\end{pmatrix}W^{\mu+},\\\
&=\frac{g}{\sqrt{2}}\begin{pmatrix}\overline{u}_{L}&\overline{c}_{L}&\overline{t}_{L}&\overline{U}_{4L}&\overline{\widetilde{U}}_{4L}\end{pmatrix}\gamma_{\mu}V_{L}^{u}\begin{pmatrix}1&0&0&0&0\\\
0&1&0&0&0\\\ 0&0&1&0&0\\\ 0&0&0&1&0\\\
0&0&0&0&0\end{pmatrix}V_{L}^{d\dagger}\begin{pmatrix}d_{L}\\\ s_{L}\\\
b_{L}\\\ D_{4L}\\\ \widetilde{D}_{4L}\end{pmatrix}W^{\mu+},\end{split}$ (66)

where the zero in the matrix between the mixing matrix $V_{L}^{u}$ and
$V_{L}^{d}$ arises from the fact that the vectorlike singlet quarks
$\widetilde{U}_{4}$ and $\widetilde{D}_{4}$ do not interact with the SM $W$
gauge boson, therefore featuring non-unitarity for the CKM mixing matrix in
this BSM model. The non-unitarity $\Delta$ reads Branco:2021vhs ;
Belfatto:2021jhf  as follows:

$\Delta=1-V_{ud}^{2}-V_{us}^{2}-V_{ub}^{2}$ (67)

What we found via the numerical scan in the previous work
CarcamoHernandez:2021yev  is the non-unitarity $\Delta$ is quite small
($\Delta=0.00035$ CarcamoHernandez:2021yev ) compared to its experimental
upper-limit ($\Delta_{\mathop{\mathrm{e}xp}}=0.04$ Belfatto:2019swo ). The
most fitted CKM prediction $(\chi_{\mathop{\mathrm{C}KM}}^{2}=956.828)$
corresponding to the non-unitarity $0.00035$ has the form of
CarcamoHernandez:2021yev :

$V_{\mathop{\mathrm{C}KM}}^{\mathop{\mathrm{p}red}}=\begin{pmatrix}0.97409&0.22602&0.00799&-6.38471\times
10^{-6}&-0.00036\\\\[2.15277pt] 0.22615&-0.97372&-0.02697&3.80102\times
10^{-5}&0.00147\\\\[2.15277pt]
0.00166&0.02815&-0.99880&-0.00766&0.00874\\\\[2.15277pt]
0.00003&0.00019&-0.00916&0.99919&-0.01773\\\\[2.15277pt]
-0.00057&0.00112&0.03812&0.03539&-0.00096\end{pmatrix}.$ (68)

Restricting our attention up to the upper-left $3\times 3$ block of Equation
68, the partial block can be compared to the experimental CKM mixing matrix
featuring non-unitarity Branco:2021vhs ; ParticleDataGroup:2020ssz :

$\lvert K_{\mathop{\mathrm{C}KM}}\rvert=\begin{pmatrix}0.97370\pm
0.00014&0.22450\pm 0.00080&0.00382\pm 0.00024\\\\[6.45831pt] 0.22100\pm
0.00400&0.98700\pm 0.01100&0.04100\pm 0.00140\\\\[6.45831pt] 0.00800\pm
0.00030&0.03880\pm 0.00110&1.01300\pm 0.03000\end{pmatrix},$ (69)

and the $31,32,13,23$ element of Equation 68 can not be fitted to their
experimental bounds at $3\sigma$ by a small difference. However, considering
this CKM prediction is an effective approach to the CKM mixing matrix since we
only have 23 mixing angle in the up-quark sector, the prediction is a good
approximation. For the purpose of making the numerical analysis for the quark
sector which will be carried out soon economical, we reuse the conditions
imposed in the previous work CarcamoHernandez:2021yev  that we collect
benchmark points which satisfy $\chi_{\mathop{\mathrm{C}KM}}^{2}<980$ and the
predicted $m_{c,s,b}$ and $m_{t}$ should be put between $\left[1\pm
0.1\right]\times m_{c,s,b}^{\mathop{\mathrm{e}xp}}$ and $\left[1\pm
0.01\right]\times m_{t}^{\mathop{\mathrm{e}xp}}$ respectively.

###  VI.5 Numerical analysis in the quark sector

A main goal of this numerical scan for the quark sector is to constrain the
$Z^{\prime}$ coupling constant to $b$ and $s$ quarks, in order to determine a
numerical $Z^{\prime}$ mass range from the $R_{K^{*}}$ anomaly (Equation
LABEL:eqn:RK_const_theoretical_clean_fit and LABEL:eqn:RK_const_global_fit)
and from the $B_{s}$ meson oscillation experimental constraints (Equation 56,
57 and 59), and then we try to find an overlapped region among the numerical
$Z^{\prime}$ mass ranges. For the task, we reuse the best fitted benchmark
point ($\chi_{\mathop{\mathrm{C}KM}}^{2}=956.828$) given in our previous work
CarcamoHernandez:2021yev  instead of setting numerical scan up from scratch
since this BSM model can also gain access to the CKM mixing matrix in the same
way as done in CarcamoHernandez:2021yev

$\begin{gathered}M^{u}=\begin{pmatrix}0&0&0&0&0\\\\[2.15277pt]
0&0&0&0&14.474\\\\[2.15277pt] 0&0&0&1206.340&277.563\\\\[2.15277pt]
0&0&273.503&-1775.200&0\\\\[2.15277pt]
0&550.990&434.462&0&-5624.050\end{pmatrix}\quad
M^{d}=\begin{pmatrix}0&0&0&0&-0.938\\\\[2.15277pt]
0&0&0&0&-4.041\\\\[2.15277pt] 0&0&0&1206.340&-27.427\\\\[2.15277pt]
0&0&-5.636&-1775.200&0\\\\[2.15277pt]
0&72.915&-75.760&0&2623.620\end{pmatrix}\\\\[8.61108pt]
M_{\mathop{\mathrm{d}iag}}^{u}=\begin{pmatrix}0&0&0&0&0\\\\[2.15277pt]
0&1.255&0&0&0\\\\[2.15277pt] 0&0&171.303&0&0\\\\[2.15277pt]
0&0&0&2155.890&0\\\\[2.15277pt] 0&0&0&0&5674.840\end{pmatrix}\quad
M_{\mathop{\mathrm{d}iag}}^{d}=\begin{pmatrix}0&0&0&0&0\\\\[2.15277pt]
0&0.094&0&0&0\\\\[2.15277pt] 0&0&3.875&0&0\\\\[2.15277pt]
0&0&0&2146.190&0\\\\[2.15277pt] 0&0&0&0&2625.960\end{pmatrix},\end{gathered}$
(70)

which gives rise to the best fitted CKM prediction of Equation 68. We vary
mass parameters in the best fitted benchmark point of Equation 70 by a factor
of $\left[1\pm\kappa\right]$ where $\kappa=0.1$ in order to find the mixing
matrices $V_{L,R}^{u}$ and $V_{L,R}^{d}$, required to derive the $Z^{\prime}$
coupling constant to $b$ and $s$ quarks, and then collect the benchmark points
which satisfy $\chi_{\mathop{\mathrm{C}KM}}^{2}<980$ as done in
CarcamoHernandez:2021yev . The varying process is shown in Table 8:

Mass parameter | Scanned Region($\mathop{\mathrm{G}eV}$)
---|---
$y_{24}^{u}v_{u}=m_{24}^{u}$ | $\left[1\pm\kappa\right]\times m_{24r}^{u}$
$y_{34}^{u}v_{u}=m_{34}^{u}$ | $\left[1\pm\kappa\right]\times m_{34r}^{u}$
$y_{43}^{u}v_{u}=m_{43}^{u}$ | $\left[1\pm\kappa\right]\times m_{43r}^{u}$
$x_{34}^{Q}v_{\phi}=m_{35}^{Q}$ | $m_{35}^{Q}$
$x_{42}^{u}v_{\phi}=m_{52}^{u}$ | $\left[1\pm\kappa\right]\times m_{42r}^{u}$
$x_{43}^{u}v_{\phi}=m_{53}^{u}$ | $\left[1\pm\kappa\right]\times m_{43r}^{u}$
$M_{45}^{Q}$ | $M_{45}^{Q}$
$M_{54}^{u}$ | $\left[1\pm\kappa\right]\times M_{54r}^{u}$
$y_{14}^{d}v_{d}=m_{14}^{d}$ | $\left[1\pm\kappa\right]\times m_{14r}^{d}$
$y_{24}^{d}v_{d}=m_{24}^{d}$ | $\left[1\pm\kappa\right]\times m_{24r}^{d}$
$y_{34}^{d}v_{d}=m_{34}^{d}$ | $\left[1\pm\kappa\right]\times m_{34r}^{d}$
$y_{43}^{d}v_{d}=m_{43}^{d}$ | $\left[1\pm\kappa\right]\times m_{43r}^{d}$
$x_{34}^{Q}v_{\phi}=m_{35}^{Q}$ | $\left[1\pm\kappa\right]\times m_{35r}^{Q}$
$x_{42}^{d}v_{\phi}=m_{52}^{d}$ | $\left[1\pm\kappa\right]\times m_{42r}^{d}$
$x_{43}^{d}v_{\phi}=m_{53}^{d}$ | $\left[1\pm\kappa\right]\times m_{43r}^{d}$
$M_{45}^{Q}$ | $\left[1\pm\kappa\right]\times M_{45r}^{Q}$
$M_{54}^{d}$ | $\left[1\pm\kappa\right]\times M_{54r}^{d}$
$\kappa$ | $0.1$

Table 8: Parameter setup, based on the best fitted benchmark point of Equation
70, to find the mixing matrices $V_{L,R}^{u}$ and $V_{L,R}^{d}$ required to
derive the $Z^{\prime}$ coupling constant to $b$ and $s$ quarks. The imposed
conditions to fit the mass parameters to the SM are the predicted $c,s,b$ and
$t$ quark masses are put between $\left[1\pm 0.1\right]\times
m_{c,s,b}^{\mathop{\mathrm{e}xp}}$ and between $\left[1\pm 0.01\right]\times
m_{t}^{\mathop{\mathrm{e}xp}}$ and the CKM mixing matrix
CarcamoHernandez:2021yev .

####  VI.5.1 Numerical scan result and ranges of the $Z^{\prime}$ coupling
constant to $b$ and $s$ quarks

The scanned vectorlike quark masses are shown in Figure 13:

Figure 13: Scanned mass ranges of vectorlike quarks. $U_{4}$ and $D_{4}$ are
vectorlike doublet up- and down-type quarks respectively and
$\widetilde{U}_{4}$ and $\widetilde{D}_{4}$ are vectorlike singlet up- and
down-type quarks respectively. The two above plots are exactly consistent with
ones obtained in our previous work CarcamoHernandez:2021yev . The vectorlike
doublet quarks $M_{U_{4}}$ and $M_{D_{4}}$ shows strong correlation between
their masses, which can be directly seen from the fifth column of mass
matrices of Equation 18, and this correlation plays an important role in
suppressing fermion contributions to the oblique parameters $T,S,U$ as we will
see soon. Here, the $Z^{\prime}$ mass is still free parameter.

The upper-left plot in Figure 13 shows mass ranges for the vectorlike doublet
up-type quark $M_{U_{4}}$ against the vectorlike singlet up-type quark
$M_{\widetilde{U}_{4}}$. The upper-right plot shows mass ranges for the
vectorlike doublet down-type quark $M_{D_{4}}$ against the vectorlike singlet
down-type quark $M_{\widetilde{D}_{4}}$. The reported vectorlike doublet and
singlet quark masses in our quark sector simulation are free from the most
recent experimental vectorlike mass constraints at $95\%$ C.L. given by the
ATLAS experiments ATLAS:2021ibc :

$\begin{split}M_{U_{4}}&>1.46\mathop{\mathrm{T}eV},\\\
M_{\widetilde{U}_{4}}&>1.27\mathop{\mathrm{T}eV},\\\
M_{D_{4}}&>1.32\mathop{\mathrm{T}eV},\\\
M_{\widetilde{D}_{4}}&>1.20\mathop{\mathrm{T}eV}.\end{split}$ (71)

Considering 100$\%$ branching ratio for $U\rightarrow Zt$ and for
$D\rightarrow Zb$, the upper-limits for the vectorlike up- and down-type quark
masses increase up to $1.60\mathop{\mathrm{T}eV}$ and
$1.42\mathop{\mathrm{T}eV}$ at $95\%$ C.L. respectively ATLAS:2021ibc :

$\begin{split}M_{U_{4}},M_{\widetilde{U}_{4}}&>1.60\mathop{\mathrm{T}eV},\\\
M_{D_{4}},M_{\widetilde{D}_{4}}&>1.42\mathop{\mathrm{T}eV}.\\\ \end{split}$
(72)

All plots in Figure 13 show no correlation, displaying only the allowed range
of the vectorlike quark masses, excepting for the lower-left plot. The lower-
left plot shows that the vectorlike doublet up-type and down type quarks are
degenerate and this property can be directly seen from the quark mass matrices
of Equation 18. The equal masses for the vectorlike up-type and down type
quarks doublets are quite important to fit the experimental oblique $S,T,U$
parameters. In other words, if their mass difference is large, it will yield
maximal violation of the custodial symmetry thus implying unacceptably large
contributions to these oblique parameters. Thus, a large up and down type
vector like doublet quark mass splitting will not allow to successfully
accommodate the experimental bounds for the oblique parameters. Therefore, the
fact that the structure of the quark mass matrices can provide the degenerate
masses for the vectorlike doublet quarks in an analytic way makes this BSM
theory more convincible. Going back to scanning the $Z^{\prime}$ coupling
constant to $b$ and $s$ quarks, they are given by (we also included the
derived $Z^{\prime}$ coupling constant to muon pair together for convenience):

Coupling constant | Scanned range
---|---
$g_{\mu\mu}^{L}$ | $4.888\times 10^{-6}<g_{\mu\mu}^{L}<2.113\times 10^{-3}$
$g_{\mu\mu}^{R}$ | $-1.195\times 10^{-2}<g_{\mu\mu}^{R}<-3.756\times 10^{-4}$
$g_{bs}^{L}$ | $-1.439\times 10^{-2}<g_{bs}^{L}<1.404\times 10^{-2}$
$g_{bs}^{R}$ | $-9.773\times 10^{-4}<g_{bs}^{R}<9.723\times 10^{-4}$
$g_{\mu\mu}^{L}+g_{\mu\mu}^{R}$ | $-1.175\times 10^{-2}<g_{\mu\mu}^{L}+g_{\mu\mu}^{R}<1.179\times 10^{-3}$

Table 9: Scanned ranges of the $Z^{\prime}$ coupling constant to $b$ and $s$
quarks as well as to muon pair, required to derive mass range of $Z^{\prime}$
gauge boson from the $R_{K^{*}}$ anomaly and from the $B_{s}$ meson
oscillation.

####  VI.5.2 $M_{Z^{\prime}}$ mass ranges resulting from the $R_{K^{*}}$
anomaly and $B_{s}$ meson oscillation

In the model under consideration, the most relevant mass ranges of
$Z^{\prime}$ come from the $R_{K^{*}}$ anomaly and the $B_{s}$ meson
oscillation. We start from the $R_{K^{*}}$ anomaly first.

$R_{K^{*}}$ anomaly

Reminding the experimental bounds for the Wilson coefficients
$G_{bs\mu}^{L,R}$ in both “theoretically clean fit” (Equation
LABEL:eqn:RK_const_tcfit) and “global fit” (Equation LABEL:eqn:RK_const_gfit),
they read respectively:

$\begin{split}\frac{1}{\left(44.44\mathop{\mathrm{T}eV}\right)^{2}}&<G_{bs\mu}^{L}=-\frac{g_{bs}^{L}g_{\mu\mu}^{L}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(40\mathop{\mathrm{T}eV}\right)^{2}},\\\
-\frac{1}{\left(35.9\mathop{\mathrm{T}eV}\right)^{2}}&<G_{bs\mu}^{R}=-\frac{g_{bs}^{L}g_{\mu\mu}^{R}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(205\mathop{\mathrm{T}eV}\right)^{2}},\\\
\end{split}$ (73)
$\begin{split}\frac{1}{\left(39.75\mathop{\mathrm{T}eV}\right)^{2}}&<G_{bs\mu}^{L}=-\frac{g_{bs}^{L}g_{\mu\mu}^{L}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(36.67\mathop{\mathrm{T}eV}\right)^{2}},\\\
-\frac{1}{\left(83.8\mathop{\mathrm{T}eV}\right)^{2}}&<G_{bs\mu}^{R}=-\frac{g_{bs}^{L}g_{\mu\mu}^{R}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(44.44\mathop{\mathrm{T}eV}\right)^{2}}.\\\
\end{split}$ (74)

In order to confirm $M_{Z^{\prime}}$ ranges, the first task is to determine
numerical ranges of the combined coupling constants $g_{bs}g_{\mu\mu}$ using
the varied range of each coupling constant in Table 9 as given in Table 10:

Combination of the coupling constants | Range of the combination
---|---
$g_{bs}^{L}g_{\mu\mu}^{L}$ | $-3.041\times 10^{-5}<g_{bs}^{L}g_{\mu\mu}^{L}<2.967\times 10^{-5}$
$g_{bs}^{L}g_{\mu\mu}^{R}$ | $-1.678\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{R}<1.720\times 10^{-4}$

Table 10: Combined coupling constants, required to confirm mass ranges of
$Z^{\prime}$ from the $R_{K^{*}}$ anomaly.

Using each boundary value of the combined coupling constants
$g_{bs}g_{\mu\mu}$, we can confirm mass ranges of $Z^{\prime}$ when both the
CMS and $\tau\rightarrow\mu\gamma$ experimental constraints are considered
together as given in Table 11:

Fit | Combination of the coupling constants | Mass range of $Z^{\prime}$ ($\mathop{\mathrm{G}eV}$)
---|---|---
Theoretically clean fit | $g_{bs}^{L}g_{\mu\mu}^{L}=-3.041\times 10^{-5}$ | $220.581<M_{Z^{\prime}}<245.066$
$g_{bs}^{L}g_{\mu\mu}^{L}=2.967\times 10^{-5}$ | no result
(CMS + $\left(\tau\rightarrow\mu\gamma\right)$) | $g_{bs}^{L}g_{\mu\mu}^{R}=-1.678\times 10^{-4}$ | $M_{Z^{\prime}}>2655.522$
| $g_{bs}^{L}g_{\mu\mu}^{R}=1.720\times 10^{-4}$ | $M_{Z^{\prime}}>470.824$
Global fit | $g_{bs}^{L}g_{\mu\mu}^{L}=-3.041\times 10^{-5}$ | $202.218<M_{Z^{\prime}}<219.202$
$g_{bs}^{L}g_{\mu\mu}^{L}=2.967\times 10^{-5}$ | no result
(CMS + $\left(\tau\rightarrow\mu\gamma\right)$) | $g_{bs}^{L}g_{\mu\mu}^{R}=-1.678\times 10^{-4}$ | $M_{Z^{\prime}}>575.665$
| $g_{bs}^{L}g_{\mu\mu}^{R}=1.720\times 10^{-4}$ | $M_{Z^{\prime}}>1099.027$
Table 11: Scanned ranges of $M_{Z^{\prime}}$, based on the combined coupling
constants, from the theoretically clean and global fit in the $R_{K^{*}}$
anomaly. The parenthesis below each fit means the constraints considered. “No
result” means there is no possible range due to overall sign difference.

Simplifying the given results in Table 11 by contracting each combined
coupling constants, it can be rearranged as follows:

Fit | Combination of the coupling constants | Mass range of $Z^{\prime}$ ($\mathop{\mathrm{G}eV}$)
---|---|---
Theoretically clean fit | $-3.041\times 10^{-5}<g_{bs}^{L}g_{\mu\mu}^{L}<2.967\times 10^{-5}$ | $220.581<M_{Z^{\prime}}<245.066$
(CMS + $\left(\tau\rightarrow\mu\gamma\right)$) | $-1.678\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{R}<1.720\times 10^{-4}$ | $M_{Z^{\prime}}>2655.522$
Theoretically clean fit | $-1.927\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{L}<1.880\times 10^{-4}$ | $555.266<M_{Z^{\prime}}<616.900$
$\left(\tau\rightarrow\mu\gamma\right)$ | $-1.726\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{R}<1.769\times 10^{-4}$ | $M_{Z^{\prime}}>2693.235$
Global fit | $-3.041\times 10^{-5}<g_{bs}^{L}g_{\mu\mu}^{L}<2.967\times 10^{-5}$ | $202.218<M_{Z^{\prime}}<219.202$
(CMS + $\left(\tau\rightarrow\mu\gamma\right)$) | $-1.678\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{R}<1.720\times 10^{-4}$ | $M_{Z^{\prime}}>1099.027$
Global fit | $-1.927\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{L}<1.880\times 10^{-4}$ | $509.040<M_{Z^{\prime}}<551.795$
$\left(\tau\rightarrow\mu\gamma\right)$ | $-1.726\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{R}<1.769\times 10^{-4}$ | $M_{Z^{\prime}}>1114.570$
Table 12: Ranges of $M_{Z^{\prime}}$ simplified. Apart from the given results
in Table 11, we also include $Z^{\prime}$ mass ranges derived from the
$R_{K^{*}}$ anomaly when the CLFV $\tau\rightarrow\mu\gamma$ experimental
constraint is considered alone to show that how the derived $Z^{\prime}$ mass
ranges get affected and that the $Z^{\prime}$ and CP-odd scalar can be
connected via the scalar potential in this case only.

What we found via this numerical study for the $R_{K^{*}}$ anomaly is three
features: the first is there is no overlapped region between $M_{Z^{\prime}}$
ranges from $G_{bs\mu}^{L}$ and those from $G_{bs\mu}^{R}$ in both
theoretically clean and global fit, which indicates that the $R_{K^{*}}$
anomaly is quite likely to be explained by only left-handed interactions, and
the second is there is no overlapped region for $G_{bs\mu}^{L}$ even in
theoretically clean fit and in global fit. Lastly, we include the derived
$Z^{\prime}$ mass ranges from the $R_{K^{*}}$ anomaly when the CLFV
$\tau\rightarrow\mu\gamma$ constraint is considered alone to show that how the
$Z^{\prime}$ mass ranges get affected and that the $Z^{\prime}$ and the CP-odd
scalar can be connected in this case only via the scalar potential. Next we
discuss the $B_{s}$ meson oscillation.

$B_{s}$ meson oscillation

Reminding the two theoretical constraints of Equation 56, 57 and the
experimental constraint of Equation 59, they read respectively:

$\begin{gathered}G_{bs}=\frac{g_{bs}^{L2}+g_{bs}^{L}g_{bs}^{R*}+g_{bs}^{R2}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(330\mathop{\mathrm{T}eV}\right)^{2}},\quad\text{from
FLAG${}^{\prime}$19}\\\
G_{bs}=\frac{g_{bs}^{L2}+g_{bs}^{L}g_{bs}^{R*}+g_{bs}^{R2}}{M_{Z^{\prime}}^{2}}<\frac{1}{\left(220\mathop{\mathrm{T}eV}\right)^{2}},\quad\text{from
Average${}^{\prime}$19}\\\ -1.07237\times
10^{-13}\mathop{\mathrm{G}eV}<\frac{g_{bs}^{L2}+g_{bs}^{L}g_{bs}^{R*}+g_{bs}^{R2}}{M_{Z^{\prime}}^{2}}\times
0.188\mathop{\mathrm{G}eV}^{3}<9.125\times
10^{-13}\mathop{\mathrm{G}eV}.\end{gathered}$ (75)

A numerical range of summing over the flavor violating coupling constants
$g_{bs}^{L,R}$, referring to the coupling constants in Table 9, reads:

$1.833\times
10^{-4}<g_{bs}^{2}=g_{bs}^{L2}+g_{bs}^{L}g_{bs}^{R*}+g_{bs}^{R2}<2.106\times
10^{-4}$ (76)

Then, $M_{Z^{\prime}}$ can be determined as in the $R_{K^{*}}$ anomaly and the
final result for the $B_{s}$ meson oscillation is given in Table 13:

Fit | Combination of the coupling constants | Mass range of $Z^{\prime}$ ($\mathop{\mathrm{G}eV}$)
---|---|---
Theoretical fit (FLAG) | $1.833\times 10^{-4}<g_{bs}^{2}=g_{bs}^{L2}+g_{bs}^{L}g_{bs}^{R*}+g_{bs}^{R2}<2.106\times 10^{-4}$ | $M_{Z^{\prime}}>4788.980$
Theoretical fit (Average) | $M_{Z^{\prime}}>3192.650$
Experimental result | $M_{Z^{\prime}}>6591.740$
Table 13: Ranges of $M_{Z^{\prime}}$ derived from the $B_{s}$ meson
oscillation experiments. They are derived from two theoretical fits, FLAG and
Average, and from one experimental result respectively.

The confirmed $Z^{\prime}$ mass ranges from the $B_{s}$ meson oscillation fits
are too much tightened when compared to those derived from the $R_{K^{*}}$
anomaly as well as from the neutrino trident production. For comparison, it is
good to look at the derived results from the diverse constraints at one sight
in Table 14:

Observable | Fit | Coupling constants | Mass range of $Z^{\prime}$ ($\mathop{\mathrm{G}eV}$)
---|---|---|---
$\nu$ Trident |  |  | $M_{Z^{\prime}}>0.920\mathop{\mathrm{G}eV}$
$R_{K^{*}}$ anomaly | Theoretically clean fit | $-3.041\times 10^{-5}<g_{bs}^{L}g_{\mu\mu}^{L}<2.967\times 10^{-5}$ | $220.581<M_{Z^{\prime}}<245.066$
(CMS + $\left(\tau\rightarrow\mu\gamma\right)$) | $-1.678\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{R}<1.720\times 10^{-4}$ | $M_{Z^{\prime}}>2655.522$
Theoretically clean fit | $-1.927\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{L}<1.880\times 10^{-4}$ | $555.266<M_{Z^{\prime}}<616.900$
$\left(\tau\rightarrow\mu\gamma\right)$ | $-1.726\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{R}<1.769\times 10^{-4}$ | $M_{Z^{\prime}}>2693.235$
Global fit | $-3.041\times 10^{-5}<g_{bs}^{L}g_{\mu\mu}^{L}<2.967\times 10^{-5}$ | $202.218<M_{Z^{\prime}}<219.202$
(CMS + $\left(\tau\rightarrow\mu\gamma\right)$) | $-1.678\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{R}<1.720\times 10^{-4}$ | $M_{Z^{\prime}}>1099.027$
Global fit | $-1.927\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{L}<1.880\times 10^{-4}$ | $509.040<M_{Z^{\prime}}<551.795$
$\left(\tau\rightarrow\mu\gamma\right)$ | $-1.726\times 10^{-4}<g_{bs}^{L}g_{\mu\mu}^{R}<1.769\times 10^{-4}$ | $M_{Z^{\prime}}>1114.570$
$B_{s}-\overline{B}_{s}$ | Theoretical fit (FLAG) | $1.833\times 10^{-4}<g_{bs}^{2}<2.106\times 10^{-4}$ | $M_{Z^{\prime}}>3192.650$
Theoretical fit (Average) | $M_{Z^{\prime}}>4788.980$
Experimental result | $M_{Z^{\prime}}>6591.740$
Table 14: Ranges of $M_{Z^{\prime}}$ derived from the neutrino trident
production ($\nu$ Trident), the $R_{K^{*}}$ anomaly and finally the $B_{s}$
meson oscillation.

The first constraint we discuss in Table 14 is the neutrino trident production
and it can not constrain $Z^{\prime}$ mass compared to ones from the
$R_{K^{*}}$ anomaly and $B_{s}$ meson oscillation. Next we discuss constraints
resulting from the $R_{K^{*}}$ anomaly and we consider two cases where the
first is to consider both the CMS and $\tau\rightarrow\mu\gamma$ experimental
constraints and the other is to consider only the $\tau\rightarrow\mu\gamma$
experimental constraint. The reasons why we consider both the cases are as
follows:

  1. 1.

The first reason is to show that how the $Z^{\prime}$ mass range is affected
when the left-handed muon coupling constant $g_{L}^{\mu\mu}$ is suppressed by
the CMS experimental limit.

  2. 2.

The second reason is the $R_{K^{*}}$ anomaly global fit with the experimental
CLFV $\tau\rightarrow\mu\gamma$ decay alone is the only possible region where
non-SM scalars for the muon $g-2$ and the $Z^{\prime}$ gauge boson for the
$R_{K^{*}}$ anomaly can be connected via the scalar potential under
consideration, while fitting the oblique parameters $T,S,U$ as well as the $W$
mass anomaly and not changing out conclusion that the muon $g-2$ and
$R_{K^{*}}$ anomaly can not be explained by the same new physics. So we call
the region “theoretically interesting $Z^{\prime}$ mass range”.

Taking the reasons listed into account, we consider the $R_{K^{*}}$ anomaly
constraint with both the CMS and $\tau\rightarrow\mu\gamma$ experiments first.
The first important feature is there is no overlapped region between the two
combination of coupling constants in both the theoretically clean fit and
global fit and what this implements is the $R_{K^{*}}$ anomaly is likely to be
explained by the only left-handed interactions. This tendency is not changed
for the constraints with the $\tau\rightarrow\mu\gamma$ experimental limit
alone. The second feature is the $Z^{\prime}$ mass range from the $R_{K^{*}}$
anomaly global fit in Table 14 with both the experimental bounds, CMS and
$\tau\rightarrow\mu\gamma$, is not compatible to the lightest $Z^{\prime}$
mass, $261\mathop{\mathrm{G}eV}$, derived in the squared gauge mass matrix of
Equation 106 under this BSM model when we turn off the vev $v_{3}$ (even
though we consider the $R_{K^{*}}$ anomaly constraint at $3\sigma$, giving
rise to
$186.942\mathop{\mathrm{G}eV}<M_{Z^{\prime}}<237.786\mathop{\mathrm{G}eV}$,
its tension is not relaxed enough to cover the lightest $Z^{\prime}$ mass). On
top of that, the vev $v_{3}$ is somewhat constrained by both the CKM mixing
matrix and the muon $g-2$. Looking at the up-type mass matrix of Equation 70,
giving rise to the best fitted CKM prediction of Equation 68 in the model
under consideration, one can know that the vev $v_{3}$ must be heavier than
$340\mathop{\mathrm{G}eV}$, which is given by diving the
$m_{35}=1206.340\mathop{\mathrm{G}eV}$ by the purturbative limit of a Yukawa
constant $\sqrt{4\pi}=3.54$. Plus, the vev $v_{3}$ is also constrained by the
muon $g-2$ prediction of Equation 101. Referring to our previous work
Hernandez:2021tii , one of the Yukawa constant, $y_{24}^{e}$, in the muon
prediction with scalar exchange remains as a free parameter whereas the other,
$x_{42}^{e}$, is governed by the vev $v_{3}$ locating in the denominator.
Therefore, the vev $v_{3}$ can not be heavy, for example, up to
$\mathop{\mathrm{T}eV}$ scale in the model under consideration, otherwise the
muon $g-2$ prediction will have a smaller order compared to its experimental
order, therefore not possible to fit its experimental bound at all. What this
implements is the $R_{K^{*}}$ anomaly and muon $g-2$ can not be explained by
the same new physics and their new physics sources, $Z^{\prime}$ for the
$R_{K^{*}}$ anomaly and non-SM scalars for the muon $g-2$ can not be connected
to each other in the scalar potential under consideration, therefore
concluding the $R_{K^{*}}$ anomaly and muon $g-2$ should be considered as an
independent new physics signal in this BSM model. Considering the
“theoretically interesting $Z^{\prime}$ mass range” is less interesting for
the reason that the CMS experimental bound is not considered, however we will
show that the $Z^{\prime}$ for the $R_{K^{*}}$ anomaly and the non-SM scalars
for the muon $g-2$ can be connected via the scalar potential while fitting
oblique parameters $T,S,U$ and $W$ mass anomaly at their $1$ or $2\sigma$
constraints in the range. Lastly, since the $Z^{\prime}$ mass ranges derived
from the $B_{s}$ meson oscillations are too heavy when compared to ones from
the neutrino trident and $R_{K^{*}}$ anomaly, we consider them to be explained
by some other heavier new physics sources, not considered in this work.

##  VII Muon anomalous magnetic moment $g-2$ revisited with scalar exchange

We have seen that the muon $g-2$ mediated by $Z^{\prime}$ at one loop level is
not an answer since it can not explain the experimental muon $g-2$ bound at
all by yielding only negative contributions. In order to cure this problem, it
requires another approach to muon $g-2$ and we have shown that scalar exchange
at one loop level can be an answer for muon $g-2$ Hernandez:2021tii . The
corresponding muon $g-2$ diagram with scalar exchange at one loop level can be
drawn by closing the scalars of the diagrams in Figure 3:

Figure 14: Feynman diagram contributing to the muon $g-2$ mediated by the
vectorlike singlet charged leptons and by the CP-even scalars $H_{1,2}$ and
the CP-odd scalar $A_{1}$,

where it is worth mentioning that the muon $g-2$ mediated by the vectorlike
doublet charged leptons does not work in this BSM model as there is no seesaw
operator mediated by the vectorlike doublet charged lepton in the fully
rotated mass matrix of Equation 12 for the purpose of imposing the SM
hierarchy dynamically. In order to setup an analytic prediction for the muon
$g-2$ by scalar exchange, it requires to discuss the scalar potential and its
relevant features and this process is interconnected with the $M_{Z^{\prime}}$
as we will see soon.

###  VII.1 Scalar potential

In this BSM model, the Higgs alignments in the interaction basis read
Hernandez:2021tii  ($v_{u}=v_{1},v_{d}=v_{2},v_{\phi}=v_{3}$):

$\begin{split}H_{u}&=\begin{pmatrix}H_{u}^{+}\\\
v_{u}+\frac{1}{\sqrt{2}}\left(\mathop{\mathrm{R}e}H_{u}^{0}+i\mathop{\mathrm{I}m}H_{u}^{0}\right)\end{pmatrix},\\\
H_{d}&=\begin{pmatrix}v_{d}+\frac{1}{\sqrt{2}}\left(\mathop{\mathrm{R}e}H_{d}^{0}+i\mathop{\mathrm{I}m}H_{d}^{0}\right)\\\
,H_{d}^{-}\end{pmatrix},\\\
\phi&=\frac{1}{\sqrt{2}}\left(v_{\phi}+\mathop{\mathrm{R}e}\phi+i\mathop{\mathrm{I}m}\phi\right).\end{split}$
(77)

The scalar potential takes the form of Hernandez:2021tii :

$\begin{split}V&=\mu_{1}^{2}(H_{u}^{\dagger}H_{u})+\mu_{2}^{2}(H_{d}^{\dagger}H_{d})+\mu_{3}^{2}(\phi^{*}\phi)\\\
&+\lambda_{1}(H_{u}^{\dagger}H_{u})^{2}+\lambda_{2}(H_{d}^{\dagger}H_{d})^{2}+\lambda_{3}(H_{u}^{\dagger}H_{u})(H_{d}^{\dagger}H_{d})+\lambda_{4}(H_{u}^{\dagger}H_{d})(H_{d}^{\dagger}H_{u})\\\
&+\lambda_{5}(\epsilon_{ij}H_{u}^{i}H_{d}^{j}\phi^{2}+\mathop{\mathrm{H}.c.})+\lambda_{6}(\phi^{*}\phi)^{2}+\lambda_{7}(\phi^{*}\phi)(H_{u}^{\dagger}H_{u})+\lambda_{8}(\phi^{*}\phi)(H_{d}^{\dagger}H_{d}),\end{split}$
(78)

where the $\mu_{1,2,3}^{2}$ are dimensionful mass parameters and the
$\lambda_{1,2,\cdots,8}$ are dimensionless quartic coupling constants. The
minimization conditions of the scalar potential allow the dimensionful mass
parameters to be written in terms of the input parameters of the scalar
potential as follows:

$\begin{split}\mu_{1}^{2}&=-2\lambda_{1}v_{1}^{2}-\lambda_{3}v_{2}^{2}+\frac{\lambda_{5}v_{2}v_{3}^{2}}{2v_{1}}-\lambda_{7}v_{3}^{2},\\\
\mu_{2}^{2}&=-2\lambda_{2}v_{2}^{2}-\lambda_{3}v_{1}^{2}+\frac{\lambda_{5}v_{1}v_{3}^{2}}{2v_{2}}-\lambda_{8}v_{3}^{2},\\\
\mu_{3}^{2}&=-2\lambda_{6}v_{3}^{2}-\lambda_{7}v_{1}^{2}-\lambda_{8}v_{2}^{2}+2\lambda_{5}v_{1}v_{2}.\end{split}$
(79)

The dimensionful mass parameters are important in exploring vacuum stability
of this BSM model and will be discussed in the next subsection.

###  VII.2 Mass matrix for CP-even, CP-odd neutral and charged scalars

The squared CP-even mass matrix in the interaction basis of
$\left(\mathop{\mathrm{R}e}H_{u}^{0},\mathop{\mathrm{R}e}H_{d}^{0},\mathop{\mathrm{R}e}\phi\right)$
reads:

$\mathbf{M}_{\mathop{\mathrm{C}P-even}}^{2}=\left(\begin{array}[]{ccc}4\lambda_{1}v_{1}^{2}+\frac{v_{2}v_{3}^{2}\lambda_{5}}{2v_{1}}&2v_{1}v_{2}\lambda_{3}-\frac{1}{2}v_{3}^{2}\lambda_{5}&\sqrt{2}v_{3}\left(v_{1}\lambda_{7}-v_{2}\lambda_{5}\right)\\\
2v_{1}v_{2}\lambda_{3}-\frac{1}{2}v_{3}^{2}\lambda_{5}&4\lambda_{2}v_{2}^{2}+\frac{v_{1}v_{3}^{2}\lambda_{5}}{2v_{2}}&\sqrt{2}v_{3}\left(v_{2}\lambda_{8}-v_{1}\lambda_{5}\right)\\\
\sqrt{2}v_{3}\left(v_{1}\lambda_{7}-v_{2}\lambda_{5}\right)&\sqrt{2}v_{3}\left(v_{2}\lambda_{8}-v_{1}\lambda_{5}\right)&2v_{3}^{2}\lambda_{6}\\\
\end{array}\right).$ (80)

In order to make our analysis for the scalar sector simpler, we apply
decoupling limit to the scalar sector in this BSM model. Under the decoupling
limit, the SM Higgs $h$ has no mixing with CP-even scalars $H_{1}$ and $H_{2}$
and it can be achieved by turning off $(12)$, $(13)$, $(21)$ and $(31)$
elements of the CP-even mass matrix of Equation 80. The decoupling limit can
be implemented by setting a couple of relations between the relevant quartic
coupling constants as follows:

$\begin{split}\lambda_{5}&=\frac{4v_{1}v_{2}}{v_{3}^{2}}\lambda_{3},\\\
\lambda_{7}&=\frac{v_{2}}{v_{1}}\lambda_{5}=\frac{4v_{2}^{2}}{v_{3}^{2}}\lambda_{3}.\end{split}$
(81)

Then, the squared CP-even mass matrix after implementing the decoupling limit
is reduced to:

$\textbf{M}_{\mathop{\mathrm{C}P-even}}^{2}=\left(\begin{array}[]{ccc}4\lambda_{1}v_{1}^{2}+2v_{2}^{2}\lambda_{3}&0&0\\\
0&2\lambda_{3}v_{1}^{2}+4v_{2}^{2}\lambda_{2}&\sqrt{2}v_{3}\left(v_{2}\lambda_{8}-\frac{4v_{1}^{2}v_{2}\lambda_{3}}{v_{3}^{2}}\right)\\\
0&\sqrt{2}v_{3}\left(v_{2}\lambda_{8}-\frac{4v_{1}^{2}v_{2}\lambda_{3}}{v_{3}^{2}}\right)&2v_{3}^{2}\lambda_{6}\\\
\end{array}\right).$ (82)

Diagonalizing the squared CP-even mass matrix, it reveals masses of the
physical SM Higgs and non-SM scalars $H_{1,2}$ in the physical basis of
$\left(h,H_{1},H_{2}\right)$:

$R_{\mathop{\mathrm{C}P-even}}\mathbf{M}_{\mathop{\mathrm{C}P-even}}^{2}R_{\mathop{\mathrm{C}P-even}}^{\dagger}=\mathop{\mathrm{d}iag}\left(m_{h}^{2},m_{H_{1}}^{2},m_{H_{2}}^{2}\right).$
(83)

The squared CP-odd mass matrix in the interaction basis of
$\left(\mathop{\mathrm{I}m}H_{u}^{0},\mathop{\mathrm{I}m}H_{d}^{0},\mathop{\mathrm{I}m}\phi\right)$
features:

$\mathbf{M}_{\mathop{\mathrm{C}P-odd}}^{2}=\left(\begin{array}[]{ccc}\frac{v_{2}v_{3}^{2}\lambda_{5}}{2v_{1}}&\frac{1}{2}v_{3}^{2}\lambda_{5}&\sqrt{2}v_{2}v_{3}\lambda_{5}\\\
\frac{1}{2}v_{3}^{2}\lambda_{5}&\frac{v_{1}v_{3}^{2}\lambda_{5}}{2v_{2}}&\sqrt{2}v_{1}v_{3}\lambda_{5}\\\
\sqrt{2}v_{2}v_{3}\lambda_{5}&\sqrt{2}v_{1}v_{3}\lambda_{5}&4v_{1}v_{2}\lambda_{5}\\\
\end{array}\right)$ (84)

After diagonalizing the squared CP-odd mass matrix, it reveals masses of one
CP-odd scalar $A_{1}$ and two Goldstone bosons $G_{Z}$ and $G_{Z^{\prime}}$ in
the physical basis of $\left(G_{Z},A_{1},G_{Z^{\prime}}\right)$, which will be
the SM $Z$ gauge boson and a massive neutral $Z^{\prime}$ gauge boson as
follows:

$R_{\mathop{\mathrm{C}P-odd}}\mathbf{M}_{\mathop{\mathrm{C}P-odd}}^{2}R_{\mathop{\mathrm{C}P-odd}}^{\dagger}=\mathop{\mathrm{d}iag}\left(0,m_{A_{1}},0\right).$
(85)

The squared charged mass matrix in the interaction basis of
$\left(H_{1}^{\pm},H_{2}^{\pm}\right)$ has the form of:

$\textbf{M}_{\mathop{\mathrm{c}harged}}^{2}=\left(\begin{array}[]{cc}\lambda_{4}v_{2}^{2}+\frac{v_{3}^{2}\lambda_{5}v_{2}}{2v_{1}}&\frac{1}{2}\lambda_{5}v_{3}^{2}+v_{1}v_{2}\lambda_{4}\\\
\frac{1}{2}\lambda_{5}v_{3}^{2}+v_{1}v_{2}\lambda_{4}&\lambda_{4}v_{1}^{2}+\frac{v_{3}^{2}\lambda_{5}v_{1}}{2v_{2}}\\\
\end{array}\right)$ (86)

After diagonalizing the squared charged mass matrix, it reveals masses of one
non-SM charged scalar and one Goldstone boson for the SM $W^{\pm}$ gauge
boson. The squared physical charged mass matrix can be written in the basis of
$\left(G_{W^{\pm}},H^{\pm}\right)$:

$R_{\mathop{\mathrm{c}harged}}\mathbf{M}_{\mathop{\mathrm{c}harged}}^{2}R_{\mathop{\mathrm{c}harged}}^{\dagger}=\mathop{\mathrm{d}iag}\left(0,m_{H^{\pm}}\right).$
(87)

###  VII.3 SM Higgs diphoton signal strength

As the BSM model under consideration generates the scalars which are the SM
Higgs and the CP-even and -odd scalars, the SM Higgs diphoton signal strength
should be discussed and we will see that the Higgs diphoton signal strength is
quite close to the SM prediction under the decoupling limit. The decay rate
for the $h\rightarrow\gamma\gamma$ reads Hernandez:2021tii :

$\Gamma\left(h\rightarrow\gamma\gamma\right)=\frac{\alpha_{\mathop{\mathrm{e}m}}^{2}m_{h}^{3}}{256\pi^{3}v^{2}}\lvert\sum_{f}a_{hff}N_{C}Q_{f}^{2}F_{1/2}(\rho_{f})+a_{hWW}F_{1}(\rho_{W})+\frac{C_{hH^{\pm}H^{\mp}}v}{2m_{H^{\pm}}^{2}}F_{0}(\rho_{H_{k}^{\pm}})\rvert^{2},$
(88)

where $\alpha_{\mathop{\mathrm{e}m}}$ is the fine structure constant,
$\rho_{i}$ is the mass ratio $\rho_{i}=m_{h}^{2}/(4M_{i}^{2}$ with
$M_{i}=m_{f},M_{W}$, $N_{C}$ is the number of colors, and $Q_{f}$ is the
electric charge of the fermion running in the loop, $C_{hH^{\pm}H^{\mp}}$ is
the coupling constant between the SM Higgs $h$ and a charged Higgs pair
$H^{\pm}$, $a_{htt}$ and $a_{hWW}$ are the deviation factors

$\begin{split}a_{htt}&\simeq 1,\\\
a_{hWW}&=\frac{1}{\sqrt{v_{1}^{2}+v_{2}^{2}}}\frac{\partial}{\partial
h}\sum_{i,j=1,2,3}v_{i}\left(R_{\mathop{\mathrm{C}P-even}}^{T}\right)_{ij}\left(h,H_{1},H_{2}\right)_{j}^{T}=\frac{v_{1}}{\sqrt{v_{1}^{2}+v_{2}^{2}}},\end{split}$
(89)

which are very close to unity in the SM, and lastly $F_{0,1/2,1}$ are loop
functions for scalar, fermion and boson respectively. The loop functions read
Hernandez:2021tii :

$\begin{split}F_{0}(z)&=-\left(z-f(z)\right)z^{-2},\\\
F_{1/2}(z)&=2\left(z+(z-1)f(z)\right)z^{-2},\\\
F_{1}(z)&=-2\left(2z^{2}+3z+3(2z-1)f(z)\right)z^{-2}.\end{split}$ (90)

where

$f(z)=\left\\{\begin{array}[]{lcc}\arcsin^{2}\sqrt{2}&\text{for}&z\leq 1\\\
&&\\\
-\frac{1}{4}\left(\ln\left(\frac{1+\sqrt{1-z^{-1}}}{1-\sqrt{1-z^{-1}}-i\pi}\right)^{2}\right)&\text{for}&z>1\end{array}\right.$
(91)

The defined decay width for $h\rightarrow\gamma\gamma$ is appeared in the
Higgs diphoton signal strength $R_{\gamma\gamma}$, which takes the form of:

$R_{\gamma\gamma}=\frac{\sigma\left(pp\rightarrow
h\right)\Gamma\left(h\rightarrow\gamma\gamma\right)}{\sigma\left(pp\rightarrow
h\right)_{\mathop{\mathrm{S}M}}\Gamma\left(h\rightarrow\gamma\gamma\right)_{\mathop{\mathrm{S}M}}}\simeq
a_{htt}^{2}\frac{\Gamma\left(h\rightarrow\gamma\gamma\right)}{\Gamma\left(h\rightarrow\gamma\gamma\right)_{\mathop{\mathrm{S}M}}}.$
(92)

The experimental $R_{\gamma\gamma}$ is given by ATLAS and CMS collaborations
ATLAS:2019nkf ; CMS:2018piu :

$R_{\gamma\gamma}^{\mathop{\mathrm{A}TLAS}}=0.96\pm 0.14.,\quad
R_{\gamma\gamma}^{\mathop{\mathrm{C}MS}}=1.18_{-0.14}^{+0.17}$ (93)

###  VII.4 Vacuum stability

In order for all the scalars appearing in this BSM model to be physical, it is
important to secure vacuum stability of the scalar potential. We generally
follow the argument for vacuum stability covered in our previous work
Hernandez:2021tii , however there is a big difference between the previous
work and the current work: the SM gauge symmetry is extended by global
$U(1)^{\prime}$ symmetry in the previous work so there is one more CP-odd
scalar, whereas the freedom for the CP-odd scalar is replaced by $Z^{\prime}$
gauge boson as the SM gauge symmetry is extended by local $U(1)^{\prime}$
symmetry in the current work. The first condition for vacuum stability is the
mass parameters $\mu_{1,2,3}^{2}$ given by the minimization condition should
be negative otherwise it would yield massless scalar particles

$\begin{split}\mu_{1}^{2}&=-2\lambda_{1}v_{1}^{2}-\lambda_{3}v_{2}^{2}+\frac{\lambda_{5}v_{2}v_{3}^{2}}{2v_{1}}-\lambda_{7}v_{3}^{2}=-2\lambda_{1}v_{1}^{2}-\lambda_{3}v_{2}^{2}<0,\\\
\mu_{2}^{2}&=-2\lambda_{2}v_{2}^{2}-\lambda_{3}v_{1}^{2}+\frac{\lambda_{5}v_{1}v_{3}^{2}}{2v_{2}}-\lambda_{8}v_{3}^{2}=-2\lambda_{2}v_{2}^{2}+\lambda_{3}v_{1}^{2}-\lambda_{8}v_{3}^{2}<0,\\\
\mu_{3}^{2}&=-2\lambda_{6}v_{3}^{2}-\lambda_{7}v_{1}^{2}-\lambda_{8}v_{2}^{2}+2\lambda_{5}v_{1}v_{2}=\frac{4\lambda_{3}v_{1}^{2}v_{2}^{2}}{v_{3}^{2}}-\lambda_{6}v_{3}^{2}-\lambda_{8}v_{2}^{2}<0,\end{split}$
(94)

where at the second equality of each equation the decoupling limit of Equation
81 was used. From the squared mass parameter $\mu_{2}^{2}$, one can know that
the quartic coupling constant $\lambda_{3}$ must be negative, taking into
account hierarchy among vevs
$(\mathcal{O}(v_{2})\ll\mathcal{O}(v_{1})\approx\mathcal{O}(v_{3}))$ as we
will see soon in the following numerical discussion. The negative sign of
$\lambda_{3}$ also leads to negative sign of $\lambda_{5}$ and $\lambda_{7}$
due to the decoupling limit. As order of the vev $v_{1}$ is compatible to that
of $v_{3}$, it is difficult to derive analytically further conditions for
vacuum stability at the current stage as the other parameters depend on their