CTP-SCU/2021021, APCTP Pre2021-013

# Quark and lepton flavor model with leptoquarks in a modular $A_{4}$ symmetry

Takaaki Nomura<EMAIL_ADDRESS>College of Physics, Sichuan University,
Chengdu 610065, China  Hiroshi Okada<EMAIL_ADDRESS>Asia Pacific
Center for Theoretical Physics, Pohang 37673, Republic of Korea  Department of
Physics, Pohang University of Science and Technology, Pohang 37673, Republic
of Korea  Yuta Orikasa<EMAIL_ADDRESS>Institute of Experimental
and Applied Physics, Czech Technical University in Prague, 110 00 Prague 1,
Czech Republic

###### Abstract

We propose a quark-lepton model via leptoquarks and modular $A_{4}$ symmetry.
Since the neutrino mass is induced at one-loop level mediated by down quarks
as well as leptoquarks, we have to explain lepton and quark masses and mixings
with a single modulus $\tau$. Here, we find predictions for lepton and quark
sectors with unified modulus $\tau$, and show several constraints originating
from leptoquarks.

##  I Introduction

Since lepto-quark(LQ) bosons connect lepton and quark sectors, these models
potentially explain several new physics beyond the standard model (SM); e.g.,
lepton(muon or electron) anomalous magnetic dipole moment ($g-2$)
Bauer:2015knc ; Chen:2016dip ; Chen:2017hir ; Nomura:2021oeu ;
ColuccioLeskow:2016dox , $B$ meson decays such as $b\to s\mu\bar{\mu}$
Sahoo:2015wya ; Becirevic:2016yqi ; Chen:2016dip ; Chen:2017hir ;
Nomura:2021oeu ; Cai:2017wry  and $b\to
c\ell\bar{\nu}_{\ell}$($\ell=e,\mu,\tau$) Sakaki:2013bfa ; Becirevic:2016yqi ;
Chen:2017hir ; Cai:2017wry  111The anomaly of $b\to c\ell\bar{\nu}_{\ell}$
processes are observed in experiments Huschle:2015rga ; Lees:2012xj ;
Lees:2013uzd ; Hirose:2016wfn ; Abdesselam:2016cgx ; Aaij:2015yra ;
Aaij:2017deq , and LQ model is one of the most promising explanations on this
anomaly., and nonzero neutrino masses AristizabalSierra:2007nf ;
Cheung:2016fjo ; Cai:2017wry . Especially, muon $g-2$ anomaly is recently
reported by E989 experiment at Fermilab combining BNL result Abi:2021gix , and
its value is deviated from the SM by 4.2$\sigma$ as follows:

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

Also, the LHCb collaboration Aaij:2021vac  recently reported anomaly of rare
$B$ meson decays of $b\to s\mu\bar{\mu}$ that is understood as violation of
lepton universality. The updated result is given by

$\displaystyle\frac{BR(B^{+}\to K^{+}\mu^{-}\mu^{+})}{BR(B^{+}\to
K^{+}e^{-}e^{+})}=0.846_{-0.039-0.012}^{+0.042+0.013}\quad(1.1{\rm
GeV}^{2}<q^{2}<6{\rm GeV}^{2}),$ (I.2)

where first(second) uncertainty is statistical(systematic) one and $q^{2}$ is
the invariant mass squared for dilepton. In addition to the above
phenomenologies, interestingly, the nonzero Majorana neutrino mass at one-loop
can be realized without any additional symmetries by introducing appropriate
LQs Cheung:2016fjo . This may be natural realization of tiny neutrino mass
model due to loop suppression.

Considering above issues, one finds that Yukawa flavor structure is also very
important to explain them. Recently, powerful symmetries to restrict the
number of parameters in Yukawa couplings, so called ”modular flavor
symmetries”, were proposed by authors in refs. Feruglio:2017spp ;
deAdelhartToorop:2011re , in which they have applied modular originated non-
Abelian discrete flavor symmetries to quark and lepton sectors. One remarkable
advantage of applying this symmetries is that dimensionless couplings of model
can be transformed into non-trivial representations under those symmetries,
and all the dimensionless values are uniquely fixed once modulus is determined
in fundamental region. We then do not need the scalar fields to obtain a
predictive mass matrix. Along the line of this idea, a vast reference has
recently appeared in the literature, e.g., $A_{4}$ Feruglio:2017spp ;
Criado:2018thu ; Kobayashi:2018scp ; Okada:2018yrn ; Nomura:2019jxj ;
Okada:2019uoy ; deAnda:2018ecu ; Novichkov:2018yse ; Nomura:2019yft ;
Okada:2019mjf ; Ding:2019zxk ; Nomura:2019lnr ; Kobayashi:2019xvz ;
Asaka:2019vev ; Zhang:2019ngf ; Gui-JunDing:2019wap ; Kobayashi:2019gtp ;
Nomura:2019xsb ; Wang:2019xbo ; Okada:2020dmb ; Okada:2020rjb ; Behera:2020lpd
; Behera:2020sfe ; Nomura:2020opk ; Nomura:2020cog ; Asaka:2020tmo ;
Okada:2020ukr ; Nagao:2020snm ; Okada:2020brs ; Yao:2020qyy ; Chen:2021zty ;
Kobayashi:2021bgy ; Kashav:2021zir ; Okada:2021qdf , $S_{3}$ Kobayashi:2018vbk
; Kobayashi:2018wkl ; Kobayashi:2019rzp ; Okada:2019xqk ; Mishra:2020gxg ;
Du:2020ylx , $S_{4}$ Penedo:2018nmg ; Novichkov:2018ovf ; Kobayashi:2019mna ;
King:2019vhv ; Okada:2019lzv ; Criado:2019tzk ; Wang:2019ovr ; Zhao:2021jxg ;
King:2021fhl ; Ding:2021zbg ; Zhang:2021olk ; gui-jun , $A_{5}$
Novichkov:2018nkm ; Ding:2019xna ; Criado:2019tzk , double covering of $A_{5}$
Wang:2020lxk ; Yao:2020zml ; Wang:2021mkw , larger groups Baur:2019kwi ,
multiple modular symmetries deMedeirosVarzielas:2019cyj , and double covering
of $A_{4}$ Liu:2019khw ; Chen:2020udk , $S_{4}$ Novichkov:2020eep ;
Liu:2020akv , and the other types of groups Kikuchi:2020nxn ; Almumin:2021fbk
; Ding:2021iqp ; Feruglio:2021dte ; Kikuchi:2021ogn ; Novichkov:2021evw  in
which masses, mixing, and CP phases for the quark and/or lepton have been
predicted 222For interest readers, we provide some literature reviews, which
are useful to understand the non-Abelian group and its applications to flavor
structure Altarelli:2010gt ; Ishimori:2010au ; Ishimori:2012zz ;
Hernandez:2012ra ; King:2013eh ; King:2014nza ; King:2017guk ; Petcov:2017ggy
.. Moreover, a systematic approach to understand the origin of CP
transformations has been discussed in Ref. Baur:2019iai , and CP/flavor
violation in models with modular symmetry was discussed in Refs.
Kobayashi:2019uyt ; Novichkov:2019sqv ; Kobayashi:2021bgy ; 1869542 , and a
possible correction from Kähler potential was discussed in Ref. Chen:2019ewa .
Furthermore, systematic analysis of the fixed points (stabilizers) has been
discussed in Ref. deMedeirosVarzielas:2020kji . A very recent paper of Ref.
Ishiguro:2020tmo  finds a favorable fixed point $\tau=\omega$ among three
fixed points, which are the fundamental domain of PSL$(2,Z)$, by
systematically analyzing the stabilized moduli values in the possible
configurations of flux compactifications as well as investigating the
probabilities of moduli values. It is then interesting to discuss a LQ model
under the framework of modular flavor symmetry since we are motivated to
consider lepton and quark sector together as a LQ connect these sectors and
some predictions in both sector can be expected.

In this paper, we focus on the quark and lepton masses and mixings based on a
LQ model in ref. Cheung:2016fjo , introducing modular $A_{4}$ symmetry to
reduce free parameters of Yukawa couplings. Since the quark sector connects to
the lepton sector via LQ, charge assignments for quarks(leptons) directly
affect the leptons(quarks). In this sense, it would be a good motivation
towards unification of quark and lepton flavor in $A_{4}$ modular symmetry.

This paper is organized as follows. In Sec. II, we review our model of quark
and lepton. In Sec. III, we have numerical analysis and show several results
for normal and inverted hierarchies. We conclude in Sec. IV. In appendix, we
summarize several features of modular $A_{4}$ symmetry.

##  II Model

|  $\eta$ |  $\Delta$
---|---|---
$SU(3)_{C}$ | $\bm{3}$ | $\bar{\bm{3}}$
$SU(2)_{L}$ | $\bm{2}$ | $\bm{3}$
$U(1)_{Y}$ | $\frac{1}{6}$ | $\frac{1}{3}$
$A_{4}$ | $\bm{1}$ | $\bm{1}$
$-k_{I}$ | $-2$ | $-2$
Table 1: Charge assignments of the LQ bosons $\eta$ and $\Delta$ under $SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}\times A_{4}$ where $k_{I}$ is the number of modular weight. | Fermions
---|---
|  $Q_{L}$ |  $\bar{u}_{R}$ |  $\bar{d}_{R}$ |  $L_{L}$ |  $\bar{e}_{R}$
$SU(3)_{C}$ | $\bm{3}$ | $\bar{\bm{3}}$ | $\bar{\bm{3}}$ | $\bm{1}$ | $\bm{1}$
$SU(2)_{L}$ | $\bm{2}$ | $\bm{1}$ | $\bm{1}$ | $\bm{2}$ | $\bm{1}$
$U(1)_{Y}$ | $\frac{1}{6}$ | $-\frac{2}{3}$ | $\frac{1}{3}$ | $-\frac{1}{2}$ | $1$
$A_{4}$ | $\bm{3}$ | $\bm{1,1^{\prime\prime},1^{\prime}}$ | $\bm{1,1^{\prime\prime},1^{\prime}}$ | $\bm{3}$ | $\bm{1,1^{\prime\prime},1^{\prime}}$
$-k_{I}$ | $-2$ | $-4$ | $0$ | $-2$ | $0$
Table 2: Charge assignments of the SM fermions under $SU(3)_{C}\times
SU(2)_{L}\times U(1)_{Y}\times A_{4}$ where $k_{I}$ is the number of modular
weight.

In this section, we review our model. It is known that introducing proper
leptoquarks lead us to a radiative seesaw model without any additional
symmetries such as $Z_{2}$. Here, we introduce two types of leptoquarks $\eta$
and $\Delta$ based on Ref. Cheung:2016fjo . The color-triplet $\eta$ has
$SU(2)_{L}$ doublet with $1/6$ hypercharge, and the color-antitriplet $\Delta$
has $SU(2)_{L}$ triplet with $1/3$ hypercharge, where these new bosons and
their charges are summarized in Table 1. Then, the valid Lagrangian to induce
the quark and lepton mass matrices is given by

$\displaystyle-\mathcal{L}_{Y}^{q}$
$\displaystyle=y^{u}_{ij}\bar{u}_{R_{i}}(i\sigma_{2})H^{*}Q_{L_{j}}+y^{d}_{ij}\bar{d}_{R_{i}}HQ_{L_{j}}+{\rm
h.c.},$ (II.1) $\displaystyle-\mathcal{L}_{Y}^{\ell}$
$\displaystyle=h_{ij}\bar{e}_{R_{i}}HL_{L_{j}}+{\rm h.c.}.$ (II.2)

The Lagrangian for the mixing between the quark and lepton and nontrivial
potential are given by

$\displaystyle-\mathcal{L}_{Y}^{mix}$
$\displaystyle=f_{ij}\overline{d_{R_{i}}}\eta^{T}(i\sigma_{2})L_{L_{j}}+g_{ij}\overline{Q^{c}_{L_{i}}}(i\sigma_{2})\Delta
L_{L_{j}}+{\rm h.c.},$ (II.3) $\displaystyle\mathcal{V}$
$\displaystyle\supset-\mu H^{\dagger}\Delta\eta+{\rm h.c.},$ (II.4)

where $(i,j)=1-3$ are family indices, $\sigma_{2}$ is the second Pauli matrix,
and $H$ is the SM Higgs field that develops a nonzero VEV, which is symbolized
by $\langle H\rangle\equiv v/\sqrt{2}\approx 246/\sqrt{2}$ GeV, and $H$ has
nonzero modular weight. Here, we parameterize components of the scalars as
follows:

$\displaystyle H=\left[\begin{array}[]{c}w^{+}\\\
\frac{v+\phi+iz}{\sqrt{2}}\end{array}\right],\quad\eta=\left[\begin{array}[]{c}\eta_{2/3}\\\
\eta_{-1/3}\end{array}\right],\quad\Delta=\left[\begin{array}[]{cc}\frac{\delta_{1/3}}{\sqrt{2}}&\delta_{4/3}\\\
\delta_{-2/3}&-\frac{\delta_{1/3}}{\sqrt{2}}\end{array}\right],$ (II.11)

where the subscript of the fields represents the electric charge, and $w^{+}$
and $z$ are absorbed by the longitudinal component of the $W^{+}$ and $Z$
bosons, respectively. Due to the $\mu$ term in Eq. (II.4), the charged
components with $1/3$ and $2/3$ electric charges mix each other. Here, we
parametrize their mixing matrices and mass eigenstates as follows:

$\displaystyle\left[\begin{array}[]{c}\eta_{i/3}\\\
\delta_{i/3}\end{array}\right]=O_{i}\left[\begin{array}[]{c}A_{i}\\\
B_{i}\end{array}\right],\quad
O_{i}\equiv\left[\begin{array}[]{cc}c_{a_{i}}&s_{a_{i}}\\\
-s_{a_{i}}&c_{a_{i}}\end{array}\right],\quad(i=1,2),$ (II.18)

where their masses are denoted as $m_{A_{i}}$ and $m_{B_{i}}$ respectively.
The interactions in terms of the mass eigenstates can be written as

$\displaystyle-\mathcal{L}_{Y}^{q}\approx$ $\displaystyle\
m_{u_{ij}}\bar{u}_{R_{i}}u_{L_{j}}+m_{d_{ij}}\bar{d}_{R_{i}}d_{L_{j}}+{\rm
h.c.},$ (II.19) $\displaystyle-\mathcal{L}_{Y}^{\ell}\approx$ $\displaystyle\
m_{\ell_{ij}}\bar{e}_{R_{i}}\ell_{L_{j}}+{\rm h.c.},$ (II.20)
$\displaystyle-\mathcal{L}_{Y}^{mix}\approx$ $\displaystyle\
f_{ij}\overline{d_{R_{i}}}\nu_{L_{j}}(c_{a_{1}}A_{1}^{*}+s_{a_{1}}B_{1}^{*})-\frac{g_{ij}}{\sqrt{2}}\overline{d_{L_{i}}^{c}}\nu_{L_{j}}(-s_{a_{1}}A_{1}+c_{a_{1}}B_{1})$
(II.21) $\displaystyle-
f_{ij}\overline{d_{R_{i}}}\ell_{L_{j}}(c_{a_{2}}A_{2}+s_{a_{2}}B_{2})-\frac{g_{ij}}{\sqrt{2}}\overline{u_{L_{i}}^{c}}\ell_{L_{j}}(-s_{a_{1}}A_{1}+c_{a_{1}}B_{1})$
(II.22)
$\displaystyle-{g_{ij}}\overline{d_{L_{i}}^{c}}\ell_{L_{j}}\delta_{4/3}\;+{g_{ij}}\overline{u_{L_{i}}^{c}}\nu_{L_{j}}(-s_{a_{2}}A_{2}^{*}+c_{a_{2}}B_{2}^{*})+{\rm
h.c.},$ (II.23)

where we define $m_{u_{ij}}\equiv\frac{vy^{u}_{ij}}{\sqrt{2}}$,
$m_{d_{ij}}\equiv\frac{vy^{d}_{ij}}{\sqrt{2}}$, and
$m_{\ell_{ij}}\equiv\frac{vh_{ij}}{\sqrt{2}}$.

The next task is to determine the matrices of $y^{u},y^{d},h,f,g$ via modular
$A_{4}$ symmetry. In the quark sector, we assign $Q_{L}$ to be ${\bf 3}$ and
$-2$, $\bar{u}_{R}$ to be $\\{\bf 1,1^{\prime\prime},1^{\prime}\\}$ and $-4$,
and $\bar{d}_{R}$ to be $\\{{\bf 1,1^{\prime\prime},1^{\prime}}\\}$ and $0$
under $A_{4}$ and $-k$, respectively. This assignment is the same as the one
in ref. Okada:2020rjb , and it is already known that allowed region
Okada:2019uoy . Thus, we will work on the same $\tau$ region of the lepton
sector in our numerical analysis. The up-type quark mass matrix is written as:

$\displaystyle\begin{aligned} y^{u}=\begin{pmatrix}a_{u}&0&0\\\ 0&a_{c}&0\\\
0&0&a_{t}\end{pmatrix}\left[\begin{pmatrix}f_{1}&f_{3}&f_{2}\\\
f_{2}&f_{1}&f_{3}\\\
f_{3}&f_{2}&f_{1}\end{pmatrix}+\begin{pmatrix}g_{u1}&0&0\\\ 0&g_{u2}&0\\\
0&0&g_{u3}\end{pmatrix}\begin{pmatrix}f^{\prime}_{1}&f^{\prime}_{3}&f^{\prime}_{2}\\\
f^{\prime}_{2}&f^{\prime}_{1}&f^{\prime}_{3}\\\
f^{\prime}_{3}&f^{\prime}_{2}&f^{\prime}_{1}\end{pmatrix}\right],\end{aligned}$
(II.24)

where $Y^{(6)}_{3}\equiv[f_{1},f_{2},f_{3}]^{T}$ and
$Y^{(6)}_{3^{\prime}}\equiv[f^{\prime}_{1},f^{\prime}_{2},f^{\prime}_{3}]^{T}$,
$g_{u1}=\alpha^{\prime}_{u}/\alpha_{u}$, $g_{u2}=\beta^{\prime}_{u}/\beta_{u}$
and $g_{u3}=\gamma^{\prime}_{u}/\gamma_{u}$ are complex parameters, and
$a_{u}$, $a_{c}$ and $a_{t}$ can be used to fit the masses of up quarks. The
explicit forms of $f_{i}$ and $f^{\prime}_{i}$ are summarized in Appendix.
Then $m_{u}$ is diagonalized by two unitary matrices as
$D_{u}=V_{u_{R}}^{\dagger}m_{u}V_{u_{L}}$, where $D_{u}\equiv{\rm
diag}(m_{u},m_{c},m_{t})$ is mass eigenvalues. Therefore, we find
$|D_{u}|^{2}=V^{\dagger}_{u_{L}}m_{u}^{\dagger}m_{u}V_{u_{L}}$. On the other
hand, the down-type quark mass matrix is given as:

$\displaystyle\begin{aligned} y^{d}=\begin{pmatrix}a_{d}&0&0\\\ 0&a_{s}&0\\\
0&0&a_{b}\end{pmatrix}\begin{pmatrix}y_{1}&y_{3}&y_{2}\\\ y_{2}&y_{1}&y_{3}\\\
y_{3}&y_{2}&y_{1}\end{pmatrix}\,,\end{aligned}$ (II.25)

where and $a_{d}$, $a_{s}$ and $a_{b}$ can be used to fit the masses of down
quarks, and $Y^{(2)}_{3}\equiv[y_{1},y_{2},y_{3}]^{T}$ in Appendix. Then
$m_{d}$ is diagonalized by two unitary matrices as
$D_{d}=V_{d_{R}}^{\dagger}m_{d}V_{d_{L}}$, where $D_{d}\equiv{\rm
diag}(m_{d},m_{s},m_{b})$ is mass eigenvalues. Therefore, we find
$|D_{d}|^{2}=V^{\dagger}_{d_{L}}m_{d}^{\dagger}m_{d}V_{d_{L}}$. Finally, we
get the observable mixing matrix $V_{CKM}$ as follows:

$\displaystyle V_{CKM}=V^{\dagger}_{u_{L}}V_{d_{L}}.$ (II.26)

###  II.1 Lepton sector

Now let us move on the lepton sector. We assign $L_{L}$ to be ${\bf 3}$ and
$-2$ and $\bar{e}_{R}$ to be $\\{\bf 1,1^{\prime\prime},1^{\prime}\\}$ and $0$
under $A_{4}$ and $-k$, respectively. Here, both of the leptoquark scalars are
assigned to be true $A_{4}$ singlets with $-2$ modular weight. The assignments
of $A_{4}$ and $-k$ are also summarized in Tables 1 and 2. Under these
assignments, we can write down the concrete matrices as follows:

$\displaystyle h$ $\displaystyle=\left[\begin{array}[]{ccc}a_{\ell}&0&0\\\
0&b_{\ell}&0\\\ 0&0&c_{\ell}\\\
\end{array}\right]\left[\begin{array}[]{ccc}y_{1}&y_{3}&y_{2}\\\
y_{2}&y_{1}&y_{3}\\\ y_{3}&y_{2}&y_{1}\\\ \end{array}\right],$ (II.33)
$\displaystyle f$ $\displaystyle=\left[\begin{array}[]{ccc}a_{\eta}&0&0\\\
0&b_{\eta}&0\\\ 0&0&c_{\eta}\\\
\end{array}\right]\left[\begin{array}[]{ccc}y^{\prime}_{1}&y^{\prime}_{3}&y^{\prime}_{2}\\\
y^{\prime}_{2}&y^{\prime}_{1}&y^{\prime}_{3}\\\
y^{\prime}_{3}&y^{\prime}_{2}&y^{\prime}_{1}\\\ \end{array}\right],$ (II.40)
$\displaystyle g$
$\displaystyle=aY^{(6)}_{1}\left[\begin{array}[]{ccc}1&0&0\\\ 0&0&1\\\
0&1&0\\\
\end{array}\right]+\frac{b}{3}\left[\begin{array}[]{ccc}2f_{1}&-f_{3}&-f_{2}\\\
-f_{3}&2f_{2}&-f_{1}\\\ -f_{2}&-f_{1}&2f_{3}\\\
\end{array}\right]+\frac{c}{2}\left[\begin{array}[]{ccc}0&f_{3}&-f_{2}\\\
-f_{3}&0&f_{1}\\\ f_{2}&-f_{1}&0\\\ \end{array}\right]$ (II.50)
$\displaystyle+\frac{b^{{}^{\prime}}}{3}\left[\begin{array}[]{ccc}2f^{{}^{\prime}}_{1}&-f^{{}^{\prime}}_{3}&-f^{{}^{\prime}}_{2}\\\
-f^{{}^{\prime}}_{3}&2f^{{}^{\prime}}_{2}&-f^{{}^{\prime}}_{1}\\\
-f^{{}^{\prime}}_{2}&-f^{{}^{\prime}}_{1}&2f^{{}^{\prime}}_{3}\\\
\end{array}\right]+\frac{c^{{}^{\prime}}}{2}\left[\begin{array}[]{ccc}0&f^{{}^{\prime}}_{3}&-f^{{}^{\prime}}_{2}\\\
-f^{{}^{\prime}}_{3}&0&f^{{}^{\prime}}_{1}\\\
f^{{}^{\prime}}_{2}&-f^{{}^{\prime}}_{1}&0\\\ \end{array}\right]$ (II.57)

where $Y^{(4)}_{3}\equiv[y^{\prime}_{1},y^{\prime}_{2},y^{\prime}_{3}]^{T}$ in
Appendix.

$\nu_{L_{i}}$$\nu_{L_{j}}$$A_{1}/B_{1}$$d_{a}$ Figure 1:  One-loop diagrams
for generating the neutrino mass matrix.

The charged-lepton sector after spontaneous symmetry breaking is given by

$\displaystyle-\mathcal{L}_{Y}^{\ell}=v\frac{h_{ij}}{\sqrt{2}}\ell_{L_{i}}e_{R_{j}}+{\rm
h.c.}.$ (II.58)

Then $m_{\ell}(\equiv vh/\sqrt{2})$ is diagonalized by two unitary matrices as
$D_{\ell}=V_{\ell_{R}}^{\dagger}m_{\ell}V_{\ell_{L}}$, where
$D_{\ell}\equiv{\rm diag}(m_{e},m_{\mu},m_{\tau})$ is mass eigenvalues.
Therefore, we find
$|D_{\ell}|^{2}=V^{\dagger}_{\ell_{L}}m_{\ell}^{\dagger}m_{\ell}V_{\ell_{L}}$.
The active neutrino mass matrix $m_{\nu}$ is given at one-loop level through
the following interactions:

$\displaystyle-{L}_{Y}^{\ell}=F_{aj}\overline{d^{\prime}_{R_{a}}}\nu_{L_{j}}(c_{a_{1}}A_{1}+s_{a_{1}}B_{1})-G_{aj}\overline{d^{\prime
c}_{L_{a}}}\nu_{L_{j}}(-s_{a_{1}}A_{1}+c_{a_{1}}B_{1}),$ (II.59)

where $F\equiv V^{\dagger}_{d_{R}}f$ and $G\equiv V^{T}_{d_{L}}g$ and
$d^{\prime}$ is mass eigenstate. Then, the neutrino masss matrix in Fig. 1 is
given at one-loop level as follows:

$\displaystyle(m_{\nu})_{ij}=s_{2a_{1}}\frac{3}{4(4\pi)^{2}}\left[1-\frac{m^{2}_{A_{1}}}{m^{2}_{B_{1}}}\right]\sum_{i=1}^{3}\left[F^{T}_{ia}D_{d_{a}}G_{aj}+G^{T}_{ia}D_{d_{a}}F_{aj}\right]F_{I}(r_{A_{1}},r_{D_{d_{i}}}),$
(II.60) $\displaystyle F_{I}(r_{1},r_{2})=\frac{r_{1}(r_{2}-1)\ln
r_{1}-r_{2}(r_{1}-1)\ln
r_{2}}{(r_{1}-1)(r_{2}-1)(r_{1}-r_{2})},\quad(r_{1}\neq 1),$ (II.61)

where we define $r_{A_{1}}\equiv(m_{A_{1}}/m_{B_{1}})^{2}$ and
$r_{D_{d_{i}}}\equiv(D_{d_{i}}/m_{B_{1}})^{2}$. $m_{\nu}$ is diagonalzied by a
unitary matrix $V_{\nu}$; $D_{\nu}\equiv V_{\nu}^{T}m_{\nu}V_{\nu}$. Here, we
define a modified neutrino mass matrix as $\tilde{m}_{\nu}\equiv
m_{\nu}/s_{2a_{1}}$. Then, we rewrite this diagonalization in terms of the
modified form $\tilde{D}_{\nu}\equiv V_{\nu}^{T}\tilde{m}_{\nu}V_{\nu}$. Thus,
we fix $s_{2a_{1}}$ by

$\displaystyle({\rm NH}):\ s_{2a_{1}}^{2}=\frac{|\Delta m_{\rm
atm}^{2}|}{\tilde{D}_{\nu_{3}}^{2}-\tilde{D}_{\nu_{1}}^{2}},\quad({\rm IH}):\
s_{2a_{1}}^{2}=\frac{|\Delta m_{\rm
atm}^{2}|}{\tilde{D}_{\nu_{2}}^{2}-\tilde{D}_{\nu_{3}}^{2}},$ (II.62)

where $\tilde{m}_{\nu}$ is diagonalized by
$V^{\dagger}_{\nu}(\tilde{m}_{\nu}^{\dagger}\tilde{m}_{\nu})V_{\nu}=(\tilde{D}_{\nu_{1}}^{2},\tilde{D}_{\nu_{2}}^{2},\tilde{D}_{\nu_{3}}^{2})$
and $\Delta m_{\rm atm}^{2}$ is the atmospheric neutrino mass-squared
difference. Here, NH and IH respectively stand for the normal hierarchy and
the inverted hierarchy. Subsequently, the solar neutrino mass-squared
difference is depicted in terms of $s_{2a_{1}}$ as follows:

$\displaystyle\Delta m_{\rm
sol}^{2}={s_{2a_{1}}^{2}}({\tilde{D}_{\nu_{2}}^{2}-\tilde{D}_{\nu_{1}}^{2}}).$
(II.63)

This should be within the experimental value, where we adopt NuFit 5.0
Esteban:2020cvm  to our numerical analysis later. The neutrinoless double beta
decay is also given by

$\displaystyle\langle
m_{ee}\rangle=s_{2a_{1}}|\tilde{D}_{\nu_{1}}\cos^{2}\theta_{12}\cos^{2}\theta_{13}+\tilde{D}_{\nu_{2}}\sin^{2}\theta_{12}\cos^{2}\theta_{13}e^{i\alpha_{2}}+\tilde{D}_{\nu_{3}}\sin^{2}\theta_{13}e^{i(\alpha_{3}-2\delta_{CP})}|,$
(II.64)

which may be able to observed by KamLAND-Zen in future KamLAND-Zen:2016pfg .
The observed mixing matrix of lepton sector Maki:1962mu  is given by $V_{\rm
PMNS}\equiv V^{\dagger}_{\ell_{L}}V_{\nu}$.

##  III Numerical analysis

Here, we perform numerical analysis. Before searching for allowed region, we
fix some mass parameters as $m_{A_{2}}=m_{A_{1}}$ and
$m_{B_{2}}=m_{\delta}=m_{B_{1}}$, where we require degenerate masses for the
components of $\eta$ and $\Delta$ to suppress the oblique parameters $\Delta
S$ and $\Delta T$. Notice here that our theoretical parameters
$a_{u,c,t},a_{d,s,b},a_{\ell},b_{\ell},c_{\ell}$ are used to determined the
experimental masses for quarks and charged-leptons. Thus, only the following
input parameters are randomly selected in the range of

$\displaystyle(m_{A_{1}},m_{B_{1}})\in[1\,,100\,]\ \text{TeV},$
$\displaystyle|g_{u1,u2,u3}|\in\left[10^{-5},\
1.5\,\right],\quad(|a_{\eta}|,|b_{\eta}|,|c_{\eta}|,|a|,|b|,|c|,|b^{\prime}|,|c^{\prime}|)\in\left[10^{-5},\
10\,\right].$ (III.1)

Above the range, we have numerical analysis in cases for quark and lepton,
where experimental data in the quark sector should be within the range at
3$\sigma$. While the one in the lepton sector is discussed in the range within
3$\sigma$(yellow dots) and 5$\sigma$(red dots) applying $\chi^{2}$ analysis in
Nufit 5.0.

###  III.1 NH

For NH case, we show our results of lepton sector in Figs. 5, 5, 5, 5. In Fig.
5, allowed value of $\tau$ is shown where yellow(red) points present the
values within 3(5)$\sigma$. One finds that allowed space is rather localized.
Especially, the region at nearby $\tau\sim 1.75i$ would be interesting since
it is close to the fixed point that has a remnant $Z_{3}$ symmetry. In Fig. 5,
we demonstrate allowed region of $\delta_{CP}$ in terms of $\sum m_{i}$. $\sum
m_{i}$ is rather localized at $0.06-0.08$ eV while whole the region is allowed
for $\delta_{CP}$. Moreover, almost all the points are within the cosmological
constraint $\sim 0.12$ eV pdg . In Fig. 5, we present allowed region of
neutrinoless double beta decay $\langle m_{ee}\rangle$ in terms of the
lightest neutrino mass $m_{1}$. $\langle m_{ee}\rangle$ is allowed up to
$0.025$ eV while $m_{1}$ is allowed up to $0.0025$ eV. Moreover, allowed
region of $m_{1}$ is localized at around $10^{-6}$ eV indicating tiny mass of
the lightest neutrino mass. In Fig. 5, we depict allowed region of Majorana
phases $\alpha_{21}$ and $\alpha_{31}$. Both are allowed for whole the region
but there would be tendency that $\alpha_{21}$ is localized at around 180∘.

Figure 2: Allowed value of $\tau$. Yellow and red points present the values
within 3 and 5$\sigma$.

Figure 3: Allowed region of $\delta_{CP}$ in terms of $\sum m_{i}$.

Figure 4: Allowed region of the mass of neutrinoless double beta decay in
terms of the lightest neutrino mass.

Figure 5: Allowed region of Majorana phases $\alpha_{21}$ in terms of
$\alpha_{31}$.

In addition to the lepton sector, we will search for our allowed region of
quark sector in Figs. 8, 8, 8. Here, the dotted red line at 3$\sigma$ interval
while the black line is best fit value. And the yellow(red) points correspond
to 3(5)$\sigma$ of the lepton sector where $\tau$ is commonly used.

In Fig. 8, we show the CP phase of quark $\delta$ in term of (1, 3) component
of CKM matrix; $|V_{ub}|$, and find whole the region is allowed at 3$\sigma$
interval. In Fig. 8, we show $|V_{ub}|$ and $|V_{td}|$, and found that there
is a weak linearly correlation between them. In Fig. 8, we show $|V_{cb}|$ and
$|V_{ub}|$, and find that there is also a weak linear correlation between
them.

Figure 6: The CP phase of quark $\delta$ versus (1, 3) component of CKM
matrix. The red dashed lines represent 3$\sigma$ experimental bounds.

Figure 7: $|V_{ub}|$ versus $|V_{td}|$.

Figure 8: $|V_{cb}|$ versus $|V_{ub}|$.

Bench mark point for NH: We also give a benchmark point to satisfy the quark
and lepton masses and mixings as well as phases in the left sides of Tables 3
and 4, where we extracted a value at nearby $\tau=1.75i$. The corresponding
lepton and neutrino mixings are given by

$\displaystyle V_{\ell_{L}}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.75+0.00017i&0.35-0.000026i&-0.56\\\
0.20-0.000012i&-0.68+0.000074i&-0.70-0.000040i\\\
-0.63+0.00013i&-0.64+0.000077i&0.44+0.000044i\\\ \end{array}\right],$ (III.5)
$\displaystyle V_{\nu}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.91-0.20i&0.26-0.075i&-0.22-0.027i\\\
0.23+0.018i&-0.045-0.28i&-0.87-0.33i\\\ -0.13+0.24i&-0.25+0.88i&-0.25-0.15i\\\
\end{array}\right].$ (III.9)

And the quark mixings are given by

$\displaystyle V_{u_{L}}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.75-0.057i&0.41+0.24i&-0.46+0.0034i\\\
-0.62+0.017i&-0.31-0.068i&0.71-0.0045i\\\
0.19-0.071i&0.76+0.30i&0.53-0.0015i\\\ \end{array}\right],$ (III.13)
$\displaystyle V_{d_{L}}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.64+0.000077i&-0.63+0.00013i&-0.44-0.000044i\\\
-0.68+0.000074i&0.20-0.000012i&0.70+0.000040i\\\
0.35-0.000026i&-0.75+0.00017i&0.56\\\ \end{array}\right].$ (III.17)

Lepton | NH($\tau\approx 1.75i$) | IH($\tau\approx 1.06i$) | IH($\tau\approx 1.76i$)
---|---|---|---
$\tau$ | $-0.0000945+1.75i$ | $-0.000689+1.06i$ | $-0.000829+1.76i$
$a_{\eta}$ | $-0.23-1.4i$ | $-0.31+0.013i$ | $-4.1+4.3i$
$b_{\eta}$ | $-0.38+1.3i$ | $-0.045-0.027i$ | $-0.0014+0.0032i$
$c_{\eta}$ | $0.0077-0.031i$ | $0.0014-0.000047i$ | $-0.0023+0.0035i$
$a$ | $0.00016+0.00011i$ | $3.0+0.98i$ | $0.0085+0.025i$
$b$ | $-0.017+0.0013i$ | $0.00014+0.00015i$ | $-0.096+0.041i$
$c$ | $-0.00030-0.00010i$ | $0.0056-0.0021i$ | $-1.8+4.5i$
$b^{\prime}$ | $0.00014-0.000016i$ | $0.000024+0.000016i$ | $-0.00011+0.00011i$
$c^{\prime}$ | $-0.21-0.27i$ | $-0.15+0.031i$ | $-0.0000034-0.000060i$
$[\alpha_{e},\beta_{e},\gamma_{e}]$ | $[0.0002,9.3\times 10^{-7},0.003]$ | $[0.0005,10^{-5},0.007]$ | $[9.1\times 10^{-7},1.9\times 10^{-4},0.003]$
$\sin^{2}\theta_{12}$ | $0.32$ | $0.28$ | $0.33$
$\sin^{2}\theta_{23}$ | $0.56$ | $0.46$ | $0.58$
$\sin^{2}\theta_{13}$ | $0.024$ | $0.024$ | $0.022$
$\delta_{CP}^{\ell}$ | $328^{\circ}$ | $170^{\circ}$ | $335^{\circ}$
$[\alpha_{21},\,\alpha_{31}]$ | $[169^{\circ},\,336^{\circ}]$ | $[167^{\circ},\,159^{\circ}]$ | $[157^{\circ},\,130^{\circ}]$
$\sum m_{i}$ |  $0.071$ eV |  $0.11$ eV |  $0.11$ eV
$s_{2a_{1}}$ | $4.6\times 10^{-9}$ | $5.7\times 10^{-5}$ | $1.9\times 10^{-9}$
$\langle m_{ee}\rangle$ |  $3.2$ meV |  $21$ meV |  $19$ meV
$[m_{A_{1}},m_{B_{1}}]$ |  $[18,6.0]$ TeV |  $[37,37]$ TeV |  $[33,39]$ TeV
$\sqrt{\chi^{2}}$ | $2.9$ | $4.8$ | $4.5$
Table 3: Numerical values of parameters and observables at the best fit points of NH and IH. Quark | NH($\tau\approx 1.75i$) | IH($\tau\approx 1.06i$) | IH($\tau\approx 1.76i$)
---|---|---|---
$\tau$ | $-0.0000945+1.75i$ | $-0.000689+1.06i$ | $-0.000829+1.76i$
$a_{u}$ | $1.3\times 10^{-7}$ | $8.2\times 10^{-8}$ | $1.1\times 10^{-5}$
$a_{c}$ | $4.6\times 10^{-5}$ | $4.2\times 10^{-5}$ | $0.0007$
$a_{t}$ | $0.017$ | $0.017$ | $0.22$
$g_{u1}$ | $0.00066+0.0016i$ | $-0.00013+0.013i$ | $0.57-0.35i$
$g_{u2}$ | $0.053+0.30i$ | $0.094+0.24i$ | $-0.060-0.41i$
$g_{u3}$ | $0.11+0.0061i$ | $0.091+0.018i$ | $0.046+0.021i$
$a_{d}$ | $0.00018$ | $0.00014$ | $0.00047$
$a_{s}$ | $1.2\times 10^{-5}$ | $8.7\times 10^{-6}$ | $5.2\times 10^{-5}$
$a_{b}$ | $0.011$ | $0.011$ | $0.025$
$|V_{us}|$ | $0.23$ | $0.22$ | $0.22$
$|V_{cb}|$ | $0.033$ | $0.027$ | $0.042$
$|V_{ub}|$ | $0.0031$ | $0.0020$ | $0.0039$
$\delta_{CP}$ | $58^{\circ}$ | $51^{\circ}$ | $83^{\circ}$
Table 4: Numerical values of parameters and observables at the best fit points
of NH and IH.

###  III.2 IH

In case of IH, we obtain less allowed parameter points compared to the case of
$NH$ since it is more difficult to fit the data. Since there are no points
within 3$\sigma$ region but few points within 5$\sigma$ region, we will
explain the tendency instead of showing scattering plots. The value of $\tau$
is interestingly localized at nearby two fixed points $i$ and $1.74i$, each of
which has remnant symmetry of $Z_{2}$ and $Z_{3}$. $\sum m_{i}$ is localized
at $0.10-0.12$ eV while $\delta_{CP}$ is allowed for the range of
$150^{\circ}-360^{\circ}$. Moreover, almost all the points are within the
cosmological constraint $\sim 0.12$ eV that is similar to the NH case.
$\langle m_{ee}\rangle$ is localized at around $0.016-0.024$ eV while $m_{1}$
is allowed up to $1.2\times 10^{-4}$ eV. Moreover, $m_{1}$ is also localized
at around $10^{-6}$ eV . $\alpha_{21}$ is localized at around 180∘,while
$\alpha_{21}$ is allowed for the range of $100^{\circ}-360^{\circ}$.

In addition to the lepton sector, we discuss our allowed region of quark
sector. Even though the allowed points are not so many, we might say something
from our analysis as follows. As for the CP phase of quark $\delta$ in term of
(1, 3) component of CKM matrix; $|V_{ub}|$, we found whole the region is
allowed at 3$\sigma$ interval. As for $|V_{ub}|$ and $|V_{td}|$, we find that
there is a weak linearly correlation between them. As for $|V_{cb}|$ and
$|V_{ub}|$, we find that there is also a weak linear correlation between them.

Bench mark point for IH: We give two interesting benchmark points;
$\tau\approx 1.06i,\ 1.76i$ to satisfy the quark and lepton masses and mixings
as well as phases in the center and right sides of Tables 3 and 4. The lepton
and neutrino mixings are given by

$\displaystyle\tau\approx 1.06i:$ $\displaystyle V_{\ell_{L}}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.65+0.0068i&0.72+0.00024i&-0.25+0.00061i\\\
-0.47+0.0067i&-0.64+0.0012i&-0.61+0.00072i\\\
-0.60+0.0068i&-0.28+0.00098i&0.75+0.000042i\\\ \end{array}\right],$ (III.21)
$\displaystyle V_{\nu}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.63+0.12i&0.053-0.029i&-0.13-0.75i\\\
0.090+0.11i&0.14+0.98i&0.015-0.089i\\\ -0.75+0.10i&0.12+0.097i&0.068+0.63i\\\
\end{array}\right],$ (III.25) $\displaystyle\tau\approx 1.76i:$ $\displaystyle
V_{\ell_{L}}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.21-0.000031i&0.80-0.0015i&0.56+0.000052i\\\
0.63-0.00065i&-0.33+0.00026i&0.70+0.00032i\\\
0.75-0.00091i&0.50-0.0010i&-0.44-0.00036i\\\ \end{array}\right],$ (III.30)
$\displaystyle V_{\nu}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.010+0.016i&0.063-0.098i&-0.87+0.48i\\\
0.28-0.037i&-0.77+0.56i&-0.12+0.014i\\\ 0.53+0.80i&0.25+0.13i&0.016+0.0094i\\\
\end{array}\right].$ (III.34)

The quark mixings are given by

$\displaystyle\tau\approx 1.06i:$ $\displaystyle V_{u_{L}}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.59+0.26i&-0.18-0.067i&-0.71-0.23i\\\
-0.55+0.21i&-0.50-0.030i&0.61+0.18i\\\ -0.41+0.28i&0.84-0.062i&0.20+0.080i\\\
\end{array}\right],$ (III.38) $\displaystyle V_{d_{L}}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.60+0.0068i&-0.28+0.00096i&-0.75-0.000044i\\\
-0.47+0.0067i&-0.64+0.0011i&0.61-0.00072i\\\
-0.65+0.0068i&0.72+0.00022i&0.25-0.00061i\\\ \end{array}\right],$ (III.42)
$\displaystyle\tau\approx 1.76i:$ $\displaystyle V_{u_{L}}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.76-0.042i&0.43+0.199&-0.45+0.0188i\\\
-0.62+0.015i&-0.33-0.052i&0.71-0.026i\\\ 0.19-0.054i&0.78+0.23i&0.54-0.014i\\\
\end{array}\right],$ (III.46) $\displaystyle V_{d_{L}}$
$\displaystyle=\left[\begin{array}[]{ccc}-0.64+0.00066i&-0.63+0.0011i&-0.44-0.00039i\\\
-0.68+0.00064i&0.21-0.00013i&0.70+0.00037i\\\
0.35-0.00023i&-0.75+0.0015i&0.56+0.000089i\\\ \end{array}\right].$ (III.50)

##  IV Conclusions

We have proposed a LQ model to explain the masses and mixings for quark and
lepton, introducing modular $A_{4}$ symmetry. Due to nature of LQ model that
lepton(quark) directly connects to the quark(lepton) via LQ, a single modulus
number has to be applied that leads to a good motivation towards unification
of quark and lepton flavor in $A_{4}$ modular symmetry. After giving an
assignment for quark sector to reproduce the experimental results at 3$\sigma$
interval, we have also constructed the lepton sector, where the neutrino mass
matrix is induced at one-loop level running down quark sector, unified value
of $\tau$ is used for quark and lepton.

Then, we have performed numerical analysis to search for allowed region
satisfying experimental measurements for both quark and lepton sector,
depending on NH and IH. In case of NH, we have found rather wide allowed space
within 3$\sigma$ interval and obtained tendency of observables for quark and
lepton. Especially, we have found allowed region at nearby $\tau=1.75i$ that
is close to the fixed point of $\tau=i\infty$. Thus, we have also shown a
promising bench mark point at around the solution.

In case of IH, we would not found the allowed region within 3$\sigma$
interval, but found within 5$\sigma$ interval. Although the number of allowed
point is few, we have found all the allowed regions are localized at nearby
$\tau=i,\ 1.76i$, both of which are nearby fixed points. We have shown them as
benchmark points. These would be tested near future.

## Acknowledgments

This research was supported by an appointment to the JRG Program at the APCTP
through the Science and Technology Promotion Fund and Lottery Fund of the
Korean Government. This was also supported by the Korean Local Governments -
Gyeongsangbuk-do Province and Pohang City (H.O.), European Regional
Development Fund-Project Engineering Applications of Microworld Physics (Grant
No. CZ.02.1.01/0.0/0.0/16_019/0000766) (Y.O.). H. O. is sincerely grateful for
the KIAS member.

## Appendix

The modular forms of weight 2, $Y^{(2)}_{\bf 3}=[y_{1},y_{2},y_{3}]^{T}$,
transforming as a triplet of $A_{4}$ is written in terms of Dedekind eta-
function $\eta(\tau)$ and its derivative:

$\displaystyle y_{1}(\tau)$ $\displaystyle=$
$\displaystyle\frac{i}{2\pi}\left(\frac{\eta^{\prime}(\tau/3)}{\eta(\tau/3)}+\frac{\eta^{\prime}((\tau+1)/3)}{\eta((\tau+1)/3)}+\frac{\eta^{\prime}((\tau+2)/3)}{\eta((\tau+2)/3)}-\frac{27\eta^{\prime}(3\tau)}{\eta(3\tau)}\right),$
$\displaystyle y_{2}(\tau)$ $\displaystyle=$
$\displaystyle\frac{-i}{\pi}\left(\frac{\eta^{\prime}(\tau/3)}{\eta(\tau/3)}+\omega^{2}\frac{\eta^{\prime}((\tau+1)/3)}{\eta((\tau+1)/3)}+\omega\frac{\eta^{\prime}((\tau+2)/3)}{\eta((\tau+2)/3)}\right),$
(IV.1) $\displaystyle y_{3}(\tau)$ $\displaystyle=$
$\displaystyle\frac{-i}{\pi}\left(\frac{\eta^{\prime}(\tau/3)}{\eta(\tau/3)}+\omega\frac{\eta^{\prime}((\tau+1)/3)}{\eta((\tau+1)/3)}+\omega^{2}\frac{\eta^{\prime}((\tau+2)/3)}{\eta((\tau+2)/3)}\right),$
$\displaystyle\eta(\tau)$ $\displaystyle=$ $\displaystyle
q^{1/24}\Pi_{n=1}^{\infty}(1-q^{n}),\quad q=e^{2\pi i\tau},\quad\omega=e^{2\pi
i/3}$

Then, any multiplets of higher weight are constructed by multiplication rules
of $A_{4}$, and one finds the following :

$\displaystyle Y^{(4)}_{\bf 1}=y^{2}_{1}+2y_{2}y_{3},\quad Y^{(4)}_{\bf
3}\equiv\left[\begin{array}[]{c}y^{\prime}_{1}\\\ y^{\prime}_{2}\\\
y^{\prime}_{3}\\\
\end{array}\right]=\left[\begin{array}[]{c}y^{2}_{1}-y_{2}y_{3}\\\
y^{2}_{3}-y_{1}y_{2}\\\ y^{2}_{2}-y_{1}y_{3}\\\ \end{array}\right],\quad
Y^{(6)}_{\bf 1}=y_{1}^{2}+y_{2}^{2}+y_{3}^{2}-3y_{1}y_{2}y_{3},$ (IV.8)
$\displaystyle Y^{(6)}_{\bf 3}\equiv\left[\begin{array}[]{c}f_{1}\\\ f_{2}\\\
f_{3}\\\
\end{array}\right]=\left[\begin{array}[]{c}y^{3}_{1}+2y_{1}y_{2}y_{3}\\\
y^{2}_{1}y_{2}+2y_{2}^{2}y_{3}\\\ y^{2}_{1}y_{3}+2y_{3}^{2}y_{2}\\\
\end{array}\right],\quad Y^{(6)}_{\bf
3^{\prime}}\equiv\left[\begin{array}[]{c}f^{\prime}_{1}\\\ f^{\prime}_{2}\\\
f^{\prime}_{3}\\\
\end{array}\right]=\left[\begin{array}[]{c}y^{3}_{3}+2y_{1}y_{2}y_{3}\\\
y^{2}_{3}y_{1}+2y_{1}^{2}y_{2}\\\ y^{2}_{3}y_{2}+2y_{2}^{2}y_{1}\\\
\end{array}\right].$ (IV.21)

## References

  * (1) M. Bauer and M. Neubert, Phys. Rev. Lett. 116, no. 14, 141802 (2016) doi:10.1103/PhysRevLett.116.141802 [arXiv:1511.01900 [hep-ph]].
  * (2) C. H. Chen, T. Nomura and H. Okada, Phys. Rev. D 94 (2016) no.11, 115005 doi:10.1103/PhysRevD.94.115005 [arXiv:1607.04857 [hep-ph]].
  * (3) E. Coluccio Leskow, G. D’Ambrosio, A. Crivellin and D. Müller, Phys. Rev. D 95, no.5, 055018 (2017) doi:10.1103/PhysRevD.95.055018 [arXiv:1612.06858 [hep-ph]].
  * (4) C. H. Chen, T. Nomura and H. Okada, Phys. Lett. B 774 (2017), 456-464 doi:10.1016/j.physletb.2017.10.005 [arXiv:1703.03251 [hep-ph]].
  * (5) T. Nomura and H. Okada, [arXiv:2104.03248 [hep-ph]].
  * (6) S. Sahoo and R. Mohanta, Phys. Rev. D 91, no. 9, 094019 (2015) doi:10.1103/PhysRevD.91.094019 [arXiv:1501.05193 [hep-ph]].
  * (7) D. Becirevic, S. Fajfer, N. Kosnik, and O. Sumensari, Phys. Rev. D 94, no. 11, 115021 (2016) doi:10.1103/PhysRevD.94.115021 [arXiv:1608.08501 [hep-ph]].
  * (8) Y. Cai, J. Gargalionis, M. A. Schmidt and R. R. Volkas, JHEP 10 (2017), 047 doi:10.1007/JHEP10(2017)047 [arXiv:1704.05849 [hep-ph]].
  * (9) Y. Sakaki, M. Tanaka, A. Tayduganov, and R. Watanabe, Phys. Rev. D 88, no. 9, 094012 (2013) doi:10.1103/PhysRevD.88.094012 [arXiv:1309.0301 [hep-ph]].
  * (10) D. Aristizabal Sierra, M. Hirsch and S. G. Kovalenko, Phys. Rev. D 77 (2008), 055011 doi:10.1103/PhysRevD.77.055011 [arXiv:0710.5699 [hep-ph]].
  * (11) K. Cheung, T. Nomura and H. Okada, Phys. Rev. D 94, no.11, 115024 (2016) doi:10.1103/PhysRevD.94.115024 [arXiv:1610.02322 [hep-ph]].
  * (12) M. Huschle et al. [Belle Collaboration], Phys. Rev. D 92, no. 7, 072014 (2015) doi:10.1103/PhysRevD.92.072014 [arXiv:1507.03233 [hep-ex]].
  * (13) J. P. Lees et al. [BaBar Collaboration], Phys. Rev. Lett. 109, 101802 (2012) doi:10.1103/PhysRevLett.109.101802 [arXiv:1205.5442 [hep-ex]].
  * (14) J. P. Lees et al. [BaBar Collaboration], Phys. Rev. D 88, no. 7, 072012 (2013) doi:10.1103/PhysRevD.88.072012 [arXiv:1303.0571 [hep-ex]].
  * (15) A. Abdesselam et al. [Belle Collaboration], arXiv:1603.06711 [hep-ex].
  * (16) S. Hirose et al. [Belle], Phys. Rev. Lett. 118, no.21, 211801 (2017) doi:10.1103/PhysRevLett.118.211801 [arXiv:1612.00529 [hep-ex]].
  * (17) R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 115, no. 11, 111803 (2015) Addendum: [Phys. Rev. Lett. 115, no. 15, 159901 (2015)] doi:10.1103/PhysRevLett.115.159901, 10.1103/PhysRevLett.115.111803 [arXiv:1506.08614 [hep-ex]].
  * (18) R. Aaij et al. [LHCb], Phys. Rev. D 97, no.7, 072013 (2018) doi:10.1103/PhysRevD.97.072013 [arXiv:1711.02505 [hep-ex]].
  * (19) B. Abi et al. [Muon g-2], Phys. Rev. Lett. 126 (2021) no.14, 141801 doi:10.1103/PhysRevLett.126.141801 [arXiv:2104.03281 [hep-ex]].
  * (20) R. Aaij et al. [LHCb], [arXiv:2103.11769 [hep-ex]].
  * (21) F. Feruglio, doi:10.1142/9789813238053_0012 [arXiv:1706.08749 [hep-ph]].
  * (22) R. de Adelhart Toorop, F. Feruglio and C. Hagedorn, Nucl. Phys. B 858, 437-467 (2012) doi:10.1016/j.nuclphysb.2012.01.017 [arXiv:1112.1340 [hep-ph]].
  * (23) J. C. Criado and F. Feruglio, SciPost Phys. 5, no.5, 042 (2018) doi:10.21468/SciPostPhys.5.5.042 [arXiv:1807.01125 [hep-ph]].
  * (24) T. Kobayashi, N. Omoto, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, JHEP 11, 196 (2018) doi:10.1007/JHEP11(2018)196 [arXiv:1808.03012 [hep-ph]].
  * (25) H. Okada and M. Tanimoto, Phys. Lett. B 791, 54-61 (2019) doi:10.1016/j.physletb.2019.02.028 [arXiv:1812.09677 [hep-ph]].
  * (26) T. Nomura and H. Okada, Phys. Lett. B 797, 134799 (2019) doi:10.1016/j.physletb.2019.134799 [arXiv:1904.03937 [hep-ph]].
  * (27) H. Okada and M. Tanimoto, Eur. Phys. J. C 81, no.1, 52 (2021) doi:10.1140/epjc/s10052-021-08845-y [arXiv:1905.13421 [hep-ph]].
  * (28) F. J. de Anda, S. F. King and E. Perdomo, Phys. Rev. D 101, no.1, 015028 (2020) doi:10.1103/PhysRevD.101.015028 [arXiv:1812.05620 [hep-ph]].
  * (29) P. P. Novichkov, S. T. Petcov and M. Tanimoto, Phys. Lett. B 793, 247-258 (2019) doi:10.1016/j.physletb.2019.04.043 [arXiv:1812.11289 [hep-ph]].
  * (30) T. Nomura and H. Okada, Nucl. Phys. B 966, 115372 (2021) doi:10.1016/j.nuclphysb.2021.115372 [arXiv:1906.03927 [hep-ph]].
  * (31) H. Okada and Y. Orikasa, [arXiv:1907.13520 [hep-ph]].
  * (32) G. J. Ding, S. F. King and X. G. Liu, JHEP 09, 074 (2019) doi:10.1007/JHEP09(2019)074 [arXiv:1907.11714 [hep-ph]].
  * (33) T. Nomura, H. Okada and O. Popov, Phys. Lett. B 803, 135294 (2020) doi:10.1016/j.physletb.2020.135294 [arXiv:1908.07457 [hep-ph]].
  * (34) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, Phys. Rev. D 100, no.11, 115045 (2019) [erratum: Phys. Rev. D 101, no.3, 039904 (2020)] doi:10.1103/PhysRevD.100.115045 [arXiv:1909.05139 [hep-ph]].
  * (35) T. Asaka, Y. Heo, T. H. Tatsuishi and T. Yoshida, JHEP 01, 144 (2020) doi:10.1007/JHEP01(2020)144 [arXiv:1909.06520 [hep-ph]].
  * (36) D. Zhang, Nucl. Phys. B 952, 114935 (2020) doi:10.1016/j.nuclphysb.2020.114935 [arXiv:1910.07869 [hep-ph]].
  * (37) G. J. Ding, S. F. King, X. G. Liu and J. N. Lu, JHEP 12, 030 (2019) doi:10.1007/JHEP12(2019)030 [arXiv:1910.03460 [hep-ph]].
  * (38) T. Kobayashi, T. Nomura and T. Shimomura, Phys. Rev. D 102, no.3, 035019 (2020) doi:10.1103/PhysRevD.102.035019 [arXiv:1912.00637 [hep-ph]].
  * (39) T. Nomura, H. Okada and S. Patra, Nucl. Phys. B 967, 115395 (2021) doi:10.1016/j.nuclphysb.2021.115395 [arXiv:1912.00379 [hep-ph]].
  * (40) X. Wang, Nucl. Phys. B 957, 115105 (2020) doi:10.1016/j.nuclphysb.2020.115105 [arXiv:1912.13284 [hep-ph]].
  * (41) H. Okada and Y. Shoji, Nucl. Phys. B 961, 115216 (2020) doi:10.1016/j.nuclphysb.2020.115216 [arXiv:2003.13219 [hep-ph]].
  * (42) H. Okada and M. Tanimoto, [arXiv:2005.00775 [hep-ph]].
  * (43) M. K. Behera, S. Singirala, S. Mishra and R. Mohanta, [arXiv:2009.01806 [hep-ph]].
  * (44) M. K. Behera, S. Mishra, S. Singirala and R. Mohanta, [arXiv:2007.00545 [hep-ph]].
  * (45) T. Nomura and H. Okada, [arXiv:2007.04801 [hep-ph]].
  * (46) T. Nomura and H. Okada, [arXiv:2007.15459 [hep-ph]].
  * (47) T. Asaka, Y. Heo and T. Yoshida, Phys. Lett. B 811, 135956 (2020) doi:10.1016/j.physletb.2020.135956 [arXiv:2009.12120 [hep-ph]].
  * (48) H. Okada and M. Tanimoto, Phys. Rev. D 103, no.1, 015005 (2021) doi:10.1103/PhysRevD.103.015005 [arXiv:2009.14242 [hep-ph]].
  * (49) K. I. Nagao and H. Okada, [arXiv:2010.03348 [hep-ph]].
  * (50) H. Okada and M. Tanimoto, JHEP 03, 010 (2021) doi:10.1007/JHEP03(2021)010 [arXiv:2012.01688 [hep-ph]].
  * (51) C. Y. Yao, J. N. Lu and G. J. Ding, JHEP 05 (2021), 102 doi:10.1007/JHEP05(2021)102 [arXiv:2012.13390 [hep-ph]].
  * (52) P. Chen, G. J. Ding and S. F. King, JHEP 04 (2021), 239 doi:10.1007/JHEP04(2021)239 [arXiv:2101.12724 [hep-ph]].
  * (53) T. Kobayashi, T. Shimomura and M. Tanimoto, [arXiv:2102.10425 [hep-ph]].
  * (54) M. Kashav and S. Verma, [arXiv:2103.07207 [hep-ph]].
  * (55) H. Okada, Y. Shimizu, M. Tanimoto and T. Yoshida, [arXiv:2105.14292 [hep-ph]].
  * (56) T. Kobayashi, K. Tanaka and T. H. Tatsuishi, Phys. Rev. D 98, no.1, 016004 (2018) doi:10.1103/PhysRevD.98.016004 [arXiv:1803.10391 [hep-ph]].
  * (57) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto, T. H. Tatsuishi and H. Uchida, Phys. Lett. B 794, 114-121 (2019) doi:10.1016/j.physletb.2019.05.034 [arXiv:1812.11072 [hep-ph]].
  * (58) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, PTEP 2020, no.5, 053B05 (2020) doi:10.1093/ptep/ptaa055 [arXiv:1906.10341 [hep-ph]].
  * (59) H. Okada and Y. Orikasa, Phys. Rev. D 100, no.11, 115037 (2019) doi:10.1103/PhysRevD.100.115037 [arXiv:1907.04716 [hep-ph]].
  * (60) S. Mishra, [arXiv:2008.02095 [hep-ph]].
  * (61) X. Du and F. Wang, JHEP 02, 221 (2021) doi:10.1007/JHEP02(2021)221 [arXiv:2012.01397 [hep-ph]].
  * (62) J. T. Penedo and S. T. Petcov, Nucl. Phys. B 939, 292-307 (2019) doi:10.1016/j.nuclphysb.2018.12.016 [arXiv:1806.11040 [hep-ph]].
  * (63) P. P. Novichkov, J. T. Penedo, S. T. Petcov and A. V. Titov, JHEP 04, 005 (2019) doi:10.1007/JHEP04(2019)005 [arXiv:1811.04933 [hep-ph]].
  * (64) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, JHEP 02, 097 (2020) doi:10.1007/JHEP02(2020)097 [arXiv:1907.09141 [hep-ph]].
  * (65) S. F. King and Y. L. Zhou, Phys. Rev. D 101, no.1, 015001 (2020) doi:10.1103/PhysRevD.101.015001 [arXiv:1908.02770 [hep-ph]].
  * (66) H. Okada and Y. Orikasa, [arXiv:1908.08409 [hep-ph]].
  * (67) J. C. Criado, F. Feruglio and S. J. D. King, JHEP 02, 001 (2020) doi:10.1007/JHEP02(2020)001 [arXiv:1908.11867 [hep-ph]].
  * (68) X. Wang and S. Zhou, JHEP 05, 017 (2020) doi:10.1007/JHEP05(2020)017 [arXiv:1910.09473 [hep-ph]].
  * (69) Y. Zhao and H. H. Zhang, JHEP 03 (2021), 002 doi:10.1007/JHEP03(2021)002 [arXiv:2101.02266 [hep-ph]].
  * (70) S. F. King and Y. L. Zhou, JHEP 04 (2021), 291 doi:10.1007/JHEP04(2021)291 [arXiv:2103.02633 [hep-ph]].
  * (71) G. J. Ding, S. F. King and C. Y. Yao, [arXiv:2103.16311 [hep-ph]].
  * (72) X. Zhang and S. Zhou, [arXiv:2106.03433 [hep-ph]].
  * (73) Bu-Yao Qu, Xiang-Gan Liu, Ping-Tao Chen, Gui-Jun Ding [arXiv:2106.11659 [hep-ph]].
  * (74) P. P. Novichkov, J. T. Penedo, S. T. Petcov and A. V. Titov, JHEP 04, 174 (2019) doi:10.1007/JHEP04(2019)174 [arXiv:1812.02158 [hep-ph]].
  * (75) G. J. Ding, S. F. King and X. G. Liu, Phys. Rev. D 100, no.11, 115005 (2019) doi:10.1103/PhysRevD.100.115005 [arXiv:1903.12588 [hep-ph]].
  * (76) X. Wang, B. Yu and S. Zhou, Phys. Rev. D 103, no.7, 076005 (2021) doi:10.1103/PhysRevD.103.076005 [arXiv:2010.10159 [hep-ph]].
  * (77) C. Y. Yao, X. G. Liu and G. J. Ding, Phys. Rev. D 103, no.9, 095013 (2021) doi:10.1103/PhysRevD.103.095013 [arXiv:2011.03501 [hep-ph]].
  * (78) X. Wang and S. Zhou, [arXiv:2102.04358 [hep-ph]].
  * (79) A. Baur, H. P. Nilles, A. Trautner and P. K. S. Vaudrevange, Phys. Lett. B 795, 7-14 (2019) doi:10.1016/j.physletb.2019.03.066 [arXiv:1901.03251 [hep-th]].
  * (80) I. de Medeiros Varzielas, S. F. King and Y. L. Zhou, Phys. Rev. D 101, no.5, 055033 (2020) doi:10.1103/PhysRevD.101.055033 [arXiv:1906.02208 [hep-ph]].
  * (81) X. G. Liu and G. J. Ding, JHEP 08, 134 (2019) doi:10.1007/JHEP08(2019)134 [arXiv:1907.01488 [hep-ph]].
  * (82) P. Chen, G. J. Ding, J. N. Lu and J. W. F. Valle, Phys. Rev. D 102, no.9, 095014 (2020) doi:10.1103/PhysRevD.102.095014 [arXiv:2003.02734 [hep-ph]].
  * (83) P. P. Novichkov, J. T. Penedo and S. T. Petcov, Nucl. Phys. B 963, 115301 (2021) doi:10.1016/j.nuclphysb.2020.115301 [arXiv:2006.03058 [hep-ph]].
  * (84) X. G. Liu, C. Y. Yao and G. J. Ding, Phys. Rev. D 103, no.5, 056013 (2021) doi:10.1103/PhysRevD.103.056013 [arXiv:2006.10722 [hep-ph]].
  * (85) S. Kikuchi, T. Kobayashi, H. Otsuka, S. Takada and H. Uchida, JHEP 11, 101 (2020) doi:10.1007/JHEP11(2020)101 [arXiv:2007.06188 [hep-th]].
  * (86) Y. Almumin, M. C. Chen, V. Knapp-Pérez, S. Ramos-Sánchez, M. Ratz and S. Shukla, JHEP 05 (2021), 078 doi:10.1007/JHEP05(2021)078 [arXiv:2102.11286 [hep-th]].
  * (87) G. J. Ding, F. Feruglio and X. G. Liu, SciPost Phys. 10 (2021), 133 doi:10.21468/SciPostPhys.10.6.133 [arXiv:2102.06716 [hep-ph]].
  * (88) F. Feruglio, V. Gherardi, A. Romanino and A. Titov, JHEP 05 (2021), 242 doi:10.1007/JHEP05(2021)242 [arXiv:2101.08718 [hep-ph]].
  * (89) S. Kikuchi, T. Kobayashi and H. Uchida, [arXiv:2101.00826 [hep-th]].
  * (90) P. P. Novichkov, J. T. Penedo and S. T. Petcov, JHEP 04 (2021), 206 doi:10.1007/JHEP04(2021)206 [arXiv:2102.07488 [hep-ph]].
  * (91) G. Altarelli and F. Feruglio, Rev. Mod. Phys. 82, 2701-2729 (2010) doi:10.1103/RevModPhys.82.2701 [arXiv:1002.0211 [hep-ph]].
  * (92) H. Ishimori, T. Kobayashi, H. Ohki, Y. Shimizu, H. Okada and M. Tanimoto, Prog. Theor. Phys. Suppl. 183, 1-163 (2010) doi:10.1143/PTPS.183.1 [arXiv:1003.3552 [hep-th]].
  * (93) H. Ishimori, T. Kobayashi, H. Ohki, H. Okada, Y. Shimizu and M. Tanimoto, Lect. Notes Phys. 858, 1-227 (2012) doi:10.1007/978-3-642-30805-5
  * (94) D. Hernandez and A. Y. Smirnov, Phys. Rev. D 86, 053014 (2012) doi:10.1103/PhysRevD.86.053014 [arXiv:1204.0445 [hep-ph]].
  * (95) S. F. King and C. Luhn, Rept. Prog. Phys. 76, 056201 (2013) doi:10.1088/0034-4885/76/5/056201 [arXiv:1301.1340 [hep-ph]].
  * (96) S. F. King, A. Merle, S. Morisi, Y. Shimizu and M. Tanimoto, New J. Phys. 16, 045018 (2014) doi:10.1088/1367-2630/16/4/045018 [arXiv:1402.4271 [hep-ph]].
  * (97) S. F. King, Prog. Part. Nucl. Phys. 94, 217-256 (2017) doi:10.1016/j.ppnp.2017.01.003 [arXiv:1701.04413 [hep-ph]].
  * (98) S. T. Petcov, Eur. Phys. J. C 78, no.9, 709 (2018) doi:10.1140/epjc/s10052-018-6158-5 [arXiv:1711.10806 [hep-ph]].
  * (99) A. Baur, H. P. Nilles, A. Trautner and P. K. S. Vaudrevange, Nucl. Phys. B 947, 114737 (2019) doi:10.1016/j.nuclphysb.2019.114737 [arXiv:1908.00805 [hep-th]].
  * (100) T. Kobayashi, Y. Shimizu, K. Takagi, M. Tanimoto, T. H. Tatsuishi and H. Uchida, Phys. Rev. D 101, no.5, 055046 (2020) doi:10.1103/PhysRevD.101.055046 [arXiv:1910.11553 [hep-ph]].
  * (101) P. P. Novichkov, J. T. Penedo, S. T. Petcov and A. V. Titov, JHEP 07, 165 (2019) doi:10.1007/JHEP07(2019)165 [arXiv:1905.11970 [hep-ph]].
  * (102) T. Kobayashi, T. Shimomura and M. Tanimoto, [arXiv:2102.10425 [hep-ph]].
  * (103) M. Tanimoto and K. Yamamoto, [arXiv:2106.10919 [hep-ph]].
  * (104) M. C. Chen, S. Ramos-Sánchez and M. Ratz, Phys. Lett. B 801, 135153 (2020) doi:10.1016/j.physletb.2019.135153 [arXiv:1909.06910 [hep-ph]].
  * (105) I. de Medeiros Varzielas, M. Levy and Y. L. Zhou, JHEP 11, 085 (2020) doi:10.1007/JHEP11(2020)085 [arXiv:2008.05329 [hep-ph]].
  * (106) K. Ishiguro, T. Kobayashi and H. Otsuka, JHEP 03, 161 (2021) doi:10.1007/JHEP03(2021)161 [arXiv:2011.09154 [hep-ph]].
  * (107) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, JHEP 09, 178 (2020) doi:10.1007/JHEP09(2020)178 [arXiv:2007.14792 [hep-ph]].
  * (108) A. Gando et al. [KamLAND-Zen], Phys. Rev. Lett. 117, no.8, 082503 (2016) doi:10.1103/PhysRevLett.117.082503 [arXiv:1605.02889 [hep-ex]].
  * (109) Z. Maki, M. Nakagawa and S. Sakata, Prog. Theor. Phys. 28, 870-880 (1962) doi:10.1143/PTP.28.870
  * (110) P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020).


<Ref>(a) suggests that any value $a \gtrsim 0.85$ works well for the random BMC considered in these experiments.
Fixing $a=0.9$, <Ref>(b) suggests a strictly decreasing relationship between the relative accuracy and the exponent $b$.
Fixing $b=0.1$, <Ref>(c) exhibits a clear peak of the relative accuracy as we vary the singular value threshold's exponent $c$ at around $0.75$.

§.§.§ Effect of embedding dimension on cluster density

<Ref> shows a plot of the average intracluster distance $\bar{d}_{\text{intra}}$, and the average intercluster distance $\bar{d}_{\text{inter}}$, both as a function of the embedding rank $r$ for some BMC sampled from the same ensemble as was used in <Ref>.
We define these distances as

\begin{equation}
\bar{d}_{\textnormal{intra}}
\eqdef
\frac{1}{K}\sum_{k=1}^K\frac{1}{|\mcV_k|}\sum_{x\in\mcV_k}\lVert \hat{X}_{x,\cdot} - \mu_{k}\rVert_2
\quad
\bar{d}_{\textnormal{inter}}
\eqdef
\frac{1}{2K(K-1)}\sum_{k=1}^K\sum_{k'\neq k}\lVert\hat{\mu}_k - \hat{\mu}_{k'}\rVert_2
\end{equation}

where $\hat{\mu}_k\eqdef (1/|\mcV_k|)\sum_{x\in \mcV_k} \hat{X}_{x,\cdot}$ denote the cluster centroids, and $\lVert\cdot\rVert_2$ is the usual Euclidean norm.
Note that $\bar{d}_{\textnormal{inter}}$ is defined as half the average distance between cluster centroids.

Observe that for each $\beta$, the intra- and intercluster distances both increase as $r$ increases.
However, for $r \gtrsim K$, the $\bar{d}_{\textnormal{inter}}$ plateaus whereas $\bar{d}_{\textnormal{inter}}$ continues to increase.
Also, observe that $\bar{d}_{\textnormal{intra}}$ grows slower for $r \lesssim K$ than for $r \gtrsim K$.

§.§.§ Singular values in sparse regimes

<Ref> shows how the singular values $\{ \sigma_i \}_{i=1, 2, 3}$ of $\hat{N}_{\Gamma}$ depend on the parameters $p_0$ and $\beta$ from <Ref>.

Observe in <Ref> that $\sigma_1$ and $\sigma_3$ are insensitive to $p_0$; while on the contrary, $\sigma_2$ depends approximately linearly on $p_0$.
A brief computation shows that the expectation $N=\mbbE[\hat{N}]$ has nonzero singular values $(n/\ell_n)\sigma_{1,2}(N)=(1,2p_0-1)$.
The proof of <Ref> in <Ref> furthermore shows that $(n/\ell_n)\sigma_{1,2}(\hat{N}_{\Gamma})$ converge in probability to $(n/\ell_n)\sigma_{1,2}(N)=(1,2p_0-1)$ as $n \rightarrow \infty$.
However, to understand this phenomena completely, a more detailed analysis would be needed to describe the concentration of $\sigma_{1,2}(\hat{N}_{\Gamma})$ at finite $n$.

This caveat is underscored by <Ref>.
These figures demonstrate how $\sigma_{2}(\hat{N}_{\Gamma})$ deviates from the asymptotic scaling predicted by (<ref>).
Specifically, since $\gamma_n=(\ell_n/n)^{0.95}=o(\ell_n/n)$, (<ref>) predicts that $\sigma_2(\hat{N}_{\Gamma})/\gamma_n=\omega_{\mbbP}(1)$; but on the contrary, observe that $\sigma_2(\hat{N}_{\Gamma})/\gamma_n$ is a decreasing function of $\ell_n$ and $n$.
The same is to a lesser extent true for $\sigma_1(\hat{N}_{\Gamma})$.

§.§.§ Confidence intervals for Table 1

<Ref> reports the number of repetitions $N_{\text{rep}}$ and the width of the estimated $95\%$-condfidence interval for the sample means $\mu_{\hat{K}}$ reported in <Ref>.
Note that all but three confidence intervals are at most $0.15$.

$l_n$         $K$     Test     2cALG2     2cHDBS              2cMEGH          2cLLSC     2cLLCI     2cCAIC

$N_{\text{rep}}$     MoE     $N_{\text{rep}}$     MoE     $N_{\text{rep}}$     MoE     $N_{\text{rep}}$     MoE     $N_{\text{rep}}$     MoE     $N_{\text{rep}}$     MoE

4*$n\ln(n)^2$     4*$3$
1     $250$     $ 0.00 $     $250$     $ 0.00 $     $250$     $ 0.04 $     $700$     $ 0.15 $     $250$     $ 0.00 $     $350$     $ 0.15 $

2     $300$     $ 0.14 $     $300$     $ 0.14 $     $250$     $ 0.14 $     $1550$     $ 0.15 $     $350$     $ 0.15 $     $350$     $ 0.15 $

3     $250$     $ 0.10 $     $250$     $ 0.07 $     $250$     $ 0.14 $     $900$     $ 0.14 $     $450$     $ 0.15 $     $450$     $ 0.15 $

4     $250$     $ 0.14 $     $250$     $ 0.14 $     $350$     $ 0.14 $     $750$     $ 0.15 $     $250$     $ 0.12 $     $350$     $ 0.14 $

4*$n\ln(n)^2$     4*$6$
1     $250$     $ 0.04 $     $250$     $ 0.03 $     $250$     $ 0.15 $     $4450$     $ 0.15 $     $850$     $ 0.15 $     $250$     $ 0.14 $

2     $450$     $ 0.15 $     $400$     $ 0.14 $     $250$     $ 0.08 $     $2950$     $ 0.15 $     $900$     $ 0.14 $     $500$     $ 0.15 $

3     $250$     $ 0.13 $     $250$     $ 0.13 $     $250$     $ 0.13 $     $1800$     $ 0.15 $     $1250$     $ 0.15 $     $500$     $ 0.14 $

4     $300$     $ 0.14 $     $250$     $ 0.13 $     $250$     $ 0.09 $     $2500$     $ 0.15 $     $850$     $ 0.15 $     $500$     $ 0.15 $

4*$n\ln(n)^2$     4*$10$
1     $250$     $ 0.03 $     $250$     $ 0.01 $     $250$     $ 0.07 $     $1550$     $ 0.15 $     $1850$     $ 0.15 $     $250$     $ 0.08 $

2     $250$     $ 0.11 $     $250$     $ 0.07 $     $250$     $ 0.08 $     $2500$     $ 0.15 $     $2550$     $ 0.15 $     $650$     $ 0.14 $

3     $350$     $ 0.14 $     $300$     $ 0.14 $     $250$     $ 0.14 $     $3150$     $ 0.15 $     $1700$     $ 0.15 $     $950$     $ 0.15 $

4     $250$     $ 0.11 $     $250$     $ 0.10 $     $250$     $ 0.13 $     $4000$     $ 0.15 $     $2500$     $ 0.15 $     $700$     $ 0.15 $

4*$n^2$     4*$3$
1     $250$     $ 0.00 $     $250$     $ 0.01 $     $250$     $ 0.04 $     $250$     $ 0.00 $     $250$     $ 0.00 $     $350$     $ 0.14 $

2     $2600$     $ 0.15 $     $9850$     $ 0.38^* $     $250$     $ 0.04 $     $250$     $ 0.10 $     $250$     $ 0.09 $     $350$     $ 0.15 $

3     $250$     $ 0.07 $     $9650$     $ 0.47^* $     $350$     $ 0.14 $     $250$     $ 0.05 $     $250$     $ 0.01 $     $450$     $ 0.15 $

4     $2450$     $ 0.15 $     $9850$     $ 0.91^* $     $1750$     $ 0.15 $     $550$     $ 0.15 $     $250$     $ 0.04 $     $300$     $ 0.14 $

4*$n^2$     4*$6$
1     $250$     $ 0.01 $     $250$     $ 0.01 $     $250$     $ 0.12 $     $250$     $ 0.01 $     $250$     $ 0.03 $     $250$     $ 0.10 $

2     $2450$     $ 0.15 $     $250$     $ 0.04 $     $250$     $ 0.06 $     $2000$     $ 0.15 $     $350$     $ 0.14 $     $450$     $ 0.14 $

3     $250$     $ 0.11 $     $250$     $ 0.01 $     $250$     $ 0.07 $     $250$     $ 0.01 $     $250$     $ 0.00 $     $250$     $ 0.10 $

4     $2250$     $ 0.15 $     $12000$     $ 0.15 $     $250$     $ 0.08 $     $2350$     $ 0.15 $     $400$     $ 0.14 $     $450$     $ 0.15 $

4*$n^2$     4*$10$
1     $250$     $ 0.00 $     $250$     $ 0.01 $     $250$     $ 0.08 $     $250$     $ 0.13 $     $350$     $ 0.14 $     $250$     $ 0.08 $

2     $3000$     $ 0.15 $     $250$     $ 0.06 $     $250$     $ 0.09 $     $1150$     $ 0.15 $     $2650$     $ 0.15 $     $1200$     $ 0.15 $

3     $250$     $ 0.12 $     $250$     $ 0.00 $     $250$     $ 0.09 $     $450$     $ 0.15 $     $250$     $ 0.04 $     $250$     $ 0.05 $

4     $3150$     $ 0.15 $     $250$     $ 0.11 $     $250$     $ 0.07 $     $1600$     $ 0.15 $     $2450$     $ 0.15 $     $1350$     $ 0.15 $

This table shows the margin of error (MoE) of the estimated $95\%$-confidence intervals for the sample means of the different estimators tested in <Ref>, as well as the number of replications $N_{\text{rep}}$ that were used in each case.

§ FURTHER ALGORITHMIC DETAILS

§.§ Singular value thresholding via Sylvester's law

<Ref> combines Sylvester's law of inertia with a shift property for eigenvalues of matrices to compute the number of singular values of $\hat{N}_\Gamma$ that lie above the threshold $\gamma_n$.

The argument as to why this works is as follows.
Recall first that for $A \in \mathbb{C}^{n \times n}$, $\gamma \in \mathbb{C}$, if $\sigma^2$ is an eigenvalue of $A^\intercal A$, then $\sigma^2 - \gamma^2$ is an eigenvalue of $A^\intercal A - \gamma^2 I_{n \times n}$.
Recall second that Sylvester's law of inertia says that for any nonsingular matrix $S$, the inertia of a symmetric matrix such as $B := \hat{N}_\Gamma^\intercal \hat{N}_\Gamma - \gamma_n^2 I_{n \times n}$ equals the inertia of $C := S B S^\intercal$.
Inertia, here, is defined as the triple of number of positive, zero, and negative eigenvalues.
Note finally that the number of positive eigenvalues of $C$ coincides with the number of singular values of $\hat{N}_\Gamma$ that are larger than $\gamma_n$.
Now, all of this is especially useful to us because the $LDL^\intercal$-decomposition happens to be such a candidate transformation.
Moreover, the inertia of the diagonal matrix in the $LDL^\intercal$-decomposition is easy to determine (one just needs to look at the signs of $D$'s diagonal entries), and the $LDL^\intercal$-decomposition is computationally easier than e.g. computing a full SVD.

§.§ Alternative algorithms, applied off–the–shelve, to try and estimate K

This section present details on alternative methods that may be used to try and detect the number of clusters $K$ of a BMC.
These were studied numerically in <ref>.
Contrary to <ref>, these algorithms do not come with consistency guarantees when applied to a trajectory of a BMC.

§.§.§ Via HDBSCAN

Given the embedding $\hat{X}$, with dimension set by <Ref>, we also use the popular HDBSCAN algorithm to detect the number of clusters via a density-based clustering of the embedding $\hat{X}$.
For our simulations we use the scikit-learn package [52], whose implementation of HDBSCAN is based on [46].
HDBSCAN is an extension of DBSCAN which we briefly describe here.

DBSCAN takes two parameterrs which are analogous to the parameters used in <Ref>.
In particular, DBSCAN first identifies states with dense neighborhoods by calculating the distance to each state's $\rho_n$th nearest neighbor.
If this distance is less than $h_n$, this state is identified as a core point.
All non-core points are labeled as noise.
To highlight the connection with our method, the set of core points can equivalently be defined as

\begin{equation}
\mcC(\rho_n,h_n)\eqdef \{x\in[n]: |\mcN_x|>\rho_n\},
\end{equation}

where $\mcN_x$ is defined as in (<ref>).
A graph is then constructed with vertices $\mcC(\rho_n,h_n)$, and edges between all pairs $x,y\in\mcC(\rho_n,h_n)$ such that $\lVert \hat{X}_{x,\cdot}-\hat{X}_{y,\cdot}\rVert_2 < h_n$.
Clusters are defined as the connected components of this graph.

HDBSCAN extends DBSCAN by performing a sweep over all possible values of $h_n$, and identifying clusters which are stable across the largest range of values of $h_n$.
Such a sweep allows to identify regions of different density, corresponding to clusters which are stable across different ranges of $h_n$.
In our experiments, the remaining parameter $\rho_n$, corresponding to the minimum cluster size, is set to the same value as in <Ref> (see (<ref>)).
We refer to [46] for further details, as well as a statistical interpretation of this method.

Pseudocode for HDBSCAN
A transition count matrix $\hat{N}=(\hat{N}_{x,y})_{x,y\in[n]}$, and neighborhood size threshold $\rho_n$
An estimated number of clusters $\hat{K}_{\mathrm{HDB}}$
Calculate SVD $U\Sigma V^\intercal$ of $\hat{N}_{\Gamma}$;

Construct spectral embedding $\hat{X}\gets [\tilde{U},\tilde{V}]$;

Estimate $\hat{K}_{\mathrm{HDB}}$ using HDBSCAN with data matrix $\hat{X}$ and parameter $\rho_n$;

§.§.§ Via the maximum eigengap

The next method likewise computes eigenvalues of a matrix associated to $\hat{N}$.
A common heuristic for estimating the number of clusters in data is to choose $\hat{K}$ as the location of the maximum gap between successive eigenvalues of some appropriate similarity matrix of the data points.
We implement this idea following the approach of [30].

Define the ingoing and outgoing degrees of state $x\in[n]$ as $d_x^{\mathrm{in}}\eqdef\sum_{y\in[n]}\hat{N}_{y,x}$ and $d_x^{\mathrm{out}}\eqdef \sum_{y\in[n]}\hat{N}_{x,y}$, corresponding to the number of transitions going into and out of state $x$, respectively.
Note that $|d_x^{\mathrm{in}}-d_x^{\mathrm{out}}|\leq 1$ because each transition into a state also has to leave it, except for the first and last time steps.
Then, define the modularity matrix $M$ elementwise as

\begin{equation}
\hat{M}_{x,y}\eqdef\hat{N}_{x,y}-\frac{d_x^{\mathrm{in}} d_y^{\mathrm{out}}}{\ell}.
\end{equation}

Intuitively, the elements of $\hat{M}_{x,y}$ compare the observed transitions between $x$ and $y$ to the expected number of transitions if each transition was sampled indpendently with the probability of observing a transition given by $d_x^{\mathrm{in}}d_y^{\mathrm{out}}/\ell^2$.
The latter gives an expected in-degree equal to $d_x^{\mathrm{in}}$ and an expected out-degree equal to $d_x^{\mathrm{out}}$.

After $\hat{M}$ has been calculated we can define the method to estimate the number of clusters in a BMC.
Let $\lambda_1(\hat{M}),\allowbreak\lambda_2(\hat{M}),\allowbreak\ldots,\allowbreak\lambda_n(\hat{M})\allowbreak\in\mbb{C}$ be the eigenvalues of $\hat{M}$ such that $|\lambda_i(\hat{M})|\geq|\lambda_{i+1}(\hat{M})|$ for all $i\in\{1,\ldots,n-1\}$.
Then we define the estimated number of clusters $\hat{K}_M$ as

\begin{equation}
\hat{K}_M=\argmax_{i=2,\ldots,n}\{|\lambda_{i-1}(\hat{M})|-|\lambda_{i}(\hat{M})|\}
\end{equation}

Pseudocode for maximum eigengap
A transition count matrix $\hat{N}=(\hat{N}_{x,y})_{x,y\in[n]}$
An estimate number of clusters $\hat{K}_M$
Compute in- and out-degrees $d_x^{\mathrm{in}}\gets \sum_{y\in[n]}\hat{N}_{y,x}$, $d_x^{\mathrm{out}}\gets \sum_{y\in[n]}\hat{N}_{x,y}$ for $x\in[n]$;

Build modularity matrix $\hat{M}_{x,y}\gets \hat{N}_{x,y}-d_x^{\mathrm{in}}d_x^{\mathrm{out}}/\ell_n$ for $x,y\in[n]$;

Calculate sorted eigenvalues $|\lambda_1(\hat{M})|,|\lambda_2(\hat{M})|,\ldots,|\lambda_n(\hat{M})|$;

Calculate maximum eigengap $\hat{K}_M\gets \argmax_{i=2,\ldots,n}\{|\lambda_{i-1}(\hat{M})-|\lambda_i(\hat{M})|\}$;

§.§.§ Via a crossvalidated log-likelihood

Stepping away from spectral methods, a natural estimator for the number of clusters is to use a crossvalidated log-likelihood.
We consider BMC with different numbers of clusters as competing model classes and try to select the model which best explains the data.
The procedure is as follows.

First, split the observed sequence $X_0,\ldots,X_{\ell}$ into two parts, $X_{\text{train}}\eqdef X_0,\ldots,X_{\floor{\ell/2}}$, and $X_{\text{val}}\eqdef X_{\floor{\ell/2}+1},\ldots,X_{\ell}$.
For each $k=1,2,\ldots,K_{\max}$, where $K_{\max}\in\mbbN_+$ is some threshold value, we use the training data $X_{\text{train}}$ to estimate a model BMC for the data with $k$ clusters.
Specifically, we use the BMC clustering algorithm of [7] to estimate a cluster assignment.
Let $\hat{Y}^{(k)}_t\eqdef \hat{\sigma}^{(k)}(X_t)$ for $t=0,\ldots,\ell$ denote the sequence of clusters under the estimated cluster assignment $\hat{\sigma}^{(k)}$.
We subsequently compute the maximum-likelihood estimator $\hat{p}^{(k)}=\{\hat{p}^{(k)}_{i,j}\}_{i,j \in [k]}$ of the cluster transition probabilities for the given cluster assignment $\hat{\sigma}^{(k)}$ as

\begin{equation}
\hat{p}^{(k)}_{i,j}\eqdef\frac{\sum_{t=0}^{\floor{\ell/2}-1}\mathbbm{1}\{\hat{Y}^{(k)}_t=i, \hat{Y}^{(k)}_{t+1}=j\}}{\sum_{t=0}^{\floor{\ell/2}-1}\mathbbm{1}\{\hat{Y}^{(k)}_t=i\}}.
\end{equation}

Recall that we consider two estimators for $\sigma$, namely the spectral clustering algorithm <cit.>, and the refined estimate obtained from <cit.>.
The improvement algorithm in particular uses a local optimization of a log-likelihood function, starting from the initial estimate provided by <cit.>.
The combined model $(\hat{\sigma}^{(k)},\hat{p}^{(k)})$ is therefore an approximation of the maximum-likelihood estimator among all BMC with $k$ clusters.

To determine the appropriate value of $k=1,\ldots,K_{\max}$ we estimate for each candidate model $(\hat{\sigma}^{(k)},\hat{p}^{(k)})$ the maximum log-likelihood using the validation data $X_{\text{val}}$ as follows

\begin{equation}
\hat{\mcL}^{(k)}\eqdef \frac{1}{\floor{\ell/2}}\sum_{t=\floor{\ell/2}+1}^{\ell-1}\ln \frac{\hat{p}^{(k)}_{\hat{Y}^{(k)}_t,\hat{Y}^{(k)}_{t+1}}}{|\hat{\mcV}^{(k)}_{\hat{Y}^{(k)}_{t+1}}|},
\end{equation}

Here, $\hat{\mcV}^{(k)}_{i}\eqdef \lb \hat{\sigma}^{(k)}\rb^{-1}(i)\subset[n]$ denotes the $i$'th cluster according to the estimate assignment $\hat{\sigma}^{(k)}$.
We then define the estimator $\hat{K}_{\mcL}$ the value of $k$ that maximizes log-likelihood ratio:

\begin{equation}
\hat{K}_{\mcL}\eqdef\argmax_{k=1,\ldots,K_{\max}}\hat{\mcL}^{(k)}.
\end{equation}

By evaluating $\hat{\mcL}^{(k)}$ on a (nearly) independent set of validation data, we limit the bias of $\hat{K}_{\mcL}$, and reduce overfitting.

Pseudocode for crossvalidated log-likelihood
Training data $X_{\mathrm{train}}=X_0,\ldots,X_{\floor{\ell/2}}$, validation data $X_{\mathrm{val}}=X_{\floor{\ell/2}+1},\ldots,X_{\ell}$, maximum number of clusters $K_{\max}$
Estimated number of clusters $\hat{K}_{\mcL}$
$k\in 1,\ldots,K_{\max}$
Estimate BMC $(\hat{\sigma}^{(k)},\hat{p}^{(k)})$ by running <cit.> with $X_{\mathrm{train}}$;

Estimate cluster sequence $\hat{Y}_t^{(k)}\gets \hat{\sigma}^{(k)}(X_t)$ for $t=\floor{\ell/2}+1,\ldots,\ell$;

Estimate log-likelihood $\hat{\mcL}^{(k)}\gets (1/\floor{\ell/2})\sum_{t=\floor{\ell/2}+1}^{\ell-1}\ln\bigl( \hat{p}^{(k)}_{\hat{Y}^{(k)}_t,\hat{Y}^{(k)}_{t+1}}/|\hat{\mcV}^{(k)}_{\hat{Y}^{(k)}_{t+1}}|\bigr)$;

Select maximum log-likelihood $\hat{K}_{\mcL}\gets \argmax_{k=1,\ldots,K_{\max}}\hat{\mcL}^{(k)}$;

§.§.§ Via a CAIC

The holdout method used to minimize the bias of $\hat{K}_{\mcL}$ means that we can only use half of our data to cluster.
Information criterion methods attempt to avoid this issue by adding a correction term to the objective function, which penalizes model complexity and reduces overfitting.
Here we use the CAIC [48].
Let $(\hat{\sigma}^{(k)},\hat{p}^{(k)})$ be as in <ref>, except that we now use all of the data $X_0,\ldots,X_{\ell}$.
Let also $\hat{Y}_t^{(k)}=\hat{\sigma}^{(k)}(X_t)$.
Then, the CAIC is defined as

\begin{equation}
\text{CAIC}(k)\eqdef -2\sum_{t=0}^{\ell-1}\ln \frac{\hat{p}^{(k)}_{\hat{Y}_t^{(k)},\hat{Y}_{t+1}^{(k)}}}{|\hat{\mcV}^{(k)}_{\hat{Y}_{t+1}^{(k)}}|} + \text{DF}(k)(\ln {\ell} - 1),
\end{equation}

where $\text{DF}(k)$ denotes the number of degrees of freedom of the model with $k$ clusters.
Here, $\text{DF}(k) = DF(\sigma) + \text{DF}(p)=n + k(k-1)$, with $\text{DF}(\sigma)$ the degrees of freedom of the cluster assignment $\sigma$, and $\text{DF}(p)$ the degrees of freedom of the cluster transition matrix $p$.
Note that $p$ has only $k(k-1)$ degrees of freedom, because each row is constrained to sum to one.
We finally define $\hat{K}_{\mathrm{CAIC}}$ as the value that minimizes the CAIC

\begin{equation}
\hat{K}_{\mathrm{CAIC}}\eqdef \argmin_{k=1,\ldots,K_{\max}}\mathrm{CAIC}(k).
\end{equation}

Pseudocode for CAIC
A trajectory $X=X_0,X_1,\ldots,X_{\ell}$, and a maximum number of clusters $K_{\max}$
An estimate number of clusters $\hat{K}_{\mcL}$
$k\in 1,\ldots,K_{\max}$
Estimate BMC $(\hat{\sigma}^{(k)},\hat{p}^{(k)})$ by running <cit.> with $X$;

Estimate cluster sequence $\hat{Y}_t^{(k)}\gets \hat{\sigma}^{(k)}(X_t)$ for $t=0,\ldots,\ell$;

Calculate $\mathrm{CAIC}(k)\gets -2\sum_{t=0}^{\ell-1}\ln \bigl(\hat{p}^{(k)}_{\hat{Y}_t^{(k)},\hat{Y}_{t+1}^{(k)}}/|\hat{\mcV}^{(k)}_{\hat{Y}_{t+1}^{(k)}}|\bigr) + \text{DF}(k)(\ln {\ell} - 1)$;

Select minimum CAIC $\hat{K}_{\mathrm{CAIC}}\gets \argmin_{k=1,\ldots,K_{\max}}\mathrm{CAIC}(k)$;