$SO(10)$ models with $A_4$ modular symmetry

We combine $SO(10)$ Grand Unified Theories (GUTs) with $A_4$ modular symmetry and present a comprehensive analysis of the resulting quark and lepton mass matrices for all the simplest cases. We focus on the case where the three fermion families in the 16 dimensional spinor representation form a triplet of $\Gamma_3\simeq A_4$, with a Higgs sector comprising a single Higgs multiplet $H$ in the ${\mathbf{10}}$ fundamental representation and one Higgs field $\overline{\Delta}$ in the ${\mathbf{\overline{126}}}$ for the minimal models, plus and one Higgs field $\Sigma$ in the ${\mathbf{120}}$ for the non-minimal models, all with specified modular weights. The neutrino masses are generated by the type-I and/or type II seesaw mechanisms and results are presented for each model following an intensive numerical analysis where we have optimized the free parameters of the models in order to match the experimental data. For the phenomenologically successful models, we present the best fit results in numerical tabular form as well as showing the most interesting graphical correlations between parameters, including leptonic CP phases and neutrinoless double beta decay, which have yet to be measured, leading to definite predictions for each of the models.


Introduction
Understanding the pattern of masses and mixings of quarks and leptons (including neutrinos) is one of the greatest challenges in particle physics. The quark and charged lepton masses are described by different interaction strengths with the Higgs doublets within the SM, but their values cannot be predicted. Much effort has been devoted to addressing the flavour puzzle, and symmetry has been an important guiding principle, where in particular non-Abelian discrete flavor symmetry is particularly suitable to predict the large lepton mixing angles [1]. However, in conventional flavour symmetry models, a complicated vacuum alignment is typically required, in which symmetry breaking Higgs fields (flavons) are introduced with their vacuum expectation values (VEVs) oriented along certain directions in flavour space.
Recently modular invariance, inspired by superstring theory with compactified extra dimensions, has been suggested to provide an origin of the flavour symmetry, especially as applied to the neutrino sector [2]. The finite discrete non-Abelian flavour symmetry groups emerge from the quotient group of the modular group SL(2, Z) over the principal congruence subgroups. The quark and lepton fields transform nontrivially under the finite modular groups and are assigned various modular weights. The modular invariance of Yukawa couplings, whose modular weights do not sum to zero, are necessarily modular forms which are holomorphic functions of the complex modulus τ . In such an approach, flavon fields other than the modulus may not be needed and the flavour symmetry can be entirely broken by the VEV of the complex modulus field τ . In principle, all higher dimensional operators in the superpotential are completely fixed by modular invariance. Modular forms are characterised by a positive integer level N and an integer weight k which can be arranged into modular multiplets of the homogeneous finite modular group Γ N ≡ Γ/Γ(N ) [3]. They are also organised into modular multiplets of the inhomogeneous finite modular group Γ N ≡ Γ/Γ(N ) if k is an even number [2]. The inhomogeneous finite modular group Γ N of the smaller levels N = 2 [4][5][6][7], N = 3 [2,4,5,, N = 4 [21,[34][35][36][37][38][39][40][41][42], N = 5 [39,43,44] and N = 7 [45] have been analysed and realistic models have been constructed. All the modular forms with integer weights can be generated from the tensor products of weight one modular forms and the odd weight modular forms which are in representations of Γ N with ρ r (S 2 ) = −1. The modular symmetry can not only help to explain the observed lepton mixing patter, but also possibly account for the hierarchical charged lepton masses which can arise from a small deviation from the modular symmetry fixed points [33,46,47].
These ideas have been generalised in several different directions. The homogeneous finite modular groups Γ N provide an even richer structure of modular forms applicable to flavor model building, with the smaller groups Γ 3 ∼ = T [3,48], Γ 4 ∼ = S 4 [49,50], Γ 5 ∼ = A 5 [51,52] and Γ 6 ∼ = S 3 × T [53] having been studied. If the modular weight k of the operator is not an integer, (cτ + d) k is not the automorphy factor anymore and some multiplier is needed, consequently the modular group should be extended to its metaplectic covering [54]. In this regard, the formalism of modular invariance has been extended to include the modular forms of rational weights k/2 [54] and k/5 [52]. Furthermore, generalized CP symmetry can be consistently imposed in the context of modular symmetry, where the modulus is determined to transform as τ → −τ * under the action of CP [55][56][57][58]. The CP transformation matrix is completely fixed by the consistency condition up to an overall phase, it has been shown to coincide with canonical CP transformation in the symmetric basis so that all the coupling constants would be real if the Clebsch-Gordan coefficients in the symmetric basis are real [55,58]. Finally, the underlying string theory typically involves a compact space with more than one modulus parametrizing its shape, and motivated by this, the modular invariance approach has been extended to incorporate several factorizable [36] and non-factorizable moduli [59]. Grand unified theories (GUTs) are amongst the most well motivated theories beyond the SM, realising the elegant aspiration to unify the three gauge interactions of the SM into a simple gauge group [60]. In GUTs, the quark and lepton fields are embedded into one or few gauge multiplets, resulting in quark and lepton mass matrices whose masses and mixing parameters are related. There is good motivation for imposing a family symmetry together with GUTs in order to address the problem of quark and lepton masses and mixing and especially large lepton mixing [61]. Among the many possible choices of family symmetry, A 4 is the minimal choice which admits triplet representations [62], while SU (5) [60] is the minimal GUT choice, being the smallest simple group which can accommodate the SM gauge symmetry. However, combining A 4 family symmetry with SU (5) GUTs [63], for example, also requires vacuum alignment of the flavons in order to break the A 4 , and so there is a strong motivation for introducing modular symmetry in such frameworks. Indeed modular symmetry has been combined with SU (5) GUTs in an (Γ 3 A 4 ) × SU (5) model in [10,64]. Other modular SU (5) GUT models have also been subsequently constructed based on (Γ 2 S 3 ) × SU (5) [6,65], and (Γ 4 S 4 ) × SU (5) [66,67].
While SU (5) GUTs is the most minimal choice, it is well known that it does not require neutrino masses, although they can easily be accommodated as singlet representations of the GUT group. This motivates the study of the SO(10) GUT group where all known chiral fermions of one generation plus one additional right-handed neutrino fit into a single 16 dimensional spinor representation of SO(10) [68], making neutrino mass inevitable. Just as in SU (5), so in SO(10) GUTs, there is good motivation for introducing a family symmetry [61,69]. However, just as before, the conventional approach of combining a family symmetry such as A 4 with SO(10) GUTs [70] also involves the introduction of flavons and vacuum alignment, leading to a complication which can be removed by the use of modular symmetry. It is perhaps surprising, therefore, that even though there is good motivation for studying modular symmetry as applied to SO(10) GUTs, there seems be no literature on this so far. The purpose of this paper is to remedy this deficiency.
In this paper, we shall perform a comprehensive study of the Γ 3 A 4 modular symmetry in the framework of supersymmetric (SUSY) SO(10) GUTs. The SO(10) gauge symmetry is spontaneously broken down to the SM gauge group SU (3) c ×SU (2)×U (1) Y by the VEV of Higgs fields ∆ in the 126 and/or Σ in the 120. In the minimal SO(10) model, one Higgs multiplet H in the 10 fundamental representation is required to further break the SM gauge symmetry into SU (3) c × U (1) EM , and we shall also assume such a Higgs. The neutrino masses are generated by the type-I and/or type II seesaw mechanisms [71][72][73][74][75] arising from the SU (2) L singlet and/or triplet components of the Higgs ∆ in the 126 representation. The purpose of this work is to find the simplest phenomenologically viable SO(10) GUT models based on A 4 modular symmetry. Models we consider are simply defined by organising the three families of fermions in the 16 dimensional spinor representation into a triplet of Γ 3 A 4 with a specified modular weight, plus a single Higgs multiplet H in the 10 fundamental representation, whose modular weight can be taken to be zero without loss of generality, supplemented by one Higgs field ∆ in the 126 and/or one Higgs field Σ in the 120, with specified modular weights. Once these modular weights are specified, the Yukawa couplings are then determined up to a number of overall dimensionless complex coefficients, and the value of the single complex modulus field τ , which is the only flavon in the theory. All models with sums of modular weights up to 10 are considered, and results are presented for each model following an intensive numerical analysis where we have optimized the free parameters of the models in order to match the experimental data. We present results only for the simplest phenomenologically viable models. We find that the minimal models containing only the Higgs fields in the 10 and the 126 are the least viable, requiring both type I and type II seesaw mechanisms to be present simultaneously and also necessitating at least some sums of modular weights of 10 (if the sums of modular weights were restricted to 8 then no viable such models were found). On the other hand, non-minimal models involving in addition to the Higgs fields in the 10 and the 126, also a Higgs field in the 120, proved to be more successful, with many such models being found with sums of modular weights of up to 8 or less, and with the type I seesaw as well as the combined type I+II seesaw also being viable. For the phenomenologically successful models, we present the best fit results in numerical tabular form as well as showing the most interesting graphical correlations between parameters, including leptonic CP phases and the effective mass in neutrinoless double beta decay, which have yet to be measured, leading to definite predictions for each of the models. This paper is organized as follows. In section 2 we briefly review modular symmetry and modular forms of level N = 3. In section 3 we discuss fermion masses from modular forms in SO (10) GUTs, based on the simplest models with three families of fermions in the 16 dimensional spinor representation in a triplet of Γ 3 A 4 , plus a single Higgs multiplet H in the 10 fundamental representation and one Higgs field ∆ in the 126 for the minimal models, and plus one Higgs field Σ in the 120 for the non-minimal models, with specified modular weights. In section 4 we present our numerical analysis for the phenomenologically successful models, presenting the best fit results in numerical tabular form as well as showing the most interesting graphical correlations between parameters. Section 5 concludes the paper.

Modular flavor symmetry
In modular-invariant supersymmetric theories, the action is generally of the form where K(τ,τ , Φ I ,Φ I ) is the Kähler potential, and W(τ, Φ I ) is the superpotential. The supermultiplets are divided into several sectors Φ I . The action is required to be modular invariant. The modular group Γ acts on the complex modulus τ as linear fraction transformation, where a, b, c and d are integers satisfying ad − bc = 1. The modular group Γ has an infinite number of elements and it can be generated by two generators S and T : which obey the relations S 2 = (ST ) 3 = 1. Under the action of Γ, the supermultiplets Φ I transform as [2] where −k I is the modular weight of Φ I , and ρ I is a unitary representation of the finite modular group Γ N = Γ/Γ(N ). Γ(N ) is the principal congruence subgroups of the modular group, and the level N is fixed in a concrete model. The Kähler potential is taken to the minimal form, and it gives rise to kinetic terms of Φ I and τ [2]. The superpotential W(τ, Φ I ) is strongly constrained by modular invariance, and it can be expanded in power series of Φ I as follows Modular invariance of W requires that Y I 1 ...In (τ ) should be a modular forms of weight k Y and level N transforming in the representation ρ Y of Γ N , i.e., The modular weights and the representations should satisfy the conditions where 1 denotes the invariant singlet of Γ N .

Modular group and modular forms of level 3
In the present work, we intend to impose modular flavor symmetry in SO(10) GUT to explain the flavor structure of quarks and leptons. We are concerned with the finite modular group Γ 3 which is isomorphic to A 4 . Γ 3 ∼ = A 4 is the even permutation group of four objects, and it can be generated by two generators S and T obeying the previous relations plus the additional one T 3 = 1, The A 4 group has four inequivalent irreducible representations including three one-dimensional representations 1, 1 , 1 and a three-dimensional representation 3. We utilize the symmetric basis in which the generators S and T are represented by unitary and symmetric matrices, i.e., 1 : S = 1, T = 1 , The multiplication rules for the tensor products of the A 4 representations are where 3 S and 3 A denote symmetric and antisymmetric contractions respectively. In terms of the components of the two triplets α = (α 1 , α 2 , α 3 ) T and β = (β 1 , β 2 , β 3 ) T , in this working basis we have The modular forms of level N and weight k span a linear space M k (Γ(N )). The space M k (Γ(3)) has dimension k + 1, and it can be expressed in terms of the Dedekind eta functions as follows [3,76] M k (Γ(3)) = a+b=k, a,b≥0 where the η(τ ) function is The three linearly independent weight 2 modular forms of level 3 can be arranged into a triplet of A 4 [3]: with Obviously we see that the constraint Y 2 2 + 2Y 1 Y 3 = 0 is fulfilled. Notice that ϑ(τ ) and ε(τ ) span the linear space of weight 1 modular forms of level 3 [3]. The expressions of the q−expansion of Y 1,2,3 (τ ) are given by Y 1 (τ ) = 1 + 12q + 36q 2 + 12q 3 + 84q 4 + 72q 5 + 36q 6 + 96q 7 + 180q 8 + 12q 9 + 216q 10 + . . . , Y 2 (τ ) = −6q 1/3 1 + 7q + 8q 2 + 18q 3 + 14q 4 + 31q 5 + 20q 6 + 36q 7 + 31q 8 + 56q 9 + 32q 10 + . . . , Y 3 (τ ) = −18q 2/3 1 + 2q + 5q 2 + 4q 3 + 8q 4 + 6q 5 + 14q 6 + 8q 7 + 14q 8 + 10q 9 + 21q 10 + . . . . (16) It agrees with the q−expansion derived in [2], where the modular forms Y 1,2,3 (τ ) are constructed in terms of η(τ ) and its derivative. The whole ring of even weight modular forms can be generated by the modular forms Y 1,2,3 (τ ) of weight 2. At weight 4, the tensor product of Y gives rise to three independent modular multiplets, There are seven linearly independent modular forms of level 3 and weight 6 and they decompose as The weight 8 modular forms can be arranged into three singlets 1, 1 , 1 and two triplets 3 of A 4 , The weight 10 modular forms of level 3 decompose as 1 ⊕ 1 ⊕ 3 ⊕ 3 ⊕ 3 under A 4 , and they are We summarize the even weight modular forms of level 3 and their decomposition under A 4 in table 1.
3II , Y (10) 3III Table 1: Modular forms of level 3 up to weight 10, where the superscripts indicate the modular weights and the subscripts denote how they transform under the A4 modular symmetry. Here Y 3I and Y 3II stand for the two linearly independent weight 6 modular forms in the triplet 3 of A4. Similar notations are adopted for Y 3III .

Fermion masses from modular forms in SO(10) GUTs
The SO(10) GUT theory embeds all SM fermions of a generation plus a right-handed neutrino into a single spinor representation denoted as 16. As a consequence, the gauge sector and the fermionic matter sector are generally quite simple. However, the same is not true of the Higgs sector. Since the larger GUT symmetry group SO(10) needs to be broken down to the Standard Model gauge group SU (3) C × SU (2) L × U (1) Y , generally one needs to introduce a large number of Higgs multiplets, with different symmetry properties under gauge transformations. There are two options when constructing SO(10) GUT models, one can use Higgs multiplets in low-dimensional representations of SO(10) then the Lagrangian contains non-renormalizable terms, or the Higgs multiplet fields are in high-dimensional representations and the Lagrangian is renormalizable. In this work, we shall stick to the renormalizable Supersymmetric (SUSY) SO(10) models for simplicity. The fermion masses are generated by the Yukawa couplings of fermion bilinears in the spinor representation 16 with the Higgs fields multiplets. From the following tensor product of SO(10) group where the subscripts S and A stand for the symmetric and antisymmetric parts of the tensor products respectively in flavor space, we see that the Higgs fields in the SO(10) representations 10, 120 and 126 can have renormalizable Yukawa couplings. Hence the Yukawa superpotential in renormalizable SO(10) models can be generally written as follows: where a, b = 1, 2, 3 are indices of generation, ψ refers to the matter fields in the 16 dimensional representation of the SO(10), H, ∆ and Σ denote the Higgs fields in the representations 10, 126 and 120 respectively. The Yukawa coupling matrices Y 10 and Y 126 are 3 × 3 complex symmetric in flavor space, while Y 120 is a complex antisymmetric matrix, i.e.
After the SO ( where the components (15, 2, 2) and (1, 2, 2) both contain a pair of the Y = ±1 SU (2) L Higgs doublets, whose neutral components give masses to the fermions. The component (10, 1, 3) contained in 126 gives Majorana masses of the right-handed neutrinos and the Majorana masses of the lefthanded neutrinos are generated due to the (10, 3, 1) component of 126. To be more specific, the SO(10) Higgs fields H, ∆ and Σ have SM Higgs components, and the decomposition of theses fields under the SM gauge symmetry SU Decomposing the Yukawa superpotential in Eq. (22), we find the quark and lepton mass matrices are of the following form where the VEVs are defined as Therefore both type-I and type-II seesaw mechanisms contribute to the neutrino masses, and the effective mass matrix of light neutrinos is given by where the first and the second terms denote the type-II and a type-I seesaw contributions respectively, the two contributions to neutrino mass depend on two different parameters v L and v R . It is possible to have a symmetry breaking pattern in SO (10) such that the first contribution (the type-II term) dominates over the type-I term. It is convenient to redefine the parameters as follows where v u and v d are the VEVs of the MSSM Higgs pair H u and H d . The parameters r a (a = 1, 2, 3 and c b (b = e, ν) are the mixing parameters which relate the H u,d to the doublets in the various GUT multiplets [77,78]. Notice that Y 10 , Y 126 and Y 120 are proportional to the Yukawa matrices u vu can be absorbed into the coupling constants α i , β i and γ i given in Eq. (34). Then the mass matrices of the quarks and leptons can be written as The effective neutrino mass matrix is still given by Eq. (31).

Combining SO(10) with A 4 modular symmetry
The three generation of fermions ψ 1,2,3 are assumed to transform as a triplet 3 under A 4 modular symmetry, and its modular weight is denoted as k F . All the Higgs multiplets H, ∆ and Σ are assigned to trivial A 4 singlet with modular weights k 10 , k 126 and k 120 respectively. Without loss of generality we can set k 10 = 0, since k 10 can be absorbed by shifting the modular weight of the matter fields. Modular invariance requires the Yukawa couplings Y 10 ab , Y 126 ab and Y 120 ab in Eq. (22) are just modular forms of level 3. The most general Yukawa superpotential invariant under both SO(10) and modular symmetry is of the following form where the representations r a,b,c and r a,b,c fulfill r a ⊗r a = r b ⊗r b = r c ⊗r c = 1, and the allowed values of r a,b,c depend on the modular weights of matter fields and Higgs multiplets, as summarized in table 1.
Note that for simplicity we have written α a ≡ α . Using the group theory contraction rules of A 4 in Eq. (11), the Yukawa matrices can be straightforwardly derived in terms of the complex coefficients α a , β b , γ c and the modular forms Y where S (k) r are symmetric 3 × 3 matrices defined below. The Yukawa matrix Y 126 is of a similar form with α i replaced by γ i . By comparison Y 120 is antisymmetric, it only receives contribution from the antisymmetric contraction (ψψ) 3 A , and its concrete forms crucially depends on the modular weight k = 2k F + k 120 with the antisymmetric matrices A (k) r : Parameters 0.02948 ± 0.00087 Table 2: The best fit values µi and 1σ uncertainties of the quark and lepton parameters when evolve to the GUT scale as calculated in [79], with the SUSY breaking scale MSUSY = 500 GeV and tan β = 10, where the error widths represent 1σ intervals. The parameter r ≡ ∆m 2 21 /∆m 2 31 is the ratio of neutrino mass-squared differences. The values of lepton mixing angles, leptonic Dirac CP violation phases δ l CP and the neutrino mass squared difference are taken from NuFIT 5.0 [80].
For convenience we have defined the symmetric and antisymmetric 3 × 3 matrices and which depend on the modular forms Y 3I,II and S (10) 3I,II,III can be read out from Eq. (38) with the components of triplet modular forms in Eqs. (18,19,20), and similarly for A (6) 3I,II , A 3I,II and A (10) 3I,II,III .

Minimal models
If only one of the Higgs fields H, ∆, Σ is employed, all fermion mass matrices are proportional to the same Yukawa matrix which can be diagonalized by a basis transformation. As a consequence, there would be no flavor mixing between up type and down type quarks or between charged leptons and neutrinos. Hence at least two Higgs fields are necessary for realistic fermion spectrum and flavor mixing. One of the Higgs fields must be a ∆ in the SO(10) representation 126 in order to generate tiny neutrino masses. The Yukawa couplings Y 126 ab of ∆ to the fermions provide Majorana masses to both right-handed and left-handed neutrinos and the seesaw mechanism is naturally realized, as can be seen from Eq. (33). In the so-called minimal SO(10) GUT, only the Higgs fields H and ∆ couple to the fermion sector. It is remarkable that the Yukawa superpotential is completely fixed by the summation of the modular weights of the matter and Higgs fields. In the following we present three benchmark models.
Then we have k F = 5 and k 126 = −4. Modular invariance requires that the modular forms of weight 6 and weight 10 should be present in the Yukawa interactions, and the most general form of the superpotential reads as • Minimal model 2: (2k F + k 10 , 2k F + k 126 ) = (10,8) For k 10 = 0, we have k F = 5 and k 126 = −2. Different from minimal model 1, the weight 8 modular forms of level 3 enter into the Yukawa couplings of ∆, and the superpotenital is given by • Minimal model 3: (2k F + k 10 , 2k F + k 126 ) = (10, 10) Both H and ∆ couple with the weight 10 modular forms, and the modular weights of fields can be taken as k F = 5, k 10 = k 126 = 0. Thus we can read off the superpotential as follows, For each term of each model the explicit Yukawa matrices can be constructed in terms of the complex coefficients α a , γ c and the modular forms Y (k) r (τ ) as described in the previous subsection. The total mass matrices will also depend on the mixing parameters r 1 , r 2 and the VEVs v u , v d plus the neutrino parameters as in Eq. (33). The above three minimal SO(10) models with A 4 modular symmetry will be confronted with the experimental data, and the results are listed in table 3. We see that the experimental data of quark and lepton masses and mixing parameters can be accommodated except the minimal model 1, the details would be discussed in section 4.

Non-minimal models
It is known that the minimal model is highly constrained in explaining the fermion masses and mixings. The Higgs field Σ which is a 120 dimensional multiplet of SO(10) make it easy to account for the different masses and mixing patterns of the quarks and leptons, because the antisymmetric Yukawa coupling matrix Y 120 have different coefficients in the Dirac mass matrices of the up-type quarks, down-type quarks, charged leptons and neutrinos, as shown in Eq. (33). As a consequence, the experimental data can be accommodated by using less terms than the minimal SO(10) models. The superpotential is strongly constrained by modular symmetry and it is fixed by the modular weights of matter fields and Higgs fields. Since only the combination of the modular weights of matter fields and Higgs fields is relevant, the modular weight k 10 can be taken to be vanishing without loss of generality. For illustration, we present three typical models with small number of free parameters in the following.
• Non-minimal model 2: (2k F + k 10 , 2k F + k 120 , 2k F + k 126 ) = (4, 2, 4) The modular weights are determined to be k F = 2, k 120 = −2, k 126 = 0 for k 10 = 0. Modular invariance requires that the weight 2 and weight 4 modular forms of level 3 should be involved in the Yukawa coupling, and the superpotenital is given by • Non-minimal model 3: (2k F + k 10 , 2k F + k 120 , 2k F + k 126 ) = (4, 6, 4) This model differs from the non-minimal model 2 in the modular weight of Σ, and the Yukawa superpotential reads as Similar to the previous minimal models, for each term of each model the explicit Yukawa matrices can be constructed in terms of the complex coefficients α a , β b , γ c and the modular forms Y (k) r (τ ). The total mass matrices will also depend on the mixing parameters r 1 , r 2 , r 3 and the VEVs v u , v d plus the neutrino parameters as in Eq. (33). A comprehensive numerical analysis is performed for the above three models, as described in section 4. It can be seen that excellent agreement with experimental data can be achieved for certain values of the free parameters.

Numerical analysis
A priori we cannot know which type of seesaw dominates or if they are of the same order of magnitude. In this section, we shall perform a detailed analysis of the predictions of the benchmark SO(10) models given in section 3, and we shall consider the type I seesaw dominant scenario, type II seesaw dominant scenario and the mixture of both seesaw mechanisms in the neutrino sector. For any given values of the coupling constants and the complex modulus τ , one can numerically diagonalize the fermion masses and then the mass eigenvalues, mixing angles and CP violation phases of both quarks and leptons can be extracted. The variation of the model parameters will dynamically affect the values of these experimental observables. In order to quantitatively estimate how well the above SO(10) modular models can reproduce the fermion masses and mixing patterns at the GUT scale, we perform a χ 2 analysis to find out the best fit values of the free parameters and the corresponding predictions for fermion masses and mixing parameters. Our full set of observables to which the models are fitted are the masses of quarks and charged leptons, the solar and atmospheric neutrino mass-squared differences, mixing angles of quarks and leptons, the CP violation phases δ q CP and δ l CP in the CKM matrix and lepton mixing matrix. The χ 2 function is defined as where P i is the theoretical predictions as functions of the parameter set x = {τ, α i , β i , γ i , r 1 , r 2 , r 3 , c e , c ν , v L , v R }, µ i are the experimental values extrapolated to the GUT scale, and σ i denote the 1σ errors in µ i . Note that we have redefined the coupling constants α i , β i and γ i to absorb the coefficients v 10 u vu respectively. The χ 2 sum runs over the seven mass ratios between charged fermions, ratio of the solar to atmospheric mass squared differences, the mixing angles and Dirac CP phases of quarks and leptons. We collect the values of µ i and σ i underlying our analysis in table 2, where the charged fermion masses and the quark mixing parameters measured at the electroweak scale have been evolved up to the corresponding ones at the GUT scale. The normal ordering neutrino masses are favored at 2.7σ level if the atmospheric neutrino data of Super-Kamiokande is taken into account [80]. Therefore we assume normal ordering spectrum, the neutrino masses and mixing that we use are the low scale values taken from NuFIT 5.0 [80]. and we have ignored the effects of the evolution from the low energy scale to the GUT scale. The effects of evolution on the neutrino mass ratios and on the mixing angles are known to be negligible to a good approximation for the normal hierarchical spectrum.
In order to fit the model parameters to the observables, we numerically minimize χ 2 with respect to free parameter vector x. We use the numerical minimization algorithm TMinuit [81] developed by CERN to numerically minimize the χ 2 function to determine the best fit values of the input parameters. The mechanism of moduli stabilization is still an open question at present, consequently we treat the modulus τ as a free parameter to adjust the agreement with the data, and it randomly varies in the fundamental domain {τ |Reτ | ≤ 1/2, Imτ > 0, |τ | ≥ 1}. After performing minimization of the χ 2 -function, we evaluate the fermion masses and mixing parameters with the best fit values of the free parameters. The overall scale α 1 v u of the up quark mass matrix is fixed by the top quark mass. The down quarks and charged lepton matrices share the same factor α 1 r 1 v d which is fixed by the measured value of the bottom quark mass. The scale factor of the light neutrino mass matrix is α 1 v 2 u /v R for type I seesaw and α 1 v L for type II seesaw and its value is determined by the solar neutrino mass splitting ∆m 2 21 . Because the parameter space is high dimensional and the χ 2 is a non-linear and complex function of the free parameters, generally many local minima exist. It is impossible to numerically determine whether a minimum is the global minimum of the χ 2 function under consideration. We have repeated the minimization procedure many times with different initial values of the free parameters, and then choose the lowest one out of the many local minima found by the program. However, it is difficult to rule out the existence of still lower minima and predictions may be improved if they exist.
In order to quantitatively measure the degree of fine-tuning needed in the models, we use the fine-tuning factor d F T which was firstly introduced in [82]. The parameter d F T is defined as where par i is the best fit value of i th real input parameter and err i is the corresponding offset of par i which changes χ 2 by one unit with all other parameters fixed at their best fit values. Notice that err i is not the uncertainty of the experimental data given in table 2. If some |err i /par i | are very small, then a small variation of the corresponding parameters would make a large difference on the χ 2 . Hence d F T can roughly measure the amount of fine-tuning involved in the fit. Similar to Ref. [82], one can understand the significance of the fine-tuning by comparing d F T with a similar parameter d Data derived from the data: where µ i and σ i are the central values and the 1σ of the i th observables respectively. We have d Data = 531.986 for the set of the experimental data in table 2. As shown in Eq. (31), generally both type I and type II seesaw can contribute to the light neutrino masses in SO(10) GUTs. A priori, we do not know which type of seesaw dominates or if they are of the same order of magnitude. In the following, we shall perform a detailed analysis of the predictions of the benchmark SO(10) models given in section 3, and we shall consider three scenarios: the type I seesaw dominance, the type II seesaw dominance and the mixture of both seesaw mechanisms in the neutrino sector.

Numerical results of the minimal models
In this section, we consider the numerical results of the models with the minimal Higgs content H and ∆ which are in the SO(10) representations 10 and 126 respectively. As a result, the term Y 120 is absent in the fermion mass matrices of Eq. (33). From Eqs. (35,36) we see that the forms of Yukawa couplings are determined by the modular weights 2k F + k 10 and 2k F + k 126 , and they depend on a number of coupling constants and the complex modulus τ . In this work, we will be concerned with the modular forms of level 3 up to weight 10, and higher weight modular forms can be discussed in a similar way although more modular invariant contractions as well as couplings would be involved. Thus the possible values of 2k F + k 10 and 2k F + k 126 are 0, 2, 4, 6, 8, 10 up to weight 10 modular forms, and consequently we can obtain 6 × 6 = 36 minimal SO(10) GUT models.
We numerically analyze all these models and minimize the corresponding χ 2 function with the TMinuit package [81]. We find that only two out of the 36 models can give a good fit to the experimental data for certain values of input parameters. Nevertheless some minimal models with less parameters can achieve a relatively small χ 2 min , and only the quark mass ratio m d /m s lies outside the 3σ allowed region. It's reasonable to regard these models as a good leading order approximation, we would like to give such an example model named as minimal model 1. This model is specified by the sum of modular weights (2k F + k 10 , 2k F + k 126 ) = (10, 6) and the corresponding superpotential is given in Eq. (39). We find that if both types of seesaw are present, almost all fermion observables can be correctly reproduced except the m d /m s ratio. The fitting results of the minimal model 1 is summarized in table 3, where we split the χ 2 into three parts χ 2 = χ 2 l + χ 2 q + χ 2 bτ and χ 2 l , χ 2 q and χ 2 bτ denote the pieces of χ 2 function arising from the deviations of the lepton sector observables, the quark sector observables and m b /m τ from their central values respectively. We see that m d /m s is predicted to be small and it is about 5σ away from the experimental best fit value while all other observables are within the 3σ ranges of the experimental data.
The good agreement with data can be achieved with higher modular weight 2k F + k 126 but at the price of introducing more free parameters in the model. The two phenomenologically viable minimal models are named as minimal model 2 and minimal model 3 with modular weights (2k F + k 10 , 2k F + k 126 ) equal to (10,8) and (10,10) respectively. The corresponding forms of the Yukawa superpotentials are given in Eq. (40) and Eq. (41). Comprehensive numerical analysis shows that minimal model 2 and minimal model 3 can accommodate the experimental data if both type-I and type-II seesaw mechanisms contribute to the light neutrino masses. There are total 27 free real parameters in fermion mass matrices in both models. Our fits yield χ 2 = 8.59476 for the minimal model 2 and χ 2 = 3.60588 for the minimal model 3. We tabulate the best fit values of the input parameters and the predictions for the fermion observables in table 3. It is notable that almost all the measured observables fall within the 1σ experimentally allowed ranges for these two viable SO(10) minimal models. Moreover, we give the model predictions for several as yet unmeasured observables including the Majorana CP violation phases α 21 , α 31 , the lightest neutrino mass m 1 , the masses of the heavy right-handed neutrinos M 1,2,3 and the effective mass m ββ in neutrinoless double beta decay, they can be understood as predictions of the models. The effective Majorana neutrino mass m ββ are defined as m ββ = m 1 cos 2 θ l 12 cos 2 θ l 13 + m 2 sin 2 θ l 12 cos 2 θ l 13 e iα 21 + m 3 sin 2 θ l 13 e i(α 31 −2δ l CP ) .
The predicted value of m ββ at the best-fit point in the minimal model 2 and minimal model 3 are m ββ = 0.738522 meV and m ββ = 0.667982 meV respectively, which are far below the sensitivity of future tonne-scale neutrinoless double beta decay experiments. We can use the measured value of the solar neutrino mass squared ∆m 2 21 = 7.42×10 −5 eV 2 [80] to fix the overall mass scale α 1 v 2 u /v R of neutrino mass matrix, subsequently the absolute values of light neutrino masses can be determined as follows, The seesaw scale is of order 10 10 ∼ 10 11 GeV in the two models.
Furthermore, we use the widely-used sampling program MultiNest [83,84] to scan the parameter space around the best fit points, and the predictions for the fermion masses and mixing parameters are required to be compatible with data at 3σ level. Notice that different predictions could possibly be obtained around other local minima. The allowed values of the modulus τ and the correlations between observables for the minimal model 2 and minimal model 3 are plotted in figure 1 and  figure 2 respectively. We can see that the regions of τ compatible with data are very narrow in the two minimal models, the effective Majorana mass m ββ characterizing the neutrinoless double beta decay amplitude lies in the narrow intervals around 0.7 meV which is too small to be detected. Although both minimal model 2 and minimal model 3 can give good accommodation to the experimental results, we note that a substantial amount of fine tuning of the free parameters is needed. The fine-tuning parameter d F T defined in Eq. (46) is found to be of order d F T ≈ 10 6 in both minimal model 2 and minimal model 3, to be compared with d Data = 531.986.

Numerical results of the non-minimal models
In this section, we proceed to consider the numerical results of the non-minimal models with the full Higgs content H, ∆ and Σ which are 10, 126 and 120 dimensional multiplets of SO (10). For the non-minimal SO(10) models, the general form of the Yukawa matrices for different fermions and the right-handed neutrino mass matrix are given in Eq. (33). In comparison with the minimal SO(10) models, three more Higgs mixing parameters r 3 , c e and c ν accompanying the Yukawa coupling Y 120 of Σ are present. The SO(10) non-minimal models are completely specified by the modular weights 2k F + k 10 , 2k F + k 120 and 2k F + k 126 . Hence there are total 6 × 6 × 5 = 180 possible non-minimal models if one concerns with the modular forms up to weight 10.
We scan the parameter space of the above non-minimal models one by one. In contrast with the minimal models, there are plenty of non-minimal models which are compatible with the experimental data of fermion masses and mixings even in either type-I or type-II dominated seesaw mechanisms. It is too lengthy to list all the viable models here. In the following, we present three benchmark non-minimal models which are characterized by Non-minimal model 1: Non-minimal model 2: (2k F + k 10 , 2k F + k 120 , 2k F + k 126 ) = (4, 2, 4) , Non-minimal model 3: (2k F + k 10 , 2k F + k 120 , 2k F + k 126 ) = (4, 6, 4) .
The superpotentials can be straightforwardly read out as given in Eq. (42), Eq. (43) and Eq. (44) respectively. We report the fitting results of the above three non-minimal models in table 4, table 5  and table 6 respectively. The non-minimal model 1 can give a reasonably good fit to all observables if the neutrino masses are generated by the type-I seesaw and the mixture of type-I and type-II seesaw which are denoted as SS-I and SS-I+II respectively. It turns out that this model is unable to reproduce the fermion masses and mixing for the pure type-II seesaw case, and accordingly the minima of χ 2 from the numerical optimization algorithm is quite large. In the type-I seesaw dominant scenario, 23 real free parameters are involved and the main contribution to χ 2 min comes from the mass ratio m b /m τ which is about 3σ larger than its central value, and all other observables are reproduced correctly with small deviations as shown in the table 4. The light neutrino masses for the type-I and type-I+II cases are predicted to be for the non-minimal model 1.
The allowed values of τ and the correlations between observables are displayed in figure 3. The combination of type-I and type-II seesaw mechanism can give a better fit of the data but at the price of introducing one more complex free parameter v L . The smallest χ 2 is found to be χ 2 min = 9.98822 and among all observables the largest contributions to χ 2 is χ 2 bτ = 3.10314 from m b /m τ . The corresponding correlations between predictions are plotted in figure 4. We can see from figure 4 that the allowed regions of lepton mixing angles θ l 12 , θ l 23 and CP-violation phases δ l CP , α 21 as well as α 31 are quite narrow. The effective Majorana mass m ββ is predicted to be around 1.9 meV and it is far below the sensitivities of next generation neutrinoless double beta decay experiments.
The numerical fitting results of the non-minimal model 2 are summarized in table 5, and this model can match the measured values of observables given in table 2 for all the three types of neutrino mass generation mechanism. The mass ratio m b /m τ is about 3.5σ away from its mean value in the type-II seesaw dominant scenario, while all other observables can be explained well with small derivations. The issue about the prediction of m b /m τ can be well resolved if the type-I seesaw mechanism dominates in neutrino sector, and the minimum of χ 2 is found to be χ 2 min = 2.71438 with χ 2 bτ = 0.32689. The best fits values of fermion masses and the mixing parameters are in good agreement with the corresponding inputted reference values. The fitting results can be slightly improved in the mixture of type-I and type-II seesaw mechanisms. The overall scale of the neutrino masses is fixed by the solar neutrino mass difference ∆m Notice that the neutrino mass spectrum is strongly hierarchical in the case of type-I and type-I+II seesaw, while it is quasi-degenerate in the pure type-II seesaw case. For the pure type-II seesaw case, the light neutrino masses are quasi-degenerate and the neutrino mass sum is i m i = 218.004 meV which is above the aggressive upper bound 120meV but still below the conservative upper limit 600 meV given by the Planck collaboration [85]. Notice that the cosmological bound on the neutrino masses significantly depend on the data sets which are combined to break the degeneracies of the many cosmological parameters, and the upper limit on neutrino mass becomes weaker when one departs from the framework of ΛCDM plus neutrino mass to frameworks with more cosmological parameters. Three heavy right-handed neutrinos are introduced in the type-I seesaw mechanism, and the scale of the right-handed neutrino mass are predicted to be 10 10 ∼ 10 13 GeV in the type-I and type-I+II cases. We show the allowed values of τ and the correlations between observables for SS-I, SS-II and SS-I+II in figure 5, figure 6 and figure 7 respectively. Notice that the possible values of certain observables lie in rather small regions.
The non-minimal model 3 differs from the non-minimal model 2 in the value of 2k F + k 120 . From Eq. (36) we know that there are two independent A 4 invariant terms in the antisymmetric Yukawa coupling Y 120 for 2k F + k 120 = 6 in the non-minimal model 3 while only a single term is allowed for 2k F + k 120 = 2 as in non-minimal model 2. As a consequence, the non-minimal model 3 will have one more complex coupling than the non-minimal model 2. The fitting results are given in table 6, we see that all three types of seesaw mechanisms can give quite good fits to the fermion observables. The minimal χ 2 min = 5.61610 in pure type-II case is smaller than the χ 2 min = 7.02930 in the type-I seesaw dominant case. Therefore the type-II seesaw can explain the data a bit better than the type-I seesaw. The predictions of light neutrino masses for three seesaw cases are Analogous to the non-minimal model 2, the light neutrino masses are also quasi-degenerate for type-II seesaw dominance and they are consistent with the conservative limit on neutrino masses from Planck [85].
We display the allowed values of τ and the correlations among the fermion masses and mixing parameters for the non-minimal model 3 with type-I, type-II and type-I+II seesaw in figure 8, figure 9 and figure 10 respectively. From figure 8 we see that the values of the complex modulus τ scatter in a small region close to the boundary of the fundamental region. The atmospheric mixing angle θ l 23 varies in the second octant and the allowed region of δ l CP is about [π, 1.5π]. The effective Majorana neutrino mass m ββ is found to be about 1 meV with the lightest neutrino mass m 1 around 6 meV. If the contribution of the type-II seesaw dominates over the other, θ l 23 is predicted to lie in the experimental favored 1σ region [80], as shown in figure 9. In the scenario of the mixed type-I and type-II seesaw mechanisms, from figure 10 we can see that the allowed region of τ is close to the residual symmetry preserved point −1/2 + √ 3i/2, the atmospheric mixing angle θ l 23 is less constrained and its whole 3σ range can be covered, and the predicted value of m ββ is around 1 meV.

Conclusion
In this paper we have combined SO(10) Grand Unified Theories (GUTs) with Γ 3 A 4 modular symmetry and presented a comprehensive analysis of the resulting quark and lepton mass matrices for all the simplest cases. We have focussed on the case where the three fermion families in the 16 dimensional spinor representation form a triplet of Γ 3 A 4 , with a Higgs sector comprising a single Higgs multiplet H in the 10 fundamental representation and one Higgs field ∆ in the 126 for the minimal models, plus an additional Higgs field Σ in the 120 for the non-minimal models, all with specified modular weights. The models are completely specified by the summation of the modular weights of matter fields and Higgs fields. The neutrino masses are generated by the type-I and/or type II seesaw mechanisms and results are presented for each model following an intensive numerical analysis where we have optimized the free parameters of the models in order to match the experimental data. For the phenomenologically successful models, we present the best fit results in numerical tabular form as well as showing the most interesting graphical correlations between parameters, including leptonic CP phases and neutrinoless double beta decay, which have yet to be measured, leading to definite predictions for each of the models.
Once the modular weights are specified, the Yukawa couplings are determined up to a number of overall dimensionless complex coefficients, and the value of the single complex modulus field τ , which is the only flavon in the theory. All models with sums of modular weights up to 10 were considered, and we presented results only for the simplest phenomenologically viable models. We found that the minimal models containing only the Higgs fields in the 10 and the 126 are the least viable, requiring both type I and type II seesaw mechanisms to be present simultaneously and also necessitating at least some sums of modular weights of 10 (if the sums of modular weights were restricted to 8 then no viable such models were found). On the other hand, non-minimal models involving in addition to the Higgs fields in the 10 and the 126, also a Higgs field in the 120, proved to be more successful, with many such models being found with sums of modular weights of up to 8 or less, and with the type I seesaw as well as the combined type I+II seesaw also being viable. In the minimal models, more free parameters in the Yukawa couplings are necessary accommodate the experimental data of fermion masses and mixing, while fewer parameters are required in the Yukawa superpotential of non-minimal models, but the Higgs sector is more complicated due to the presence of Σ in the representation 120. The Higgs potential required for the two light MSSM Higgs doublets to emerge from components of the 10, 126 and 120 SO(10) multiplets, leaving the other components heavy (the doublet-triplet splitting problem) was not considered.
In conclusion, we have successfully combined SO(10) GUTs with Γ 3 A 4 modular symmetry and analysed the simplest models of this kind, presenting the best fit results in numerical tabular form as well as showing the most interesting graphical correlations between parameters, leading to definite predictions for each of the models. The right-handed neutrino masses are predicted to be in the typical range for leptogenesis, and it would be interesting to study this in a future publication.      Table 6: The best fit values of the free parameters and the corresponding predictions for the masses and mixing parameters of lepton and quark mixing at the best fit point in the SO(10) non-minimal model 3 specified by the modular weights (2kF + k10, 2kF + k120, 2kF + k 126 ) = (4, 6, 4). Figure 1: The values of the complex modulus τ compatible with experimental data and the correlations between the neutrino mixing angles, CP violation phases, quark mass ratios and mixing parameters in the minimal model 2 with SS-I+II. The vertical and horizontal dashed lines are the 3σ bounds taken from [80]. In the panel of m ββ with respect to the lightest neutrino mass mmin, the blue (red) dashed lines stand for the most general allowed regions for normal ordering (inverted ordering) neutrino mass spectrum respectively, where the neutrino oscillation parameter are varied within their 3σ ranges. The present upper limit m ββ < (61 − 165) meV from KamLAND-Zen [86] is shown by horizontal grey band. The vertical grey exclusion band represents the most aggressive bound i mi < 0.120eV from Planck [85].         with SS-I+II. The vertical and horizontal dashed lines are the 3σ bounds taken from [80].