Simple Modular invariant model
for Quark, Lepton, and flavored-QCD axion

Y. H. Ahn¹¹¹1Email: [email protected], Sin Kyu Kang^1,2²²2Email: [email protected] ¹Institute for Convergence Fundamental Study, Seoul National University of Science and Technology,
Seoul 139-743, Korea ²School of Liberal Arts, Seoul National University of Science and Technology, Seoul 139-743, Korea

Abstract

We propose a minimal extension of the Standard Model by incorporating sterile neutrinos and a QCD axion to account for the mass and mixing hierarchies of quarks and leptons and to solve the strong CP problem, and by introducing $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ symmetry. We demonstrate that the Kähler transformation corrects the weight of modular forms in the superpotential and show that the model is consistent with the modular and $U(1)_{X}$ anomaly-free conditions. This enables a simple construction of a modular-independent superpotential for scalar potential. Using minimal supermultiplets, we demonstrate a level 3 modular form-induced superpotential. Sterile neutrinos explain small active neutrino masses via the seesaw mechanism and provide a well-motivated $U(1)_{X}$ breaking scale, whereas gauge singlet scalar fields play crucial roles in generating the QCD axion, heavy neutrino mass, and fermion mass hierarchy. The model predicts a range for the $U(1)_{X}$ breaking scale from $10^{13}$ GeV to $10^{15}$ GeV for $1\,\mbox{TeV}<m_{3/2}<10^{6}\,\mbox{TeV}$ . In the supersymmetric limit, all Yukawa coefficients in the superpotential are given by complex numbers with an absolute value of unity, implying a democratic distribution. Performing numerical analysis, we study how model parameters are constrained by current experimental results. In particular, the model predicts that the value of the quark Dirac CP phase falls between $38^{\circ}$ to $87^{\circ}$ , which is consistent with experimental data, and the favored value of the neutrino Dirac CP phase is around $250^{\circ}$ . Furthermore, the model can be tested by ongoing and future experiments on axion searches, neutrino oscillations, and $0\nu\beta\beta$ -decay.

I Introduction

Despite being theoretically self-consistent and successfully demonstrating experimental results in low-energy experiments so far, the Standard Model (SM) of particle physics leaves unanswered questions in theoretical and cosmological issues and fails to explain some physical phenomena such as neutrino oscillations, muon $g-2$ , etc. Various attempts have been made to extend the SM in order to address these questions and account for experimental results that cannot be explained within the SM. For instance, the canonical seesaw mechanism Minkowski:1977sc has been proposed to explain the tiny masses of neutrinos by introducing new heavy neutral fermions alongside the SM particles. Additionally, the Peccei-Quinn (PQ) mechanism Peccei-Quinn has been suggested to solve the strong CP problem in quantum chromodynamics (QCD) by extending the SM to include an anomalous $U(1)_{X}$ symmetry.

Recently, Feruglio Feruglio:2017spp proposed a new idea regarding the origin of the structure of lepton mixing. He applied modular invariance ³³3Modular invariance was analyzed for supersymmetric field theories in Ref.Ferrara:1989bc . under the modular group to determine the flavor structure of leptons without introducing a number of scalar fields. This approach requires the Yukawa couplings among twisted states to be modular forms. It is a string-derived mechanism that naturally restricts the possible variations in the flavor structure of quarks and leptons, which are unconstrained by the SM gauge invariance. However, explaining the hierarchies of the masses and mixing in the quark and lepton sectors remains a challenge. As studied in most references modular , the Yukawa coefficients are assumed to be free parameters ⁴⁴4In Ref.modular2 , the Yukawa coefficients are assumed to be of order unity. which can be determined by matching them with experimental data on fermion mass and mixing hierarchies. This approach is not significantly different from that in the SM, except for the introduction of modular forms. Alternatively, it is also possible to take the Yukawa coefficient to be of order unity, accommodating the hierarchies of the fermion masses and mixing. Recently, Ref.Feruglio:2023uof has demonstrated that the vanishing QCD angle, a large CKM phase, and the reproduction of quark and lepton masses and mixings can be achieved by using coefficients up to order one, see also Ref.Kobayashi:2020oji .

To incorporate sterile neutrinos and a QCD axion into the SM and provide a natural explanation for the mass and mixing hierarchies of quarks and leptons, we propose an extension of a modular invariant model based on the four dimensional (4D) effective action derived from superstring theory with $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ symmetry. The non-Abelian discrete symmetry $\Gamma_{N}$ with $N=2,3,4,5$ plays a role of modular invariance, and may originate from superstring theory in compactified extra dimensions, where it acts as a finite subgroup of the modular group deAdelhartToorop:2011re . To ensure the validity of a modular invariant model with $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ , we take the followings into account:

(i)

T-duality relates one type of superstring theory to another, and it also appears in the 4D low-energy effective field theory derived from superstring theory (for a review, see Ref.string_book ). In particular, 4D low-energy effective field theory of type IIA string theory with a certain compactification is invariant under the modular transformation of the modulus $\tau$ ,

\displaystyle\tau\rightarrow\gamma\tau=\frac{a\tau+b}{c\tau+d}\,,\quad(a,b,c,d\in{Z},~{}~{}ad-bc=1)\,.

(1)

So, the 4D action we conisder is required to be invariant under the modular transformation and gauged $U(1)_{X}$ symmetry, as well as the Kähler transformation (refer to Eq.(9)). This is necessary to cancel out the modular anomaly (see Ref.Derendinger:1991hq ) associated with the modular transformation Eq.(1) under the non-local modular group $\Gamma_{N}$ and the gauged $U(1)_{X}$ anomaly, at the quantum level.

(ii)

While type II string theory allows for low axion decay constant models via D-branes, leading to the gauged $U(1)_{X}$ that becomes a global PQ symmetry when the $U(1)_{X}$ gauge boson is decoupled Dine:1987xk , heterotic string theory typically has a $U(1)_{X}$ breaking scale with a decay constant close to the string scale. The broken $U(1)_{X}$ gauge symmetry leaves behind a protected global $U(1)_{X}$ that is immune to quantum-gravitational effects, achieved via the Green-Schwarz (GS) mechanism Green:1984sg . The PQ breaking scale, or the low axion decay constant, can be determined by taking into account both supersymmetry (SUSY)-breaking effects Burgess:2003ic and supersymmetric next-leading-order Planck-suppressed terms Nanopoulos:1983sp ; Ahn:2016hbn ; Ahn:2017dpf .

The model features a minimal set of fields that transform based on representations of $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ , and includes modular forms of level $N$ . These modular forms act as Yukawa couplings and transform under the modular group $\Gamma_{N}$ . It should be noted that the Kähler transformation (refer to Eq.(9)) corrects the weight of modular forms in the superpotential due to the modular invariance of both the superpotential and Kähler potential, see Eq.(20). This enables a simple construction of a $\tau$ -independent superpotential for scalar potential. The so-called flavored-PQ symmetry $U(1)_{X}$ guarantees the absence of bare mass terms Ahn:2014gva . We minimally extend the model by incorporating three right-handed neutrinos $N^{c}$ and SM gauge singlet scalar fields $\chi(\tilde{\chi})$ . The scalar fields with a modular weight of zero and charged by $+(-)$ under $U(1)_{X}$ play a crucial role in generating the QCD axion, heavy neutrino mass, and fermion mass hierarchy. Then the complex scalar field ${\cal F}=\chi(\tilde{\chi})$ with modular weight zero acts on dimension-four(three) operators well-sewed by $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ and modular invariance with different orders, which generate the effective interactions for the SM and the right-handed neutrinos as follows,

\displaystyle\tilde{c}_{1}\,{\cal O}_{3}\,({\cal F})^{1}+{\cal O}_{4}\sum_{n=0}^{\rm finite}c_{n}\,\left(\frac{{\cal F}}{\Lambda}\right)^{n}+...

(2)

Here, $\Lambda$ is the scale of flavor dynamics above which unknown physics exists as a UV cutoff, and Yukawa coefficients $c_{n}(\tilde{c}_{1})$ are all complex numbers assumed to have a unit absolute value ( $|\tilde{c}_{1}|,|c_{n}|=1$ ). The dimension- $4(3)$ operators ${\cal O}_{4(3)}$ are determined by $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ and modular invariance in the supersymmetric limit. These operators include modular forms of level $N$ , which transform according to the representation of $\Gamma_{N}$ Feruglio:2017spp . We will demonstrate that any additive finite correction terms, which could potentially be generated by higher weight modular forms, are prohibited due to the modular weight of the $\chi(\tilde{\chi})$ fields being zero. Note that there exist the infinite series of higher-dimensional operators induced solely by the combination of $\chi\tilde{\chi}$ in the supersymmetric limit. These operators, represented by dots in Eq.(2), can be absorbed into the finite leading order terms and effectively modify the coefficients $\tilde{c}_{1}$ and $c_{n}$ at the leading order. Furthermore, to avoid the breaking effects of the axionic shift symmetry caused by gravity that spoil the axion solution to the strong CP problem Kamionkowski:1992mf , we imposed a $U(1)_{X}$ -mixed gravitational anomaly-free condition Ahn:2016hbn ; Ahn:2018cau ; Ahn:2021ndu .

The rest of this paper is organized as follows: The next section discusses modular and $U(1)_{X}$ anomaly-free conditions under $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ symmetry, along with the modular forms of superpotential corrected by Kähler transformation. Sec.III presents an example of a superpotential induced by level 3 modular forms. We introduce minimal supermultiplets to address the challenges of tiny neutrino masses, the strong CP problem, and the hierarchies of SM fermion mass and mixing. For our purpose, we show how to derive Yukawa superpotentials and a modular-independent superpotential for the scalar potential and determining the relevant $U(1)_{X}$ PQ symmetry breaking scale (or seesaw scale). Additionally, we provide comments on the modular invariant model. In Sec.IV, we visually demonstrate the interconnections between quarks, leptons, and flavored-QCD axion. In Sec.V, we present numerical values of physical parameters that satisfy the current experimental data on flavor mixing and mass for quarks and leptons while also favoring the assumption in Eq.(2). The study predicts the Dirac CP phases of quarks and leptons, as well as the mass of the flavored-QCD axion and its coupling to photons and electrons. The final section provides a summary of our work.

II Modular and U(1) anomaly-free

T-duality relates different types of superstring theory and is also present in the 4D low-energy effective field theory derived from superstring theory (see Ref.string_book for a review). In particular, type IIA intersecting D-brane models are related to magnetized D-brane models through T-duality string_book . The group $\Gamma(N)$ acts on the complex variable $\tau$ , varying in the upper-half complex plane ${\rm Im}(\tau)>0$ , as the modular transformation Eq.(1). Then the low-energy effective field theory of type IIA intersecting D-brane models must have the symmetry under the modular transformation Eq.(1). First, we shortly review the modular symmetry. The infinite groups $\Gamma(N)$ , called principal congruence subgroups of level $N=1,2,3,...$ , are defined by

\displaystyle\Gamma(N)=\Big{\{}\begin{pmatrix}a&b\\ c&d\end{pmatrix}\in SL(2,Z),\begin{pmatrix}a&b\\ c&d\end{pmatrix}=\begin{pmatrix}1&0\\ 0&1\end{pmatrix}\quad(\text{mod}~{}N)\Big{\}}\,,

(3)

which are normal subgroups of homogeneous modular group $\Gamma\equiv\Gamma(1)\simeq SL(2,Z)$ , where $SL(2,Z)$ is the group of $2\times 2$ matrices with integer entries and determinant equal to one. The projective principal congruence subgroups are defined as $\bar{\Gamma}(N)=\Gamma(N)/\{\pm{I}\}$ for $N=1,2$ . For $N\geq 3$ we have $\bar{\Gamma}(N)=\Gamma(N)$ because the elements $-I$ does not belong to $\Gamma(N)$ . The modular group $\bar{\Gamma}\equiv\Gamma/\{\pm{I}\}$ is generated by two elements $S$ and $T$ ,

\displaystyle S:\tau\rightarrow-\frac{1}{\tau}\,,\qquad\qquad T:\tau\rightarrow\tau+1\,,

(4)

satisfying

\displaystyle S^{2}=(ST)^{3}=(TS)^{3}={\bf 1}\,.

(5)

They can be represented by the $PSL(2,Z)$ matrices

\displaystyle S=\begin{pmatrix}0&1\\ -1&0\end{pmatrix}\,,\qquad T=\begin{pmatrix}1&1\\ 0&1\end{pmatrix}\,.

(6)

The groups $\Gamma_{N}$ are finite modular groups obtained by imposing the condition $T^{N}={\bf 1}$ in addition to Eq.(5), where $\Gamma_{N}\equiv\bar{\Gamma}/\bar{\Gamma}(N)$ . The groups $\Gamma_{N}$ are isomorphic to the permutation groups $S_{3}$ , $A_{4}$ , $S_{4}$ , and $A_{5}$ for $N=2,3,4,5$ , respectively deAdelhartToorop:2011re .

We work in the 4D ${\cal N}=1$ string-derived supergravity framework defined by a general Kähler potential $G(\Phi,\bar{\Phi})$ of the chiral superfields $\Phi$ and their conjugates,

\displaystyle G(\Phi,\bar{\Phi})=\frac{K(\Phi,\bar{\Phi})}{M^{2}_{P}}+\ln\frac{|W(\Phi)|^{2}}{M^{6}_{P}}\,,

(7)

and by an analytic gauge kinetic function $f(\Phi)$ of the chiral superfields $\Phi$ , where $M_{P}=(8\pi G_{N})^{-1/2}=2.436\times 10^{18}$ GeV is the reduced Planck mass with Newton’s gravitational constant $G_{N}$ , $K(\Phi,\bar{\Phi})$ is a real gauge-invariant function of $\Phi$ and $\bar{\Phi}$ , and $W(\Phi)$ is a holomorphic gauge-invariant function of $\Phi$ . Based on the 4D effective field theory derived from type IIA intersecting D-brane models, we build a modular invariant model with minimal chiral superfields transforming according to representations of $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ . Here we assume that the non-Abelian discrete symmetry $\Gamma_{N}$ as a finite subgroup of the modular group Feruglio:2017spp , and the anomalous gauged $U(1)_{X}$ including the SM gauge symmetry $G_{\rm SM}$ may arise from several stacks on D-brane models string_book . In the 4D global supersymmetry the most general form of the action can be written as

\displaystyle{\cal S}=\int d^{4}xd^{2}\theta d^{2}\bar{\theta}K(\Phi,\bar{\Phi}e^{2A})+\Big{\{}\int d^{4}xd^{2}\theta\Big{(}W(\Phi)+\frac{f_{ab}(\Phi)}{4}{\cal W}^{\alpha a}{\cal W}^{b}_{\alpha}\Big{)}+{\rm h.c.}\Big{\}}\,,

(8)

where $A\equiv A^{a}T^{a}$ is the gauge multiplet containing Yang-Mills multiplet and $T^{a}$ are the gauge group generators, and ${\cal W}_{\alpha}$ is a gauge-invariant chiral spinor superfield containing the Yang-Mills field strength. The chiral superfields $\Phi$ denote all chiral supermultiplets with Kähler moduli, complex structure moduli, axio-dilaton, and matter superfields, transforming under $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ . We assume that the low-energy Kähler potential $K$ , superpotential $W$ , and gauge kinetic function $f$ for moduli and matter superfields are given at a scale where Kähler moduli and complex structure moduli are stabilized through fluxes (see Ref.Gukov:1999ya ; Giddings:2001yu ; Dasgupta:1999ss ) leading to a consistent low-energy SM gauge theory.

Under the modular transformation Eq.(1) and the gauged $U(1)_{X}$ symmetry, the action (8) should be invariant with the transformations ⁵⁵5The upper two shifts in Eq.(9) of the Kähler potential and superpotential are known as the ‘Kähler transformation’ with reference to Eq.(7).

	$\displaystyle K(\Phi,\bar{\Phi}e^{2A})\rightarrow K(\Phi,\bar{\Phi}e^{2A})+\big{(}g(\Phi)+g(\bar{\Phi})\big{)}M^{2}_{P}\,,$
	$\displaystyle W(\Phi)\rightarrow W(\Phi)e^{-g(\Phi)}\,,$
	$\displaystyle f(\Phi){\cal W}^{\alpha}{\cal W}_{\alpha}\rightarrow f(\Phi){\cal W}^{\alpha}{\cal W}_{\alpha}\,,$		(9)

where $g(\Phi)$ is a function of modulus $\tau$ . Then the given symmetry $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ can be violated at the quantum level by (i) an anomalous triangle graph associated with modular transformation Eq.(1) under the non-local modular group $\Gamma_{N}$ , and (ii) anomalous triangle graphs with external states $A^{a}_{\nu}A^{b}_{\rho}V_{X\mu}$ where $A^{a}_{\nu}$ , $A^{b}_{\rho}$ are gauge bosons of the SM gauge group $G_{\rm SM}$ and $V^{\mu}_{X}$ is the connection associated with the gauged $U(1)_{X}$ . These anomalies can be cancelled by the GS mechanism Green:1984sg .

II.1 Modular anomaly-free and modular forms of level $N$

To demonstrate the invariance of $K(\Phi,\bar{\Phi}e^{2A})$ and $f(\Phi){\cal W}^{\alpha}{\cal W}_{\alpha}$ of Eq.(8) under the finite modular group $\Gamma_{N}$ and the gauged $U(1)_{X}$ , we consider a low-energy Kähler potential ⁶⁶6It is similar to the one-loop Kähler potential presented in Ref.Derendinger:1991hq .:

	$\displaystyle K$	$\displaystyle=$	$\displaystyle-M^{2}_{P}\ln\Big{\{}\Big{(}S+\bar{S}-3\tilde{c}\ln(-i\tau+i\bar{\tau})\Big{)}\Big{(}U_{X}+\bar{U}_{X}-\frac{\delta^{\rm GS}_{X}}{16\pi^{2}}V_{X}\Big{)}\prod^{2}_{i=1}({\cal U}_{i}+\bar{\cal U}_{i})\Big{\}}$		(10)
		$\displaystyle-$	$\displaystyle M^{2}_{P}\ln\big{(}-i\tau+i\bar{\tau}\big{)}^{3}+(-i\tau+i\bar{\tau})^{-k}\|\varphi\|^{2}+Z_{X}\varphi^{\dagger}_{X}e^{-XV_{X}}\varphi_{X}+...\,,$		(10)

where $-k$ is the modular weight, $Z_{X}$ is the normalization factor, $S$ denotes the axio-dilaton field, $\tau$ represents the overall Kähler modulus, and $U_{X}$ and ${\cal U}_{i}$ correspond to the complex structure moduli. The dots in Eq.(10) denote the contributions of non-renormalizable terms scaled by an UV cutoff $M_{P}$ invariant under $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ . We note that the matter fields $\varphi_{X}$ with $U(1)_{X}$ charge, complex structure modulus $U_{X}$ and the vector superfield $V_{X}$ of the gauged $U(1)_{X}$ including the gauge field $A^{\mu}_{X}$ participate in the 4D GS mechanism. We take the holomorphic gauge kinetic function to be linear in the complex structure moduli $U_{X}$ and ${\cal U}_{i}$ , $f_{ab}(\Phi)\supset\delta_{ab}(S+U_{X}+{\cal U}_{i})$ . The complex structure moduli are associated with the SM gauge theory, which we will not be focusing on. The GS parameter $\delta^{\rm GS}_{X}$ characterizes the coupling of the anomalous gauge boson to the axion $\theta_{X}$ . The matter superfields in $K$ consist of all scalar fields that are not moduli and do not have Planck-sized vacuum expectation values (VEVs). The scalar components of $\varphi$ and $\varphi_{X}$ are neutral under the $U(1)_{X}$ symmetry and the modular group $\Gamma_{N}$ , respectively.

Calculating $K_{IJ}=\partial_{I}\partial_{J}$ from the Kähler potential Eq.(10), we obtain the kinetic terms for the scalar components of the supermultiplets which are approximated well for $M_{P}\gg\langle\varphi\rangle,\langle\varphi_{X}\rangle$ and $V_{X}=0$ as follows,

	$\displaystyle{\cal L}_{\rm kinetic}$	$\displaystyle\simeq$	$\displaystyle\frac{3M^{2}_{P}}{\langle-i\tau+i\bar{\tau}\rangle^{2}}\partial_{\mu}\bar{\tau}\partial^{\mu}\tau+\frac{M^{2}_{P}}{\langle U_{X}+\bar{U}_{X}\rangle^{2}}\partial_{\mu}\bar{U}_{X}\partial^{\mu}U_{X}$		(11)
		$\displaystyle+$	$\displaystyle\frac{M^{2}_{P}}{\langle S+\bar{S}-3\tilde{c}\ln(-i\tau+i\bar{\tau})\rangle^{2}}\partial_{\mu}\bar{S}\partial^{\mu}S+K_{\varphi\bar{\varphi}}\partial_{\mu}\bar{\varphi}\partial^{\mu}\varphi+K_{\varphi_{X}\bar{\varphi}_{X}}\partial_{\mu}\bar{\varphi}_{X}\partial^{\mu}\varphi_{X},$		(11)

where $K_{\varphi\bar{\varphi}}=K_{\varphi_{X}\bar{\varphi}_{X}}=1$ for canonically normalized scalar fields achieved by rescaling the fields $\varphi$ and $\varphi_{X}$ for given values of the VEVs of $\tau$ and $U_{X}$ . The $U(1)_{X}$ charged modulus $U_{X}$ and scalar field $\varphi_{X}$ can be decomposed as:

\displaystyle U_{X}=\frac{\rho_{X}}{2}+i\theta_{X},\qquad\varphi_{X}\big{|}_{\theta=\bar{\theta}=0}=\frac{1}{\sqrt{2}}e^{i\frac{A_{X}}{v_{X}}}(v_{X}+h_{X}),

(12)

where $Re[U_{X}]=\rho_{X}/2=1/g^{2}_{X}$ with $g_{X}$ being 4D $U(1)_{X}$ gauge coupling, $A_{X}$ , $v_{X}$ , and $h_{X}$ are the axion, VEV, and Higgs boson of scalar components, respectively. Due to the axionic shift symmetry, the kinetic terms of Eq.(11) for the axionic and size part of $U_{X}$ do not mix in perturbation theory, where any non-perturbative violations are small enough to be irrelevant, and the same goes for the axion and Higgs boson of the scalar components of $\varphi_{X}$ for $v_{X}\rightarrow\infty$ .

Since the matter superfields $\varphi$ and axio-dilaton $S$ transform as,

\displaystyle\varphi\rightarrow(c\tau+d)^{-k}\rho(\gamma)\varphi\,,\qquad S\rightarrow S-3\tilde{c}\ln(c\tau+d)\,,

(13)

where $\rho(\gamma)$ is the unitary representation of the modular group $\Gamma_{N}$ and $\tilde{c}$ is a constant, the transformation of the Kähler potential $K$ given in Eq.(9) leads us to

\displaystyle g(\tau)=\ln(c\tau+d)^{3}\,.

(14)

Generically, the transformation of $K$ in Eq.(9) incorporating Eq.(14) gives rise to a modular anomaly arising from $\delta S=-\tilde{c}\frac{1}{4}\int d^{4}xd^{2}\theta{\cal W}^{\alpha}{\cal W}_{\alpha}\,g(\tau)+{\rm h.c.}$ Derendinger:1991hq ,

\displaystyle-\frac{1}{8}\tilde{c}\Big{\{}\big{(}g(\tau)+g(\bar{\tau})\big{)}Q^{\mu\nu}Q_{\mu\nu}+i\big{(}g(\tau)-g(\bar{\tau})\big{)}Q^{\mu\nu}\tilde{Q}_{\mu\nu}\Big{\}}

(15)

where $\tilde{Q}_{\mu\nu}=\frac{1}{2}\varepsilon_{\mu\nu\rho\sigma}Q^{\rho\sigma}$ with associated gauge field strengths $Q$ , and the first term in the bracket represents the kinetic term for gauge bosons and the second term is the CP-odd term. After receiving a correction due to the modular transformation of $S$ in Eq.(13), The gauge kinetic function $f_{ab}$ is given at leading order by,

\displaystyle f^{\rm 1-loop}(\Phi)=\delta_{ab}(S+U_{X})-\tilde{c}\ln(c\tau+d)^{3}\,.

(16)

where the second term in the right hand side is the correction. It is worthwhile to notice that this correction cancels the modular anomaly Eq.(15) generated by $g(\tau),g(\bar{\tau})$ .

The modular invariance $W(\Phi)$ under the modular group $\Gamma_{N}$ ( $N\geq 2$ ) provides a strong restriction on the flavor structure Feruglio:2017spp . The superpotential $W(\Phi)$ can be expanded in power series of the multiplets $\varphi$ which are separated into brane sectors $\varphi_{(I)}$

\displaystyle W(\Phi)=\sum_{n}Y_{I_{1}...I_{n}}(\tau)\varphi_{(I_{1})}\cdot\cdot\cdot\varphi_{(I_{n})}\,,

(17)

where the functions $Y_{I_{1}...I_{n}}(\tau)$ are generically ⁷⁷7In Type II string orientifold compactifications, the Yukawa couplings have modular properties Cremades:2003qj . $\tau$ -dependent in type IIA intersecting D-brane models string_book ; Cremades:2004wa . The superpotential $W(\Phi)$ must have modular invariance under the transformation $W(\Phi)\rightarrow W(\Phi)e^{-g(\tau)}$ , where $g(\tau)$ is given by Eq.(14). To ensure this, we need to satisfy two conditions: (i) the matter superfields $\varphi_{I_{i}}$ of the brane sector $I_{i}$ should transform

\displaystyle\varphi_{(I_{i})}\rightarrow(c\tau+d)^{-k_{I_{i}}}\rho_{(I_{i})}(\gamma)\varphi_{(I_{i})}

(18)

in a representation $\rho_{(I_{i})}(\gamma)$ of the modular group $\Gamma_{N}$ , where $-k_{I_{i}}$ is the modular weight of sector $I_{i}$ , and (ii) the functions $Y_{I_{1}...I_{n}}(\tau)$ should be modular forms of weight $k_{Y}(n)$ transforming in the representation $\rho(\gamma)$ of $\Gamma_{N}$ ,

\displaystyle Y_{I_{1}...I_{n}}(\gamma\tau)=(c\tau+d)^{k_{Y}(n)}\rho(\gamma)Y_{I_{1}...I_{n}}(\tau)\,,

(19)

with the requirements

	$\displaystyle k_{Y}(n)-3=k_{I_{1}}+...+k_{I_{n}}\,,$
	$\displaystyle\rho(\gamma)\otimes\rho_{(I_{1})}\otimes\cdot\cdot\cdot\otimes\rho_{(I_{n})}\ni{\bf 1}\,.$		(20)

The weight of modular forms in the superpotential is corrected by the Kähler transformation in Eq.(9) due to the modular invariance of both the superpotential and Kähler potential. For example, a $\tau$ -independent superpotential for scalar potential can be simply constructed by the matter supermultiplets that belong to the untwisted sector in the orbifold compactification of type II string theory (see Eq.(37)). We will show an explicit example of the superpotential induced by the modular forms of level 3 in the section III.

II.2 Gauged U(1) anomaly-free

The 4D action given by Eq.(8) should also be $U(1)_{X}$ gauge invariant. Under the $U(1)_{X}$ gauge transformation $V_{X}\rightarrow V_{X}+i(\Lambda_{X}-\bar{\Lambda}_{X})$ , the matter superfields $\Phi_{X}$ and complex structure modulus $U_{X}$ transform as,

\displaystyle\Phi_{X}\rightarrow e^{iX\Lambda_{X}}\Phi_{X}\,,\qquad U_{X}\rightarrow U_{X}+i\frac{\delta^{\rm GS}_{X}}{16\pi^{2}}\Lambda_{X}\,,

(21)

where $\Lambda_{X}(\bar{\Lambda}_{X})$ are (anti-) chiral superfields parameterizing $U(1)_{X}$ transformation on the superspace. So the axionic modulus $\theta_{X}$ and axion $a_{X}$ have shift symmetries

\displaystyle\theta_{X}\rightarrow\theta_{X}-\frac{\delta^{\rm GS}_{X}}{16\pi^{2}}\xi_{X}\,,\qquad a_{X}\rightarrow a_{X}+\frac{\delta^{\rm GS}_{X}}{\delta^{Q}_{X}}f_{X}\xi_{X}\,,

(22)

where $\xi_{X}=-{\rm Re}\Lambda_{X}\big{|}_{\theta=\bar{\theta}=0}$ , $f_{X}=Xv_{X}$ is the axion decay constant, and $\delta^{Q}_{X}$ are anomaly coefficients defined in Eq.(25). Then, the $U(1)_{X}$ gauge field $A^{\mu}_{X}$ transforms as

\displaystyle A^{\mu}_{X}\rightarrow A^{\mu}_{X}-\partial^{\mu}\xi_{X}\,.

(23)

Since the gauged $U(1)_{X}$ is anomalous, the axion $a_{X}$ and axionic modulus $\theta_{X}$ couple to the (non)-Abelian Chern-Pontryagin densities for the SM gauge group in the compactification. In type II string vacuum, the $U(1)_{X}$ anomalies should be cancelled by appropriate shifts of Ramond-Ramond axions in the bulk Sagnotti:1992qw ; Ibanez:1998qp ; Poppitz:1998dj ; Lalak:1999bk . The 4D effective action of the axions, $\theta_{X}$ and $a_{X}$ , and its corresponding gauge field $A^{\mu}_{X}$ contains the followings Ahn:2017dpf ; Ahn:2016typ :

	$\displaystyle K_{U_{X}\bar{U}_{X}}\Big{(}\partial^{\mu}\theta_{X}-\frac{\delta^{\rm GS}_{X}}{16\pi^{2}}A^{\mu}_{X}\Big{)}^{2}-\frac{1}{4g^{2}_{X}}F^{\mu\nu}_{X}F_{X\mu\nu}-g_{X}\xi^{\rm FI}_{X}D_{X}+D_{X}g_{X}X\|\varphi_{X}\|^{2}$
	$\displaystyle+\|D_{\mu}\varphi_{X}\|^{2}+\theta_{X}{\rm Tr}(Q^{\mu\nu}\tilde{Q}_{\mu\nu})+\frac{a_{X}}{f_{X}}\frac{\delta^{Q}_{X}}{16\pi^{2}}{\rm Tr}(Q^{\mu\nu}\tilde{Q}_{\mu\nu})\,,$		(24)

where the gauge field strengths $Q=G,W,Y$ for $SU(3)_{C}$ , $SU(2)_{L}$ , and $U(1)_{Y}$ , respectively, and their gauge couplings are absorbed into their corresponding gauge field strengths. $F^{\mu\nu}_{X}$ is the $U(1)_{X}$ gauge field strength defined by $F^{\mu\nu}_{X}=\partial^{\mu}A^{\nu}_{X}-\partial^{\nu}A^{\mu}_{X}$ . In $|D_{\mu}\varphi_{X}|^{2}$ the scalar components of $\varphi_{X}$ couple to the $U(1)_{X}$ gauge boson, where the gauge coupling $g_{X}$ is absorbed into the gauge boson $A^{\mu}_{X}$ in the $U(1)_{X}$ gauge covariant derivative $D^{\mu}=\partial^{\mu}-iXA^{\mu}_{X}$ . The coefficients of the mixed $U(1)_{X}\times[SU(3)_{C}]^{2}$ , $U(1)_{X}\times[SU(2)_{L}]^{2}$ , and $U(1)_{X}\times[U(1)_{Y}]^{2}$ anomalies are given respectively by

\displaystyle\delta^{G}_{X}=2{\rm Tr}[XT^{2}_{SU(3)}]\,,\quad\delta^{W}_{X}=2{\rm Tr}[XT^{2}_{SU(2)}]\,,\quad\delta^{Y}_{X}=2{\rm Tr}[XY^{2}]\,.

(25)

Here $U(n)$ generators ( $n\geq 2$ ) are normalized according to ${\rm Tr}[T^{a}T^{b}]=\delta_{ab}/2$ , and for convenience, $\delta^{Y}_{X}=2{\rm Tr}[XY^{2}]$ is defined for hypercharge. The FI term ${\cal L}^{\rm FI}_{X}=-\xi^{\rm FI}_{X}\int d^{2}\theta V_{X}=-\xi^{\rm FI}_{X}g_{X}D_{X}$ with $D_{X}=g_{X}(\xi^{\rm FI}_{X}-X|\varphi_{X}|^{2})$ leads to D-term potential for the anomalous $U(1)_{X}$ ,

\displaystyle V_{D}=\frac{1}{U_{X}+\bar{U}_{X}}\big{(}-\xi^{\rm FI}_{X}+X|\varphi_{X}|^{2}\big{)}^{2}\,,

(26)

where $\xi^{\rm FI}_{X}$ is FI factor produced by expanding the Kähler potential Eq.(10) in components linear in $V_{X}$ and depends on the closed string modulus ${\rm Re}[U_{X}]=\rho_{X}/2$ . Since the FI term is controlled by the string coupling it can not be zero. The re-stabilization of VEVs by $\varphi_{X}$ necessarily implies spontaneous breaking of the anomalous $U(1)_{X}$ , which will be shown later.

The first, third, fourth, and fifth terms in Eq.(24) result from expanding the Kähler potential of Eq.(10). The first and sixth terms together, and the fifth and seventh terms in Eq.(24), are gauge invariant under the anomalous $U(1)_{X}$ gauge transformations of Eqs.(22, 23). The gauge invariant interaction Lagrangian is given by

\displaystyle{\cal L}^{\rm int}_{A\theta}=-A^{\mu}_{X}J^{\theta}_{\mu}+\theta_{X}{\rm Tr}(Q^{\mu\nu}\tilde{Q}_{\mu\nu})-A^{\mu}_{X}J^{X}_{\mu}+\frac{a_{X}}{f_{X}}\frac{\delta^{Q}_{X}}{16\pi^{2}}{\rm Tr}(Q^{\mu\nu}\tilde{Q}_{\mu\nu})

(27)

where the anomalous currents $J^{X}_{\mu}$ and $J^{\theta}_{\mu}$ coupling to the gauge boson $A^{\mu}_{X}$ (that is, $\partial_{\mu}J^{\mu}_{X}=\frac{\delta^{\rm GS}_{X}}{16\pi^{2}}{\rm Tr}(Q^{\mu\nu}\tilde{Q}_{\mu\nu})=-\partial_{\mu}J^{\mu}_{\theta}$ with $\delta^{\rm GS}_{X}=\alpha^{Q}_{X}\delta^{Q}_{X}$ ) are represented by $J^{\theta}_{\mu}=K_{U_{X}\bar{U}_{X}}\frac{\delta^{\rm GS}_{X}}{8\pi^{2}}\partial_{\mu}\theta_{X}$ and $J^{X}_{\mu}=-iX\varphi_{X}^{\dagger}\overleftrightarrow{\partial_{\mu}}\varphi_{X}$ .

Expanding Eq.(24) and setting $\theta_{X}=a_{\theta}/8\pi^{2}f_{\theta}$ with $f_{\theta}=\sqrt{\frac{2K_{U_{X}\bar{U}_{X}}}{(8\pi^{2})^{2}}}$ to canonically normalize, ${\cal L}^{\rm int}_{A\theta}$ becomes

	$\displaystyle\frac{1}{2}(\partial^{\mu}a_{\theta})^{2}+\frac{a_{\theta}}{8\pi^{2}f_{\theta}}{\rm Tr}(Q^{\mu\nu}\tilde{Q}_{\mu\nu})+\frac{1}{2}(\partial^{\mu}A_{X})^{2}+\frac{A_{X}}{f_{X}}\frac{\delta^{Q}_{X}}{16\pi^{2}}{\rm Tr}(Q^{\mu\nu}\tilde{Q}_{\mu\nu})$
	$\displaystyle-A^{\mu}_{X}(J^{X}_{\mu}+J^{\theta}_{\mu})+\frac{1}{2g^{2}_{X}}m^{2}_{X}A^{\mu}_{X}A_{X\mu}-\frac{1}{4g^{2}_{X}}F^{\mu\nu}_{X}F_{X\mu\nu}-\frac{g^{2}_{X}}{2}(\xi^{\rm FI}_{X}-X\|\varphi_{X}\|^{2})^{2}\,,$		(28)

where the gauge boson mass $m_{X}$ obtained by the super-Higgs mechanism is given by $m_{X}=\sqrt{2K_{U_{X}\bar{U}_{X}}(\delta^{\rm GS}_{X}/16\pi^{2})^{2}+2f^{2}_{X}}$ . Then the open string axion $a_{X}$ (decay constant $f_{X}$ ) is mixed linearly with the closed string $a_{\theta}$ (decay constant $f_{\theta}$ ):

\displaystyle\tilde{A}=\frac{a_{X}\frac{\delta^{\rm GS}_{X}}{2}f_{\theta}-a_{\theta}f_{X}}{\sqrt{f^{2}_{X}+(\frac{\delta^{\rm GS}_{X}}{2}f_{\theta})^{2}}}\approx a_{X}\,,\qquad G=\frac{a_{\theta}\frac{\delta^{\rm GS}_{X}}{2}f_{\theta}+A_{X}f_{X}}{\sqrt{f^{2}_{X}+(\frac{\delta^{\rm GS}_{X}}{2}f_{\theta})^{2}}}\approx a_{\theta}\,.

(29)

where the approximations are valid under the assumption that $f_{\theta}$ is much larger than $f_{X}$ . The gauged $U(1)_{X}$ absorbs one linear combination of $a_{X}$ and $a_{\theta}$ , denoted $G$ , giving it a string scale mass through the $U(1)_{X}$ gauge boson, while the other combination, $\tilde{A}\approx a_{X}$ , remains at low energies and contributes to the QCD axion. At energies below the scale $m_{X}$ , the gauge boson decouples, leaving behind an anomalous global $U(1)_{X}$ symmetry.

III Minimal model set-up

For our purpose, we take into account $\Gamma(3)$ modular symmetry, which gives the modular forms of level 3. The group $\Gamma_{3}$ is isomorphic to $A_{4}$ which is the symmetry group of the tetrahedron and the finite groups of the even permutation of four objects having four irreducible representations. Its irreducible representations are three singlets ${\bf 1},{\bf 1}^{\prime}$ , and ${\bf 1}^{\prime\prime}$ and one triplet ${\bf 3}$ with the multiplication rules ${\mathbf{3}}\otimes{\mathbf{3}}={\mathbf{3}}_{s}\oplus{\mathbf{3}}_{a}\oplus{\mathbf{1}}\oplus{\mathbf{1}}^{\prime}\oplus{\mathbf{1}}^{\prime\prime}$ and ${\mathbf{1}}^{\prime}\oplus{\mathbf{1}}^{\prime}={\mathbf{1}}^{\prime\prime}$ , where the subscripts $s$ and $a$ denote symmetric and antisymmetric combinations respectively. Let $(a_{1},a_{2},a_{3})$ and $(b_{1},b_{2},b_{3})$ denote the basis vectors for two ${\mathbf{3}}$ ’s. Then we have

$\displaystyle(a\otimes b)_{{\mathbf{3}}_{s}}$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{3}}\big{(}2a_{1}b_{1}-a_{2}b_{3}-a_{3}b_{2},~{}2a_{3}b_{3}-a_{1}b_{2}-a_{2}b_{1},~{}2a_{2}b_{2}-a_{3}b_{1}-a_{1}b_{3}\big{)}\,,$
$\displaystyle(a\otimes b)_{{\mathbf{3}}_{a}}$	$\displaystyle=$	$\displaystyle\big{(}a_{2}b_{3}-a_{3}b_{2},~{}a_{1}b_{2}-a_{2}b_{1},~{}a_{3}b_{1}-a_{1}b_{3}\big{)}\,,$
$\displaystyle(a\otimes b)_{\mathbf{1}}$	$\displaystyle=$	$\displaystyle a_{1}b_{1}+a_{2}b_{3}+a_{3}b_{2}\,,$
$\displaystyle(a\otimes b)_{{\mathbf{1}}^{\prime}}$	$\displaystyle=$	$\displaystyle a_{3}b_{3}+a_{1}b_{2}+a_{2}b_{1}\,,$
$\displaystyle(a\otimes b)_{{\mathbf{1}}^{\prime\prime}}$	$\displaystyle=$	$\displaystyle a_{2}b_{2}+a_{3}b_{1}+a_{1}b_{3}\,.$	(30)

The details of the $A_{4}$ group are shown in Appendix A. The modular forms $f(\tau)$ of level $3$ and weight $k$ , such as Eq.(19), are holomorphic functions of the complex variable $\tau$ with well-defined transformation properties

\displaystyle f(\gamma\tau)=(c\tau+d)^{k}f(\tau)\quad\gamma=\begin{pmatrix}a&b\\ c&d\end{pmatrix}\in\Gamma_{3}

(31)

with an integer $k\geq 0$ , under the group $\Gamma_{3}$ . The three linearly independent weight 2 and level-3 modular forms are given by Feruglio:2017spp

	$\displaystyle Y_{1}(\tau)=\frac{i}{2\pi}\Big{[}\frac{\eta^{\prime}(\frac{\tau}{3})}{\eta(\frac{\tau}{3})}+\frac{\eta^{\prime}(\frac{\tau+1}{3})}{\eta(\frac{\tau+1}{3})}+\frac{\eta^{\prime}(\frac{\tau+2}{3})}{\eta(\frac{\tau+2}{3})}-\frac{27\eta^{\prime}(3\tau)}{\eta(3\tau)}\Big{]}\,,$
	$\displaystyle Y_{2}(\tau)=\frac{-i}{\pi}\Big{[}\frac{\eta^{\prime}(\frac{\tau}{3})}{\eta(\frac{\tau}{3})}+\omega^{2}\frac{\eta^{\prime}(\frac{\tau+1}{3})}{\eta(\frac{\tau+1}{3})}+\omega\frac{\eta^{\prime}(\frac{\tau+2}{3})}{\eta(\frac{\tau+2}{3})}\Big{]}\,,$
	$\displaystyle Y_{3}(\tau)=\frac{-i}{\pi}\Big{[}\frac{\eta^{\prime}(\frac{\tau}{3})}{\eta(\frac{\tau}{3})}+\omega\frac{\eta^{\prime}(\frac{\tau+1}{3})}{\eta(\frac{\tau+1}{3})}+\omega^{2}\frac{\eta^{\prime}(\frac{\tau+2}{3})}{\eta(\frac{\tau+2}{3})}\Big{]}\,,$		(32)

where $\omega=-1/2+i\sqrt{3}/2$ and $\eta(\tau)$ is the Dedekind eta-function defined by

\displaystyle\eta(\tau)=q^{1/24}\prod^{\infty}_{n=1}(1-q^{n})\quad\text{with}~{}q\equiv e^{i2\pi\tau}~{}\text{and}~{}{\rm Im}(\tau)>0\,.

(33)

The Dedekind eta-function satisfies

\displaystyle\eta(-1/\tau)=\sqrt{-i\tau}\,\eta(\tau)\,,\qquad\eta(\tau+1)=e^{i\pi/12}\,\eta(\tau)\,.

(34)

The three linear independent modular functions transform as a triplet of $A_{4}$ , i.e. $Y^{(2)}_{\bf 3}=(Y_{1},Y_{2},Y_{3})$ . The $q$ -expansion of $Y_{i}(\tau)$ reads

	$\displaystyle Y_{1}(\tau)=1+12q+36q^{2}+12q^{3}+...$
	$\displaystyle Y_{2}(\tau)=-6q^{1/3}(17q+8q^{2}+...)$
	$\displaystyle Y_{3}(\tau)=-18q^{2/3}(1+2q+5q^{2}+...)\,.$		(35)

$Y^{(2)}_{\bf 3}$ is constrained by the relation,

\displaystyle(Y^{(2)}_{\bf 3}Y^{(2)}_{\bf 3})_{{\bf 1}^{\prime\prime}}=Y^{2}_{2}+2Y_{1}Y_{3}=0\,.

(36)

III.1 Modular invariant supersymmetric potential and a Nambu-Goldstone mode

Using Eqs.(17-20), we construct unique supersymmetric and modular invariant scalar potential by introducing minimal supermultiplets. Those include SM singlet fields $\chi_{0}$ ⁸⁸8The field $\chi_{0}$ can act as an inflaton Ahn:2017dpf . with modular weight three and $\chi(\tilde{\chi})$ with modular weight zero. Additionally, we have the usual two Higgs doublets $H_{u,d}$ with modular weight zero, which are responsible for electroweak (EW) symmetry breaking. The fields $\chi$ and $\tilde{\chi}$ are charged by $+1$ and $-1$ , respectively, and are ensured by the extended $U(1)_{X}$ symmetry due to the holomorphy of the superpotential. (If the seesaw mechanism Minkowski:1977sc is implemented, the field $\chi$ or $\tilde{\chi}$ may be responsible for the heavy neutrino mass term).

Under $k_{I}\times A_{4}\times U(1)_{X}$ with the modular weights $k_{I}$ according to Eq.(20), we assign the two Higgs doublets $H_{u,d}$ to be $(0,{\mathbf{1}},0)$ and three SM gauge singlets $\chi$ , $\tilde{\chi}$ , $\chi_{0}$ to be $(0,{\bf 1},+1)$ , $(0,{\bf 1},-1)$ , $(3,{\bf 1},0)$ , respectively ⁹⁹9As a consequence of $k_{I}\times A_{4}\times U(1)_{X}$ the other superpotential term $\kappa_{\alpha}L_{\alpha}H_{u}$ and the terms violating the lepton and baryon number symmetries are not allowed. Besides, dimension 6 supersymmetric operators like $Q_{i}Q_{j}Q_{k}L_{l}$ ( $i,j,k$ must not all be the same) are not allowed either, and stabilizing proton.. The $A_{4}$ -singlet $\chi_{0}$ field with modular weight three ensures that the functions $Y_{I_{1}...I_{n}}(\tau)$ are independent of $\tau$ . Then the supersymmetric scalar potential invariant under $G_{\rm SM}\times U(1)_{X}\times A_{4}$ is given at leading order by

\displaystyle W_{v}

\displaystyle=

\displaystyle g_{\chi_{0}}\chi_{0}\,H_{u}H_{d}+\chi_{0}(g_{\chi}\chi\tilde{\chi}-\mu^{2}_{\chi})\,\,,

(37)

where dimensionless coupling constants $g_{\chi_{0}}$ and $g_{\chi}$ are assumed to be equal to one, but are modified to Eq.(60) by considering all higher order terms induced by $\chi\tilde{\chi}$ combinations. Note that the PQ breaking parameter $\mu_{\chi}$ corresponds to the scale of the spontaneous symmetry breaking.

In the global SUSY limit, i.e. $M_{P}\rightarrow\infty$ , the scalar potential obtained by the $F$ - and $D$ -term of all fields is required to vanish. Then the relevant $F$ -term from Eq.(37) and $D$ -term of the scalar potential given by Eq.(26) reads

\displaystyle V_{F}^{\rm global}

\displaystyle=

\displaystyle\big{|}g_{\chi}\chi\tilde{\chi}-\mu^{2}_{\chi}\big{|}^{2}\,,\qquad V_{D}^{\rm global}=\frac{|X|^{2}g^{2}_{X}}{2}\Big{(}-\frac{\xi^{\rm FI}_{X}}{|X|}+|\chi|^{2}-|\tilde{\chi}|^{2}\Big{)}^{2}\,\,.

(38)

The scalar fields $\chi$ and $\tilde{\chi}$ have $X$ -charges $+1$ and $-1$ , respectively, i.e.,

\displaystyle\chi\rightarrow e^{+i\xi}\chi\,,\qquad\tilde{\chi}\rightarrow e^{-i\xi}\tilde{\chi}\,,

(39)

with a constant $\xi$ . So the potential $V_{\rm SUSY}$ has $U(1)_{X}$ symmetry. Since SUSY is preserved after the spontaneous breaking of $U(1)_{X}$ , the scalar potential in the limit of $M_{P}\rightarrow\infty$ vanishes at its ground states, i.e. $\langle V_{F}^{\rm global}\rangle=0$ and $\langle V_{D}^{\rm global}\rangle=0$ vanishing $F$ -term must have also vanishing $D$ -term. From the minimization of the $F$ -term scalar potential we obtain

\displaystyle\langle\chi\rangle=\langle\tilde{\chi}\rangle=\frac{v_{\chi}}{\sqrt{2}}\qquad\text{with}~{}\mu_{\chi}=v_{\chi}\sqrt{\frac{g_{\chi}}{2}}\,.

(40)

where we have assumed $\langle\chi\rangle,\langle\tilde{\chi}\rangle\gg\langle H_{u,d}\rangle$ . The above supersymmetric solution is taken by the $D$ -flatness condition for Ahn:2016hbn ; Ahn:2017dpf

\displaystyle\xi^{\rm FI}_{X}=0\,,\qquad\langle\chi\rangle=\langle\tilde{\chi}\rangle\,.

(41)

The tension between $\langle\chi\rangle=\langle\tilde{\chi}\rangle$ and $\xi^{\rm FI}_{X}\neq 0$ arises because the FI term cannot be cancelled, unless the VEV of flux in the FI term is below the string scale Achucarro:2006zf ; Burgess:2003ic . The FI term acts as an uplifting potential,

\displaystyle\xi^{\rm FI}_{X}=M^{2}_{P}\frac{\delta^{\rm GS}_{X}}{16\pi^{2}}\frac{\Delta\rho}{\rho_{0}}\,,

(42)

where $\Delta\rho=\rho_{X}-\rho_{0}$ , which raises the Anti-de Sitter minimum to the de Sitter minimum Burgess:2003ic . To achieve this, the $F$ -term must necessarily break SUSY for the $D$ -term to act as an uplifting potential. The PQ scale $\mu_{\chi}$ can be determined by taking into account both the SUSY-breaking effect, which lifts up the flat direction, and supersymmetric next-leading-order Planck-suppressed terms Nanopoulos:1983sp ; Ahn:2016hbn ; Ahn:2017dpf . The supersymmetric next-to-leading order term invariant under $A_{4}\times U(1)_{X}$ satisfying Eq.(20) are given by

\displaystyle\Delta W_{v}\simeq\frac{\alpha}{M^{2}_{P}}\chi_{0}(\chi\tilde{\chi})^{2}\,,

(43)

where $\alpha$ is assumed to be a real-valued constant being of unity. Since soft SUSY-breaking terms are already present at the scale relevant to flavor dynamics, the scalar potential for $\chi$ , $\tilde{\chi}$ at leading order read

\displaystyle V(\chi,\tilde{\chi})\simeq-\alpha_{1}m^{2}_{3/2}|\chi|^{2}-\alpha_{2}m^{2}_{3/2}|\tilde{\chi}|^{2}+\alpha^{2}\frac{|\chi|^{4}|\tilde{\chi}|^{4}}{M^{4}_{P}}\,,

(44)

where $m_{3/2}$ represents soft SUSY-breaking mass, and $\alpha_{1},\alpha_{2}$ are real-valued constants. It leads to the PQ breaking scale (equivalently, the seesaw scale),

\displaystyle\mu_{\chi}\simeq\big{(}\frac{g^{6}_{\chi}\alpha_{1}\alpha_{2}}{16\alpha^{4}}\big{)}^{\frac{1}{12}}\big{(}m_{3/2}M^{2}_{P}\big{)}^{\frac{1}{3}}\,,

(45)

indicating that $\mu_{\chi}$ lies within the range of approximately $1.2\times 10^{13}$ to $1.7\times 10^{14}$ GeV (or $2.6\times 10^{13}$ to $1.2\times 10^{15}$ GeV) for $m_{3/2}$ values ranging from 1 TeV to $10^{3}$ TeV (or from 10 TeV to $10^{6}$ TeV) for $\alpha_{1}$ and $\alpha_{2}$ of order unity.

The model includes the SM gauge singlet scalar fields $\chi$ and $\tilde{\chi}$ charged under $U(1)_{X}$ , which have interactions invariant under $G_{\rm SM}\times U(1)_{X}\times A_{4}$ with the transformations Eq.(9). These interactions result in a chiral symmetry, which is reflected in the form of the kinetic and Yukawa terms, as well as the scalar potential $V_{\rm SUSY}$ in the SUSY limit:

\displaystyle{\cal L}

\displaystyle\supset

\displaystyle\partial_{\mu}\chi^{\ast}\partial^{\mu}\chi+\partial_{\mu}\tilde{\chi}^{\ast}\partial^{\mu}\tilde{\chi}+{\cal L}_{Y}-V_{\rm SUSY}+{\cal L}_{\vartheta}+\overline{\psi}\,i\!\!\not\!\partial\psi+\frac{1}{2}\overline{N}\,i\!\!\not\!\partial N+\frac{1}{2}\overline{\nu}\,i\!\!\not\!\partial\nu\,,

(46)

where $\psi$ denotes Dirac fermions, and $V_{\rm SUSY}$ is replaced by $V_{\rm total}$ when SUSY breaking effects are considered. The above kinetic terms for $\chi(\tilde{\chi})$ are canonically normalized from the Kähler potential Eq.(10). Here four component Majorana spinors ( $N^{c}=N$ and $\nu^{c}=\nu$ ) are used. The global $U(1)_{X}$ PQ symmetry guarantees the absence of bare mass term in the Yukawa Lagrangian ${\cal L}_{Y}$ in Eq.(46). The QCD Lagrangian has a CP-violating term

\displaystyle{\cal L}_{\vartheta}

\displaystyle=

\displaystyle\vartheta_{\rm QCD}\,\frac{g^{2}_{s}}{32\pi^{2}}\,G^{a\mu\nu}\tilde{G}^{a}_{\mu\nu}

(47)

where $g_{s}$ stands for the gauge coupling constant of $SU(3)_{C}$ , and $G^{a\mu\nu}$ is the color field strength tensor and its dual $\tilde{G}^{a}_{\mu\nu}=\frac{1}{2}\varepsilon_{\mu\nu\rho\sigma}G^{a\mu\nu}$ (here $a$ is an $SU(3)$ -adjoint index), coming from the strong interaction. After obtaining VEV $\langle\chi\rangle\neq 0$ , which generates the heavy neutrino masses given by Eq.(53), the PQ $U(1)_{X}$ symmetry breaks spontaneously at a much higher scale than EW scale. This is manifested through the existence of the Nambu-Goldstone (NG) mode $A_{X}$ , which interacts with ordinary quarks and leptons via Yukawa interactions, see Eqs.(81, 94, 121). To extract the associated boson resulting from spontaneous breaking of $U(1)_{X}$ , we set the decomposition of complex scalar fields Ahn:2014gva ; Ahn:2016hbn ; Ahn:2018cau as follows

\displaystyle\chi=\frac{v_{\chi}}{\sqrt{2}}e^{i\frac{A_{X}}{u_{\chi}}}\left(1+\frac{h_{\chi}}{u_{\chi}}\right)\,,\quad\,\tilde{\chi}=\frac{v_{\tilde{\chi}}}{\sqrt{2}}e^{-i\frac{A_{X}}{u_{\chi}}}\left(1+\frac{h_{\tilde{\chi}}}{u_{\chi}}\right)\qquad\text{with}~{}u_{\chi}=\sqrt{v^{2}_{\chi}+v^{2}_{\tilde{\chi}}}\,,

(48)

in which $A_{X}$ is the NG mode and we set $v_{\chi}=v_{\tilde{\chi}}$ and $h_{\chi}=h_{\tilde{\chi}}$ in the supersymmetric limit. The derivative coupling of NG boson $A_{X}$ arises from the kinetic term

\displaystyle\partial_{\mu}\chi^{\ast}\partial^{\mu}\chi+\partial_{\mu}\tilde{\chi}^{\ast}\partial^{\mu}\tilde{\chi}=\frac{1}{2}(\partial_{\mu}A_{X})^{2}\Big{(}1+\frac{h_{\chi}}{u_{\chi}}\Big{)}^{2}+\frac{1}{2}(\partial_{\mu}h_{\chi})^{2}\,.

(49)

Performing $u_{\chi}\rightarrow\infty$ , the NG mode $A_{X}$ , whose interaction is determined by symmetry, is distinguished from the radial mode $h_{\chi}$ , which is invariant under the symmetry $U(1)_{X}$ .

III.2 Modular invariant Yukawa superpotentials and Anomaly coefficients

By introducing just two $A_{4}$ -singlet fields, $\chi$ and $\tilde{\chi}$ , with modular weight zero and charged under $U(1)_{X}$ by $+1$ and $-1$ , respectively, and using economic weight modular forms, we construct Yukawa superpotentials that are invariant under $G_{\rm SM}\times U(1)_{X}\times A_{4}$ satisfying Eq.(20). This approach can explain the observed hierarchy of fermion masses and mixing given by the Cabibbo-Kobayashi-Maskawa (CKM) matrix for quarks as well as by Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix for leptons. Furthermore, the approach provides a solution to the strong CP problem by breaking down the $U(1)_{X}$ flavor symmetry. Since the modular weights of the fields $\chi(\tilde{\chi})$ are zero, any additive correction terms induced by higher weight modular forms are forbidden in the superpotential (see Eqs.(37, 50, 53)). However, higher-order corrections arising from the combination $\chi\tilde{\chi}$ are allowed, but they do not modify the leading-order flavor structure.

Now let us assign $A_{4}\times U(1)_{X}$ representations and quantum numbers as well as modular weights $k_{I}$ to the SM quarks and leptons including SM gauge singlet Majorana neutrinos as presented in Table-1 ¹⁰¹⁰10All fields appearing in Table-1 are left-handed particles/antiparticles.. Here, three quark $SU(2)_{L}$ doublets and three up-type quark singlets are denoted as $Q_{i(=1,2,3)}$ and $(u^{c},c^{c},t^{c})$ , respectively. ${\cal D}^{c}=\{d^{c},s^{c},b^{c}\}$ represents the down-type quark singlets.

Table 1: Representations and quantum numbers of the quark fields under

G_{\rm SM}\times A_{4}\times U(1)_{X}

and modular weight

k_{I}

according to Eq.(20). In

({\cal Q}_{1},{\cal Q}_{2})_{Y}

G_{\rm SM}

{\cal Q}_{1}

and

{\cal Q}_{2}

are the representations under

SU(3)_{C}

and

SU(2)_{L}

respectively, and the script

Y

denotes the

U(1)

hypercharge.

Field	$Q_{1}$	$Q_{2}$	$Q_{3}$	${\cal D}^{c}$	$u^{c}$	$c^{c}$	$t^{c}$
$G_{\rm SM}$	$(3,2)_{1/6}$	$(3,2)_{1/6}$	$(3,2)_{1/6}$	$(3,1)_{1/3}$	$(3,1)_{-2/3}$	$(3,1)_{-2/3}$	$(3,1)_{-2/3}$
$A_{4}$	$\mathbf{1}$	$\mathbf{1}^{\prime\prime}$	$\mathbf{1^{\prime}}$	$\mathbf{3}$	$\mathbf{1}$	$\mathbf{1}^{\prime}$	$\mathbf{1^{\prime\prime}}$
$k_{I}$	$0$	$0$	$0$	$-3$	$-3$	$-3$	$3$
$U(1)_{X}$	$f_{b}-f_{d}$	$f_{b}-f_{s}$	$0$	$-f_{b}$	$f_{d}-f_{b}-f_{u}$	$f_{s}-f_{b}-f_{c}$	$0$

Then the quark Yukawa superpotential invariant under $G_{\rm SM}\times A_{4}\times U(1)_{X}$ with modular forms is sewed with ${\cal F}=\{\chi$ or $\tilde{\chi}\}$ through Eq.(2) as

$\displaystyle W_{q}$	$\displaystyle=$	$\displaystyle\alpha_{t}^{(0)}\,t^{c}Q_{3}H_{u}+\alpha_{c}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{c}\|}Y^{(6)}_{\bf 1}c^{c}Q_{2}H_{u}+\alpha_{u}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{u}\|}Y^{(6)}_{\bf 1}u^{c}Q_{1}H_{u}$	(50)
	$\displaystyle+$	$\displaystyle\alpha_{b}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{b}\|}(Y^{(6)}_{\bf 3}{\cal D}^{c})_{{\bf 1}^{\prime\prime}}Q_{3}H_{d}+\alpha_{s}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{s}\|}(Y^{(6)}_{\bf 3}{\cal D}^{c})_{{\bf 1}^{\prime}}Q_{2}H_{d}$
	$\displaystyle+$	$\displaystyle\alpha_{d}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{d}\|}(Y^{(6)}_{\bf 3}{\cal D}^{c})_{{\bf 1}}Q_{1}H_{d}+W^{(h)}_{q}\,,$

where $\alpha_{i}^{(0)}$ denotes coefficient at leading order and $W^{(h)}_{q}$ stand for higher order contributions, which are simply contructed by the leadging order operators in Eq.(50) multiplied by $\sum^{\infty}_{n=1}\Big{(}\frac{\chi\tilde{\chi}}{\Lambda^{2}}\Big{)}^{n}$ . Note that all Yukawa coefficients in the above suprpotential, $\alpha_{i}^{(0)}$ , are assumed to be complex numbers with an absolute value of unity. Since it is hard to reproduce the experimental data of fermion masses and mixing with Yukawa terms constructed with modular forms of weight 4 in quark and charged-lepton sectors in this model, we take into account Yukawa terms with modular forms of weight 6 which are decomposed as ${\bf 1}\oplus{\bf 3}\oplus{\bf 3}$ under $A_{4}$ given explicitly by Feruglio:2017spp

	$\displaystyle Y^{(6)}_{\bf 1}=Y^{3}_{1}+Y^{3}_{2}+Y^{3}_{3}-3Y_{1}Y_{2}Y_{3}\,$
	$\displaystyle Y^{(6)}_{{\bf 3},1}=(Y^{3}_{1}+2Y_{1}Y_{2}Y_{3},Y^{2}_{1}Y_{2}+2Y^{2}_{2}Y_{3},Y^{2}_{1}Y_{3}+2Y^{2}_{3}Y_{2})$
	$\displaystyle Y^{(6)}_{{\bf 3},2}=(Y^{3}_{3}+2Y_{1}Y_{2}Y_{3},Y^{2}_{2}Y_{1}+2Y^{2}_{1}Y_{2},Y^{2}_{3}Y_{2}+2Y^{2}_{2}Y_{1})\,.$		(51)

In the above superpotential only the top quark operator is renormalizable and does not contain a modular form, leading to the top quark mass as the pole mass, while the other quark operators driven by $\chi$ (or $\tilde{\chi})$ are dependent on modular forms. Using modular forms of weight 6, $Y^{(6)}_{\bf 1}$ and $Y^{(6)}_{\bf 3}$ , with the quark fields charged under $A_{4}\times U(1)_{X}$ , which does not allow mixing among up-type quarks, the off-diagonal entries in the up-type quark mass matrix are forbidden, as indicated in Eq.(68). From the above superpotential the effective Yukawa couplings of quarks can be visualized as functions of the SM gauge-singlet fields $\chi(\tilde{\chi})$ and modular forms $Y^{(6)}_{\bf 1(3)}$ , except for the top Yukawa coupling (see the details given in Sec. IV).

According to the quantum numbers of the quark sectors as in Table-1, the color anomaly coefficient of $U(1)_{X}\times[SU(3)_{C}]^{2}$ defined as $N_{C}\equiv 2{\rm Tr}[X_{\psi}T^{2}_{SU(3)_{C}}]$ reads

\displaystyle N_{C}=-(f_{u}+f_{c}+f_{d}+f_{s}+f_{b})\,.

(52)

Note that $U(n)$ generators ( $n\geq 2$ ) are normalized according to ${\rm Tr}[T^{a}T^{b}]=\delta^{ab}/2$ . The $U(1)_{X}$ is broken down to its discrete subgroup $Z_{N_{\rm DW}}$ in the backgrounds of QCD instanton, and the quantity $N_{C}$ (non-zero integer) is given by the axionic domain-wall number $N_{\rm DW}$ . At the QCD phase transition, each axionic string becomes the edge to $N_{\rm DW}$ domain-walls, and the process of axion radiation stops. To avoid the domain-wall problem one should consider $N_{\rm DW}=1$ or the PQ phase transition occurred during (or before) inflation for $N_{\rm DW}>1$ .

Next, we turn to the lepton sector, where the fields are charged under $G_{\rm SM}\times A_{4}\times U(1)_{X}$ with modular weight $k_{I}$ . Remark that the sterile neutrinos $N^{c}$ (which interact with gravity) are introduced (i) to solve the anomaly-free condition of $U(1)\times[gravity]^{2}$ , (ii) to explain the small active neutrino masses via the seesaw mechanism, and (iii) to provide a theoretically well-motivated PQ symmetry breaking scale.

Table 2: Representations and quantum numbers of the lepton fields under

G_{\rm SM}\times A_{4}\times U(1)_{X}

and modular weight

k_{I}

determined according to Eq.(20).

Field	$L_{e}$	$L_{\mu}$	$L_{\tau}$	$e^{c}$	$\mu^{c}$	$\tau^{c}$	$N^{c}$
$G_{\rm SM}$	$(1,2)_{-1/2}$	$(1,2)_{-1/2}$	$(1,2)_{-1/2}$	$(1,1)_{1}$	$(1,1)_{1}$	$(1,1)_{1}$	$(1,1)_{0}$
$A_{4}$	$\mathbf{1}$	$\mathbf{1}^{\prime}$	$\mathbf{1}^{\prime\prime}$	$\mathbf{1}$	$\mathbf{1^{\prime\prime}}$	$\mathbf{1^{\prime}}$	$\mathbf{3}$
$k_{I}$	$\frac{5}{2}$	$\frac{5}{2}$	$\frac{5}{2}$	$-\frac{11}{2}$	$-\frac{11}{2}$	$-\frac{11}{2}$	$-\frac{3}{2}$
$U(1)_{X}$	$\frac{1}{2}-g_{e}$	$\frac{1}{2}-g_{\mu}$	$\frac{1}{2}-g_{\tau}$	$g_{e}-\frac{1}{2}-f_{e}$	$g_{\mu}-\frac{1}{2}-f_{\mu}$	$g_{\tau}-\frac{1}{2}-f_{\tau}$	$-\frac{1}{2}$

In Table 2, the representations and quantum numbers of the lepton fields as well as modular weight $k_{f}$ determined along with Eq.(20) are presented. Here, $L_{e}$ , $L_{\mu}$ , and $L_{\tau}$ denote $SU(2)_{L}$ lepton doublets and $e^{c}$ , $\mu^{c}$ , and $\tau^{c}$ are three charged lepton singlets. The field $N^{c}$ represents the right-handed $SU(2)_{L}$ singlet neutrino, which is introduced to generate active neutrino masses via canonical seesaw mechanism Minkowski:1977sc .

We note that mixing between different charged-leptons does not occur when lepton Yukawa superpotential is economically constructed with modular forms of weight 6 prevents, which leads to the diagonal form of the charged lepton mass matrix as can be seen in Eq.(91)). In contrast, modular forms $Y^{(2)}_{\bf 3}$ , $Y^{(6)}_{\bf 1}$ , and $Y^{(6)}_{\bf 3}$ are used to construct neutrino mass matrices ¹¹¹¹11By selecting appropriate modular weight of particle contents, lower weight modular forms can be used, such as (i) $Y^{(2)}_{\bf 3}$ in the Dirac neutrino sector, and $Y^{(4)}_{\bf 1(1^{\prime},1^{\prime\prime})}$ and $Y^{(4)}_{\bf 3}$ in the Majorana neutrino sector. However, this leads to additional interactions, including $\frac{1}{2}Y^{(4)}_{\bf 1}(N^{c}N^{c})_{\bf 1}\chi$ , $\frac{1}{2}Y^{(4)}_{{\bf 1}^{\prime}}(N^{c}N^{c})_{{\bf 1}^{\prime\prime}}\chi$ , $\frac{1}{2}Y^{(4)}_{{\bf 1}^{\prime\prime}}(N^{c}N^{c})_{{\bf 1}^{\prime}}\chi$ , and $\frac{1}{2}Y^{(4)}_{\bf 3}(N^{c}N^{c})_{\bf 3}\chi$ . (ii) Another option is to use $Y^{(2)}_{\bf 3}$ in the Dirac neutrino sector and no modular form in the Majorana neutrino sector, which results in only $\frac{1}{2}(N^{c}N^{c})_{\bf 1}\chi$ and degenerate heavy Majorana neutrino mass states at the seesaw scale. However, we have found that this approach is difficult to reconcile with experimental neutrino data. Then the Yukawa superpotential for lepton invariant under $G_{\rm SM}\times A_{4}\times U(1)_{X}$ with economic modular forms are sewed with ${\cal F}=\{\chi$ or $\tilde{\chi}\}$ through Eq.(2), respectively, as

$\displaystyle W_{\ell\nu}$	$\displaystyle=$	$\displaystyle\alpha_{\tau}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{\tau}\|}Y^{(6)}_{\bf 1}\tau^{c}L_{\tau}H_{d}+\alpha_{\mu}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{\mu}\|}Y^{(6)}_{\bf 1}\mu^{c}L_{\mu}H_{d}+\alpha_{e}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{e}\|}Y^{(6)}_{\bf 1}e^{c}L_{e}H_{d}$	(53)
	$\displaystyle+$	$\displaystyle\beta_{1}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|g_{e}\|}(Y^{(2)}_{\bf 3}N^{c})_{{\bf 1}}L_{e}H_{u}+\beta_{2}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|g_{\mu}\|}(Y^{(2)}_{\bf 3}N^{c})_{{\bf 1}^{\prime\prime}}L_{\mu}H_{u}$
	$\displaystyle+$	$\displaystyle\beta_{3}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|g_{\tau}\|}(Y^{(2)}_{\bf 3}N^{c})_{{\bf 1}^{\prime}}L_{\tau}H_{u}+\gamma_{1}^{(0)}\frac{1}{2}Y^{(6)}_{\bf 1}(N^{c}N^{c})_{{\bf 1}}\chi+\gamma_{2}^{(0)}\frac{1}{2}Y^{(6)}_{\bf 3}(N^{c}N^{c})_{{\bf 3}}\chi+W^{(h)}_{l\nu},$

where $\alpha_{i}^{(0)},\beta_{i}^{(0)},\gamma_{i}^{(0)}$ denote coefficients at leading order and $W^{(h)}_{l\nu}$ stands for higher order contributions triggered by the combination $\chi\tilde{\chi}$ . Like in the quark sector, the Yukawa coefficients in the above superpotential, such as $\alpha_{i}^{(0)},\beta_{i}^{(0)},\gamma_{i}^{(0)}$ , are assumed to be complex numbers with an absolute value of unity. In the above superpotential, the charged-lepton and Dirac neutrino parts have three distinct Yukawa terms each, with their common modular forms being $Y^{(6)}_{\bf 1}$ and $Y^{(2)}_{\bf 3}$ , respectively. Each term involves ${\cal F}/\Lambda$ to the power of an appropriate $U(1)_{X}$ quantum number. The flavored $U(1)_{X}$ PQ symmetry allows for two renormalizable terms for the right-handed neutrino $N^{c}$ , which implement the seesaw mechanism Minkowski:1977sc by making the VEV $\langle\chi\rangle$ large. The details on how the active neutrino masses and mixing are predicted will be presented in Sec. IV.3.

Nonperturbative quantum gravitational anomaly effects Kamionkowski:1992mf violate the conservation of the corresponding current, $\partial_{\mu}J^{\mu}_{X}\propto R\tilde{R}$ , where $R$ is the Riemann tensor and $\tilde{R}$ is its dual, and make the axion solution to the strong CP problem problematic. To consistently couple gravity to matter charged under $U(1)_{X}$ , the mixed-gravitational anomaly $U(1)_{X}\times[{\rm gravity}]^{2}$ (related to the color anomaly $U(1)_{X}\times[SU(3)_{C}]^{2}$ ) must be cancelled, as shown in Refs. Ahn:2016hbn ; Ahn:2018cau ; Ahn:2021ndu , which leads to the relation,

\displaystyle 3N_{C}=f_{e}+f_{\mu}+f_{\tau}+g_{e}+g_{\mu}+g_{\tau}\,.

(54)

Thus the choice of $U(1)_{X}$ charge for ordinary quarks and leptons is strictly restricted.

Below the $U(1)_{X}$ symmetry breaking scale (here, equivalent to the seesaw scale) the effective interactions of QCD axion with the weak and hypercharge gauge bosons and with the photon are expressed through the chiral rotation of Eq.(62), respectively, as

	$\displaystyle{\cal L}^{WY}_{A}$	$\displaystyle=$	$\displaystyle\frac{A_{X}}{f_{A}}\frac{1}{32\pi^{2}}\Big{\{}g^{2}_{W}\,N_{W}\,W^{\mu\nu}\tilde{W}_{\mu\nu}+g^{2}_{Y}\,N_{Y}\,Y^{\mu\nu}\tilde{Y}_{\mu\nu}\Big{\}}\,,$		(55)
	$\displaystyle{\cal L}^{\gamma}_{A}$	$\displaystyle=$	$\displaystyle\frac{A_{X}}{f_{A}}\frac{e^{2}}{32\pi^{2}}\,E\,F^{\mu\nu}\tilde{F}_{\mu\nu}\,,$		(56)

where $g_{W}$ , $g_{Y}$ , and $e$ stand for the gauge coupling constant of $SU(2)_{L}$ , $U(1)_{Y}$ , and $U(1)_{EM}$ , respectively; their corresponding gauge field strengths $W^{\mu\nu},Y^{\mu\nu}$ , and $F^{\mu\nu}$ with their dual forms $\tilde{W}_{\mu\nu},\tilde{Y}_{\mu\nu}$ , and $\tilde{F}_{\mu\nu}$ , respectively. Here $N_{W}\equiv 2{\rm Tr}[X_{\psi_{f}}T^{2}_{SU(2)}]$ and $N_{Y}\equiv 2{\rm Tr}[X_{\psi_{f}}(Q_{f}^{Y})^{2}]$ are the anomaly coefficients of $U(1)_{X}\times[SU(2)_{L}]^{2}$ and $U(1)_{X}\times[U(1)_{Y}]^{2}$ , respectively. And the electromagnetic anomaly coefficient $E$ of $U(1)_{X}\times[U(1)_{EM}]^{2}$ defined by $E=2\sum_{\psi_{f}}X_{\psi_{f}}(Q^{\rm em}_{\psi_{f}})^{2}$ with $Q^{\rm em}_{\psi_{f}}$ being the $U(1)_{\rm EM}$ charge of field $\psi_{f}$ is expressed as

\displaystyle E

\displaystyle=

\displaystyle N_{W}+N_{Y}=-2(f_{e}+f_{\mu}+f_{\tau})-\frac{2}{3}(4f_{u}+4f_{c}+f_{d}+f_{s}+f_{b})\,.

(57)

The physical quantities of QCD axion, such as axion mass $m_{a}$ and axion-photon coupling $g_{a\gamma\gamma}$ , are dependent on the ratio of electromagnetic anomaly coefficient $E$ to color one $N_{C}$ . The value of $E/N_{C}$ is determined in terms of the $X$ -charges for quarks and leptons by the relation,

	$\displaystyle\frac{E}{N_{C}}$	$\displaystyle=$	$\displaystyle\frac{2(f_{e}+f_{\mu}+f_{\tau})+\frac{2}{3}(4f_{u}+4f_{c}+f_{d}+f_{s}+f_{b})}{f_{u}+f_{c}+f_{d}+f_{s}+f_{b}}\,$		(58)
		$\displaystyle=$	$\displaystyle\frac{6(f_{e}+f_{\mu}+f_{\tau})+2(4f_{u}+4f_{c}+f_{d}+f_{s}+f_{b})}{-f_{e}-f_{\mu}-f_{\tau}-g_{e}-g_{\mu}-g_{\tau}}\,,$		(58)

where the first and second equality follow from Eqs.(52) and (54), respectively. Our model with a specific value of $E/N_{C}$ can be tested by ongoing experiments such as KLASH Alesini:2019nzq and FLASH Gatti:2021cel ( see Eq.(86) and Fig.1 and 2) by considering the scale of $U(1)_{X}$ breakdown induced by Eq.(45).

Compared to conventional $A_{4}$ symmetry models resulting in tribimaximal Ma:2001dn or nearly tribimaximal Ahn:2012cg mixing in the neutrino sector, the modular invariant model leads to neutrino mixing without the need for special breaking patterns and the introduction of multiple scalar fields. Our model can be uniquely realized for quark sector by assigning $A_{4}\times U(1)_{X}$ quantum numbers to matter fields with appropriate modular forms based on Eq.(2). Some comments are worth noting. (i) By selecting the appropriate modular weight for the right-handed down-type quark fields, it is possible to construct down-type quark Yukawa superpotential with lower modular weight forms $Y^{(2)}_{\bf 3}$ or $Y^{(4)}_{\bf 3}$ while keeping the same up-type quark Yukawa superpotential given in Eq.(50). However, it is hard to reproduce the experimental data for quark masses and mixing hierarchies in this way due to the limited number of parameters. (ii) Unlike the case in Table-1, the quark $SU(2)_{L}$ doublets and singlets can be assigned to $A_{4}$ triplets and singlets by choosing appropriate modular weight forms and $U(1)_{X}$ quantum numbers, respectively. In this case, the quark mass hierarchies can be realized in the limit of $\langle\tau\rangle=i\infty$ (i.e. $Y_{1}\rightarrow 1$ , $Y_{2}\rightarrow 0$ , and $Y_{3}\rightarrow 0$ ), whereas it is hard to reproduce the CKM mixing angles since additive correction terms induced by higher weight modular forms are forbidden by the modular weight zero of $\chi(\tilde{\chi})$ fields. (iii) In the opposite scenario where the quark $SU(2)_{L}$ doublets and singlets are assigned to $A_{4}$ singlets and triplets, respectively, it is not possible to account for the observed quark mass hierarchy due to the charge assignment of $U(1)_{X}$ . (iv) For leptons, unlike the case in Table-2, the left-handed charged lepton $SU(2)_{L}$ doublets $L$ can be assigned to the $A_{4}$ triplet and their $U(1)_{X}$ quantum numbers are taken to be $\frac{1}{2}-g_{l}$ , whereas $SU(2)_{L}$ singlets ( $e^{c}$ , $\mu^{c}$ , $\tau^{c}$ ) are assigned to the $A_{4}$ singlets (1,1 ${}^{{}^{\prime\prime}}$ ,1 ${}^{{}^{\prime}}$ ) and $U(1)_{X}$ quantum numbers are taken to be $(g_{l}-f_{e}-\frac{1}{2},g_{l}-f_{\mu}-\frac{1}{2},g_{l}-f_{\tau}-\frac{1}{2})$ . To generate neutrino mass through the seesaw mechanism, $N^{c}$ is assigned to the $A_{4}$ triplet, and $U(1)_{X}$ quantum number is taken to be $-\frac{1}{2}$ . In this case, we have the freedom to select the weights. For instance, we can choose the following weights: $k_{L}=\frac{5}{2}$ , $k_{e^{c}}=k_{\mu^{c}}=k_{\tau^{c}}=-\frac{3}{2}$ , and $k_{N^{c}}=\frac{1}{2}$ . Then the lepton Yukawa superpotential reads

	$\displaystyle W_{\ell\nu}$	$\displaystyle=$	$\displaystyle\Big{[}\alpha_{\tau}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{\tau}\|}(Y^{(2)}_{\bf 3}L)_{{\bf 1}^{\prime\prime}}\tau^{c}+\alpha_{\mu}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{\mu}\|}(Y^{(2)}_{\bf 3}L)_{{\bf 1}^{\prime}}\mu^{c}+\alpha_{e}^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|f_{e}\|}(Y^{(2)}_{\bf 3}L)_{\bf 1}e^{c}\Big{]}H_{d}$		(59)
		$\displaystyle+$	$\displaystyle\beta^{(0)}\Big{(}\frac{\cal F}{\Lambda}\Big{)}^{\|g_{l}\|}(N^{c}L)_{{\bf 1}}H_{u}+\gamma^{(0)}\frac{1}{2}(Y^{(2)}_{\bf 3}N^{c}N^{c})_{{\bf 1}}\chi+...\,,$		(59)

where dots stand for higher order contributions triggered by the combination $\chi\tilde{\chi}$ . It is worth noting that the above superpotential enables mixing between different charged leptons, analogous to the down-type quark sector. Additionally, the Dirac neutrino Yukawa matrix, denoted as $m_{D}$ , exhibits a proportional relationship to $m_{D}^{\dagger}m_{D}\propto(1,1,1)$ , and the heavy Majorana neutrino mass term follows a similar form, as found in Ref.Feruglio:2017spp . By selecting another specific weights, namely $k_{L}=\frac{9}{2}$ , $k_{e^{c}}=k_{\mu^{c}}=k_{\tau^{c}}=-\frac{7}{2}$ , and $k_{N^{c}}=-\frac{3}{2}$ , a notable change occurs in the modular form of the Majorana neutrino operator. Specifically, the term $Y^{(2)}_{\bf 3}$ transforms into $Y^{(6)}_{\bf 1(3)}$ , resulting in an expression that aligns with the form presented in Eq.(53). While these cases show potential for reproducing lepton mass and mixing, further investigation is necessary to confirm its viability. (v) However, the assignment of the right-handed charged lepton $E^{c}=(e^{c},\mu^{c},\tau^{c})$ to the $A_{4}$ triplet, the left-handed charged lepton denoted as $L_{e,\mu,\tau}$ to $A_{4}$ singlets (1,1 ${}^{{}^{\prime\prime}}$ ,1 ${}^{{}^{\prime}}$ ), and $N^{c}$ to the $A_{4}$ triplet may not provide an explanation for the observed charged-lepton mass hierarchy. This difficulty arises from the charge assignment of $U(1)_{X}$ .

IV Quark and Lepton interactions with QCD axion

Let us discuss how quark and lepton masses and mixings are derived from Yukawa interactions within a framework based on $A_{4}\times U(1)_{X}$ symmetries with modular invariance. Non-zero VEVs of scalar fields $\chi(\tilde{\chi})$ spontaneously break the flavor symmetry $U(1)_{X}$ ¹²¹²12If the symmetry $U(1)_{X}$ is broken spontaneously, the massless mode $A_{X}$ of the scalar $\chi$ appear as a phase. at high energies above EW scale and create a heavy Majorana neutrino mass term. Then, the effective Yukawa structures in the low-energy limit depend on a small dimensionless parameter $\langle{\cal F}\rangle/\Lambda\equiv\Delta_{\cal F}$ . The higher order contributions of superpotentials $W^{(h)}_{q(l\nu)}$ become $\sum^{\infty}_{n=1}\tilde{c}_{i}\,\Delta^{2n}_{\chi}\cdot$ (leading order operators) with $\tilde{c}_{i}=e^{i\tilde{\theta}_{i}}$ , which make the Yukawa coefficients of the leading order terms in the superpotentials given in Eqs.(37, 50, 53) shifted. Denoting the effective Yukawa coefficients shifted by higher order contributions as $\alpha_{i},\beta_{i},\gamma_{i}$ , we see that they are constrained as

\displaystyle 1-\frac{\Delta^{2}_{\chi}}{1-\Delta^{2}_{\chi}}\leq|\alpha_{i}|,|\beta_{|},|\gamma_{i}|\leq 1+\frac{\Delta^{2}_{\chi}}{1-\Delta^{2}_{\chi}}\qquad\text{with}~{}\Delta_{\chi}\equiv\frac{v_{\chi}}{\sqrt{2}\,\Lambda},

(60)

where the lower (upper) limit corresponds to the sum of higher order terms with $\tilde{\theta}_{i}=\pi(0)$ . When $H_{u(d)}$ acquire non-zero VEVs, all quarks and leptons obtain masses. The relevant quark and lepton interactions with their chiral fermions are given by

	$\displaystyle-{\cal L}$	$\displaystyle\supset$	$\displaystyle\overline{q^{u}_{R}}\,\mathcal{M}_{u}\,q^{u}_{L}+\overline{q^{d}_{R}}\,\mathcal{M}_{d}\,q^{d}_{L}+\frac{g}{\sqrt{2}}W^{+}_{\mu}\overline{q^{u}_{L}}\gamma^{\mu}\,q^{d}_{L}$		(61)
		$\displaystyle+$	$\displaystyle\overline{\ell_{R}}\,{\cal M}_{\ell}\,\ell_{L}+\frac{1}{2}\begin{pmatrix}\overline{\nu^{c}_{L}}&\overline{N_{R}}\end{pmatrix}\begin{pmatrix}0&m^{T}_{D}\\ m_{D}&~{}~{}~{}e^{i\frac{A_{X}}{u_{\chi}}}\,M_{R}\end{pmatrix}\begin{pmatrix}\nu_{L}\\ N^{c}_{R}\end{pmatrix}+\frac{g}{\sqrt{2}}W^{-}_{\mu}\overline{\ell_{L}}\gamma^{\mu}\,\nu_{L}+\text{h.c.}\,,$		(61)

where $g$ is the $SU(2)_{L}$ coupling constant, $q^{u}=(u,c,t)$ , $q^{d}=(d,s,b)$ , $\ell=(e,\mu,\tau)$ , $\nu=(\nu_{e},\nu_{\mu},\nu_{\tau})$ , and $N=(N_{1},N_{2},N_{3})$ . $M_{R}$ contains a VEV of $\chi$ presented by Eq.(48). The explicit forms of $\mathcal{M}_{u,d,l}$ will be given later. The above Lagrangian of the fermions, including their kinetic terms of Eq.(46), should be invaruant under $U(1)_{X}$ :

\displaystyle\psi_{f}\rightarrow e^{if_{\psi_{f}}\frac{\gamma_{5}}{2}\alpha}\psi_{f}\,,\quad t=\text{invariant}\,,\quad N\rightarrow e^{i\frac{\gamma_{5}}{2}\alpha}N

(62)

where $\psi_{f}=\{u,c,d,s,b,e,\mu,\tau\}$ and $\alpha$ is a transformation constant parameter.

IV.1 Quark and flavored-QCD axion

As axion models, the axion-Yukawa coupling matrices and quark mass matrices in our model can be aresimultaneously diagonalized. The quark mass matrices are diagonalized through bi-unitary transformations: $V^{\psi}_{R}{\cal M}_{\psi}V^{\psi{\dagger}}_{L}=\hat{{\cal M}}_{\psi}$ (diagonal form) and the mass eigenstates are $\psi^{\prime}_{R}=V^{\psi}_{R}\,\psi_{R}$ and $\psi^{\prime}_{L}=V^{\psi}_{L}\,\psi_{L}$ . These transformation include, in particular, the chiral transformation of Eq.(62) that necessarily makes ${\cal M}_{u,d,\ell}$ real and positive. This induces a contribution to the QCD vacuum angle in Eq.(46), i.e.,

\displaystyle\vartheta_{\rm QCD}\rightarrow\vartheta_{\rm eff}=\vartheta_{\rm QCD}+\arg\{\det({\cal M}_{u})\det({\cal M}_{d})\}

(63)

with $-\pi\leq\vartheta_{\rm eff}\leq\pi$ . Then one obtains the vanishing QCD anomaly term

\displaystyle{\cal L}_{\vartheta}=\Big{(}\vartheta_{\rm eff}+\frac{A_{X}}{F_{a}}\Big{)}\frac{\alpha_{s}^{\prime}}{8\pi}G^{a\mu\nu}\tilde{G}^{a}_{\mu\nu}\qquad\text{with}~{}F_{a}=\frac{f_{a}}{N_{C}}\,,

(64)

where $\alpha_{s}^{\prime}=g^{2}_{s}/4\pi$ and the axion decay constant $F_{a}$ with $f_{a}=u_{\chi}$ of Eq.(48). At low energies $A_{X}$ will get a VEV, $\langle A_{X}\rangle=-F_{a}\vartheta_{\rm eff}$ , eliminating the constant $\vartheta_{\rm eff}$ term. The QCD axion then is the excitation of the $A_{X}$ field, $a=A_{X}-\langle A_{X}\rangle$ .

Substituting the VEV of Eq.(40) into the superpotential Eq.(50), the mass matrices ${\cal M}_{u}$ and ${\cal M}_{d}$ for up- and down-type quarks given in the Lagrangian (61) are derived as,

	$\displaystyle{\cal M}_{u}={\left(\begin{array}[]{ccc}\alpha_{u}\Delta^{\|f_{u}\|}_{\chi}\,Y^{(6)}_{\bf 1}e^{if_{u}\frac{A_{X}}{u_{\chi}}}&0&0\\ 0&\alpha_{c}\Delta^{\|f_{c}\|}_{\chi}\,Y^{(6)}_{\bf 1}e^{if_{c}\frac{A_{X}}{u_{\chi}}}&0\\ 0&0&\alpha_{t}\end{array}\right)}\,v_{u}\,,$		(68)
	$\displaystyle{\cal M}_{d}=\Big{[}{\left(\begin{array}[]{ccc}\alpha_{d}\Delta^{\|f_{d}\|}_{\chi}&\alpha_{s}\,x\Delta^{\|f_{s}\|}_{\chi}&\alpha_{b}\,y\Delta^{\|f_{b}\|}_{\chi}\\ \alpha_{d}\,y\Delta^{\|f_{d}\|}_{\chi}&\alpha_{s}\Delta^{\|f_{s}\|}_{\chi}&\alpha_{b}\,x\Delta^{\|f_{b}\|}_{\chi}\\ \alpha_{d}\,x\Delta^{\|f_{d}\|}_{\chi}&\alpha_{s}\,y\Delta^{\|f_{s}\|}_{\chi}&\alpha_{b}\Delta^{\|f_{b}\|}_{\chi}\end{array}\right)}(1+2xy)$		(72)
	$\displaystyle\qquad\qquad\qquad\qquad+{\left(\begin{array}[]{ccc}\tilde{\alpha}_{d}\,y\Delta^{\|f_{d}\|}_{\chi}&\tilde{\alpha}_{s}\Delta^{\|f_{s}\|}_{\chi}&\tilde{\alpha}_{b}\,x\Delta^{\|f_{b}\|}_{\chi}\\ \tilde{\alpha}_{d}\,x\Delta^{\|f_{d}\|}_{\chi}&\tilde{\alpha}_{s}\,y\Delta^{\|f_{s}\|}_{\chi}&\tilde{\alpha}_{b}\Delta^{\|f_{b}\|}_{\chi}\\ \tilde{\alpha}_{d}\Delta^{\|f_{d}\|}_{\chi}&\tilde{\alpha}_{s}\,x\Delta^{\|f_{s}\|}_{\chi}&\tilde{\alpha}_{b}\,y\Delta^{\|f_{b}\|}_{\chi}\end{array}\right)}(y^{2}+2x)\Big{]}\tilde{C}\,Y^{3}_{1}\,v_{d}\,,$		(76)

where $v_{d}\equiv\langle H_{d}\rangle=v\cos\beta/\sqrt{2}$ , $v_{u}\equiv\langle H_{u}\rangle=v\sin\beta/\sqrt{2}$ with $v\simeq 246$ GeV, and

\displaystyle\tilde{C}={\rm diag.}\big{(}e^{if_{d}\frac{A_{X}}{u_{\chi}}},e^{if_{s}\frac{A_{X}}{u_{\chi}}},e^{if_{b}\frac{A_{X}}{u_{\chi}}}\big{)}\,,\quad x=\frac{Y_{2}}{Y_{1}}\,,\quad y=\frac{Y_{3}}{Y_{1}}\,.

(77)

The terms with $\alpha_{d,s,b}$ in Eq.(76) generate by taking the modular form $Y^{(6)}_{{\bf 3},1}$ given in Eq.(51), whereas the terms with $\tilde{\alpha}_{d,s,b}$ in Eq.(76) generate by taking $Y^{(6)}_{{\bf 3},2}$ .

The quark mass matrices ${\cal M}_{u}$ in Eq.(68) and ${\cal M}_{d}$ in Eq.(76) generate the up- and down-type quark masses:

\displaystyle\hat{\mathcal{M}}_{u}=V^{u}_{R}\,{\cal M}_{u}\,V^{u{\dagger}}_{L}={\rm diag}(m_{u},m_{c},m_{t})\,,\quad\hat{\mathcal{M}}_{d}=V^{d}_{R}\,{\cal M}_{d}\,V^{d{\dagger}}_{L}={\rm diag}(m_{d},m_{s},m_{b})\,.

(78)

Diagonalizing the matrices ${\cal M}_{f}^{\dagger}{\cal M}_{f}$ and ${\cal M}_{f}{\cal M}_{f}^{\dagger}$ ( $f=u,d$ ) determine the mixing matrices $V_{L}^{f}$ and $V_{R}^{f}$ , respectively Ahn:2011yj . The left-handed quark mixing matrices $V_{L}^{u}$ and $V_{L}^{d}$ are components of the CKM matrix $V_{\rm CKM}=V_{L}^{u}V_{L}^{d\dagger}$ , which is generated from the down-type quark matrix in Eq.(76) due to the diagonal form of the up-type quark mass matrix in Eq.(68). The CKM matrix is parameterized by the Wolfenstein parametrization Wolfenstein:1983yz , see Eq.(139), and has been determined with high precision Ahn:2011fg . The current best-fit values of the CKM mixing angles in the standard parameterization Chau:1984fp read in the $3\sigma$ range ckm

\displaystyle\theta^{q}_{23}[^{\circ}]=2.376^{+0.054}_{-0.070}\,,\quad\theta^{q}_{13}[^{\circ}]=0.210^{+0.016}_{-0.010}\,,\quad\theta^{q}_{12}[^{\circ}]=13.003^{+0.048}_{-0.036}\,,\quad\delta^{q}_{CP}[^{\circ}]=65.5^{+3.1}_{-4.9}\,.

(79)

The physical structure of the up- and down-type quark Lagrangian should match up with the empirical results calculated from the Particle Data Group (PDG) PDG :

	$\displaystyle m_{d}$	$\displaystyle=$	$\displaystyle 4.67^{+0.48}_{-0.17}\,{\rm MeV}\,,\qquad~{}m_{s}=93^{+11}_{-5}\,{\rm MeV}\,,\qquad\qquad~{}~{}m_{b}=4.18^{+0.03}_{-0.02}\,{\rm GeV}\,,$
	$\displaystyle m_{u}$	$\displaystyle=$	$\displaystyle 2.16^{+0.49}_{-0.29}\,{\rm MeV}\,,\qquad~{}m_{c}=1.27\pm 0.02\,{\rm GeV}\,,\qquad~{}m_{t}=173.1\pm 0.9\,{\rm GeV}\,,$		(80)

where $t$ -quark mass is the pole mass, $c$ - and $b$ -quark masses are the running masses in the $\overline{\rm MS}$ scheme, and the light $u$ -, $d$ -, $s$ -quark masses are the current quark masses in the $\overline{\rm MS}$ scheme at the momentum scale $\mu\approx 2$ GeV. Below the scale of spontaneous $SU(2)_{L}\times U(1)_{Y}$ gauge symmetry breaking, the running masses of $c$ - and $b$ -quark receive corrections from QCD and QED loops PDG . The top quark mass at scales below the pole mass is unphysical since the $t$ -quark decouples at its scale, and its mass is determined more directly by experiments PDG .

After diagonalizing the mass matrices of Eqs.(68, 76), the flavored-QCD axion to quark interactions are written at leading order as

$\displaystyle-{\cal L}^{aq}$	$\displaystyle\simeq$	$\displaystyle\frac{\partial_{\mu}a}{2u_{\chi}}\Big{\{}f_{u}\,\bar{u}\gamma^{\mu}\gamma_{5}u+f_{c}\,\bar{c}\gamma^{\mu}\gamma_{5}c+f_{d}\,\bar{d}\gamma^{\mu}\gamma_{5}d+f_{s}\,\bar{s}\gamma^{\mu}\gamma_{5}s+f_{b}\,\bar{b}\gamma^{\mu}\gamma_{5}b\Big{\}}$	(81)
	$\displaystyle+$	$\displaystyle\frac{\partial_{\mu}a}{2u_{\chi}}\Big{\{}(f_{d}-f_{s})\lambda\Big{(}1-\frac{\lambda^{2}}{2}\Big{)}\bar{d}\gamma^{\mu}\gamma_{5}s+(f_{s}-f_{b})A_{d}\lambda^{2}\,\bar{s}\gamma^{\mu}\gamma_{5}b$
		$\displaystyle\qquad\qquad+A_{d}\lambda^{3}\Big{(}f_{d}(\rho+i\eta)-f_{s}+f_{b}(1-\rho-i\eta)\Big{)}\bar{b}\gamma^{\mu}\gamma_{5}d+{\rm h.c.}\Big{\}}$
	$\displaystyle+$	$\displaystyle i\frac{a}{2u_{\chi}}\Big{\{}(f_{d}-f_{s})\lambda\Big{(}1-\frac{\lambda^{2}}{2}\Big{)}(m_{d}-m_{s})\,\bar{d}s+(f_{s}-f_{b})A_{d}\lambda^{2}(m_{s}-m_{b})\,\bar{s}b$
		$\displaystyle\qquad\qquad+A_{d}\lambda^{3}\Big{(}f_{d}(\rho+i\eta)-f_{s}+f_{b}(1-\rho-i\eta)\Big{)}(m_{b}-m_{d})\,\bar{b}d+{\rm h.c.}\Big{\}}$
	$\displaystyle+$	$\displaystyle m_{u}\,\bar{u}u+m_{c}\,\bar{c}c+m_{t}\,\bar{t}t+m_{d}\,\bar{d}d+m_{s}\,\bar{s}s+m_{b}\,\bar{b}b-\bar{q}i\!\!\not\!D\,q\,,$

where $V^{d{\dagger}}_{L}=V_{\rm CKM}$ of Eq.(139) is used by rotating the phases in ${\cal M}_{u}$ away, which is the result of a direct interaction of the SM gauge singlet scalar field $\chi$ with the SM quarks charged under $U(1)_{X}$ . The flavored-QCD axion $a$ is produced by flavor-changing neutral Yukawa interactions in Eq.(81), which leads to induced rare flavor-changing processes. The strongest bound on the QCD axion decay constant is from the flavor-changing process $K^{+}\rightarrow\pi^{+}+a$ Wilczek:1982rv ; Bolton:1988af ; Artamonov:2008qb ; raredecay ; Berezhiani:1989fp , induced by the flavored-QCD axion $a$ . From Eq.(81) the flavored-QCD axion interactions with the flavor violating coupling to the $s$ - and $d$ -quark are given by

\displaystyle-{\cal L}^{asd}_{Y}\simeq\frac{i}{2}(f_{d}-f_{s})\frac{a}{N_{C}F_{a}}\bar{s}d\,(m_{d}-m_{s})\lambda\Big{(}1-\frac{\lambda^{2}}{2}\Big{)}\,.

(82)

Then the decay width of $K^{+}\rightarrow\pi^{+}+a$ is given by

\displaystyle\Gamma(K^{+}\rightarrow\pi^{+}+a)=\frac{m^{3}_{K}}{16\pi}\Big{(}1-\frac{m^{2}_{\pi}}{m^{2}_{K}}\Big{)}^{3}\Big{|}\frac{f_{d}-f_{s}}{2\,F_{a}N_{C}}\lambda\Big{(}1-\frac{\lambda^{2}}{2}\Big{)}\Big{|}^{2}\,,

(83)

where $m_{K^{\pm}}=493.677\pm 0.013$ MeV, $m_{\pi^{\pm}}=139.57061\pm 0.00024$ MeV PDG . From the present experimental upper bound ${\rm Br}(K^{+}\rightarrow\pi^{+}a)<(3-6)\times 10^{-11}(1\times 10^{-11})$ for $m_{a}=0-110$ (160-260) MeV at $90\%$ CL with ${\rm Br}(K^{+}\rightarrow\pi^{+}\nu\bar{\nu})=(10.6^{+4.0}_{-3.4}|_{\rm stat}\pm 0.9_{\rm syst})\times 10^{-11}$ at $68\%$ CL NA62:2021zjw , we obtain the lower limit on the QCD axion decay constant,

\displaystyle F_{a}\frac{2|N_{C}|}{|f_{d}-f_{s}|}\gtrsim(0.86-1.90)\times 10^{11}\,{\rm GeV}\,.

(84)

The QCD axion mass $m_{a}$ in terms of the pion mass and pion decay constant reads Ahn:2014gva ; Ahn:2016hbn

\displaystyle m^{2}_{a}F^{2}_{a}=m^{2}_{\pi^{0}}f^{2}_{\pi}F(z,w)\,,

(85)

where $f_{\pi}\simeq 92.1$ MeV PDG and $F(z,w)=z/(1+z)(1+z+w)$ with $\omega=0.315\,z$ . Here the Weinberg value lies in $z\equiv m^{\overline{\rm MS}}_{u}(2\,{\rm GeV})/m^{\overline{\rm MS}}_{d}(2\,{\rm GeV})=0.47^{+0.06}_{-0.07}$ PDG . After integrating out the heavy $\pi^{0}$ and $\eta$ at low energies, there is an effective low energy Lagrangian with an axion-photon coupling $g_{a\gamma\gamma}$ : ${\cal L}_{a\gamma\gamma}=-g_{a\gamma\gamma}\,a\,\vec{E}\cdot\vec{B}$ where $\vec{E}$ and $\vec{B}$ are the electromagnetic field components. The axion-photon coupling is expressed in terms of the QCD axion mass, pion mass, pion decay constant, $z$ and $w$ ,

\displaystyle g_{a\gamma\gamma}=\frac{\alpha_{\rm em}}{2\pi}\frac{m_{a}}{f_{\pi}m_{\pi^{0}}}\frac{1}{\sqrt{F(z,w)}}\left(\frac{E}{N_{C}}-\frac{2}{3}\,\frac{4+z+w}{1+z+w}\right)\,.

(86)

The upper bound on the axion-photon coupling, derived from the recent analysis of the horizontal branch stars in galactic globular clusters Ayala:2014pea , can be translated to

\displaystyle|g_{a\gamma\gamma}|<6.6\times 10^{-11}\,{\rm GeV}^{-1}\,(95\%\,{\rm CL})\Leftrightarrow F_{a}\gtrsim 2.525\times 10^{7}\Big{|}\frac{E}{N_{C}}-1.903\Big{|}\,{\rm GeV}\,,

(87)

where $z=0.47$ is used.

IV.2 Charged-lepton and flavored-QCD axion

Substituting the VEV ofEq.(40) into the superpotential (53), the charged-lepton mass matrix given in the Lagrangian (61) is derived as,

\displaystyle{\cal M}_{\ell}

\displaystyle=

\displaystyle{\left(\begin{array}[]{ccc}\alpha_{e}\Delta^{|f_{e}|}_{\chi}\,e^{if_{e}\frac{A_{X}}{u_{\chi}}}&0&0\\ 0&\alpha_{\mu}\,\Delta^{|f_{\mu}|}_{\chi}\,e^{if_{\mu}\frac{A_{X}}{u_{\chi}}}&0\\ 0&0&\alpha_{\tau}\,\Delta^{|f_{\tau}|}_{\chi}\,e^{if_{\tau}\frac{A_{X}}{u_{\chi}}}\end{array}\right)}Y^{(6)}_{\bf 1}\,v_{d}\,.

(91)

Recall that the coefficients $\alpha_{i}$ are complex numbers with an effective absolute value satisfying Eq.(60). Then, the corresponding charged-lepton masses are given by

\displaystyle m_{e}

\displaystyle=

\displaystyle\alpha_{e}\,Y^{(6)}_{\bf 1}\Delta^{|f_{e}|}_{\chi}\,v_{d}\,,\quad m_{\mu}=\alpha_{\mu}\,Y^{(6)}_{\bf 1}\Delta^{|f_{\mu}|}_{\chi}\,v_{d}\,,\quad m_{\tau}=\alpha_{\tau}\,Y^{(6)}_{\bf 1}\Delta^{|f_{\tau}|}_{\chi}\,v_{d}\,,

(92)

where $Y^{(6)}_{\bf 1}$ is given in Eq.(51) and the phases in each term can be absorbed into $(l_{i})_{R}$ . They are matched with the empirical values from the PDG PDG given by

\displaystyle m_{e}

\displaystyle=

\displaystyle 0.511\,{\rm MeV}\,,\qquad m_{\mu}=105.658\,{\rm MeV}\,,\qquad m_{\tau}=1776.86\pm 0.12\,{\rm MeV}\,.

(93)

Flavored axions typically interact with charged leptons (electrons, muons, taus) Ahn:2014gva ; Ahn:2016hbn ; Ahn:2018cau ; Ahn:2021ndu and can be emitted through atomic axio-recombination, axio-deexcitation, axio-bremsstrahlung in electron-ion or electron-electron collisions, and Compton scatterings Redondo:2013wwa . Then the flavored-QCD axion to charged-lepton interactions read

\displaystyle-{\cal L}^{a\ell}

\displaystyle\simeq

\displaystyle\frac{\partial_{\mu}a}{2u_{\chi}}\Big{(}f_{e}\,\bar{e}\gamma^{\mu}\gamma_{5}e+f_{\mu}\,\bar{\mu}\gamma^{\mu}\gamma_{5}\mu+f_{\tau}\,\bar{\tau}\gamma^{\mu}\gamma_{5}\tau\Big{)}+\sum_{\ell=e,\mu,\tau}\big{(}m_{\ell}\bar{\ell}\ell-\bar{\ell}i\!\!\not\!\partial\,\ell\big{)}\,.

(94)

Like rare neutral flavor-changing decays in particle physics, the interaction of the flavored-QCD axion $a$ with leptons makes it possible to search for the QCD axion in astroparticle physics through stellar evolution. The flavored-QCD axion coupling to electrons reads

\displaystyle g_{aee}

\displaystyle=

\displaystyle|f_{e}|\,\frac{m_{e}}{u_{\chi}}\,.

(95)

Stars in the red giant branch (RGB) of color-magnitude diagrams in globular clusters provide a strict constraint on axion-electron couplings which leads to a lower bound on the axion decay constant. This constraint is expressed as Viaux:2013lha

\displaystyle|g_{aee}|<4.3\times 10^{-13}\quad(95\%\,{\rm CL})~{}\Leftrightarrow N_{C}F_{a}\gtrsim 1.19|f_{e}|\times 10^{9}\,{\rm GeV}\,.

(96)

Bremsstrahlung off electrons $e+Ze\rightarrow Ze+e+a$ in white dwarfs is an effective process for detecting axions as the Primakoff and Compton processes are suppressed due to the large plasma frequency. Comparing the theoretical and observed WD luminosity functions (WDLFs) provides a way to place limits ¹³¹³13Note that Refs. WD01 ; Bertolami:2014wua have pointed out features in some WDLFs DeGennaro:2007yw ; Rowell:2011wp that could imply axion-electron couplings in the range $7.2\times 10^{-14}\lesssim|g_{aee}|\lesssim 2.2\times 10^{-13}$ . on $|g_{aee}|$ Raffelt:1985nj . Recent analyses of WDLFs, using detailed WD cooling treatment and new data on the WDLF of the Galactic disk, suggest electron couplings $|g_{aee}|\lesssim 2.8\times 10^{-13}$ Bertolami:2014wua . However, these results come with large theoretical and observational uncertainties.

We note that the entries of the quark and charged lepton mass matrices given in Eqs.(68), (76), and (91) except for the entry corresponding to the top quark are expressed as a combination of $\Delta_{\chi}^{|f_{\alpha}|}$ and modular forms for each component. Accurate determination of the values of $\Delta_{\chi}$ , its power, and the value of $\tau$ is crucial to reproduce the observed CKM mixing angles given in Eq.(79) and quark masses in Eq.(80). The values of those parameters are also closely linked to those in the lepton sector, and they should necessarily be determined in order to reproduce the observed values of charged lepton masses and to predict light active neutrinos derived from Eqs.(104) and (112). The $U(1)_{X}$ PQ scale, which corresponds to the seesaw scale (as shown n Eq.(105)), can be estimated as $\mu_{\chi}\sim 6\times 10^{13}$ GeV from Eq.(45) for $m_{3/2}\sim 100$ TeV.

IV.3 Neutrino

Similar to the case of charged lepton mass matrix, the heavy Majorana mass matrix given in the Lagrangian (61) is derived from the superpotential (53) as

	$\displaystyle M_{R}$	$\displaystyle=$	$\displaystyle M{\left(\begin{array}[]{ccc}1+y^{3}-5xy&0&0\\ 0&0&1+y^{3}-5xy\\ 0&1+y^{3}-5xy&0\end{array}\right)}$		(100)
		$\displaystyle+$	$\displaystyle M{\left(\begin{array}[]{ccc}\frac{2}{\sqrt{3}}(\gamma\,p+\gamma^{\prime}\,yr)&-\frac{1}{\sqrt{3}}(\gamma\,yp+\gamma^{\prime}\,xr)&-\frac{1}{\sqrt{3}}(\gamma\,xp+\gamma^{\prime}\,r)\\ -\frac{1}{\sqrt{3}}(\gamma\,yp+\gamma^{\prime}\,xr)&\frac{2}{\sqrt{3}}(\gamma\,xp+\gamma^{\prime}\,r)&-\frac{1}{\sqrt{3}}(\gamma\,p+\gamma^{\prime}\,yr)\\ -\frac{1}{\sqrt{3}}(\gamma\,xp+\gamma^{\prime}\,r)&-\frac{1}{\sqrt{3}}(\gamma\,p+\gamma^{\prime}\,yr)&\frac{2}{\sqrt{3}}(\gamma\,yp+\gamma^{\prime}\,xr)\end{array}\right)}\,,$		(104)

where $p=1+2xy$ , $r=y^{2}+2xy$ , $\gamma=\gamma_{2}/\gamma_{1}$ , $\gamma^{\prime}=\gamma^{\prime}_{2}/\gamma_{1}$ , and the common factor $M$ can be replaced by the QCD axion decay constant $F_{a}$ ,

\displaystyle M\equiv|\gamma_{1}Y^{3}_{1}\langle\chi\rangle|=\frac{|\gamma_{1}|}{2}F_{a}|N_{C}Y^{3}_{1}|\,.

(105)

The terms with $\gamma_{2}$ in Eq.(104) are derived by taking the modular form $Y^{(6)}_{{\bf 3},1}$ safisfying Eq.(51), whereas the terms with $\gamma^{\prime}_{2}$ are derived by taking $Y^{(6)}_{{\bf 3},2}$ . Eq.(104) has three unknown complex parameters, $\gamma$ $\gamma^{\prime}$ , and $\gamma_{1}$ , where the phase of $\gamma_{1}$ contributes as an overall factor after seesawing. Other variables such as $x,y,Y_{1}$ are determined from the analysis for the quark and charged-lepton sectors, and $\langle\chi\rangle$ is fixed from the seesaw formula Eq.(113) whose scale is given by PQ scale Eq.(45). The Dirac mass term in the Lagrangian (61) reads

\displaystyle m_{D}={\left(\begin{array}[]{ccc}\beta_{1}\,\Delta^{|g_{e}|}_{\chi}&\beta_{2}\,y\Delta^{|g_{\mu}|}_{\chi}&\beta_{3}\,x\Delta^{|g_{\tau}|}_{\chi}\\ \beta_{1}\,y\Delta^{|g_{e}|}_{\chi}&\beta_{2}\,x\Delta^{|g_{\mu}|}_{\chi}&\beta_{3}\,\Delta^{|g_{\tau}|}_{\chi}\\ \beta_{1}\,x\Delta^{|g_{e}|}_{\chi}&\beta_{2}\,\Delta^{|g_{\mu}|}_{\chi}&\beta_{3}\,y\Delta^{|g_{\tau}|}_{\chi}\end{array}\right)}{\left(\begin{array}[]{ccc}e^{ig_{e}\frac{A_{X}}{u_{\chi}}}&0&0\\ 0&e^{ig_{\mu}\frac{A_{X}}{u_{\chi}}}&0\\ 0&0&e^{ig_{\tau}\frac{A_{X}}{u_{\chi}}}\end{array}\right)}\,Y_{1}\,v_{u}\,.

(112)

The coefficients $\beta_{i}$ and $\gamma_{i}$ in the neutrino sector, like in the quark and charged-lepton sectors, are complex numbers corrected by higher-dimensional operators, resulting in an effective absolute value satisfying Eq.(60). Eq.(112) contains three complex parameters ( $\beta_{1},\beta_{2}$ , and $\beta_{3}$ ), where one of the phases can be removable as an overall factor after seesawing. As shwon before, the parameter $\Delta_{\chi}$ can be determined from quark and charged-lepton sectors. In addition, its $U(1)_{X}$ quantum number $g_{\alpha}$ can be determined from the numerical analysis for the neutrino sector with the help of the condition of $U(1)_{X}$ -mixed gravitational anomaly-free given in in Eq.(54).

After integrating out the right-handed heavy Majorana neutrinos, the effective neutrino mass matrix ${\cal M}_{\nu}$ is given at leading order by

\displaystyle{\cal M}_{\nu}\simeq-m^{T}_{D}M^{-1}_{R}m_{D}=U^{\ast}_{\nu}{\rm diag.}(m_{\nu_{1}},m_{\nu_{2}},m_{\nu_{3}})U^{{\dagger}}_{\nu}\,,

(113)

where $U_{\nu}$ is the rotation matrix diagonalizing ${\cal M}_{\nu}$ and $m_{\nu_{i}}$ ( $i=1,2,3$ ) are the light neutrino masses. Then the PMNS mixing matrix becomes

\displaystyle U_{\rm PMNS}=U_{\nu}.

(114)

The matrix $U_{\rm PMNS}$ is expressed in terms of three mixing angles, $\theta_{12},\theta_{13},\theta_{23}$ , and a Dirac type CP violaitng phase $\delta_{CP}$ and two additional CP violating phases $\varphi_{1,2}$ if light neutrinos are Majorana particle as PDG

\displaystyle U_{\rm PMNS}={\left(\begin{array}[]{ccc}c_{13}c_{12}&c_{13}s_{12}&s_{13}e^{-i\delta_{CP}}\\ -c_{23}s_{12}-s_{23}c_{12}s_{13}e^{i\delta_{CP}}&c_{23}c_{12}-s_{23}s_{12}s_{13}e^{i\delta_{CP}}&s_{23}c_{13}\\ s_{23}s_{12}-c_{23}c_{12}s_{13}e^{i\delta_{CP}}&-s_{23}c_{12}-c_{23}s_{12}s_{13}e^{i\delta_{CP}}&c_{23}c_{13}\end{array}\right)}Q_{\nu}\,,

(118)

where $s_{ij}\equiv\sin\theta_{ij}$ , $c_{ij}\equiv\cos\theta_{ij}$ and $Q_{\nu}={\rm Diag}(e^{-i\varphi_{1}/2},e^{-i\varphi_{2}/2},1)$ . Then the neutrino masses are obtained by the transformation

\displaystyle U^{T}_{\rm PMNS}\,{\cal M}_{\nu}\,U_{\rm PMNS}={\rm Diag}.(m_{\nu_{1}},m_{\nu_{2}},m_{\nu_{3}})\,.

(119)

Here $m_{\nu_{i}}$ ( $i=1,2,3$ ) are the light neutrino masses. The observed hierarchy $|\Delta m^{2}_{\rm Atm}|=|m^{2}_{\nu_{3}}-(m^{2}_{\nu_{1}}+m^{2}_{\nu_{2}})/2|\gg\Delta m^{2}_{\rm Sol}\equiv m^{2}_{\nu_{2}}-m^{2}_{\nu_{1}}>0$ and the requirement of a Mikheyev-Smirnov-Wolfenstein resonance Wolfenstein:1977ue for solar neutrinos lead to two possible neutrino mass spectra: normal mass ordering (NO) $m^{2}_{\nu_{1}}<m^{2}_{\nu_{2}}<m^{2}_{\nu_{3}}$ and inverted mass ordering (IO) $m^{2}_{\nu_{3}}<m^{2}_{\nu_{1}}<m^{2}_{\nu_{2}}$ . Nine physical observables can be derived from Eqs.(118) and (119): $\theta_{23}$ , $\theta_{13}$ , $\theta_{12}$ , $\delta_{CP}$ , $\varphi_{1}$ , $\varphi_{2}$ , $m_{\nu_{1}}$ , $m_{\nu_{2}}$ , and $m_{\nu_{3}}$ .

Table 3: The global fit of three-flavor oscillation parameters at the best-fit and

3\sigma

level with Super-Kamiokande atmospheric data Esteban:2020cvm . NO = normal neutrino mass ordering; IO = inverted mass ordering. And

\Delta m^{2}_{\rm Sol}\equiv m^{2}_{\nu_{2}}-m^{2}_{\nu_{1}}

\Delta m^{2}_{\rm Atm}\equiv m^{2}_{\nu_{3}}-m^{2}_{\nu_{1}}

for NO, and

\Delta m^{2}_{\rm Atm}\equiv m^{2}_{\nu_{2}}-m^{2}_{\nu_{3}}

for IO.

	$\theta_{13}[^{\circ}]$	$\delta_{CP}[^{\circ}]$	$\theta_{12}[^{\circ}]$	$\theta_{23}[^{\circ}]$	$\Delta m^{2}_{\rm Sol}[10^{-5}{\rm eV}^{2}]$	$\Delta m^{2}_{\rm Atm}[10^{-3}{\rm eV}^{2}]$
$\begin{array}[]{ll}\hbox{NO}\\ \hbox{IO}\end{array}$	$\begin{array}[]{ll}8.58^{+0.33}_{-0.35}\\ 8.57^{+0.37}_{-0.34}\end{array}$	$\begin{array}[]{ll}232^{+118}_{-88}\\ 276^{+68}_{-82}\end{array}$	$33.41^{+2.33}_{-2.10}$	$\begin{array}[]{ll}42.2^{+8.8}_{-2.5}\\ 49.0^{+2.5}_{-9.1}\end{array}$	$7.41^{+0.62}_{-0.59}$	$\begin{array}[]{ll}2.507^{+0.083}_{-0.080}\\ 2.486^{+0.084}_{-0.080}\end{array}$

Recent global fits Esteban:2018azc ; deSalas:2017kay ; Capozzi:2018ubv of neutrino oscillations have enabled a more precise determination of the mixing angles and mass squared differences, with large uncertainties remaining for $\theta_{23}$ and $\delta_{CP}$ at 3 $\sigma$ . The most recent analysis Esteban:2020cvm lists global fit values and $3\sigma$ intervals for these parameters in Table-3. Furthermore, recent constraints on the rate of $0\nu\beta\beta$ decay have added to these findings. Specifically, the most tight upper bounds for the effective Majorana mass ( ${\cal M}_{\nu})_{ee}$ , which is the modulus of the $ee$ -entry of the effective neutrino mass matrix, are given by

\displaystyle({\cal M}_{\nu})_{ee}<0.036-0.156\,{\rm eV}\,~{}(^{136}\text{Xe-based experiment\,\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{KamLAND-Zen:2022tow}{\@@citephrase{(}}{\@@citephrase{)}}}})

(120)

at $90\%$ confidence level. $0\nu\beta\beta$ decay is a low-energy probe of lepton number violation and its measurement could provide the strongest evidence for lepton number violation at high energy. Its discovery would suggest the Majorana nature of neutrinos and, consequently, the existence of heavy Majorana neutrinos via the seesaw mechanism Minkowski:1977sc .

Transforming the neutrino fields by chiral rotations of Eq.(62) under $U(1)_{X}$ we obtain the flavored-QCD axion interactions to neutrinos

\displaystyle-{\cal L}^{a\nu}

\displaystyle\simeq

\displaystyle\frac{a}{2u_{\chi}}\sum_{i\neq j}\overline{\nu_{i}}\big{\{}(m_{\nu_{i}}-m_{\nu_{j}}){\rm Im}[V]_{ij}-i\gamma_{5}(m_{\nu_{i}}+m_{\nu_{j}}){\rm Re}[V]_{ij}\big{\}}\nu_{j}+\frac{\partial a}{4u_{\chi}}\overline{N}\gamma_{\mu}\gamma_{5}N^{c}

(121)

where $m_{\nu_{i}}$ are real and positive, ${\cal Q}[V]_{ij}=(g_{e}-\frac{1}{2}){\cal Q}[U_{1i}U^{\ast}_{1j}]+(g_{\mu}-\frac{1}{2}){\cal Q}[U_{2i}U^{\ast}_{2j}]+(g_{\tau}-\frac{1}{2}){\cal Q}[U_{3i}U^{\ast}_{3j}]$ with ${\cal Q}={\rm Re}$ or ${\rm Im}$ , and ${\rm Im}[V]_{ij}=-{\rm Im}[V]_{ji}$ with $U\equiv U_{\rm PMNS}$ . Since the light neutrino mass is less than 0.1 eV, the coupling between the flavored-QCD axion and light neutrinos is subject to a stringent constraint given by Eq.(84), which significantly suppresses the interaction. Therefore, we will not take it into consideration. Ref.Kharusi:2021jez provides the latest experimental constraints on Majoron-neutrino coupling, which are below the range of $(0.4-0.9)\times 10^{-5}$ .

Table 4:

U(1)_{X}

charges linked to seesaw scale.

	$\|g_{e}\|$	$\|g_{\mu}\|$	$\|g_{\tau}\|$	$\langle\chi\rangle$ /GeV
NO	$6$	$4$	$5$	$10^{13}$
	$5$	$3$	$4$	$5\times 10^{13}$
	$4$	$2$	$3$	$10^{14}$
	$3$	$1$	$2$	$5\times 10^{14}$
	$2$	$0$	$1$	$10^{15}$
IO	$2$	$3$	3	$5\times 10^{13}$
	$1$	2	2	$10^{14}$
	$0$	$1$	1	$5\times 10^{14}$

Once the lepton $U(1)_{X}$ quantum numbers are fixed, the seesaw scale $M\sim v_{\chi}$ of Eq.(105) comparable to the PQ scale of Eq.(45) can be roughly determined using the seesaw formula Eq.(113). By putting Eqs.(104, 112) into the seesaw formula Eq.(113), we obtain numerically a range of values for $\langle\chi\rangle$ . For instance, see Table-4: it implies that for normal mass ordering, the maximum scale should be below $\sim 10^{15}$ GeV, and for inverted mass ordering, the maximum scale should be below $\sim 5\times 10^{14}$ GeV. Refer to Table-5 and -6 for NO, and Table-7 for IO, with $F_{a}=2|\langle\chi\rangle/N_{C}|$ . By using the seesaw formula Eq.(113), one can set the scale $\langle\chi\rangle$ along with the $U(1)_{X}$ quantum numbers $g_{\alpha}$ without loss of generality. Doing so results in an effective mass matrix with nine physical degrees of freedom

\displaystyle m_{0}\equiv\Big{|}\frac{\beta^{2}_{1}}{\gamma_{1}Y_{\bf 1}}\Big{|}\frac{v^{2}_{u}}{\langle\chi\rangle}\,,~{}|\gamma|\,,~{}|\gamma^{\prime}|\,,~{}|\tilde{\beta}_{2}|\,,~{}|\tilde{\beta}_{3}|\,,~{}\arg(\gamma)\,,~{}\arg(\gamma^{\prime})\,,~{}\arg(\tilde{\beta}_{2})\,,~{}\arg(\tilde{\beta}_{3})\,,

(122)

in which $m_{0}$ is an overall factor of Eq.(113), $\tilde{\beta}_{2}\equiv\beta_{2}/\beta_{1},\tilde{\beta}_{3}\equiv\beta_{3}/\beta_{1}$ , and $1-2\Delta^{2}_{\chi}\leq|\tilde{\beta}_{2(3)}|,|\gamma|,|\gamma^{\prime}|\leq\frac{1}{1-2\Delta^{2}_{\chi}}$ . Out of the nine observables corresponding to Eq.(122), the five measured quantities ( $\theta_{12}$ , $\theta_{23}$ , $\theta_{13}$ , $\Delta m^{2}_{\rm Sol}$ , $\Delta m^{2}_{\rm Atm}$ ) can be used as constraints. The remaining four degrees of freedom correspond to four measurable quantities ( $\delta_{CP}$ , $\varphi_{1,2}$ , and the $0\nu\beta\beta$ -decay rate), which can be determined through measurements.

V Numerical analysis for quark, lepton, and a QCD axion

To simulate and match experimental results for quarks and leptons (Eqs.(79, 80) and Table-3), we use linear algebra tools from Ref.Antusch:2005gp . By analyzing experimental data for quarks and charged leptons, we determined the $U(1)_{X}$ quantum numbers listed in Tables-5 and -6 for normal neutrino mass ordering and Table-7 for inverted neutrino mass ordering. We also ensured the $U(1)_{X}$ -mixed gravitational anomaly-free condition of Eq.(54) and consistency of the seesaw scale discussed above Eq.(122) with the PQ breaking scale of Eq.(45).

Notably, in our model, the flavored-QCD axion mass (and its associated PQ breaking or seesaw scale) is closely linked to the soft SUSY-breaking mass. Our analysis covers the PQ scale $\langle\chi\rangle$ from roughly $10^{13}$ GeV to $10^{15}$ GeV due to Table-4, corresponding to $m_{3/2}$ values of 1 TeV to $10^{6}$ TeV, by considering Eq.(45) and Table-4. The given $U(1)_{X}$ quantum numbers can then be used to predict the branching ratio of $K^{+}\rightarrow\pi^{+}+a$ (Eq.(83)), as well as the axion coupling to photon (Eq.(86)) and electron (Eq.(95)). See Table-5, -6, and -7 for more details. Our proposed model’s predictions can be tested by current axion search experiments. KLASH Alesini:2019nzq is sensitive to the mass range of $0.27-0.93\,\mu{\rm eV}$ , while FLASH Gatti:2021cel covers the mass range of $0.5-1.5\,\mu{\rm eV}$ .

Refer to caption — Figure 1: Plots for axion-photon coupling $|g_{a\gamma\gamma}|$ as a function of the flavored-QCD axion mass $m_{a}$ for NO and IO. Orange shaded region and vertical red lines indicate the conventional QCD axion predictions and the exclusion region of various axion search experiments, respectively, see Ref.PDG .

Fig.1 illustrates plots of the axion-photon coupling $|g_{a\gamma\gamma}|$ as a function of the flavored-QCD axion mass $m_{a}$ for NO (left) and IO (right), respectively. Each plotted point corresponds to values listed in Tables-5, -6, and -7, which are consistent with the experimental constraints described in Eqs.(84), (87), and (96). And Fig.1 illustrates that certain data points in Table-5 (I-c, I-d, I-e) have been fully excluded by the ADMX experiment ADMX , while another data point (II) in Table-5 has been marginally excluded by the same experiment. Fig.2 shows plots for axion-electron coupling $|g_{aee}|$ as a function of the flavored-QCD axion mass $m_{a}$ for NO (left) and IO (right).

Table 5: The quark and lepton

U(1)_{X}

quantum numbers that satisfy experimental results, including Eqs.(79, 80, 84) and Table-3, as well as the

U(1)

-mixed gravitational anomaly-free condition for normal neutrino mass ordering. And QCD anomaly coefficient

N_{C}

, QCD axion decay constant

F_{a}=2|\langle\chi\rangle/N_{C}|

(with PQ breakdown or seesaw scale

\langle\chi\rangle

), the ratio of QCD and QED anomaly coefficients

E/N_{C}

, axion-electron coupling

g_{ae}

, and for

z=0.47

axion-photon coupling

g_{a\gamma\gamma}

, axion mass

m_{a}

, and branching ratio

{\rm Br}(K^{+}\rightarrow\pi^{+}+a)\equiv{\rm Br}(K\pi a)

U(1)_{X}

f_{u}

f_{c}

f_{d}

f_{s}

f_{b}

f_{e}

f_{\mu}

f_{\tau}

g_{e}

g_{\mu}

g_{\tau}

N_{C}

\frac{F_{a}}{\rm GeV}

\frac{E}{N_{C}}

\frac{g_{ae}}{10^{-17}}

\frac{|g_{a\gamma\gamma}|}{10^{-17}\,{\rm GeV}^{-1}}

\frac{m_{a}}{10^{-7}\,{\rm eV}}

{\rm Br}(K\pi a)

20

8

14

11

5

20

9

4

6

4

5

4.2^{+1.9}_{-1.7}\times 10^{-15}

I-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm 4

5\times 10^{12}

-\frac{5}{6}

36.13

66.00

10.89

I-b

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm

\pm

\mp

\pm 4

5\times 10^{12}

\frac{19}{6}

36.13

26.94

10.89

I-c

\mp

\pm

\mp

\pm

\pm

\pm

\mp

\pm

\pm

\pm

\pm

\pm 10

2\times 10^{12}

\frac{1}{15}

36.13

112.71

27.21

I-d

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\mp

\pm

\pm

\mp

\pm 10

2\times 10^{12}

-\frac{29}{15}

36.13

228.88

27.21

I-e

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\pm

\mp

\mp

\pm 10

2\times 10^{12}

-\frac{53}{15}

36.13

321.82

27.21

21

8

14

11

5

20

9

4

5

3

4

4.2^{+1.9}_{-1.7}\times 10^{-15}

\mp

\pm

\pm

\mp

\mp

\pm

\pm

\pm

\pm

\pm

\pm

\pm 15

1.3\times 10^{12}

-2

37.06

358.09

41.87

20

8

14

11

5

19

9

4

5

3

4

1.7^{+0.8}_{-0.7}\times 10^{-16}

III-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm

\mp

\pm

\pm 4

2.5\times 10^{13}

\frac{11}{3}

6.87

7.71

2.18

III-b

\mp

\pm

\pm

\mp

\pm

\pm

\pm

\mp

\mp

\mp

\mp

\pm 4

2.5\times 10^{13}

-\frac{16}{3}

6.87

34.11

2.18

III-c

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm 4

2.5\times 10^{13}

-\frac{1}{3}

6.87

10.88

2.18

III-d

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\pm

\mp

\mp

\pm 10

10^{13}

-\frac{10}{3}

6.87

62.04

5.44

III-e

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\mp

\pm

\mp

\pm

\pm 10

10^{13}

-\frac{26}{15}

6.87

43.45

5.44

20

8

14

11

5

20

9

4

4

2

3

4.2^{+1.9}_{-1.7}\times 10^{-17}

IV-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\mp

\mp

\pm

\pm 4

5\times 10^{13}

-\frac{5}{6}

3.61

6.60

1.09

IV-b

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm

\mp

\pm

\pm 4

5\times 10^{13}

\frac{19}{6}

3.61

2.69

1.09

IV-c

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\mp

\pm

\mp

\pm

\pm 10

2\times 10^{13}

-\frac{29}{15}

3.61

22.89

2.72

IV-d

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\mp

\mp

\pm

\pm 10

2\times 10^{13}

-\frac{53}{15}

3.61

32.18

2.72

IV-e

\mp

\pm

\pm

\mp

\mp

\pm

\pm

\pm

\pm

\pm

\pm

\pm 14

1.43\times 10^{13}

-\frac{7}{3}

3.61

35.30

3.81

Table 6: The same as in Table 5.

U(1)_{X}

f_{u}

f_{c}

f_{d}

f_{s}

f_{b}

f_{e}

f_{\mu}

f_{\tau}

g_{e}

g_{\mu}

g_{\tau}

N_{C}

\frac{F_{a}}{\rm GeV}

\frac{E}{N_{C}}

\frac{g_{ae}}{10^{-17}}

\frac{|g_{a\gamma\gamma}|}{10^{-17}\,{\rm GeV}^{-1}}

\frac{m_{a}}{10^{-7}\,{\rm eV}}

{\rm Br}(K\pi a)

21

8

14

11

5

19

9

4

4

2

3

4.2^{+1.9}_{-1.7}\times 10^{-17}

V-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\mp

\pm

\pm

\pm 5

4\times 10^{13}

\frac{4}{15}

3.43

5.05

1.36

V-b

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm

\pm

\pm

\pm 5

4\times 10^{13}

\frac{52}{15}

3.43

4.24

1.36

V-c

\mp

\pm

\pm

\mp

\pm

\pm

\pm

\mp

\mp

\mp

\mp

\pm 5

4\times 10^{13}

-\frac{56}{15}

3.43

16.67

1.36

V-d

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\mp

\pm

\pm

\pm 11

1.8\times 10^{13}

-\frac{92}{33}

3.43

30.64

2.99

V-e

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\mp

\pm

\pm

\pm

\pm 11

1.8\times 10^{13}

-\frac{4}{3}

3.43

21.34

2.99

20

8

14

11

5

19

9

4

3

1

2

1.7^{+0.8}_{-0.7}\times 10^{-18}

VI-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\mp

\mp

\pm

\pm 4

2.5\times 10^{14}

-\frac{1}{3}

0.69

1.09

0.22

VI-b

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm

\pm

\pm

\pm 4

2.5\times 10^{14}

\frac{11}{3}

0.69

0.77

0.22

VI-c

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\mp

\mp

\pm

\pm 10

10^{14}

-\frac{10}{3}

0.69

6.20

0.54

VI-d

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\mp

\pm

\pm

\pm

\pm 10

10^{14}

-\frac{26}{15}

0.69

4.35

0.54

21

8

14

11

5

20

9

4

3

1

2

1.7^{+0.8}_{-0.7}\times 10^{-18}

VII-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\mp

\pm

\pm

\pm 5

2\times 10^{14}

-\frac{2}{15}

0.72

1.24

0.27

VII-b

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\pm

\mp

\mp

\pm 11

9\times 10^{13}

-\frac{98}{33}

0.72

6.36

0.60

21

8

14

11

5

19

9

4

2

0

1

4.2^{+1.9}_{-1.7}\times 10^{-19}

VIII-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\pm

0

\mp

\pm 5

4\times 10^{14}

\frac{4}{15}

0.34

0.51

0.14

VIII-b

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\pm

0

\mp

\pm 11

1.8\times 10^{14}

-\frac{92}{33}

0.34

3.06

0.30

20

8

14

11

5

20

9

4

2

0

1

4.2^{+1.9}_{-1.7}\times 10^{-19}

IX-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\mp

0

\mp

\pm 4

5\times 10^{14}

-\frac{5}{6}

0.36

0.66

0.11

IX-b

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\mp

0

\mp

\pm 10

2\times 10^{14}

-\frac{53}{15}

0.36

3.22

0.27

V.1 Quark and charged-lepton

The Yukawa matrices for charged fermions in the SM, as given in Eqs.(68, 76, 91), are taken at the scale of $U(1)_{X}$ symmetry breakdown. Hence their masses are subject to quantum corrections. Subsequently, these matrices are run down to $m_{t}$ and diagonalized. We assume that the Yukawa matrices at the scale of $U(1)_{X}$ breakdown are the same as those at the scale $m_{t}$ , since the one-loop renormalization group running effect on observables for hierarchical mass spectra is expected to be negligible. The low-energy Yukawa couplings required for experimental values are obtained from the physical masses and mixing angles compiled by the PDG PDG and CKMfitter ckm .

We have thirteen physical observables in the quark and charged-lepton sector: $m_{d},m_{s},m_{b}$ , $m_{u},m_{c},m_{t},m_{e},m_{\mu},m_{\tau}$ , and $\theta^{q}_{12},\theta^{q}_{23},\theta^{q}_{13},\delta^{q}_{CP}$ . These observables are used to determine thirteen effective model parameters: 21 parameters ( $|\alpha_{d}|$ , $|\alpha_{s}|$ , $|\alpha_{b}|$ , $|\tilde{\alpha}_{d}|$ , $|\tilde{\alpha}_{s}|$ , $|\tilde{\alpha}_{b}|$ , $\alpha_{e}$ , $\alpha_{\mu}$ , $\alpha_{\tau}$ , $\alpha_{t}$ , $\alpha_{c}$ , $\alpha_{u}$ ; $\arg(\alpha_{d}),\arg(\alpha_{s}),\arg(\tilde{\alpha}_{d}),\arg(\tilde{\alpha}_{s}),\arg(\tilde{\alpha}_{b})$ ; $\Delta_{\chi}$ , $\tan\beta$ ; ${\rm Re[\tau]},{\rm Im}[\tau]$ ) among whihc 8 parameters are fixed by quantum numbers ( $f_{e,\mu,\tau},f_{d,s,b},f_{u,c}$ ). Using highly precise data as constraints for both quarks and charged leptons, with the exception of the quark Dirac CP phase, as described in Eqs.(79), (80), and (93), we scanned all parameter ranges and determined that

	$\displaystyle\Delta_{\chi}=[0.596,0.602]\,,\qquad\tan\beta=[6.8,7.3]\,,$
	$\displaystyle\tau=(0.0001\sim 0.046)+(1.0906\sim 1.1086)i\,.$		(123)

${\rm Re}(\tau)$ contributes to the phase of the Yukawa coupling, while ${\rm Im}(\tau)$ influences the magnitude of Yukawa coupling, as demonstrated in Eqs.(68, 76, 91, 105, 112). When $\langle\tau\rangle=i$ , resulting in real values of $x,y$ with $Y^{(2)}_{\bf 3}=Y_{1}(i)(1,1-\sqrt{3},-2+\sqrt{3})$ , it becomes apparent that the empirical results of quark masses and CKM mixing angles cannot be satisfied due to the overall factors in Eq.(76) being real. Therefore, a departure of $\tau$ from $i$ is necessary. Fig.3 shows how the quark Dirac CP phase $\delta^{q}_{CP}$ behaves based on certain constrained parameters. Our model predicts that $\delta^{q}_{CP}$ falls between $38^{\circ}$ and $87^{\circ}$ , which aligns well with experimental data. The horizontal black-dotted lines in Fig.3 represent the $3\sigma$ experimental bound for $\delta^{q}_{CP}$ . Notably, the effective Yukawa coefficients satisfying the experimental data fall well within the bound specified in Eq.(60), as shown in the top left panel of Fig.3. This reflects that these coefficients have a natural size of unity, as stated in Eq.(2).

We choose reference values, for example, that satisfy the experimental data :

\displaystyle\Delta_{\chi}=0.597\,,\qquad\tau=0.0074+1.0997i\,,\qquad\tan\beta=6.8

(124)

which result in effective Yukawa coefficients from Eq.(60) satisfying $0.45\lesssim|\alpha_{i}|\lesssim 1.55$ . With the inputs

	$\displaystyle\arg(\alpha_{d})=1.007,~{}\arg(\alpha_{s})=2.232,~{}\arg(\tilde{\alpha}_{d})=4.723,~{}\arg(\tilde{\alpha}_{s})=3.388,~{}\arg(\tilde{\alpha}_{b})=1.164,$
	$\displaystyle\alpha_{u}=1.320~{}\text{for $\|f_{u}\|=21$}(\alpha_{u}=0.788~{}\text{for $\|f_{u}\|=20$}),~{}\alpha_{c}=0.950,~{}\alpha_{t}=1.006,$
	$\displaystyle\|\alpha_{d}\|=1.039,~{}\|\alpha_{s}\|=1.218,~{}\|\alpha_{b}\|=0.790,~{}\|\tilde{\alpha}_{d}\|=0.896,~{}\|\tilde{\alpha}_{s}\|=0.822,~{}\|\tilde{\alpha}_{b}\|=1.158,$		(125)

we obtain the mixing angles and Dirac CP phase $\theta^{q}_{12}=12.980^{\circ}$ , $\theta^{q}_{23}=2.320^{\circ}$ , $\theta^{q}_{13}=0.218^{\circ}$ , $\delta^{q}_{CP}=64.216^{\circ}$ compatible with the $3\sigma$ Global fit of CKMfitter ckm , see Eq.(79); the quark masses $m_{d}=4.593$ MeV, $m_{s}=103.819$ MeV, $m_{b}=4.206$ GeV, $m_{u}=2.164$ MeV, $m_{c}=1.271$ GeV, and $m_{t}=173.1$ GeV compatible with the values in PDG PDG , see Eq.(80). Here, without loss of generality, the up-type quark masses $m_{u}$ , $m_{c}$ , and $m_{t}$ are a one-to-one correspondence with $\alpha_{u}$ , $\alpha_{c}$ , and $\alpha_{t}$ , which have been taken real, and we have set $\arg(\alpha_{b})=0$ .

The masses of the charged leptons $m_{e}$ , $m_{\mu}$ , and $m_{\tau}$ are in a one-to-one correspondence with the real parameters $\alpha_{e}$ , $\alpha_{\mu}$ , and $\alpha_{\tau}$ from Eq.(92). Using the numerical results of Eq.(124) from the quark sector, with the inputs

\displaystyle\alpha_{e}=1.268~{}\text{for}\,|f_{e}|=20\,(\alpha_{e}=0.757~{}\text{for}\,|f_{e}|=19),\alpha_{\mu}=0.900,\alpha_{\tau}=1.148\,,

(126)

we obtain the charged lepton masses, which well agree with the empirical values of Eq.(93).

Table 7: The same as in Table 5, except for the inverted neutrino mass ordering.

U(1)_{X}

f_{u}

f_{c}

f_{d}

f_{s}

f_{b}

f_{e}

f_{\mu}

f_{\tau}

g_{e}

g_{\mu}

g_{\tau}

N_{C}

\frac{F_{a}}{\rm GeV}

\frac{E}{N_{C}}

\frac{g_{ae}}{10^{-17}}

\frac{|g_{a\gamma\gamma}|}{10^{-17}\,{\rm GeV}^{-1}}

\frac{m_{a}}{10^{-7}\,{\rm eV}}

{\rm Br}(K\pi a)

19

8

14

11

5

20

9

4

2

3

3

1.7^{+0.8}_{-0.7}\times 10^{-16}

I-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm

\begin{array}[]{ll}\mp\\ \pm\end{array}

\begin{array}[]{ll}\pm\\ \mp\end{array}

\pm 3

3.33\times 10^{13}

\frac{10}{3}

7.23

4.62

1.63

I-b

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\mp

\pm

\begin{array}[]{ll}\mp\\ \pm\end{array}

\begin{array}[]{ll}\pm\\ \mp\end{array}

\pm 9

1.11\times 10^{13}

-\frac{22}{9}

7.23

46.54

4.90

20

8

14

11

5

19

9

4

2

3

3

1.7^{+0.8}_{-0.7}\times 10^{-16}

II-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\mp

\begin{array}[]{ll}\mp\\ \pm\end{array}

\begin{array}[]{ll}\pm\\ \mp\end{array}

\pm 4

2.5\times 10^{13}

-\frac{1}{3}

6.87

10.88

2.18

II-b

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\mp

\begin{array}[]{ll}\mp\\ \pm\end{array}

\begin{array}[]{ll}\pm\\ \mp\end{array}

\pm 10

10^{13}

-\frac{10}{3}

6.87

62.04

5.44

20

8

14

11

5

20

9

4

1

2

2

4.2^{+1.9}_{-1.7}\times 10^{-17}

III-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm

\pm

\pm

\pm 4

5\times 10^{13}

\frac{19}{6}

3.61

6.60

1.09

III-b

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\pm 4

5\times 10^{13}

-\frac{5}{6}

3.61

2.69

1.09

III-c

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\mp

\pm

\pm

\pm

\pm 10

2\times 10^{13}

-\frac{29}{15}

3.61

22.89

2.72

III-d

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\pm

\mp

\mp

\pm 10

2\times 10^{13}

-\frac{53}{15}

3.61

32.18

2.72

19

8

14

11

5

19

9

4

1

2

2

4.2^{+1.9}_{-1.7}\times 10^{-17}

IV-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

\mp

\pm

\pm

\pm 3

6.67\times 10^{13}

4

3.43

3.47

8.16

IV-b

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\pm

\mp

\mp

\mp

\pm 3

6.67\times 10^{13}

-\frac{4}{3}

3.43

5.82

8.16

IV-c

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\mp

\mp

\pm

\pm

\pm 9

2.22\times 10^{13}

-\frac{20}{9}

3.43

22.11

2.45

IV-d

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\pm

\mp

\mp

\mp

\pm 9

2.22\times 10^{13}

-4

3.43

31.40

2.45

19

8

14

11

5

20

9

4

0

1

1

1.7^{+0.8}_{-0.7}\times 10^{-18}

V-a

\mp

\pm

\pm

\mp

\pm

\pm

\mp

\mp

0

\pm

\pm

\pm 3

3.33\times 10^{14}

\frac{10}{3}

0.72

0.46

0.16

V-b

\mp

\pm

\mp

\pm

\pm

\pm

\pm

\mp

0

\pm

\pm

\pm 9

1.11\times 10^{14}

-\frac{22}{9}

0.72

4.65

0.49

V.2 Neutrino

The seesaw mechanism in Eq.(113) operates at the $U(1)_{X}$ symmetry breakdown scale, while its implications are measured by experiments below EW scale. Therefore, quantum corrections to neutrino masses and mixing angles can be crucial, especially for degenerate neutrino masses Antusch:2005gp . However, based on our observation that the neutrino mass spectra exhibit hierarchy at the scale of $U(1)_{X}$ breakdown (as depicted in Fig.4), we can safely assume that the renormalization group running effect on observables can be ignored.

Neutrino oscillation experiments currently aim to make precise measurements of the Dirac CP-violating phase $\delta_{CP}$ and atmospheric mixing angle $\theta_{23}$ . Using our model, we investigate which values of $\delta_{CP}$ and $\theta_{23}$ can predict the mass hierarchy of neutrinos (NO or IO) and identify observables that can be tested in current and next-generation experiments. To explore the parameter spaces, we scan the precision constraints $\{\theta_{13}$ , $\theta_{23}$ , $\theta_{12}$ , $\Delta m^{2}_{\rm Sol}$ , $\Delta m^{2}_{\rm Atm}\}$ at $3\sigma$ from Table-3. Using the reference values from Eq.(124) in the quark and charged-lepton sectors, we determine the input parameter spaces of Eq.(122) for both NO and IO at the $U(1)_{X}$ breaking scale : for example, taking $\langle\chi\rangle=5\times 10^{13}$ GeV (see, Table-6, 6, and -7) and for NO

	$\displaystyle\tilde{\beta}_{2}=[0.48,1.10],\quad\tilde{\beta}_{3}=[0.61,1.18]\,,\quad\arg(\tilde{\beta}_{2(3)})=[0,2\pi]$
	$\displaystyle m_{0}\langle\chi\rangle/v^{2}_{u}=[0.59,1.84],\quad\gamma=[0.37,0.84],\quad\gamma^{\prime}=[0.37,0.65]$
	$\displaystyle\arg(\gamma)=[4.42,5.55],\quad\arg(\gamma^{\prime})=[1.30,2.19];$		(127)

where $\beta_{1}=[0.97,1.47]$ , $\beta_{2}=[0.61,1.11]$ , and $\beta_{3}=[0.69,1.15]$ ; for IO

	$\displaystyle\tilde{\beta}_{2}=[0.55,0.87],\quad\tilde{\beta}_{3}=[0.56,0.86]\,,\quad\arg(\tilde{\beta}_{2(3)})=[0,2\pi]$
	$\displaystyle m_{0}\langle\chi\rangle/v^{2}_{u}=[0.81,1.39],\quad\gamma=[0.758,1.154],\quad\gamma^{\prime}=[0.350,0.618]$
	$\displaystyle\arg(\gamma)=[2.92,4.19],\quad\arg(\gamma^{\prime})=[0.1,1.29]\,\&\,[5.66,6.26]\,.$		(128)

where $\beta_{1}=[0.98,1.29]$ , $\beta_{2}=[0.69,0.95]$ , and $\beta_{3}=[0.68,0.93]$ . For these parameter regions, we investigate how $0\nu\beta\beta$ -decay rate and Dirac CP phase can be determined for the normal and inverted mass ordering. These predictions are represented by crosses and X-marks for NO and IO, respectively, in Fig.5. Referring to the two-dimensional allowed regions at $3\sigma$ presented in Ref.Esteban:2020cvm , we note that the most favored regions correspond to $\delta_{CP}\sim 250^{\circ}$ , whereas there are no favored regions with respect to $\theta_{23}$ . Ongoing experiments like DUNE DUNE:2018tke , as well as proposed next-generation experiments such as Hyper-K Hyper-Kamiokande:2018ofw , are poised to greatly reduce uncertainties in the values of $\theta_{23}$ and $\delta_{CP}$ , providing a rigorous test for our proposed model. Furthermore, ongoing and future experiments on $0\nu\beta\beta$ -decay like NEXT NEXT:2020amj , SNO $+$ SNO:2022trz , KamLAND-Zen KamLAND-Zen:2022tow , Theia Theia:2019non , SuperNEMO Arnold:2004xq may soon reach a sensitivity to exclude the inverted mass ordering of our model. Cosmological and astrophysical measurements provide powerful constraints on the sum of neutrino masses. The upper bound on the sum of the three active neutrino masses can be summarized as $\sum m_{\nu}=m_{\nu_{1}}+m_{\nu_{2}}+m_{\nu_{3}}<0.120$ eV at $95\%$ CL for TT, TE, EE+lowE+lensing+BAO Planck:2018vyg . Fig.6 illustrates that the sum of neutrino masses lies in the range of $0.0582$ to $0.0605$ eV for the NO when $0.0003397\lesssim m_{\rm lightest}=m_{\nu_{1}}\lesssim 0.0007945$ eV. Conversely, for the IO, the sum of neutrino masses falls within the range of $0.1003$ to $0.1038$ eV when $0.001125\lesssim m_{\rm lightest}=m_{\nu_{3}}\lesssim 0.001647$ eV.

VI conclusion

We proposed a minimal extension of a modular invariant model that incorporates sterile neutrinos and a QCD axion (as a strong candidate of dark matter) into the SM to account for the mass and mixing hierarchies of quarks and leptons, as well as the strong CP problem. Our model, based on the 4D effective action, features the $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ symmetry. To ensure the reliability of our model, we have examined the modular forms of the superpotential, corrected by Kähler transformation, under the $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ symmetry, while also considering the modular and $U(1)_{X}$ anomaly-free conditions. The model features a minimal set of fields that transform based on representations of $G_{\rm SM}\times\Gamma_{N}\times U(1)_{X}$ , and includes modular forms of level $N$ . These modular forms act as Yukawa couplings and transform under the modular group $\Gamma_{N}$ . Our numerical analysis guarantees that, in the supersymmetric limit, all Yukawa coefficients in the superpotential are complex numbers with a unit absolute value, implying a democratic distribution.

We demonstrated, as an explicit example, a level 3 modular form-induced superpotential by introducing minimal supermultiplets. The extension includes right-handed neutrinos ( $N^{c}$ ) and SM gauge singlet scalar fields ( $\chi$ and $\tilde{\chi}$ ) with zero modular weight and ( $+$ and $-$ ) charge under $U(1)_{X}$ . These scalar fields are crucial in generating the QCD axion, heavy neutrino mass, and fermion mass hierarchy. Modular invariance of both the superpotential and Kähler potential allows for Kähler transformation to correct modular form weight in the superpotential, enabling a $\tau$ -independent superpotential for the scalar potential. The sterile neutrinos are introduced to satisfy the $U(1)_{X}$ -mixed gravitational anomaly-free condition, explain small active neutrino masses via the seesaw mechanism, and provide a well-motivated PQ symmetry breaking scale. As the fields $\chi(\tilde{\chi})$ have modular weights of zero, any additional correction terms arising from higher weight modular forms are not permitted in the superpotentials. However, the combination $\chi\tilde{\chi}$ can trigger higher-order corrections that are permissible and do not modify the leading-order flavor structures. Taking into account both SUSY-breaking effects and supersymmetric next-leading-order Planck-suppressed terms, we have determined the low axion decay constant (or seesaw scale). This leads to an approximate range for the PQ scale $\langle\chi\rangle$ (equivalently, the seesaw scale) of $10^{13}$ GeV to $10^{15}$ GeV for $m_{3/2}$ values between 1 TeV and $10^{6}$ TeV, see Table-5, -6 and -7. Interestingly enough, in our model, the PQ breaking scale (or axion mass) is closely linked to the seesaw scale and the soft SUSY breaking mass. Our model with $E/N_{C}$ could be tested by ongoing experiments such as KLASH Alesini:2019nzq and FLASH Gatti:2021cel , see Fig.1 and 2, by considering the scale of $U(1)_{X}$ breakdown.

We explored numerical values of physical parameters that satisfy the highly precise data on the mass of quarks and charged-leptons, as well as the quark mixing angles, except for the quark Dirac CP phase. Our model predicts that the value of $\delta^{q}_{CP}$ falls within the range of $38^{\circ}$ to $87^{\circ}$ , which is consistent with experimental data. Notably, the effective Yukawa coefficients satisfying the experimental data fall well within the bound specified in Eq.(60), as shown in the top left panel of Fig.3. This suggests that our assumption, as stated in Eq.(2), that the Yukawa coefficients have a natural size of unity is plausible. Using precise neutrino oscillation data as constraints, we investigated how the $0\nu\beta\beta$ -decay rate and Dirac CP phase could be determined for the normal and inverted mass ordering in the neutrino sector. Referring to the $3\sigma$ allowed regions in Ref.Esteban:2020cvm , we note that the most favored regions for our proposed model are $\delta_{CP}\sim 250^{\circ}$ , with no favored regions with respect to $\theta_{23}$ . Ongoing experiments, such as DUNE DUNE:2018tke and proposed next-generation experiments, such as Hyper-K Hyper-Kamiokande:2018ofw , are expected to greatly reduce uncertainties in the values of $\theta_{23}$ and $\delta_{CP}$ , providing a rigorous test for our model. Additionally, ongoing and future experiments on $0\nu\beta\beta$ -decay, such as NEXT NEXT:2020amj , SNO $+$ SNO:2022trz , KamLAND-Zen KamLAND-Zen:2022tow , Theia Theia:2019non , and SuperNEMO Arnold:2004xq , may soon have the sensitivity to exclude the inverted mass ordering in our model.

Acknowledgements.

YHA was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIT) (No.2020R1A2C1010617). SKK was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIT) (No. 2019R1A2C1088953, No.2020K1A3A7A09080135, No.2023R1A2C1006091).

Appendix A The group $A_{4}$

The group $A_{4}$ is the symmetry group of the tetrahedron, isomorphic to the finite group of the even permutations of four objects. The group $A_{4}$ has two generators, denoted $S$ and $T$ , satisfying the relations $S^{2}=T^{3}=(ST)^{3}={\bf 1}$ . In the three-dimensional complex representation, $S$ and $T$ are given by

\displaystyle S=\frac{1}{3}{\left(\begin{array}[]{ccc}-1&2&2\\ 2&-1&2\\ 2&2&-1\end{array}\right)}\,,\qquad T={\left(\begin{array}[]{ccc}1&0&0\\ 0&\omega&0\\ 0&0&\omega^{2}\end{array}\right)}\,,

(135)

where $\omega=e^{i2\pi/3}=-1/2+i\sqrt{3}/2$ is a complex cubic-root of unity. $A_{4}$ has four irreducible representations: three singlets ${\bf 1},{\bf 1}^{\prime}$ , and ${\bf 1}^{\prime\prime}$ and one triplet ${\bf 3}$ . An $A_{4}$ singlet ${\bf a}$ is invariant under the action of $S$ ( $S{\bf a}={\bf a}$ ), while the action of $T$ produces $T{\bf a}={\bf a}$ for ${\bf 1}$ , $T{\bf a}=\omega{\bf a}$ for ${\bf 1}^{\prime}$ , and $T{\bf a}=\omega^{2}{\bf a}$ for ${\bf 1}^{\prime\prime}$ . Products of two $A_{4}$ representations decompose into irreducible representations according to the following multiplication rules: ${\mathbf{3}}\otimes{\mathbf{3}}={\mathbf{3}}_{s}\oplus{\mathbf{3}}_{a}\oplus{\mathbf{1}}\oplus{\mathbf{1}}^{\prime}\oplus{\mathbf{1}}^{\prime\prime}$ , ${\mathbf{1}}^{\prime}\oplus{\mathbf{1}}^{\prime}={\mathbf{1}}^{\prime\prime}$ and ${\mathbf{1}}^{\prime\prime}\oplus{\mathbf{1}}^{\prime\prime}={\mathbf{1}}^{\prime}$ . Explicitly, if $(a_{1},a_{2},a_{3})$ and $(b_{1},b_{2},b_{3})$ denote two $A_{4}$ triplets, then we have Eq.(30).

Appendix B The CKM mixing matrix

The CKM mixing matrix is given in the Wolfenstein parametrization Wolfenstein:1983yz by

\displaystyle V_{\rm CKM}={\left(\begin{array}[]{ccc}1-\frac{1}{2}\lambda^{2}&\lambda&A_{d}\lambda^{3}(\rho-i\eta)\\ -\lambda&1-\frac{1}{2}\lambda^{2}&A_{d}\lambda^{2}\\ A_{d}\lambda^{3}(1-\rho-i\eta)&-A_{d}\lambda^{2}&1\end{array}\right)}+{\cal O}(\lambda^{4})\,.

(139)

where $\lambda=0.22500^{+0.00082}_{-0.00063}$ , $A_{d}=0.813^{+0.026}_{-0.018}$ , $\bar{\rho}=\rho/(1-\lambda^{2}/2)=0.157^{+0.036}_{-0.018}$ , and $\bar{\eta}=\eta/(1-\lambda^{2}/2)=0.347^{+0.030}_{-0.020}$ with $3\sigma$ errors ckm .

References

(1) P. Minkowski, Phys. Lett. B 67, 421 (1977); T. Yanagida, in Proc. of the Workshop on Unified Theories and Baryon Number in the Universe, ed.O. Sawada and A. Sugamoto, 95 (KEK, Japan, 1979); M. Gell-Mann, P. Ramond and R. Slansky, in Supergravity, ed. P. Nieuwenhuizen and D. Freeman (North Holland, Amsterdam, 1979); R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, (1980) 912.
(2) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977); R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977).
(3) F. Feruglio, doi:10.1142/9789813238053_0012 [arXiv:1706.08749 [hep-ph]].
(4) S. Ferrara, D. Lust, A. D. Shapere and S. Theisen, Phys. Lett. B 225 (1989), 363; S. Ferrara, D. Lust and S. Theisen, Phys. Lett. B 233 (1989), 147-152.
(5) T. Kobayashi, N. Omoto, Y. Shimizu, K. Takagi, M. Tanimoto and T. H. Tatsuishi, JHEP 11 (2018), 196 [arXiv:1808.03012 [hep-ph]]; S. J. D. King and S. F. King, JHEP 09 (2020), 043 [arXiv:2002.00969 [hep-ph]]; F. J. de Anda, S. F. King and E. Perdomo, Phys. Rev. D 101 (2020) no.1, 015028 [arXiv:1812.05620 [hep-ph]]; C. Y. Yao, J. N. Lu and G. J. Ding, JHEP 05 (2021), 102 [arXiv:2012.13390 [hep-ph]]; H. Okada and M. Tanimoto, Eur. Phys. J. C 81 (2021) no.1, 52 [arXiv:1905.13421 [hep-ph]]; C. Y. Yao, J. N. Lu and G. J. Ding, JHEP 05 (2021), 102 [arXiv:2012.13390 [hep-ph]]; S. Kikuchi, T. Kobayashi, H. Otsuka, M. Tanimoto, H. Uchida and K. Yamamoto, Phys. Rev. D 106, no.3, 035001 (2022) [arXiv:2201.04505 [hep-ph]]; K. Hoshiya, S. Kikuchi, T. Kobayashi and H. Uchida, Phys. Rev. D 106 (2022) no.11, 115003 [arXiv:2209.07249 [hep-ph]]; A. Baur, H. P. Nilles, S. Ramos-Sanchez, A. Trautner and P. K. S. Vaudrevange, JHEP 09 (2022), 224 [arXiv:2207.10677 [hep-ph]]; X. K. Du and F. Wang, JHEP 01 (2023), 036 [arXiv:2209.08796 [hep-ph]]; Y. Abe, T. Higaki, J. Kawamura and T. Kobayashi, [arXiv:2302.11183 [hep-ph]]; H. Okada and M. Tanimoto, Phys. Dark Univ. 40 (2023), 101204 [arXiv:2005.00775 [hep-ph]]; S. T. Petcov and M. Tanimoto, [arXiv:2306.05730 [hep-ph]]; S. T. Petcov and M. Tanimoto, [arXiv:2212.13336 [hep-ph]].
(6) G. J. Ding, S. F. King and X. G. Liu, JHEP 09 (2019), 074 [arXiv:1907.11714 [hep-ph]]; Y. H. Ahn, S. K. Kang, R. Ramos and M. Tanimoto, Phys. Rev. D 106 (2022) no.9, 095002 [arXiv:2205.02796 [hep-ph]].
(7) F. Feruglio, A. Strumia and A. Titov, [arXiv:2305.08908 [hep-ph]].
(8) T. Kobayashi and H. Otsuka, Phys. Lett. B 807 (2020), 135554 [arXiv:2002.06931 [hep-ph]].
(9) R. de Adelhart Toorop, F. Feruglio and C. Hagedorn, Nucl. Phys. B 858, 437-467 (2012) [arXiv:1112.1340 [hep-ph]].
(10) L.E. Ibanez and A.M. Uranga, String theory and particle physics: an introduction to string phenomenology, Cambridge University Press, Cambridge U.K. (2012).
(11) J. P. Derendinger, S. Ferrara, C. Kounnas and F. Zwirner, Nucl. Phys. B 372, 145-188 (1992).
(12) M. Dine, N. Seiberg and E. Witten, Nucl. Phys. B 289, 589-598 (1987); P. Svrcek and E. Witten, JHEP 06, 051 (2006) [arXiv:hep-th/0605206 [hep-th]].
(13) M. B. Green and J. H. Schwarz, Phys. Lett. 149B, 117 (1984); M. B. Green and J. H. Schwarz, Nucl. Phys. B 255, 93-114 (1985); M. B. Green, J. H. Schwarz and P. C. West, Nucl. Phys. B 254, 327-348 (1985).
(14) C. P. Burgess, R. Kallosh and F. Quevedo, JHEP 10, 056 (2003) [arXiv:hep-th/0309187 [hep-th]].
(15) D. V. Nanopoulos, K. A. Olive, M. Srednicki and K. Tamvakis, Phys. Lett. B 124, 171 (1983); H. Murayama, H. Suzuki and T. Yanagida, Phys. Lett. B 291, 418-425 (1992); K. Choi, E. J. Chun and J. E. Kim, Phys. Lett. B 403, 209-217 (1997) [arXiv:hep-ph/9608222 [hep-ph]].
(16) Y. H. Ahn, Phys. Rev. D 100, no. 1, 015002 (2019) [arXiv:1706.09707 [hep-ph]].
(17) Y. H. Ahn, Phys. Rev. D 96, no. 1, 015022 (2017) [arXiv:1611.08359 [hep-ph]].
(18) Y. H. Ahn, Phys. Rev. D 91, 056005 (2015) [arXiv:1410.1634 [hep-ph]].
(19) M. Kamionkowski and J. March-Russell, Phys. Lett. B 282, 137 (1992) [hep-th/9202003]; R. Kallosh, A. D. Linde, D. A. Linde and L. Susskind, Phys. Rev. D 52, 912 (1995) [hep-th/9502069]; G. Dvali, hep-th/0507215.
(20) Y. H. Ahn, Phys. Rev. D 98, no. 3, 035047 (2018) [arXiv:1804.06988 [hep-ph]]; Y. H. Ahn, Nucl. Phys. B 939, 534 (2019) [arXiv:1802.05044 [hep-ph]]. Y. H. Ahn and X. Bi, Nucl. Phys. B 960, 115210 (2020) [arXiv:1912.09038 [hep-ph]].
(21) Y. H. Ahn, S. K. Kang and H. M. Lee, Phys. Rev. D 106, no.7, 075029 (2022) [arXiv:2112.13392 [hep-ph]].
(22) S. Gukov, C. Vafa and E. Witten, Nucl. Phys. B 584, 69-108 (2000) [erratum: Nucl. Phys. B 608, 477-478 (2001)] [arXiv:hep-th/9906070 [hep-th]].
(23) S. B. Giddings, S. Kachru and J. Polchinski, Phys. Rev. D 66, 106006 (2002) [arXiv:hep-th/0105097 [hep-th]].
(24) K. Dasgupta, G. Rajesh and S. Sethi, JHEP 08, 023 (1999) [arXiv:hep-th/9908088 [hep-th]].
(25) D. Cremades, L. E. Ibanez and F. Marchesano, JHEP 05, 079 (2004) [arXiv:hep-th/0404229 [hep-th]].
(26) D. Cremades, L. E. Ibanez and F. Marchesano, JHEP 07 (2003), 038 [arXiv:hep-th/0302105 [hep-th]]; R. Blumenhagen, M. Cvetic, P. Langacker and G. Shiu, Ann. Rev. Nucl. Part. Sci. 55 (2005), 71-139 [arXiv:hep-th/0502005 [hep-th]]; S. A. Abel and M. D. Goodsell, JHEP 10 (2007), 034 [arXiv:hep-th/0612110 [hep-th]]; R. Blumenhagen, B. Kors, D. Lust and S. Stieberger, Phys. Rept. 445 (2007), 1-193 [arXiv:hep-th/0610327 [hep-th]]; F. Marchesano, Fortsch. Phys. 55 (2007), 491-518 [arXiv:hep-th/0702094 [hep-th]]; I. Antoniadis, A. Kumar and B. Panda, Nucl. Phys. B 823 (2009), 116-173 [arXiv:0904.0910 [hep-th]]; T. Kobayashi, S. Nagamoto and S. Uemura, PTEP 2017 (2017) no.2, 023B02 [arXiv:1608.06129 [hep-th]].
(27) A. Sagnotti, Phys. Lett. B 294, 196-203 (1992) [arXiv:hep-th/9210127 [hep-th]].
(28) L. E. Ibanez, R. Rabadan and A. M. Uranga, Nucl. Phys. B 542, 112-138 (1999) [arXiv:hep-th/9808139 [hep-th]].
(29) E. Poppitz, Nucl. Phys. B 542, 31-44 (1999) [arXiv:hep-th/9810010 [hep-th]].
(30) Z. Lalak, S. Lavignac and H. P. Nilles, Nucl. Phys. B 559, 48-70 (1999) [arXiv:hep-th/9903160 [hep-th]].
(31) Y. H. Ahn, Phys. Rev. D 93, no. 8, 085026 (2016) [arXiv:1604.01255 [hep-ph]].
(32) A. Achucarro, B. de Carlos, J. A. Casas and L. Doplicher, JHEP 06, 014 (2006) [arXiv:hep-th/0601190 [hep-th]]; E. Dudas and Y. Mambrini, JHEP 10, 044 (2006) [arXiv:hep-th/0607077 [hep-th]].
(33) D. Alesini, D. Babusci, P. Beltrame, S.J., F. Björkeroth, F. Bossi, P. Ciambrone, G. Delle Monache, D. Di Gioacchino, P. Falferi and A. Gallo, et al. [arXiv:1911.02427 [physics.ins-det]].
(34) C. Gatti, P. Gianotti, C. Ligi, M. Raggi and P. Valente, Universe 7, no.7, 236 (2021).
(35) E. Ma and G. Rajasekaran, Phys. Rev. D 64, 113012 (2001) [arXiv:hep-ph/0106291 [hep-ph]]; G. Altarelli and F. Feruglio, Nucl. Phys. B 720, 64-88 (2005) [arXiv:hep-ph/0504165 [hep-ph]]; P. F. Harrison and W. G. Scott, Phys. Lett. B 557, 76 (2003) doi:10.1016/S0370-2693(03)00183-7 [arXiv:hep-ph/0302025 [hep-ph]]; X. G. He, Y. Y. Keum and R. R. Volkas, JHEP 04, 039 (2006) [arXiv:hep-ph/0601001 [hep-ph]].
(36) Y. H. Ahn, S. Baek and P. Gondolo, Phys. Rev. D 86, 053004 (2012) [arXiv:1207.1229 [hep-ph]]; Y. H. Ahn, S. K. Kang and C. S. Kim, Phys. Rev. D 87, no.11, 113012 (2013) [arXiv:1304.0921 [hep-ph]]; Z. z. Xing, Phys. Lett. B 533, 85-93 (2002) doi:10.1016/S0370-2693(02)01649-0 [arXiv:hep-ph/0204049 [hep-ph]]; S. K. Kang, Z. z. Xing and S. Zhou, Phys. Rev. D 73, 013001 (2006) doi:10.1103/PhysRevD.73.013001 [arXiv:hep-ph/0511157 [hep-ph]]; S. Chang, S. K. Kang and K. Siyeon, Phys. Lett. B 597, 78-88 (2004) doi:10.1016/j.physletb.2004.06.104 [arXiv:hep-ph/0404187 [hep-ph]].
(37) Y. H. Ahn, H. Y. Cheng and S. Oh, Phys. Rev. D 83, 076012 (2011) [arXiv:1102.0879 [hep-ph]].
(38) L. Wolfenstein, Phys. Rev. Lett. 51, 1945 (1983).
(39) Y. H. Ahn, H. Y. Cheng and S. Oh, Phys. Lett. B 703, 571 (2011) [arXiv:1106.0935 [hep-ph]].
(40) L. L. Chau and W. Y. Keung, Phys. Rev. Lett. 53, 1802 (1984).
(41) http://ckmfitter.in2p3.fr.
(42) R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022).
(43) R. D. Bolton et al., Phys. Rev. D 38, 2077 (1988).
(44) A. V. Artamonov et al. [E949 Collaboration], Phys. Rev. Lett. 101, 191802 (2008) [arXiv:0808.2459 [hep-ex]].
(45) A. Davidson and K. C. Wali, Phys. Rev. Lett. 48 (1982), 11; F. Wilczek, Phys. Rev. Lett. 49, 1549 (1982).
(46) J. L. Feng, T. Moroi, H. Murayama and E. Schnapka, Phys. Rev. D 57, 5875 (1998) [hep-ph/9709411].
(47) Z. G. Berezhiani and M. Y. Khlopov, Sov. J. Nucl. Phys. 51 (1990), 739-746; Z. G. Berezhiani and M. Y. Khlopov, Z. Phys. C 49 (1991), 73-78; A. S. Sakharov and M. Y. Khlopov, Phys. Atom. Nucl. 57 (1994), 651-658
(48) E. Cortina Gil et al. [NA62], JHEP 06 (2021), 093 [arXiv:2103.15389 [hep-ex]].
(49) A. Ayala, I. DomÃnguez, M. Giannotti, A. Mirizzi and O. Straniero, Phys. Rev. Lett. 113, no. 19, 191302 (2014) [arXiv:1406.6053 [astro-ph.SR]].
(50) J. Redondo, JCAP 12 (2013), 008 [arXiv:1310.0823 [hep-ph]].
(51) N. Viaux, M. Catelan, P. B. Stetson, G. Raffelt, J. Redondo, A. A. R. Valcarce and A. Weiss, Astron. Astrophys. 558 (2013), A12 [arXiv:1308.4627 [astro-ph.SR]]; N. Viaux, M. Catelan, P. B. Stetson, G. Raffelt, J. Redondo, A. A. R. Valcarce and A. Weiss, Phys. Rev. Lett. 111 (2013), 231301 [arXiv:1311.1669 [astro-ph.SR]].
(52) M. M. Miller Bertolami, B. E. Melendez, L. G. Althaus and J. Isern, JCAP 1410, no. 10, 069 (2014) [arXiv:1406.7712 [hep-ph]].
(53) J. Isern, E. Garcia-Berro, S. Torres and S. Catalan, Astrophys. J. 682, L109 (2008) [arXiv:0806.2807 [astro-ph]]; J. Isern, S. Catalan, E. Garcia-Berro and S. Torres, J. Phys. Conf. Ser. 172, 012005 (2009) [arXiv:0812.3043 [astro-ph]]; J. Isern, M. Hernanz and E. Garcia-Berro, Astrophys. J. 392, L23 (1992).
(54) G. G. Raffelt, Phys. Lett. B 166, 402 (1986); S. I. Blinnikov and N. V. Dunina-Barkovskaya, Mon. Not. Roy. Astron. Soc. 266, 289 (1994).
(55) S. DeGennaro, T. von Hippel, D. E. Winget, S. O. Kepler, A. Nitta, D. Koester and L. Althaus, Astron. J. 135 (2008), 1-9 [arXiv:0709.2190 [astro-ph]].
(56) N. Rowell and N. Hambly, [arXiv:1102.3193 [astro-ph.GA]].
(57) L. Wolfenstein, Phys. Rev. D 17, 2369 (1978); S. P. Mikheyev and A. Y. Smirnov, Sov. J. Nucl. Phys. 42, 913 (1985) [Yad. Fiz. 42, 1441 (1985)].
(58) I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, JHEP 1901, 106 (2019) [arXiv:1811.05487 [hep-ph]].
(59) P. F. de Salas, D. V. Forero, C. A. Ternes, M. Tortola and J. W. F. Valle, Phys. Lett. B 782, 633 (2018) [arXiv:1708.01186 [hep-ph]].
(60) F. Capozzi, E. Lisi, A. Marrone and A. Palazzo, Prog. Part. Nucl. Phys. 102, 48 (2018) [arXiv:1804.09678 [hep-ph]].
(61) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, JHEP 09, 178 (2020) [arXiv:2007.14792 [hep-ph]]; NuFIT 5.2 (2022), www.nu-fit.org.
(62) S. Abe et al. [KamLAND-Zen], Phys. Rev. Lett. 130, no.5, 051801 (2023) [arXiv:2203.02139 [hep-ex]].
(63) S. A. Kharusi, G. Anton, I. Badhrees, P. S. Barbeau, D. Beck, V. Belov, T. Bhatta, M. Breidenbach, T. Brunner and G. F. Cao, et al. Phys. Rev. D 104, no.11, 112002 (2021) [arXiv:2109.01327 [hep-ex]].
(64) S. Antusch, J. Kersten, M. Lindner, M. Ratz and M. A. Schmidt, JHEP 0503, 024 (2005) [hep-ph/0501272].
(65) S. J. Asztalos et al. [ADMX Collaboration], Phys. Rev. Lett. 104, 041301 (2010) [arXiv:0910.5914 [astro-ph.CO]].
(66) B. Abi et al. [DUNE], [arXiv:1807.10334 [physics.ins-det]].
(67) K. Abe et al. [Hyper-Kamiokande], [arXiv:1805.04163 [physics.ins-det]].
(68) C. Adams et al. [NEXT], JHEP 2021 (2021) no.08, 164 [arXiv:2005.06467 [physics.ins-det]].
(69) A. Allega et al. [SNO+], Phys. Rev. D 105 (2022) no.11, 112012 [arXiv:2205.06400 [hep-ex]].
(70) M. Askins et al. [Theia], Eur. Phys. J. C 80 (2020) no.5, 416 [arXiv:1911.03501 [physics.ins-det]].
(71) R. Arnold, C. Augier, A. M. Bakalyarov, J. Baker, A. Barabash, P. Bernaudin, M. Bouchel, V. Brudanin, A. J. Caffrey and J. Cailleret, et al. Nucl. Instrum. Meth. A 536 (2005), 79-122 [arXiv:physics/0402115 [physics]].
(72) N. Aghanim et al. [Planck], Astron. Astrophys. 641 (2020), A6 [erratum: Astron. Astrophys. 652 (2021), C4] [arXiv:1807.06209 [astro-ph.CO]].

Simple Modular invariant model for Quark, Lepton, and flavored-QCD axion