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

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{\"a}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 [1] 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 [2] 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 [3] proposed a new idea regarding the origin of the structure of lepton mixing. He applied modular invariance 1 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 [5], the Yukawa coefficients are assumed to be free parameters 2 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. [7] 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.
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 1 Modular invariance was analyzed for supersymmetric field theories in Ref. [4]. 2 In Ref. [6], the Yukawa coefficients are assumed to be of order unity. derived from superstring theory with G SM ×Γ N ×U (1) X symmetry. The non-Abelian discrete symmetry Γ 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 [8]. To ensure the validity of a modular invariant model with G SM × Γ N × 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. [9]). 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 τ , 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. [10]) associated with the modular transformation Eq.(1) under the non-local modular group Γ 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 [11], 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 [12]. The PQ breaking scale, or the low axion decay constant, can be determined by taking into account both supersymmetry (SUSY)-breaking effects [13] and supersymmetric next-leading-order Planck-suppressed terms [14][15][16].
The model features a minimal set of fields that transform based on representations of G SM × Γ N × U (1) X , and includes modular forms of level N . These modular forms act as Yukawa couplings and transform under the modular group Γ 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 τ -independent superpotential for scalar potential.
The so-called flavored-PQ symmetry U (1) X guarantees the absence of bare mass terms [17].
We minimally extend the model by incorporating three right-handed neutrinos N c and SM gauge singlet scalar fields χ(χ). 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 F = χ(χ) with modular weight zero acts on dimension-four(three) operators well-sewed by G SM × Γ N × U (1) X and modular invariance with different orders, which generate the effective interactions for the SM and the right-handed neutrinos as follows, Here, Λ is the scale of flavor dynamics above which unknown physics exists as a UV cutoff, and Yukawa coefficients c n (c 1 ) are all complex numbers assumed to have a unit absolute value (|c 1 |, |c n | = 1). The dimension-4(3) operators O 4(3) are determined by G SM ×Γ N ×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 Γ N [3]. 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 χ(χ) fields being zero. Note that there exist the infinite series of higher-dimensional operators induced solely by the combination of χχ 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 coefficientsc 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 [18], we imposed a U (1) X -mixed gravitational anomaly-free condition [16,19,20].
The rest of this paper is organized as follows: The next section discusses modular and U (1) X anomaly-free conditions under G SM × Γ N × 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 lowenergy effective field theory derived from superstring theory (see Ref. [9] for a review). In particular, type IIA intersecting D-brane models are related to magnetized D-brane models through T-duality [9]. The group Γ(N ) acts on the complex variable τ , varying in the upperhalf complex plane Im(τ ) > 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 Γ(N ), called principal congruence subgroups of level N = 1, 2, 3, ..., are defined by which are normal subgroups of homogeneous modular group Γ ≡ Γ(1) ≃ SL(2, Z), where SL(2, Z) is the group of 2×2 matrices with integer entries and determinant equal to one. The projective principal congruence subgroups are defined asΓ(N ) = Γ(N )/{±I} for N = 1, 2.
For N ≥ 3 we haveΓ(N ) = Γ(N ) because the elements −I does not belong to Γ(N ). The modular groupΓ ≡ Γ/{±I} is generated by two elements S and T , satisfying They can be represented by the P SL(2, Z) matrices The groups Γ N are finite modular groups obtained by imposing the condition T N = 1 in addition to Eq.(5), where Γ N ≡Γ/Γ(N ). The groups Γ N are isomorphic to the permutation groups S 3 , A 4 , S 4 , and A 5 for N = 2, 3, 4, 5, respectively [8].
We work in the 4D N = 1 string-derived supergravity framework defined by a general Kähler potential G(Φ,Φ) of the chiral superfields Φ and their conjugates, and by an analytic gauge kinetic function f (Φ) of the chiral superfields Φ, where M P = (8πG N ) −1/2 = 2.436 × 10 18 GeV is the reduced Planck mass with Newton's gravitational is a real gauge-invariant function of Φ andΦ, and W (Φ) is a holomorphic gauge-invariant function of Φ. 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 SM × Γ N × U (1) X . Here we assume that the non-Abelian discrete symmetry Γ N as a finite subgroup of the modular group [3], and the anomalous gauged U (1) X including the SM gauge symmetry G SM may arise from several stacks on D-brane models [9]. In the 4D global supersymmetry the most general form of the action can be written as where Then the given symmetry G SM × Γ N × 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 Γ N , and (ii) anomalous triangle graphs with external states A a ν A b ρ V Xµ where A a ν , A b ρ are gauge bosons of the SM gauge group G SM and V µ X is the connection associated with the gauged U (1) X . These anomalies can be cancelled by the GS mechanism [12].
A. Modular anomaly-free and modular forms of level N To demonstrate the invariance of K(Φ,Φe 2A ) and f (Φ)W α W α of Eq.(8) under the finite modular group Γ N and the gauged U (1) X , we consider a low-energy Kähler potential 4 : where −k is the modular weight, Z X is the nomalization factor, S denotes the axio-dilaton field, τ represents the overall Kähler modulus, and U X and U i correspond to the complex structure moduli. The dots in Eq.(10) denote the contributions of non-renomalizable terms scaled by an UV cutoff M P invariant under G SM × Γ N × U (1) X . We note that the matter fields φ 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 µ 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 U i , f ab (Φ) ⊃ δ ab (U X + U i ). These moduli are associated with the SM gauge theory, which we will not be focusing on. The GS parameter δ GS X characterizes the coupling of the 3 The 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). 4 It is similar to the one-loop Kähler potential presented in Ref. [10].
anomalous gauge boson to the axion θ 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 φ and φ X are neutral under the U (1) X symmetry and the modular group Γ N , respectively.
Calculating K IJ = ∂ I ∂ 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 ≫ ⟨φ⟩, ⟨φ X ⟩ and V X = 0 as follows, where K φφ = K φ XφX = 1 for canonically normalized scalar fields achieved by rescaling the fields φ and φ X for given values of the VEVs of τ and U X . The U (1) X charged modulus U X and scalar field φ X can be decomposed as: where g X is 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 φ X for v X → ∞.
Since the matter superfields φ and axio-dilaton S transform as, where ρ(γ) is the unitary representation of the modular group Γ N andc is a constant, the transformation of the Kähler potential K given in Eq.(9) leads us to Generically, the transformation of K in Eq.(9) incorporating Eq. (14) gives rise to a modular anomaly arising from δS = −c 1 4 d 4 xd 2 θW α W α g(τ ) + h.c. [10], whereQ µν = 1 2 ε µνρσ Q ρσ 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, 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(τ ), g(τ ).
The modular invariance W (Φ) under the modular group Γ N (N ≥ 2) provides a strong restriction on the flavor structure [3]. The superpotential W (Φ) can be expanded in power series of the multiplets φ which are separated into brane sectors φ (I) where the functions Y I 1 ...In (τ ) are generically 5 τ -dependent in type IIA intersecting Dbrane models [9,24]. The superpotential W (Φ) must have modular invariance under the , where g(τ ) is given by Eq. (14). To ensure this, we need to satisfy two conditions: (i) the matter superfields φ I i of the brane sector I i should in a representation ρ (I i ) (γ) of the modular group Γ N , where −k I i is the modular weight of sector I i , and (ii) the functions Y I 1 ...In (τ ) should be modular forms of weight k Y (n) transforming in the representation ρ(γ) of Γ N , with the requirements 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. 5 In Type II string orientifold compactifications, the Yukawa couplings have modular properties [25].
For example, a τ -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.
anomalies are given respectively by Here U (n) generators (n ≥ 2) are normalized according to Tr[T a T b ] = δ ab /2, and for leads to D-term potential for the anomalous U (1) X , where ξ 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 Re[U X ] = ρ X /2. Since the FI term is controlled by the string coupling it can not be zero. The re-stabilization of VEVs by φ 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 where the anomalous currents J X µ and J θ µ coupling to the gauge boson A µ X (that is, Expanding Eq.(24) and setting where the gauge boson mass m X obtained by the super-Higgs mechanism is given by Then the open string axion a X (decay constant f X ) is mixed linearly with the closed string a θ (decay constant f θ ): where the approximations are valid under the assumption that f θ is much larger than f X .
The gauged U (1) X absorbs one linear combination of a X and a θ , denoted G, giving it a string scale mass through the U (1) X gauge boson, while the other combination,Ã ≈ 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 Γ(3) modular symmetry, which gives the modular forms of level 3. The group Γ 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 1, 1 ′ , and 1 ′′ and one triplet 3 with the multiplication rules 3⊗3 = 3 s ⊕3 a ⊕1⊕1 ′ ⊕1 ′′ and 1 ′ ⊕1 ′ = 1 ′′ , 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 3's. Then we have The details of the A 4 group are shown in Appendix A. The modular forms f (τ ) of level 3 and weight k, such as Eq. (19), are holomorphic functions of the complex variable τ with well-defined transformation properties with an integer k ≥ 0, under the group Γ 3 . The three linearly independent weight 2 and level-3 modular forms are given by [3] where ω = −1/2 + i √ 3/2 and η(τ ) is the Dedekind eta-function defined by The Dedekind eta-function satisfies The three linear independent modular functions transform as a triplet of is constrained by the relation, A. Modular invariant supersymmetric potential and a Nambu-Goldstone mode Using Eqs. (17)(18)(19)(20), we construct unique supersymmetric and modular invariant scalar potential by introducing minimal supermultiplets. Those include SM singlet fields χ 0 6 with modular weight three and χ(χ) 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 χ andχ 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 [1] is implemented, the field χ orχ may be responsible for the heavy neutrino mass term).
Under k I × A 4 × 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, 1, 0) and three SM gauge singlets χ,χ, χ 0 to be (0, 1, +1), The A 4 -singlet χ 0 field with modular weight three ensures that the functions Y I 1 ...In (τ ) are independent of τ . Then the supersymmetric scalar potential invariant under G SM × U (1) X × A 4 is given at leading order by where dimensionless coupling constants g χ 0 and g χ are assumed to be equal to one, but are modified to Eq.(60) by considering all higher order terms induced by χχ combinations. Note that the PQ breaking parameter µ χ corresponds to the scale of the spontaneous symmetry breaking.
In the global SUSY limit, i.e. M P → ∞, the scalar potential obtained by the F -and Dterm 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 The scalar fields χ andχ have X-charges +1 and −1, respectively, i.e., with a constant ξ. So the potential V 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 → ∞ vanishes at its ground states, i.e. ⟨V global F ⟩ = 0 and ⟨V global D ⟩ = 0 vanishing F -term must have also vanishing D-term. From the minimization of the F -term scalar potential we obtain 7 As a consequence of k I × A 4 × U (1) X the other superpotential term κ α L α 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.
where we have assumed ⟨χ⟩, ⟨χ⟩ ≫ ⟨H u,d ⟩. The above supersymmetric solution is taken by the D-flatness condition for [15,16] ξ FI The tension between ⟨χ⟩ = ⟨χ⟩ and ξ FI X ̸ = 0 arises because the FI term cannot be cancelled, unless the VEV of flux in the FI term is below the string scale [13,31]. The FI term acts as an uplifting potential, where ∆ρ = ρ X − ρ 0 , which raises the Anti-de Sitter minimum to the de Sitter minimum [13]. To achieve this, the F -term must necessarily break SUSY for the D-term to act as an uplifting potential. The PQ scale µ χ can be determined by taking into account both the SUSY-breaking effect, which lifts up the flat direction, and supersymmetric nextleading-order Planck-suppressed terms [14][15][16]. The supersymmetric next-to-leading order where α 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 χ,χ at leading order read where m 3/2 represents soft SUSY-breaking mass, and α 1 , α 2 are real-valued constants. It leads to the PQ breaking scale (equivalently, the seesaw scale), indicating that µ χ lies within the range of approximately 1. 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 SUSY in the SUSY limit: where ψ denotes Dirac fermions, and V SUSY is replaced by V total when SUSY breaking effects are considered. The above kinetic terms for χ(χ) are canonically normalized from the Kähler potential Eq.(10). Here four component Majorana spinors (N c = N and ν c = ν) are used.
The global U (1) X PQ symmetry guarantees the absence of bare mass term in the Yukawa Lagrangian L Y in Eq. (46). The QCD Lagrangian has a CP-violating term where g s stands for the gauge coupling constant of SU (3) C , and G aµν is the color field strength tensor and its dualG a µν = 1 2 ε µνρσ G aµν (here a is an SU (3)-adjoint index), coming from the strong interaction. After obtaining VEV ⟨χ⟩ ̸ = 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.(71, 81, 92). To extract the associated boson resulting from spontaneous breaking of U (1) X , we set the decomposition of complex scalar fields [16,17,19] as follows in which A X is the NG mode and we set v χ = vχ and h χ = hχ in the supersymmetric limit.
The derivative coupling of NG boson A X arises from the kinetic term Performing u χ → ∞, the NG mode A X , whose interaction is determined by symmetry, is distinguished from the radial mode h χ , which is invariant under the symmetry U (1) X .

B. Modular invariant Yukawa superpotentials and Anomaly coefficients
By introducing just two A 4 -singlet fields, χ andχ, with modular weight zero and charged under U (1) X by +1 and −1, respectively, and using economic weight modular forms, we con- 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 χ(χ) 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 χχ are allowed, but they do not modify the leading-order flavor structure. Now let us assign A 4 × 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-I  Field where α  . Note that all Yukawa coefficients in the above suprpotential, α i , 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 1 ⊕ 3 ⊕ 3 under A 4 given explicitly by [3] 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 χ (orχ) are dependent on modular forms. Using modular forms of weight 6, Y  According to the quantum numbers of the quark sectors as in Table-I, the color anomaly Note that U (n) generators (n ≥ 2) are normalized according to Tr[T a T b ] = δ ab /2. The U (1) X is broken down to its discrete subgroup Z N DW in the backgrounds of QCD instanton, and the quantity N C (non-zero integer) is given by the axionic domain-wall number N DW . At the QCD phase transition, each axionic string becomes the edge to N DW domain-walls, and the process of axion radiation stops. To avoid the domain-wall problem one should consider N DW = 1 or the PQ phase transition occurred during (or before) inflation for N DW > 1.
Next, we turn to the lepton sector, where the fields are charged under G SM × A 4 × 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)×[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. In Table II, the representations and quantum and modular weight k I determined according to Eq. (20).
numbers of the lepton fields as well as modular weight k f determined along with Eq. (20) are presented. Here, L e , L µ , and L τ denote SU (2) L lepton doublets and e c , µ c , and τ 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 [1].
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.(78)). In contrast, modular forms Y 3 , Y 1 , and Y 3 are used to construct neutrino mass matrices 9 Then the Yukawa superpotential for lepton invariant under G SM × A 4 × U (1) X with economic modular forms are sewed with F = {χ orχ} through Eq.(2), respectively, as 9 By selecting appropriate modular weight of particle contents, lower weight modular forms can be used, such in the Dirac neutrino sector, and Y  lν stands for higher order contributions triggered by the combination χχ. Like in the quark sector, the Yukawa coefficients in the above superpotential, such as α i , 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 3 , respectively. Each term involves F/Λ 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 [1] by making the VEV ⟨χ⟩ large. The details on how the active neutrino masses and mixing are predicted will be presented in Sec. IV C.
Nonperturbative quantum gravitational anomaly effects [18] violate the conservation of the corresponding current, ∂ µ J µ X ∝ RR, where R is the Riemann tensor andR 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 × [gravity] 2 (related to the color anomaly U (1) X × [SU (3) C ] 2 ) must be cancelled, as shown in Refs. [16,19,20], which leads to the relation, 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 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 µν , Y µν , and F µν with their dual formsW µν ,Ỹ µν , andF µν , respectively. Here And the electromagnetic anomaly coefficient The physical quantities of QCD axion, such as axion mass m a and axion-photon coupling g aγγ , 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, 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 [32] and FLASH [33]( see Eq.(76) 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 [34] or nearly tribimaximal [35] 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 × 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 righthanded down-type quark fields, it is possible to construct down-type quark Yukawa superpotential with lower modular weight forms Y (1,1 ′′ ,1 ′ ) and U (1) X quantum numbers are taken to be (g l −f e − 1 2 , g l −f µ − 1 2 , g l −f τ − 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 − 1 2 . In this case, we have the freedom to select the weights. For instance, we can choose the following weights: k L = 5 2 , k e c = k µ c = k τ c = − 3 2 , and k N c = 1 2 . Then the lepton Yukawa superpotential reads where dots stand for higher order contributions triggered by the combination χχ. 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 m D ∝ (1, 1, 1), and the heavy Majorana neutrino mass term follows a similar form, as found in Ref. [3]. By selecting another specific weights, namely k L = 9 2 , k e c = k µ c = k τ c = − 7 2 , and k N c = − 3 2 , a notable change occurs in the modular form of the Majorana neutrino operator. Specifically, the term 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 [36]. (v) However, the assignment of the right-handed charged lepton E c = (e c , µ c , τ c ) to the A 4 triplet, the left-handed charged lepton denoted as L e,µ,τ to A 4 singlets (1,1 ′′ ,1 ′ ), 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 × U (1) X symmetries with modular invariance.
Non-zero VEVs of scalar fields χ(χ) spontaneously break the flavor symmetry U (1) X 10 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 10 If the symmetry U (1) X is broken spontaneously, the massless mode A X of the scalar χ appear as a phase.
parameter ⟨F⟩/Λ ≡ ∆ F . The higher order contributions of superpotentials W (h) q(lν) become ∞ n=1c i ∆ 2n χ ·(leading order operators) withc i = e iθ 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 α i , β i , γ i , we see that they are constrained as where the lower (upper) limit corresponds to the sum of higher order terms withθ i = π(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 where g is the SU where ψ f = {u, c, d, s, b, e, µ, τ } and α is a transformation constant parameter.

A. 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 ψ R M ψ V ψ † L =M ψ (diagonal form) and the mass eigenstates are ψ ′ R = V ψ R ψ R and ψ ′ L = V ψ L ψ L . These transformation include, in particular, the chiral transformation of Eq.(62) that necessarily makes M u,d,ℓ real and positive. This induces a contribution to the QCD vacuum angle in Eq.(46), i.e., with −π ≤ ϑ eff ≤ π. Then one obtains the vanishing QCD anomaly term where α ′ s = g 2 s /4π and the axion decay constant F a with f a = u χ of Eq. (48). At low energies A X will get a VEV, ⟨A X ⟩ = −F a ϑ eff , eliminating the constant ϑ eff term. The QCD axion then is the excitation of the A X field, a = A X − ⟨A X ⟩.
Substituting the VEV of Eq.(40) into the superpotential Eq.(50), the mass matrices M u and M d for up-and down-type quarks given in the Lagrangian (61) are derived as, where The terms with α d,s,b in Eq.(66) generate by taking the modular form Y Diagonalizing the matrices M † f M f and M f M † f (f = u, d) determine the mixing matrices V f L and V f R , respectively [37]. The left-handed quark mixing matrices V u L and V d L are components of the CKM matrix V CKM = V u L V d † L , which is generated from the down-type quark matrix in Eq.(66) due to the diagonal form of the up-type quark mass matrix in Eq. (65). The CKM matrix is parameterized by the Wolfenstein parametrization [38], see Eq.(B1), and has been determined with high precision [39]. The current best-fit values of the CKM mixing angles in the standard parameterization [40]   QED loops [42]. 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 [42].
After diagonalizing the mass matrices of Eqs. (65,66), the flavored-QCD axion to quark interactions are written at leading order as where V d † L = V CKM of Eq.(B1) is used by rotating the phases in M u away, which is the result of a direct interaction of the SM gauge singlet scalar field χ with the SM quarks charged under U (1) X . The flavored-QCD axion a is produced by flavor-changing neutral Yukawa interactions in Eq.(71), which leads to induced rare flavor-changing processes. The strongest bound on the QCD axion decay constant is from the flavor-changing process K + → π + + a [43][44][45][46], induced by the flavored-QCD axion a. From Eq.(71) the flavored-QCD axion interactions with the flavor violating coupling to the s-and d-quark are given by Then the decay width of K + → π + + a is given by where m K ± = 493.677 ± 0.013 MeV, m π ± = 139.57061 ± 0.00024 MeV [42]. From the present experimental upper bound Br(K + → π + a) < (3−6)×10 −11 (1×10 −11 ) for m a = 0−110 (160-260) MeV at 90% CL with Br(K + → π + νν) = (10.6 +4.0 −3.4 | stat ± 0.9 syst ) × 10 −11 at 68%CL [47], we obtain the lower limit on the QCD axion decay constant, The QCD axion mass m a in terms of the pion mass and pion decay constant reads [16,17] where f π ≃ 92.1 MeV [42] and F (z, w) = z/(1 + z)(1 + z + w) with ω = 0.315 z. Here the Weinberg value lies in z ≡ m MS u (2 GeV)/m MS d (2 GeV) = 0.47 +0.06 −0.07 [42]. After integrating out the heavy π 0 and η at low energies, there is an effective low energy Lagrangian with an axion-photon coupling g aγγ : L aγγ = −g aγγ a ⃗ E · ⃗ B where ⃗ E and ⃗ 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, The upper bound on the axion-photon coupling, derived from the recent analysis of the horizontal branch stars in galactic globular clusters [48], can be translated to where z = 0.47 is used.
Then the flavored-QCD axion to charged-lepton interactions read 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 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 [50] |g aee | < 4.3 × 10 −13 (95% CL) ⇔ N C F a ≳ 1.19|f e | × 10 9 GeV .
Bremsstrahlung off electrons e + Ze → 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 11 on |g aee | [53]. Recent analyses of WDLFs, using detailed WD cooling treatment and new data on the WDLF of the Galactic disk, suggest electron couplings |g aee | ≲ 2.8×10 −13 [51]. 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.

C. 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 where p = 1 + 2xy, r = y 2 + 2xy, γ = γ 2 /γ 1 , γ ′ = γ ′ 2 /γ 1 , and the common factor M can be replaced by the QCD axion decay constant F a , The terms with γ 2 in Eq.(84) are derived by taking the modular form Y 3,1 safisfying Eq.(51), whereas the terms with γ ′ 2 are derived by taking Y 3,2 . Eq.(84) has three unknown complex parameters, γ γ ′ , and γ 1 , where the phase of γ 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 ⟨χ⟩ is fixed from the seesaw formula Eq.(87) whose scale is given by PQ scale Eq. (45). The Dirac mass term in the Lagrangian (61) reads The coefficients β i and γ 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.(86) contains three complex parameters (β 1 , β 2 , and β 3 ), where one of the phases can be removable as an overall factor after seesawing. As shwon before, the parameter ∆ χ can be determined from quark and charged-lepton sectors.
In addition, its U (1) X quantum number g α 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 M ν is given at leading order by where U ν is the rotation matrix diagonalizing M ν and m ν i (i = 1, 2, 3) are the light neutrino masses. Then the PMNS mixing matrix becomes The matrix U PMNS is expressed in terms of three mixing angles, θ 12 , θ 13 , θ 23 , and a Dirac type CP violaitng phase δ CP and two additional CP violating phases φ 1,2 if light neutrinos are Majorana particle as [42] U where s ij ≡ sin θ ij , c ij ≡ cos θ ij and Q ν = Diag(e −iφ 1 /2 , e −iφ 2 /2 , 1). Then the neutrino masses are obtained by the transformation Here m ν i (i = 1, 2, 3) are the light neutrino masses. The observed hierarchy |∆m 2 for IO. have enabled a more precise determination of the mixing angles and mass squared differences, with large uncertainties remaining for θ 23 and δ CP at 3σ. The most recent analysis [60] lists global fit values and 3σ intervals for these parameters in Table-III at 90% confidence level. 0νββ 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 [1].
Transforming the neutrino fields by chiral rotations of Eq.(62) under U (1) X we obtain the flavored-QCD axion interactions to neutrinos where m ν i are real and positive, 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.(74), which significantly suppresses the interaction. Therefore, we will not take it into consideration. Ref. [62] provides the latest experimental constraints on Majoron-neutrino coupling, which are below the range of (0.4 − 0.9) × 10 −5 .   Table-IV: it implies that for normal mass ordering, the maximum scale should be below ∼ 10 15 GeV, and for inverted mass ordering, the maximum scale should be below ∼ 5 × 10 14 GeV. Refer to Table-V and -VI for NO, and   Table-VII for IO, with F a = 2|⟨χ⟩/N C |. By using the seesaw formula Eq.(87), one can set the scale ⟨χ⟩ along with the U (1) X quantum numbers g α without loss of generality. Doing so results in an effective mass matrix with nine physical degrees of freedom in which m 0 is an overall factor of Eq.(87), Out of the nine observables corresponding to Eq.(93), the five measured quantities (θ 12 , θ 23 , θ 13 , ∆m 2 Sol , ∆m 2 Atm ) can be used as constraints. The remaining four degrees of freedom correspond to four measurable quantities (δ CP , φ 1,2 , and the 0νββ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. (69,70) and Table-III), we use linear algebra tools from Ref. [63]. By analyzing experimental data for quarks and charged leptons, we determined the U (1) X quantum numbers listed in Tables-V and -VI for normal neutrino mass ordering and Table-VII 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.(93) 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 ⟨χ⟩ from roughly 10 13 GeV to 10 15 GeV due to Table-IV, corresponding to m 3/2 values of 1 TeV to 10 6 TeV, by considering Eq. (45) and Table-IV. The given U (1) X quantum numbers can then be used to predict the branching ratio of K + → π + + a (Eq.(73)), as well as the axion coupling to photon (Eq.(76)) and electron (Eq.(82)). See Table-V, -VI, and -VII for more details. Our proposed model makes predictions that can be tested by current axion search experiments like KLASH [32] and FLASH [33]. Fig.1 illustrates plots of the axion-photon coupling |g aγγ | 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-V, -VI, and -VII, which are consistent with the experimental constraints described in Eqs.(74),   (77), and (83). And Fig.1 illustrates that certain data points in Table-V (I-c, I-d, I-e) have been fully excluded by the ADMX experiment [64], while another data point (II) in Table-V 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).

A. Quark and charged-lepton
The Yukawa matrices for charged fermions in the SM, as given in Eqs. (65,66,78), 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 [42] and CKMfitter [41].
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 µ , m τ , and θ q 12 , θ q 23 , θ q 13 , δ q CP . These observables are used to determine thir-  Table-V and -VI  I  II  III  IV  V  VI  VII  VIII Table-VII  I  II  III  IV  Here, without loss of generality, the up-type quark masses m u , m c , and m t are a one-to-one correspondence with α u , α c , and α t , which have been taken real, and we have set arg(α b ) = 0.
The masses of the charged leptons m e , m µ , and m τ are in a one-to-one correspondence with the real parameters α e , α µ , and α τ from Eq.(79). Using the numerical results of Eq.(95) from the quark sector, with the inputs α e = 1.268 for |f e | = 20 (α e = 0.757 for |f e | = 19), α µ = 0.900, α τ = 1.148 , we obtain the charged lepton masses, which well agree with the empirical values of Eq.(80).

B. Neutrino
The seesaw mechanism in Eq.(87) 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 [63]. 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 δ CP and atmospheric mixing angle θ 23 . Using our model, we investigate which values of δ CP and θ 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 {θ 13 , θ 23 , θ 12 , ∆m 2 Sol , ∆m 2 Atm } at 3σ from Table-   we investigate how 0νββ-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σ presented in Ref. [60], we note that the most favored regions correspond to δ CP ∼ 250 • , whereas there are no favored regions with respect to θ 23 . Ongoing experiments like DUNE [65], as well as proposed next-generation experiments such as Hyper-K [66], are poised to greatly reduce uncertainties in the values of θ 23 and δ CP , providing a rigorous test for our proposed model. Furthermore, ongoing and future experiments on 0νββ-decay like NEXT [67], SNO+ [68], KamLAND-Zen [61], Theia [69], SuperNEMO [70] may soon reach a sensitivity to exclude the inverted mass ordering of our model.

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 SM ×Γ N ×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 SM × Γ N × 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 SM × Γ N × U (1) X , and includes modular forms of level N . These modular forms act as Yukawa couplings and transform under the modular group Γ 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 (χ andχ) 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 τ -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 χ(χ) have modular weights of zero, any additional correction terms arising from higher weight modular forms are not permitted in the superpotentials. However, the combination χχ 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 ⟨χ⟩ (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-V, -VI and -VII. 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 [32] and FLASH [33], see Fig.1 and 2, by considering the scale of U (1) X breakdown.