# A minimal axion model for mass matrices with five texture-zeros

Yithsbey Giraldo<EMAIL_ADDRESS>Departamento de Física, Universidad de
Nariño, A.A. 1175, San Juan de Pasto, Colombia  R. Martinez
<EMAIL_ADDRESS>Departamento de Física, Universidad Nacional de
Colombia, Ciudad Universitaria, K.45 No.26-85, Bogotá D.C., Colombia  Eduardo
Rojas<EMAIL_ADDRESS>Departamento de Física, Universidad de Nariño, A.A.
1175, San Juan de Pasto, Colombia  Juan C. Salazar<EMAIL_ADDRESS>Departamento de Física, Universidad de Nariño, A.A. 1175, San Juan de Pasto,
Colombia

###### Abstract

A model with fermion and scalar fields charged under a Peccei-Queen (PQ)
symmetry is proposed. The PQ charges are chosen in such a way that they can
reproduce mass matrices with five texture zeros, which can generate the
fermion masses, the CKM matrix, and the PMNS matrix of the Standard Model
(SM). To obtain this result, at least 4 Higgs doublets are needed. As we will
see in the manuscript this is a highly non-trivial result since the texture
zeros of the mass matrices impose a large number of restrictions. This model
shows a route to understand the different scales of the SM by extending it
with a multi-Higgs sector and an additional PQ symmetry. Since the PQ charges
are not universal, the model presents flavor-changing neutral currents (FCNC)
at the tree level, a feature that constitutes the main source of restrictions
on the parameter space. We report the allowed regions by lepton decays and
compare them with those coming from the semileptonic decays
$K^{\pm}\longrightarrow\pi\bar{\nu}\nu$. We also show the excluded regions and
the projected bounds of future experiments for the axion-photon coupling as a
function of the axion mass and compare it with the parameter space of our
model.

##  I Introduction

The discovery of the Higgs with a mass of 125 GeV, by the ATLAS Aad:2012tfa
and CMS Chatrchyan:2012ufa  collaborations, is very important because it
provides experimental support for spontaneous symmetry breaking, which is the
mechanism that explains the origin of the masses of fermions and gauge bosons.
Additionally, it opens up the possibility of new physics in the scalar sector,
such as the two Higgs doublet model He:2008qm ; Grzadkowski:2009iz ;
Logan:2010nw ; Boucenna:2011hy ; He:2011gc ; Bai:2012nv ; He:2013suk ;
Cai:2013zga ; Wang:2014elb ; Drozd:2014yla ; Campbell:2015fra ;
vonBuddenbrock:2016rmr ; Carena:2015moc , models with additional singlet
scalar fields Branco:2011iw , or scalar fields that could be candidates for
Dark Matter McDonald:1993ex ; LopezHonorez:2006gr ; Barbieri:2006dq ;
Honorez:2010re ; Carena:2015moc . On the other hand, in the Standard Model
(SM) Glashow:1961tr ; Weinberg:1967tq ; Salam:1968rm , symmetry breaking
generates a coupling of the Higgs to fermions, proportional to their masses,
which is consistent with experimental data. However, there are several orders
of magnitude between the fermion mass hierarchies that cannot be explained
within the context of the SM. Six masses must be defined for the up and down
quarks, three Cabibbo-Kobayashi-Maskawa (CKM) mixing angles, and a complex
phase that involves CP violation. On the other hand, in the lepton sector,
there are three masses for charged leptons, two squared mass differences for
neutrinos, three mixing angles, and a complex phase that involves CP violation
in the lepton sector. In this case, it is necessary to determine the mass of
the lightest neutrino and the character of neutrinos, whether they are Dirac
or Majorana fermions.

In the Davis experiment Davis:1968cp , which was designed to detect solar
neutrinos, a deficiency in the solar neutrino flux was first observed.
According to the results of Bahcall, only one-third of solar neutrinos would
reach the Earth Bahcall:2004cc . Neutrino oscillation was first proposed by
Pontecorvo Pontecorvo:1957qd , and the precise mechanism of solar neutrino
oscillations was proposed by Mikheyev, Smirnov, and Wolfenstein, involving a
resonant enhancement of neutrino oscillations due to matter effects
Wolfenstein:1977ue ; Mikheyev:1985zog . These observations have been confirmed
by many experiments from four different sources: solar neutrinos as in
Homestake Davis:1968cp , SAGEDavis:1968cp , GALLEX & GNOKaether:2010ag ;
Cleveland:1998nv , SNOSNO:2011hxd , BorexinoBorexino:2008fkj ; Bellini2014
and Super-KamiokandeSuper-Kamiokande:2005wtt ; Super-Kamiokande:2006jvq
experiments, atmospheric neutrinos as in IceCubeIceCube:2014flw , neutrinos
from reactors as KamLANDKamLAND:2010fvi , CHOOZCHOOZ:1999hei , Palo
VerdePiepke:2002ju , Daya BayDayaBay:2015ayh , RENOKim:2016yvm  and
SBLKopp:2013vaa , and from accelerators as in MINOSMINOS:2013utc ,
T2KT2K:2013ppw  and NO$\nu$ANOvA:2016kwd . Neutrino oscillations depend on
squared mass differences. On the other hand, the lightest neutrino mass has
not been determined yet, but from cosmological considerations, none of the
neutrino masses can exceed 0.3 eV, which implies that the neutrino masses are
much smaller than the charged fermion masses. However, unlike quarks and
charged leptons, in the SM the neutrinos are massless, which is explained by
assuming that neutrinos are left-handed. Therefore, the discovery of neutrino
masses implies new physics beyond the SM. By adding right-handed neutrinos,
the Higgs mechanism of the SM can give neutrinos the same type of mass
acquired by charged leptons and quarks. It is possible to add right-handed
neutrinos $\nu_{R}$ to the SM, as long as they do not participate in weak
interactions. With the presence of right-handed neutrinos, it would be
possible to generate Dirac masses $m_{D}$, similar to those of charged leptons
and quarks. In principle, it is also possible to give Majorana masses to left-
handed neutrinos, and similarly, right-handed neutrinos can have Majorana
masses $M_{R}$. For a very large $M_{R}$, it would give effective Majorana
masses for left-handed neutrinos as $m_{\text{eff}}\approx m_{D}^{2}/M_{R}$.
The presence of large Majorana masses allows to explain the tiny neutrino
masses compared to the charged fermion masses [38]. To explain the smallness
of neutrino masses, there are three types of seesaw mechanisms in the
literature: type I with three electroweak neutrinos and three heavy right-
handed neutrinos, type II Schechter:1980gr ; Ma:1998dx , type III Foot:1988aq
, and inverse seesaw Mohapatra:1986bd ; Gonzalez-Garcia:1988okv . One way to
explain the fermion mass hierarchies and the CKM and PMNS mixing angles is
through zeros in the Yukawa couplings of fermions (this is known as texture-
zeros or simply textures of the mass matrices, and these zeros are usually
chosen by hand). It is common in the literature to consider Fritzsch-type
textures Branco:2007nb ; Choubey:2008tb , or similar Adhikary:2011pv ;
Adhikary:2010fa ; Branco:2005jr ; Mei:2005qp ; Hagedorn:2004ba ;
Antusch:2005gp , for the neutrino and charged lepton mass matrices.

There is no theory that provides values for the entries of the Yukawa
Lagrangian, and consequently, there is no a first-principle explanation for
the masses and their large differences in the SM. This lack of naturalness is
known as the hierarchy problem Ludl:2014axa ; Ferreira:2014vna ; Ludl:2015lta
; Grimus:2004hf ; Fonseca:2014koa . In this direction, a way that has been
explored in the literature is to propose a sector with multiple scalar
doublets along with discrete symmetries Giraldo:2020hwl ; Bjorkeroth:2018ipq ,
to reduce the number of Yukawa couplings, or equivalently, by introducing
texture-zeros in the mass matrices Cheng:1987rs ; Matsuda:2006xa ;
CarcamoHernandez:2006iuj ; Langacker:2000ju ; Leroux:2001fx ; King:2005jy . It
is also possible to consider global symmetry groups that prohibit certain
Yukawas, which somehow generate the texture-zeros mentioned Ludl:2014axa ;
Ferreira:2014vna ; Ludl:2015lta ; Grimus:2004hf ; Fonseca:2014koa . Another
way of obtaining these textures is through horizontal gauge symmetries, with
the assignment of quantum numbers to the fermion sector, which can break the
universality of the SM Alvarado:2012xi ; CarcamoHernandez:2006iuj ;
Ibanez:1994ig ; Binetruy:1994ru ; Nir:1995bu ; Jain:1994hd ; Dudas:1995yu ;
Pisano:1992bxx ; Frampton:1992wt ; Foot:1994ym ; Long:1995ctv ; Long:1996rfd ;
Long:1997vbr ; Diaz:2002dx ; Diaz:2003dk ; Diaz:2004fs ; Ochoa:2005ch . This
gauge symmetry generates textures that produce flavor changes in the neutral
currents and that, in principle, could be seen in future colliders. There are
models with electroweak extensions of the SM such as $SO(14)$, $SU(9)$,
$3$-$3$-$1$, $U(1)_{X}$, etc. Barr:1979xt ; Taylor:1989tm ; Barr:1989ta ;
Weinberg:1972ws ; Mohapatra:1974wk ; Barr:1976bk ; Balakrishna:1987qd ;
Barr:2001vj ; Ferretti:2006df ; Barr:2008pn ; Frampton:2009ce ; Babu:1995hr ;
Sato:1997hv ; Irges:1998ax ; Barr:2000ka ; Haba:2000be  that attempt to
explain the flavor and the mass hierarchy problem of the SM. Another mechanism
to generate textures in the Yukawa Lagrangian is through additional discrete
or global symmetries. Some groups that have been used in the literature are:
$S_{3}$, $A_{4}$, $\Delta_{27}$, $Z_{2}$, etc. He:2006dk ; Ahn:2013mva ;
GonzalezFelipe:2013xok ; Ishimori:2012fg ; GonzalezCanales:2012blg ;
Mohapatra:2012tb ; Ding:2013hpa ; Ma:2013xqa ; Nishi:2013jqa ;
Altarelli:2005yp ; Ishimori:2010fs ; Kadosh:2010rm ; Ding:2013eca . The
simplest symmetries are of abelian type, which can be used to impose texture-
zeros in the mass matrices to make them predictive. On the other hand, given
fermion mass matrices with texture-zeros, it is possible to find an extended
scalar sector so that the texture-zeros can be generated from abelian
symmetries Giraldo:2020hwl ; Bjorkeroth:2018ipq .

Due to the fact that there are three up-type quarks and three down-type
quarks, the mass operators are $3\times 3$ complex matrices with 36 degrees of
freedom. If we consider these operators to be Hermitian Giraldo:2011ya ;
Branco:1999nb ; Giraldo:2015cpp , the number of free parameters reduces to 18,
which cannot be fully determined from the 10 available physical quantities,
namely masses and mixing angles Ramond:1993kv . This provides freedom to
reduce the number of free parameters in the matrices and search for matrix
structures with zeros that provide eigenvalues and mixing angles consistent
with the masses and mixing matrices of the fermions. One way to find zeros in
the mass matrices that is automatically consistent with experimental data is
based on weak basis transformations (WBT) for quarks and leptons
Fritzsch:1999ee ; Branco:1988iq ; Branco:1999nb ; Giraldo:2015cpp . Fritzsch
proposed an ansatz with six zeros He:1989eh ; Fritzsch:1977vd ;
Fritzsch:1977za ; Fritzsch:1977vd ; Fritzsch:1979zq ; Du:1992iy ;
Fritzsch:2002ga , but the value of $|V_{ub}/V_{cb}|\approx 0.06$ is too small
compared to the experimental value $|V_{ub}/V_{cb}|_{exp}\approx 0.09$
ParticleDataGroup:2022pth . For this reason, the use of 4 and 5 zero-textures
was proposed Fritzsch:2002ga ; Xing:2003yj ; Branco:1999nb ; Zhou:2003ji ;
Giraldo:2011ya ; Benavides:2020pjx ; Giraldo:2015ffa . Reference
Giraldo:2011ya ; Giraldo:2015cpp  showed that matrices with five zero-textures
could reproduce the mass hierarchy and mixing angles of the CKM matrix.

The strong CP problem arises from the fact that the QCD Lagrangian has a non-
perturbative term (”$\theta$-term”) that explicitly violates CP in strong
interactions. Peccei and Quinn proposed a solution to the strong CP problem
Peccei:1977hh ; Peccei:1977ur , where it is assumed that the SM has an
additional global chiral symmetry $U(1)$, which is spontaneously broken at a
large energy scale $f_{a}$. One consequence of this breaking is the existence
of a particle called the axion, which is the Goldstone boson of the broken
$U(1)_{PQ}$ symmetry Weinberg:1977ma ; Wilczek:1977pj . Due to the fact that
the PQ symmetry is not exact at the quantum level, as a result of a chiral
anomaly, the axion is massive and its mass (see Appendix D) is given by:

$m_{a}=\frac{f_{\pi}m_{\pi}}{f_{a}}\frac{\sqrt{z}}{1+z}\approx 6\mu
eV\left(\frac{10^{12}GeV}{f_{a}}\right),$ (1)

where $z=0.56$ is assumed for the up and down quark mass ratio, while
$f_{\pi}\approx 92$ MeV and $m_{\pi}=135$ MeV are the pion decay constant and
mass, respectively.

The effective couplings of axions to ordinary particles are inversely
proportional to $f_{a}$, and also depend on the model. It was originally
thought that the PQ symmetry breaking occurred at the electroweak scale, but
experiments have ruled this out. The mass of the axion and its coupling to
matter and radiation scale as $1/f_{a}$, making its direct detection extremely
difficult. The combined limits from unsuccessful searches in nuclear and
particle physics experiments and from stellar evolution imply that $f_{a}\geq
3\times 10^{9}$ GeV Kim:1979if . Furthermore, there is an upper limit of
$f_{a}\leq 10^{12}$GeV that comes from cosmology, since light axions are
produced in abundance during the QCD phase transition Turner:1989vc ;
Raffelt:1990yz ; Preskill:1982cy ; Abbott:1982af ; Dine:1982ah . Hence, these
models are generically referred to as ”invisible” axion models and remain
phenomenologically viable. There are two classes of invisible axion models in
the literature: KSVZ (Kim, Shifman, Vainshtein, and Zakharov) Kim:1979if ;
Shifman:1979if  and DFSZ (Dine, Fischler, Srednicki, and Zhitnitsky)
Dine:1981rt ; Zhitnitsky:1980tq . The main difference between KSVZ-type and
DFSZ-type axions is that the former do not couple to ordinary quarks and
leptons at tree level, but instead require an exotic quark that ensures a
nonzero QCD anomaly to generate CP violation. Depending on the assumed value
of $f_{a}$, the existence of axions could have interesting consequences in
astrophysics and cosmology. The emission of axions produced in stellar plasma
through their coupling to photons, electrons, and nucleons would provide a new
mechanism for energy loss in stars. This could accelerate the evolutionary
process of stars and, therefore, shorten their lifespan. Axions can also exist
as primordial cosmic relics produced copiously in early times and could be
candidates for dark matter. Of course, if axions are candidates for dark
matter Sikivie:2009qn ; Sikivie:2006ni ; Hannestad:2003ye ; Hannestad:2005df ;
Melchiorri:2007cd ; Hannestad:2007dd ; Hannestad:2008js , they should be very
light, with $m_{a}\leq 10\mu eV$. However, data from numerous laboratory
experiments and astrophysical observations, together with the cosmological
requirement that the contribution to the mass density of the Universe from
relic axions does not saturate the Universe, restricts the allowed values of
the axion mass to a range of $10^{-5}eV<m_{a}<10^{-2}$ eV. One source of
axions would be the Sun, which, coupled to two photons, could be produced
through the Primakoff conversion of thermal photons in the electric and
magnetic fields of the solar plasma. The limits are primarily useful for
complementing the arguments of stellar energy loss Viaux:2013lha  and the
searches for solar axions by CAST at CERN CAST:2008ixs  and the Tokyo axion
helioscope Inoue:2008zp .

The axion-photon coupling (see Appendix D) can be calculated in chiral
perturbation theory as Peccei:1977hh ; Peccei:1977ur .

$g_{a\gamma}=-\frac{\alpha}{2\pi
f_{a}}\left(\frac{E}{N}-\frac{2}{3}\frac{z+4}{z+1}\right).$ (2)

This coupling and the axion mass are related to each other through the
relation $E/N$, which depends on the model and can be tested in experiments.

The strongest limits on the axion-electron coupling are derived from
observations of stars with a dense core, where bremsstrahlung is very
effective. These conditions are realized in White Dwarfs and Red Giant Stars,
where the evolution of a White Dwarf is a cooling process by photon radiation
and neutrino emission, with the possible addition of new energy loss channels
such as axions. Current numerical analysis suggest a limit of $g_{ae}\leq
2.8\times 10^{-13}$ DiLuzio:2020wdo . In particular, using data from the Sloan
Digital Sky Survey (SDSS) and SuperCOSMOS Sky Survey (SCSS)
MillerBertolami:2014rka , they showed that the axion-electron coupling is
approximately $1.4\times 10^{-13}$. More recent analysis of the data in Ref.
MillerBertolami:2014rka  found $g_{ae}=1.5^{+0.6}_{-0.9}\times 10^{-13}$
($95\%$ CL) Giannotti:2017hny ; 10.1093/mnras/sty1162 . The two groups
studying the axion-electron coupling are M5 Viaux:2013lha  and M3
2018arXiv180210357S . Their combination yields the limit
$g_{ae}=1.6^{+0.29}_{-0.34}\times 10^{-13}$. For a recent and comprehensive
review of axion physics, see DiLuzio:2020wdo .

This document is organized as follows: In Section 2, we review the textures
for the quark and lepton mass matrices that will be used in this work. We also
write the real parameters of these matrices in terms of the masses of the SM
fermions and two free parameters. In Section 3, we present the particle
content of our model and the necessary PQ charges to generate the mass matrix
textures presented in Section 2. In Section 4, we adjust the Yukawa couplings
to obtain the masses of the charged leptons and neutrinos. It is important to
note that we cannot use the VEVs to adjust the lepton masses, as these were
already adjusted to reproduce the quark masses. It is also important to note
that by using a seesaw mechanism, we can avoid adjusting the Yukawas, however,
that is not our purpose in the present work. In Section 5, we show the
Lagrangian of our model. In Section 6, we present some constraints in the
parameter space, as well as projected constraints for upcoming experimental
results, both for experiments under construction and in the data-taking phase.

##  II The Five texture-zero mass matrices

The reason for dealing with texture zeros in the Standard Model (SM) and its
extensions is to simplify as much as possible the number of free parameters
that allow us to see relationships between masses and mixings present in these
models. The Yukawa Lagrangian is responsible for giving mass to SM fermions
after spontaneous symmetry breaking. A first simplification, without losing
generality, is to consider that the fermion mass matrices are Hermitian, so
the number of free parameters for each sector of quarks and leptons is reduced
to 18, but there is still an excess of parameters to reproduce the
experimental data provided in the literature. Due to the lack of a model to
make predictions, discrete symmetries can be used to prohibit some components
in the Yukawa matrix, generating the so-called texture zeros for the mass
matrices. In many works, instead of proposing discrete symmetries, texture
zeros are proposed as practical and direct alternatives. The advantage of this
approach is that it is possible to choose each mass matrix in an optimal way
for the analytical treatment of the problem, and at the same time adjust the
mixing angles and the masses of the fermions.

###  II.1 Quark Sector

We should keep in mind that six-zero textures in the SM have already been
discarded because their predictions are outside the experimental ranges
allowed; but, five-zero textures for quark mass matrices is a viable
possibility Verma:2017ppl ; Xing:2019vks ; Fritzsch:1999ee ; Desai:2000bu ;
Ludl:2015lta ; Ponce:2013nsa . Specifically, we chose the following five-zero
textures because they fit well with experimental quark masses and mixing
parameters Giraldo:2011ya ; Giraldo:2015cpp ; Giraldo:2018mqi :

$\begin{split}M^{U}&=\begin{pmatrix}0&0&C_{u}\\\ 0&A_{u}&B_{u}\\\
C_{u}^{*}&B_{u}^{*}&D_{u}\end{pmatrix},\\\ M^{D}&=\begin{pmatrix}0&C_{d}&0\\\
C_{d}^{*}&0&B_{d}\\\ 0&B_{d}^{*}&A_{d}\end{pmatrix}.\end{split}$ (3)

In addition, the phases in $M^{D}$ can be removed by a weak basis
transformation (WBT) Giraldo:2011ya ; Branco:1988iq ; Branco:1999nb , so that
they are absorbed by the off-diagonal terms in $M^{U}$. In this way, the mass
matrices (3) can be rewritten as:

$\begin{split}M^{U}&=\begin{pmatrix}0&0&|C_{u}|e^{i\phi_{C_{u}}}\\\
0&A_{u}&|B_{u}|e^{i\phi_{B_{u}}}\\\
|C_{u}|e^{-i\phi_{C_{u}}}&|B_{u}|e^{-i\phi_{B_{u}}}&D_{u}\end{pmatrix},\\\
M^{D}&=\begin{pmatrix}0&|C_{d}|&0\\\ |C_{d}|&0&|B_{d}|\\\
0&|B_{d}|&A_{d}\end{pmatrix},\end{split}$ (4)

By applying the trace and the determinant to the mass matrices (4), before and
after the diagonalization process, the free real parameters of $M^{U}$ and
$M^{D}$ can be written in terms of their masses:

$\displaystyle D_{u}$ $\displaystyle=m_{u}-m_{c}+m_{t}-A_{u},$ (5a)
$\displaystyle|B_{u}|$
$\displaystyle=\sqrt{\frac{(A_{u}-m_{u})(A_{u}+m_{c})(m_{t}-A_{u})}{A_{u}}},$
(5b) $\displaystyle|C_{u}|$
$\displaystyle=\sqrt{\frac{m_{u}\,m_{c}\,m_{t}}{A_{u}}},$ (5c) $\displaystyle
A_{d}$ $\displaystyle=m_{d}-m_{s}+m_{b},$ (5d) $\displaystyle|B_{d}|$
$\displaystyle=\sqrt{\frac{(m_{b}-m_{s})(m_{d}+m_{b})(m_{s}-m_{d})}{m_{d}-m_{s}+m_{b}}},$
(5e) $\displaystyle|C_{d}|$
$\displaystyle=\sqrt{\frac{m_{d}\,m_{s}\,m_{b}}{m_{d}-m_{s}+m_{b}}}.$ (5f)

A possibility that works very well is to consider the second generation of
quark masses to be negative, i.e., with eigenvalues $-m_{c}$ and $-m_{s}$. And
$A_{u}$ is a free parameter, whose value, determined by the quark mass
hierarchy, must be in the following range:

$m_{u}\leq A_{u}\leq m_{t}.$ (6)

The exact analytical procedure for diagonalizing the mass matrices (4) is
indicated in Appendix C.

###  II.2 Lepton sector

In this work, we will consider Dirac neutrinos. This is achieved, in part, by
extending the SM with right-handed neutrinos. In this way, we can carry out a
treatment similar to that of the quark sector, that is, the mass matrices of
the lepton sector can be considered Hermitian and the weak basis
transformation (WBT) can be applied Branco:1999nb ; Giraldo:2011ya . In the
literature, work has been done considering various texture-zeros for the Dirac
mass matrices of the lepton sector Hagedorn:2005kz ; Ahuja:2009jj ;
Liu:2012axa ; Verma:2013cza ; Ludl:2014axa ; Verma:2014lpa ; Fakay:2014nea ;
Cebola:2015dwa ; Gautam:2015kya ; Ahuja:2016san ; Singh:2018lao ;
Ahuja:2018fmw ; Barranco:2012ci . In our treatment, we are going to consider
the following five-zero texture model studied in the paper Benavides:2020pjx ,
which can accurately reproduce the Pontecorvo-Maki-Nakagawa-Sakata (PMNS)
mixing matrix $V_{\text{PMNS}}$ (mixing angles and the CP violating phase),
the charged lepton masses, and the squared mass differences in the normal mass
ordering.

$\begin{split}M^{N}&=\begin{pmatrix}0&|C_{\nu}|e^{ic_{\nu}}&0\\\
|C_{\nu}|e^{-ic_{\nu}}&E_{\nu}&|B_{\nu}|e^{ib_{\nu}}\\\
0&|B_{\nu}|e^{-ib_{\nu}}&A_{\nu}\end{pmatrix},\\\
M^{E}&=\begin{pmatrix}0&|C_{\ell}|&0\\\ |C_{\ell}|&0&|B_{\ell}|\\\
0&|B_{\ell}|&A_{\ell}\end{pmatrix}.\end{split}$ (7)

Without loss of generality, by using a WBT, the phases of the charged lepton
mass matrix, $M^{E}$, can be absorbed into the entries $C_{\nu}$ and $B_{\nu}$
of the neutrino mass matrix, $M^{N}$. Similarly, as was done in the case of
the quark sector, the parameters present in the mass matrices of the lepton
sector (7) can be expressed in terms of the masses of the charged leptons
$m_{e},m_{\mu}$ and $m_{\tau}$ and the masses of the neutrinos $m_{1},m_{2}$
and $m_{3}$, in the normal ordering ($m_{1}<m_{2}<m_{3}$):

$\displaystyle A_{\ell}$ $\displaystyle=m_{e}-m_{\mu}+m_{\tau},$ (8a)
$\displaystyle|B_{\ell}|$
$\displaystyle=\sqrt{\frac{(m_{\tau}-m_{\mu})(m_{e}+m_{\tau})(m_{\mu}-m_{e})}{m_{e}-m_{\mu}+m_{\tau}}},$
(8b) $\displaystyle|C_{\ell}|$
$\displaystyle=\sqrt{\frac{m_{e}m_{\mu}m_{\tau}}{m_{e}-m_{\mu}+m_{\tau}}},$
(8c) $\displaystyle E_{\nu}$ $\displaystyle=m_{1}-m_{2}+m_{3}-A_{\nu},$ (8d)
$\displaystyle|B_{\nu}|$
$\displaystyle=\sqrt{\frac{(A_{\nu}-m_{1})(A_{\nu}+m_{2})(m_{3}-A_{\nu})}{A_{\nu}}},$
(8e) $\displaystyle|C_{\nu}|$
$\displaystyle=\sqrt{\frac{m_{1}m_{2}m_{3}}{A_{\nu}}},$ (8f)

where the values of the masses and the parameter $A_{\nu}$ are given in Table
5. Furthermore, for the adjustment of the mass matrices (7) it is very
convenient to assume that the eigenvalues associated with the masses of the
second family, $-m_{2}$ and $-m_{\mu}$, are negative quantities. The exact
diagonalizing matrices of the mass matrices (7) are shown in appendix C,
equations (49) and (48).

##  III PQ symmetry and the minimal particle content

###  III.1 Yukawa Lagrangian and the PQ symmetry

The texture-zeros of the mass matrices defined in the equations (4) and (7)
can be generated by imposing a Peccei-Queen symmetry $U(1)_{PQ}$ on the
Lagrangian model, Eq. (9) Bjorkeroth:2018ipq ; Garnica:2019hvn ;
Ringwald:2015dsf . As will be explained below, the minimal Lagrangian that
allows us to implement this symmetry is given by Giraldo:2020hwl ;
Brivio:2017ije

$\displaystyle\mathcal{L}_{\text{LO}}$
$\displaystyle\supset(D_{\mu}\Phi^{\alpha})^{\dagger}D^{\mu}\Phi^{\alpha}+\sum_{\psi}i\bar{\psi}\gamma^{\mu}D_{\mu}\psi+\sum_{i=1}^{2}(D_{\mu}S_{i})^{\dagger}D^{\mu}S_{i}$
$\displaystyle-\Bigg{(}\bar{q}_{Li}y^{D\alpha}_{ij}\Phi^{\alpha}d_{Rj}+\bar{q}_{Li}y^{U\alpha}_{ij}\tilde{\Phi}^{\alpha}u_{Rj}$
$\displaystyle+\bar{\ell}_{Li}y^{E\alpha}_{ij}\Phi^{\alpha}e_{Rj}+\bar{\ell}_{Li}y^{N\alpha}_{ij}\tilde{\Phi}^{\alpha}\nu_{Rj}+\text{h.c}\Bigg{)}$
$\displaystyle+$
$\displaystyle(\lambda_{Q}\bar{Q}_{R}Q_{L}S_{2}+\text{h.c})-V(\Phi,S_{1},S_{2})\,.$
(9)

As it was shown in reference Giraldo:2020hwl , at least four Higgs doublets
are required to generate the quark mass textures, therefore $\alpha=1,2,3,4$.
In (9) $i,j$ are family indices (there is an implicit sum over repeated
indices). The superscripts $U$, $D$, $E$, $N$ refer to up-type quarks, down-
type quarks, electron-like and neutrino-like fermions, respectively; and
$D_{\mu}=\partial_{\mu}+i\Gamma_{\mu}$ is the covariant derivative in the SM.
The scalar potential $V(\Phi,S_{1},S_{2})$ is shown in appendix A (for further
details, see reference Giraldo:2020hwl ). In Eq. (9) $\psi$ stands for the SM
fermion fields plus the heavy quark $Q$ (see Tables 1 and 2). As it is shown
in Table 2 the PQ charges of the heavy quark can be chosen in such a way that
only the interaction with the scalar singlet $S_{2}$ is allowed. We assign
$Q_{\text{PQ}}$ charges for the left-handed quark doublets ($q_{L}$):
$x_{q_{i}}$, right-handed up-type quark singlets ($u_{R}$): $x_{u_{i}}$,
right-handed down-type quark singlets ($d_{R}$): $x_{d_{i}}$, left-handed
lepton doublets ($\ell_{L}$): $x_{\ell_{i}}$, right-handed charged leptons
($e_{R}$): $x_{e_{i}}$ and right-handed Dirac neutrinos ($\nu_{R}$):
$x_{\nu_{i}}$ for each family ($i=1,2,3$). We follow a similar notation for
the scalar doublets, $x_{\phi_{\alpha}}$ ($\alpha=1,2,3,4$), and the scalar
singlets $x_{{}_{S_{1,2}}}$.

In this work, the PQ charges assigned to the quark sector and the scalar
sector, as well as the VEVs assigned to the scalar doublets, will be the same
as those assigned in Giraldo:2020hwl  (Tables (1) and (2)), and we will adjust
the PQ charges of the lepton sector to reproduce the texture-zeros given in
Eq. (7). To forbid a given entry in the lepton mass matrices, the
corresponding sum of PQ charges must be equal to zero, so that we can obtain
texture-zeros by imposing the following conditions:

$\displaystyle M^{N}=\begin{pmatrix}0&x&0\\\ x&x&x\\\
0&x&x\end{pmatrix}\longrightarrow\begin{pmatrix}S_{11}^{N\alpha}\neq
0&S_{12}^{N\alpha}=0&S_{13}^{N\alpha}\neq 0\\\
S_{21}^{N\alpha}=0&S_{22}^{N\alpha}=0&S_{23}^{N\alpha}=0\\\
S_{31}^{N\alpha}\neq 0&S_{32}^{N\alpha}=0&S_{33}^{N\alpha}=0\end{pmatrix},$
(10) $\displaystyle M^{E}=\begin{pmatrix}0&x&0\\\ x&0&x\\\
0&x&x\end{pmatrix}\longrightarrow\begin{pmatrix}S_{11}^{E\alpha}\neq
0&S_{12}^{E\alpha}=0&S_{13}^{E\alpha}\neq 0\\\
S_{21}^{E\alpha}=0&S_{22}^{E\alpha}\neq 0&S_{23}^{E\alpha}=0\\\
S_{31}^{E\alpha}\neq 0&S_{32}^{E\alpha}=0&S_{33}^{E\alpha}=0\end{pmatrix},$
(11)

where $S_{ij}^{N\alpha}=(-x_{\ell_{i}}+x_{\nu_{j}}-x_{\phi_{\alpha}})$ and
$S_{ij}^{E\alpha}=(-x_{\ell_{i}}+x_{e_{j}}+x_{\phi_{\alpha}})$.

Since the PQ charges of the Higgs doublets ($\alpha=1,2,3,4$) are already
given, the possible solutions of (10) and (11) are strongly constrained. Table
1 provides a solution for the PQ charges of the lepton sector.

Particles | Spin | $SU(3)_{C}$ | $SU(2)_{L}$ | $U(1)_{Y}$ | $U_{\text{PQ}}(i=1)$ | $U_{\text{PQ}}(i=2)$ | $U_{\text{PQ}}(i=3)$ | $Q_{\text{PQ}}$
---|---|---|---|---|---|---|---|---
$q_{Li}$ | 1/2 | 3 | 2 | 1/6 | $-2s_{1}+2s_{2}+\alpha_{q}$ | $-s_{1}+s_{2}+\alpha_{q}$ | $\alpha_{q}$ | $x_{q_{i}}$
$u_{Ri}$ | 1/2 | 3 | 1 | 2/3 | $s_{1}+\alpha_{q}$ | $s_{2}+\alpha_{q}$ | $-s_{1}+2s_{2}+\alpha_{q}$ | $x_{u_{i}}$
$d_{Ri}$ | 1/2 | 3 | 1 | $-1/3$ | $2s_{1}-3s_{2}+\alpha_{q}$ | $s_{1}-2s_{2}+\alpha_{q}$ | $-s_{2}+\alpha_{q}$ | $x_{d_{i}}$
$\ell_{Li}$ | 1/2 | 1 | 2 | $-1/2$ | $-2s_{1}+2s_{2}+\alpha_{\ell}$ | $-s_{1}+s_{2}+\alpha_{\ell}$ | $\alpha_{\ell}$ | $x_{\ell_{i}}$
$e_{Ri}$ | 1/2 | 1 | 1 | $-1$ | $2s_{1}-3s_{2}+\alpha_{\ell}$ | $s_{1}-2s_{2}+\alpha_{\ell}$ | $-s_{2}+\alpha_{\ell}$ | $x_{e_{i}}$
$\nu_{Ri}$ | 1/2 | 1 | 1 | 0 | $-4s_{1}+5s_{2}+\alpha_{\ell}$ | $-s_{1}+2s_{2}+\alpha_{\ell}$ | $s_{2}+\alpha_{\ell}$ | $x_{\nu_{i}}$
Table 1: Particle content. The subindex $i=1,2,3$ stand for the family number in the interaction basis. Columns 6-8 are the Peccei-Quinn charges, $Q_{PQ}$, for each family of quarks and leptons in the SM. $s_{1},s_{2}$ and $\alpha$ are real parameters, with $s_{1}\neq s_{2}$. Particles | Spin | $SU(3)_{C}$ | $SU(2)_{L}$ | $U(1)_{Y}$ | $U_{\text{PQ}}$ | $Q_{\text{PQ}}$
---|---|---|---|---|---|---
$\Phi_{1}$ | 0 | 1 | 2 | 1/2 | $s_{1}$ | $x_{\phi_{1}}$
$\Phi_{2}$ | 0 | 1 | 2 | 1/2 | $s_{2}$ | $x_{\phi_{2}}$
$\Phi_{3}$ | 0 | 1 | 2 | 1/2 | $-s_{1}+2s_{2}$ | $x_{\phi_{3}}$
$\Phi_{4}$ | 0 | 1 | 2 | 1/2 | $-3s_{1}+4s_{2}$ | $x_{\phi_{4}}$
$Q_{L}$ | 1/2 | 3 | 0 | 0 | $x_{Q_{L}}$ | $x_{Q_{L}}$
$Q_{R}$ | 1/2 | 3 | 0 | 0 | $x_{Q_{R}}$ | $x_{Q_{R}}$
$S_{1}$ | 0 | 1 | 1 | 0 | $s_{1}-s_{2}$ | $x_{{}_{S_{1}}}$
$S_{2}$ | 0 | 1 | 1 | 0 | $x_{Q_{R}}-x_{Q_{L}}$ | $x_{{}_{S_{2}}}$
Table 2: Beyond the SM fields and their respective PQ charges. The parameters
$s_{1},s_{2}$ are reals, with $s_{1}\neq s_{2}$ and $x_{Q_{R}}\neq x_{Q_{L}}$.

In our model we include two scalar singlets $S_{1}$ and $S_{2}$ that break the
global symmetry $U(1)_{PQ}$. The QCD anomaly of the PQ charges is

$\displaystyle
N=2\sum_{i}^{3}x_{q_{i}}-\sum_{i}^{3}x_{u_{i}}-\sum_{i}^{3}x_{d_{i}}+A_{Q}\,,$
(12)

where $A_{Q}=x_{QL}-x_{QR}$ is the contribution to the anomaly of the heavy
quark $Q$, which is a singlet under the electroweak gauge group, with left
(right) Peccei-Quinn charges $x_{QL,R}$, respectively. We can write the
charges as a function of $N$ (since $N$ must be different from zero), such
that

$\displaystyle s_{1}=\frac{N}{9}\hat{s}_{1},\hskip
28.45274pts_{2}=\frac{N}{9}\left(\epsilon+\hat{s}_{1}\right),\ \
\text{with}\hskip 14.22636pt\epsilon=1-\frac{A_{Q}}{N}\,,$ (13)

where $\hat{s}_{1}$ and $\epsilon$ are arbitrary real numbers. To solve the
strong CP problem with $N\neq 0$ and simultaneously generate the texture-zeros
in the mass matrices, it is necessary to maintain
$\epsilon=\frac{9(s_{1}-s_{2})}{N}\neq 0$. With these definitions for Flavor-
Changing Neutral Currents (FCNC) observables, the relevant parameters are
$\hat{s}_{1}$ and $\epsilon$. This parameterization is quite convenient (for
those cases where the parameters $\alpha_{q}$ and $\alpha_{\ell}$ are not
relevant) because by fixing $N$ and $f_{a}$, we can vary $\hat{s}_{1}$ and
$\epsilon$ for a fixed $\Lambda_{\text{PQ}}=f_{a}N$ in such a way that the
parameter space naturally reduces to two dimensions.

##  IV Naturalness of Yukawa couplings

###  IV.1 The mass matrices in the quark sector

In reference Giraldo:2020hwl , it was shown that to generate five texture
zeros in the quark mass matrices (3), as a consequence of a PQ symmetry, it is
necessary to include at least four scalar doublets in the model. After
spontaneous symmetry breaking, the quark sector mass matrices take on the
following form:

$\displaystyle
M^{U}=\hat{v}_{\alpha}y^{U\alpha}_{ij}=\begin{pmatrix}0&0&y^{U1}_{13}\hat{v}_{1}\\\
\vspace{0.2cm}0&y^{U1}_{22}\hat{v}_{1}&y^{U2}_{23}\hat{v}_{2}\\\
y^{U1^{*}}_{13}\hat{v}_{1}&y^{U2^{*}}_{23}\hat{v}_{2}&y^{U3}_{33}\hat{v}_{3}\end{pmatrix},\hskip
14.22636ptM^{D}=\hat{v}_{\alpha}y^{D\alpha}_{ij}=\begin{pmatrix}0&|y^{D4}_{12}|\hat{v}_{4}&0\\\
|y^{D4}_{12}|\hat{v}_{4}&0&|y^{D3}_{23}|\hat{v}_{3}\\\
0&|y^{D3}_{23}|\hat{v}_{3}&y^{D2}_{33}\hat{v}_{2}\end{pmatrix},$ (14)

where the $\hat{v}_{i}$ are defined in terms of the vacuum expectation values,
$\hat{v}_{i}=v_{i}/\sqrt{2}$. In Giraldo:2020hwl  it was shown that the five-
texture zeros (4) are flexible enough to set the quark Yukawa couplings close
to 1 for most of them (except for $y^{U2}_{23}$, $y^{D3}_{23}$ and
$y_{13}^{U1}$), in this way we obtain:

$\displaystyle\hat{v}_{1}=1.71\,\text{GeV},\hskip
14.22636pt\hat{v}_{2}=2.91\,\text{GeV},\hskip
14.22636pt\hat{v}_{3}=174.085\,\text{GeV},\hskip
14.22636pt\hat{v}_{4}=13.3\,\text{MeV}.$ (15)

As we can see, the hermiticity of the mass matrices is not fully achieved, but
it is good to impose it for several reasons: (i) In the SM and its extensions,
in which the right chirality fields are singlets under $SU(2)$, the mass
matrices can be assumed Hermitian without losing generality, (ii) the previous
fact allows us to consider Hermitian mass matrices, even after imposing an
additional $PQ$ symmetry in the model, (iii) we can implement the WBT method
Giraldo:2011ya , and (iv) there is an extensive literature on physically
viable Hermitian mass matrices. It is important to noticing that the mass
matrices in Eq. (14) are Hermitian.

###  IV.2 The mass matrices in the lepton sector

We can obtain the lepton mass matrices by starting from the Yukawa Lagrangian
(9), which is invariant under the Peccei-Quinn $U(1)_{PQ}$ symmetry, and
taking into account the Yukawa parameters and expectation values (15). After
the spontaneous symmetry breaking, the mass matrices for neutral and charged
leptons are given respectively by Benavides:2020pjx :

$M^{N}=\hat{v}_{\alpha}y^{N\alpha}_{ij}=\begin{pmatrix}0&y^{N1}_{12}\hat{v}_{1}&0\\\
\vspace{0.2cm}y^{N4}_{21}\hat{v}_{4}&y^{N2}_{22}\hat{v}_{2}&y^{N1}_{23}\hat{v}_{1}\\\
0&y^{N3}_{32}\hat{v}_{3}&y^{N2}_{33}\hat{v}_{2}\end{pmatrix},\hskip
14.22636ptM^{E}=\hat{v}_{\alpha}y^{E\alpha}_{ij}=\begin{pmatrix}0&|y^{E4}_{12}|\hat{v}_{4}&0\\\
|y^{E4}_{12}|\hat{v}_{4}&0&|y^{E3}_{23}|\hat{v}_{3}\\\
0&|y^{E3}_{23}|\hat{v}_{3}&y^{E2}_{33}\hat{v}_{2}\end{pmatrix}.$ (16)

As we previously mentioned, at least four Higgs doublets are needed to obtain
the five texture-zeros for the chosen quark mass matrices. Our goal in this
work is to keep the same number of Higgs doublets and their respective $PQ$
charges to generate the mass matrices and texture zeros for the lepton sector,
Eq. (7). To get an Hermitian mass matrix $M^{N}$, it is necessary to impose
$y_{21}^{N4}/y_{12}^{N1*}=\hat{v}_{1}/\hat{v}_{4}$ and
$y_{32}^{N3}/y_{23}^{N1*}=\hat{v}_{1}/\hat{v}_{3}$, requiring that the
diagonal elements be real, i.e., $y^{N2}_{22}=y^{N2*}_{22}$ and
$y^{N2}_{33}=y^{N2*}_{33}$. On the other hand, to obtain a symmetric mass
matrix, $M^{E}$, for the charged leptons, it is sufficient to assume that the
Yukawa couplings are Hermitian. Through these choices it is possible to avoid
additional Higgs doublets.

Based on the results of Table 5, Appendix C, and the relationships established
in (8), we find the following values for the Yukawa couplings of the lepton
sector:

$\begin{split}|y^{E4}_{12}|&=0.569582,\hskip
42.67912pt|y^{E3}_{23}|=0.00248291,\hskip 31.2982pty^{E2}_{33}=0.574472,\\\
|y^{N1}_{12}|&=4.74362\times 10^{-6},\hskip
14.22636pt|y^{N4}_{21}|=0.000609894,\hskip 25.6073pty^{N2}_{22}=6.68808\times
10^{-6},\\\ |y^{N1}_{23}|&=0.0000159881,\hskip
22.76228pt|y^{N3}_{32}|=1.57047\times 10^{-7},\hskip
14.22636pty^{N2}_{33}=8.65364\times 10^{-6}.\end{split}$ (17)

To reproduce the neutrino masses quoted in Benavides:2020pjx , in the SM is
required a Yukawa coupling around $10^{-14}$. In our case, the smallest Yukawa
coupling is $10^{-7}$, which significantly reduces the fine-tuning in
comparison to that given by the SM.

##  V The Effective Lagrangian

The strongest constraints on non-universal PQ charges come from the FCNC. To
determine these constraints, we start by writing the most general next-to-
leading order (NLO) effective Lagrangian as Georgi:1986df ; Gavela:2019wzg :

$\displaystyle\mathcal{L}_{\text{NLO}}=c_{a\Phi^{\alpha}}O_{a\Phi^{\alpha}}+c_{1}\frac{\alpha_{1}}{8\pi}O_{B}+c_{2}\frac{\alpha_{2}}{8\pi}O_{W}+c_{3}\frac{\alpha_{3}}{8\pi}O_{G},$
(18)

$c_{a\Phi^{\alpha}}$ and $c_{1,2,3}$ are Wilson coefficients;
$\alpha_{1,2,3}=\frac{g_{1,2,3}^{2}}{4\pi}$, where $g_{1,2,3}$ are the
coupling strengths of the electroweak interaction in the interaction basis;
and the Wilson operators are:

$\displaystyle O_{a\Phi}=$ $\displaystyle
i\frac{\partial^{\mu}a}{\Lambda_{\text{PQ}}}\left((D_{\mu}\Phi^{\alpha})^{\dagger}\Phi^{\alpha}-\Phi^{\alpha\dagger}(D_{\mu}\Phi^{\alpha})\right),\hskip
14.22636ptO_{B}=-\frac{a}{\Lambda_{\text{PQ}}}B_{\mu\nu}\tilde{B}^{\mu\nu},$
$\displaystyle O_{W}=$
$\displaystyle-\frac{a}{\Lambda_{\text{PQ}}}W_{\mu\nu}^{a}\tilde{W}^{a\mu\nu},\hskip
99.58464ptO_{G}=-\frac{a}{\Lambda_{\text{PQ}}}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}\,,$
(19)

where $B$, $W^{a}$ and $G^{a}$ correspond to the gauge fields associated with
the SM gauge groups $U(1)_{Y}$, $SU(2)_{L}$ and $SU(3)_{C}$, respectively. $a$
is the axion field which corresponds to the CP odd component of $S_{1}$. It is
possible to redefine the fields by multiplying by a phase

$\displaystyle\Phi^{\alpha}$ $\displaystyle\longrightarrow
e^{i\frac{x_{{}_{\Phi^{\alpha}}}}{\Lambda_{\text{PQ}}}a}\Phi^{\alpha},$
$\displaystyle\psi_{L}$ $\displaystyle\longrightarrow
e^{i\frac{x_{\psi_{L}}}{\Lambda_{\text{PQ}}}a}\psi_{L},$
$\displaystyle\psi_{R}$ $\displaystyle\longrightarrow
e^{i\frac{x_{\psi_{R}}}{\Lambda_{\text{PQ}}}a}\psi_{R},$ $\displaystyle S_{i}$
$\displaystyle\longrightarrow
e^{i\frac{x_{{}_{S_{i}}}}{\Lambda_{\text{PQ}}}a}S_{i}.$ (20)

In this expression, $x_{\psi}$ corresponds to the PQ charges of the SM
fermions, i.e., $\\{x_{\psi_{L,R}}\\}=\\{x_{q_{i}}$, $x_{u_{i}}$,
$x_{d_{i}},x_{l_{i}}$, $x_{e_{i}}$, $x_{\nu_{i}}\\}$ and
$\\{x_{\Phi^{\alpha}}\\}$ are the PQ charges of the Higgs doublets
$\\{\Phi^{\alpha}\\}$. Replacing these definitions in the kinetic terms of Eq.
(9), we obtain new contributions to the effective Lagrangian Eq. (18) (the NLO
contributions in the non-derivative terms cancel out). The leading order (LO)
terms in $\Lambda^{-1}_{\text{PQ}}$ can be written as Georgi:1986df ;
Brivio:2017ije :

$\displaystyle\mathcal{L}_{\text{NLO}}\longrightarrow\mathcal{L}_{\text{NLO}}+\Delta\mathcal{L}_{\text{NLO}},$
(21)

where

$\displaystyle\Delta\mathcal{L}_{\text{NLO}}=\Delta\mathcal{L}_{K^{\Phi}}+\Delta\mathcal{L}_{K^{\psi}}+\Delta\mathcal{L}_{K^{S}}+\Delta\mathcal{L}(F_{\mu\nu}),$
(22)

with

$\displaystyle\Delta\mathcal{L}_{K^{\Phi}}=$ $\displaystyle
ix_{\Phi^{\alpha}}\frac{\partial^{\mu}a}{\Lambda_{\text{PQ}}}\left[(D_{\mu}\Phi^{\alpha})^{\dagger}\Phi^{\alpha}-\Phi^{\alpha\dagger}(D_{\mu}\Phi^{\alpha})\right],$
$\displaystyle\Delta\mathcal{L}_{K^{\psi}}=$
$\displaystyle\frac{\partial_{\mu}a}{2\Lambda_{\text{PQ}}}\sum_{\psi}(x_{\psi_{L}}-x_{\psi_{R}})\bar{\psi}\gamma^{\mu}\gamma^{5}\psi-(x_{\psi_{L}}+x_{\psi_{R}})\bar{\psi}\gamma^{\mu}\psi,$
$\displaystyle\Delta\mathcal{L}_{K^{S}}=$ $\displaystyle
ix_{{}_{S_{i}}}\frac{\partial^{\mu}a}{\Lambda_{\text{PQ}}}\left[(D_{\mu}S_{i})^{\dagger}S_{i}-S^{\dagger}_{i}(D_{\mu}S_{i})\right]+\text{h.c}\,.$
(23)

The field redefinitions (V) induce a modification in the measure of the
functional path integral whose effects can be obtained from the divergence of
the axial-vector current:
$J^{PQ5}_{\mu}=\sum_{\psi}(x_{\psi_{L}}-x_{\psi_{R}})\bar{\psi}\gamma_{\mu}\gamma^{5}\psi$
Bauer:2017ris ,

$\displaystyle\partial^{\mu}J^{PQ5}_{\mu}=$
$\displaystyle\sum_{\psi}2im_{\psi}(x_{\psi_{L}}-x_{\psi_{R}})\bar{\psi}\gamma^{5}\psi-\sum_{\psi}(x_{\psi_{L}}-x_{\psi_{R}})\frac{\alpha_{1}Y^{2}(\psi)}{2\pi}B_{\mu\nu}\tilde{B}^{\mu\nu}$
$\displaystyle-$ $\displaystyle\sum_{\text{ $SU(2)_{L}$
doublets}}x_{\psi_{L}}\frac{\alpha_{2}}{4\pi}W_{\mu\nu}^{a}\tilde{W}^{a\mu\nu}-\sum_{\text{$SU(3)$
triplets}}(x_{\psi_{L}}-x_{\psi_{R}})\frac{\alpha_{3}}{4\pi}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu},$
(24)

where the hypercharge is normalized by $Q=T_{3L}+Y$. The relation (V) is an
on-shell relation, which is consistent with the momentum of an on-shell axion.
Substituting this result into
$\mathcal{L}_{K^{\psi}}=\frac{\partial^{\mu}a}{2\Lambda_{\text{PQ}}}J^{PQ5}_{\mu}=-\frac{a}{2\Lambda_{\text{PQ}}}\partial^{\mu}J^{PQ5}_{\mu}$
we obtain new contributions to the leading-order Wilson coefficients
Salvio:2013iaa

$\displaystyle c_{1}$ $\displaystyle\longrightarrow c_{1}-\frac{1}{3}\Sigma
q+\frac{8}{3}\Sigma u+\frac{2}{3}\Sigma d-\Sigma\ell+2\Sigma e,$
$\displaystyle c_{2}$ $\displaystyle\longrightarrow c_{2}-3\Sigma
q-\Sigma\ell,$ $\displaystyle c_{3}$ $\displaystyle\longrightarrow
c_{3}-2\Sigma q+\Sigma u+\Sigma d-A_{Q},$ (25)

where $\Sigma q\equiv x_{q_{1}}+x_{q_{2}}+x_{q_{3}}$. The corresponding NLO
Lagrangian is

$\displaystyle\Delta\mathcal{L}(F_{\mu\nu})=$
$\displaystyle\frac{a}{\Lambda_{\text{PQ}}}\frac{\alpha_{1}}{8\pi}B_{\mu\nu}\tilde{B}^{\mu\nu}\left(\frac{1}{3}\Sigma
q-\frac{8}{3}\Sigma u-\frac{2}{3}\Sigma d+\Sigma\ell-2\Sigma e\right)$
$\displaystyle+$
$\displaystyle\frac{a}{\Lambda_{\text{PQ}}}\frac{\alpha_{2}}{8\pi}W_{\mu\nu}^{a}\tilde{W}^{a\mu\nu}\left(3\Sigma
q+\Sigma\ell\right)$ $\displaystyle+$
$\displaystyle\frac{a}{\Lambda_{\text{PQ}}}\frac{\alpha_{3}}{8\pi}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}\left(2\Sigma
q-\Sigma u-\Sigma d+A_{Q}\right).$ (26)

It is convenient to define $c_{3}^{\text{eff}}=c_{3}-2\Sigma q+\Sigma u+\Sigma
d-A_{Q}=-N$. In our case, $c_{i}=0$ and the only contributions to
$c_{i}^{\text{eff}}$ come from the anomaly. It is customary to define
$\Lambda_{\text{PQ}}=f_{a}\lvert c_{3}^{\text{eff}}\rvert$ to include the
factor $c_{3}^{\text{eff}}$ in the normalization of the PQ charges. From now
on, we will assume that all the PQ charges are normalized in this way, so that
$x_{\psi}$ corresponds to $x_{\psi}/c_{3}^{\text{eff}}$. For normalized
charges, $c_{3}^{\text{eff}}=1$, therefore, we still maintain the general form
despite writing all the expressions in terms of the effective scale $f_{a}$.

The scalar fields and their PQ charges are the same as in the reference
Giraldo:2020hwl , so the scalar potential $V(\Phi,S)$ is identical to that of
the mentioned reference. With the VEVs and couplings given in Giraldo:2020hwl
, the model reproduces the mass of the SM Higgs, while the masses of the
exotic scalars are above the TeV scale. This potential has the appropriate
number of Goldstone bosons to give masses to the SM gauge bosons
$Z^{0},W^{\pm}$ and has an extra field that can be identified with the axion
$a$.

##  VI Low energy constraints

###  VI.1 Flavor changing neutral currents

Due to the non-universal PQ charges in our model, a tree-level analysis of
flavor-changing neutral currents is necessary. As mentioned in reference
DiLuzio:2020wdo , the strongest limits on the axion-quark FCNC couplings come
from meson decays in light mesons and missing energy. The decays
$K^{\pm}\rightarrow\pi^{\pm}a$ provide the tightest limits (NA62 Collaboration
NA62:2021zjw ) for the axion mass DiLuzio:2020wdo . Currently the most
restrictive limits come from the semileptonic decays of kaons
$K^{\pm}\rightarrow\pi^{\pm}\bar{\nu}\nu$ and leptons
$\ell_{1}\rightarrow\ell_{2}$+missing energy. From the term
$\Delta\mathcal{L}_{K^{\psi}}$, we obtain the vector and axial couplings for a
multi-Higgs sector model, as shown in references DiLuzio:2020wdo ;
Giraldo:2020hwl

$\displaystyle\Delta\mathcal{L}_{K^{\psi}}=-\partial_{\mu}a\bar{f}_{i}\gamma^{\mu}\left(g_{af_{i}f_{j}}^{V}+\gamma^{5}g_{af_{i}f_{j}}^{A}\right)f_{j},$
(27)

where

$\displaystyle
g_{af_{i}f_{j}}^{V,A}=\frac{1}{2f_{a}c^{\text{eff}}_{3}}\Delta^{Fij}_{V,A},$
(28)

where $\Delta^{Fij}_{V,A}=\Delta^{Fij}_{RR}(d)\pm\Delta^{Fij}_{LL}(q)$ with
$\Delta^{Fij}_{LL}(q)=\left(U^{F}_{L}x_{q}~{}U_{L}^{F\dagger}\right)^{ij}$ and
$\Delta^{Fij}_{RR}(d)=\left(U^{F}_{R}x_{d}~{}U_{R}^{F\dagger}\right)^{ij}$. In
these expressions, $F$ stands for $U,D,N$ or $E$ and the $U_{L,R}^{F}$ are de
diagonalizing matrices (see Appendix C). In Eq. (28), we normalize the charges
using $c^{\text{eff}}_{3}$, as explained in the last paragraph of section V
(in other references $|c^{\text{eff}}_{3}|=|N|$ is considered, corresponding
to the $SU(3)\times U(1)_{PQ}$ anomaly). The branching ratio for lepton decays
$\ell_{i}\rightarrow\ell_{j}\,a$ is given by Bjorkeroth:2018dzu

$\text{Br}(\ell_{1}\rightarrow\ell_{2}\,a)=\frac{m_{\ell_{1}}^{3}}{16\pi\Gamma(\ell_{1})}\left(1-\frac{m_{\ell_{2}}^{2}}{m_{\ell_{1}}^{2}}\right)^{3}|g_{a\ell_{1}\ell_{2}}|^{2},$

in this expression, the vector and axial couplings contribute in the same way

$|g_{a\ell_{1}\ell_{2}}|^{2}=|g_{a\ell_{1}\ell_{2}}^{V}|^{2}+|g_{a\ell_{1}\ell_{2}}^{A}|^{2}.$

where $m$ is the mass of the leptons and $\Gamma(\ell_{i})$ is the total decay
width of the particle $\ell_{j}$.

For the lepton decay $\ell_{i}\rightarrow\ell_{j}\,a\,\gamma$, we can relate
this branching ratio to the branching ratio of the process without the photon
in the final state, according to the expression:

$\text{Br}(\ell_{1}\rightarrow\ell_{2}\,a\,\gamma)=\left(\frac{\alpha}{2\pi}\int
dx\,dyf(x,y)\right)\text{Br}(\ell_{1}\rightarrow\ell_{2}\,a),$

where $\alpha$ is the fine structure constant, and the function,
$f(x,y)=\frac{(1-x)(2-y-xy)}{y^{2}(x+y-1)}$, depends on the mass and the
energies $x=2E_{\ell_{2}}/m_{\ell_{1}}$ and $y=2E_{\gamma}/m_{\ell_{1}}$. For
the lepton decay $\mu\rightarrow e\,a\,\gamma$, the constraints come from the
Crystal Box experiment Bolton:1988af , with cut energies
$E_{\gamma},E_{e}>30\,\text{MeV}$, $\theta_{e\gamma}>140^{\circ}$, where
$\text{cos}\,\theta_{e\gamma}=1+\frac{2(1-x-y)}{xy}$, so that $\int
dx\,dyf(x,y)\approx 0.011$.

$\begin{array}[]{|l|l|}\hline\cr\hskip 14.22636pt\text{Collaboration}&\hskip
28.45274pt\text{Upper bound}\\\ \hline\cr\text{N62
Collaboration\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{NA62:2021zjw}{\@@citephrase{(}}{\@@citephrase{)}}}}&\mathcal{B}\left(K^{+}\rightarrow\pi^{+}a\right)<(10.6^{+4.0}_{-3.4}\big{|}_{stat}\pm
0.9_{syst})\times 10^{-11}\\\ \text{TRIUMF~{}\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Jodidio:1986mz}{\@@citephrase{(}}{\@@citephrase{)}}}
}&\mathcal{B}\left(\mu^{+}\rightarrow e^{+}a\right)<2.6\times 10^{-6}\\\
\text{Crystal Box~{}\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{TWIST:2014ymv}{\@@citephrase{(}}{\@@citephrase{)}}}
}&\mathcal{B}\left(\mu^{+}\rightarrow e^{+}\gamma a\right)<1.1\times
10^{-9}\\\ \text{ARGUS~{}\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Albrecht:1995ht}{\@@citephrase{(}}{\@@citephrase{)}}}
}&\mathcal{B}\left(\tau^{+}\rightarrow e^{+}a\right)<1.5\times 10^{-2}\\\
\text{ARGUS~{}\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Albrecht:1995ht}{\@@citephrase{(}}{\@@citephrase{)}}}
}&\mathcal{B}\left(\tau^{+}\rightarrow\mu^{+}a\right)<2.6\times 10^{-2}\\\
\hline\cr\end{array}$ Table 3: These inequalities come from the window for new
physics in the branching ratio uncertainty of the meson decay in a pair
$\bar{\nu}\nu$.
---
Figure 1: Allowed regions by lepton decays.

In our model, there is a natural alignment between the $\Phi_{3}$ (which is
quite similar to $H_{1}$ in the Georgi basis Georgi:1978ri ) and the standard
model Higgs boson as a consequence of the large suppression of the VEVs of the
scalar doublets $v_{i}$, with $i=1,2,4$, respect to $v_{3}$, the VEV of
$\Phi_{3}$. To some extent, this alignment avoids FCNC involving the SM Higgs
boson Georgi:1978ri ; however, after alignment, there are other sources of
FCNC associated with the additional scalar doublets, which cannot be avoided
by any means; however, as argued in Ref. Giraldo:2020hwl  they are suppressed
by a factor $1/M^{4}$ (where $M>1$TeV is the mass of the exotic scalar
doublets), and therefore, our model avoids these potential sources of FCNC in
agreement with the general argument presented in Georgi:1978ri .

From astrophysical considerations we have: bounds from black holes
superradiance and the SN 1987A upper limit on the neutron electric dipole
moment, which, when combined, impose a constraint on the axion decay constant
in the range DiLuzio:2020wdo  (see Figure 1) : $0.8\times 10^{6}\text{GeV}\leq
f_{a}\leq 2.8\times 10^{17}\text{GeV}$.

###  VI.2 Constraints on the axion-photon coupling

Figure 2:  The excluded parameter space by various experiments corresponds to
the colored regions, the dashed-lines correspond to the projected bounds of
coming experiments looking for axion signals. The gray region corresponds to
the parameter space scanned by our model.

There are several experiments designed to look for exotic particles. The
sources studied in the search for axions are: the solar axion flux
(helioscopes experiments), dark matter halo (haloscopes experiments), and
axions produced in the laboratory.

Among the experiments with the potential to search for evidence of axions in
regions that cover areas within the limits established by the parameters of
our model are: DM-Radio Chaudhuri:2014dla , KLASH Alesini:2017ifp ;
Gatti:2018ojx , ADMX Stern:2016bbw , ALPHA Lawson:2019brd , MADMAX
Beurthey:2020yuq , IAXO Armengaud:2014gea ; Shilon:2013bla  and ABRACADABRA
Ouellet:2018beu . Similarly, some experiments have already ruled out regions
established by the parameters of our model, among which are: ADMX ADMX:2021nhd
; ADMX:2019uok ; ADMX:2020hay , CAST CAST:2017uph ; CAST:2007jps ,
CAPPLee:2020cfj ; Jeong:2020cwz ; CAPP:2020utb , HAYSTACK HAYSTAC:2020kwv ;
HAYSTAC:2018rwy , Solar $\nu$ Gondolo:2008dd , Horizontal Branch Ayala:2014pea
, MUSE Regis:2020fhw  and VIMOS Grin:2006aw

##  VII Summary and conclusions

We have presented a model in which the fermion and scalar fields are charged
under a $U(1)_{PQ}$ Peccei-Quinn symmetry. A recent work Giraldo:2020hwl
showed that at least four Higgs doublets are required to generate Hermitian
mass matrices in the quark sector with five texture-zeros, reproducing the
quark masses, the mixing angles, and the CP-violating phase of the CKM mixing
matrix. In this work, we show that using the same number of Higgs doublets,
without changing the PQ charges in the quark and Higgs sectors, it is possible
to generate Hermitian mass matrices in the lepton sector that reproduce the
neutrino mass-squared differences in the normal mass ordering, the mixing
angles, and the CP-violating phase of the PMNS mixing matrix. This result is
quite non-trivial as we maintain the same four Higgs doublets required in the
quark sector to generate a different texture pattern in the lepton sector.
When compared to the SM, our model has almost all Yukawa couplings close to 1
in the quark sector. In the neutrino sector, the smallest Yukawa coupling is
of the order of $1.6\times 10^{-7}$, which is seven orders of magnitude larger
than the corresponding Yukawa coupling in the SM, so it requires less fine-
tuning than the SM.

The polar decomposition theorem VPrasolov ; Hassani  allows any matrix to be
written as the product of a Hermitian matrix and a unitary matrix. In the SM
and in theories where the right-handed fermion fields are singlets under the
gauge group, it is possible to absorb the unitary matrix into the right-handed
fields by redefining them; from this procedure, we can write any mass matrix
as a Hermitian matrix. In our work, we assume that the mass matrices are
Hermitian in the interaction space, this hypothesis has been used in previous
studies on textures Fritzsch:1999ee ; Xing:2015sva ; Gupta:2013yha ;
Ludl:2015lta ; Giraldo:2011ya ; Fusaoka:1998vc , and it is quite useful for
studying the flavor problem. In our work, we have normalized the PQ charges
with the QCD anomaly $-N$ in such a way that by keeping the parameter
$\epsilon\neq 0$, we obtain the textures of the mass matrices, addressing the
flavor and strong CP problems simultaneously.

If nature is not fine-tuned in a more fundamental high-energy theory, we
expect that, eventually, it will be possible to find a texture that allows us
to obtain all the scales of the SM from the VEVs of a Higgs sector with a
minimal scalar content without the need to adjust the Yukawa couplings.

In our analysis, we report the constraints from lepton decays and compare them
with the constraints from the search for neutrino pairs in charged Kaon decays
$K^{\pm}\longrightarrow\pi^{\pm}\bar{\nu}\nu$. The results are shown in figure
1, where the allowed region in the parameter space generated by $\epsilon$ and
the axion decay constant $f_{a}$ is displayed. This figure shows that the
strongest constraints come from the semileptonic meson decay
$K^{\pm}\rightarrow\pi\nu\bar{\nu}$. It is important to note that the lepton
decays do not further constrain the parameter space of our model (compared to
the region excluded by the meson decay). We also show the excluded regions for
the axion-photon coupling as a function of the axion mass; these results are
summarized in Figure 2; the gray region corresponds to the parameter space of
our model in the interval $-1<\epsilon<1$.

In this article, we have demonstrated that with four Higgs doublets, it is
possible to fit the textures of the mass matrices, both in the lepton and
quark sectors. These matrices generate the masses and the mixing matrices for
quarks and leptons within the experimentally reported values in the
literature. The introduction these doublets improves the fine-tuning problem
of the Yukawa couplings and shows that this approach is a viable way to tackle
the flavor problem. We hope to improve our results in future work by using the
Seesaw mechanism in the lepton sector.

## Acknowledgments

We thank Financial support from “Patrimonio Autónomo Fondo Nacional de
Financiamiento para la Ciencia, la Tecnología y la Innovación, Francisco José
de Caldas”. This research was partly supported by the “Vicerrectoría de
Investigaciones e Interacción Social VIIS de la Universidad de Nariño”,
project numbers 1928, 2172, 2686, 2693 and 2679. E.R and Y.G. acknowledge
additional financial support from Minciencias CD 82315 CT ICETEX 2021-1080.

##  Appendix A The mass operator matrices

The most general Yukawa Lagrangian for the interaction of four Higgs doublets
$\Phi_{\alpha}$ with the SM fermions is given by

$\displaystyle\mathcal{L}=-\bar{q}_{L}^{\prime
i}\Phi_{\alpha}y_{ij}^{D\alpha}d^{\prime j}_{R}-\bar{q}_{L}^{\prime
i}\tilde{\Phi}_{\alpha}y_{ij}^{U\alpha}u^{\prime j}_{R}-\bar{\ell}_{L}^{\prime
i}\Phi_{\alpha}y_{ij}^{E\alpha}e^{\prime j}_{R}-\bar{\ell}_{L}^{\prime
i}\tilde{\Phi}_{\alpha}y_{ij}^{N\alpha}\nu^{\prime j}_{R}+\text{h.c},$ (29)

where a sum is assumed on repeated indices. Here $i,j$ run over $1,2,3$ and
$\alpha$ over $1,2,3,4$. The Higgs boson doublet fields are parameterized as
follows:

$\displaystyle\Phi_{\alpha}=\begin{pmatrix}\phi_{\alpha}^{+}\\\
\frac{v_{\alpha}+h_{\alpha}+i\eta_{\alpha}}{\sqrt{2}}\end{pmatrix},\hskip
28.45274pt\tilde{\Phi}_{\alpha}=i\sigma_{2}\Phi_{\alpha}^{*}.$ (30)

Similar to the two Higgs doublet model Cardozo:2020uol  we rotate the Higgs
fields to the (generalized) Georgi basis, that is,

$\displaystyle\begin{pmatrix}H_{1}\\\ H_{2}\\\ H_{3}\\\
H_{4}\end{pmatrix}=R_{1}(\beta_{1})R_{2}(\beta_{2})R_{3}(\beta_{3})\begin{pmatrix}\Phi_{1}\\\
\Phi_{2}\\\ \Phi_{3}\\\ \Phi_{4}\end{pmatrix}\equiv
R_{\beta\alpha}\Phi_{\alpha},$ (31)

where the orthogonal matrices

$R_{1}(\beta_{1})=\begin{pmatrix}\cos\beta_{1}&\sin\beta_{1}&0&0\\\
-\sin\beta_{1}&\cos\beta_{1}&0&0\\\ 0&0&1&0\\\ 0&0&0&1\end{pmatrix},$ (32a)
$R_{2}(\beta_{2})=\begin{pmatrix}1&0&0&0\\\ 0&\cos\beta_{2}&\sin\beta_{2}&0\\\
0&-\sin\beta_{2}&\cos\beta_{2}&0\\\ 0&0&0&1\\\ \end{pmatrix},$ (32b)
$R_{3}(\beta_{3})=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\
0&0&\cos\beta_{3}&\sin\beta_{3}\\\ 0&0&-\sin\beta_{3}&\cos\beta_{3}\\\
\end{pmatrix},$ (32c)

where $\tan\beta_{1}=\frac{\sqrt{v_{2}^{2}+v_{3}^{2}+v_{4}^{2}}}{v_{1}}$,
$\tan\beta_{2}=\frac{\sqrt{v_{3}^{2}+v_{4}^{2}}}{v_{2}}$ and
$\tan\beta_{3}=\frac{v_{4}}{v_{3}}$, and
$H_{\beta}=(H_{\beta}^{+},(H_{\beta}^{0}+iH_{\beta}^{\text{odd}})/\sqrt{2})^{\text{T}}$.
This basis is chosen in such a way that only the neutral component of $H_{1}$
acquires a vacuum expectation value

$\displaystyle\langle
H_{1}^{0}\rangle=\sqrt{v_{1}^{2}+v_{2}^{2}+v_{3}^{2}+v_{4}^{2}}\equiv v,\hskip
14.22636pt$ $\displaystyle\langle H_{2}^{0}\rangle=0,\hskip 14.22636pt\langle
H_{3}^{0}\rangle=0,\hskip 14.22636pt\langle H_{4}^{0}\rangle=0.$ (33)

In this way
$\Phi_{\alpha}y_{ij}^{F\alpha}=y_{ij}^{F\alpha}R_{\alpha\beta}^{T}R_{\beta\gamma}\Phi_{\gamma}=\mathcal{Y}_{ij}^{F\beta}H_{\beta}$,
and $F=U,D,N,E$; where we have defined

$\displaystyle\mathcal{Y}^{F\beta}_{ij}=R_{\beta\alpha}y_{ij}^{F\alpha}.$ (34)

With these definitions, equation (29) becomes

$\displaystyle\mathcal{L}=-\bar{q}_{L}^{\prime
i}H_{\beta}\mathcal{Y}_{ij}^{D\beta}d^{\prime j}_{R}-\bar{q}_{L}^{\prime
i}\tilde{H}_{\beta}\mathcal{Y}_{ij}^{U\beta}u^{\prime
j}_{R}-\bar{\ell}_{L}^{\prime i}H_{\beta}\mathcal{Y}_{ij}^{E\beta}e^{\prime
j}_{R}-\bar{\ell}_{L}^{\prime
i}\tilde{H}_{\beta}\mathcal{Y}_{ij}^{N\beta}\nu^{\prime j}_{R}+\text{h.c}.$
(35)

It is necessary to rotate to the fermion mass eigenstates, i.e.,

$\displaystyle f_{L,R}=U^{F}_{L,R}f^{\prime}_{L,R},$ (36)

where the diagonalization matrices $U_{L,R}$ are defined below, in section C.
From the Lagrangian for the charged currents

$\displaystyle\mathcal{L}_{CC}=$
$\displaystyle-\frac{g}{\sqrt{2}}\bar{u}_{Li}^{\prime}\gamma^{\mu}d_{Li}^{\prime}W^{+}-\frac{g}{\sqrt{2}}\bar{e}_{Li}^{\prime}\gamma^{\mu}\nu_{Li}^{\prime}W^{-}+\text{h.c}$
$\displaystyle=$
$\displaystyle-\frac{g}{\sqrt{2}}\bar{u}_{Li}\gamma^{\mu}\left(V_{{}_{\text{CKM}}}\right)_{ij}d_{Lj}W^{+}-\frac{g}{\sqrt{2}}\bar{e}_{Li}\gamma^{\mu}\left(V_{{}_{\text{PMNS}}}\right)_{ij}\nu_{Lj}W^{-}+\text{h.c},$
(37)

it is possible to obtain the CKM
($V_{{}_{\text{CKM}}}=U^{U}_{L}U^{D\dagger}_{L}$) and PMNS
($V_{{}_{\text{PMNS}}}=U^{E}_{L}U^{\nu\dagger}_{L}$) mixing matrices by
rotating to the fermion mass eigenstates. In particular, we are interested in
the coupling of the axial neutral current to the axion in the mass
eigenstates.

$\displaystyle\mathcal{L}_{H^{0}}=$
$\displaystyle-\frac{1}{\sqrt{2}}\bar{d}_{L}^{\prime
i}H_{\beta}^{0}\mathcal{Y}_{ij}^{D\beta}d^{\prime
j}_{R}-\frac{1}{\sqrt{2}}\bar{u}_{L}^{\prime
i}H_{\beta}^{0*}\mathcal{Y}_{ij}^{U\beta}u^{\prime
j}_{R}-\frac{1}{\sqrt{2}}\bar{e}_{L}^{\prime
i}H_{\beta}^{0}\mathcal{Y}_{ij}^{E\beta}e^{\prime
j}_{R}-\frac{1}{\sqrt{2}}\bar{\nu}_{L}^{\prime
i}H_{\beta}^{0*}\mathcal{Y}_{ij}^{N\beta}\nu^{\prime j}_{R}+\text{h.c},$
$\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{2}}\bar{d}_{L}^{i}H_{\beta}^{0}Y_{ij}^{D\beta}d^{j}_{R}-\frac{1}{\sqrt{2}}\bar{u}_{L}^{i}H_{\beta}^{0*}Y_{ij}^{U\beta}u^{j}_{R}-\frac{1}{\sqrt{2}}\bar{e}_{L}^{i}H_{\beta}^{0}Y_{ij}^{E\beta}e^{j}_{R}-\frac{1}{\sqrt{2}}\bar{\nu}_{L}^{i}H_{\beta}^{0*}Y_{ij}^{N\beta}\nu^{j}_{R}+\text{h.c},$

where
$Y_{ij}^{F\beta}=\left(U^{F}_{L}\mathcal{Y}^{F\beta}U^{F\dagger}_{R}\right)_{ij}$.
In these expressions the mass functions in the interaction basis are:

$\displaystyle M^{D}_{ij}=\frac{v}{\sqrt{2}}\mathcal{Y}_{ij}^{D1},\hskip
14.22636ptM^{U}_{ij}=\frac{v}{\sqrt{2}}\mathcal{Y}_{ij}^{U1},\hskip
14.22636ptM^{E}_{ij}=\frac{v}{\sqrt{2}}\mathcal{Y}_{ij}^{E1},\hskip
14.22636ptM^{N}_{ij}=\frac{v}{\sqrt{2}}\mathcal{Y}_{ij}^{N1},$ (38)

where $v=\langle H_{1}^{0}\rangle$ is the Higgs vacuum expectation value.

##  Appendix B Scalar potential

As studied in Giraldo:2020hwl , the scalar sector requires four scalar
doublets $\phi^{\alpha}$ to reproduce the mass textures of the fermion sector
correctly, and two scalar singlets $S_{1}$ and $S_{2}$ that break the PQ
symmetry while generating a phenomenologically viable scalar mass spectrum.
The $S_{2}$ singlet also gives mass to the heavy quark. The most general
potential allowed by the PQ symmetry according to the charges established in
Table 2 is:

$\displaystyle V(\Phi,S_{i})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{4}\mu_{i}^{2}\Phi_{i}^{\dagger}\Phi_{i}+\sum_{k=1}^{2}\mu_{s_{k}}^{2}S_{k}^{*}S_{k}+\sum_{i=1}^{4}\lambda_{i}\left(\Phi_{i}^{\dagger}\Phi_{i}\right)^{2}$
(39) $\displaystyle+$
$\displaystyle\sum_{k=1}^{2}\lambda_{s_{k}}\left(S_{k}^{*}S_{k}\right)^{2}+\sum_{i=1}^{4}\sum_{k=1}^{2}\lambda_{is_{k}}\left(\Phi_{i}^{\dagger}\Phi_{i}\right)\left(S_{k}^{*}S_{k}\right)$
$\displaystyle+$
$\displaystyle\sum_{\underbrace{i,j=1}_{i<j}}^{4}\bigg{(}\lambda_{ij}\left(\Phi_{i}^{\dagger}\Phi_{i}\right)\left(\Phi_{j}^{\dagger}\Phi_{j}\right)+J_{ij}\left(\Phi_{i}^{\dagger}\Phi_{j}\right)\left(\Phi_{j}^{\dagger}\Phi_{i}\right)\bigg{)}$
$\displaystyle+$
$\displaystyle\lambda_{s_{1}s_{2}}\left(S_{1}^{*}S_{1}\right)\left(S_{2}^{*}S_{2}\right)$
$\displaystyle+$ $\displaystyle
K_{1}\left(\left(\Phi_{1}^{\dagger}\Phi_{2}\right)\left(\Phi_{3}^{\dagger}\Phi_{2}\right)+h.c.\right)$
$\displaystyle+$ $\displaystyle
K_{2}\left(\left(\Phi_{3}^{\dagger}\Phi_{4}\right)\left(\Phi_{3}^{\dagger}\Phi_{1}\right)+h.c.\right)$
$\displaystyle+$ $\displaystyle
F_{1}\left(\left(\Phi_{2}^{\dagger}\Phi_{3}\right)S_{1}+h.c.\right)$
$\displaystyle+$ $\displaystyle
F_{2}\left(\left(\Phi_{1}^{\dagger}\Phi_{2}\right)S_{1}+h.c.\right)$
$\displaystyle+$
$\displaystyle\frac{1}{2}\left(m_{\zeta_{S_{2}}}\right)^{2}_{\text{SB}}\zeta^{2}_{S_{2}}+\frac{1}{2}\left(m_{\xi_{S_{2}}}\right)^{2}_{\text{SB}}\xi^{2}_{S_{2}}.$

where the terms proportional to $F_{i}$ are allowed by the particular choice
of PQ charges and these couplings $F_{i}$ have units of mass. After
spontaneous symmetry breaking (SSB), the four Higgs doublets acquire VEVs that
give mass to all the SM particles. The scalar doublets and singlets are
written as follows:

$\displaystyle\Phi_{\alpha}=$
$\displaystyle\begin{pmatrix}\phi_{\alpha}^{+}\\\
\frac{v_{\alpha}+h_{\alpha}+i\eta_{\alpha}}{\sqrt{2}}\end{pmatrix},\hskip
28.45274pt\tilde{\Phi}_{\alpha}=i\sigma_{2}\Phi_{\alpha}^{*},\hskip
14.22636pt{\alpha=1,2,3,4,}$ $\displaystyle S_{i}=$
$\displaystyle\frac{v_{{}_{S_{i}}}+\xi_{S_{i}}+i\zeta_{S_{i}}}{\sqrt{2}};\hskip
28.45274pti=1,2,$ (40)

where the VEVs satisfy the following hierarchy: $v_{4}\ll v_{1},v_{2}\ll
v_{3}\ll v_{S_{1}}\sim v_{S_{2}}$. The scalar singlets $S_{1}$ and $S_{2}$
break the PQ symmetry at the high energy scale given by $v_{s_{1}}\approx
v_{s_{2}}$. The last two terms in equation (39) correspond to the soft-
breaking masses of the imaginary and the real parts of $S_{2}$, which are
generated at one loop in the Coleman-Weinberg potential from the interaction
term $\lambda_{Q}S_{2}\bar{Q}_{R}Q_{L}+\text{h.c.}$ Additionally, we choose
numerical values for the parameters of the potential (39) in order to obtain a
scalar sector mass spectrum consistent with the existing phenomenology. The
values of these parameters are:

$\displaystyle\lambda_{1}$ $\displaystyle=$
$\displaystyle\lambda_{2}=\lambda_{4}=\lambda_{s_{1}}=\lambda_{s_{2}}=\lambda_{s_{1}s_{2}}=1,\hskip
28.45274pt$ $\displaystyle\lambda_{3}$ $\displaystyle=$ $\displaystyle 0.463$
$\displaystyle\lambda_{ij}$ $\displaystyle=$ $\displaystyle 1\textup{ for any
}i,j,$ $\displaystyle\lambda_{js_{1}}$ $\displaystyle=$
$\displaystyle\lambda_{js_{2}}=1\textup{ for any }j,$ $\displaystyle J_{12}$
$\displaystyle=$ $\displaystyle J_{13}=J_{23}=J_{24}=-1,\ \ \textup{otherwise
}J_{ij}=1,$ $\displaystyle K_{1}$ $\displaystyle=$ $\displaystyle K_{2}=-1,$
$\displaystyle F_{1}$ $\displaystyle=$ $\displaystyle F_{2}=-1\text{GeV}.$
(41)

In particular, the value of $\lambda_{3}$ adjusts the SM Higgs mass. The
$v_{i}$ are determined from the SM fermion masses and the quark mass matrix
textures, Eq. (15). The VEV $v_{s_{1}}$ remains a free parameter; however,
this parameter is important for the axion physics due to the relationship
Giannotti:2017hny ,

$\displaystyle f_{a}=\frac{v_{s_{1}}}{2N}.$ (42)

In our calculations we took $v_{s_{1}}\approx v_{s_{2}}\approx 10^{6}$GeV. It
is important to emphasize that in our model, $f_{a}$ can take arbitrary
values; nevertheless, a small $f_{a}$ restricts $\epsilon$ (Eq. 13) to values
close to zero. Taking into account all these considerations, including Eq.
(B), the scalar mass spetrum (in GeV) is:

$\displaystyle\text{CP even}=\\{1.73\times 10^{6},1.\times 10^{6},6.54\times
10^{3},1.97\times 10^{3},$ $\displaystyle\hskip 51.21504pt1.09\times
10^{3},125\\},$ $\displaystyle\text{CP odd}=\\{6.54\times 10^{3},1.97\times
10^{3},1.09\times 10^{3},0,0,m_{\zeta_{S_{2}}}\\},$
$\displaystyle\text{Charged fields}=\\{6.54\times 10^{3},1.97\times
10^{3},1.11\times 10^{3},0\\}.$ (43)

The mass spectrum of the scalar fields is above the TeVs scale, except for the
SM Higgs, which is at 125 GeV. The pseudoscalar sector (CP odd fields) have
two massless eigenstates, the axion field and the Goldstone boson which is
absorbed by the longitudinal component of the SM $Z$ boson. A similar result
is obtained in the charged sector, where it is possible to identify the two
Goldstone bosons required to give mass to the SM $W^{\pm}$ fields.

##  Appendix C diagonalization matrices

To compare with physical quantities, it is necessary to rotate fields to the
mass eigenstates, i.e., $f_{L,R}=U^{F}_{L,R}f^{\prime}_{L,R}$, where the prime
symbol stands for the interaction basis. In our formalism the quark mass
matrices are Hermitian, so the right- and left-handed diagonalizing matrices
are identical; additionally, we establish that the eigenvalues of the second
family of quarks are negative in order to generate texture-zeros in some
diagonal terms of the mass matrices, as indicated in Branco:1999nb . This sign
is taken into account by introducing the identity matrix written as
$I_{2}I_{2}=1$ with $I_{2}=\text{diag}(1,-1,1)$, i.e.,

$\displaystyle M^{F}_{ij}=$
$\displaystyle\left(U^{F\dagger}\lambda^{F}U^{F}\right)_{ij}=\left(U^{F\dagger}_{L}m^{F}U^{F}_{R}\right)_{ij}=\frac{v}{\sqrt{2}}\mathcal{Y}^{F1}_{ij}=\frac{v}{\sqrt{2}}R_{1\alpha}y^{F\alpha}_{ij},$
(44)

where $\mathcal{Y}_{ij}^{F\beta}$ and $R_{\alpha\beta}$ were defined in
section A, $\lambda^{U,D}=\text{diag}(m_{u,d},-m_{c,s},m_{t,b})$ and
$m^{U,D}=\text{diag}(m_{u,d},m_{c,s},m_{t,b})$, with similar definitions in
the lepton sector, i.e.,
$\lambda^{N,E}=\text{diag}(m_{1,e},-m_{2,\mu},m_{3,\tau})$,
$m^{N,E}=\text{diag}(m_{1,e},m_{2,\mu},m_{3,\tau})$, and

$\displaystyle U_{L}^{F}=U^{F},\hskip 28.45274ptU_{R}^{F}=I_{2}U^{F},$ (45)

where the $U^{F}$ diagonalization matrices are defined below. It is important
to stress that the texture-zeros pattern in the matrix $\mathcal{Y}^{F1}_{ij}$
are identical to those in the original Yukawa couplings $y^{F\alpha}_{ij}$,
since the sum over $\alpha$ does not mix the $i,j$ indices. In fact, according
to equations (14) and (16),
$M^{F}=\frac{v_{\alpha}}{\sqrt{2}}y^{F\alpha}_{ij}=\frac{v}{\sqrt{2}}R_{1\alpha}y^{F\alpha}_{ij}$,
therefore $R_{1\alpha}=\frac{v_{\alpha}}{v}$. The diagonalization matrices
are:

$U^{U\dagger}=\begin{pmatrix}e^{i(\phi_{C_{u}}+\theta_{1u})}\sqrt{\frac{m_{c}m_{t}(A_{u}-m_{u})}{A_{u}(m_{c}+m_{u})(m_{t}-m_{u})}}&-e^{i(\phi_{C_{u}}+\theta_{2u})}\sqrt{\frac{(A_{u}+m_{c})m_{t}m_{u}}{A_{u}(m_{c}+m_{t})(m_{c}+m_{u})}}&e^{i(\phi_{C_{u}}+\theta_{3u})}\sqrt{\frac{m_{c}(m_{t}-A_{u})m_{u}}{A_{u}(m_{c}+m_{t})(m_{t}-m_{u})}}\\\
-e^{i(\phi_{B_{u}}+\theta_{1u})}\sqrt{\frac{(A_{u}+m_{c})(m_{t}-A_{u})m_{u}}{A_{u}(m_{c}+m_{u})(m_{t}-m_{u})}}&-e^{i(\phi_{B_{u}}+\theta_{2u})}\sqrt{\frac{m_{c}(m_{t}-A_{u})(A_{u}-m_{u})}{A_{u}(m_{c}+m_{t})(m_{c}+m_{u})}}&e^{i(\phi_{B_{u}}+\theta_{3u})}\sqrt{\frac{(A_{u}+m_{c})m_{t}(A_{u}-m_{u})}{A_{u}(m_{c}+m_{t})(m_{t}-m_{u})}}\\\
e^{i\theta_{1u}}\sqrt{\frac{m_{u}(A_{u}-m_{u})}{(m_{c}+m_{u})(m_{t}-m_{u})}}&e^{i\theta_{2u}}\sqrt{\frac{m_{c}(A_{u}+m_{c})}{(m_{c}+m_{t})(m_{c}+m_{u})}}&e^{i\theta_{3u}}\sqrt{\frac{m_{t}(m_{t}-A_{u})}{(m_{c}+m_{t})(m_{t}-m_{u})}}\end{pmatrix},$
(46)
$U^{D\dagger}=\begin{pmatrix}e^{i\theta_{1d}}\sqrt{\frac{m_{b}(m_{b}-m_{s})m_{s}}{(m_{b}-m_{d})(m_{d}+m_{s})(m_{b}+m_{d}-m_{s})}}&-e^{i\theta_{2d}}\sqrt{\frac{m_{b}(m_{b}+m_{d})m_{d}}{(m_{d}+m_{s})(m_{b}+m_{d}-m_{s})(m_{b}+m_{s})}}&\sqrt{\frac{m_{d}(m_{s}-m_{d})m_{s}}{(m_{b}-m_{d})(m_{b}+m_{d}-m_{s})(m_{b}+m_{s})}}\\\
e^{i\theta_{1d}}\sqrt{\frac{m_{d}(m_{b}-m_{s})}{(m_{b}-m_{d})(m_{d}+m_{s})}}&e^{i\theta_{2d}}\sqrt{\frac{(m_{b}+m_{d})m_{s}}{(m_{d}+m_{s})(m_{b}+m_{s})}}&\sqrt{\frac{m_{b}(m_{s}-m_{d})}{(m_{b}-m_{d})(m_{b}+m_{s})}}\\\
-e^{i\theta_{1d}}\sqrt{\frac{m_{d}(m_{b}+m_{d})(m_{s}-m_{d})}{(m_{b}-m_{d})(m_{d}+m_{s})(m_{b}+m_{d}-m_{s})}}&-e^{i\theta_{2d}}\sqrt{\frac{(m_{b}-m_{s})m_{s}(m_{s}-m_{d})}{(m_{d}+m_{s})(m_{b}+m_{d}-m_{s})(m_{b}+m_{s})}}&\sqrt{\frac{m_{b}(m_{b}+m_{d})(m_{b}-m_{s})}{(m_{b}-m_{d})(m_{b}+m_{d}-m_{s})(m_{b}+m_{s})}}\end{pmatrix},$
(47)

where $\theta_{1u},\theta_{2u},\theta_{3u},\theta_{1d}$ and $\theta_{2d}$ are
arbitrary phases (a third phase for the diagonalization matrix (47) can be
absorbed by the remaining phases) that are useful for conforming to the
$V_{\text{CKM}}=U^{U}_{L}\,U^{D{\dagger}}_{L}$ matrix convention. Taking as
input the SM parameters at the $Z$ pole, the best fit values are given in
Table 4.

$\theta_{1u}$ | $\theta_{2u}$ | $\theta_{3u}$ | $\theta_{1d}$ | $\theta_{2d}$ | $\phi_{C_{u}}$ | $\phi_{B_{u}}$
---|---|---|---|---|---|---
$-2.84403$ | 1.85606 | $-0.00461668$ | 1.93013 | $-0.976639$ | $-1.49697$ | 0.301461
$A_{u}$ | $m_{u}$ | $m_{c}$ | $m_{t}$ | $m_{d}$ | $m_{s}$ | $m_{b}$
1690.29 MeV | 1.2684 MeV | 633.197 MeV | 171268 MeV | 3.14751 MeV | 56.1169 MeV | 2910.01 MeV
Table 4: Best-fit point of the mass matrix parameters with respect to
experimental data for the masses and mixing angles of the quark sector at the
Z pole.

Similarly, in the lepton sector, the diagonalization matrices of the mass
matrices (7) are:

$\displaystyle U^{N{\dagger}}$
$\displaystyle=\begin{pmatrix}e^{i(\theta_{1\nu}+c_{\nu})}\,\sqrt{\frac{m_{2}m_{3}(A_{\nu}-m_{1})}{A_{\nu}(m_{2}+m_{1})(m_{3}-m_{1})}}&-e^{i(\theta_{2\nu}+c_{\nu})}\sqrt{\frac{m_{1}m_{3}(m_{2}+A_{\nu})}{A_{\nu}(m_{2}+m_{1})(m_{3}+m_{2})}}&e^{i(\theta_{3\nu}+c_{\nu})}\,\sqrt{\frac{m_{1}m_{2}(m_{3}-A_{\nu})}{A_{\nu}(m_{3}-m_{1})(m_{3}+m_{2})}}\\\
&&&\\\\[-5.69054pt]
e^{i\theta_{1\nu}}\sqrt{\frac{m_{1}(A_{\nu}-m_{1})}{(m_{1}+m_{2})(m_{3}-m_{1})}}&e^{i\theta_{2\nu}}\sqrt{\frac{m_{2}(A_{\nu}+m_{2})}{(m_{2}+m_{1})(m_{3}+m_{2})}}&e^{i\theta_{3\nu}}\,\sqrt{\frac{m_{3}(m_{3}-A_{\nu})}{(m_{3}-m_{1})(m_{3}+m_{2})}}\\\
&&&\\\\[-5.69054pt]
-e^{i(\theta_{1\nu}-b_{\nu})}\sqrt{\frac{m_{1}(A_{\nu}+m_{2})(m_{3}-A_{\nu})}{A_{\nu}(m_{1}+m_{2})(m_{3}-m_{1})}}&-e^{i(\theta_{2\nu}-b_{\nu})}\sqrt{\frac{m_{2}(A_{\nu}-m_{1})(m_{3}-A_{\nu})}{A_{\nu}(m_{2}+m_{1})(m_{3}+m_{2})}}&e^{i(\theta_{3\nu}-b_{\nu})}\,\sqrt{\frac{m_{3}(A_{\nu}-m_{1})(A_{\nu}+m_{2})}{A_{\nu}(m_{3}-m_{1})(m_{3}+m_{2})}}\end{pmatrix},$
(48) $\displaystyle U^{E{\dagger}}$
$\displaystyle=\begin{pmatrix}e^{i\theta_{1\ell}}\,\sqrt{\frac{m_{\mu}m_{\tau}(m_{\tau}-m_{\mu})}{(m_{e}-m_{\mu}+m_{\tau})(m_{\mu}+m_{e})(m_{\tau}-m_{e})}}&-e^{i\theta_{2\ell}}\sqrt{\frac{m_{e}m_{\tau}(m_{e}+m_{\tau})}{(m_{e}-m_{\mu}+m_{\tau})(m_{\mu}+m_{e})(m_{\tau}+m_{\mu})}}&\sqrt{\frac{m_{e}m_{\mu}(m_{\mu}-m_{e})}{(m_{e}-m_{\mu}+m_{\tau})(m_{\tau}-m_{e})(m_{\tau}+m_{\mu})}}\\\
&&&\\\\[-5.69054pt]
e^{i\theta_{1\ell}}\sqrt{\frac{m_{e}(m_{\tau}-m_{\mu})}{(m_{\mu}+m_{e})(m_{\tau}-m_{e})}}&e^{i\theta_{2\ell}}\sqrt{\frac{m_{\mu}(m_{e}+m_{\tau})}{(m_{\mu}+m_{e})(m_{\tau}+m_{\mu})}}&\sqrt{\frac{m_{\tau}(m_{\mu}-m_{e})}{(m_{\tau}-m_{e})(m_{\tau}+m_{\mu})}}\\\
&&&\\\\[-5.69054pt]
-e^{i\theta_{1\ell}}\sqrt{\frac{m_{e}(m_{e}+m_{\tau})(m_{\mu}-m_{e})}{(m_{e}-m_{\mu}+m_{\tau})(m_{\mu}+m_{e})(m_{\tau}-m_{e})}}&-e^{i\theta_{2\ell}}\sqrt{\frac{m_{\mu}(m_{\tau}-m_{\mu})(m_{\mu}-m_{e})}{(m_{e}-m_{\mu}+m_{\tau})(m_{\mu}+m_{e})(m_{\tau}+m_{\mu})}}&\sqrt{\frac{m_{\tau}(m_{\tau}-m_{\mu})(m_{e}+m_{\tau})}{(m_{e}-m_{\mu}+m_{\tau})(m_{\tau}-m_{e})(m_{\tau}+m_{\mu})}}\end{pmatrix},$
(49)

where $\theta_{1\ell}$, $\theta_{2\ell}$, $\theta_{1\nu}$, $\theta_{2\nu}$,
$\theta_{3\nu}$ are necessary phases in order to adjust to the established
convention for the PMNS mixing matrix Esteban:2018azc  111 NuFIT collaboration
(http://www.nu-fit.org/?q=node/211)(with SK atmospheric data).; and $c_{\nu}$
and $b_{\nu}$ are the phases of $C_{\nu}$ and $B_{\nu}$ in the neutral mass
matrix $M^{N}$ in Eq. (7). The best fit values for these quantities are shown
in Table 5.

$\theta_{1\ell}$ | $\theta_{2\ell}$ | $\theta_{1\nu}$ | $\theta_{2\nu}$ | $\theta_{3\nu}$ | $c_{\nu}$ | $b_{\nu}$
---|---|---|---|---|---|---
0.154895 | 2.01797 | $-0.835504$ | 2.21169 | 1.81786 | 1.01608 | 2.03726
$A_{\nu}\,(\text{eV})$ | $m_{e}\,(\text{MeV})$ | $m_{\mu}\,(\text{MeV})$ | $m_{\tau}\,(\text{MeV})$ | $m_{1}\,(\text{eV})$ | $m_{2}\,(\text{eV})$ | $m_{3}\,(\text{eV})$
0.0251821 | 0.5109989461 | 105.6583745 | 1776.86 | 0.00353647 | 0.00929552 | 0.0504034
Table 5: Best fit values.

##  Appendix D Axion decay into photons

In the SM, $B^{\mu}=\cos\theta_{W}A^{\mu}-\sin\theta_{W}Z^{\mu}$ and
$W^{3\mu}=\sin\theta_{W}A^{\mu}+\cos\theta_{W}Z^{\mu}$, where $A^{\mu}$ and
$Z^{\mu}$ are the SM fields for the photon and $Z$ gauge bosons, replacing
these expresions in Eq. (V) we obtain

$\displaystyle\mathcal{L}$ $\displaystyle\supset-
c_{1}^{\text{eff}}\frac{\alpha_{1}}{8\pi}\frac{a}{\Lambda_{\text{PQ}}}B_{\mu\nu}\tilde{B}^{\mu\nu}-c_{2}^{\text{eff}}\frac{\alpha_{2}}{8\pi}\frac{a}{\Lambda_{\text{PQ}}}W_{\mu\nu}^{3}\tilde{W}^{3\mu\nu}$
$\displaystyle=-\frac{\alpha}{8\pi}(c_{1}^{\text{eff}}+c_{2}^{\text{eff}})\frac{a}{\Lambda_{\text{PQ}}}F_{\mu\nu}\tilde{F}^{\mu\nu}$
$\displaystyle-\frac{\alpha}{8\pi
c_{W}^{2}s_{W}^{2}}(s^{4}_{W}c_{1}^{\text{eff}}+c^{4}_{W}c_{2}^{\text{eff}})\frac{a}{f_{a}}Z_{\mu\nu}\tilde{Z}^{\mu\nu}$
$\displaystyle-\frac{2\alpha}{8\pi
c_{W}s_{W}}(c^{2}_{W}c_{2}^{\text{eff}}-c^{2}_{W}c_{1}^{\text{eff}})\frac{a}{\Lambda_{\text{PQ}}}F_{\mu\nu}\tilde{Z}^{\mu\nu}$
$\displaystyle=e^{2}C_{\gamma\gamma}\frac{a}{\Lambda_{\text{PQ}}}F_{\mu\nu}\tilde{F}^{\mu\nu}+\frac{e^{2}C_{ZZ}}{c_{W}^{2}s_{W}^{2}}\frac{a}{\Lambda_{\text{PQ}}}Z_{\mu\nu}\tilde{Z}^{\mu\nu}+\frac{2e^{2}C_{\gamma
Z}}{c_{W}s_{W}}\frac{a}{\Lambda_{\text{PQ}}}F_{\mu\nu}\tilde{Z}^{\mu\nu}$ (50)
$\displaystyle c_{1}^{\text{eff}}$ $\displaystyle=c_{1}-\frac{1}{3}\Sigma
q+\frac{8}{3}\Sigma u+\frac{2}{3}\Sigma d-\Sigma l+2\Sigma e$ (51)
$\displaystyle c_{2}^{\text{eff}}$ $\displaystyle=c_{2}-3\Sigma q-\Sigma l$
(52)

where $\Sigma f\equiv f_{1}+f_{2}+f_{3}$ is the sum of the PQ charges of the
three families. There are similar definitions for the interaction of the axion
with the gluons

$\displaystyle-
c_{3}^{\text{eff}}\frac{\alpha_{3}}{8\pi}\frac{a}{\Lambda_{\text{PQ}}}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu}=g_{s}^{2}C_{GG}\frac{a}{\Lambda_{\text{PQ}}}G_{\mu\nu}^{a}\tilde{G}^{a\mu\nu},$
(53)

where $c_{3}^{\text{eff}}=c_{3}-2\Sigma q+\Sigma u+\Sigma d-A_{Q}$, in our
particular case $c_{i}=0$. In axion phenomenology, it is usual to define

$\displaystyle C_{\gamma\gamma}=$
$\displaystyle-\frac{1}{32\pi^{2}}(c_{1}^{\text{eff}}+c_{2}^{\text{eff}}),\hskip
39.83368ptC_{ZZ}=-\frac{1}{32\pi^{2}}(s_{W}^{4}c_{1}^{\text{eff}}+c_{W}^{4}c_{2}^{\text{eff}}),$
$\displaystyle C_{\gamma Z}=$
$\displaystyle-\frac{1}{32\pi^{2}}(c_{W}^{2}c_{2}^{\text{eff}}-c_{W}^{2}c_{1}^{\text{eff}}),\hskip
14.22636ptC_{GG}=-\frac{1}{32\pi^{2}}c_{3}^{\text{eff}}.$ (54)

The decay widths of an axion decaying in two photons and a $Z$ decaying in an
axion and a photon are Bauer:2017ris

$\displaystyle\Gamma(a\rightarrow\gamma\gamma)=$
$\displaystyle\frac{4\pi\alpha^{2}m_{a}^{3}}{\Lambda^{2}_{\text{PQ}}}\lvert
C_{\gamma\gamma}^{\text{eff}}\rvert^{2},$
$\displaystyle\Gamma(Z\rightarrow\gamma a)=$
$\displaystyle\frac{8\pi\alpha(m_{Z})m_{Z}^{3}}{3s_{W}^{2}c_{W}^{2}\Lambda^{2}_{\text{PQ}}}\lvert
C_{\gamma
Z}^{\text{eff}}\rvert^{2}\left(1-\frac{m_{a}^{2}}{m_{Z}^{2}}\right)^{3}.$ (55)

Another possible decay channel of the axion in two photons is due to the
mixing between the axion and the pion since the latter can decay in two
photons, this decay mode generates an additional correction that only depends
on the couplings of the axion to the gluons Alonso-Alvarez:2018irt

$\displaystyle C_{\gamma\gamma}^{\text{eff}}=$
$\displaystyle-\frac{c_{3}^{\text{eff}}}{32\pi^{2}}\left(\frac{c_{1}^{\text{eff}}+c_{2}^{\text{eff}}}{c_{3}^{\text{eff}}}-2.03\right),$
$\displaystyle C_{\gamma Z}^{\text{eff}}=$
$\displaystyle-\frac{c_{3}^{\text{eff}}}{32\pi^{2}}\left(\frac{c_{W}^{2}c_{2}^{\text{eff}}-c_{W}^{2}c_{1}^{\text{eff}}}{c_{3}^{\text{eff}}}-0.74/2\right).$
(56)

It is usual to define $\Lambda_{\text{PQ}}=\lvert c_{3}^{\text{eff}}\rvert
f_{a}$.

$\displaystyle\frac{E}{N}=\frac{c_{1}^{\text{eff}}+c_{2}^{\text{eff}}}{c_{3}^{\text{eff}}}.$
(57)

The axion-photon interaction is given by

$\displaystyle
g_{a\gamma\gamma}=\frac{4e^{2}C_{\gamma\gamma}^{\text{eff}}}{\Lambda_{PQ}}=-\frac{\alpha}{2\pi
f_{PQ}}\left(\frac{E}{N}-2.03\right)$ (58)

where $\alpha=\frac{e^{2}}{4\pi}$. Due to the gluon-axion interaction, the
axion gets a mass term, which is described at low energies as an axion-pion
interaction diCortona:2015ldu

$\displaystyle m_{a}=5.7(7)\mu
eV\left(\frac{10^{12}\text{GeV}}{f_{a}}\right).$ (59)

## References

  * (1) G. Aad et al., “Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC,” Phys. Lett. B, vol. 716, pp. 1–29, 2012.
  * (2) S. Chatrchyan et al., “Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC,” Phys. Lett. B, vol. 716, pp. 30–61, 2012.
  * (3) X.-G. He, T. Li, X.-Q. Li, J. Tandean, and H.-C. Tsai, “Constraints on Scalar Dark Matter from Direct Experimental Searches,” Phys. Rev. D, vol. 79, p. 023521, 2009.
  * (4) B. Grzadkowski and P. Osland, “Tempered Two-Higgs-Doublet Model,” Phys. Rev. D, vol. 82, p. 125026, 2010.
  * (5) H. E. Logan, “Dark matter annihilation through a lepton-specific Higgs boson,” Phys. Rev. D, vol. 83, p. 035022, 2011.
  * (6) M. S. Boucenna and S. Profumo, “Direct and Indirect Singlet Scalar Dark Matter Detection in the Lepton-Specific two-Higgs-doublet Model,” Phys. Rev. D, vol. 84, p. 055011, 2011.
  * (7) X.-G. He, B. Ren, and J. Tandean, “Hints of Standard Model Higgs Boson at the LHC and Light Dark Matter Searches,” Phys. Rev. D, vol. 85, p. 093019, 2012.
  * (8) Y. Bai, V. Barger, L. L. Everett, and G. Shaughnessy, “Two-Higgs-doublet-portal dark-matter model: LHC data and Fermi-LAT 135 GeV line,” Phys. Rev. D, vol. 88, p. 015008, 2013.
  * (9) X.-G. He and J. Tandean, “Low-Mass Dark-Matter Hint from CDMS II, Higgs Boson at the LHC, and Darkon Models,” Phys. Rev. D, vol. 88, p. 013020, 2013\.
  * (10) Y. Cai and T. Li, “Singlet dark matter in a type II two Higgs doublet model,” Phys. Rev. D, vol. 88, no. 11, p. 115004, 2013.
  * (11) L. Wang and X.-F. Han, “A simplified 2HDM with a scalar dark matter and the galactic center gamma-ray excess,” Phys. Lett. B, vol. 739, pp. 416–420, 2014.
  * (12) A. Drozd, B. Grzadkowski, J. F. Gunion, and Y. Jiang, “Extending two-Higgs-doublet models by a singlet scalar field - the Case for Dark Matter,” JHEP, vol. 11, p. 105, 2014.
  * (13) R. Campbell, S. Godfrey, H. E. Logan, A. D. Peterson, and A. Poulin, “Implications of the observation of dark matter self-interactions for singlet scalar dark matter,” Phys. Rev. D, vol. 92, no. 5, p. 055031, 2015\.  [Erratum: Phys.Rev.D 101, 039905 (2020)].
  * (14) S. von Buddenbrock, N. Chakrabarty, A. S. Cornell, D. Kar, M. Kumar, T. Mandal, B. Mellado, B. Mukhopadhyaya, R. G. Reed, and X. Ruan, “Phenomenological signatures of additional scalar bosons at the LHC,” Eur. Phys. J. C, vol. 76, no. 10, p. 580, 2016.
  * (15) M. Carena, H. E. Haber, I. Low, N. R. Shah, and C. E. M. Wagner, “Alignment limit of the NMSSM Higgs sector,” Phys. Rev. D, vol. 93, no. 3, p. 035013, 2016.
  * (16) G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, “Theory and phenomenology of two-Higgs-doublet models,” Phys. Rept., vol. 516, pp. 1–102, 2012.
  * (17) J. McDonald, “Gauge singlet scalars as cold dark matter,” Phys. Rev. D, vol. 50, pp. 3637–3649, 1994.
  * (18) L. Lopez Honorez, E. Nezri, J. F. Oliver, and M. H. G. Tytgat, “The Inert Doublet Model: An Archetype for Dark Matter,” JCAP, vol. 02, p. 028, 2007\.
  * (19) R. Barbieri, L. J. Hall, and V. S. Rychkov, “Improved naturalness with a heavy Higgs: An Alternative road to LHC physics,” Phys. Rev. D, vol. 74, p. 015007, 2006.
  * (20) L. Lopez Honorez and C. E. Yaguna, “The inert doublet model of dark matter revisited,” JHEP, vol. 09, p. 046, 2010.
  * (21) S. L. Glashow, “Partial Symmetries of Weak Interactions,” Nucl. Phys., vol. 22, pp. 579–588, 1961.
  * (22) S. Weinberg, “A Model of Leptons,” Phys. Rev. Lett., vol. 19, pp. 1264–1266, 1967.
  * (23) A. Salam, “Weak and Electromagnetic Interactions,” Conf. Proc. C, vol. 680519, pp. 367–377, 1968.
  * (24) R. Davis, Jr., D. S. Harmer, and K. C. Hoffman, “Search for neutrinos from the sun,” Phys. Rev. Lett., vol. 20, pp. 1205–1209, 1968.
  * (25) J. N. Bahcall, “Solving the mystery of the missing neutrinos,” 6 2004.
  * (26) B. Pontecorvo, “Inverse beta processes and nonconservation of lepton charge,” Zh. Eksp. Teor. Fiz., vol. 34, p. 247, 1957.
  * (27) L. Wolfenstein, “Neutrino Oscillations in Matter,” Phys. Rev. D, vol. 17, pp. 2369–2374, 1978.
  * (28) S. P. Mikheyev and A. Y. Smirnov, “Resonance Amplification of Oscillations in Matter and Spectroscopy of Solar Neutrinos,” Sov. J. Nucl. Phys., vol. 42, pp. 913–917, 1985.
  * (29) F. Kaether, W. Hampel, G. Heusser, J. Kiko, and T. Kirsten, “Reanalysis of the GALLEX solar neutrino flux and source experiments,” Phys. Lett. B, vol. 685, pp. 47–54, 2010.
  * (30) B. T. Cleveland, T. Daily, R. Davis, Jr., J. R. Distel, K. Lande, C. K. Lee, P. S. Wildenhain, and J. Ullman, “Measurement of the solar electron neutrino flux with the Homestake chlorine detector,” Astrophys. J., vol. 496, pp. 505–526, 1998.
  * (31) B. Aharmim et al., “Combined Analysis of all Three Phases of Solar Neutrino Data from the Sudbury Neutrino Observatory,” Phys. Rev. C, vol. 88, p. 025501, 2013.
  * (32) G. Bellini et al., “Measurement of the solar 8B neutrino rate with a liquid scintillator target and 3 MeV energy threshold in the Borexino detector,” Phys. Rev. D, vol. 82, p. 033006, 2010.
  * (33) G. Bellini, J. Benziger, D. Bick, G. Bonfini, D. Bravo, B. Caccianiga, L. Cadonati, F. Calaprice, A. Caminata, P. Cavalcante, A. Chavarria, A. Chepurnov, D. D’Angelo, S. Davini, A. Derbin, A. Empl, A. Etenko, K. Fomenko, D. Franco, F. Gabriele, C. Galbiati, S. Gazzana, C. Ghiano, M. Giammarchi, M. Göger-Neff, A. Goretti, M. Gromov, C. Hagner, E. Hungerford, A. Ianni, A. Ianni, V. Kobychev, D. Korablev, G. Korga, D. Kryn, M. Laubenstein, B. Lehnert, T. Lewke, E. Litvinovich, F. Lombardi, P. Lombardi, L. Ludhova, G. Lukyanchenko, I. Machulin, S. Manecki, W. Maneschg, S. Marcocci, Q. Meindl, E. Meroni, M. Meyer, L. Miramonti, M. Misiaszek, M. Montuschi, P. Mosteiro, V. Muratova, L. Oberauer, M. Obolensky, F. Ortica, K. Otis, M. Pallavicini, L. Papp, L. Perasso, A. Pocar, G. Ranucci, A. Razeto, A. Re, A. Romani, N. Rossi, R. Saldanha, C. Salvo, S. Schönert, H. Simgen, M. Skorokhvatov, O. Smirnov, A. Sotnikov, S. Sukhotin, Y. Suvorov, R. Tartaglia, G. Testera, D. Vignaud, R. B. Vogelaar, F. von Feilitzsch, H. Wang, J. Winter, M. Wojcik, A. Wright, M. Wurm, O. Zaimidoroga, S. Zavatarelli, K. Zuber, G. Zuzel, and B. Collaboration, “Neutrinos from the primary proton–proton fusion process in the sun,” Nature, vol. 512, pp. 383–386, Aug 2014.
  * (34) J. Hosaka et al., “Solar neutrino measurements in super-Kamiokande-I,” Phys. Rev. D, vol. 73, p. 112001, 2006.
  * (35) J. Hosaka et al., “Three flavor neutrino oscillation analysis of atmospheric neutrinos in Super-Kamiokande,” Phys. Rev. D, vol. 74, p. 032002, 2006.
  * (36) M. G. Aartsen et al., “Determining neutrino oscillation parameters from atmospheric muon neutrino disappearance with three years of IceCube DeepCore data,” Phys. Rev. D, vol. 91, no. 7, p. 072004, 2015.
  * (37) A. Gando et al., “Constraints on $\theta_{13}$ from A Three-Flavor Oscillation Analysis of Reactor Antineutrinos at KamLAND,” Phys. Rev. D, vol. 83, p. 052002, 2011.
  * (38) M. Apollonio et al., “Limits on neutrino oscillations from the CHOOZ experiment,” Phys. Lett. B, vol. 466, pp. 415–430, 1999.
  * (39) A. Piepke, “Final results from the Palo Verde neutrino oscillation experiment,” Prog. Part. Nucl. Phys., vol. 48, pp. 113–121, 2002.
  * (40) F. P. An et al., “New Measurement of Antineutrino Oscillation with the Full Detector Configuration at Daya Bay,” Phys. Rev. Lett., vol. 115, no. 11, p. 111802, 2015.
  * (41) S.-B. Kim, “Measurement of neutrino mixing angle $\Theta_{13}$ and mass difference $\Delta m^{2}_{ee}$ from reactor antineutrino disappearance in the RENO experiment,” Nucl. Phys. B, vol. 908, pp. 94–115, 2016.
  * (42) J. Kopp, P. A. N. Machado, M. Maltoni, and T. Schwetz, “Sterile Neutrino Oscillations: The Global Picture,” JHEP, vol. 05, p. 050, 2013.
  * (43) P. Adamson et al., “Measurement of Neutrino and Antineutrino Oscillations Using Beam and Atmospheric Data in MINOS,” Phys. Rev. Lett., vol. 110, no. 25, p. 251801, 2013.
  * (44) K. Abe et al., “Observation of Electron Neutrino Appearance in a Muon Neutrino Beam,” Phys. Rev. Lett., vol. 112, p. 061802, 2014.
  * (45) P. Adamson et al., “First measurement of electron neutrino appearance in NOvA,” Phys. Rev. Lett., vol. 116, no. 15, p. 151806, 2016.
  * (46) J. Schechter and J. W. F. Valle, “Neutrino Masses in SU(2) x U(1) Theories,” Phys. Rev. D, vol. 22, p. 2227, 1980.
  * (47) E. Ma and U. Sarkar, “Neutrino masses and leptogenesis with heavy Higgs triplets,” Phys. Rev. Lett., vol. 80, pp. 5716–5719, 1998.
  * (48) R. Foot, H. Lew, X. G. He, and G. C. Joshi, “Seesaw Neutrino Masses Induced by a Triplet of Leptons,” Z. Phys. C, vol. 44, p. 441, 1989.
  * (49) R. N. Mohapatra and J. W. F. Valle, “Neutrino Mass and Baryon Number Nonconservation in Superstring Models,” Phys. Rev. D, vol. 34, p. 1642, 1986.
  * (50) M. C. Gonzalez-Garcia and J. W. F. Valle, “Fast Decaying Neutrinos and Observable Flavor Violation in a New Class of Majoron Models,” Phys. Lett. B, vol. 216, pp. 360–366, 1989.
  * (51) G. C. Branco, D. Emmanuel-Costa, M. N. Rebelo, and P. Roy, “Four Zero Neutrino Yukawa Textures in the Minimal Seesaw Framework,” Phys. Rev. D, vol. 77, p. 053011, 2008.
  * (52) S. Choubey, W. Rodejohann, and P. Roy, “Phenomenological consequences of four zero neutrino Yukawa textures,” Nucl. Phys. B, vol. 808, pp. 272–291, 2009.  [Erratum: Nucl.Phys.B 818, 136–136 (2009)].
  * (53) B. Adhikary, A. Ghosal, and P. Roy, “Neutrino Masses, Cosmological Bound and Four Zero Yukawa Textures,” Mod. Phys. Lett. A, vol. 26, pp. 2427–2435, 2011.
  * (54) B. Adhikary, A. Ghosal, and P. Roy, “Baryon asymmetry from leptogenesis with four zero neutrino Yukawa textures,” JCAP, vol. 01, p. 025, 2011.
  * (55) G. C. Branco, M. N. Rebelo, and J. I. Silva-Marcos, “Leptogenesis, Yukawa textures and weak basis invariants,” Phys. Lett. B, vol. 633, pp. 345–354, 2006.
  * (56) J.-w. Mei, “Running neutrino masses, leptonic mixing angles and CP-violating phases: From M(Z) to Lambda(GUT),” Phys. Rev. D, vol. 71, p. 073012, 2005\.
  * (57) C. Hagedorn, J. Kersten, and M. Lindner, “Stability of texture zeros under radiative corrections in see-saw models,” Phys. Lett. B, vol. 597, pp. 63–72, 2004.
  * (58) S. Antusch, J. Kersten, M. Lindner, M. Ratz, and M. A. Schmidt, “Running neutrino mass parameters in see-saw scenarios,” JHEP, vol. 03, p. 024, 2005.
  * (59) P. O. Ludl and W. Grimus, “A complete survey of texture zeros in the lepton mass matrices,” JHEP, vol. 07, p. 090, 2014.  [Erratum: JHEP 10, 126 (2014)].
  * (60) P. M. Ferreira and L. Lavoura, “New textures for the lepton mass matrices,” Nucl. Phys. B, vol. 891, pp. 378–400, 2015.
  * (61) P. O. Ludl and W. Grimus, “A complete survey of texture zeros in general and symmetric quark mass matrices,” Phys. Lett. B, vol. 744, pp. 38–42, 2015\.
  * (62) W. Grimus, A. S. Joshipura, L. Lavoura, and M. Tanimoto, “Symmetry realization of texture zeros,” Eur. Phys. J. C, vol. 36, pp. 227–232, 2004.
  * (63) R. M. Fonseca and W. Grimus, “Classification of lepton mixing matrices from finite residual symmetries,” JHEP, vol. 09, p. 033, 2014.
  * (64) Y. Giraldo, R. Martinez, E. Rojas, and J. C. Salazar, “Flavored axions and the flavor problem,” Eur. Phys. J. C, vol. 82, no. 12, p. 1131, 2022.
  * (65) F. Björkeroth, L. Di Luzio, F. Mescia, and E. Nardi, “$U(1)$ flavour symmetries as Peccei-Quinn symmetries,” JHEP, vol. 02, p. 133, 2019.
  * (66) T. P. Cheng and M. Sher, “Mass Matrix Ansatz and Flavor Nonconservation in Models with Multiple Higgs Doublets,” Phys. Rev. D, vol. 35, p. 3484, 1987\.
  * (67) K. Matsuda and H. Nishiura, “Can four-zero-texture mass matrix model reproduce the quark and lepton mixing angles and CP violating phases?,” Phys. Rev. D, vol. 74, p. 033014, 2006.
  * (68) A. E. Carcamo Hernandez, R. Martinez, and J. A. Rodriguez, “Different kind of textures of Yukawa coupling matrices in the two Higgs doublet model type III,” Eur. Phys. J. C, vol. 50, pp. 935–948, 2007.
  * (69) P. Langacker and M. Plumacher, “Flavor changing effects in theories with a heavy $Z^{\prime}$ boson with family nonuniversal couplings,” Phys. Rev. D, vol. 62, p. 013006, 2000.
  * (70) K. Leroux and D. London, “Flavor changing neutral currents and leptophobic $Z^{\prime}$ gauge bosons,” Phys. Lett. B, vol. 526, pp. 97–103, 2002.
  * (71) S. F. King, S. Moretti, and R. Nevzorov, “Theory and phenomenology of an exceptional supersymmetric standard model,” Phys. Rev. D, vol. 73, p. 035009, 2006.
  * (72) C. Alvarado, R. Martinez, and F. Ochoa, “Quark mass hierarchy in 3-3-1 models,” Phys. Rev. D, vol. 86, p. 025027, 2012.
  * (73) L. E. Ibanez and G. G. Ross, “Fermion masses and mixing angles from gauge symmetries,” Phys. Lett. B, vol. 332, pp. 100–110, 1994.
  * (74) P. Binetruy and P. Ramond, “Yukawa textures and anomalies,” Phys. Lett. B, vol. 350, pp. 49–57, 1995.
  * (75) Y. Nir, “Gauge unification, Yukawa hierarchy and the mu problem,” Phys. Lett. B, vol. 354, pp. 107–110, 1995.
  * (76) V. Jain and R. Shrock, “Models of fermion mass matrices based on a flavor dependent and generation dependent U(1) gauge symmetry,” Phys. Lett. B, vol. 352, pp. 83–91, 1995.
  * (77) E. Dudas, S. Pokorski, and C. A. Savoy, “Yukawa matrices from a spontaneously broken Abelian symmetry,” Phys. Lett. B, vol. 356, pp. 45–55, 1995.
  * (78) F. Pisano and V. Pleitez, “An SU(3) x U(1) model for electroweak interactions,” Phys. Rev. D, vol. 46, pp. 410–417, 1992.
  * (79) P. H. Frampton, “Chiral dilepton model and the flavor question,” Phys. Rev. Lett., vol. 69, pp. 2889–2891, 1992.
  * (80) R. Foot, H. N. Long, and T. A. Tran, “$SU(3)_{L}\otimes U(1)_{N}$ and $SU(4)_{L}\otimes U(1)_{N}$ gauge models with right-handed neutrinos,” Phys. Rev. D, vol. 50, no. 1, pp. R34–R38, 1994.
  * (81) H. N. Long, “The 331 model with right handed neutrinos,” Phys. Rev. D, vol. 53, pp. 437–445, 1996.
  * (82) H. N. Long, “SU(3)-L x U(1)-N model for right-handed neutrino neutral currents,” Phys. Rev. D, vol. 54, pp. 4691–4693, 1996.
  * (83) H. N. Long, “Scalar sector of the 3 3 1 model with three Higgs triplets,” Mod. Phys. Lett. A, vol. 13, pp. 1865–1874, 1998.
  * (84) R. A. Diaz, R. Martinez, and J. A. Rodriguez, “A New supersymmetric SU(3)(L) x U(1)(X) gauge model,” Phys. Lett. B, vol. 552, pp. 287–292, 2003.
  * (85) R. A. Diaz, R. Martinez, and F. Ochoa, “The Scalar sector of the SU(3)(c) x SU(3)(L) x U(1)(X) model,” Phys. Rev. D, vol. 69, p. 095009, 2004.
  * (86) R. A. Diaz, R. Martinez, and F. Ochoa, “SU(3)(c) x SU(3)(L) x U(1)(X) models for beta arbitrary and families with mirror fermions,” Phys. Rev. D, vol. 72, p. 035018, 2005.
  * (87) F. Ochoa and R. Martinez, “Family dependence in SU(3)(c) x SU(3)(L) x U(1)(X) models,” Phys. Rev. D, vol. 72, p. 035010, 2005.
  * (88) S. M. Barr, “Light Fermion Mass Hierarchy and Grand Unification,” Phys. Rev. D, vol. 21, p. 1424, 1980.
  * (89) T. R. Taylor and G. Veneziano, “Quenching the Cosmological Constant,” Phys. Lett. B, vol. 228, pp. 311–316, 1989.
  * (90) S. M. Barr, “A Simple and Predictive Model for Quark and Lepton Masses,” Phys. Rev. Lett., vol. 64, p. 353, 1990.
  * (91) S. Weinberg, “Electromagnetic and weak masses,” Phys. Rev. Lett., vol. 29, pp. 388–392, 1972.
  * (92) R. N. Mohapatra, “Gauge Model for Chiral Symmetry Breaking and Muon electron Mass Ratio,” Phys. Rev. D, vol. 9, p. 3461, 1974.
  * (93) S. M. Barr and A. Zee, “A New Approach to the electron-Muon Mass Ratio,” Phys. Rev. D, vol. 15, p. 2652, 1977.
  * (94) B. S. Balakrishna, “Fermion Mass Hierarchy From Radiative Corrections,” Phys. Rev. Lett., vol. 60, p. 1602, 1988.
  * (95) S. M. Barr, “Flavor without flavor symmetry,” Phys. Rev. D, vol. 65, p. 096012, 2002.
  * (96) L. Ferretti, S. F. King, and A. Romanino, “Flavour from accidental symmetries,” JHEP, vol. 11, p. 078, 2006.
  * (97) S. M. Barr, “Doubly Lopsided Mass Matrices from Unitary Unification,” Phys. Rev. D, vol. 78, p. 075001, 2008.
  * (98) P. H. Frampton and T. W. Kephart, “Fermion Mixings in SU(9) Family Unification,” Phys. Lett. B, vol. 681, pp. 343–346, 2009.
  * (99) K. S. Babu and S. M. Barr, “Large neutrino mixing angles in unified theories,” Phys. Lett. B, vol. 381, pp. 202–208, 1996.
  * (100) J. Sato and T. Yanagida, “Large lepton mixing in a coset space family unification on E(7) / SU(5) x U(1)**3,” Phys. Lett. B, vol. 430, pp. 127–131, 1998.
  * (101) N. Irges, S. Lavignac, and P. Ramond, “Predictions from an anomalous U(1) model of Yukawa hierarchies,” Phys. Rev. D, vol. 58, p. 035003, 1998.
  * (102) S. M. Barr and I. Dorsner, “A General classification of three neutrino models and U(e3),” Nucl. Phys. B, vol. 585, pp. 79–104, 2000.
  * (103) N. Haba and H. Murayama, “Anarchy and hierarchy,” Phys. Rev. D, vol. 63, p. 053010, 2001.
  * (104) X.-G. He, Y.-Y. Keum, and R. R. Volkas, “A(4) flavor symmetry breaking scheme for understanding quark and neutrino mixing angles,” JHEP, vol. 04, p. 039, 2006.
  * (105) Y. H. Ahn, S. K. Kang, and C. S. Kim, “Spontaneous CP Violation in $A_{4}$ Flavor Symmetry and Leptogenesis,” Phys. Rev. D, vol. 87, no. 11, p. 113012, 2013.
  * (106) R. González Felipe, H. Serôdio, and J. P. Silva, “Models with three Higgs doublets in the triplet representations of $A_{4}$ or $S_{4}$,” Phys. Rev. D, vol. 87, no. 5, p. 055010, 2013.
  * (107) H. Ishimori and E. Ma, “New Simple $A_{4}$ Neutrino Model for Nonzero $\theta_{13}$ and Large $\delta_{CP}$,” Phys. Rev. D, vol. 86, p. 045030, 2012.
  * (108) F. Gonzalez Canales, A. Mondragon, and M. Mondragon, “The $S_{3}$ Flavour Symmetry: Neutrino Masses and Mixings,” Fortsch. Phys., vol. 61, pp. 546–570, 2013.
  * (109) R. N. Mohapatra and C. C. Nishi, “$S_{4}$ Flavored CP Symmetry for Neutrinos,” Phys. Rev. D, vol. 86, p. 073007, 2012.
  * (110) G.-J. Ding, S. F. King, C. Luhn, and A. J. Stuart, “Spontaneous CP violation from vacuum alignment in $S_{4}$ models of leptons,” JHEP, vol. 05, p. 084, 2013.
  * (111) E. Ma, “Neutrino Mixing and Geometric CP Violation with Delta(27) Symmetry,” Phys. Lett. B, vol. 723, pp. 161–163, 2013.
  * (112) C. C. Nishi, “Generalized $CP$ symmetries in $\Delta(27)$ flavor models,” Phys. Rev. D, vol. 88, no. 3, p. 033010, 2013.
  * (113) G. Altarelli and F. Feruglio, “Tri-bimaximal neutrino mixing from discrete symmetry in extra dimensions,” Nucl. Phys. B, vol. 720, pp. 64–88, 2005\.
  * (114) H. Ishimori, Y. Shimizu, M. Tanimoto, and A. Watanabe, “Neutrino masses and mixing from $S_{4}$ flavor twisting,” Phys. Rev. D, vol. 83, p. 033004, 2011.
  * (115) A. Kadosh and E. Pallante, “An A(4) flavor model for quarks and leptons in warped geometry,” JHEP, vol. 08, p. 115, 2010.
  * (116) G.-J. Ding and Y.-L. Zhou, “Dirac Neutrinos with $S_{4}$ Flavor Symmetry in Warped Extra Dimensions,” Nucl. Phys. B, vol. 876, pp. 418–452, 2013\.
  * (117) Y. Giraldo, “Texture Zeros and WB Transformations in the Quark Sector of the Standard Model,” Phys. Rev. D, vol. 86, p. 093021, 2012.
  * (118) G. C. Branco, D. Emmanuel-Costa, and R. Gonzalez Felipe, “Texture zeros and weak basis transformations,” Phys. Lett. B, vol. 477, pp. 147–155, 2000\.
  * (119) Y. Giraldo and E. Rojas, “CKM mixings from mass matrices with five texture zeros,” Phys. Rev. D, vol. 104, no. 7, p. 075009, 2021.
  * (120) P. Ramond, R. G. Roberts, and G. G. Ross, “Stitching the Yukawa quilt,” Nucl. Phys. B, vol. 406, pp. 19–42, 1993.
  * (121) H. Fritzsch and Z.-z. Xing, “Mass and flavor mixing schemes of quarks and leptons,” Prog. Part. Nucl. Phys., vol. 45, pp. 1–81, 2000.
  * (122) G. C. Branco, L. Lavoura, and F. Mota, “Nearest Neighbor Interactions and the Physical Content of Fritzsch Mass Matrices,” Phys. Rev. D, vol. 39, p. 3443, 1989.
  * (123) X.-G. He and W.-S. Hou, “Relating the Long $B$ Lifetime to a Very Heavy Top Quark,” Phys. Rev. D, vol. 41, p. 1517, 1990.
  * (124) H. Fritzsch, “Weak Interaction Mixing in the Six - Quark Theory,” Phys. Lett. B, vol. 73, pp. 317–322, 1978.
  * (125) H. Fritzsch, “Calculating the Cabibbo Angle,” Phys. Lett. B, vol. 70, pp. 436–440, 1977.
  * (126) H. Fritzsch, “Quark Masses and Flavor Mixing,” Nucl. Phys. B, vol. 155, pp. 189–207, 1979.
  * (127) D.-s. Du and Z.-z. Xing, “A Modified Fritzsch ansatz with additional first order perturbation,” Phys. Rev. D, vol. 48, pp. 2349–2352, 1993.
  * (128) H. Fritzsch and Z.-z. Xing, “Four zero texture of Hermitian quark mass matrices and current experimental tests,” Phys. Lett. B, vol. 555, pp. 63–70, 2003.
  * (129) R. L. Workman et al., “Review of Particle Physics,” PTEP, vol. 2022, p. 083C01, 2022.
  * (130) Z.-z. Xing and H. Zhang, “Complete parameter space of quark mass matrices with four texture zeros,” J. Phys. G, vol. 30, pp. 129–136, 2004.
  * (131) Y.-F. Zhou, “Textures and hierarchies in quark mass matrices with four texture zeros,” 9 2003.
  * (132) R. H. Benavides, Y. Giraldo, L. Muñoz, W. A. Ponce, and E. Rojas, “Five Texture Zeros for Dirac Neutrino Mass Matrices,” J. Phys. G, vol. 47, no. 11, p. 115002, 2020.
  * (133) Y. Giraldo, “Seeking Texture Zeros in the Quark Mass Matrix Sector of the Standard Model,” Nucl. Part. Phys. Proc., vol. 267-269, pp. 76–78, 2015\.
  * (134) R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett., vol. 38, pp. 1440–1443, 1977.
  * (135) R. D. Peccei and H. R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D, vol. 16, pp. 1791–1797, 1977.
  * (136) S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett., vol. 40, pp. 223–226, 1978.
  * (137) F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett., vol. 40, pp. 279–282, 1978.
  * (138) J. E. Kim, “Weak Interaction Singlet and Strong CP Invariance,” Phys. Rev. Lett., vol. 43, p. 103, 1979.
  * (139) M. S. Turner, “Windows on the Axion,” Phys. Rept., vol. 197, pp. 67–97, 1990.
  * (140) G. G. Raffelt, “Astrophysical methods to constrain axions and other novel particle phenomena,” Phys. Rept., vol. 198, pp. 1–113, 1990.
  * (141) J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B, vol. 120, pp. 127–132, 1983.
  * (142) L. F. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B, vol. 120, pp. 133–136, 1983.
  * (143) M. Dine and W. Fischler, “The Not So Harmless Axion,” Phys. Lett. B, vol. 120, pp. 137–141, 1983.
  * (144) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, “Can Confinement Ensure Natural CP Invariance of Strong Interactions?,” Nucl. Phys. B, vol. 166, pp. 493–506, 1980.
  * (145) M. Dine, W. Fischler, and M. Srednicki, “A Simple Solution to the Strong CP Problem with a Harmless Axion,” Phys. Lett. B, vol. 104, pp. 199–202, 1981.
  * (146) A. R. Zhitnitsky, “On Possible Suppression of the Axion Hadron Interactions. (In Russian),” Sov. J. Nucl. Phys., vol. 31, p. 260, 1980.
  * (147) P. Sikivie and Q. Yang, “Bose-Einstein Condensation of Dark Matter Axions,” Phys. Rev. Lett., vol. 103, p. 111301, 2009.
  * (148) P. Sikivie, “Axion Cosmology,” Lect. Notes Phys., vol. 741, pp. 19–50, 2008.
  * (149) S. Hannestad and G. Raffelt, “Cosmological mass limits on neutrinos, axions, and other light particles,” JCAP, vol. 04, p. 008, 2004.
  * (150) S. Hannestad, A. Mirizzi, and G. Raffelt, “New cosmological mass limit on thermal relic axions,” JCAP, vol. 07, p. 002, 2005.
  * (151) A. Melchiorri, O. Mena, and A. Slosar, “An improved cosmological bound on the thermal axion mass,” Phys. Rev. D, vol. 76, p. 041303, 2007.
  * (152) S. Hannestad, A. Mirizzi, G. G. Raffelt, and Y. Y. Y. Wong, “Cosmological constraints on neutrino plus axion hot dark matter,” JCAP, vol. 08, p. 015, 2007.
  * (153) S. Hannestad, A. Mirizzi, G. G. Raffelt, and Y. Y. Y. Wong, “Cosmological constraints on neutrino plus axion hot dark matter: Update after WMAP-5,” JCAP, vol. 04, p. 019, 2008.
  * (154) N. Viaux, M. Catelan, P. B. Stetson, G. Raffelt, J. Redondo, A. A. R. Valcarce, and A. Weiss, “Neutrino and axion bounds from the globular cluster M5 (NGC 5904),” Phys. Rev. Lett., vol. 111, p. 231301, 2013.
  * (155) E. Arik et al., “Probing eV-scale axions with CAST,” JCAP, vol. 02, p. 008, 2009.
  * (156) Y. Inoue, Y. Akimoto, R. Ohta, T. Mizumoto, A. Yamamoto, and M. Minowa, “Search for solar axions with mass around 1 eV using coherent conversion of axions into photons,” Phys. Lett. B, vol. 668, pp. 93–97, 2008.
  * (157) L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, “The landscape of QCD axion models,” Phys. Rept., vol. 870, pp. 1–117, 2020.
  * (158) M. M. Miller Bertolami, B. E. Melendez, L. G. Althaus, and J. Isern, “Revisiting the axion bounds from the Galactic white dwarf luminosity function,” JCAP, vol. 10, p. 069, 2014.
  * (159) M. Giannotti, I. G. Irastorza, J. Redondo, A. Ringwald, and K. Saikawa, “Stellar Recipes for Axion Hunters,” JCAP, vol. 10, p. 010, 2017.
  * (160) J. Isern, E. García-Berro, S. Torres, R. Cojocaru, and S. Catalán, “Axions and the luminosity function of white dwarfs: the thin and thick discs, and the halo,” Monthly Notices of the Royal Astronomical Society, vol. 478, pp. 2569–2575, 05 2018.
  * (161) O. Straniero, I. Dominguez, M. Giannotti, and A. Mirizzi, “Axion-electron coupling from the RGB tip of Globular Clusters,” arXiv e-prints, p. arXiv:1802.10357, Feb. 2018.
  * (162) R. Verma, “Exploring the predictability of symmetric texture zeros in quark mass matrices,” Phys. Rev. D, vol. 96, no. 9, p. 093010, 2017.
  * (163) Z.-z. Xing, “Flavor structures of charged fermions and massive neutrinos,” Phys. Rept., vol. 854, pp. 1–147, 2020.
  * (164) B. R. Desai and A. R. Vaucher, “Quark mass matrices with four and five texture zeroes, and the CKM matrix, in terms of mass eigenvalues,” Phys. Rev. D, vol. 63, p. 113001, 2001.
  * (165) W. A. Ponce, J. D. Gómez, and R. H. Benavides, “Five texture zeros and CP violation for the standard model quark mass matrices,” Phys. Rev. D, vol. 87, no. 5, p. 053016, 2013.
  * (166) Y. Giraldo and E. Rojas, “Five Non-Fritzsch Texture Zeros for Quarks Mass Matrices in the Standard Model,” in 38th International Symposium on Physics in Collision, 11 2018.
  * (167) C. Hagedorn and W. Rodejohann, “Minimal mass matrices for dirac neutrinos,” JHEP, vol. 07, p. 034, 2005.
  * (168) G. Ahuja, M. Gupta, M. Randhawa, and R. Verma, “Texture specific mass matrices with Dirac neutrinos and their implications,” Phys. Rev. D, vol. 79, p. 093006, 2009.
  * (169) X.-w. Liu and S. Zhou, “Texture Zeros for Dirac Neutrinos and Current Experimental Tests,” Int. J. Mod. Phys. A, vol. 28, p. 1350040, 2013.
  * (170) R. Verma, “Lower bound on neutrino mass and possible CP violation in neutrino oscillations,” Phys. Rev. D, vol. 88, p. 111301, 2013.
  * (171) R. Verma, “Lepton textures and neutrino oscillations,” Int. J. Mod. Phys. A, vol. 29, no. 21, p. 1444009, 2014.
  * (172) P. Fakay, S. Sharma, G. Ahuja, and M. Gupta, “Leptonic mixing angle $\theta_{13}$ and ruling out of minimal texture for Dirac neutrinos,” PTEP, vol. 2014, no. 2, p. 023B03, 2014.
  * (173) L. M. Cebola, D. Emmanuel-Costa, and R. G. Felipe, “Confronting predictive texture zeros in lepton mass matrices with current data,” Phys. Rev. D, vol. 92, no. 2, p. 025005, 2015.
  * (174) R. R. Gautam, M. Singh, and M. Gupta, “Neutrino mass matrices with one texture zero and a vanishing neutrino mass,” Phys. Rev. D, vol. 92, no. 1, p. 013006, 2015.
  * (175) G. Ahuja, S. Sharma, P. Fakay, and M. Gupta, “General lepton textures and their implications,” Mod. Phys. Lett. A, vol. 30, no. 34, p. 1530025, 2015\.
  * (176) M. Singh, “Texture One Zero Dirac Neutrino Mass Matrix With Vanishing Determinant or Trace Condition,” Nucl. Phys. B, vol. 931, pp. 446–468, 2018.
  * (177) G. Ahuja and M. Gupta, “Texture zero mass matrices and nature of neutrinos,” Int. J. Mod. Phys. A, vol. 33, no. 31, p. 1844032, 2018.
  * (178) J. Barranco, D. Delepine, and L. Lopez-Lozano, “Neutrino Mass Determination from a Four-Zero Texture Mass Matrix,” Phys. Rev. D, vol. 86, p. 053012, 2012.
  * (179) Y. A. Garnica, S. F. Mantilla, R. Martinez, and H. Vargas, “From Peccei Quinn symmetry to mass hierarchy problem,” J. Phys. G, vol. 48, no. 9, p. 095002, 2021.
  * (180) A. Ringwald and K. Saikawa, “Axion dark matter in the post-inflationary Peccei-Quinn symmetry breaking scenario,” Phys. Rev. D, vol. 93, no. 8, p. 085031, 2016.  [Addendum: Phys.Rev.D 94, 049908 (2016)].
  * (181) I. Brivio, M. B. Gavela, L. Merlo, K. Mimasu, J. M. No, R. del Rey, and V. Sanz, “ALPs Effective Field Theory and Collider Signatures,” Eur. Phys. J. C, vol. 77, no. 8, p. 572, 2017.
  * (182) H. Georgi, D. B. Kaplan, and L. Randall, “Manifesting the Invisible Axion at Low-energies,” Phys. Lett. B, vol. 169, pp. 73–78, 1986.
  * (183) M. B. Gavela, R. Houtz, P. Quilez, R. Del Rey, and O. Sumensari, “Flavor constraints on electroweak ALP couplings,” Eur. Phys. J. C, vol. 79, no. 5, p. 369, 2019.
  * (184) M. Bauer, M. Neubert, and A. Thamm, “Collider Probes of Axion-Like Particles,” JHEP, vol. 12, p. 044, 2017.
  * (185) A. Salvio, A. Strumia, and W. Xue, “Thermal axion production,” JCAP, vol. 01, p. 011, 2014.
  * (186) E. Cortina Gil et al., “Measurement of the very rare K${}^{+}\rightarrow{\pi}^{+}\nu\overline{\nu}$ decay,” JHEP, vol. 06, p. 093, 2021.
  * (187) F. Björkeroth, E. J. Chun, and S. F. King, “Flavourful Axion Phenomenology,” JHEP, vol. 08, p. 117, 2018.
  * (188) R. D. Bolton et al., “Search for Rare Muon Decays with the Crystal Box Detector,” Phys. Rev. D, vol. 38, p. 2077, 1988.
  * (189) A. Jodidio, B. Balke, J. Carr, G. Gidal, K. A. Shinsky, H. M. Steiner, D. P. Stoker, M. Strovink, R. D. Tripp, B. Gobbi, and C. J. Oram, “Erratum: Search for right-handed currents in muon decay [phys. rev. d 34, 1967 (1986)],” Phys. Rev. D, vol. 37, pp. 237–238, Jan 1988.
  * (190) R. Bayes et al., “Search for two body muon decay signals,” Phys. Rev. D, vol. 91, no. 5, p. 052020, 2015.
  * (191) H. Albrecht et al., “A Search for lepton flavor violating decays tau —-$>$ e alpha, tau —$>$ mu alpha,” Z. Phys. C, vol. 68, pp. 25–28, 1995.
  * (192) H. Georgi and D. V. Nanopoulos, “Suppression of Flavor Changing Effects From Neutral Spinless Meson Exchange in Gauge Theories,” Phys. Lett. B, vol. 82, pp. 95–96, 1979.
  * (193) S. Chaudhuri, P. W. Graham, K. Irwin, J. Mardon, S. Rajendran, and Y. Zhao, “Radio for hidden-photon dark matter detection,” Phys. Rev. D, vol. 92, no. 7, p. 075012, 2015.
  * (194) D. Alesini, D. Babusci, D. Di Gioacchino, C. Gatti, G. Lamanna, and C. Ligi, “The KLASH Proposal,” 7 2017.
  * (195) C. Gatti et al., “The Klash Proposal: Status and Perspectives,” in 14th Patras Workshop on Axions, WIMPs and WISPs, 11 2018.
  * (196) I. Stern, “ADMX Status,” PoS, vol. ICHEP2016, p. 198, 2016.
  * (197) M. Lawson, A. J. Millar, M. Pancaldi, E. Vitagliano, and F. Wilczek, “Tunable axion plasma haloscopes,” Phys. Rev. Lett., vol. 123, no. 14, p. 141802, 2019.
  * (198) S. Beurthey et al., “MADMAX Status Report,” 3 2020.
  * (199) E. Armengaud et al., “Conceptual Design of the International Axion Observatory (IAXO),” JINST, vol. 9, p. T05002, 2014.
  * (200) I. Shilon, A. Dudarev, H. Silva, U. Wagner, and H. H. J. ten Kate, “New Superconducting Toroidal Magnet System for IAXO, the International AXion Observatory,” AIP Conf. Proc., vol. 1573, no. 1, pp. 1559–1566, 2015\.
  * (201) J. L. Ouellet et al., “First Results from ABRACADABRA-10 cm: A Search for Sub-$\mu$eV Axion Dark Matter,” Phys. Rev. Lett., vol. 122, no. 12, p. 121802, 2019.
  * (202) C. Bartram et al., “Search for Invisible Axion Dark Matter in the 3.3–4.2 $\mu$eV Mass Range,” Phys. Rev. Lett., vol. 127, no. 26, p. 261803, 2021.
  * (203) T. Braine et al., “Extended Search for the Invisible Axion with the Axion Dark Matter Experiment,” Phys. Rev. Lett., vol. 124, no. 10, p. 101303, 2020.
  * (204) C. Bartram et al., “Axion dark matter experiment: Run 1B analysis details,” Phys. Rev. D, vol. 103, no. 3, p. 032002, 2021.
  * (205) V. Anastassopoulos et al., “New CAST Limit on the Axion-Photon Interaction,” Nature Phys., vol. 13, pp. 584–590, 2017.
  * (206) S. Andriamonje et al., “An Improved limit on the axion-photon coupling from the CAST experiment,” JCAP, vol. 04, p. 010, 2007.
  * (207) S. Lee, S. Ahn, J. Choi, B. R. Ko, and Y. K. Semertzidis, “Axion Dark Matter Search around 6.7 $\mu$eV,” Phys. Rev. Lett., vol. 124, no. 10, p. 101802, 2020.
  * (208) J. Jeong, S. Youn, S. Bae, J. Kim, T. Seong, J. E. Kim, and Y. K. Semertzidis, “Search for Invisible Axion Dark Matter with a Multiple-Cell Haloscope,” Phys. Rev. Lett., vol. 125, no. 22, p. 221302, 2020.
  * (209) O. Kwon et al., “First Results from an Axion Haloscope at CAPP around 10.7 $\mu$eV,” Phys. Rev. Lett., vol. 126, no. 19, p. 191802, 2021.
  * (210) K. M. Backes et al., “A quantum-enhanced search for dark matter axions,” Nature, vol. 590, no. 7845, pp. 238–242, 2021.
  * (211) L. Zhong et al., “Results from phase 1 of the HAYSTAC microwave cavity axion experiment,” Phys. Rev. D, vol. 97, no. 9, p. 092001, 2018.
  * (212) P. Gondolo and G. G. Raffelt, “Solar neutrino limit on axions and keV-mass bosons,” Phys. Rev. D, vol. 79, p. 107301, 2009.
  * (213) A. Ayala, I. Domínguez, M. Giannotti, A. Mirizzi, and O. Straniero, “Revisiting the bound on axion-photon coupling from Globular Clusters,” Phys. Rev. Lett., vol. 113, no. 19, p. 191302, 2014.
  * (214) M. Regis, M. Taoso, D. Vaz, J. Brinchmann, S. L. Zoutendijk, N. F. Bouché, and M. Steinmetz, “Searching for light in the darkness: Bounds on ALP dark matter with the optical MUSE-faint survey,” Phys. Lett. B, vol. 814, p. 136075, 2021.
  * (215) D. Grin, G. Covone, J.-P. Kneib, M. Kamionkowski, A. Blain, and E. Jullo, “A Telescope Search for Decaying Relic Axions,” Phys. Rev. D, vol. 75, p. 105018, 2007.
  * (216) V. Prasolov, Problems and theorems in linear algebra.  USA: American Mathematical Society, 1944.
  * (217) S. Hassani, Mathematical Physics: A Modern Introduction to Its Foundations.  USA: Springer International Publishing, 2013.
  * (218) Z.-z. Xing and Z.-h. Zhao, “On the four-zero texture of quark mass matrices and its stability,” Nucl. Phys. B, vol. 897, pp. 302–325, 2015.
  * (219) M. Gupta and G. Ahuja, “Flavor mixings and textures of the fermion mass matrices,” Int. J. Mod. Phys. A, vol. 27, p. 1230033, 2012.
  * (220) H. Fusaoka and Y. Koide, “Updated estimate of running quark masses,” Phys. Rev. D, vol. 57, pp. 3986–4001, 1998.
  * (221) J. Cardozo, J. H. Muñoz, N. Quintero, and E. Rojas, “Analysing the charged scalar boson contribution to the charged-current $B$ meson anomalies,” J. Phys. G, vol. 48, no. 3, p. 035001, 2021.
  * (222) I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni, and T. Schwetz, “Global analysis of three-flavour neutrino oscillations: synergies and tensions in the determination of $\theta_{23}$, $\delta_{CP}$, and the mass ordering,” JHEP, vol. 01, p. 106, 2019.
  * (223) G. Alonso-Álvarez, M. B. Gavela, and P. Quilez, “Axion couplings to electroweak gauge bosons,” Eur. Phys. J. C, vol. 79, no. 3, p. 223, 2019\.
  * (224) G. Grilli di Cortona, E. Hardy, J. Pardo Vega, and G. Villadoro, “The QCD axion, precisely,” JHEP, vol. 01, p. 034, 2016.